home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 21 / IOPROG_21.ISO / SOFT / LIBMAT.ZIP / MATSPEC.CPP < prev    next >
Encoding:
C/C++ Source or Header  |  1992-10-16  |  17.0 KB  |  651 lines
  1. /**************************************************/
  2. /*    matspec.c source for matSpecFunc class      */
  3. /**************************************************/
  4.  
  5.  
  6. /**************************************************/
  7. /*            MatClass Source File                */
  8. /*       Copyright of C. R. Birchenhall           */
  9. /*       University of Manchester, UK.            */
  10. /*   MatClass is freeware. This file should be    */
  11. /* made freely available to users of any software */
  12. /* whose creation is wholly or partly dependent   */
  13. /*                on this file.                   */
  14. /**************************************************/
  15.  
  16. #include "matspec.hpp"
  17.  
  18. #ifdef __cplusplus
  19.  
  20. const REAL matSpecFunc::zero  =  0.0 ;
  21. const REAL matSpecFunc::half  =  0.5 ;
  22. const REAL matSpecFunc::one   =  1.0 ;
  23. const REAL matSpecFunc::two   =  2.0 ;
  24. const REAL matSpecFunc::three =  3.0 ;
  25. const REAL matSpecFunc::four  =  4.0 ;
  26. const REAL matSpecFunc::five  =  5.0 ;
  27. const REAL matSpecFunc::six   =  6.0 ;
  28. const REAL matSpecFunc::seven =  7.0 ;
  29. const REAL matSpecFunc::eight =  8.0 ;
  30. const REAL matSpecFunc::nine  =  9.0 ;
  31. const REAL matSpecFunc::ten   = 10.0 ;
  32. const REAL matSpecFunc::root2 =  1.414213562 ;
  33. const REAL matSpecFunc::pi    =  3.141592654 ;
  34.  
  35. matSpecFunc::matSpecFunc( void ) : matObject()
  36. {
  37.    error = NOERROR ;
  38.    maxIter = 100 ;
  39.    eps = 1.0E-7 ;
  40. } // matSpecFunc
  41.  
  42. #else
  43.  
  44. static void matSpecInitial( void )
  45. {
  46.  if ( matSpecFunc::one == 1.0 )
  47.    return ;
  48.  matSpecFunc::zero  =  0.0 ;
  49.  matSpecFunc::half  =  0.5 ;
  50.  matSpecFunc::one   =  1.0 ;
  51.  matSpecFunc::two   =  2.0 ;
  52.  matSpecFunc::three =  3.0 ;
  53.  matSpecFunc::four  =  4.0 ;
  54.  matSpecFunc::five  =  5.0 ;
  55.  matSpecFunc::six   =  6.0 ;
  56.  matSpecFunc::seven =  7.0 ;
  57.  matSpecFunc::eight =  8.0 ;
  58.  matSpecFunc::nine  =  9.0 ;
  59.  matSpecFunc::ten   = 10.0 ;
  60.  matSpecFunc::root2 =  1.414213562 ;
  61.  matSpecFunc::pi    =  3.141592654 ;
  62. } // matSpecialInitial
  63.  
  64.  
  65. matSpecFunc::matSpecFunc( void ) : ()
  66. {
  67.    error = 0 ;
  68.    maxIter = 100 ;
  69.    eps = 1.0E-7 ;
  70.    matSpecInitial() ;
  71. } // matSpecFunc
  72.  
  73. #endif
  74.  
  75. matSpecFunc::matSpecFunc( const matSpecFunc& fn )
  76. {
  77.    error = fn.error ;
  78.    maxIter = fn.maxIter ;
  79.    eps = fn.eps ;
  80. } // matSpecFunc
  81.  
  82. matSpecFunc::~matSpecFunc( void ) {}
  83.  
  84. outFile& matSpecFunc::info( outFile& f ) M_CONST
  85. {
  86.     errorExit( "matSpec::info", NIMPL ) ;
  87.     return f ;
  88. } // matSpecFunc info
  89.  
  90. outFile& matSpecFunc::specInfo( outFile& f, const char* fName )
  91.          M_CONST
  92. {
  93.    if ( matListCtrl() > 2 )
  94.       return f ;
  95.    objectInfo( f ) ;
  96.    putName( "", f ) ;
  97.    putName( fName, f ) ;
  98.    f.newLine() ;
  99.    return f ;
  100. } // matSpecFunc specInfo
  101.  
  102. REAL matSpecFunc::value( REAL r )
  103. {
  104.    errorExit( "value", NIMPL ) ;
  105.    return r ;
  106. } // matSpecFunc value
  107.  
  108. REAL matSpecFunc::inv( REAL r )
  109. {
  110.    errorExit( "inv", NIMPL ) ;
  111.    return r ;
  112. } // matSpecFunc inv
  113.  
  114. REAL matSpecFunc::operator() ( REAL r )
  115. {
  116.    REAL v = value(r) ;
  117.    if ( error )
  118.       errorExit( "op()", error ) ;
  119.    return v ;
  120. } // matSpecFunc ()
  121.  
  122. matrix& matSpecFunc::transform( matrix& z, const matrix& x )
  123. {
  124.    static char *mName = "transform" ;
  125.    matFunc func( mName ) ; debugInfo( func ) ;
  126.    INDEX nr = x.nRows(), nc = x.nCols(), i, j ;
  127.    z.reset( nr, nc ) ;
  128.    error = NOERROR ;
  129.    for ( j = 1 ; !error && j <= nc ; j++ )
  130.       for ( i = 1 ; !error && i <= nr ; i++ ) {
  131.          z(i,j) = value( x(i,j) ) ;
  132.          if ( error ) {
  133.             x.errorij( i, j ) ;
  134.             errorExit( mName, error ) ;
  135.          } // if
  136.       } // for i
  137.    return z ;
  138. } // matSpecfunc transform
  139.  
  140. matrix& matSpecFunc::invTrans( matrix& z, const matrix& x )
  141. {
  142.    static char *mName = "invTrans" ;
  143.    matFunc func( mName ) ; debugInfo( func ) ;
  144.    INDEX nr = x.nRows(), nc = x.nCols(), i, j ;
  145.    z.reset( nr, nc ) ;
  146.    error = NOERROR ;
  147.    for ( j = 1 ; !error && j <= nc ; j++ )
  148.       for ( i = 1 ; !error && i <= nr ; i++ ) {
  149.          z(i,j) = inv( x(i,j) ) ;
  150.          if ( error ) {
  151.             x.errorij( i, j ) ;
  152.             errorExit( mName, error ) ;
  153.          } // if
  154.       } // for i
  155.    return z ;
  156. } // matSpecfunc invTrans
  157.  
  158. matrix matSpecFunc::operator () ( const matrix& x )
  159. {
  160.    static char *mName = "matSpecFunc" ;
  161.    matFunc func( mName ) ; debugInfo( func ) ;
  162.    INDEX nr = x.nRows(), nc = x.nCols() ;
  163.    matrix z( nr, nc ) ;
  164.    transform( z, x ) ;
  165.    return z ;
  166. } // matSpecFunc ()
  167.  
  168. matrix matSpecFunc::inv( const matrix& x )
  169. {
  170.    static char *mName = "matSpecFunc" ;
  171.    matFunc func( mName ) ; debugInfo( func ) ;
  172.    INDEX nr = x.nRows(), nc = x.nCols() ;
  173.    matrix z( nr, nc ) ;
  174.    invTrans( z, x ) ;
  175.    return z ;
  176. } // matSpecFunc ()
  177.  
  178. INDEX matSpecFunc::setIter( INDEX newMax )
  179. {
  180.    INDEX oldIter = maxIter ;
  181.    if ( newMax > 0 )
  182.       maxIter = newMax ;
  183.    return oldIter ;
  184. } // incomBetaFunc setIter
  185.  
  186. REAL matSpecFunc::setEps( REAL newEps )
  187. {
  188.    REAL oldEps = eps ;
  189.    if ( newEps >= 0.0 )
  190.      eps = newEps ;
  191.    return oldEps ;
  192. } // incomBeatFunc
  193.  
  194. /************************************************************/
  195. /*                        logGamma class                    */
  196. /************************************************************/
  197.  
  198.  
  199. #ifdef __cplusplus
  200.  
  201. logGammaFunc::logGammaFunc( void ) : matSpecFunc()
  202. {} // logGammaFunc
  203.  
  204. #else
  205.  
  206. logGammaFunc::logGammaFunc( void ) : ()
  207. {} // logGammaFunc
  208.  
  209. #endif
  210.  
  211. logGammaFunc::~logGammaFunc( void ) {}
  212.  
  213. outFile& logGammaFunc::info( outFile& f ) M_CONST
  214. {
  215.    return specInfo( f, "logGam" ) ;
  216. } // info
  217.  
  218. REAL logGammaFunc::value( REAL r )
  219. {
  220.    static REAL coeff[6] = { 76.18009173, -86.50532033,
  221.                               24.01409822, -1.231739516,
  222.                               0.120858003e-2, -0.536382e-5 } ;
  223.    static char *mName = "logGammaFunc" ;
  224.    matFunc func( mName ) ; debugInfo( func ) ;
  225.    REAL x, tmp, ser ;
  226.    INDEX i ;
  227.    error = NOERROR ;
  228.    if ( r < zero ) {
  229.       error = BDARG ;
  230.       return zero ;
  231.    } else if ( r == one )
  232.       return zero ;
  233.    else if ( r < one )
  234.       x = one - r ;
  235.    else // r > one
  236.       x = r - one ;
  237.    tmp = x + 5.5 ;
  238.    tmp -= ( x + half ) * log( tmp ) ;
  239.    ser = one ;
  240.    for ( i = 0 ; i <= 5 ; i++ ) {
  241.       x += one ;
  242.       ser += coeff[i] / x ;
  243.    } // for
  244.    ser = -tmp + log( 2.50662827465 * ser ) ;
  245.    if ( r > one )
  246.       return ser ;
  247.    else { // r < one
  248.       x = one - r ;
  249.       return log( ( pi * x ) / sin( pi * x ) ) - ser ;
  250.    } // else
  251. } // logGammaFunc value
  252.  
  253. logGammaFunc logGamma ;
  254.  
  255. /***********************************************************/
  256. /*                   logBeta Function                      */
  257. /***********************************************************/
  258.  
  259. REAL logBeta( REAL z, REAL w )
  260. {
  261.    static char *mName = "logBeta" ;
  262.    matFunc func( mName ) ;
  263.    return logGamma(z) + logGamma(w) - logGamma(z+w) ;
  264. } // logBeta
  265.  
  266. /************************************************************/
  267. /*                        incBeta class                    */
  268. /************************************************************/
  269.  
  270.  
  271. #ifdef __cplusplus
  272. incBetaFunc::incBetaFunc( void ) : matSpecFunc()
  273. #else
  274. incBetaFunc::incBetaFunc( void ) : ()
  275. #endif
  276. {
  277.    a = b = - 1.0 ;
  278.    c = 0 ;
  279. } // incBetaFunc
  280.  
  281. #ifdef __cplusplus
  282. incBetaFunc::incBetaFunc( incBetaFunc& fn )
  283.               : matSpecFunc(fn)
  284. #else
  285. incBetaFunc::incBetaFunc( incBetaFunc& fn )
  286.               : ( fn )
  287. #endif
  288. {
  289.    a = fn.a ;
  290.    b = fn.b ;
  291.    c = fn.c ;
  292. } // incBetaFunc
  293.  
  294. void incBetaFunc::operator = ( incBetaFunc& fn )
  295. {
  296.    matSpecFunc::operator = ( fn ) ;
  297.    a = fn.a ;
  298.    b = fn.b ;
  299.    c = fn.c ;
  300. } // incBetaFunc
  301.  
  302. incBetaFunc::~incBetaFunc( void ) {}
  303.  
  304. outFile& incBetaFunc::info( outFile& f ) M_CONST
  305. {
  306.    return specInfo( f, "inBeta" ) ;
  307. } // info
  308.  
  309. REAL incBetaFunc::fraction( REAL x, REAL a, REAL b, REAL c )
  310. /************************************************************
  311.    See Press et al NRC p.180
  312.  
  313.    Assumes c = logBeta( a, b )
  314.  
  315.    Calculate incomplete beta function from continued fraction
  316.  
  317.    incomplete_beta(a,b,x) = x^{a}.(1-x)^{b}.fraction /
  318.                                    a * beta(a,b),
  319.    where
  320.  
  321.                1    d1   d2   d3   d4
  322.    fraction = ---  ---- ---- ---- ---- ....
  323.                1+   1+   1+   1+   1+
  324.  
  325.    Use recurrence  fraction(n) = A(n) / B(n)
  326.  
  327.    where
  328.          A(n) = ( s(n) * A(n-1) + r(n) * A(n-2) ) * factor
  329.          B(n) = ( s(n) * B(n-1) + r(n) * B(n-2) ) * factor
  330.    and
  331.          A(-1) = 1, B(-1) = 0, A(0) = s(0), B(0) = 1.
  332.  
  333.    Here s(0) = 0 and for n >= 1 s(n) = 1, while r(1) = 1
  334.    and subsequently, i >= 2
  335.  
  336.          r(i) =  m(b-m)x / (a+i-1)(a+i)       if i = 2m,
  337.          r(i) = -(a+m)(a+b+m)x / (a+i-1)(a+i) if i = 2m+1,
  338.  
  339.    and factor is some scaling factor to avoid overflow.
  340.  
  341.    Hence A(0) = 0 , B(0) = 1,
  342.          r(1) = -(a+b).x / (a+1)
  343.          A(1) = A(0) + r(1).A(-1) = r(1) = 1
  344.          B(1) = B(0) + r(1).B(-1) = 1
  345. ****************************************************************/
  346. {
  347.    REAL old_bta = zero , factor = one;
  348.    REAL A0 = zero, A1 = one, B0 = one, B1 = one ;
  349.    REAL bta = one, am = a, ai = a ;
  350.    REAL m = zero, r ;
  351.    char *mName = "fraction" ;
  352.    matFunc func( mName ) ; debugInfo( func ) ;
  353.  
  354.    error = NOERROR ;
  355.    if ( x < zero || x > one ) {
  356.       error = BDARG ;
  357.       return zero ;
  358.    } else if ( x == zero )
  359.       return zero ;
  360.    else if ( x == one )
  361.       return one ;
  362.    else do {
  363.       // start with i = 1, m = 0, subsequent loops i odd
  364.       ai += one ;
  365.       r = - am * ( am + b ) * x / ( ( ai - one ) * ai ) ;
  366.       // two steps of recurrence replace A's & B's
  367.       A0 = ( A1 + r * A0 ) * factor ; // i odd
  368.       B0 = ( B1 + r * B0 ) * factor ;
  369.       // start here with i = 2, m = 1, subsequently i even
  370.       am += one ;
  371.       m  += one ;
  372.       ai += one ;
  373.       r = m * ( b - m ) * x * factor / ( ( ai - one ) * ai ) ;
  374.       A1 = A0 + r * A1 ; // i even, A0 & B0 already scaled
  375.       B1 = B0 + r * B1 ;
  376.       old_bta = bta ;
  377.       factor = one / B1 ;
  378.       bta = A1 * factor ;
  379.    } while ( m < maxIter && fabs( bta - old_bta )
  380.                                > eps * fabs( bta ) ) ;
  381.    if ( m >= maxIter ) {
  382.       error = NCONV ;
  383.       return zero ;
  384.    } // if
  385.    return bta * exp( a * log(x) + b * log(one-x) - c ) / a ;
  386. } // incBetaFunc fraction
  387.  
  388. REAL incBetaFunc::value( REAL r )
  389. /**********************************************
  390.    Approximation of Incomplete Beta, based on
  391.    Press et al NRC p. 179.   Uses betaFraction.
  392.    Assumes c = logBeta(a,b) ;
  393. ************************************************/
  394. {
  395.    static char *mName = "value" ;
  396.    matFunc func( mName ) ; debugInfo( func ) ;
  397.    error = NOERROR ;
  398.    if ( a < 0 ) {
  399.       error = NGPAR ;
  400.       return zero ;
  401.    } else if ( r < zero || r > one ) {
  402.       error = BDARG ;
  403.       return zero ;
  404.    } else if ( r == zero )
  405.       return zero ;
  406.    else if ( r == one )
  407.       return one ;
  408.    else if ( r < ( a + one ) / ( a + b + two ) )
  409.       return fraction( r, a, b, c ) ;
  410.    return one - fraction( one - r, b, a, c ) ;
  411. } // incBetaFunc value
  412.  
  413. incBetaFunc& incBetaFunc::par( REAL newA, REAL newB )
  414. {
  415.    static char *mName = "par" ;
  416.    matFunc func( mName ) ; debugInfo( func ) ;
  417.    if ( newA <= zero || newB <= zero )
  418.       errorExit( "par", NGPAR ) ;
  419.    a = newA ;
  420.    b = newB ;
  421.    c = logBeta( a, b ) ;
  422.    return *this ;
  423. } // incBetaFunc par
  424.  
  425. matrix incBetaFunc::operator () ( const matrix& x, REAL newA, REAL newB )
  426. {
  427.    static char *mName = "incBeta(m)" ;
  428.    matFunc func( mName ) ; debugInfo( func ) ;
  429.    par( newA, newB ) ;
  430.    INDEX nr = x.nRows(), nc = x.nCols() ;
  431.    matrix z(nr,nc) ;
  432.    transform(z,x) ;
  433.    return z ;
  434. } // incBeta op ()
  435.  
  436. REAL incBetaFunc::operator () ( REAL r, REAL newA, REAL newB )
  437. {
  438.    static char *mName = "incBeta(r)" ;
  439.    matFunc func( mName ) ; debugInfo( func ) ;
  440.    par( newA, newB ) ;
  441.    REAL v = value(r) ;
  442.    if ( error )
  443.       errorExit( mName, error ) ;
  444.    return v ;
  445. } // incBeta op ()
  446.  
  447. incBetaFunc incBeta ;
  448.  
  449.  
  450. /************************************************************/
  451. /*                      incGamma class                    */
  452. /************************************************************/
  453.  
  454.  
  455. #ifdef __cplusplus
  456. incGammaFunc::incGammaFunc( void ) : matSpecFunc()
  457. #else
  458. incGammaFunc::incGammaFunc( void ) : ()
  459. #endif
  460. {
  461.    a = - 1.0 ;
  462.    c = 0 ;
  463. } // incGammaFunc
  464.  
  465. #ifdef __cplusplus
  466. incGammaFunc::incGammaFunc( incGammaFunc& fn )
  467.               : matSpecFunc(fn)
  468. #else
  469. incGammaFunc::incGammaFunc( incGammaFunc& fn )
  470.               : ( fn )
  471. #endif
  472. {
  473.    a = fn.a ;
  474.    c = fn.c ;
  475. } // incGammaFunc
  476.  
  477. void incGammaFunc::operator = ( incGammaFunc& fn )
  478. {
  479.    matSpecFunc::operator = ( fn ) ;
  480.    a = fn.a ;
  481.    c = fn.c ;
  482. } // incGammaFunc
  483.  
  484. incGammaFunc::~incGammaFunc( void ) {}
  485.  
  486. outFile& incGammaFunc::info( outFile& f ) M_CONST
  487. {
  488.    return specInfo( f, "incGam" ) ;
  489. } // info
  490.  
  491. REAL incGammaFunc::series( REAL r )
  492. /*******************************************************
  493.   return series representation incomplete gamma P(a,r).
  494.   Assumes c = logGamma(a).
  495. *******************************************************/
  496. {
  497.    static char *mName = "series" ;
  498.    matFunc func( mName ) ; debugInfo( func ) ;
  499.    REAL ser, del, aplus ;
  500.    int i ;
  501.    error = NOERROR ;
  502.    if ( r < zero ) {
  503.       error = BDARG ;
  504.       return zero ;
  505.    } // if
  506.    if ( r <= zero )
  507.       return zero ;
  508.    aplus = a ;
  509.    del = ser = one / a ;
  510.    i = 1 ;
  511.    do {
  512.       aplus += one ;
  513.       del   *= r / aplus ;
  514.       ser   += del ;
  515.    } while ( ( i <= maxIter ) && ( fabs(del) >= fabs(ser) * eps ) ) ;
  516.    if ( i >= maxIter ) {
  517.       error= NCONV ;
  518.       return zero ;
  519.    } // if
  520.    return ser *= exp( - r + a * log(r) - c ) ;
  521. } // incGammaFunc series
  522.  
  523. REAL incGammaFunc::fraction( REAL r )
  524. /**********************************************************
  525.    Use continued fraction expression for complement of
  526.    incomplete gamma function. See Press et al p.171.
  527.  
  528.    gammaComplement(a,r)=exp(-r).r^{a}.fraction/logGamma(a),
  529.  
  530.    where
  531.  
  532.                1    1-a   1   2-a   2
  533.    fraction = ---  ----- --- ----- --- ....
  534.                r+    1+   r+   1+   r+
  535.  
  536.    Use recurrence  fraction(n) = A(n) / B(n)
  537.  
  538.    where
  539.          A(n) = ( s(n) * A(n-1) + r(n) * A(n-2) ) * factor
  540.          B(n) = ( s(n) * B(n-1) + r(n) * B(n-2) ) * factor
  541.    and
  542.          A(-1) = 1, B(-1) = 0, A(0) = s(0), B(0) = 1.
  543.  
  544.    Here
  545.  
  546.          s(0) = 0, s(1) = r, r(0) = 0, r(1) = 1,
  547.  
  548.    so that
  549.  
  550.          A(1) = one * factor, B(1) = r * factor
  551.  
  552.    subsequently
  553.  
  554.          r(i) = k - a  if i = 2k,   k > 0
  555.          r(i) = k      if i = 2k+1,
  556.          s(i) = 1      if i = 2k,
  557.          s(i) = x      if i = 2k+1
  558.  
  559.    and factor is some scaling factor to avoid overflow
  560.  
  561.    Assumes c = logGamma(a).
  562.  
  563. *************************************************************/
  564. {
  565.    static char *mName = "fraction" ;
  566.    matFunc func( mName ) ; debugInfo( func ) ;
  567.    REAL old_gam = zero , factor = one;
  568.    REAL A0 = zero, A1 = one, B0 = one, B1 = r ;
  569.    REAL gam = one / r, z = zero, ma = zero - a, rfact ;
  570.    REAL maxZ = REAL(maxIter) ;
  571.  
  572.    do {
  573.       z  += one ;
  574.       ma += one ;
  575.       // two steps of recurrence replace A's & B's
  576.       A0 = ( A1 + ma * A0 ) * factor ; // i even
  577.       B0 = ( B1 + ma * B0 ) * factor ;
  578.       rfact = z * factor ;
  579.       A1 = r * A0 + rfact * A1 ; // i odd, A0 already rescaled
  580.       B1 = r * B0 + rfact * B1 ;
  581.       if ( B1 ) {
  582.          factor = one / B1 ;
  583.          old_gam = gam ;
  584.          gam = A1 * factor ;
  585.       } // if
  586.    } while ( z < maxZ && fabs( gam - old_gam ) > eps * fabs( gam ) ) ;
  587.    if ( z >= maxZ ) {
  588.       error = NCONV ;
  589.       return zero ;
  590.    } // if
  591.    return exp( -r + a * log(r) - c ) * gam ;
  592. } // incGammaFunc::fraction
  593.  
  594. REAL incGammaFunc::value( REAL r )
  595. /**********************************************
  596.    Approximation of Incomplete gamma, based on
  597.    Press et al NRC p. 172.
  598. ************************************************/
  599. {
  600.    static char *mName = "value" ;
  601.    matFunc func( mName ) ; debugInfo( func ) ;
  602.    if ( a < 0 ) {
  603.       error =  NGPAR ;
  604.       return zero ;
  605.    } else if ( r < zero ) {
  606.       error = BDARG ;
  607.       return zero ;
  608.    } else if ( r == zero )
  609.       return zero ;
  610.    else if ( r < ( a + one ) )
  611.       return series( r ) ;
  612.    return one - fraction( r ) ;
  613. } // incGammaFunc value
  614.  
  615. incGammaFunc& incGammaFunc::par( REAL newA )
  616. {
  617.    static char *mName = "par" ;
  618.    matFunc func( mName ) ; debugInfo( func ) ;
  619.    if ( newA <= zero )
  620.       errorExit( "par", NGPAR ) ;
  621.    a = newA ;
  622.    c = logGamma( a ) ;
  623.    return *this ;
  624. } // incGammaFunc par
  625.  
  626. matrix incGammaFunc::operator () ( const matrix& x, REAL newA )
  627. {
  628.    static char *mName = "incGamma(m)" ;
  629.    matFunc func( mName ) ; debugInfo( func ) ;
  630.    par( newA ) ;
  631.    INDEX nr = x.nRows(), nc = x.nCols() ;
  632.    matrix z(nr,nc) ;
  633.    transform(z,x) ;
  634.    return z ;
  635. } // incGamma op ()
  636.  
  637. REAL incGammaFunc::operator () ( REAL r, REAL newA )
  638. {
  639.    static char *mName = "incGamma(r)" ;
  640.    matFunc func( mName ) ; debugInfo( func ) ;
  641.    par( newA ) ;
  642.    REAL v = value(r) ;
  643.    if ( error )
  644.       errorExit( mName, error ) ;
  645.    return v ;
  646. } // incGamma op ()
  647.  
  648. incGammaFunc incGamma ;
  649.  
  650.  
  651.