Skip to content

Instantly share code, notes, and snippets.

@cab1729

cab1729/RiemannSiegel.java

Last active September 10, 2017 16:56
Show Gist options
  • Save cab1729/ee1a25fa08aa03378bae8d9b6cdc86de to your computer and use it in GitHub Desktop.
Save cab1729/ee1a25fa08aa03378bae8d9b6cdc86de to your computer and use it in GitHub Desktop.
Download ZIP
Code Samples - Core Java/Math: Riemann Hypothesis research
Raw
RiemannSiegel.java
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));
}
}
@cab1729
Copy link
Author

cab1729 commented Aug 18, 2017

Replaced lost source file in previous gist version

@cab1729
Copy link
Author

cab1729 commented Aug 18, 2017

For functions.java source see https://gist.github.com/cab1729/1318030

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment