home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / hyperg.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  8.6 KB  |  285 lines
  1. /* specfunc/hyperg.c
  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 2 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., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. /* Miscellaneous implementations of use
  23.  * for evaluation of hypergeometric functions.
  24.  */
  25. #include <config.h>
  26. #include <gsl/gsl_math.h>
  27. #include <gsl/gsl_errno.h>
  28. #include <gsl/gsl_sf_exp.h>
  29. #include <gsl/gsl_sf_gamma.h>
  30.  
  31. #include "error.h"
  32. #include "hyperg.h"
  33.  
  34. #define SUM_LARGE  (1.0e-5*GSL_DBL_MAX)
  35.  
  36.  
  37. int
  38. gsl_sf_hyperg_1F1_series_e(const double a, const double b, const double x,
  39.                               gsl_sf_result * result
  40.                               )
  41. {
  42.   double an  = a;
  43.   double bn  = b;
  44.   double n   = 1.0;
  45.   double del = 1.0;
  46.   double abs_del = 1.0;
  47.   double max_abs_del = 1.0;
  48.   double sum_val = 1.0;
  49.   double sum_err = 0.0;
  50.  
  51.   while(abs_del/fabs(sum_val) > GSL_DBL_EPSILON) {
  52.     double u, abs_u;
  53.  
  54.     if(bn == 0.0) {
  55.       DOMAIN_ERROR(result);
  56.     }
  57.     if(an == 0.0 || n > 1000.0) {
  58.       result->val  = sum_val;
  59.       result->err  = sum_err;
  60.       result->err += 2.0 * GSL_DBL_EPSILON * n * fabs(sum_val);
  61.       return GSL_SUCCESS;
  62.     }
  63.  
  64.     u = x * (an/(bn*n));
  65.     abs_u = fabs(u);
  66.     if(abs_u > 1.0 && max_abs_del > GSL_DBL_MAX/abs_u) {
  67.       result->val = sum_val;
  68.       result->err = fabs(sum_val);
  69.       GSL_ERROR ("overflow", GSL_EOVRFLW);
  70.     }
  71.     del *= u;
  72.     sum_val += del;
  73.     if(fabs(sum_val) > SUM_LARGE) {
  74.       result->val = sum_val;
  75.       result->err = fabs(sum_val);
  76.       GSL_ERROR ("overflow", GSL_EOVRFLW);
  77.     }
  78.  
  79.     abs_del = fabs(del);
  80.     max_abs_del = GSL_MAX_DBL(abs_del, max_abs_del);
  81.     sum_err += 2.0*GSL_DBL_EPSILON*abs_del;
  82.  
  83.     an += 1.0;
  84.     bn += 1.0;
  85.     n  += 1.0;
  86.   }
  87.  
  88.   result->val  = sum_val;
  89.   result->err  = sum_err;
  90.   result->err += abs_del;
  91.   result->err += 2.0 * GSL_DBL_EPSILON * n * fabs(sum_val);
  92.  
  93.   return GSL_SUCCESS;
  94. }
  95.  
  96.  
  97. int
  98. gsl_sf_hyperg_1F1_large_b_e(const double a, const double b, const double x, gsl_sf_result * result)
  99. {
  100.   if(fabs(x/b) < 1.0) {
  101.     const double u = x/b;
  102.     const double v = 1.0/(1.0-u);
  103.     const double pre = pow(v,a);
  104.     const double uv  = u*v;
  105.     const double uv2 = uv*uv;
  106.     const double t1  = a*(a+1.0)/(2.0*b)*uv2;
  107.     const double t2a = a*(a+1.0)/(24.0*b*b)*uv2;
  108.     const double t2b = 12.0 + 16.0*(a+2.0)*uv + 3.0*(a+2.0)*(a+3.0)*uv2;
  109.     const double t2  = t2a*t2b;
  110.     result->val  = pre * (1.0 - t1 + t2);
  111.     result->err  = pre * GSL_DBL_EPSILON * (1.0 + fabs(t1) + fabs(t2));
  112.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  113.     return GSL_SUCCESS;
  114.   }
  115.   else {
  116.     DOMAIN_ERROR(result);
  117.   }
  118. }
  119.  
  120.  
  121. int
  122. gsl_sf_hyperg_U_large_b_e(const double a, const double b, const double x,
  123.                              gsl_sf_result * result,
  124.                  double * ln_multiplier
  125.                  )
  126. {
  127.   double N   = floor(b);  /* b = N + eps */
  128.   double eps = b - N;
  129.   
  130.   if(fabs(eps) < GSL_SQRT_DBL_EPSILON) {
  131.     double lnpre_val;
  132.     double lnpre_err;
  133.     gsl_sf_result M;
  134.     if(b > 1.0) {
  135.       double tmp = (1.0-b)*log(x);
  136.       gsl_sf_result lg_bm1;
  137.       gsl_sf_result lg_a;
  138.       gsl_sf_lngamma_e(b-1.0, &lg_bm1);
  139.       gsl_sf_lngamma_e(a, &lg_a);
  140.       lnpre_val = tmp + x + lg_bm1.val - lg_a.val;
  141.       lnpre_err = lg_bm1.err + lg_a.err + GSL_DBL_EPSILON * (fabs(x) + fabs(tmp));
  142.       gsl_sf_hyperg_1F1_large_b_e(1.0-a, 2.0-b, -x, &M);
  143.     }
  144.     else {
  145.       gsl_sf_result lg_1mb;
  146.       gsl_sf_result lg_1pamb;
  147.       gsl_sf_lngamma_e(1.0-b, &lg_1mb);
  148.       gsl_sf_lngamma_e(1.0+a-b, &lg_1pamb);
  149.       lnpre_val = lg_1mb.val - lg_1pamb.val;
  150.       lnpre_err = lg_1mb.err + lg_1pamb.err;
  151.       gsl_sf_hyperg_1F1_large_b_e(a, b, x, &M);
  152.     }
  153.  
  154.     if(lnpre_val > GSL_LOG_DBL_MAX-10.0) {
  155.       result->val  = M.val;
  156.       result->err  = M.err;
  157.       *ln_multiplier = lnpre_val;
  158.       GSL_ERROR ("overflow", GSL_EOVRFLW);
  159.     }
  160.     else {
  161.       gsl_sf_result epre;
  162.       int stat_e = gsl_sf_exp_err_e(lnpre_val, lnpre_err, &epre);
  163.       result->val  = epre.val * M.val;
  164.       result->err  = epre.val * M.err + epre.err * fabs(M.val);
  165.       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  166.       *ln_multiplier = 0.0;
  167.       return stat_e;
  168.     }
  169.   }
  170.   else {
  171.     double omb_lnx = (1.0-b)*log(x);
  172.     gsl_sf_result lg_1mb;    double sgn_1mb;
  173.     gsl_sf_result lg_1pamb;  double sgn_1pamb;
  174.     gsl_sf_result lg_bm1;    double sgn_bm1;
  175.     gsl_sf_result lg_a;      double sgn_a;
  176.     gsl_sf_result M1, M2;
  177.     double lnpre1_val, lnpre2_val;
  178.     double lnpre1_err, lnpre2_err;
  179.     double sgpre1, sgpre2;
  180.     gsl_sf_hyperg_1F1_large_b_e(    a,     b, x, &M1);
  181.     gsl_sf_hyperg_1F1_large_b_e(1.0-a, 2.0-b, x, &M2);
  182.  
  183.     gsl_sf_lngamma_sgn_e(1.0-b,   &lg_1mb,   &sgn_1mb);
  184.     gsl_sf_lngamma_sgn_e(1.0+a-b, &lg_1pamb, &sgn_1pamb);
  185.  
  186.     gsl_sf_lngamma_sgn_e(b-1.0, &lg_bm1, &sgn_bm1);
  187.     gsl_sf_lngamma_sgn_e(a,     &lg_a,   &sgn_a);
  188.  
  189.     lnpre1_val = lg_1mb.val - lg_1pamb.val;
  190.     lnpre1_err = lg_1mb.err + lg_1pamb.err;
  191.     lnpre2_val = lg_bm1.val - lg_a.val - omb_lnx - x;
  192.     lnpre2_err = lg_bm1.err + lg_a.err + GSL_DBL_EPSILON * (fabs(omb_lnx)+fabs(x));
  193.     sgpre1 = sgn_1mb * sgn_1pamb;
  194.     sgpre2 = sgn_bm1 * sgn_a;
  195.  
  196.     if(lnpre1_val > GSL_LOG_DBL_MAX-10.0 || lnpre2_val > GSL_LOG_DBL_MAX-10.0) {
  197.       double max_lnpre_val = GSL_MAX(lnpre1_val,lnpre2_val);
  198.       double max_lnpre_err = GSL_MAX(lnpre1_err,lnpre2_err);
  199.       double lp1 = lnpre1_val - max_lnpre_val;
  200.       double lp2 = lnpre2_val - max_lnpre_val;
  201.       double t1  = sgpre1*exp(lp1);
  202.       double t2  = sgpre2*exp(lp2);
  203.       result->val  = t1*M1.val + t2*M2.val;
  204.       result->err  = fabs(t1)*M1.err + fabs(t2)*M2.err;
  205.       result->err += GSL_DBL_EPSILON * exp(max_lnpre_err) * (fabs(t1*M1.val) + fabs(t2*M2.val));
  206.       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  207.       *ln_multiplier = max_lnpre_val;
  208.       GSL_ERROR ("overflow", GSL_EOVRFLW);
  209.     }
  210.     else {
  211.       double t1 = sgpre1*exp(lnpre1_val);
  212.       double t2 = sgpre2*exp(lnpre2_val);
  213.       result->val  = t1*M1.val + t2*M2.val;
  214.       result->err  = fabs(t1) * M1.err + fabs(t2)*M2.err;
  215.       result->err += GSL_DBL_EPSILON * (exp(lnpre1_err)*fabs(t1*M1.val) + exp(lnpre2_err)*fabs(t2*M2.val));
  216.       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  217.       *ln_multiplier = 0.0;
  218.       return GSL_SUCCESS;
  219.     }
  220.   }
  221. }
  222.  
  223.  
  224.  
  225. /* [Carlson, p.109] says the error in truncating this asymptotic series
  226.  * is less than the absolute value of the first neglected term.
  227.  *
  228.  * A termination argument is provided, so that the series will
  229.  * be summed at most up to n=n_trunc. If n_trunc is set negative,
  230.  * then the series is summed until it appears to start diverging.
  231.  */
  232. int
  233. gsl_sf_hyperg_2F0_series_e(const double a, const double b, const double x,
  234.                               int n_trunc,
  235.                               gsl_sf_result * result
  236.                               )
  237. {
  238.   const int maxiter = 2000;
  239.   double an = a;
  240.   double bn = b;  
  241.   double n   = 1.0;
  242.   double sum = 1.0;
  243.   double del = 1.0;
  244.   double abs_del = 1.0;
  245.   double max_abs_del = 1.0;
  246.   double last_abs_del = 1.0;
  247.   
  248.   while(abs_del/fabs(sum) > GSL_DBL_EPSILON && n < maxiter) {
  249.  
  250.     double u = an * (bn/n * x);
  251.     double abs_u = fabs(u);
  252.  
  253.     if(abs_u > 1.0 && (max_abs_del > GSL_DBL_MAX/abs_u)) {
  254.       result->val = sum;
  255.       result->err = fabs(sum);
  256.       GSL_ERROR ("overflow", GSL_EOVRFLW);
  257.     }
  258.  
  259.     del *= u;
  260.     sum += del;
  261.  
  262.     abs_del = fabs(del);
  263.  
  264.     if(abs_del > last_abs_del) break; /* series is probably starting to grow */
  265.  
  266.     last_abs_del = abs_del;
  267.     max_abs_del  = GSL_MAX(abs_del, max_abs_del);
  268.  
  269.     an += 1.0;
  270.     bn += 1.0;
  271.     n  += 1.0;
  272.     
  273.     if(an == 0.0 || bn == 0.0) break;        /* series terminated */
  274.     
  275.     if(n_trunc >= 0 && n >= n_trunc) break;  /* reached requested timeout */
  276.   }
  277.  
  278.   result->val = sum;
  279.   result->err = GSL_DBL_EPSILON * n + abs_del;
  280.   if(n >= maxiter)
  281.     GSL_ERROR ("error", GSL_EMAXITER);
  282.   else
  283.     return GSL_SUCCESS;
  284. }
  285.