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