karbon
vglobal.cc00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include <math.h>
00023
00024 #include "vglobal.h"
00025
00026
00027 int
00028 VGlobal::binomialCoeff( unsigned n, unsigned k )
00029 {
00030 return
00031 static_cast<int>(
00032 0.5 +
00033 exp(
00034 factorialLn( n ) -
00035 factorialLn( k ) -
00036 factorialLn( n - k ) ) );
00037 }
00038
00039 double
00040 VGlobal::factorialLn( unsigned n )
00041 {
00042 const unsigned cacheSize = 100;
00043
00044
00045 static double cache[ cacheSize ];
00046
00047
00048 if( n <= 1 )
00049 return 0.0;
00050
00051 if( n <= cacheSize - 1 )
00052 {
00053 return cache[ n ]
00054 ? cache[ n ]
00055 : ( cache[ n ] = gammaLn( n + 1.0 ) );
00056 }
00057 else
00058 {
00059 return gammaLn( n + 1.0 );
00060 }
00061 }
00062
00063 double
00064 VGlobal::gammaLn( double x )
00065 {
00066 static double coeff[ 6 ] =
00067 {
00068 76.18009172947146,
00069 -86.50532032941677,
00070 24.01409824083091,
00071 -1.231739572450155,
00072 0.1208650973866179e-2,
00073 -0.5395239384953e-5
00074 };
00075
00076 double y = x;
00077
00078 double tmp = x + 5.5;
00079 tmp -= ( x + 0.5 ) * log( tmp );
00080
00081 double ser = 1.000000000190015;
00082
00083 for( int i = 0; i < 5; ++i )
00084 {
00085 ser += coeff[ i ] / ++y;
00086 }
00087
00088 return -tmp + log( 2.5066282746310005 * ser / x );
00089 }
00090
|