/* zeta.cpp, an attempt to compute the Riemann Zeta function Copyright (C) 2000-2008 Robert P. Munafo This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. If you got ries.c from the website www.mrob.com, the GNU General Public License may be retrieved from the following URL: http://www.mrob.com/ries/COPYING.txt To compile, I currently use: gcc -lstdc++ zeta.cpp -lm -o zeta On the Intel Core 2 architecture, a speedup of about 30% is achieved using the flags -O3 and -ffast-math. If given an argument T, it prints a graph of the Zeta function on the critical line (z = 0.5 + i t) from t = T-100 to t = T. If not given an argument, it looks for a file called "cdda.wav" (which should be 44100-Hz, PCM, 16-bit, stereo) and overwrites it with the Zeta function, samples spaced 0.01 apart on the complex plane. Thus, each second of audio will cover 441.0 units on the imaginary axis. This program does not "know" enough about the WAV file format to be able to create a WAV file from scratch -- however I suspect that would be fairly easy to do. I have a description of the AIFF format in "AIFF-1.3.pdf" (in another project) but do not currently have a WAV format description. */ /* Use the following #define to make it evaluate a pre-defined number of samples, for benchmarking purposes. #define BENCH_NSAMP 10000 */ /* REVISION HISTORY 20000302 Initial version 20000303 Post query to sci.math 20000304 Add notes from a reply from sci.math 20050804 Lots of fixes to make it compatible with official ISO implementation of __complex__. Make little-endian-ness of WAV format explicit (so it works on a PowerPC). Run it on a 2-GHz G5, the speed is truly impressive! Increase the number of terms from 1000 to 10000. 20101017 Improve output formatting slightly. 20120102 Add Lanczos Gamma function 20120109 Add Pugh coefficients for Lanczos 20120110 Start using f107_o, add test_f107. 20120111 Convert Lanczos gamma functions to f107. 20120116 u64/s64 need to be 'long long' 20120117 Add gosper_gamma_z_small 20120118 Add print_macsyma_verify, integer_factorial and integer_double_factorial TO DO Try to find a better way to compute the Zeta function. Here are some candidates (compiled on 20100528): 1. The Knopp/Hasse series (formula 4 in "Sondow 1994 Analytic.pdf"): zeta(s) = 1/(1-2^(1-s)) * SIGMAv{n=0..inf} [ 1/(2^(n+1)) SIGMAv{k=0..n} [ -1^k (n k) (k+1)^-s ] ] 2. Algorithm 2 in "Borwein 2000 Efficient.pdf". Select an N, and use: zeta(s) = -1/(d(N)*(1-2^(1-s))) * SIGMAv{k=0..N-1} [ -1^k (d(k)-d(N))/(k+1)^s ] + error d(k) = N * SIGMAv{i=0..k} [ (N+i+1)! * 4^i / ((N-1)!*(2i)!) ] and note that |error| decreases like 1/(3+sqrt(8))^N for Re(s) >= 1/2. If many zeta(s) values are being computed with the same N, the d(k) can be precomputed. 3. Formula 34 in "Borwein 1999 Computational.pdf": zeta(s) = limv{N->inf} 1/(2^(N-s+1)-2^N) * SIGMAv{k=0..2N-1} [ -1^k/(k+1)^s SIGMAv{m=0..k-N} [(N m)-2^N] ] 4. Algorithm 3 in "Borwein 2000 Efficient.pdf". Select an N, and use: zeta(s) = -1/(2^N(1-2^(1-s))) * SIGMAv{j=0..2N-1} [e_j/(j+1)^s] + error e_j = -1^j * SIGMAv{k=0..j-N} [(N k) - 2^N] and note that |error| decreases like 1/8^N for Re(s)>0 (and a somewhat more complex error bound for Re(s)<=0). If many zeta(s) values are being computed with the same N, the e_j can be precomputed. Borwein notes that this method is "even simpler, though not quite as fast" as the above "Algorithm 2". LINKS http://web.viu.ca/pughg/RiemannZeta/RiemannZetaLong.html Nice Java applet and some fairly comprehensible explanations SEE ALSO There is other info about special functions in: .../proj/zeta/ken-takusagawa .../proj/ries/function-sources.txt */ #include #include #include #include /* Bring in the C++ objects for higher-precision maths */ #define F107_APPROX 0 #include "f107_o.cpp"; typedef signed char s8; typedef unsigned char u8; typedef short s16; typedef unsigned short u16; typedef int s32; typedef unsigned int u32; /* For the 64-bit types we must use 'long long' to work in the GCC compiler on 32-bit targets */ typedef long long s64; typedef unsigned long long u64; // Size of wav file header #define WAV_START 44 #define LOUDNESS 1000.0 /* FREQUENCY is the number of units on the imaginary axis that will be heard each second. Imagine that the Zeta function has been carved into a phonograph record and we're playing the record. This constant tells how many units on the real axis the needle moves each second. */ #define FREQUENCY 441.0 /* INTERVAL is the spacing between values of the argument we pass to the Zeta function to produce the audio samples. Note the constant 44100 -- we assume that our output sound file is 44100 samples per second. */ #define INTERVAL (FREQUENCY / 44100.0) /* Function prototypes */ extern "C" { void check_sizes(void); double logaddexp(double a, double b); double cabssqr(__complex__ double z); double cabs(__complex__ double z); __complex__ double clog(__complex__ double z); __complex__ double cexp(__complex__ double z); __complex__ double cpow(double b, __complex__ double e); void cprint(__complex__ double z); void test_f107(void); __complex__ double gamma1(__complex__ double z); f107_o integer_factorial(double z); f107_o pugh_ln_gamma(f107_o z); f107_o lanczos_ln_gamma(f107_o z); f107_o lanczos_gamma(f107_o z); f107_o gosper_gamma_z_small(f107_o zz); f107_o gosper_gamma_z_ge_14(f107_o zz); __complex__ double zetacrit1(double zi, s32 terms); void print_macsyma_verify(const char * s1, f107_o n1, const char * s2, f107_o n2, const char * s3); int main(int, char **); } void check_sizes(void) { int s; s = (int) sizeof(u16); if (s != 2) { printf("sizeof(u16) == %d, expected size was 2\n", s); exit(-1); } s = (int) sizeof(u32); if (s != 4) { printf("sizeof(u32) == %d, expected size was 4\n", s); exit(-1); } s = (int) sizeof(u64); if (s != 8) { printf("sizeof(u64) == %d, expected size was 8\n", s); exit(-1); } } double g_pi = 3.1415926535897932; /* logaddexp returns ln(exp(a) + exp(b)), but works even when both exp's would overflow and/or underflow */ double logaddexp(double a, double b) { double l, s; /* l is the larger, s is the smaller */ /* Given a and b, we want to calculate ln(exp(a)+exp(b)). Define eA=exp(a) and eB = exp(b). Then we want ln(eA+eB). Note that if eB is really small compared to eA, then roundoff reduces this to ln(eA) which is just a. ln(eA+eB) = ln(eA) + ln(1 + eB/eA) = a + ln(1 + exp(b - a)) We start by putting the larger (i.e. higher) of a and b into "l" and the other into "s". */ if (a > b) { l = a; s = b; } else { l = b; s = a; } if (s-l < -50.0) { /* exp(s-l) is less than 10^-20, so log(1.0+exp(s-l)) rounds off to 0 */ return l; } else { return (l + log(1.0 + exp(s-l))); } } double cabssqr(__complex__ double z) { return((__real__ z)*(__real__ z) + (__imag__ z)*(__imag__ z)); } double cabs(__complex__ double z) { return(sqrt(cabssqr(z))); } __complex__ double clog(__complex__ double z) { return(log(cabs(z)) + (atan2(__imag__ z, __real__ z) * 1.0fi)); } /* complex exponential function */ __complex__ double cexp(__complex__ double z) { __complex__ double rv; rv = exp(__real__ z) * (cos(__imag__ z) + (1.0fi * sin(__imag__ z))); return(rv); } /* real base to a complex power */ __complex__ double cpow(double b, __complex__ double e) { __complex__ double rv; //printf("base %10.6lg log(base) %10.6lg\n", b, log(b)); //printf(" exponent %10.6lg + %10.6lg i\n", (__real__ e), (__imag__ e)); rv = cexp(log(b) * e); return(rv); } void cprint(__complex__ double z) { printf("%14.10lg + %14.10lg i", __real__ z, __imag__ z); } /*^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ | | GAMMA FUNCTION | | Following are several routines for computing numerical approzimations of | the Gamma function or its logarithm. For our purposes the leading candidates | are Stirling series, Spouge's approximation and the Lanczos approximation. | A good comparison of these is found in: | | Glendon Ralph Pugh, An Analysis of the Lanczos gamma approximation | (PhD thesis), University of British Columbia, 2004. | PDF at: http://web.viu.ca/pughg/phdThesis/phdThesis.pdf | | In general, if you want to compute the value many times with the same | precision (relative error bound), Lanczos is the best option because it | requires the fewest operations for each Gamma evaluation. However, if you | want to change precisions often, Lanczos is not good because each precision | requires a different set of coefficients and these are hard to compute. | Spouge's approximation is much easier for this purpose. | Finally, the Stirling series offers the advantage of not requiring more | than one or two division operations (one of which is a reciprocal of Z^2, | which is evaluated one and used multiple times while adding the terms | of the series sum). | |___________________________________________________________________________*/ /* - - - - - - - - - - - - - - - - F107 - - - - - - - - - - - - - - - - - The f107 "double-double" type is needed to get enough precision for some of the Lanczos approximations. */ #define F107_K(x) (f107_sscan((char *) x)) /* f107 actually only provides about 32.2 digits, not the 34 you get from IEEE binary128 */ #define QUAD_DIGITS 34 void test_f107(void) { char out[QUAD_DIGITS+10]; char instr[20]; f107_o x, a; int i, j, k, l; /* Test string output */ k = QUAD_DIGITS; for(i=0; i<2; i++) { x = 1234.0; x = x / 1.0e10; x = x / 1.0e10; if (i) { x = 0.0 - x; } for(j=0; j<52; j++) { k = (k + 7) % QUAD_DIGITS; /* Make sure we get errors */ out[0] = 0; f107_spfg(out, k+3, 0, k, x); printf(" %d.%d: %s\n", k+3, k, out); x = 27.0 * x; if (fabs(x) > 1.0e6) { x = 143.0 * x; } } } /* Test string input and output */ for(i=0; i<10; i++) { for(l=0; l<2; l++) { for(j=0; j<6; j++) { strncpy(out, "0000012345", sizeof(out)); strncpy(instr, out+i, sizeof(instr)); /* Shift digits one space to make room for '.' */ for(k=sizeof(instr)-1; k>j; k--) { instr[k] = instr[k-1]; } instr[j] = '.'; if (l == 0) { printf(" %-11s", instr); } else { x = f107_sscan(instr); f107_spfg(out, sizeof(out), 0, 5, x); printf(" %-11s", out); } } printf("\n"); } printf("\n"); } /* Test divide */ for(a=2.0; a<14.0; a += 1.0) { x = 1.0 / a; f107_spfg(out, sizeof(out), 0, 3, a); printf("1/%s = ", out); f107_spfg(out, sizeof(out), 0, 33, x); printf("%s\n", out); } /* Test square root */ for(a=2.0; a<14.0; a += 1.0) { x = sqrt(a); f107_spfg(out, sizeof(out), 0, 3, a); printf("sqrt(%s) = ", out); f107_spfg(out, sizeof(out), 0, 33, x); printf("%s\n", out); } /* Test sin and cos */ for(a=-40.0; a<40.0; a += 2.0) { x = sin(a/10.0); f107_spfg(out, sizeof(out), 0, 3, a/10.0); printf("sin(%s) = ", out); f107_spfg(out, sizeof(out), 0, 33, x); printf("%s", out); x = cos(a/10.0); f107_spfg(out, sizeof(out), 0, 3, a/10.0); printf("; cos(%s) = ", out); f107_spfg(out, sizeof(out), 0, 33, x); printf("%s\n", out); } /* Test exp and log */ for(a=-40.0; a<40.0; a += 2.0) { x = exp(a/10.0); f107_spfg(out, sizeof(out), 0, 3, a/10.0); printf("exp(%s) = ", out); f107_spfg(out, sizeof(out), 0, 33, x); printf("%s", out); x = log(a/10.0); f107_spfg(out, sizeof(out), 0, 3, a/10.0); printf("; log(%s) = ", out); f107_spfg(out, sizeof(out), 0, 33, x); printf("%s\n", out); } for(a=0.0; a<20.0; a += 1.0) { x = tanh(a/10.0); f107_spfg(out, sizeof(out), 0, 3, a/10.0); printf("tanh(%s) = ", out); f107_spfg(out, sizeof(out), 0, 33, x); printf("%s", out); /* Check consistency of tanh to exp by the identity formula: exp(x) = (1+t)/(1-t) where t=tanh(x/2) */ x = tanh(a/20.0); x = (1.0 + x) / (1.0 - x); f107_spfg(out, sizeof(out), 0, 3, a/10.0); printf(" : exp(%s) ~ ", out); f107_spfg(out, sizeof(out), 0, 33, x); printf("%s\n", out); } // 0.693 is slightly less than ln(2) for(a = 69300.0; a < 69313.0; a = a + 1.0) { x = a / 100000.0; f107_spfg(out, sizeof(out), 0, 33, x); printf("exp(%s) = ", out); f107_spfg(out, sizeof(out), 0, 33, exp(x)); printf("%s\n", out); // a = a + 1.0; } } /* - - - - - - - - - - - - - STIRLING GAMMA - - - - - - - - - - - - - - - - - This is an attempt to apply analytic extension to the Stirling series for the Gamma function. I haven't verified its accuracy yet. I am using the Stirling series (generalized to analytic complex functions) and the identity: gamma(z+1) = z gamma(z) gamma(z) = gamma(z+1) / z and the (similar but not identical) identity given on the Weisstein page: ln(gamma(x + iy + 1)) = ln(x^2 + y^2) + i arctan(y/x) + ln(gamma(x + iy)) ln(gamma(x + iy)) = ln(gamma(x + iy + 1)) - ln(x^2 + y^2) - i arctan(y/x) this second identity would be identical to the first, if only they had "ln(sqrt(x^2 + y^2))" instead of "ln(x^2 + y^2)". I suspect it's an error. The whole thing is computed in logs and the answer computed with exp() at the end. Thus, it is readily converted to a log(gamma(z)) function. */ __complex__ double gamma1(__complex__ double z) { __complex__ double scale; __complex__ double l; __complex__ double zp; /* need to test for negative integers and zero */ scale = 0.0i + 0.0; /* On each identity transformation, set its limit to "0.0" to disable, around 5.0-10.0 to enable */ while(cabs(z) < 0.0) { /* strange identity, disabled because it's suspect */ scale = scale + log(cabssqr(z)) + (1.0fi * atan2(__imag__ z, __real__ z)); z = z + 1.0; } while (((__real__ z) <= 0.0) || (cabs(z) < 8.0)) { scale = scale + clog(z); z = z + 1.0; } l = 0.0 - scale; /* The series I use for gamma is actually the Stirling series adapted for computing factorials, that is, x! = gamma(x+1). Rather than change the series back to a real gamma series, I just subtract 1 from z. Works just fine (-: */ z = z - 1.0; l = l + 0.5 * log(2.0 * g_pi); l = l + (z + 0.5) * clog(z); l = l - z; l = l + 1.0 / (12.0 * z); //printf("gamma 1 "); cprint(cexp(l)); printf("\n"); zp = z * z * z; /* z^3 */ l = l - 1.0 / (360.0 * zp); //printf("gamma 2 "); cprint(cexp(l)); printf("\n"); zp = zp * z * z; /* z^5 */ l = l + 1.0 / (1260.0 * zp); //printf("gamma 3 "); cprint(cexp(l)); printf("\n"); zp = zp * z * z; /* z^7 */ l = l - 1.0 / (1680.0 * zp); //printf("gamma 4 "); cprint(cexp(l)); printf("\n"); zp = zp * z * z; /* z^9 */ l = l + 1.0 / (1188.0 * zp); printf("gamma 5 "); cprint(cexp(l)); printf("\n"); return(cexp(l)); } /* - - - - - - - - - - - - - - - LANCZOS GAMMA - - - - - - - - - - - - - - - - Following are several sets of Lanczos coefficients and two evaluation techniques. The Pugh technique is from the Pugh paper cited above. */ f107_o f107_pi = "3.1415926535 8979323846 2643383279 5028841971"; f107_o ln_2_sqrt_eopi = "0.6207822376 3524522234 5518445781 6472122518"; f107_o f107_e = "2.7182818284 5904523536 0287471352 6624977572"; /* Helper functions for exact_gamma. In addition to these we also use f107_intpow and f107_pi */ f107_o integer_factorial(double z) { f107_o x; if (z != floor(z)) { x.c0 = x.c1 = K107_NaN; return x; } if (z < 0.0) { x.c0 = x.c1 = K107_NaN; return x; } if (z > 170.0) { x.c0 = x.c1 = K107_infinity; return x; } x = 1.0; while (z > 1.0) { x = x * z; z = z - 1.0; } return x; } f107_o integer_double_factorial(double z) { f107_o x; if (z != floor(z)) { x.c0 = x.c1 = K107_NaN; return x; } if (z < 0.0) { x.c0 = x.c1 = K107_NaN; return x; } if (z > 300.0) { x.c0 = x.c1 = K107_infinity; return x; } x = 1.0; while (z > 1.0) { x = x * z; z = z - 2.0; } return x; } #if 0 /* This gives about 11.5 digits of accuracy */ /* 4, 0, 4, 3, 7, 16 */ # define USE_PUGH_METHOD 0 # define LG_g "5.0" /* Lanczos parameter "g" */ # define LG_N 6 /* Range of coefficients i=[0..N] */ const char * lct[LG_N+1] = { "1.000000000190015", "76.18009172947146", "-86.50532032941677", "24.01409824083091", "-1.231739572450155", "0.1208650973866179e-2", "-0.5395239384953e-5" }; #elif 0 /* This gives about 15.5 digits accuracy */ /* 4, 9, 13, 3, 25, 21 */ # define USE_PUGH_METHOD 0 # define LG_g "7.0" /* Lanczos parameter "g" */ # define LG_N 8 /* Range of coefficients i=[0..N] */ const char * lct[LG_N+1] = { "0.99999999999980993227684700473478", "676.520368121885098567009190444019", "-1259.13921672240287047156078755283", "771.3234287776530788486528258894", "-176.61502916214059906584551354", "12.507343278686904814458936853", "-0.13857109526572011689554707", "9.984369578019570859563e-6", "1.50563273514931155834e-7" }; #elif 0 /* This gives about 19.5 digits accuracy */ /* 4, 0, 4, 3, 8, 6 */ # define USE_PUGH_METHOD 0 # define LG_g "4.7421875" /* Lanczos parameter "g" */ # define LG_N 14 /* Range of coefficients i=[0..N] */ const char * lct[LG_N+1] = { "0.99999999999999709182", "57.156235665862923517", "-59.597960355475491248", "14.136097974741747174", "-0.49191381609762019978", ".33994649984811888699e-4", ".46523628927048575665e-4", "-.98374475304879564677e-4", ".15808870322491248884e-3", "-.21026444172410488319e-3", ".21743961811521264320e-3", "-.16431810653676389022e-3", ".84418223983852743293e-4", "-.26190838401581408670e-4", ".36899182659531622704e-5" }; #else /* This set of coefficients has the potential for 32 digits of accuracy, but f107 gives us only about 24.5 digits. It uses a slightly different formula. See "Pugh 2004 Analysis.pdf" */ # define USE_PUGH_METHOD 1 # define LG_g "22.618910" /* Lanczos parameter "g" */ # define LG_N 21 /* Range of coefficients i=[0..N] */ const char * lct[LG_N+1] = { "+2.0240434640140357514731512432760e-10", "+1.5333183020199267370932516012553e0", "-1.1640274608858812982567477805332e1", "+4.0053698000222503376927701573076e1", "-8.2667863469173479039227422723581e1", "+1.1414465885256804336106748692495e2", "-1.1135645608449754488425056563075e2", "+7.9037451549298877731413453151252e1", "-4.1415428804507353801947558814560e1", "+1.6094742170165161102085734210327e1", "-4.6223809979028638614212851576524e0", "+9.7030884294357827423006360746167e-1", "-1.4607332380456449418243363858893e-1", "+1.5330325530769204955496334450658e-2", "-1.0773862404547660506042948153734e-3", "+4.7911128916072940196391032755132e-5", "-1.2437781042887028450811158692678e-6", "+1.6751019107496606112103160490729e-8", "-9.7674656970897286097939311684868e-11", "+1.8326577220560509759575892664132e-13", "-6.4508377189118502115673823719605e-17", "+1.3382662604773700632782310392171e-21" }; #endif f107_o pugh_ln_gamma(f107_o z) { f107_o sum; f107_o base; f107_o rv; f107_o g; int i; //char fmt[40]; if (z < 0.5) { /* Use Euler's reflection formula: Gamma(z) = Pi / [Sin[Pi*z] * Gamma[1-z]]; */ return log(f107_pi / sin(f107_pi * z)) - pugh_ln_gamma(1.0 - z); } g = F107_K(LG_g); z = z - 1.0; base = z + g + 0.5; // Base of the Lanczos exponential //f107_spfg(fmt, sizeof(fmt), 0, 33, base); //printf("base: %s\n", fmt); /* We start with the terms that have the smallest coefficients and largest denominator. */ sum = F107_K(lct[0]); //f107_spfg(fmt, sizeof(fmt), 0, 33, sum); //printf("sum: %s\n", fmt); for(i=LG_N; i>=1; i--) { sum += F107_K(lct[i]) / (z + ((f107_o) i)); //f107_spfg(fmt, sizeof(fmt), 0, 33, sum); //printf("sum: %s\n", fmt); } /* Gamma[z] = 2 sqrt(e/pi) * (base/e)^(z+0.5) * sum */ return (ln_2_sqrt_eopi + (log(base)-1.0)*(z+0.5) + log(sum) ); } f107_o lanczos_ln_gamma(f107_o z) { f107_o ln_sqrt_2_pi = "0.9189385332 0467274178 0329736405 6176398613"; f107_o sum; f107_o base; f107_o rv; int i; if (z < 0.5) { /* Use Euler's reflection formula: Gamma(z) = Pi / [Sin[Pi*z] * Gamma[1-z]]; */ return log(f107_pi / sin(f107_pi * z)) - lanczos_ln_gamma(1.0 - z); } z = z - 1.0; base = z + F107_K(LG_g) + 0.5; // Base of the Lanczos exponential sum = 0.0; sum += F107_K(lct[0]); /* We start with the terms that have the smallest coefficients and largest denominator. */ for(i=LG_N; i>=1; i--) { sum += F107_K(lct[i]) / (z + ((f107_o) i)); } // printf("ls2p %7g l(b^e) %7g -b %7g l(s) %7g\n", ln_sqrt_2_pi, // log(base)*(z+0.5), -base, log(sum)); /* Gamma[z] = Sqrt(2*Pi) * sum * base^[z + 0.5] / E^base */ return ln_sqrt_2_pi + log(sum) - base + log(base)*(z+0.5); } f107_o lanczos_gamma(f107_o z) { f107_o lng; if (USE_PUGH_METHOD) { lng = pugh_ln_gamma(z); } else { lng = lanczos_ln_gamma(z); } return(exp(lng)); } /* - - - - - - - - - - - - - - - GOSPER GAMMA - - - - - - - - - - - - - - - - Bill Gosper has produced many formulas related to the Gamma function. This is a special approximation formula he cooked up for me in January 2012. As originally formulated, it is a factorial approximation formula. Given coefficients a1..a5 and b1..b5, the value of z! is approximated by (pseudo-code): ' numer = z^2*(a1*z^6+a2*z^4+a3*z^2+a4+z^8)+a5 ' = z^10+a1*z^8+a2*z^6+a3*z^4+a4*z^2+a5 ` = ((((z^2+a1)*z^2+a2)*z^2+a3)*z^2+a4)*z^2+a5 ' denom = 810*z^5*(z^2*(z^2*(b1*z^4+b2*z^2+b3+z^6)+b4)+b5) ' = 810*(z^15+b1*z^13+b2*z^11+b3*z^9+b4*z^7+b5*z^5 ` = 810*(((((z^2+b1)*z^2+b2)*z^2+b3)*z^2+b4)*z^2+b5)*z^5 ' s1 = sqrt(e^(numer/denom - 2*z)/z) ' sinh = (e^(1/z) - e^(-(1/z)))/2 ' z! ~= sqrt(2*pi)*z^(3*z/2+1)*sinh^(z/2)*s1 This gives about 33 or 34 digits of accuracy for arguments z>=12. In order to cover the full range we turn Gosper's formula into an lngamma formula: ' to compute gamma(zz): z = zz + 1 numer = ... ' same as above denom = ... ' same as above lns1 = numer/2*denom - z - ln(z)/2 ez = e^(1/z) sinh = (ez - 1/ez)/2 lngamma = ln(sqrt(2*pi)) + ln(z)*(3*z/2+1) + ln(sinh)*z/2 + lns1 ' then if desired gamma = exp(lngamma) */ f107_o gosper_gamma_z_ge_14(f107_o zz) { /* %%% In release version, these need to be set up once per runtime */ f107_o a1 = "35.84203843993396256574668265810393153"; f107_o a2 = "388.5030350779587569420887020038442966"; f107_o a3 = "1460.896036726577552102819199055099138"; f107_o a4 = "1523.289051544707117362888206714560808"; f107_o a5 = "-5.609999545496366734012531925182500548"; f107_o b1 = "36.78489558279110542277098460981101565"; f107_o b2 = "421.8240318655427519890428941274746232"; f107_o b3 = "1811.624561248534982160437726096585511"; f107_o b4 = "2760.786439188096048960781896174485476"; f107_o b5 = "1106.322477561403375269124002624801732"; f107_o z, zs, scale; f107_o numer, denom; f107_o gamma; z = zz - 1.0; char fmt[40]; /* Reflection formula for negative arguments will go here */ scale = 1.0; /* Use the induction formula z! = (z+1)! / (z+1) */ while (z < 14.0) { z = z + 1.0; scale = scale * z; } f107_o log_z = log(z); //f107_spfg(fmt, sizeof(fmt), 0, 33, log_z); printf("log_z=%s\n", fmt); zs = f107_o::f107_native_square(z); //f107_spfg(fmt, sizeof(fmt), 0, 33, zs); printf("zs=%s\n", fmt); numer = ((((zs + a1)*zs + a2)*zs + a3)*zs + a4)*zs + a5; //f107_spfg(fmt, sizeof(fmt), 0, 33, numer); printf("numer=%s\n", fmt); denom = 810.0 * (((((zs + b1)*zs + b2)*zs + b3)*zs + b4)*zs + b5)*zs*zs*z; //f107_spfg(fmt, sizeof(fmt), 0, 33, denom); printf("denom=%s\n", fmt); f107_o sinh2 = sinh(1.0/z); #if 0 // Compute Log[Gamma[z]] (for larger range of arguments) f107_o ln_sqrt_2_pi = "0.9189385332 0467274178 0329736405 6176398613"; f107_o lns1 = (z - 0.5*numer/denom) + 0.5 * log_z; f107_o lngamma = ln_sqrt_2_pi + log_z * (1.5*z + 1.0) + log(sinh2) * z * 0.5 - lns1 - log(scale); /* Correction for negative argument will go here; NaN if lngamma would be complex */ gamma = exp(lngamma); #else // Compute Gamma[z] f107_o sqrt_2_pi = "2.5066282746 3100050241 5765284811 0452530069"; gamma = sqrt_2_pi * exp(log_z*(1.5*z + 0.5) + 0.5*numer/denom) * exp(log(sinh2)*z*0.5 - z) / scale; /* Correction for negative argument goes here */ return gamma; #endif } f107_o gosper_gamma_z_small(f107_o zz) { /* Parameters or zeta-based approximation, shown in farquad18.png For -1/4 Gamma[z-1] = Gamma[z]/(z-1) */ if (zz < 0.75) { scale = scale / zz; zz = zz + 1.0; } /* Gamma[z+1] = z * Gamma[z] */ while (zz > 1.75) { // printf("induct: will compute gamma[%g] as gamma[%g] * %g\n", // zz.c0, zz.c0 - 1.0, zz.c0 - 1.0); zz = zz - 1.0; scale = scale * zz; } /* Gosper's formula is for factorial, so we will compute x! = Gamma[zz]. Thus x = zz-1 */ f107_o x = zz - 1.0; // printf("Will compute gamma[%g] as (%g)!\n", zz.c0, x.c0); // printf("then scale by %g and recip? %d\n", scale.c0, recip); f107_o sum; f107_o gamma; f107_o lngamma; f107_o xpi = x * f107_pi; int k; if (x < 0.25) { // printf("Z sum for x == %g\n", x.c0); /* Use zeta-based approximation, "farquad18.png" */ f107_o px = x; // Power of x during sum */ f107_o xs = f107_o::f107_native_square(px); f107_o term; sum = 0.0; for(k=1; k<=18; k++) { term = (px * Z[k]) / ((f107_o) ((2.0 * (double) k) - 1.0)); // printf("Z: k == %d, term == %g", k, term.c0); sum = sum + term; px = px * xs; // printf(", sum == %g\n", sum.c0); } //printf("xpi = %g, sin(xpi) = %g, log(sin(xpi)/xpi) = %g\n", // xpi.c0, sin(xpi.c0), log(sin(xpi.c0) / (xpi.c0))); if (x == 0.0) { /* sin(pi x) / pi x == 1.0, so we have log(1.0) - sum */ lngamma = 0.0 - sum; } else { lngamma = f107_o::f107_native_scale(log(sin(xpi) / xpi), -0.5) - sum; } // printf("log((%g)!) = %g\n", x.c0, lngamma.c0); } else { // printf("P sum for x == %g\n", x.c0); /* Use psi-based approximation, "farquad2.png" */ f107_o px = x - 0.5; // Power of (x-1/2) during sum */ f107_o xm2s = f107_o::f107_native_square(px); f107_o fact; f107_o term; sum = 0.0; fact = 1.0; for(k=0; k<=14; k++) { fact = fact * ((f107_o) ((2.0 * (double) k) + 1.0)); term = P[k]; // printf("P: k == %d, P[k] == %g", k, term.c0); term = (px * term) / fact; // printf(", fact == %g, term == %g", fact.c0, term.c0); sum = sum + term; px = px * xm2s; // printf(", sum == %g\n", sum.c0); fact = fact * ((f107_o) ((2.0 * (double) k) + 2.0)); } xm2s = 1.0 - x; xm2s = f107_pi * xm2s * x / sin(xpi); lngamma = sum + f107_o::f107_native_scale(log(xm2s), 0.5); // printf("log((%g)!) = %g\n", x.c0, lngamma.c0); } gamma = exp(lngamma); if (recip) { gamma = 1.0 / (scale * gamma); } else { gamma = scale * gamma; } return gamma; } /*^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ | | COMPLEX ZETA FUNCTION | |___________________________________________________________________________*/ /* Here are some formulas to try, but some are probably only valid for real values, and most probably won't be useful: zeta(z) = 1/(1-2^(-z)) PRODUCT(p prime) 1/(1-p^(-z)) don't know if "p prime" includes p=2 -- if not, the first term could be moved inside the product zeta(1-z) = (2/(2 pi)^z) cos(pi z/2) gamma(z) zeta(z) zeta(1/2 - it) = (2/(2 pi)^(1/2 + it)) cos(pi (1/2 + it)/2) gamma(1/2 + it) zeta(1/2 + it) The "main term" of the Riemann-Seigel Z function: Z(t) ~= 2 SIGMA[k=1..v(t)] cos(theta(t) - t ln k) / sqrt(k) theta(t) = arg gamma(1/4 + it/2) - (t ln pi)/2 v(t) = floor(sqrt(t/(2 pi))) */ /* This is a possibly-futile attempt to compute the Riemann-Siegel z(t) function, that is, |zeta(1/2 + i t)| For complex s with real(s) > 0, we have: Z(s) = 1/(1-2^(1-s)) SIGMA[n=1..inf] -1^(n-1) n^(-s) = 1/(1-2^(1-s)) { SIGMA[odd n>0] n^(-s) - SIGMA[even n>0] n^(-s) } The terms converge at a rate of 1/sqrt(n), which is very slowly. They're smaller if you add consecutive terms, but it still converges just as slowly. */ __complex__ double zetacrit1(double zi, s32 terms) { s32 i; __complex__ double z; __complex__ double m; /* the first multiplier, 1/(1-2^(1-s)) */ __complex__ double t; /* one term in the sum */ double n; __complex__ double sum; __complex__ double t1; __complex__ double rv; /* Compute the value of z on the critical line */ z = 0.5 + (zi * 1.0fi); /* compute the first multiplier */ m = 1.0 / (1.0 - cpow(2.0, (1.0 - z))); if (0) { printf("m = %10.6lg + %10.6lg i\n", (__real__ m), (__imag__ m)); } /* Compute first term */ n = 1.0; t = cpow(n, 0.0 - z); sum = t; /* Add some terms */ i = 0; while (i < terms) { /* a negative term */ n += 1.0; t1 = cpow(n, 0.0 - z); /* and a positive term */ n += 1.0; t = cpow(n, 0.0 - z); t = t - t1; sum += t; if(0) { printf("n=%7lg 2 terms = %10.6lg + %10.6lg i", n, (__real__ t), (__imag__ t)); printf(" sum = %10.6lg + %10.6lg i\n", (__real__ sum), (__imag__ sum)); } i++; } // printf(" sum = %10.6lg + %10.6lg i\n", (__real__ sum), (__imag__ sum)); rv = m * sum; // printf(" rv = %10.6lg + %10.6lg i\n", (__real__ rv), (__imag__ rv)); return(rv); } /* From: Raymond Manzoni Date: 03/04/2000 sci.math Hi, Thank you for the take-off sound, For which country? I don't know! Here are some of the papers that could interest you: A nice historical introduction may be found in "Computational Number Theory at CWI in 1970-1994" but the paper (Riele-Lune.ps) doesn't seem available online any longer. At CECM you'll find : "Computational strategies for the Riemann zeta function" at http://www.cecm.sfu.ca/preprints/1998pp.html#98:118 consult too Andrew Odlyzko's "Papers on Zeros of the Riemann Zeta Function and Related Topics" at http://www.research.att.com/~amo/doc/zeta.html (for example "Fast algorithms for multiple evaluations of the Riemann zeta function") and Andrew's "Tables of zeros of the Riemann zeta function" at http://www.research.att.com/~amo/zeta_tables/index.html For more papers concerning Riemann zeta look at : http://www.archetipo.web66.com/zeta/briefintro.htm Or forget all that and use simply the quick pari/gp for evaluation of zeta in the critical strip (available with C source at http://www.parigp-home.de/ ) ftp://megrez.math.u-bordeaux.fr/pub/pari/unix/ Hope it helped a little, Raymond Manzoni I added too a resume of my earlier attempts using Euler Mac-Laurin : I used the following formula for that purpose: zeta(x) ~ sum(1/k^x, k=1..N) + 1/((x-1)*N^(x-1)) -1/(2*N^x) +x/(12*N^(x+1)) -x*(x+1)*(x+2)/(720*N^(x+3)) +x*(x+1)*(x+2)*(x+3)*(x+4)/(30240*N^(x+5)) +... The greater the N the better the result (but even N=10 is good) The accuracy is rather fine when |Im(x)| < 2*PI*N and becomes very bad after that ! You may find more terms using the more general Euler Mac-Laurin formula : This allows to explore better the negative values of x (in fact you will get the exact values for x=-1,-2, until -2 times the number of Bernoulli terms!). With more terms the precision in the imaginary plane is better under 2*PI*N and worse over it! f(x) ~ sum(f(k), k=1..N) + int(f(u),u=N..oo) -f(N)/2 -B1/2!*f'(N) +B2/4!*f'''(N) -... where Bn means the nth Bernoulli number This formula was very fruitful for me even on my little TI 92. (N=100 was enough) */ /* Print a string which can be put into maxima to verify a calculation. The output looks like: gamma(-1.5b0) - 2.3632718012073547030642233111215b0; To get the output in this example, call the function like this: print_macsyma_verify("gamma(", z, ")", x); The "b0" after each number, the " - " after "gamma(-1.5b0)" and the ";" at the end are generated automatically. */ void print_macsyma_verify(const char * s1, f107_o n1, const char * s2, f107_o n2) { char tmp[40]; /* Format the argument */ f107_spfg(tmp, sizeof(tmp), 0, 4, n1); if (f107_isinf(n2)) { printf("/* %s%s%s is infinite */\n", s1, tmp, s2); } else if (f107_isnan(n2)) { printf("/* %s%s%s is a NaN */\n", s1, tmp, s2); } else { printf("%s", s1); printf("%s", tmp); printf("%s%s%s", "b0", s2, " / "); f107_spfg(tmp, sizeof(tmp), 0, 34, n2); printf("%s", tmp); printf("%s%s", "b0 - 1.0b0", ";\n"); } } int main(int nargs, char ** argv) { char ** av; int args; char * arg; double zi; s16 got_zi; /* got z? */ check_sizes(); got_zi = 0; if (nargs > 1) { av = argv; args = nargs; av++; args--; /* skip our program name/path */ while(args) { arg = *av; av++; args--; /* Test for number */ if (((arg[0] >= '0') && (arg[0] <= '9')) || (arg[0] == '-') || (arg[0] == '.')) { sscanf(arg, "%lf", &zi); got_zi = 1; } } } if (0) { printf("zetacrit1(0.5 + b i) = %14.10lg\n", cabs(zetacrit1(zi, 100))); printf(" 1000 %14.10lg\n", cabs(zetacrit1(zi, 1000))); printf(" 10000 %14.10lg\n", cabs(zetacrit1(zi, 10000))); printf(" 100000 %14.10lg\n", cabs(zetacrit1(zi, 100000))); printf(" 1000000 %14.10lg\n", cabs(zetacrit1(zi, 1000000))); exit(0); } if(0) { gamma1(0.25 + 0.0i); gamma1(3.0 + 0.0i); /* expect 2.0 */ gamma1(5.0 + 0.0i); /* expect 24.0 */ gamma1(0.25 + (1.0fi * zi)); gamma1(0.25 + (1.0fi * (zi + 1.0))); gamma1(0.25 + (1.0fi * (zi + 2.0))); exit(0); } if (1) { /* Test the Lanczos gamma function. The correct values would be: x Gamma(x) x Gamma(x) x Gamma(x) 0.1 9.5135076986687 1.1 0.9513507698668 2.5 1.3293403881791 0.2 4.5908437119988 1.2 0.9181687423997 3 2 0.3 2.9915689876876 1.3 0.8974706963062 3.5 3.3233509704478 0.4 2.2181595437577 1.4 0.887263817503 4 6 0.5 1.7724538509055 1.5 0.8862269254527 4.5 11.631728396567 0.6 1.4891922488128 1.6 0.8935153492876 5 24 0.7 1.2980553326475 1.7 0.9086387328532 5.5 52.342777784554 0.8 1.1642297137253 1.8 0.9313837709802 6 120 0.9 1.0686287021193 1.9 0.9617658319073 6.5 287.88527781504 1 1 2 1 7 720 */ //gosper_gamma_z_ge_14(26.0); f107_stats(); exit(0); f107_o x10, x, k0_1, expect, gam; // for(i=1; i<=70; i++) { expect = 1.0; k0_1 = "0.1"; for (x10=-20.0; x10 <= 270.0;) { char fmt[40]; x = x10 / 10.0; if (x < 10.0) { gam = gosper_gamma_z_small(x); } else { gam = gosper_gamma_z_ge_14(x); } // gam = lanczos_gamma(x); #if 0 f107_spfg(fmt, sizeof(fmt), 0, 3, x); printf(" gamma(%s)", fmt); f107_spfg(fmt, sizeof(fmt), 0, 32, gam); printf(" == %s", fmt); if ((x > 1.0) && (x == floor(x))) { f107_o err = (gam / expect) - 1.0; printf(" err = %10.5g", err.c0); expect = expect * x; } printf("\n"); #else print_macsyma_verify("gamma(", x, ")", gam); #endif x10 = (x10 > 19.5) ? (x10 + 5.0) : (x10 + 1.0); } } if(1) { test_f107(); exit(0); } if(got_zi) { #define TRANSFORM(x) ((int) (((x) * 10.0) + 40.0)) __complex__ double z; double t, z1, z2r, z2i; s32 terms; terms = 500; t = zi - 100.0; if (t < 0.0) { t = 0.0; } for(; t lim) { lim = i2i; } if (i2r > lim) { lim = i2r; } if (lim < 40) { lim = 40; } if(lim > 80) { lim = 80; } for(i=0; i<=lim; i++) { if (i == i2i) { printf("i"); } else if (i == i2r) { printf("r"); } else if (i == i1) { printf("."); } else if (i == 40) { printf("|"); } else { printf(" "); } } printf("\n"); } } else { __complex__ double z; double t, z1r, z1i; s32 nsamp, ctr, cs; s16 samp[2]; s8 dbuf[4]; FILE * wav; /* Measure the size of the file, which tells us how many samples we will want to compute */ wav = fopen("cdda.wav", "r+"); fseek(wav, 0, SEEK_END); nsamp = ftell(wav); nsamp -= WAV_START; nsamp /= 4; nsamp--; #ifdef BENCH_NSAMP if (nsamp > BENCH_NSAMP) { nsamp = BENCH_NSAMP; } #endif printf("WAV file length: %d samples\n", nsamp); printf("Evaluating Zeta(1/2 + i t) for t in (0.0 .. %g)\n", ((double) nsamp) * INTERVAL); printf("(each dot is 441 samples = 0.01 sec.)\n"); fseek(wav, WAV_START, SEEK_SET); printf("t=%10.6lg", 0.0); fflush(stdout); t = 0.0; ctr = 0; cs = 0; while(nsamp > 0) { z = zetacrit1(t, 10000); z1r = __real__ z; z1i = __imag__ z; samp[0] = ((int) (z1i * LOUDNESS)); samp[1] = ((int) (z1r * LOUDNESS)); dbuf[0] = samp[0] & 0xff; dbuf[1] = samp[0] >> 8; dbuf[2] = samp[1] & 0xff; dbuf[3] = samp[1] >> 8; fwrite((void *) dbuf, 2, 2, wav); ctr++; if (ctr >= 441) { printf("."); fflush(stdout); ctr = 0; cs++; if (cs >= 66) { printf("\rt=%10.6lg", t); fflush(stdout); cs = 0; } } nsamp--; t += INTERVAL; } fclose(wav); printf("\nDone!\n"); } return(0); }