home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 7 / 07.iso / c / c129 / 1.ddi / SPECFUNC.C < prev    next >
Encoding:
C/C++ Source or Header  |  1989-04-11  |  18.0 KB  |  843 lines
  1.  
  2.  
  3. #include <math.h>
  4. #include <stdio.h>
  5. #include "miscio.h"
  6.  
  7. # define NearZero 1.0E-20
  8. # define Tolerance 1.0E-10
  9.  
  10. void DisplayErrorMessage(int ErrorFlag)
  11. {
  12.  
  13.    switch ( ErrorFlag ) {
  14.       case 0:
  15.          printf("No Errors Found.");
  16.          break;
  17.       case 1:
  18.          printf( "ERROR!!  Argument must be a nonzero number");
  19.          break;
  20.       case 2:
  21.          printf("ERROR!!  Argument must not be a negative");
  22.          break;
  23.       case 3:
  24.          printf( "ERROR!!  Order of polynomial must be greater than zero");
  25.          break;
  26.       case 4:
  27.          printf( "ERROR!!  Argument must be greater than -1");
  28.          break;
  29.    }
  30. }
  31.  
  32. float CalcPoly(float XIn,float CoefVector[],int order)
  33. {
  34.    float _fret = 4;
  35.    int I;
  36.    float YOut;
  37.  
  38.    YOut = CoefVector[order];
  39.    for ( I = order - 1; I >= 0; --I ) {
  40.       YOut = YOut * XIn + CoefVector[I];
  41.    }
  42.    _fret = YOut;
  43.    return(_fret);
  44. }
  45.  
  46. float CalcPower(float x,float a)
  47. {
  48.    float _fret;
  49.  
  50.    _fret = exp(a * log(x));
  51.    return(_fret);
  52. }
  53.  
  54. float LogGamma(float x)
  55. {
  56.    float _fret;
  57.    float coef[16];
  58.    float ser;
  59.    float temp;
  60.    float xx;
  61.    float stp;
  62.    int I;
  63.    coef[0] = 7.618009173E1;
  64.    coef[1] = -8.650532033E1;
  65.    coef[2] = 2.401409822E1;
  66.    coef[3] = -1.231739516;
  67.    coef[4] =  1.20858003E-3;
  68.    coef[5] = -5.36382E-6;
  69.    stp     = 2.50662827565;
  70.    xx = x - 1.0;
  71.    temp = xx + 5.5;
  72.    temp = (xx + 0.5) * log(temp) - temp;
  73.    ser = 1.0;
  74.    for ( I = 0; I <= 5; ++I ) {
  75.       xx = xx + 1.0;
  76.       ser = ser + coef[I] / xx;
  77.    }
  78.    _fret = temp + log(stp * ser);
  79.  
  80.    return(_fret);
  81. }
  82.  
  83. float StandGamma(float y)
  84. {
  85.    float _fret;
  86.    typedef float BigArray[21];
  87.    float temp;
  88.    float f;
  89.    int I;
  90.    float BigArry[21];
  91.     BigArry[0] =   0.5772156649015329;
  92.     BigArry[1] =  -0.6558780715202538;
  93.     BigArry[2] =  -0.0420026350340952;
  94.     BigArry[3] =   0.1665386113822915;
  95.     BigArry[4] =  -0.0421977345555443;
  96.     BigArry[5] =  -0.009621971527877;
  97.     BigArry[6] =   0.0072189432466630;
  98.     BigArry[7] =  -0.0011651675918591;
  99.     BigArry[8] =  -0.0002152416741149;
  100.     BigArry[9] =   0.0001280502823882;
  101.     BigArry[10] = -0.0000201348547807;
  102.     BigArry[11] = -0.0000012504934821;
  103.     BigArry[12] =  0.0000011330272320;
  104.     BigArry[13] = -0.0000002056338417;
  105.     BigArry[14] =  0.0000000061160950;
  106.     BigArry[15] =  0.0000000050020075;
  107.     BigArry[16] = -0.0000000011812746;
  108.     BigArry[17] =  0.0000000001043427;
  109.     BigArry[18] =  0.0000000000077823;
  110.     BigArry[19] = -0.0000000000036968;
  111.     BigArry[20] =  0.00000000000005100;
  112.    f = y;
  113.    temp = f;
  114.    for ( I = 0; I <= 20; ++I ) {
  115.       f = f * y;
  116.       temp = temp + f * BigArry[I];
  117.    }
  118.    _fret = 1.0 / temp;
  119.    return(_fret);
  120. }
  121.  
  122. float Gamma(float x,int *ErrorFlag)
  123. {
  124.    float _fret;
  125.    float y;
  126.    float f;
  127.    if ( (fabs(x) < NearZero) ) {
  128.       (*ErrorFlag) = 1;
  129.    }
  130.    else {
  131.       if ( ((x < 0) && (fabs(x - ceil(0.5+x)) < NearZero)) ) {
  132.          (*ErrorFlag) = 2;
  133.       }
  134.       else {
  135.          if ( fabs(x-1.0) < NearZero ) {
  136.             _fret = 1.0;
  137.          }
  138.          else {
  139.             if ( x < 0.0 ) {
  140.                _fret = -(M_PI / (x * sin(M_PI * x) * Gamma(-x,ErrorFlag)));
  141.             }
  142.             else {
  143.                y = x;
  144.                f = 1.0;
  145.                while ( y > 1.0 ) {
  146.  
  147.                y = y - 1.0;
  148.                   f = f * y;
  149.                }
  150.                if ( fabs(y-1.0) < NearZero ) {
  151.                   _fret = f;
  152.                }
  153.                else {
  154.                   if ( y > (1.0 / 2.0) ) {
  155.                      _fret = f * M_PI / (sin(M_PI * y)* Gamma(1.0 - y,ErrorFlag));
  156.                   }
  157.                   else {
  158.                      _fret = f * StandGamma(y);
  159.                   }
  160.                }
  161.             }
  162.          }
  163.       }
  164.    }
  165.    return(_fret);
  166. }
  167.  
  168. float IncGamma(float a,float x)
  169. {
  170.    float _fret;
  171.    float ir;
  172.    float dif;
  173.    float s1;
  174.    float s2;
  175.    float gln;
  176.    float result;
  177.  
  178.    if ( (a < 0.0) || (x <= 0.0) ) {
  179.       result = 0.0;
  180.    }
  181.    else {
  182.       gln = LogGamma(a);
  183.       ir = a;
  184.       s1 = 1.0 / a;
  185.       dif = s1;
  186.       while ( (fabs(dif) > fabs(s1) * Tolerance) ) {
  187.          ir = ir + 1.0;
  188.          dif = dif * x / ir;
  189.          s1 = s1 + dif;
  190.       }
  191.       result = s1 * exp(-x + a * log(x) - gln);
  192.    }
  193.    _fret = result;
  194.    return(_fret);
  195. }
  196.  
  197. float IncGammaComp(float a,float x)
  198. {
  199.    float _fret;
  200.  
  201.    _fret = 1.0 - IncGamma(a,x);
  202.    return(_fret);
  203. }
  204.  
  205. float Bessel(float order,float x,int *ErrorFlag)
  206. {
  207.    float _fret;
  208.    float j;
  209.    float d;
  210.    float f;
  211.    float sum;
  212.    int c;
  213.  
  214.    if ( order < NearZero ) {
  215.       (*ErrorFlag) = 3;
  216.    }
  217.    else {
  218.       if ( ((x < 0.0) && (fabs(x - ceil(0.5 + x)) < NearZero)) ) {
  219.          (*ErrorFlag) = 2;
  220.       }
  221.       else {
  222.          if ( fabs(x) < NearZero ) {
  223.             if ( fabs(order) < NearZero ) {
  224.                _fret = 1.0;
  225.             }
  226.             else {
  227.                _fret = 0;
  228.             }
  229.          }
  230.          else {
  231.             if ( x <= 15.0 ) {
  232.                f = 1.0 / Gamma(order + 1,ErrorFlag);
  233.                j = -((x*x) / 4.0);
  234.                d = order + 1.0;
  235.                sum = f;
  236.                c = 1;
  237.                do {
  238.                   f = f * j / (c * d);
  239.                   d = d + 1.0;
  240.                   c = c + 1;
  241.                   sum = sum + f;
  242.                } while ( ! ((fabs(f) < Tolerance)) );
  243.                if ( x > 0.0 ) {
  244.                   _fret = sum * exp(order * log(x / 2.0));
  245.                }
  246.                else {
  247.                   _fret = (1.0 - 2.0 * (fmod(floor(order),2.0))
  248.                            * sum * exp(order * log(fabs(x / 2.0) )));
  249.                }
  250.             }
  251.             else {
  252.                _fret = sqrt(2.0 / (M_PI * x)) *
  253.                        cos(x - M_PI / 4.0 - order * M_PI * 2.0);
  254.             }
  255.          }
  256.       }
  257.    }
  258.    return(_fret);
  259. }
  260.  
  261. float Tan(float x)
  262. {
  263.    float _fret;
  264.  
  265.    _fret = sin(x) / cos(x);
  266.    return(_fret);
  267. }
  268.  
  269. float Cosh(float z)
  270. {
  271.    float _fret;
  272.  
  273.    _fret = (exp(z) + exp(-z)) / 2.0;
  274.    return(_fret);
  275. }
  276.  
  277. float Sinh(float z)
  278. {
  279.    float _fret;
  280.  
  281.    _fret = (exp(z) - exp(-z)) / 2.0;
  282.    return(_fret);
  283. }
  284.  
  285. float Sech(float z)
  286. {
  287.    float _fret;
  288.  
  289.    _fret = 1.0 / Cosh(z);
  290.    return(_fret);
  291. }
  292.  
  293. float ArcTanh(float x)
  294. {
  295.    float _fret;
  296.  
  297.    _fret = 0.5 * log((1.0 + x) / (1.0 - x));
  298.    return(_fret);
  299. }
  300.  
  301. float ErrFuncIter(float x)
  302. {
  303.    float _fret;
  304.    float a;
  305.    float s1;
  306.    float s2;
  307.    float t;
  308.    float result;
  309.    int I;
  310.    int j;
  311.  
  312.    a = x*x;
  313.    s1 = x;
  314.    t = x;
  315.    I = 1;
  316.    while ( t > (Tolerance * s1) ) {
  317.       s2 = s1;
  318.       t = 2.0 * t * (a / (1.0 + 2.0 * I));
  319.       s1 = t + s2;
  320.       I = I + 1;
  321.    }
  322.    result = 2.0 * s1 * (exp(-a)) / sqrt(M_PI);
  323.    _fret = result;
  324.    return(_fret);
  325. }
  326.  
  327. float ErrFunc(float x,int *ErrorFlag)
  328. {
  329.    float _fret;
  330.    float PolyCoef[16];
  331.    float result;
  332.    float temp;
  333.    float PosNeg;
  334.    float p;
  335.    float ANS;
  336.    int I;
  337.  
  338.    if ( fabs(x) < NearZero ) {
  339.       (*ErrorFlag) = 1;
  340.    }
  341.    else {
  342.       if ( x > 4.0 ) {
  343.          result = 1.0;
  344.       }
  345.       else {
  346.          if ( x > 1.5 ) {
  347.             result = ErrFuncIter(x);
  348.          }
  349.          else {
  350.             p = x * x;
  351.             PolyCoef[0] = 0.0;
  352.             PolyCoef[0] = 0.0;
  353.             PolyCoef[1] = 0.6666667;
  354.             PolyCoef[2] = 0.2666667;
  355.             PolyCoef[3] = 0.07619048;
  356.             PolyCoef[4] = 0.01693122;
  357.             PolyCoef[5] = 3.078403E-3;
  358.             PolyCoef[6] = 4.736005E-4;
  359.             PolyCoef[7] = 6.314673E-5;
  360.             PolyCoef[8] = 7.429027E-6;
  361.             PolyCoef[9] = 7.820028E-7;
  362.             PolyCoef[10] = 7.447646E-8;
  363.             PolyCoef[11] = 6.476214E-9;
  364.         temp = x + x * CalcPoly(p,PolyCoef,11);
  365.             result = 2.0 * temp * exp(-p) / sqrt(M_PI);
  366.          }
  367.       }
  368.    }
  369.    _fret = result;
  370.    return(_fret);
  371. }
  372.  
  373. float ErrFuncComp(float x,int *ErrorFlag)
  374. {
  375.    float _fret;
  376.  
  377.    _fret = 1.0 - ErrFunc(x,ErrorFlag);
  378.    return(_fret);
  379. }
  380.  
  381. float ErrFuncR(float x,float y, int *ErrorFlag)
  382. {
  383.    float _fret;
  384.    float jy;
  385.    float xy2;
  386.    float x2;
  387.    float sum;
  388.    float part;
  389.    float f;
  390.    int j;
  391.  
  392.    x2 = 2.0 * x;
  393.    xy2 = x2 * y;
  394.    if ( fabs(x) > 0.0e-8 ) {
  395.       part = exp(-(x*x)) / x * (1.0 / (2.0 * M_PI)) * (1.0 - cos(xy2));
  396.    }
  397.    else {
  398.       part = 0.0;
  399.    }
  400.    sum = ErrFunc(x,ErrorFlag    );
  401.    sum = sum + part;
  402.    part = 0.0;
  403.    for ( j = 1; j <= 10; ++j ) {
  404.       jy = j * y;
  405.       f = x2 - x2 * cos(xy2) * Cosh(jy) + j * Sinh(jy) * sin(xy2);
  406.       part = part + exp(-((j*j) / 4.0)) / ((j*j) + 4.0 * (x*x)) * f;
  407.    }
  408.    sum = sum + part * 2.0 / M_PI * exp(-(x*x));
  409.    _fret = sum;
  410.    return(_fret);
  411. }
  412.  
  413. float ErrFuncI(float x,float y)
  414. {
  415.    float _fret;
  416.    float jy;
  417.    float x2;
  418.    float xy2;
  419.    float sum;
  420.    float part;
  421.    float g;
  422.    int j;
  423.  
  424.    x2 = 2.0 * x;
  425.    xy2 = x2 * y;
  426.    sum = 0;
  427.    if ( fabs(x) < 0.0000001 ) {
  428.       part = y / M_PI;
  429.    }
  430.    else {
  431.       part = exp(-(x*x)) / (M_PI * x2) * sin(xy2);
  432.    }
  433.    sum = sum + part;
  434.    part = 0.0;
  435.    for ( j = 1; j <= 20; ++j ) {
  436.       jy = j * y;
  437.       g = x2 * Cosh(jy) * sin(xy2) + j * Sinh(jy) * cos(xy2);
  438.       part = part * exp(-((j*j) / 4.0)) / ((j*j) + (x2*x2)) * g;
  439.    }
  440.    sum = sum + part * exp(-(x*x)) * 2.0 / M_PI;
  441.    _fret = sum;
  442.    return(_fret);
  443. }
  444.  
  445. float GaussHyper(float a,float b,float c,float z)
  446. {
  447.    float _fret;
  448.    float g;
  449.    float sum;
  450.    int j;
  451.  
  452.    g = 1.0;
  453.    j = 1;
  454.    sum = 0.0;
  455.    do {
  456.       sum = sum + g;
  457.       g = g * a * b * z / (j * c);
  458.       j = j + 1;
  459.       a = a + 1.0;
  460.       b = b + 1.0;
  461.       c = c + 1.0;
  462.    } while ( ! ((j > 80) || (fabs(g) < Tolerance)) );
  463.    _fret = sum;
  464.    return(_fret);
  465. }
  466.  
  467. float KumrConf(float a,float b,float z)
  468. {
  469.    float _fret;
  470.    float g;
  471.    float sum;
  472.    int j;
  473.  
  474.    sum = 0.0;
  475.    g = 1.0;
  476.    j = 1;
  477.    do {
  478.       sum = sum + g;
  479.       g = g * a * z / (b * j);
  480.       j = j + 1;
  481.       a = a + 1.0;
  482.       b = b + 1.0;
  483.    } while ( ! ((j > 80) || (fabs(g) < Tolerance)) );
  484.    _fret = sum;
  485.    return(_fret);
  486. }
  487.  
  488. float AssocConf(float a,float b,float z, int *ErrorFlag)
  489. {
  490.    float _fret;
  491.    float x;
  492.    float PiCalc;
  493.    float Gam1;
  494.    float Kum1;
  495.    float exp1;
  496.    float kum2;
  497.    float gam2;
  498.    float Calc1;
  499.    float calc2;
  500.    float calc3;
  501.    float subtotal;
  502.  
  503.    if ( fabs(z) < 1E-10 ) {
  504.       _fret = 0.0;
  505.    }
  506.    else {
  507.       if ( (fabs(b - floor(b)) < 0.0001) || (fabs(b - ceil(0.5+b)) < 0.0001) ) {
  508.          _fret = (AssocConf(a,b - 0.0005,z,ErrorFlag) + AssocConf(a,b + 0.0005,z,ErrorFlag)) / 2.0;
  509.       }
  510.       else {
  511.          PiCalc = M_PI / sin(M_PI * b);
  512.          Kum1 = KumrConf(a,b,z);
  513.          Gam1 = Gamma(1.0 + a - b,ErrorFlag) * Gamma(b,ErrorFlag);
  514.          Gam1 = Gamma(1.0 + a - b,ErrorFlag) * Gamma(b,ErrorFlag);
  515.      exp1 = CalcPower(z,1.0 - b);
  516.          kum2 = KumrConf(1.0 + a - b,2.0 - b,z);
  517.          gam2 = Gamma(a,ErrorFlag) * Gamma(2.0 - b,ErrorFlag);
  518.          Calc1 = Kum1 / Gam1;
  519.          calc2 = exp1 * kum2;
  520.          calc3 = calc2 / gam2;
  521.          subtotal = Calc1 - calc3;
  522.          _fret = PiCalc * subtotal;
  523.          scanf("");
  524.       }
  525.    }
  526.    return(_fret);
  527. }
  528.  
  529.  
  530. float Hermite(int n,float x,int *ErrorFlag)
  531. {
  532.    float _fret;
  533.  
  534.    if ( n < 0 ) {
  535.       (*ErrorFlag) = 3;
  536.    }
  537.    else {
  538.       switch ( n ) {
  539.          case 0:
  540.             _fret = 1.0;
  541.             break;
  542.          case 1:
  543.             _fret = 2.0 * x;
  544.             break;
  545.          case 2:
  546.             _fret = 4.0 * (x*x) - 2.0;
  547.             break;
  548.          case 3:
  549.         _fret = 8.0 * CalcPower(x,3.0) - (12.0 * x);
  550.             break;
  551.          case 4:
  552.         _fret = 16.0 * CalcPower(x,4.0) - 48.0 * (x*x) + 12.0;
  553.             break;
  554.          case 5:
  555.         _fret = 32.0 * CalcPower(x,5.0) - 160.0 * CalcPower(x,3.0) + 120.0 * x;
  556.             break;
  557.          default:
  558.             _fret = x * exp(n * log(2)) * AssocConf(0.5 - n / 2.0,3.0 / 2.0,(x*x),ErrorFlag);
  559.             break;
  560.       }
  561.    }
  562.    return(_fret);
  563. }
  564.  
  565. float Legend(int n,float x,int *ErrorFlag)
  566. {
  567.    float _fret = 5.0;
  568.  
  569.    if ( n < 0 ) {
  570.       (*ErrorFlag) = 3;
  571.    }
  572.    else {
  573.       switch ( n ) {
  574.          case 0:
  575.             _fret = 1.0;
  576.             break;
  577.          case 1:
  578.             _fret = x;
  579.             break;
  580.          case 2:
  581.             _fret = 0.5 * (3.0 * x * x - 1.0);
  582.             break;
  583.          case 3:
  584.         _fret = 0.5 * (5.0 * CalcPower(x,3.0) - (3.0 * x));
  585.             break;
  586.          case 4:
  587.         _fret = 1.0 / 8.0 * (35.0 * CalcPower(x,4.0) - 30.0 * 4.0 + 3.0);
  588.             break;
  589.          case 5:
  590.         _fret = 1.0 / 8.0 * (63.0 * CalcPower(x,5.0) - 70.0 * CalcPower(x,3.0) + 15.0 * x);
  591.             break;
  592.          case 6:
  593.         _fret = 1.0 / 16.0 * (231.0 * CalcPower(x,6.0) - 315.0 * CalcPower(x,4.0) + 105.0 * 4.0 - 5.0);
  594.             break;
  595.          default:
  596.             _fret = GaussHyper(-n,n + 1.0,1.0,(1.0 - x) / 2.0);
  597.             break;
  598.       }
  599.    }
  600.    return(_fret);
  601. }
  602.  
  603. float Laguerre(int n,float a,float x,int *ErrorFlag)
  604. {
  605.    float _fret = 6.0;
  606.    float num;
  607.  
  608.    if ( n < 0 ) {
  609.       (*ErrorFlag) = 3;
  610.    }
  611.    else {
  612.       if ( a <= -1 ) {
  613.          (*ErrorFlag) = 4;
  614.       }
  615.       else {
  616.          if ( n % 2 == 0 ) {
  617.             num = 1.0;
  618.          }
  619.          else {
  620.             num = -1.0;
  621.          }
  622.          switch ( n ) {
  623.             case 0:
  624.                _fret = 1.0;
  625.                break;
  626.             case 1:
  627.                _fret = -x + 1.0;
  628.                break;
  629.             case 2:
  630.                _fret = 0.5 * ((x*x) - 4.0 * x + 2.0);
  631.                break;
  632.             case 3:
  633.            _fret = 1.0 / 6.0 * (-CalcPower(x,3.0) + 9.0 * (x*x) -
  634.                        18.0 * x + 6.0);
  635.                break;
  636.             case 4:
  637.            _fret = 1.0 / 24.0 * (CalcPower(x,4.0) - 16 *
  638.                CalcPower(x,3.0) + 72.0 * (x*x) - 96.0 * x + 24.0);
  639.                break;
  640.             default:
  641.                _fret = num * AssocConf(-n,a + 1.0,x,ErrorFlag) / Gamma(n + 1.0,ErrorFlag);
  642.                break;
  643.          }
  644.       }
  645.    }
  646.    return(_fret);
  647. }
  648.  
  649. float Combin(float n,int m)
  650. {
  651.    float _fret = 4.0;
  652.  
  653.    if ( m == 1 ) {
  654.       _fret = n;
  655.    }
  656.    else {
  657.       _fret = Combin(n - 1.0,m - 1) * n / (1.0*m);
  658.    }
  659.    return(_fret);
  660. }
  661.  
  662. float Jacobi(int n,float a,float b,float x,int *ErrorFlag)
  663. {
  664.    float _fret;
  665.  
  666.    if ( n < 0 ) {
  667.       (*ErrorFlag) = 3;
  668.    }
  669.    else {
  670.       if ( ((a <= -1.0) || (b <= -1.0)) ) {
  671.          (*ErrorFlag) = 4;
  672.       }
  673.       else {
  674.          _fret = Combin(n + a,n) *
  675.                  GaussHyper(-n,n + a + b + 1.0,a + 1.0,(1.0 - x) / 2.0);
  676.       }
  677.    }
  678.    return(_fret);
  679. }
  680.  
  681. float Tcheb(int n,float x,int *ErrorFlag)
  682. {
  683.    float _fret;
  684.  
  685.    if ( n < 0 ) {
  686.       (*ErrorFlag) = 3;
  687.    }
  688.    else {
  689.       switch ( n ) {
  690.          case 0:
  691.             _fret = 1.0;
  692.             break;
  693.          case 1:
  694.             _fret = x;
  695.             break;
  696.          case 2:
  697.             _fret = 2.0 * (x*x) - 1.0;
  698.             break;
  699.          case 3:
  700.         _fret = 4.0 * CalcPower(x,3.0) - 3.0 * x;
  701.             break;
  702.          case 4:
  703.         _fret = 8.0 * CalcPower(x,4.0) - 8.0 * (x*x) + 1.0;
  704.             break;
  705.          case 5:
  706.         _fret = 16.0 * CalcPower(x,5.0) - 20.0 * CalcPower(x,3.0) + 5.0 * x;
  707.             break;
  708.          case 6:
  709.         _fret = 32.0 * CalcPower(x,6.0) - 48.0 * CalcPower(x,4.0) + 18.0 * (x*x) - 1.0;
  710.             break;
  711.          default:
  712.             _fret = GaussHyper(-n,n,0.5,(1.0 - x) / 2.0);
  713.             break;
  714.       }
  715.    }
  716.    return(_fret);
  717. }
  718.  
  719. float ModBesselI(float n,float x, int *ErrorFlag)
  720. {
  721.    float _fret = 6.0;
  722.  
  723.    if ( fabs(x) < NearZero ) {
  724.       if ( fabs(n) < NearZero ) {
  725.          _fret = 1.0;
  726.       }
  727.       else {
  728.          _fret = 0.0;
  729.       }
  730.    }
  731.    else {
  732.       _fret = exp(-x) * exp(n * log(x / 2.0)) *
  733.               KumrConf(n + 0.5,2.0 * n + 1.0,2.0 * x) / Gamma(n + 1.0,ErrorFlag);
  734.    }
  735.    return(_fret);
  736. }
  737.  
  738. float ModBesselK(float n,float x, int *ErrorFlag)
  739. {
  740.    float _fret;
  741.  
  742.    if ( fabs(x) < NearZero ) {
  743.       if ( fabs(n) < NearZero ) {
  744.          _fret = 1.0;
  745.       }
  746.       else {
  747.          _fret = 0.0;
  748.       }
  749.    }
  750.    else {
  751.       _fret = exp(n * log(2.0 * x)) * exp(-x) * sqrt(M_PI) *
  752.               AssocConf(n + 0.5,2 * n + 1.0,2.0 * x, ErrorFlag);
  753.    }
  754.    return(_fret);
  755. }
  756.  
  757. float Beta(float z,float w,int *ErrorFlag)
  758. {
  759.    float _fret;
  760.    float b;
  761.  
  762.    b = LogGamma(z) + LogGamma(w) - LogGamma(z + w);
  763.    _fret = exp(b);
  764.    (*ErrorFlag) = 0;
  765.    return(_fret);
  766. }
  767.  
  768. float IBeta(float a,float b,float x)
  769. {
  770.    float _fret;
  771.    float m;
  772.    float am;
  773.    float bp;
  774.    float bm;
  775.    float az;
  776.    float qab;
  777.    float qap;
  778.    float qam;
  779.    float bz;
  780.    float em;
  781.    float tem;
  782.    float ap;
  783.    float d;
  784.    float app;
  785.    float bpp;
  786.    float aold;
  787.  
  788.    am = 1.0;
  789.    bm = 1.0;
  790.    az = 1.0;
  791.    qab = a + b;
  792.    qap = a + 1.0;
  793.    qam = a - 1.0;
  794.    bz = 1.0 - qab * x / qap;
  795.    m = 1.0;
  796.    do {
  797.       em = m;
  798.       tem = 2.0 * em;
  799.       d = em * (b - m) * x / ((qam + tem) * (a + tem));
  800.       ap = az + d * am;
  801.       bp = bz + d * bm;
  802.       d = -(a + em) * (qab + em) * x / ((a + tem) * (qap + tem));
  803.       app = ap + d * az;
  804.       bpp = bp + d * bz;
  805.       aold = az;
  806.       am = ap / bpp;
  807.       bm = bp / bpp;
  808.       az = app / bpp;
  809.       bz = 1.0;
  810.       m = m + 1.0;
  811.    } while ( ! ((fabs(az - aold) < (Tolerance * fabs(az)))) );
  812.    _fret = az;
  813.    return(_fret);
  814. }
  815.  
  816. float IncBeta(float a,float b,float x)
  817. {
  818.    float _fret;
  819.    float result;
  820.    float bt;
  821.  
  822.    if ( (x < 0.0) || (x > 1.0) ) {
  823.       printf("%s\n", "Improper Value passed IncBeta");
  824.    }
  825.    else {
  826.       if ( (fabs(x) < NearZero) || ((1.0 - x) < NearZero) ) {
  827.          bt = 0.0;
  828.       }
  829.       else {
  830.          bt = LogGamma(a + b) + a * log(x) + b * log(1.0 - x) - LogGamma(a) - LogGamma(b);
  831.       }
  832.       bt = exp(bt);
  833.       if ( x < ((a + 1.0) / (a + b + 2.0)) ) {
  834.          result = bt * (IBeta(a,b,x) / a);
  835.       }
  836.       else {
  837.          result = 1 - bt * IBeta(b,a,1.0 - x) / b;
  838.       }
  839.       _fret = result;
  840.    }
  841.    return(_fret);
  842. }
  843.