/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:32:02 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dbinom.h" #include #include #include /* PARAMETER translations */ #define MAXL 12 #define MAXN 150 /* end of PARAMETER translations */ double /*FUNCTION*/ dbinom( long n, long k) { long int _l0, i, ierr, k1, k2, k3, ki, ni, nmk; static long int lbnd[MAXL]; double dbinom_v, fni, tp, tp1; static double big, biglog, bigtst, exerr, fac[MAXN]; static long nc = 0; static long kbnd = 1; static char ermsg[2][33]={"Bad value for N or K, ","Result would overflow, "}; static long primes[33-(3)+1]={13,17,19,23,29,31,37,41,43,47,53, 59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139, 149,139}; /* OFFSET Vectors w/subscript range: 1 to dimension */ double *const Fac = &fac[0] - 1; long *const Lbnd = &lbnd[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 1998-01-23 DBINOM Krogh Work around for HP compiler bug. *>> 1997-04-14 DBINOM Krogh Added external statement, aint => anint. *>> 1995-12-14 DBINOM Krogh Initial Code *--D replaces "?": ?BINOM, ?LGAMA * * ( N } * Computes the binomial coefficient ( }. Require 0 .le. K .le. N. * ( K ) * * **************** Variable Definitions ******************************** * * D1MACH D1MACH returns the largest floating point number * DBINOM The result returned for the Binomial Coefficient. If an error * returns, this is set to -1. * DLGAMA MATH77 library for computing the log of the gamma function. * BIG Used in testing for possible overflow. * BIGTST BIG / MAXN, used to test for overflow when computing factorials * BIGLOG Used to checking for overflow before exponentiating. * EXERR If the result is bigger than EXERR some extra work is done * to eliminate the effect of small round off errors. * FAC Array of factorials. For N .le. LBND(1), N! = FAC(N). For * LBND(1) .lt. N .le. LBND(2), N! = FAC(LBND(1)) * FAC(N). For * LBND(j) .lt. N .le. LBND(j+1), N! = FAC(LBND(1)) * ... * * FAC(LBND(j)) * FAC(N).] * FNI Floating value for NI. * I Temporary index. * IERR Index for error message. * K Formal argumtent, see description above. * K1,K2,K3 Indices used in untangling factorials. K1 is originally for * NI, K2 for NMK, and K3 for KI. * KI Intenal value for K = min (K, N - K) * LBND Array used to obtain large factorials, see FAC above. * KBND Last index used in LBND. * MAXL Size of the array LBND. This is a little larger than needed * to handle N = MAXN in single precision IEEE arithmetic. May * need to be bigger if MAXN is increased. * MAXN Size of the array FAC. This is probably larger than needed for * single precision in most cases. * N Formal argumtent, see description above. * NC Highest factorial that has been computed. * NI Internal value for N. * NMK NI - KI * PRIMES A table of primes. The last entry must be equal to the third * to the last entry. * TP Used for temporary storage. * TP1 Used for temporary storage. * * **************** Variable Declarations ******************************* * */ /* 151,157,163,167, 173,179,181,191,193,197,199 * * **************** Start of Executable Code **************************** * */ ni = n; if (ni < 0) goto L_300; ki = min( k, ni - k ); if (ki <= 2) { if (ki < 1) { tp = 1.e0; if (ki != 0) goto L_300; } else { fni = ni; if (ki == 1) { tp = fni; } else { tp = .5e0*(fni*(fni - 1.e0)); } } } else { nmk = ni - ki; if (ni > nc) { if (nc == 0) { big = DBL_MAX; bigtst = big/(double)( MAXN ); biglog = log( big ) - .1e0; exerr = .01e0/DBL_EPSILON; Fac[1] = 1.e0; Fac[2] = 2.e0; nc = 2; } if (ni > MAXN) { tp = dlgama( (double)( ni + 1 ) ) - dlgama( (double)( ki + 1 ) ) - dlgama( (double)( nmk + 1 ) ); if (tp > biglog) goto L_320; tp = exp( tp ); goto L_250; } tp = nc; for (i = nc; i <= (ni - 1); i++) { tp += 1.e0; if (Fac[i] >= bigtst) { Lbnd[kbnd] = i; kbnd += 1; Fac[i + 1] = tp; } else { Fac[i + 1] = tp*Fac[i]; } } nc = ni; Lbnd[kbnd] = ni; } if (ni <= Lbnd[1]) { tp = Fac[ni]/(Fac[ki]*Fac[nmk]); } else if (nmk <= Lbnd[1]) { tp = Fac[ni]*(Fac[Lbnd[1]]/(Fac[ki])/Fac[nmk]); if (ni > Lbnd[2]) tp *= Fac[Lbnd[2]]; } else { k2 = 1; k3 = 1; i = 1; for (k1 = 1; k1 <= MAXL; k1++) { if (ni <= Lbnd[k1]) goto L_140; if (nmk > Lbnd[k1]) { k2 += 1; if (ki > Lbnd[k1]) k3 += 1; } } L_140: if (k2 == k1) { tp = Fac[ni]/(Fac[nmk]*Fac[ki]); } else { tp = Fac[Lbnd[k2]]/Fac[nmk]; L_150: k2 += 1; L_160: if (tp > 1.e0) { if (i < k3) { tp /= Fac[Lbnd[i]]; i += 1; goto L_160; } tp /= Fac[ki]; for (i = k2; i <= (k1 - 1); i++) { if (tp > big/Fac[Lbnd[i]]) goto L_320; tp *= Fac[Lbnd[i]]; } if (tp > big/Fac[ni]) goto L_320; tp *= Fac[ni]; } else if (k2 != k1) { tp *= Fac[Lbnd[k2]]; goto L_150; } else { tp = tp*Fac[ni]/Fac[ki]; } } } if (tp >= exerr) { i = .2199e0*(double)( ni ); L_200: if (primes[i-(3)] <= nmk) { i += 1; goto L_200; } if (primes[i-(3)] <= ni) { tp1 = primes[i-(3)]; if (primes[i + 1-(3)] <= ni) tp1 *= primes[i + 1-(3)]; /* anint(TP) avoided below due to bugs in recent HP Fortran compilers. * TP = TP1 * anint(TP / TP1) */ tp = .5e0 + (tp/tp1); tp = tp1*(tp - fmod( tp, 1.e0 )); } } } L_250: ; /* TP = anint(TP) */ tp = (tp + .5e0) - fmod( tp + .5e0, 1.e0 ); dbinom_v = tp; return( dbinom_v ); /* Error processing -- Default is to stop in IERV1. */ L_300: ierr = 1; goto L_340; L_320: ierr = 2; L_340: ierm1( "DBINOM", ierr, 2, (char*)ermsg[ierr - 1], " N ", n, ',' ); ierv1( " K ", k, '.' ); dbinom_v = -1.e0; return( dbinom_v ); } /* end of function */