Changeset 785

Timestamp:

01/16/07 00:03:11 (2 years ago)

Author:

peter

Message:

source/thirdparty_bevington.cpp:
source/thirdparty.cpp:
source/helike_cs.cpp:

• Delete and replace software with incompatible or unclear licenses.
  Replacement routines are taken from the Cephes library.
  
• remove gammafun and replace with standard mathlib function tgamma.
  

source/hydro_bauman.cpp:

• Move factorial and log10_factorial into thirdparty.cpp to make
  them globally available.
  

source/service.cpp:

• Replace e2 with better version using recursion relation. The old
  version had dreadful accuracy.
  

source/sanity_check.cpp:

• Make check on e2 stricter now that it is better.
  

source/atom_pop5.cpp:
source/ion_solver.cpp:
source/conv_eden_ioniz.cpp:
source/cont_gaunt.cpp:
source/thirdparty_lapack.cpp:
source/hydrolevelpop.cpp:
source/atom_seq_beryllium.cpp:
source/helike_level.cpp:
source/prt_lines_helium.cpp:
source/atom_leveln.cpp:
source/mole_h_step.cpp:
source/atom_oi.cpp:
source/atom_feii.cpp:
source/mole_co_solve.cpp:
source/helike_recom.cpp:

• Update includes.
  
• Call replaced routines from thirdparty.cpp.
  

source/bevington.h:
source/lapack.h:

• Deleted and merged into thirdparty.h.
  

source/thirdparty.h:

• Included the above, added prototypes for new routines from Cephes.
  

Location:

trunk/source

Files:

• atom_feii.cpp (modified) (1 diff)
• atom_leveln.cpp (modified) (1 diff)
• atom_oi.cpp (modified) (1 diff)
• atom_pop5.cpp (modified) (1 diff)
• atom_seq_beryllium.cpp (modified) (1 diff)
• bevington.h (deleted)
• cont_gaunt.cpp (modified) (1 diff)
• conv_eden_ioniz.cpp (modified) (1 diff)
• helike_cs.cpp (modified) (9 diffs)
• helike_level.cpp (modified) (1 diff)
• helike_recom.cpp (modified) (3 diffs)
• hydrolevelpop.cpp (modified) (1 diff)
• hydro_bauman.cpp (modified) (19 diffs)
• ion_solver.cpp (modified) (1 diff)
• lapack.h (deleted)
• mole_co_solve.cpp (modified) (1 diff)
• mole_h_step.cpp (modified) (1 diff)
• prt_lines_helium.cpp (modified) (1 diff)
• sanity_check.cpp (modified) (4 diffs)
• service.cpp (modified) (1 diff)
• thirdparty.cpp (modified) (16 diffs)
• thirdparty.h (modified) (1 diff)
• thirdparty_bevington.cpp (modified) (3 diffs)
• thirdparty_lapack.cpp (modified) (1 diff)

trunk/source/atom_feii.cpp

@@ -51 +51 @@
 #include "lines_service.h"
 #include "ipoint.h"
-#include "lapack.h"
+#include "thirdparty.h"
 #include "hydrogenic.h"
 #include "lines.h"

trunk/source/atom_leveln.cpp

@@ -8 +8 @@
 #include "dense.h"
 #include "trace.h"
-#include "lapack.h"
+#include "thirdparty.h"
 #include "atoms.h"
 

trunk/source/atom_oi.cpp

@@ -14 +14 @@
 #include "phycon.h"
 #include "lines_service.h"
-#include "lapack.h"
+#include "thirdparty.h"
 #include "atoms.h"
 

trunk/source/atom_pop5.cpp

@@ -6 +6 @@
 #include "phycon.h"
 #include "dense.h"
-#include "lapack.h"
+#include "thirdparty.h"
 #include "atoms.h"
 

trunk/source/atom_seq_beryllium.cpp

@@ -5 +5 @@
 #include "phycon.h"
 #include "lines_service.h"
-#include "lapack.h"
+#include "thirdparty.h"
 #include "dense.h"
 #include "cooling.h"

trunk/source/cont_gaunt.cpp

@@ -349 +349 @@
 
                 /* This is the coefficient in equation 8 in Sutherland  */
-                /* Karzas and Latter give gammafun(2.*i+2.), Sutherland gives gammafun(2.*i+1.) */
-                ICoef = 0.25*pow(-chi, (double)i+1.)*exp( 1.5708*fabs(etai-etaf) )/gammafun(2.*i+1.);
+                /* Karzas and Latter give tgamma(2.*i+2.), Sutherland gives tgamma(2.*i+1.)     */
+                ICoef = 0.25*pow(-chi, (double)i+1.)*exp( 1.5708*fabs(etai-etaf) )/factorial(2*i);
                 GammaProduct = cdgamma(complex<double>(i+1.,etai))*cdgamma(complex<double>(i+1.,etaf));
                 ICoef *= abs(GammaProduct);

trunk/source/conv_eden_ioniz.cpp

@@ -7 +7 @@
 #include "trace.h"
 #include "thermal.h"
-#include "bevington.h"
+#include "thirdparty.h"
 #include "phycon.h"
 #include "conv.h"

trunk/source/helike_cs.cpp

@@ -25 +25 @@
 
 /* returns thermally-averaged Seaton 62 collision strength. */
-static double S62_Therm_ave_coll_str( double EProjectile_eV );
+STATIC double S62_Therm_ave_coll_str( double EProjectile_eV );
 
 /* all of these are used in the calculation of Stark collision strengths
  * following the algorithms in Vrinceanu & Flannery (2001). */
-static double Therm_ave_coll_str_int_VF01( double EProjectileRyd);
-static double collision_strength_VF01( double velOrEner, bool lgParamIsRedVel );
-static double L_mix_integrand_VF01( double alpha );
-static double StarkCollTransProb_VF01( long int n, long int l, long int lp, double alpha, double deltaPhi);
-
-/* returns the value of elliptical integral K */
-static double ellip_int_K( double k );
+STATIC double Therm_ave_coll_str_int_VF01( double EProjectileRyd);
+STATIC double collision_strength_VF01( double velOrEner, bool lgParamIsRedVel );
+STATIC double L_mix_integrand_VF01( double alpha );
+STATIC double StarkCollTransProb_VF01( long int n, long int l, long int lp, double alpha, double deltaPhi);
 
 static long     int global_n, global_l, global_l_prime, global_s, global_z, global_Collider;
@@ -1238 +1235 @@
 
 /* The integrand for calculating the thermal average of collision strengths */
-static double S62_Therm_ave_coll_str( double proj_energy )
+STATIC double S62_Therm_ave_coll_str( double proj_energy )
 {
 
@@ -1335 +1332 @@
         zeta_of_betaone = zOverB2 * POW2(betaone);
         
-        /* cs1 = betanot * gsl_sf_bessel_K0_e(betanot) * gsl_sf_bessel_K1_e(betanot); */
-        cs2 = 0.5*zeta_of_betaone + betaone *
-                gsl_sf_bessel_K0_e(betaone) * gsl_sf_bessel_K1_e(betaone);
+        /* cs1 = betanot * bessel_k0(betanot) * bessel_k1(betanot); */
+        cs2 = 0.5*zeta_of_betaone + betaone * bessel_k0(betaone) * bessel_k1(betaone);
                 
         /* cross_section = MIN2(cs1, cs2); */
@@ -1511 +1507 @@
 
 /* The integrand for calculating the thermal average of collision strengths */
-static double Therm_ave_coll_str_int_VF01( double EOverKT )
+STATIC double Therm_ave_coll_str_int_VF01( double EOverKT )
 {
         double integrand;
@@ -1524 +1520 @@
 }
 
-static double collision_strength_VF01( double velOrEner, bool lgParamIsRedVel )
+STATIC double collision_strength_VF01( double velOrEner, bool lgParamIsRedVel )
 {
 
@@ -1674 +1670 @@
 }
 
-static double L_mix_integrand_VF01( double alpha )
+STATIC double L_mix_integrand_VF01( double alpha )
 {
         double integrand, deltaPhi, b;
@@ -1703 +1699 @@
 }
 
-static double StarkCollTransProb_VF01( long n, long l, long lp, double alpha, double deltaPhi)
+STATIC double StarkCollTransProb_VF01( long n, long l, long lp, double alpha, double deltaPhi)
 {
         double probability;
@@ -1776 +1772 @@
                                 if( A < 0. )
                                 {
-                                        probability *= ellip_int_K( sqrt(B/(B-A)) );
+                                        probability *= ellpk( -A/(B-A) );
                                         probability /= sqrt( B - A );
                                 }
                                 else
                                 {
-                                        probability *= ellip_int_K( sqrt((B-A)/B) );
+                                        probability *= ellpk( A/B );
                                         probability /= sqrt( B );
                                 }
@@ -1794 +1790 @@
 
 }
-
-
-/* This is a modified version of two routines
- * in the GNU science library, version 1.4:
- *              gsl_sf_ellint_Kcomp_e
- *              gsl_sf_ellint_RF_e
- * The routines include references:
- *              [Carlson, Numer. Math. 33 (1979) 1, (4.5)]
- *              [Abramowitz+Stegun, 17.3.33]
- */
-static double ellip_int_K( double k )
-{
-        double result;
-
-        DEBUG_ENTRY( "ellip_int_K()" );
-
-        ASSERT( k*k < 1.0 );
-  
-        if(k*k >= 1.0 - sqrt(DBL_EPSILON))
-        {
-                /* [Abramowitz+Stegun, 17.3.33] */
-                const double y = 1.0 - k*k;
-                const double a[] = { 1.38629436112, 0.09666344259, 0.03590092383 };
-                const double b[] = { 0.5, 0.12498593597, 0.06880248576 };
-                const double ta = a[0] + y*(a[1] + y*a[2]);
-                const double tb = -log(y) * (b[0] * y*(b[1] + y*b[2]));
-                result = ta + tb;
-        }
-        else
-        {
-                double y = 1.0 - k*k;
-                
-                /*const double lolim = 5.0 * DBL_MIN;
-                const double uplim = 0.2 * DBL_MAX;*/
-
-                const double c1 = 1.0 / 24.0;
-                const double c2 = 3.0 / 44.0;
-                const double c3 = 1.0 / 14.0;
-                double xn = 0.0;
-                double yn_ = y;
-                double zn = 1.0;
-                double mu =0., xndev=0., yndev=0., zndev=0., e2_, e3, s;
-                /* yn and e2 hide VS.net internal symbol yn and cddefines symbol e2, respectively. 
-                 * An underscore at the end prevents this collision. */
-                
-                ASSERT( y >= 5.0 * DBL_MIN );
-                ASSERT( y < 0.2 * DBL_MAX );
-
-                while(1)
-                {
-                        double epslon, lamda;
-                        double xnroot, ynroot, znroot;
-                        mu = (xn + yn_ + zn) / 3.0;
-                        ASSERT( mu > 0. );
-                        xndev = 2.0 - (mu + xn) / mu;
-                        yndev = 2.0 - (mu + yn_) / mu;
-                        zndev = 2.0 - (mu + zn) / mu;
-                        epslon = MAX3(fabs(xndev), fabs(yndev), fabs(zndev));
-                        if(epslon < 0.001) break;
-                        xnroot = sqrt(xn);
-                        ynroot = sqrt(yn_);
-                        znroot = sqrt(zn);
-                        lamda = xnroot * (ynroot + znroot) + ynroot * znroot;
-                        xn = (xn + lamda) * 0.25;
-                        yn_ = (yn_ + lamda) * 0.25;
-                        zn = (zn + lamda) * 0.25;
-                        ASSERT( xn>0. && yn_>0. && zn>0. );
-                }
-                e2_ = xndev * yndev - zndev * zndev;
-                e3 = xndev * yndev * zndev;
-                s = 1.0 + (c1 * e2_ - 0.1 - c2 * e3) * e2_ + c3 * e3;
-                result = s / sqrt(mu);
-        }
-
-        DEBUG_EXIT( "ellip_int_K()" );
-
-        return result;
-
-}
 /*lint +e661 possible access of out-of-bounds pointer */
-

trunk/source/helike_level.cpp

@@ -4 +4 @@
 /*HeLikeError fills uncertainty arrays */
 #include "cddefines.h" 
-#include "lapack.h"
+#include "thirdparty.h"
 #include "trace.h"
 #include "secondaries.h"

trunk/source/helike_recom.cpp

@@ -35 +35 @@
 #include "ionbal.h"
 #include "atmdat.h"
-#include "bevington.h"
-/* This will be needed if the Peach algorithms are turned on. */
-/* #include "thirdparty.h" */
+#include "thirdparty.h"
 
 /*lint -e662 creation of  out of bound pointer */
@@ -1981 +1979 @@
         /* I'm pretty sure Peach has this wrong...and the exponent
          * should be -1/2 instead of 1/2. */
-        K = 1./( nu * sqrt( zeta * gammafun(nu+l+1.) * gammafun(nu-l) ) );
+        K = 1./( nu * sqrt( zeta * tgamma(nu+l+1.) * tgamma(nu-l) ) );
         Sum = 1.;
         b = 1.;
@@ -2047 +2045 @@
         energy = POW2(kprime/(double)nelem);
         
-        Fkl = gsl_sf_bessel_Jn_e(2*lp+1,bfa)+ (energy/12)*
-                ( lp*(lp+1)*(2*lp+1)*gsl_sf_bessel_Jn_e(2*lp+1,bfa)-
-                3*(lp+1)*(2*rho)*gsl_sf_bessel_Jn_e(2*lp+3,bfa)+
-                (pow(2.*rho, 1.5))*gsl_sf_bessel_Jn_e(2*lp+4,bfa) );
+        Fkl = bessel_jn(2*lp+1,bfa)+ (energy/12)*
+                ( lp*(lp+1)*(2*lp+1)*bessel_jn(2*lp+1,bfa)-
+                3*(lp+1)*(2*rho)*bessel_jn(2*lp+3,bfa)+
+                (pow(2.*rho, 1.5))*bessel_jn(2*lp+4,bfa) );
         Fkl *= sqrt(PI*rho);
 
-        Hkl = gsl_sf_bessel_Yn_e(2*lp+1,bfa)+ (energy/12)*
-                ( lp*(lp+1)*(2*lp+1)*gsl_sf_bessel_Yn_e(2*lp+1,bfa)-
-                3*(lp+1)*(2*rho)*gsl_sf_bessel_Yn_e(2*lp+3,bfa)+
-                (pow(2.*rho, 1.5))*gsl_sf_bessel_Yn_e(2*lp+4,bfa) );
+        Hkl = bessel_yn(2*lp+1,bfa)+ (energy/12)*
+                ( lp*(lp+1)*(2*lp+1)*bessel_yn(2*lp+1,bfa)-
+                3*(lp+1)*(2*rho)*bessel_yn(2*lp+3,bfa)+
+                (pow(2.*rho, 1.5))*bessel_yn(2*lp+4,bfa) );
         Hkl *= pow(1.-sexp(tau_l*rho), 2.*lp+1.) * sqrt(PI*rho);
 

trunk/source/hydrolevelpop.cpp

@@ -13 +13 @@
 #include "dense.h"
 #include "trace.h"
-#include "lapack.h"
+#include "thirdparty.h"
 #include "hydrogenic.h"
 #include "dynamics.h"

trunk/source/hydro_bauman.cpp

@@ -15 +15 @@
 #include "cddefines.h"
 #include "physconst.h"
+#include "thirdparty.h"
 #include "hydro_bauman.h"
 
@@ -527 +528 @@
 
 inline double local_product( double K , long int lp );
-inline double log10_factorial( long int n );
 inline double log10_prodxx( long int lp, double Ksqrd );
 
@@ -560 +560 @@
 
 static const double PHYSICAL_CONSTANT_TWO = 4./3.*PI*FINE_STRUCTURE*pow2(BohrRadiusCM);
-
-/************************************************************************************************/
-/* pre-calculated factorials                                                                    */
-/************************************************************************************************/
-
-static const int NPRE_FACTORIAL = 171;
-static const double pre_factorial[NPRE_FACTORIAL] =
-{
-        1.00000000000000000000e+00,
-        1.00000000000000000000e+00,
-        2.00000000000000000000e+00,
-        6.00000000000000000000e+00,
-        2.40000000000000000000e+01,
-        1.20000000000000000000e+02,
-        7.20000000000000000000e+02,
-        5.04000000000000000000e+03,
-        4.03200000000000000000e+04,
-        3.62880000000000000000e+05,
-        3.62880000000000000000e+06,
-        3.99168000000000000000e+07,
-        4.79001600000000000000e+08,
-        6.22702080000000000000e+09,
-        8.71782912000000000000e+10,
-        1.30767436800000000000e+12,
-        2.09227898880000000000e+13,
-        3.55687428096000000000e+14,
-        6.40237370572800000000e+15,
-        1.21645100408832000000e+17,
-        2.43290200817664000000e+18,
-        5.10909421717094400000e+19,
-        1.12400072777760768000e+21,
-        2.58520167388849766400e+22,
-        6.20448401733239439360e+23,
-        1.55112100433309859840e+25,
-        4.03291461126605635592e+26,
-        1.08888694504183521614e+28,
-        3.04888344611713860511e+29,
-        8.84176199373970195470e+30,
-        2.65252859812191058647e+32,
-        8.22283865417792281807e+33,
-        2.63130836933693530178e+35,
-        8.68331761881188649615e+36,
-        2.95232799039604140861e+38,
-        1.03331479663861449300e+40,
-        3.71993326789901217463e+41,
-        1.37637530912263450457e+43,
-        5.23022617466601111726e+44,
-        2.03978820811974433568e+46,
-        8.15915283247897734264e+47,
-        3.34525266131638071044e+49,
-        1.40500611775287989839e+51,
-        6.04152630633738356321e+52,
-        2.65827157478844876773e+54,
-        1.19622220865480194551e+56,
-        5.50262215981208894940e+57,
-        2.58623241511168180614e+59,
-        1.24139155925360726691e+61,
-        6.08281864034267560801e+62,
-        3.04140932017133780398e+64,
-        1.55111875328738227999e+66,
-        8.06581751709438785585e+67,
-        4.27488328406002556374e+69,
-        2.30843697339241380441e+71,
-        1.26964033536582759243e+73,
-        7.10998587804863451749e+74,
-        4.05269195048772167487e+76,
-        2.35056133128287857145e+78,
-        1.38683118545689835713e+80,
-        8.32098711274139014271e+81,
-        5.07580213877224798711e+83,
-        3.14699732603879375200e+85,
-        1.98260831540444006372e+87,
-        1.26886932185884164078e+89,
-        8.24765059208247066512e+90,
-        5.44344939077443063905e+92,
-        3.64711109181886852801e+94,
-        2.48003554243683059915e+96,
-        1.71122452428141311337e+98,
-        1.19785716699698917933e+100,
-        8.50478588567862317347e+101,
-        6.12344583768860868500e+103,
-        4.47011546151268434004e+105,
-        3.30788544151938641157e+107,
-        2.48091408113953980872e+109,
-        1.88549470166605025458e+111,
-        1.45183092028285869606e+113,
-        1.13242811782062978295e+115,
-        8.94618213078297528506e+116,
-        7.15694570462638022794e+118,
-        5.79712602074736798470e+120,
-        4.75364333701284174746e+122,
-        3.94552396972065865030e+124,
-        3.31424013456535326627e+126,
-        2.81710411438055027626e+128,
-        2.42270953836727323750e+130,
-        2.10775729837952771662e+132,
-        1.85482642257398439069e+134,
-        1.65079551609084610774e+136,
-        1.48571596448176149700e+138,
-        1.35200152767840296226e+140,
-        1.24384140546413072522e+142,
-        1.15677250708164157442e+144,
-        1.08736615665674307994e+146,
-        1.03299784882390592592e+148,
-        9.91677934870949688836e+149,
-        9.61927596824821198159e+151,
-        9.42689044888324774164e+153,
-        9.33262154439441526381e+155,
-        9.33262154439441526388e+157,
-        9.42594775983835941673e+159,
-        9.61446671503512660515e+161,
-        9.90290071648618040340e+163,
-        1.02990167451456276198e+166,
-        1.08139675824029090008e+168,
-        1.14628056373470835406e+170,
-        1.22652020319613793888e+172,
-        1.32464181945182897396e+174,
-        1.44385958320249358163e+176,
-        1.58824554152274293982e+178,
-        1.76295255109024466316e+180,
-        1.97450685722107402277e+182,
-        2.23119274865981364576e+184,
-        2.54355973347218755612e+186,
-        2.92509369349301568964e+188,
-        3.39310868445189820004e+190,
-        3.96993716080872089396e+192,
-        4.68452584975429065488e+194,
-        5.57458576120760587943e+196,
-        6.68950291344912705515e+198,
-        8.09429852527344373681e+200,
-        9.87504420083360135884e+202,
-        1.21463043670253296712e+205,
-        1.50614174151114087918e+207,
-        1.88267717688892609901e+209,
-        2.37217324288004688470e+211,
-        3.01266001845765954361e+213,
-        3.85620482362580421582e+215,
-        4.97450422247728743840e+217,
-        6.46685548922047366972e+219,
-        8.47158069087882050755e+221,
-        1.11824865119600430699e+224,
-        1.48727070609068572828e+226,
-        1.99294274616151887582e+228,
-        2.69047270731805048244e+230,
-        3.65904288195254865604e+232,
-        5.01288874827499165889e+234,
-        6.91778647261948848943e+236,
-        9.61572319694108900019e+238,
-        1.34620124757175246000e+241,
-        1.89814375907617096864e+243,
-        2.69536413788816277557e+245,
-        3.85437071718007276916e+247,
-        5.55029383273930478744e+249,
-        8.04792605747199194159e+251,
-        1.17499720439091082343e+254,
-        1.72724589045463891049e+256,
-        2.55632391787286558753e+258,
-        3.80892263763056972532e+260,
-        5.71338395644585458806e+262,
-        8.62720977423324042775e+264,
-        1.31133588568345254503e+267,
-        2.00634390509568239384e+269,
-        3.08976961384735088657e+271,
-        4.78914290146339387432e+273,
-        7.47106292628289444390e+275,
-        1.17295687942641442768e+278,
-        1.85327186949373479574e+280,
-        2.94670227249503832518e+282,
-        4.71472363599206132029e+284,
-        7.59070505394721872577e+286,
-        1.22969421873944943358e+289,
-        2.00440157654530257674e+291,
-        3.28721858553429622598e+293,
-        5.42391066613158877297e+295,
-        9.00369170577843736335e+297,
-        1.50361651486499903974e+300,
-        2.52607574497319838672e+302,
-        4.26906800900470527345e+304,
-        7.25741561530799896496e+306
-        /*,1.24101807021766782298e+309*/
-};
 
 /************************Start of program***************************/
@@ -1221 +1040 @@
         double Ksqrd = K * K;
 
-        double ld1 = pre_factorial[ (2*n) - 1 ];
+        double ld1 = factorial( 2*n - 1 );
         double ld2 = powi((double)(4*n), n);
         double ld3 = exp(-(double)(2 * n));
@@ -1356 +1175 @@
         /******************************/
 
-        ld1 = log10_factorial( ( (2 * n) - 1) );
+        ld1 = lfactorial( 2*n - 1 );
         ASSERT( ld1 >= 0. );
 
@@ -1482 +1301 @@
         {
                 printf( "BadMagic: l and l' do NOT satisfy dipole requirements.\n\n" );
-                puts( "[Stop in log10_factorial]" );
+                puts( "[Stop in bhG_mx]" );
                 cdEXIT(EXIT_FAILURE);
         }
@@ -2294 +2113 @@
         DEBUG_ENTRY( "bhg()" );
 
-        double ld1 = pre_factorial[ (n + l) ];
-        double ld2 = pre_factorial[ ((n - l) - 1) ];
+        double ld1 = factorial( n + l );
+        double ld2 = factorial( n - l - 1 );
         double ld3 = (ld1 / ld2);
 
@@ -2377 +2196 @@
         /**************************************************************************************/
 
-        double ld1 = log10_factorial( (n + l) );
-        double ld2 = log10_factorial( ((n - l) - 1) );
+        double ld1 = lfactorial( n + l );
+        double ld2 = lfactorial( n - l - 1 );
 
         /**********************************************************************/
@@ -3053 +2872 @@
                         }
                         i2 = (2*l - 1);
-                        d3 = pre_factorial[ i2 ];
+                        d3 = factorial( i2 );
                         d4 = d2/d3;
                         d4 = sqrt( d4 );
@@ -3395 +3214 @@
                         /**************************************/
 
-                        ld2 = log10_factorial( 2*l - 1 );
+                        ld2 = lfactorial( 2*l - 1 );
 
                         ASSERT( ((2*l)+1) >= 0);
@@ -3771 +3590 @@
         /****************************************************************/
 
-        d0 = pre_factorial[ i1 ];
+        d0 = factorial( i1 );
         d1 = 4. * d0;
         d2 = 1. / d1;
@@ -3825 +3644 @@
                 cdEXIT(EXIT_FAILURE);
         }
-        d1 = pre_factorial[ i1 ];
+        d1 = factorial( i1 );
 
         i2 = np + l - 1;
@@ -3834 +3653 @@
                 cdEXIT(EXIT_FAILURE);
         }
-        d2 = pre_factorial[ i2 ];
+        d2 = factorial( i2 );
 
         i3 = n - l - 1;
@@ -3843 +3662 @@
                 cdEXIT(EXIT_FAILURE);
         }
-        d3 = pre_factorial[ i3 ];
+        d3 = factorial( i3 );