home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 41 / IOPROG_41.ISO / soft / c++ / NUMCPP11.ZIP / mathx.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  1999-05-09  |  9.7 KB  |  375 lines
  1. //===================================================================
  2. // mathx.h
  3. //
  4. // Version 1.1
  5. //
  6. // Written by:
  7. //   Brent Worden
  8. //   WordenWare
  9. //   email:  Brent@Worden.org
  10. //
  11. // Copyright (c) 1998-1999 WordenWare
  12. //
  13. // Created:  August 28, 1998
  14. // Revised:  April 10, 1999
  15. //===================================================================
  16.  
  17. #include <cmath>
  18.  
  19. #include "mathx.h"
  20. #include "numerror.h"
  21. #include "utility.hpp"
  22.  
  23. NUM_BEGIN
  24.  
  25. NUMERICS_EXPORT double acosh(double x)
  26. {
  27.     return ((x <= 1.0) ? 0.0 :
  28.             ((x > 1.0e10) ? 0.69314718055995 + log(x) :
  29.              log(x + sqrt((x-1.0)*(x+1.0)))));
  30. }
  31.  
  32. NUMERICS_EXPORT double asinh(double x)
  33. {
  34.     double ax = fabs(x);
  35.     if(ax > 1.0e10){
  36.         return ((x > 0.0) ? 0.69314718055995 + log(ax) :
  37.                 -.69314718055995 + log(ax));
  38.     } else {
  39.         double y = x*x;
  40.         return ((x == 0.0) ? 0.0 :
  41.                 ((x > 0.0) ? log1x(ax+y/(1.0+sqrt(1.0+y))) :
  42.                  -log1x(ax+y/(1.0+sqrt(1.0+y)))));
  43.     }
  44. }
  45.  
  46. NUMERICS_EXPORT double atanh(double x)
  47. {
  48.     double ax = fabs(x);
  49.     
  50.     if(ax >= 1.0){
  51.         return ((x > 0.0) ? FLT_MAX : -FLT_MAX);
  52.     } else {
  53.         return ((x == 0.0) ? 0.0 :
  54.                 ((x > 0.0) ? .0*log1x(2.0*ax/(1.0-ax)) :
  55.                  -.5*log1x(2.0*ax/(1.0-ax))));
  56.     }
  57. }
  58.  
  59. NUMERICS_EXPORT double log1x(double x)
  60. {
  61.     if(x == 0.0){
  62.         return 0.0;
  63.     } else if(x < -.2928 || x > .4142){
  64.         return log(1.0+x);
  65.     } else {
  66.         double z = x/(x+2.0);
  67.         double y = z*z;
  68.         return (z*(2.0+y*(0.666666666663366+y*(0.400000001206045+y*
  69.             (0.285714091590488+y*(0.22223823332791+y*
  70.             (0.1811136267967+y*0.16948212488)))))));
  71.     }
  72. }
  73.  
  74. NUMERICS_EXPORT double betacf(double a, double b, double x)
  75. {
  76.     int m,m2;
  77.     double aa,c=1.0,del,h,qab=a+b,qam=a-1.0,qap=a+1.0, d=1.0-qab*x/qap;
  78.     
  79.     if(areEqual(d, 0.0)) d=NUMERICS_FLOAT_MIN;
  80.     d=1.0/d;
  81.     h=d;
  82.     for(m=1; m <= NUMERICS_ITMAX; m++){
  83.         m2=2*m;
  84.         aa=m*(b-m)*x/((qam+m2)*(a+m2));
  85.         d=1.0+aa*d;
  86.         if(areEqual(d, 0.0)) d = NUMERICS_FLOAT_MIN;
  87.         c=1.0+aa/c;
  88.         if(areEqual(c, 0.0)) c = NUMERICS_FLOAT_MIN;
  89.         d=1.0/d;
  90.         h*=d*c;
  91.         aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
  92.         d=1.0+aa*d;
  93.         if(areEqual(d, 0.0)) d = NUMERICS_FLOAT_MIN;
  94.         c=1.0+aa/c;
  95.         if(areEqual(c, 0.0)) c = NUMERICS_FLOAT_MIN;
  96.         d=1.0/d;
  97.         del=d*c;
  98.         h*=del;
  99.         if(areEqual(del, 1.0)) break;
  100.     }
  101.     if(m > NUMERICS_ITMAX){
  102.         throw Exception("betacf", "a or b too big, or NUMERICS_ITMAX too small");
  103.     }
  104.     return h;
  105. }
  106.  
  107. void gser(double *gamser, double a, double x, double *gln)
  108. {
  109.     int n;
  110.     double sum, del, ap;
  111.     
  112.     *gln = gammln(a);
  113.     if(x <= 0.0){
  114.         if(x < 0.0){
  115.             throw Exception("gser", "x less than 0");
  116.         }
  117.         *gamser=0.0;
  118.         return;
  119.     } else {
  120.         ap = a;
  121.         del = sum = 1.0/a;
  122.         for(n = 1; n <=NUMERICS_ITMAX; n++){
  123.             ++ap;
  124.             del *= x / ap;
  125.             sum += del;
  126.             if(fabs(del) < fabs(sum)*NUMERICS_MAX_ERROR){
  127.                 *gamser = sum * exp(-x+a*log(x)-(*gln));
  128.                 return;
  129.             }
  130.         }
  131.     }
  132.     throw Exception("gser", "a too large, NUMERICS_ITMAX too small");
  133.     return;
  134. }
  135.  
  136. void gcf(double *gammcf, double a, double x, double *gln)
  137. {
  138.     int i;
  139.     double an, b, c, d, del, h;
  140.     
  141.     *gln = gammln(a);
  142.     b = x + 1.0 - a;
  143.     c = 1.0 / NUMERICS_FLOAT_MIN;
  144.     d = 1.0 / b;
  145.     h = d;
  146.     for(i = 1; i <= NUMERICS_ITMAX; i++){
  147.         an = -i*(i-a);
  148.         b += 2.0;
  149.         d = an*d+b;
  150.         if(areEqual(d, 0.0)) d = NUMERICS_FLOAT_MIN;
  151.         c = b + an / c;
  152.         if(areEqual(c, 0.0)) c = NUMERICS_FLOAT_MIN;
  153.         d = 1.0 / d;
  154.         del = d*c;
  155.         h *= del;
  156.         if(areEqual(del, 1.0)) break;
  157.     }
  158.     if(i > NUMERICS_ITMAX){
  159.         throw Exception("gcf", "a too large, NUMERICS_ITMAX too small");
  160.     }
  161.     *gammcf = exp(-x+a*log(x)-(*gln))*h;
  162. }
  163.  
  164. NUMERICS_EXPORT double factrl(int n)
  165. {
  166.     static int ntop = 4;
  167.     static double a[33] = {1.0, 1.0, 2.0, 6.0, 24.0};
  168.     int j;
  169.     
  170.     if(n < 0){
  171.         throw Exception("factrl", "Negative factorial");
  172.         return 0.0;
  173.     }
  174.     if(n > 32) return exp(gammln(n + 1.0));
  175.     while(ntop < n){
  176.         j = ntop++;
  177.         a[ntop] = a[j] * ntop;
  178.     }
  179.     return a[n];
  180. }
  181.  
  182. NUMERICS_EXPORT double betai(double a, double b, double x)
  183. {
  184.     double bt;
  185.     
  186.     if(x < 0.0 || x > 1.0){
  187.         throw Exception("betai", "Bad x");
  188.         return 0.0;
  189.     }
  190.     if(x == 0.0 || x == 1.0) bt = 0.0;
  191.     else bt = exp(gammln(a + b) - gammln(a) - gammln(b) + a * log(x) + b *
  192.                   log(1.0 - x));
  193.     if(x < (a + 1.0) / (a + b + 2.0)) return bt * betacf(a, b, x) / a;
  194.     else return 1.0 - bt * betacf(b, a, 1.0 - x) / b;
  195. }
  196.  
  197. NUMERICS_EXPORT double gammln(double xx)
  198. {
  199.     double x, y, tmp, ser;
  200.     static double cof[6]={76.18009172947146, -86.50532032941677,
  201.         24.01409824083091, -1.231739572460166, 0.1208650973866179e-2,
  202.         -0.5395239384953e-5};
  203.     int j;
  204.  
  205.     y = x = xx;
  206.     tmp = x + 5.5;
  207.     tmp -= (x + 0.5) * log(tmp);
  208.     ser = 1.000000000190015;
  209.     for(j = 0; j <= 5; j++) ser += cof[j] / ++y;
  210.     return -tmp + log(2.5066282746310005 * ser / x);
  211. }
  212.  
  213. NUMERICS_EXPORT double gammp(double a, double x)
  214. {
  215.     double gamser, gammcf, gln;
  216.     
  217.     if(x < 0.0 || a <= 0.0){
  218.         throw Exception("gammp", "Invalid arguments");
  219.         return 0.0;
  220.     }
  221.     if(x < (a+1.0)){
  222.         gser(&gamser, a, x, &gln);
  223.         return gamser;
  224.     } else {
  225.         gcf(&gammcf, a, x, &gln);
  226.         return 1.0 - gammcf;
  227.     }
  228. }
  229.  
  230. NUMERICS_EXPORT double gammq(double a, double x)
  231. {
  232.     double gamser, gammcf, gln;
  233.     
  234.     if(x < 0.0 || a <= 0.0){
  235.         throw Exception("gammq", "Invalid arguments");
  236.         return 0.0;
  237.     }
  238.     if(x < (a+1.0)){
  239.         gser(&gamser, a, x, &gln);
  240.         return 1.0 - gamser;
  241.     } else {
  242.         gcf(&gammcf, a, x, &gln);
  243.         return gammcf;
  244.     }
  245. }
  246.  
  247. NUMERICS_EXPORT double beta(double z, double w)
  248. {
  249.     return exp(gammln(z) + gammln(w) - gammln(z+w));
  250. }
  251.  
  252. NUMERICS_EXPORT double erff(double x)
  253. {
  254.     return x < 0.0 ? -gammp(.5, x * x) : gammp(.5, x * x);
  255. }
  256.  
  257. NUMERICS_EXPORT double bico(int n, int k)
  258. {
  259.     return floor(0.5 + exp(factln(n) - factln(k) - factln(n - k)));
  260. }
  261.  
  262. NUMERICS_EXPORT double bicoln(int n, int k)
  263. {
  264.     return factln(n) - factln(k) - factln(n - k);
  265. }
  266.  
  267. NUMERICS_EXPORT double factln(int n)
  268. {
  269.     static double a[101];
  270.  
  271.     if(n < 0){
  272.         throw Exception("factrl", "Negative factorial");
  273.         return 0.0;
  274.     }
  275.     if(n <= 1) return 0.0;
  276.     if(n <= 100) return a[n] ? a[n] : (a[n] = gammln(n + 1.0));
  277.     else return gammln(n + 1.0);
  278. }
  279.  
  280. NUMERICS_EXPORT double erffc(double x)
  281. {
  282.     return x < 0.0 ? 1.0 + gammp(0.5, x*x) : gammq(0.5, x*x);
  283. }
  284.  
  285. NUMERICS_EXPORT double erfcc(double x)
  286. {
  287.     double t,z, ans;
  288.  
  289.     z = fabs(x);
  290.     t=1.0/(1.0+.5*z);
  291.     ans=t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*
  292.         (.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+t*
  293.         (1.48851587+t*(-.82215223+t*.17087277)))))))));
  294.     
  295.     return x >= 0.0 ? ans : 2.0-ans;
  296. }
  297.  
  298. NUMERICS_EXPORT double expint(int n, double x)
  299. {
  300.     int i, ii, nm1;
  301.     double a,b,c,d,del,fact,h,psi,ans;
  302.     
  303.     nm1 = n-1;
  304.     if(n<0||x<0.0||(x==0.0 &&(n==0||n==1))) return 0.0; // Error condition
  305.     else {
  306.         if(n==0)ans=exp(-x)/x;
  307.         else {
  308.             if(x==0.0) ans=1.0/nm1;
  309.             else {
  310.                 if(x > 1.0){
  311.                     b=x+n;
  312.                     c=1.0/NUMERICS_FLOAT_MIN;
  313.                     d=1.0/b;
  314.                     h=d;
  315.                     for(i=1;i<=NUMERICS_ITMAX;i++){
  316.                         a = -i*(nm1+1);
  317.                         b+=2.0;
  318.                         d=1.0/(a*d+b);
  319.                         c=b+a/c;
  320.                         del =c*d;
  321.                         h*=del;
  322.                         if(areEqual(del, 1.0)){
  323.                             ans=h*exp(-x);
  324.                             return ans;
  325.                         }
  326.                     }
  327.                     return ans; // Error condition
  328.                 } else {
  329.                     ans =(nm1!=0 ? 1.0/nm1 : -log(x)-NUMERICS_EULER);
  330.                     fact=1.0;
  331.                     for(i=1;i<=NUMERICS_ITMAX;i++){
  332.                         fact*=-x/i;
  333.                         if(i!=nm1) del=-fact/(i-nm1);
  334.                         else {
  335.                             psi = -NUMERICS_EULER;
  336.                             for(ii=1;ii<=nm1;ii++) psi += 1.0/ii;
  337.                             del=fact*(-log(x)+psi);
  338.                         }
  339.                         ans+=del;
  340.                         if(fabs(del) < fabs(ans)*NUMERICS_MAX_ERROR) return ans;
  341.                     }
  342.                     return ans; // Error condition
  343.                 }
  344.             }
  345.         }
  346.     }
  347.     return ans;
  348. }
  349.  
  350. NUMERICS_EXPORT double pythag(double a, double b)
  351. {
  352.     double aa = fabs(a), ab = fabs(b);
  353.     
  354.     if(aa > ab) return (aa * sqrt(1.0 + (ab/aa * ab/aa)));
  355.     else return (ab == 0.0 ? 0.0 : ab * sqrt(1.0 + (aa/ab * aa/ab)));
  356. }
  357.  
  358. NUMERICS_EXPORT double sign(double x)
  359. {
  360.     if(x < 0.0) return -1.0;
  361.     if(x > 0.0) return 1.0;
  362.  
  363.     return 0.0;
  364. }
  365.  
  366. NUM_END
  367.  
  368. //====================================================================
  369. // Revision History
  370. //
  371. // Version 1.0 - 08/28/1998 - New.
  372. // Version 1.1 - 04/10/1999 - Added Numerics namespace.
  373. //                            Added use of Exception class.
  374. //====================================================================
  375.