recurse.h

开放gsl矩阵运算

/* 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., 675 Mass Ave, Cambridge, MA 02139, 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_ */