2
我写了一个计算Riemann-Siegel Z（t）函数的基本程序。我很好奇，如果有更好的方法来近似剩余期限。我现在使用的方法使用了来自Haselgrove的可怕表格近似值。在Java中实现Riemann Siegel公式，对改进余数项好奇
Further information about the Riemann Siegel formula。这可能稍微先进一点，但有关这方面的更多详细信息，请参见this thesis。另外，在Edwards book。
我知道我的for循环不是最优的，我只是用它来进行测试。我可以使用不同的方法来近似零。
下面是我写的实现：
import java.util.*;
public class Main {
public static void main(String[] args) {
double s[] = new double[10];
s[0] = 2;
for (double i = 0; i < 500; i += 0.0001) {
if (RiemennZ(i, 4) < 0.0001 && RiemennZ(i, 4) > -1*0.0001)
System.out.println("Found a zero at " + i + ", the value of Zeta(s) is "
+ RiemennZ(i, 4));
}
//System.out.println(4);
//System.out.println("Value of the Zeta Function " + Arrays.toString(Riemann.zeta(s)));
System.out.println("The function you wrote is- " + RiemennZ(16, 4));
System.out.println(fAbs(1.3) -1.0);
//System.out.println(theta(25));
}
// Riemann-Siegel theta function using the approximation by the Stirling series
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.
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
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;
// Add the remainder terms
double remainder;
double frac = val - n;
int k = 0;
double R = 0.0;
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
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)
+.000000* 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));
}
}
其余的术语定义如下：(cos[2pi(p^2-p-1/(16))])/(cos(2pip))
做Wolfram Alpha的内这个函数的导数多是一塌糊涂。有没有人曾经遇到过这种问题？
为了使用多余项，我需要计算多个衍生为：(cos[2pi(p^2-p-1/(16))])/(cos(2pip))
有一些解决这个方式，可以用Java来实现？
一种方法是使用有限差分方法。这不是一个很好的解决方案，但是我首先想到的。
// Derivation of the first C term using first order central difference
public static double firstDerivative(double p) {
double epsilon = 0.0000000001;
double d1, d2;
double dx = 0.00001;
double diff = 1.0;
d1 = (function(p + dx) - function(p - dx))/(2 * dx);
while (diff > epsilon) {
dx /= 2;
d2 = (function(p + dx) - function(p - dx))/(2 * dx);
diff = Math.abs(d2 - d1);
d1 = d2;
}
return d1;
}
// Derivation of the second C term using second order central difference
public static double secondDerivative(double p) {
double epsilon = 0.0000000001;
double d1, d2;
double dx = 0.00001;
double diff = 1.0;
d1 = (function(p + dx) - 2.0 * function(p) + function(p - dx))/Math.pow(dx, 2);
while (diff > epsilon) {
dx /= 2;
d2 = (function(p + dx) - 2.0 * function(p) + function(p - dx))/Math.pow(dx, 2);
diff = Math.abs(d2 - d1);
d1 = d2;
}
return d1;
}
其余的前8个C项将具有大量的准确性。这些上升到二十四个派生。让我回顾一下你提出的建议，我需要看看如何使用Java计算派生。 – Axion004
有关Riemann Siegel配方的另一个项目可在http://web.mit.edu/ kenta /www/six/parallel/2-Final-Report.html#Evaluating中找到。作者计算多达46个余数项。由于余数项变得非常小，因为t - >无穷大，所以没有必要全部参考46。 – Axion004