|
Packit |
67cb25 |
/* specfunc/recurse.h
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Author: G. Jungman */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#ifndef _RECURSE_H_
|
|
Packit |
67cb25 |
#define _RECURSE_H_
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#define CONCAT(a,b) a ## _ ## b
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* n_max >= n_min + 2
|
|
Packit |
67cb25 |
* f[n+1] + a[n] f[n] + b[n] f[n-1] = 0
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Trivial forward recurrence.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
#define GEN_RECURSE_FORWARD_SIMPLE(func) \
|
|
Packit |
67cb25 |
int CONCAT(recurse_forward_simple, func) ( \
|
|
Packit |
67cb25 |
const int n_max, const int n_min, \
|
|
Packit |
67cb25 |
const double parameters[], \
|
|
Packit |
67cb25 |
const double f_n_min, \
|
|
Packit |
67cb25 |
const double f_n_min_p1, \
|
|
Packit |
67cb25 |
double * f, \
|
|
Packit |
67cb25 |
double * f_n_max \
|
|
Packit |
67cb25 |
) \
|
|
Packit |
67cb25 |
{ \
|
|
Packit |
67cb25 |
int n; \
|
|
Packit |
67cb25 |
\
|
|
Packit |
67cb25 |
if(f == 0) { \
|
|
Packit |
67cb25 |
double f2 = f_n_min; \
|
|
Packit |
67cb25 |
double f1 = f_n_min_p1; \
|
|
Packit |
67cb25 |
double f0; \
|
|
Packit |
67cb25 |
for(n=n_min+2; n<=n_max; n++) { \
|
|
Packit |
67cb25 |
f0 = -REC_COEFF_A(n-1,parameters) * f1 - REC_COEFF_B(n-1, parameters) * f2; \
|
|
Packit |
67cb25 |
f2 = f1; \
|
|
Packit |
67cb25 |
f1 = f0; \
|
|
Packit |
67cb25 |
} \
|
|
Packit |
67cb25 |
*f_n_max = f0; \
|
|
Packit |
67cb25 |
} \
|
|
Packit |
67cb25 |
else { \
|
|
Packit |
67cb25 |
f[n_min] = f_n_min; \
|
|
Packit |
67cb25 |
f[n_min + 1] = f_n_min_p1; \
|
|
Packit |
67cb25 |
for(n=n_min+2; n<=n_max; n++) { \
|
|
Packit |
67cb25 |
f[n] = -REC_COEFF_A(n-1,parameters) * f[n-1] - REC_COEFF_B(n-1, parameters) * f[n-2]; \
|
|
Packit |
67cb25 |
} \
|
|
Packit |
67cb25 |
*f_n_max = f[n_max]; \
|
|
Packit |
67cb25 |
} \
|
|
Packit |
67cb25 |
\
|
|
Packit |
67cb25 |
return GSL_SUCCESS; \
|
|
Packit |
67cb25 |
} \
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* n_start >= n_max >= n_min
|
|
Packit |
67cb25 |
* f[n+1] + a[n] f[n] + b[n] f[n-1] = 0
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Generate the minimal solution of the above recursion relation,
|
|
Packit |
67cb25 |
* with the simplest form of the normalization condition, f[n_min] given.
|
|
Packit |
67cb25 |
* [Gautschi, SIAM Rev. 9, 24 (1967); (3.9) with s[n]=0]
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
#define GEN_RECURSE_BACKWARD_MINIMAL_SIMPLE(func) \
|
|
Packit |
67cb25 |
int CONCAT(recurse_backward_minimal_simple, func) ( \
|
|
Packit |
67cb25 |
const int n_start, \
|
|
Packit |
67cb25 |
const int n_max, const int n_min, \
|
|
Packit |
67cb25 |
const double parameters[], \
|
|
Packit |
67cb25 |
const double f_n_min, \
|
|
Packit |
67cb25 |
double * f, \
|
|
Packit |
67cb25 |
double * f_n_max \
|
|
Packit |
67cb25 |
) \
|
|
Packit |
67cb25 |
{ \
|
|
Packit |
67cb25 |
int n; \
|
|
Packit |
67cb25 |
double r_n = 0.; \
|
|
Packit |
67cb25 |
double r_nm1; \
|
|
Packit |
67cb25 |
double ratio; \
|
|
Packit |
67cb25 |
\
|
|
Packit |
67cb25 |
for(n=n_start; n > n_max; n--) { \
|
|
Packit |
67cb25 |
r_nm1 = -REC_COEFF_B(n, parameters) / (REC_COEFF_A(n, parameters) + r_n); \
|
|
Packit |
67cb25 |
r_n = r_nm1; \
|
|
Packit |
67cb25 |
} \
|
|
Packit |
67cb25 |
\
|
|
Packit |
67cb25 |
if(f != 0) { \
|
|
Packit |
67cb25 |
f[n_max] = 10.*DBL_MIN; \
|
|
Packit |
67cb25 |
for(n=n_max; n > n_min; n--) { \
|
|
Packit |
67cb25 |
r_nm1 = -REC_COEFF_B(n, parameters) / (REC_COEFF_A(n, parameters) + r_n); \
|
|
Packit |
67cb25 |
f[n-1] = f[n] / r_nm1; \
|
|
Packit |
67cb25 |
r_n = r_nm1; \
|
|
Packit |
67cb25 |
} \
|
|
Packit |
67cb25 |
ratio = f_n_min / f[n_min]; \
|
|
Packit |
67cb25 |
for(n=n_min; n<=n_max; n++) { \
|
|
Packit |
67cb25 |
f[n] *= ratio; \
|
|
Packit |
67cb25 |
} \
|
|
Packit |
67cb25 |
} \
|
|
Packit |
67cb25 |
else { \
|
|
Packit |
67cb25 |
double f_nm1; \
|
|
Packit |
67cb25 |
double f_n = 10.*DBL_MIN; \
|
|
Packit |
67cb25 |
*f_n_max = f_n; \
|
|
Packit |
67cb25 |
for(n=n_max; n > n_min; n--) { \
|
|
Packit |
67cb25 |
r_nm1 = -REC_COEFF_B(n, parameters) / (REC_COEFF_A(n, parameters) + r_n); \
|
|
Packit |
67cb25 |
f_nm1 = f_n / r_nm1; \
|
|
Packit |
67cb25 |
r_n = r_nm1; \
|
|
Packit |
67cb25 |
} \
|
|
Packit |
67cb25 |
ratio = f_n_min / f_nm1; \
|
|
Packit |
67cb25 |
*f_n_max *= ratio; \
|
|
Packit |
67cb25 |
} \
|
|
Packit |
67cb25 |
\
|
|
Packit |
67cb25 |
return GSL_SUCCESS; \
|
|
Packit |
67cb25 |
} \
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#endif /* !_RECURSE_H_ */
|