cab1729: math geek with no social life, Code Samples: Riemann-Siegel formula

Java implementation of the Riemann-Siegel Formula used to find non-trivial zeros of the Riemann Zeta Function

import java.io.BufferedReader; import java.io.FileNotFoundException; import java.io.FileReader; import java.io.IOException; public class RiemannSiegel { private static final double twopi = 2.0 * StrictMath.PI; private RiemannSiegel () { // prevent instantiation } /** * Reimann-Siegel Z(t) function with additional parameter * for number of R terms to add * @param t * @param r * @return */ public static double Z (double t, int r) { // t/2pi^1/2 double t2 = StrictMath.sqrt(t/twopi); double N = functions.fabs(t2); double p = t2 - N; int k = 0; double R = 0.0; double tt = theta(t); double sum = 0.0; for (int n = 1; n <= N; n++) { sum += (StrictMath.pow(n, -0.5) * StrictMath.cos(tt - (t * StrictMath.log(n)))); } sum = 2.0 * sum; // add remainder R here double piot = StrictMath.PI/t; while (k <= r) { R = R + C(k,2.0*p-1.0) * StrictMath.pow(2.0 * piot, ((double) k)*0.5); ++k; /* */ } /* */ R = functions.even((int) N-1) * StrictMath.pow(2.0 * StrictMath.PI / t,0.25) * R; /* */ return sum + R; } /** * C terms for Riemann-Siegel * coefficients of remainder terms * @param n * @param z * @return */ public static double C (int n, double z) { if (n==0) return(.38268343236508977173 * StrictMath.pow(z, 0.0) +.43724046807752044936 * StrictMath.pow(z, 2.0) +.13237657548034352332 * StrictMath.pow(z, 4.0) -.01360502604767418865 * StrictMath.pow(z, 6.0) -.01356762197010358089 * StrictMath.pow(z, 8.0) -.00162372532314446528 * StrictMath.pow(z,10.0) +.00029705353733379691 * StrictMath.pow(z,12.0) +.00007943300879521470 * StrictMath.pow(z,14.0) +.00000046556124614505 * StrictMath.pow(z,16.0) -.00000143272516309551 * StrictMath.pow(z,18.0) -.00000010354847112313 * StrictMath.pow(z,20.0) +.00000001235792708386 * StrictMath.pow(z,22.0) +.00000000178810838580 * StrictMath.pow(z,24.0) -.00000000003391414390 * StrictMath.pow(z,26.0) -.00000000001632663390 * StrictMath.pow(z,28.0) -.00000000000037851093 * StrictMath.pow(z,30.0) +.00000000000009327423 * StrictMath.pow(z,32.0) +.00000000000000522184 * StrictMath.pow(z,34.0) -.00000000000000033507 * StrictMath.pow(z,36.0) -.00000000000000003412 * StrictMath.pow(z,38.0) +.00000000000000000058 * StrictMath.pow(z,40.0) +.00000000000000000015 * StrictMath.pow(z,42.0)); else if (n==1) return(-.02682510262837534703 * StrictMath.pow(z, 1.0) +.01378477342635185305 * StrictMath.pow(z, 3.0) +.03849125048223508223 * StrictMath.pow(z, 5.0) +.00987106629906207647 * StrictMath.pow(z, 7.0) -.00331075976085840433 * StrictMath.pow(z, 9.0) -.00146478085779541508 * StrictMath.pow(z,11.0) -.00001320794062487696 * StrictMath.pow(z,13.0) +.00005922748701847141 * StrictMath.pow(z,15.0) +.00000598024258537345 * StrictMath.pow(z,17.0) -.00000096413224561698 * StrictMath.pow(z,19.0) -.00000018334733722714 * StrictMath.pow(z,21.0) +.00000000446708756272 * StrictMath.pow(z,23.0) +.00000000270963508218 * StrictMath.pow(z,25.0) +.00000000007785288654 * StrictMath.pow(z,27.0) -.00000000002343762601 * StrictMath.pow(z,29.0) -.00000000000158301728 * StrictMath.pow(z,31.0) +.00000000000012119942 * StrictMath.pow(z,33.0) +.00000000000001458378 * StrictMath.pow(z,35.0) -.00000000000000028786 * StrictMath.pow(z,37.0) -.00000000000000008663 * StrictMath.pow(z,39.0) -.00000000000000000084 * StrictMath.pow(z,41.0) +.00000000000000000036 * StrictMath.pow(z,43.0) +.00000000000000000001 * StrictMath.pow(z,45.0)); else if (n==2) return(+.00518854283029316849 * StrictMath.pow(z, 0.0) +.00030946583880634746 * StrictMath.pow(z, 2.0) -.01133594107822937338 * StrictMath.pow(z, 4.0) +.00223304574195814477 * StrictMath.pow(z, 6.0) +.00519663740886233021 * StrictMath.pow(z, 8.0) +.00034399144076208337 * StrictMath.pow(z,10.0) -.00059106484274705828 * StrictMath.pow(z,12.0) -.00010229972547935857 * StrictMath.pow(z,14.0) +.00002088839221699276 * StrictMath.pow(z,16.0) +.00000592766549309654 * StrictMath.pow(z,18.0) -.00000016423838362436 * StrictMath.pow(z,20.0) -.00000015161199700941 * StrictMath.pow(z,22.0) -.00000000590780369821 * StrictMath.pow(z,24.0) +.00000000209115148595 * StrictMath.pow(z,26.0) +.00000000017815649583 * StrictMath.pow(z,28.0) -.00000000001616407246 * StrictMath.pow(z,30.0) -.00000000000238069625 * StrictMath.pow(z,32.0) +.00000000000005398265 * StrictMath.pow(z,34.0) +.00000000000001975014 * StrictMath.pow(z,36.0) +.00000000000000023333 * StrictMath.pow(z,38.0) -.00000000000000011188 * StrictMath.pow(z,40.0) -.00000000000000000416 * StrictMath.pow(z,42.0) +.00000000000000000044 * StrictMath.pow(z,44.0) +.00000000000000000003 * StrictMath.pow(z,46.0)); else if (n==3) return(-.00133971609071945690 * StrictMath.pow(z, 1.0) +.00374421513637939370 * StrictMath.pow(z, 3.0) -.00133031789193214681 * StrictMath.pow(z, 5.0) -.00226546607654717871 * StrictMath.pow(z, 7.0) +.00095484999985067304 * StrictMath.pow(z, 9.0) +.00060100384589636039 * StrictMath.pow(z,11.0) -.00010128858286776622 * StrictMath.pow(z,13.0) -.00006865733449299826 * StrictMath.pow(z,15.0) +.00000059853667915386 * StrictMath.pow(z,17.0) +.00000333165985123995 * StrictMath.pow(z,19.0) +.00000021919289102435 * StrictMath.pow(z,21.0) -.00000007890884245681 * StrictMath.pow(z,23.0) -.00000000941468508130 * StrictMath.pow(z,25.0) +.00000000095701162109 * StrictMath.pow(z,27.0) +.00000000018763137453 * StrictMath.pow(z,29.0) -.00000000000443783768 * StrictMath.pow(z,31.0) -.00000000000224267385 * StrictMath.pow(z,33.0) -.00000000000003627687 * StrictMath.pow(z,35.0) +.00000000000001763981 * StrictMath.pow(z,37.0) +.00000000000000079608 * StrictMath.pow(z,39.0) -.00000000000000009420 * StrictMath.pow(z,41.0) -.00000000000000000713 * StrictMath.pow(z,43.0) +.00000000000000000033 * StrictMath.pow(z,45.0) +.00000000000000000004 * StrictMath.pow(z,47.0)); else return(+.00046483389361763382 * StrictMath.pow(z, 0.0) -.00100566073653404708 * StrictMath.pow(z, 2.0) +.00024044856573725793 * StrictMath.pow(z, 4.0) +.00102830861497023219 * StrictMath.pow(z, 6.0) -.00076578610717556442 * StrictMath.pow(z, 8.0) -.00020365286803084818 * StrictMath.pow(z,10.0) +.00023212290491068728 * StrictMath.pow(z,12.0) +.00003260214424386520 * StrictMath.pow(z,14.0) -.00002557906251794953 * StrictMath.pow(z,16.0) -.00000410746443891574 * StrictMath.pow(z,18.0) +.00000117811136403713 * StrictMath.pow(z,20.0) +.00000024456561422485 * StrictMath.pow(z,22.0) -.00000002391582476734 * StrictMath.pow(z,24.0) -.00000000750521420704 * StrictMath.pow(z,26.0) +.00000000013312279416 * StrictMath.pow(z,28.0) +.00000000013440626754 * StrictMath.pow(z,30.0) +.00000000000351377004 * StrictMath.pow(z,32.0) -.00000000000151915445 * StrictMath.pow(z,34.0) -.00000000000008915418 * StrictMath.pow(z,36.0) +.00000000000001119589 * StrictMath.pow(z,38.0) +.00000000000000105160 * StrictMath.pow(z,40.0) -.00000000000000005179 * StrictMath.pow(z,42.0) -.00000000000000000807 * StrictMath.pow(z,44.0) +.00000000000000000011 * StrictMath.pow(z,46.0) +.00000000000000000004 * StrictMath.pow(z,48.0)); } /** * Riemann-Siegel theta function * theta(t) approximation using Stirling series * @param t * @return */ public static double theta (double t) { return (t/2.0 * StrictMath.log(t/2.0/StrictMath.PI) - t/2.0 - StrictMath.PI/8.0 + 1.0/48.0/t + 7.0/5760.0/t/t/t); } }

ref:

H. M. Edwards, Riemann’s Zeta Function, Academic Press, 1974.

E. C. Titchmarsh, The theory of the Riemann Zeta-function, Oxford Science publications, second edition, revised by D. R. Heath-Brown (1986).

ZetaGrid project code