htswanalysis/MACS/lib/gsl/gsl-1.11/specfunc/recurse.h

   1 /* specfunc/recurse.h
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* Author: G. Jungman */
  21
  22 #ifndef _RECURSE_H_
  23 #define _RECURSE_H_
  24
  25 #define CONCAT(a,b) a ## _ ## b
  26
  27
  28 /* n_max >= n_min + 2
  29  * f[n+1] + a[n] f[n] + b[n] f[n-1] = 0
  30  *
  31  * Trivial forward recurrence.
  32  */
  33 #define GEN_RECURSE_FORWARD_SIMPLE(func)                                      \
  34 int CONCAT(recurse_forward_simple, func) (                                    \
  35                                const int n_max, const int n_min,              \
  36                                const double parameters[],                     \
  37                                const double f_n_min,                          \
  38                                const double f_n_min_p1,                       \
  39                                double * f,                                    \
  40                                double * f_n_max                               \
  41                                )                                              \
  42 {                                                                             \
  43   int n;                                                                      \
  44                                                                               \
  45   if(f == 0) {                                                                \
  46     double f2 = f_n_min;                                                      \
  47     double f1 = f_n_min_p1;                                                   \
  48     double f0;                                                                \
  49     for(n=n_min+2; n<=n_max; n++) {                                           \
  50       f0 = -REC_COEFF_A(n-1,parameters) * f1 - REC_COEFF_B(n-1, parameters) * f2; \
  51       f2 = f1;                                                                \
  52       f1 = f0;                                                                \
  53     }                                                                         \
  54     *f_n_max = f0;                                                            \
  55   }                                                                           \
  56   else {                                                                      \
  57     f[n_min]     = f_n_min;                                                   \
  58     f[n_min + 1] = f_n_min_p1;                                                \
  59     for(n=n_min+2; n<=n_max; n++) {                                           \
  60       f[n] = -REC_COEFF_A(n-1,parameters) * f[n-1] - REC_COEFF_B(n-1, parameters) * f[n-2]; \
  61     }                                                                         \
  62     *f_n_max = f[n_max];                                                      \
  63   }                                                                           \
  64                                                                               \
  65   return GSL_SUCCESS;                                                         \
  66 }                                                                             \
  67
  68
  69 /* n_start >= n_max >= n_min
  70  * f[n+1] + a[n] f[n] + b[n] f[n-1] = 0
  71  *
  72  * Generate the minimal solution of the above recursion relation,
  73  * with the simplest form of the normalization condition, f[n_min] given.
  74  * [Gautschi, SIAM Rev. 9, 24 (1967); (3.9) with s[n]=0]
  75  */
  76 #define GEN_RECURSE_BACKWARD_MINIMAL_SIMPLE(func)                             \
  77 int CONCAT(recurse_backward_minimal_simple, func) (                           \
  78                                const int n_start,                             \
  79                                const int n_max, const int n_min,              \
  80                                const double parameters[],                     \
  81                                const double f_n_min,                          \
  82                                double * f,                                    \
  83                                double * f_n_max                               \
  84                                )                                              \
  85 {                                                                             \
  86   int n;                                                                      \
  87   double r_n = 0.;                                                            \
  88   double r_nm1;                                                               \
  89   double ratio;                                                               \
  90                                                                               \
  91   for(n=n_start; n > n_max; n--) {                                            \
  92     r_nm1 = -REC_COEFF_B(n, parameters) / (REC_COEFF_A(n, parameters) + r_n); \
  93     r_n = r_nm1;                                                              \
  94   }                                                                           \
  95                                                                               \
  96   if(f != 0) {                                                                \
  97     f[n_max] = 10.*DBL_MIN;                                                      \
  98     for(n=n_max; n > n_min; n--) {                                               \
  99       r_nm1  = -REC_COEFF_B(n, parameters) / (REC_COEFF_A(n, parameters) + r_n); \
 100       f[n-1] = f[n] / r_nm1;                                                     \
 101       r_n = r_nm1;                                                               \
 102     }                                                                         \
 103     ratio = f_n_min / f[n_min];                                               \
 104     for(n=n_min; n<=n_max; n++) {                                             \
 105       f[n] *= ratio;                                                          \
 106     }                                                                         \
 107   }                                                                           \
 108   else {                                                                      \
 109     double f_nm1;                                                             \
 110     double f_n = 10.*DBL_MIN;                                                 \
 111     *f_n_max = f_n;                                                           \
 112     for(n=n_max; n > n_min; n--) {                                               \
 113       r_nm1 = -REC_COEFF_B(n, parameters) / (REC_COEFF_A(n, parameters) + r_n);  \
 114       f_nm1 = f_n / r_nm1;                                                       \
 115       r_n = r_nm1;                                                               \
 116     }                                                                         \
 117     ratio = f_n_min / f_nm1;                                                  \
 118     *f_n_max *= ratio;                                                        \
 119   }                                                                           \
 120                                                                               \
 121   return GSL_SUCCESS;                                                         \
 122 }                                                                             \
 123
 124
 125 #endif /* !_RECURSE_H_ */