Edit: Expected output is posted at the end.
I wrote a program that calculates zeroes for the Riemann Zeta Function by the Riemann-Siegel formula. I am interesting in adjusting one of the methods in the program to observe what is known as Lehmer's Phenomenon.
Inside the program, I am looking to adjust the end of the findRoots() method.
The full program references
/**************************************************************************
**
** Riemann-Siegel Formula for roots of Zeta(s) on critical line.
**
**************************************************************************
** Axion004
** 07/31/2015
**
** This program finds the roots of Zeta(s) using the well known Riemann-
** Siegel formula. The Riemann–Siegel theta function is approximated
** using Stirling's approximation. It also uses an interpolation method to
** locate zeroes. The coefficients for R(t) are handled by the Taylor
** Series approximation originally listed by Haselgrove in 1960. It is
** necessary to use these coefficients in order to increase computational
** speed.
**************************************************************************/
import java.util.Iterator;
import java.util.LinkedHashSet;
//These two imports are from the Apache Commons Math library
import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.analysis.solvers.BracketingNthOrderBrentSolver;
public class RiemannSiegelTwo{
public static void main(String[] args){
SiegelMain();
}
// Main method
public static void SiegelMain() {
System.out.println("Zeroes inside the critical line for " +
"Zeta(1/2 + it). The t values are referenced below.");
System.out.println();
findRoots();
}
/**
* The sign of a calculated double value.
* @param x - the double value.
* @return the sign in -1, 1, or 0 format.
*/
private static int sign(double x) {
if (x < 0.0)
return -1;
else if (x > 0.0)
return 1;
else
return 0;
}
/**
* Finds the roots of a specified function through the
* BracketingNthOrderBrentSolver from the Apache Commons Math library.
* See http://commons.apache.org/proper/commons-math/apidocs/org/
* apache/commons/math3/analysis/solvers/BracketingNthOrderBrentSolver
* .html
* The zeroes inside the interval of 0 < t < 10000 are printed from
* a TreeSet.
*/
public static void findRoots() {
BracketingNthOrderBrentSolver f = new
BracketingNthOrderBrentSolver();
UnivariateFunction func = (double x) -> RiemennZ(x, 4);
LinkedHashSet<Double> set = new LinkedHashSet<>();
double i = 1.0;
while (i < 1000) {
i+= 0.1;
if(sign(func.value(i)) != sign(func.value(i+0.1))) {
set.add(f.solve(1000, func, i, i+0.1));
}
}
Iterator<Double> itr = set.iterator();
while(itr.hasNext()) {
System.out.println(itr.next());
}
}
/**
* Riemann-Siegel theta function using the approximation by the
* Stirling series.
* @param t - the value of t inside the Z(t) function.
* @return Stirling's approximation for theta(t).
*/
public static double theta (double t) {
return (t/2.0 * Math.log(t/(2.0*Math.PI)) - t/2.0 - Math.PI/8.0
+ 1.0/(48.0*Math.pow(t, 1)) + 7.0/(5760*Math.pow(t, 3)));
}
/**
* Computes Math.Floor of the absolute value term passed in as t.
* @param t - the value of t inside the Z(t) function.
* @return Math.floor of the absolute value of t.
*/
public static double fAbs(double t) {
return Math.floor(Math.abs(t));
}
/**
* Riemann-Siegel Z(t) function implemented per the Riemenn Siegel
* formula. See http://mathworld.wolfram.com/Riemann-SiegelFormula.html
* for details
* @param t - the value of t inside the Z(t) function.
* @param r - referenced for calculating the remainder terms by the
* Taylor series approximations.
* @return the approximate value of Z(t) through the Riemann-Siegel
* formula
*/
public static double RiemennZ(double t, int r) {
double twopi = Math.PI * 2.0;
double val = Math.sqrt(t/twopi);
double n = fAbs(val);
double sum = 0.0;
for (int i = 1; i <= n; i++) {
sum += (Math.cos(theta(t) - t * Math.log(i))) / Math.sqrt(i);
}
sum = 2.0 * sum;
double remainder;
double frac = val - n;
int k = 0;
double R = 0.0;
// Necessary to individually calculate each remainder term by using
// Taylor Series co-efficients. These coefficients are defined below.
while (k <= r) {
R = R + C(k, 2.0*frac-1.0) * Math.pow(t / twopi,
((double) k) * -0.5);
k++;
}
remainder = Math.pow(-1, (int)n-1) * Math.pow(t / twopi, -0.25) * R;
return sum + remainder;
}
/**
* C terms for the Riemann-Siegel formula. See
* https://web.viu.ca/pughg/thesis.d/masters.thesis.pdf for details.
* Calculates the Taylor Series coefficients for C0, C1, C2, C3,
* and C4.
* @param n - the number of coefficient terms to use.
* @param z - referenced per the Taylor series calculations.
* @return the Taylor series approximation of the remainder terms.
*/
public static double C (int n, double z) {
if (n==0)
return(.38268343236508977173 * Math.pow(z, 0.0)
+.43724046807752044936 * Math.pow(z, 2.0)
+.13237657548034352332 * Math.pow(z, 4.0)
-.01360502604767418865 * Math.pow(z, 6.0)
-.01356762197010358089 * Math.pow(z, 8.0)
-.00162372532314446528 * Math.pow(z,10.0)
+.00029705353733379691 * Math.pow(z,12.0)
+.00007943300879521470 * Math.pow(z,14.0)
+.00000046556124614505 * Math.pow(z,16.0)
-.00000143272516309551 * Math.pow(z,18.0)
-.00000010354847112313 * Math.pow(z,20.0)
+.00000001235792708386 * Math.pow(z,22.0)
+.00000000178810838580 * Math.pow(z,24.0)
-.00000000003391414390 * Math.pow(z,26.0)
-.00000000001632663390 * Math.pow(z,28.0)
-.00000000000037851093 * Math.pow(z,30.0)
+.00000000000009327423 * Math.pow(z,32.0)
+.00000000000000522184 * Math.pow(z,34.0)
-.00000000000000033507 * Math.pow(z,36.0)
-.00000000000000003412 * Math.pow(z,38.0)
+.00000000000000000058 * Math.pow(z,40.0)
+.00000000000000000015 * Math.pow(z,42.0));
else if (n==1)
return(-.02682510262837534703 * Math.pow(z, 1.0)
+.01378477342635185305 * Math.pow(z, 3.0)
+.03849125048223508223 * Math.pow(z, 5.0)
+.00987106629906207647 * Math.pow(z, 7.0)
-.00331075976085840433 * Math.pow(z, 9.0)
-.00146478085779541508 * Math.pow(z,11.0)
-.00001320794062487696 * Math.pow(z,13.0)
+.00005922748701847141 * Math.pow(z,15.0)
+.00000598024258537345 * Math.pow(z,17.0)
-.00000096413224561698 * Math.pow(z,19.0)
-.00000018334733722714 * Math.pow(z,21.0)
+.00000000446708756272 * Math.pow(z,23.0)
+.00000000270963508218 * Math.pow(z,25.0)
+.00000000007785288654 * Math.pow(z,27.0)
-.00000000002343762601 * Math.pow(z,29.0)
-.00000000000158301728 * Math.pow(z,31.0)
+.00000000000012119942 * Math.pow(z,33.0)
+.00000000000001458378 * Math.pow(z,35.0)
-.00000000000000028786 * Math.pow(z,37.0)
-.00000000000000008663 * Math.pow(z,39.0)
-.00000000000000000084 * Math.pow(z,41.0)
+.00000000000000000036 * Math.pow(z,43.0)
+.00000000000000000001 * Math.pow(z,45.0));
else if (n==2)
return(+.00518854283029316849 * Math.pow(z, 0.0)
+.00030946583880634746 * Math.pow(z, 2.0)
-.01133594107822937338 * Math.pow(z, 4.0)
+.00223304574195814477 * Math.pow(z, 6.0)
+.00519663740886233021 * Math.pow(z, 8.0)
+.00034399144076208337 * Math.pow(z,10.0)
-.00059106484274705828 * Math.pow(z,12.0)
-.00010229972547935857 * Math.pow(z,14.0)
+.00002088839221699276 * Math.pow(z,16.0)
+.00000592766549309654 * Math.pow(z,18.0)
-.00000016423838362436 * Math.pow(z,20.0)
-.00000015161199700941 * Math.pow(z,22.0)
-.00000000590780369821 * Math.pow(z,24.0)
+.00000000209115148595 * Math.pow(z,26.0)
+.00000000017815649583 * Math.pow(z,28.0)
-.00000000001616407246 * Math.pow(z,30.0)
-.00000000000238069625 * Math.pow(z,32.0)
+.00000000000005398265 * Math.pow(z,34.0)
+.00000000000001975014 * Math.pow(z,36.0)
+.00000000000000023333 * Math.pow(z,38.0)
-.00000000000000011188 * Math.pow(z,40.0)
-.00000000000000000416 * Math.pow(z,42.0)
+.00000000000000000044 * Math.pow(z,44.0)
+.00000000000000000003 * Math.pow(z,46.0));
else if (n==3)
return(-.00133971609071945690 * Math.pow(z, 1.0)
+.00374421513637939370 * Math.pow(z, 3.0)
-.00133031789193214681 * Math.pow(z, 5.0)
-.00226546607654717871 * Math.pow(z, 7.0)
+.00095484999985067304 * Math.pow(z, 9.0)
+.00060100384589636039 * Math.pow(z,11.0)
-.00010128858286776622 * Math.pow(z,13.0)
-.00006865733449299826 * Math.pow(z,15.0)
+.00000059853667915386 * Math.pow(z,17.0)
+.00000333165985123995 * Math.pow(z,19.0)
+.00000021919289102435 * Math.pow(z,21.0)
-.00000007890884245681 * Math.pow(z,23.0)
-.00000000941468508130 * Math.pow(z,25.0)
+.00000000095701162109 * Math.pow(z,27.0)
+.00000000018763137453 * Math.pow(z,29.0)
-.00000000000443783768 * Math.pow(z,31.0)
-.00000000000224267385 * Math.pow(z,33.0)
-.00000000000003627687 * Math.pow(z,35.0)
+.00000000000001763981 * Math.pow(z,37.0)
+.00000000000000079608 * Math.pow(z,39.0)
-.00000000000000009420 * Math.pow(z,41.0)
-.00000000000000000713 * Math.pow(z,43.0)
+.00000000000000000033 * Math.pow(z,45.0)
+.00000000000000000004 * Math.pow(z,47.0));
else
return(+.00046483389361763382 * Math.pow(z, 0.0)
-.00100566073653404708 * Math.pow(z, 2.0)
+.00024044856573725793 * Math.pow(z, 4.0)
+.00102830861497023219 * Math.pow(z, 6.0)
-.00076578610717556442 * Math.pow(z, 8.0)
-.00020365286803084818 * Math.pow(z,10.0)
+.00023212290491068728 * Math.pow(z,12.0)
+.00003260214424386520 * Math.pow(z,14.0)
-.00002557906251794953 * Math.pow(z,16.0)
-.00000410746443891574 * Math.pow(z,18.0)
+.00000117811136403713 * Math.pow(z,20.0)
+.00000024456561422485 * Math.pow(z,22.0)
-.00000002391582476734 * Math.pow(z,24.0)
-.00000000750521420704 * Math.pow(z,26.0)
+.00000000013312279416 * Math.pow(z,28.0)
+.00000000013440626754 * Math.pow(z,30.0)
+.00000000000351377004 * Math.pow(z,32.0)
-.00000000000151915445 * Math.pow(z,34.0)
-.00000000000008915418 * Math.pow(z,36.0)
+.00000000000001119589 * Math.pow(z,38.0)
+.00000000000000105160 * Math.pow(z,40.0)
-.00000000000000005179 * Math.pow(z,42.0)
-.00000000000000000807 * Math.pow(z,44.0)
+.00000000000000000011 * Math.pow(z,46.0)
+.00000000000000000004 * Math.pow(z,48.0));
}
}
The zeroes which are printed include
14.134728277620736
21.022037047686375
25.010858848857314
30.42487545102874
32.93506194020564
37.5861781599404
40.91871891395697
43.32707339344611
48.0051507135618
49.773832659813834
52.97032138572673
56.446247773482625
59.34704389126199
60.83177861205533
65.11254400260545
67.07981057423619
69.54640168845049
72.067157682421
75.70469070082693
77.14484006152179
79.33737502795745
82.91038084149456
84.73549299832806
87.42527458836074
88.809111229007
92.49189925490495
94.65134407077679
95.87063420239032
98.83119422921574
101.31785099451231
103.72553805583749
105.44662303529545
107.1686111943489
111.02953552739187
111.87465919250342
114.32022090805314
116.22668032564269
118.79078286343636
121.37012500478542
122.9468292915669
124.25681855486741
127.51668388015422
129.5787041982112
131.08768853318705
133.49773719981772
134.7565097562393
138.1160420512064
139.73620895846443
141.1237073980395
143.11184581101554
146.00098248274014
147.4227653471707
150.05352041306668
150.92525762000432
153.02469380843965
156.11290929776234
157.5975918115767
158.8499881762274
161.18896413483722
163.0307096895234
165.5370691851963
167.18443998174794
169.09451541024382
169.91197648339372
173.41153651776858
174.7541915257848
176.441434295885
178.37740777747973
179.91648401935672
182.2070784847937
184.87446784765135
185.59878367831283
187.22892258339354
189.41615865594937
192.02665636108685
193.07972660281158
195.26539668016733
196.87648183995336
198.0153096770279
201.26475194277046
202.49359451564763
204.18967180130994
205.39469720336314
207.90625888685815
209.57650971791733
211.69086259398296
213.34791936192354
214.5470447814087
216.16953850912452
219.06759634795878
220.71491884220737
221.43070555208672
224.0070002561315
224.98332466649444
... and much more. To observe Lehmer's Phenomenon, I am looking to adjust the following
Iterator<Double> itr = set.iterator();
while(itr.hasNext()) {
System.out.println(itr.next());
}
To something similar to
// Not workable code
Iterator<Double> itr = set.iterator();
while(itr.hasNext()) {
double EPSILON = .001;
if(itr.nextVal - itr.currentVal < EPSILON )
System.out.println(itr.next());
}
So that zeroes are only printed if the two values are close together. Can I do this with the Iterator? Do I need to use a ListIterator instead?
I tried an enhanced for loop and ran into the same issue. Is there a better or more efficient way to compare the values that are printed from the LinkedHashSet set?
(Note: I originally created a TreeSet instead of a LinkedHashSet. I later changed this to a LinkedHashSet when I realized that LinkedHashSet is much faster than a TreeSet)
Expected output (changing i < 50000 inside the findRoots() method)
public static void findRoots() {
BracketingNthOrderBrentSolver f = new
BracketingNthOrderBrentSolver();
UnivariateFunction func = (double x) -> RiemennZ(x, 4);
LinkedHashSet<Double> set = new LinkedHashSet<>();
double i = 1.0;
while (i < 50000) {
i+= 0.1;
if(sign(func.value(i)) != sign(func.value(i+0.1))) {
set.add(f.solve(1000, func, i, i+0.1));
}
}
double EPSILON = .05;
Iterator<Double> itr = set.iterator();
Double prevVal = null;
while(itr.hasNext()) {
Double currentVal = itr.next();
if (prevVal != null) {
if(currentVal - prevVal < EPSILON ) {
System.out.println(prevVal);
System.out.println(currentVal);
}
}
prevVal = currentVal;
}
}
Output. The two zeroes are within EPSILON (0.05 of each other). The reason why this is important deals with the mathematics behind the Z(t) function. EPSILON can be adjusted to find zeroes which are closer together.
Zeroes inside the critical line for Zeta(1/2 + it). The t values are referenced below.
5229.198557200015
5229.2418112597425
7005.062866174953
7005.100564672575
17143.786536183896
17143.82184350522
33179.36529436468
33179.40157549113
42525.79593689609
42525.835168267266
I think the best way to do this would be to simply store the previous value, and skip the first entry. So it would look something like the following:
double EPSILON = .001;
Iterator<Double> itr = set.iterator();
Double prevVal = null;
while(itr.hasNext()) {
Double currentVal = itr.next();
if (prevVal != null) {
if(currentVal - prevVal < EPSILON ) {
System.out.println(currentVal);
}
}
prevVal = currentVal;
}
This way you do not need to do a lot of manipulation with the iterator
Furthermore, I moved EPSILON outside of the while loop. Since it is a constant, it is best practice to have it outside of any loops and methods.