Last active
September 10, 2017 16:56
-
-
Save cab1729/ee1a25fa08aa03378bae8d9b6cdc86de to your computer and use it in GitHub Desktop.
Code Samples - Core Java/Math: Riemann Hypothesis research
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 | |
} | |
/** | |
* Riemann-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); | |
} | |
/** | |
* @param args | |
*/ | |
public static void main(String[] args) { | |
// TODO Auto-generated method stub | |
//double t = Double.parseDouble(args[0]); | |
int r = Integer.parseInt(args[0]); | |
//System.out.println("Z(t) for " + args[0]+ "=" + RiemannSiegel.Z(t,r)); | |
// test values - Haselgrove Table II Z(t) | |
//System.out.println("Z(t) for " + 102.5+ "=" + RiemannSiegel.Z(102.5,r)); | |
//System.out.println("Z(t) for " + 108.5+ "=" + RiemannSiegel.Z(108.5,r)); | |
//System.out.println("Z(t) for " + 114.5+ "=" + RiemannSiegel.Z(114.5,r)); | |
//System.out.println("Z(t) for " + 124.5+ "=" + RiemannSiegel.Z(124.5,r)); | |
//System.out.println("Z(t) for " + 131.0+ "=" + RiemannSiegel.Z(131.0,r)); | |
/* | |
System.out.println("Z(t) for " + 14.134725142+ "=" + RiemannSiegel.Z(14.134725142,r)); | |
System.out.println("Z(t) for " + 21.022039639+ "=" + RiemannSiegel.Z(21.022039639,r)); | |
System.out.println("Z(t) for " + 25.010857580+ "=" + RiemannSiegel.Z(25.010857580,r)); | |
System.out.println("Z(t) for " + 30.424876126+ "=" + RiemannSiegel.Z(30.424876126,r)); | |
System.out.println("Z(t) for " + 32.935061588+ "=" + RiemannSiegel.Z(32.935061588,r)); | |
*/ | |
long start = System.currentTimeMillis(); | |
try { | |
BufferedReader tf = | |
new BufferedReader(new FileReader("zeros3p.txt")); | |
double c = 267653395647.0; | |
String line = null; | |
while ((line = tf.readLine()) != null) { | |
double t = Double.parseDouble(line.trim()); | |
System.out.println("Z(t) for 267653395647.0+" + (t) + "=" + RiemannSiegel.Z(c+t,r)); | |
} | |
} catch (FileNotFoundException e) { | |
// TODO Auto-generated catch block | |
e.printStackTrace(); | |
} catch (IOException e) { | |
// TODO Auto-generated catch block | |
e.printStackTrace(); | |
} | |
long stop = System.currentTimeMillis(); | |
System.out.println("elapsed time in msec: " + (stop-start)); | |
} | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
For functions.java source see https://gist.github.com/cab1729/1318030