diff --git a/miscfunc.cc b/miscfunc.cc index 33f13fe..d92edce 100644 --- a/miscfunc.cc +++ b/miscfunc.cc @@ -18,6 +18,7 @@ #include #include +#include #include #include "miscfunc.h" @@ -57,6 +58,26 @@ return a[n]; } +//n! in double precision for the cases where the exact factorial() would overflow +#define MAXDOUBLEFACT 170 +double double_factorial(const int n) +{ +static int ntop=MAXFACT; +static double a[MAXDOUBLEFACT+1]; + +int j; +if (n < 0) laerror("negative argument of double_factorial"); +if (n > MAXDOUBLEFACT) laerror("overflow in double_factorial"); +if (n <= MAXFACT) return (double)factorial(n); +if (ntop == MAXFACT) a[MAXFACT]=(double)factorial(MAXFACT); + while (ntopd;n--) p *= n; -value=p/factorial(k); +/* value = C(n-k+i,i) after step i is exact; divide by the common factor first to avoid overflow of intermediates */ +value=1; +for(int i=1; i<=k; ++i) + { + unsigned long long m=n-k+i; + unsigned long long g=ullgcd(m,i); + m /= g; + value /= (i/g); /*exact, since value*m is divisible by i and m/g is coprime to i/g*/ + if(value > std::numeric_limits::max()/m) laerror("overflow in binom"); + value *= m; + } if(n<=MAXBINOM) ibitab[ind]=value; return value; } @@ -201,436 +234,30 @@ return y ; -static void splintx(double ary[][3], int n,double x,double *y) +/* Euler-Maclaurin summation: zeta(x)= sum_{k=1}^{N-1} k^-x + N^(1-x)/(x-1) + N^-x/2 + sum_{j=1}^m B_2j/(2j)! x(x+1)...(x+2j-2) N^(-x-2j+1) */ +/* accurate to about 1e-15 for -1.02<=x<9.1 (the pole term 1/(x-1) is exact) */ +static double zeta_euler_maclaurin(double x) { - int klo,khi,k; - double h,b,a; +const int N=10, m=10; +double s,nx,t; +int j,k; - klo=1; - khi=n; - while (khi-klo > 1) { - k=(khi+klo) >> 1; - if (ary[k-1][0] > x) khi=k; - else klo=k; +s=0.0; +for(k=1;k0.96 && x<1.04 && x!=1.0) /*use laurent series -1 and 0 term and spline */ - { - double l,e; - l= -std::log10(std::fabs(x-1.0)); - if(l>=6.5) std::cerr <<"Argument of zeta too near to the pole, result will be very imprecise!\n"; - if(l>12.0) e=0.0; - else - { - if(x>1.0) splintx(zet3,nzet3,l,&e); - else splintx(zet4,nzet4,l,&e); - } - return(c_euler+e+1.0/(x-1.0)); - } - if(x==std::floor(x) && std::fabs(x)<=(((unsigned int) 1)<<(4*sizeof(int)-1)) ) { int n; @@ -643,13 +270,11 @@ if(x==std::floor(x) && std::fabs(x)<=(((unsigned int) 1)<<(4*sizeof(int)-1)) ) else return(0.0); } if(n==2||n==4||n==6||n==8) /*precalculated few bernoulli numbers */ - return(std::fabs(bernoulli_number(n))*pow(2*M_PI,(double)n)/factorial(n)/2.0); + return(std::fabs(bernoulli_number(n))*pow(2*M_PI,(double)n)/double_factorial(n)/2.0); /*otherwise continue like if general real argument*/ } -if(x<-1.02 || (x>0.58 && x<1.94)) return (2.0*pow(2.0*M_PI,x-1.0)*sin(M_PI*x/2.0)*gamma(1.0-x)*zeta(1.0-x)); +if(x<-1.02) return (2.0*pow(2.0*M_PI,x-1.0)*sin(M_PI*x/2.0)*gamma(1.0-x)*zeta(1.0-x)); - -/* make different approximations in different intervals - some small non-smoothness is introduced */ if(x>=9.1) /* in this region the summation per definition is best and quite exact*/ { double s; @@ -659,8 +284,7 @@ if(x>=9.1) /* in this region the summation per definition is best and quite exac return(1.0+s); } -if(x>=1.94) { splintx(zet2,nzet2,x,&hlp); return(hlp);} -splintx(zet1,nzet1,x,&hlp); return(hlp); +return zeta_euler_maclaurin(x); } //end zeta @@ -672,10 +296,10 @@ double gamma(double x) bool f; double p; -if(x>MAXFACT) laerror("too big argument in gamma"); +if(x>MAXDOUBLEFACT+1) laerror("too big argument in gamma"); f= (std::fabs(x-std::floor(x+0.5))= MAXBERN) /* use zeta function for it */ {int m; m=2*n; - return(zeta((double)m)*(n&1?2.0:-2.0)*factorial(m)/pow(2.0*M_PI,(double)m)); + return(zeta((double)m)*(n&1?2.0:-2.0)*double_factorial(m)/pow(2.0*M_PI,(double)m)); } if(n<=top) return(bern[n]); @@ -751,7 +375,7 @@ for(i=top+1;i<=n;i++) for(i=top+1;i<=n;i++) {int m; m=2*i; - bern[i]=zeta((double)m)*(i&1?2.0:-2.0)*factorial(m)/pow(2.0*M_PI,(double)m); + bern[i]=zeta((double)m)*(i&1?2.0:-2.0)*double_factorial(m)/pow(2.0*M_PI,(double)m); } top=n; return(bern[n]); diff --git a/miscfunc.h b/miscfunc.h index 9b97e09..7c1ea2d 100644 --- a/miscfunc.h +++ b/miscfunc.h @@ -27,7 +27,8 @@ namespace LA { const double c_euler= 0.5772156649015328606065120900824024310422; -extern unsigned long long factorial(const int n); +extern unsigned long long factorial(const int n); //exact, n<=20 +extern double double_factorial(const int n); //n! (not n!!) in double precision, n<=170 extern double factorialln(const int n); extern unsigned long long binom(int n, int k); extern unsigned long long longpow(unsigned long long x, int i);