【发布时间】：2015-12-26 10:53:17
【问题描述】：
编辑：预期的输出在最后发布。
我编写了一个程序，通过Riemann-Siegel formula 计算Riemann Zeta Function 的零。我有兴趣调整程序中的一种方法来观察所谓的Lehmer's Phenomenon。
在程序内部，我正在寻找调整findRoots()方法的结尾。
完整的程序参考
/**************************************************************************
**
** 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));
}
}
打印的零包括
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
...等等。为了观察 Lehmer 的现象，我希望调整以下内容
Iterator<Double> itr = set.iterator();
while(itr.hasNext()) {
System.out.println(itr.next());
}
类似于
// 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());
}
只有当两个值接近时才会打印零。我可以用迭代器做到这一点吗？我需要改用 ListIterator 吗？
我尝试了增强的 for 循环并遇到了同样的问题。是否有更好或更有效的方法来比较从 LinkedHashSet 集中打印的值？
（注意：我最初创建的是 TreeSet 而不是 LinkedHashSet。后来当我意识到 LinkedHashSet 比 TreeSet 快得多时，我将其更改为 LinkedHashSet）
预期输出（在 findRoots() 方法中更改 i < 50000）
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;
}
}
输出。这两个零在 EPSILON 之内（彼此相差 0.05）。这很重要的原因在于 Z(t) 函数背后的数学。可以调整 EPSILON 以找到更靠近的零点。
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
【问题讨论】：
-
为什么需要Set？可以重复吗？
-
设置是去重复。很难对重复的数据进行排序。
-
你能发布预期的输出吗？
标签： java iterator treeset linkedhashset