我编写了一个 Java 程序来计算 the Riemann Zeta Function 的值。在程序内部，我制作了一个库来计算必要的复杂函数，例如 atan、cos 等。两个程序中的所有内容都通过 double 和 BigDecimal 数据类型访问。这在评估 Zeta 函数的大值时会产生重大问题。
Zeta 函数引用的数值近似
当 s 具有大型复杂形式(例如 s = (230+30i) )时，直接以高值评估此近似值会产生问题。我非常感谢获得有关此 here 的信息。 S2.minus(S1) 的计算会产生错误，因为我在 adaptiveQuad 方法中写错了一些东西。
例如，Zeta(2+3i) 通过这个程序生成
Calculation of the Riemann Zeta Function in the form Zeta(s) = a + ib.
Enter the value of [a] inside the Riemann Zeta Function: 2
Enter the value of [b] inside the Riemann Zeta Function: 3
The value for Zeta(s) is 7.980219851133409E-1 - 1.137443081631288E-1*i
Total time taken is 0.469 seconds.
这是 correct 。
Zeta(100+0i) 生成Calculation of the Riemann Zeta Function in the form Zeta(s) = a + ib.
Enter the value of [a] inside the Riemann Zeta Function: 100
Enter the value of [b] inside the Riemann Zeta Function: 0
The value for Zeta(s) is 1.000000000153236E0
Total time taken is 0.672 seconds.
与 Wolfram 相比，这也是正确的。问题是由标记为
adaptiveQuad 的方法内部的某些内容引起的。Zeta(230+30i) 生成Calculation of the Riemann Zeta Function in the form Zeta(s) = a + ib.
Enter the value of [a] inside the Riemann Zeta Function: 230
Enter the value of [b] inside the Riemann Zeta Function: 30
The value for Zeta(s) is 0.999999999999093108519845391615339162047254997503854254342793916541606842461539820124897870147977114468145672577664412128509813042591501204781683860384769321084473925620572315416715721728082468412672467499199310913504362891199180150973087384370909918493750428733837552915328069343498987460727711606978118652477860450744628906250 - 38.005428584222228490409289204403133867487950535704812764806874887805043029499897666636162309572126423385487374863788363786029170239477119910868455777891701471328505006916099918492113970510619110472506796418206225648616641319533972054228283869713393805956289770456519729094756021581247296126093715429306030273437500E-15*i
Total time taken is 1.746 seconds.
与 Wolfram 相比，虚部有点偏离。
评估积分的算法称为 Adaptive Quadrature ，找到
double Java 实现 here 。自适应四边形方法应用以下内容// adaptive quadrature
public static double adaptive(double a, double b) {
double h = b - a;
double c = (a + b) / 2.0;
double d = (a + c) / 2.0;
double e = (b + c) / 2.0;
double Q1 = h/6 * (f(a) + 4*f(c) + f(b));
double Q2 = h/12 * (f(a) + 4*f(d) + 2*f(c) + 4*f(e) + f(b));
if (Math.abs(Q2 - Q1) <= EPSILON)
return Q2 + (Q2 - Q1) / 15;
else
return adaptive(a, c) + adaptive(c, b);
}
这是我第四次尝试编写程序
/**************************************************************************
**
** Abel-Plana Formula for the Zeta Function
**
**************************************************************************
** Axion004
** 08/16/2015
**
** This program computes the value for Zeta(z) using a definite integral
** approximation through the Abel-Plana formula. The Abel-Plana formula
** can be shown to approximate the value for Zeta(s) through a definite
** integral. The integral approximation is handled through the Composite
** Simpson's Rule known as Adaptive Quadrature.
**************************************************************************/
import java.util.*;
import java.math.*;
public class AbelMain5 extends Complex {
private static MathContext MC = new MathContext(512,
RoundingMode.HALF_EVEN);
public static void main(String[] args) {
AbelMain();
}
// Main method
public static void AbelMain() {
double re = 0, im = 0;
double start, stop, totalTime;
Scanner scan = new Scanner(System.in);
System.out.println("Calculation of the Riemann Zeta " +
"Function in the form Zeta(s) = a + ib.");
System.out.println();
System.out.print("Enter the value of [a] inside the Riemann Zeta " +
"Function: ");
try {
re = scan.nextDouble();
}
catch (Exception e) {
System.out.println("Please enter a valid number for a.");
}
System.out.print("Enter the value of [b] inside the Riemann Zeta " +
"Function: ");
try {
im = scan.nextDouble();
}
catch (Exception e) {
System.out.println("Please enter a valid number for b.");
}
start = System.currentTimeMillis();
Complex z = new Complex(new BigDecimal(re), new BigDecimal(im));
System.out.println("The value for Zeta(s) is " + AbelPlana(z));
stop = System.currentTimeMillis();
totalTime = (double) (stop-start) / 1000.0;
System.out.println("Total time taken is " + totalTime + " seconds.");
}
/**
* The definite integral for Zeta(z) in the Abel-Plana formula.
* <br> Numerator = Sin(z * arctan(t))
* <br> Denominator = (1 + t^2)^(z/2) * (e^(2*pi*t) - 1)
* @param t - the value of t passed into the integrand.
* @param z - The complex value of z = a + i*b
* @return the value of the complex function.
*/
public static Complex f(double t, Complex z) {
Complex num = (z.multiply(Math.atan(t))).sin();
Complex D1 = new Complex(1 + t*t).pow(z.divide(TWO));
Complex D2 = new Complex(Math.pow(Math.E, 2.0*Math.PI*t) - 1.0);
Complex den = D1.multiply(D2);
return num.divide(den, MC);
}
/**
* Adaptive quadrature - See http://www.mathworks.com/moler/quad.pdf
* @param a - the lower bound of integration.
* @param b - the upper bound of integration.
* @param z - The complex value of z = a + i*b
* @return the approximate numerical value of the integral.
*/
public static Complex adaptiveQuad(double a, double b, Complex z) {
double EPSILON = 1E-10;
double step = b - a;
double c = (a + b) / 2.0;
double d = (a + c) / 2.0;
double e = (b + c) / 2.0;
Complex S1 = (f(a, z).add(f(c, z).multiply(FOUR)).add(f(b, z))).
multiply(step / 6.0);
Complex S2 = (f(a, z).add(f(d, z).multiply(FOUR)).add(f(c, z).multiply
(TWO)).add(f(e, z).multiply(FOUR)).add(f(b, z))).multiply
(step / 12.0);
Complex result = (S2.subtract(S1)).divide(FIFTEEN, MC);
if(S2.subtract(S1).mod() <= EPSILON)
return S2.add(result);
else
return adaptiveQuad(a, c, z).add(adaptiveQuad(c, b, z));
}
/**
* The definite integral for Zeta(z) in the Abel-Plana formula.
* <br> value = 1/2 + 1/(z-1) + 2 * Integral
* @param z - The complex value of z = a + i*b
* @return the value of Zeta(z) through value and the
* quadrature approximation.
*/
public static Complex AbelPlana(Complex z) {
Complex C1 = ONEHALF.add(ONE.divide(z.subtract(ONE), MC));
Complex C2 = TWO.multiply(adaptiveQuad(1E-16, 100.0, z));
if ( z.real().doubleValue() == 0 && z.imag().doubleValue() == 0)
return new Complex(0.0, 0.0);
else
return C1.add(C2);
}
}
复数 (
BigDecimal )/**************************************************************************
**
** Complex Numbers
**
**************************************************************************
** Axion004
** 08/20/2015
**
** This class is necessary as a helper class for the calculation of
** imaginary numbers. The calculation of Zeta(z) inside AbelMain is in
** the form of z = a + i*b.
**************************************************************************/
import java.math.BigDecimal;
import java.math.MathContext;
import java.text.DecimalFormat;
import java.text.NumberFormat;
public class Complex extends Object{
private BigDecimal re;
private BigDecimal im;
/**
BigDecimal constant for zero
*/
final static Complex ZERO = new Complex(BigDecimal.ZERO) ;
/**
BigDecimal constant for one half
*/
final static Complex ONEHALF = new Complex(new BigDecimal(0.5));
/**
BigDecimal constant for one
*/
final static Complex ONE = new Complex(BigDecimal.ONE);
/**
BigDecimal constant for two
*/
final static Complex TWO = new Complex(new BigDecimal(2.0));
/**
BigDecimal constant for four
*/
final static Complex FOUR = new Complex(new BigDecimal(4.0)) ;
/**
BigDecimal constant for fifteen
*/
final static Complex FIFTEEN = new Complex(new BigDecimal(15.0)) ;
/**
Default constructor equivalent to zero
*/
public Complex() {
re = BigDecimal.ZERO;
im = BigDecimal.ZERO;
}
/**
Constructor with real part only
@param x Real part, BigDecimal
*/
public Complex(BigDecimal x) {
re = x;
im = BigDecimal.ZERO;
}
/**
Constructor with real part only
@param x Real part, double
*/
public Complex(double x) {
re = new BigDecimal(x);
im = BigDecimal.ZERO;
}
/**
Constructor with real and imaginary parts in double format.
@param x Real part
@param y Imaginary part
*/
public Complex(double x, double y) {
re= new BigDecimal(x);
im= new BigDecimal(y);
}
/**
Constructor for the complex number z = a + i*b
@param re Real part
@param im Imaginary part
*/
public Complex (BigDecimal re, BigDecimal im) {
this.re = re;
this.im = im;
}
/**
Real part of the Complex number
@return Re[z] where z = a + i*b.
*/
public BigDecimal real() {
return re;
}
/**
Imaginary part of the Complex number
@return Im[z] where z = a + i*b.
*/
public BigDecimal imag() {
return im;
}
/**
Complex conjugate of the Complex number
in which the conjugate of z is z-bar.
@return z-bar where z = a + i*b and z-bar = a - i*b
*/
public Complex conjugate() {
return new Complex(re, im.negate());
}
/**
* Returns the sum of this and the parameter.
@param augend the number to add
@param mc the context to use
@return this + augend
*/
public Complex add(Complex augend,MathContext mc)
{
//(a+bi)+(c+di) = (a + c) + (b + d)i
return new Complex(
re.add(augend.re,mc),
im.add(augend.im,mc));
}
/**
Equivalent to add(augend, MathContext.UNLIMITED)
@param augend the number to add
@return this + augend
*/
public Complex add(Complex augend)
{
return add(augend, MathContext.UNLIMITED);
}
/**
Addition of Complex number and a double.
@param d is the number to add.
@return z+d where z = a+i*b and d = double
*/
public Complex add(double d){
BigDecimal augend = new BigDecimal(d);
return new Complex(this.re.add(augend, MathContext.UNLIMITED),
this.im);
}
/**
* Returns the difference of this and the parameter.
@param subtrahend the number to subtract
@param mc the context to use
@return this - subtrahend
*/
public Complex subtract(Complex subtrahend, MathContext mc)
{
//(a+bi)-(c+di) = (a - c) + (b - d)i
return new Complex(
re.subtract(subtrahend.re,mc),
im.subtract(subtrahend.im,mc));
}
/**
* Equivalent to subtract(subtrahend, MathContext.UNLIMITED)
@param subtrahend the number to subtract
@return this - subtrahend
*/
public Complex subtract(Complex subtrahend)
{
return subtract(subtrahend,MathContext.UNLIMITED);
}
/**
Subtraction of Complex number and a double.
@param d is the number to subtract.
@return z-d where z = a+i*b and d = double
*/
public Complex subtract(double d){
BigDecimal subtrahend = new BigDecimal(d);
return new Complex(this.re.subtract(subtrahend, MathContext.UNLIMITED),
this.im);
}
/**
* Returns the product of this and the parameter.
@param multiplicand the number to multiply by
@param mc the context to use
@return this * multiplicand
*/
public Complex multiply(Complex multiplicand, MathContext mc)
{
//(a+bi)(c+di) = (ac - bd) + (ad + bc)i
return new Complex(
re.multiply(multiplicand.re,mc).subtract(im.multiply
(multiplicand.im,mc),mc),
re.multiply(multiplicand.im,mc).add(im.multiply
(multiplicand.re,mc),mc));
}
/**
Equivalent to multiply(multiplicand, MathContext.UNLIMITED)
@param multiplicand the number to multiply by
@return this * multiplicand
*/
public Complex multiply(Complex multiplicand)
{
return multiply(multiplicand,MathContext.UNLIMITED);
}
/**
Complex multiplication by a double.
@param d is the double to multiply by.
@return z*d where z = a+i*b and d = double
*/
public Complex multiply(double d){
BigDecimal multiplicand = new BigDecimal(d);
return new Complex(this.re.multiply(multiplicand, MathContext.UNLIMITED)
,this.im.multiply(multiplicand, MathContext.UNLIMITED));
}
/**
Modulus of a Complex number or the distance from the origin in
* the polar coordinate plane.
@return |z| where z = a + i*b.
*/
public double mod() {
if ( re.doubleValue() != 0.0 || im.doubleValue() != 0.0)
return Math.sqrt(re.multiply(re).add(im.multiply(im))
.doubleValue());
else
return 0.0;
}
/**
* Modulus of a Complex number squared
* @param z = a + i*b
* @return |z|^2 where z = a + i*b
*/
public double abs(Complex z) {
double doubleRe = re.doubleValue();
double doubleIm = im.doubleValue();
return doubleRe * doubleRe + doubleIm * doubleIm;
}
public Complex divide(Complex divisor)
{
return divide(divisor,MathContext.UNLIMITED);
}
/**
* The absolute value squared.
* @return The sum of the squares of real and imaginary parts.
* This is the square of Complex.abs() .
*/
public BigDecimal norm()
{
return re.multiply(re).add(im.multiply(im)) ;
}
/**
* The absolute value of a BigDecimal.
* @param mc amount of precision
* @return BigDecimal.abs()
*/
public BigDecimal abs(MathContext mc)
{
return BigDecimalMath.sqrt(norm(),mc) ;
}
/** The inverse of the the Complex number.
@param mc amount of precision
@return 1/this
*/
public Complex inverse(MathContext mc)
{
final BigDecimal hyp = norm() ;
/* 1/(x+iy)= (x-iy)/(x^2+y^2 */
return new Complex( re.divide(hyp,mc), im.divide(hyp,mc)
.negate() ) ;
}
/** Divide through another BigComplex number.
@param oth the other complex number
@param mc amount of precision
@return this/other
*/
public Complex divide(Complex oth, MathContext mc)
{
/* implementation: (x+iy)/(a+ib)= (x+iy)* 1/(a+ib) */
return multiply(oth.inverse(mc),mc) ;
}
/**
Division of Complex number by a double.
@param d is the double to divide
@return new Complex number z/d where z = a+i*b
*/
public Complex divide(double d){
BigDecimal divisor = new BigDecimal(d);
return new Complex(this.re.divide(divisor, MathContext.UNLIMITED),
this.im.divide(divisor, MathContext.UNLIMITED));
}
/**
Exponential of a complex number (z is unchanged).
<br> e^(a+i*b) = e^a * e^(i*b) = e^a * (cos(b) + i*sin(b))
@return exp(z) where z = a+i*b
*/
public Complex exp () {
return new Complex(Math.exp(re.doubleValue()) * Math.cos(im.
doubleValue()), Math.exp(re.doubleValue()) *
Math.sin(im.doubleValue()));
}
/**
The Argument of a Complex number or the angle in radians
with respect to polar coordinates.
<br> Tan(theta) = b / a, theta = Arctan(b / a)
<br> a is the real part on the horizontal axis
<br> b is the imaginary part of the vertical axis
@return arg(z) where z = a+i*b.
*/
public double arg() {
return Math.atan2(im.doubleValue(), re.doubleValue());
}
/**
The log or principal branch of a Complex number (z is unchanged).
<br> Log(a+i*b) = ln|a+i*b| + i*Arg(z) = ln(sqrt(a^2+b^2))
* + i*Arg(z) = ln (mod(z)) + i*Arctan(b/a)
@return log(z) where z = a+i*b
*/
public Complex log() {
return new Complex(Math.log(this.mod()), this.arg());
}
/**
The square root of a Complex number (z is unchanged).
Returns the principal branch of the square root.
<br> z = e^(i*theta) = r*cos(theta) + i*r*sin(theta)
<br> r = sqrt(a^2+b^2)
<br> cos(theta) = a / r, sin(theta) = b / r
<br> By De Moivre's Theorem, sqrt(z) = sqrt(a+i*b) =
* e^(i*theta / 2) = r(cos(theta/2) + i*sin(theta/2))
@return sqrt(z) where z = a+i*b
*/
public Complex sqrt() {
double r = this.mod();
double halfTheta = this.arg() / 2;
return new Complex(Math.sqrt(r) * Math.cos(halfTheta), Math.sqrt(r) *
Math.sin(halfTheta));
}
/**
The real cosh function for Complex numbers.
<br> cosh(theta) = (e^(theta) + e^(-theta)) / 2
@return cosh(theta)
*/
private double cosh(double theta) {
return (Math.exp(theta) + Math.exp(-theta)) / 2;
}
/**
The real sinh function for Complex numbers.
<br> sinh(theta) = (e^(theta) - e^(-theta)) / 2
@return sinh(theta)
*/
private double sinh(double theta) {
return (Math.exp(theta) - Math.exp(-theta)) / 2;
}
/**
The sin function for the Complex number (z is unchanged).
<br> sin(a+i*b) = cosh(b)*sin(a) + i*(sinh(b)*cos(a))
@return sin(z) where z = a+i*b
*/
public Complex sin() {
return new Complex(cosh(im.doubleValue()) * Math.sin(re.doubleValue()),
sinh(im.doubleValue())* Math.cos(re.doubleValue()));
}
/**
The cos function for the Complex number (z is unchanged).
<br> cos(a +i*b) = cosh(b)*cos(a) + i*(-sinh(b)*sin(a))
@return cos(z) where z = a+i*b
*/
public Complex cos() {
return new Complex(cosh(im.doubleValue()) * Math.cos(re.doubleValue()),
-sinh(im.doubleValue()) * Math.sin(re.doubleValue()));
}
/**
The hyperbolic sin of the Complex number (z is unchanged).
<br> sinh(a+i*b) = sinh(a)*cos(b) + i*(cosh(a)*sin(b))
@return sinh(z) where z = a+i*b
*/
public Complex sinh() {
return new Complex(sinh(re.doubleValue()) * Math.cos(im.doubleValue()),
cosh(re.doubleValue()) * Math.sin(im.doubleValue()));
}
/**
The hyperbolic cosine of the Complex number (z is unchanged).
<br> cosh(a+i*b) = cosh(a)*cos(b) + i*(sinh(a)*sin(b))
@return cosh(z) where z = a+i*b
*/
public Complex cosh() {
return new Complex(cosh(re.doubleValue()) *Math.cos(im.doubleValue()),
sinh(re.doubleValue()) * Math.sin(im.doubleValue()));
}
/**
The tan of the Complex number (z is unchanged).
<br> tan (a+i*b) = sin(a+i*b) / cos(a+i*b)
@return tan(z) where z = a+i*b
*/
public Complex tan() {
return (this.sin()).divide(this.cos());
}
/**
The arctan of the Complex number (z is unchanged).
<br> tan^(-1)(a+i*b) = 1/2 i*(log(1-i*(a+b*i))-log(1+i*(a+b*i))) =
<br> -1/2 i*(log(i*a - b+1)-log(-i*a + b+1))
@return arctan(z) where z = a+i*b
*/
public Complex atan(){
Complex ima = new Complex(0.0,-1.0); //multiply by negative i
Complex num = new Complex(this.re.doubleValue(),this.im.doubleValue()
-1.0);
Complex den = new Complex(this.re.negate().doubleValue(),this.im
.negate().doubleValue()-1.0);
Complex two = new Complex(2.0, 0.0); // divide by 2
return ima.multiply(num.divide(den).log()).divide(two);
}
/**
* The Math.pow equivalent of two Complex numbers.
* @param z - the complex base in the form z = a + i*b
* @return z^y where z = a + i*b and y = c + i*d
*/
public Complex pow(Complex z){
Complex a = z.multiply(this.log(), MathContext.UNLIMITED);
return a.exp();
}
/**
* The Math.pow equivalent of a Complex number to the power
* of a double.
* @param d - the double to be taken as the power.
* @return z^d where z = a + i*b and d = double
*/
public Complex pow(double d){
Complex a=(this.log()).multiply(d);
return a.exp();
}
/**
Override the .toString() method to generate complex numbers, the
* string representation is now a literal Complex number.
@return a+i*b, a-i*b, a, or i*b as desired.
*/
public String toString() {
NumberFormat formatter = new DecimalFormat();
formatter = new DecimalFormat("#.###############E0");
if (re.doubleValue() != 0.0 && im.doubleValue() > 0.0) {
return formatter.format(re) + " + " + formatter.format(im)
+"*i";
}
if (re.doubleValue() !=0.0 && im.doubleValue() < 0.0) {
return formatter.format(re) + " - "+ formatter.format(im.negate())
+ "*i";
}
if (im.doubleValue() == 0.0) {
return formatter.format(re);
}
if (re.doubleValue() == 0.0) {
return formatter.format(im) + "*i";
}
return formatter.format(re) + " + i*" + formatter.format(im);
}
}
我正在查看下面的答案。
一个问题可能是由于
Complex num = (z.multiply(Math.atan(t))).sin();
Complex D1 = new Complex(1 + t*t).pow(z.divide(TWO));
Complex D2 = new Complex(Math.pow(Math.E, 2.0*Math.PI*t) - 1.0);
Complex den = D1.multiply(D2, MathContext.UNLIMITED);
我没有申请
BigDecimal.pow(BigDecimal) 。虽然，我不认为这是导致浮点运算产生差异的直接问题。编辑 :我尝试了 Zeta 函数的新积分近似。最终，我将开发一种新方法来计算
BigDecimal.pow(BigDecimal) 。
最佳答案
警告 我同意 @laune's answer 中的所有评论，但我的印象是你可能希望无论如何都要追求这个。尤其要确保您真的了解 1) 以及这对您意味着什么 - 很容易进行大量繁重的计算以产生无意义的结果。
Java中的任意精度浮点函数
重申一下，我认为您的问题确实与您选择的数学和数值方法有关，但这里有一个使用 Apfloat library 的实现。我强烈建议您使用现成的任意精度库(或类似的库)，因为它可以避免您“推出自己的”任意精度数学函数(例如 pow 、 exp 、 sin 、 atan 等)。你说
Ultimately, I will develop a new method to calculate BigDecimal.pow(BigDecimal)
真的很难做到 。
您还需要注意常量的精度 - 请注意，我使用 Apfloat 示例实现将
PI 计算为大量(对于大的某些定义!)的 sig figs。我在某种程度上相信 Apfloat 库在幂运算中使用了适当精确的 e 值 - 如果您想检查，可以使用源代码。计算zeta的不同积分公式
您在其中一个编辑中提出了三种不同的基于集成的方法:
标记为 25.5.12 的那个是您当前在问题中遇到的那个(尽管可以很容易地将其计算为零)，但由于 2) 在 @laune's answer 中很难使用。我在代码中将 25.5.12 实现为
integrand1() - 我敦促您使用不同 t 的 s = a + 0i 范围绘制它，并了解它的行为方式。或者查看 Wolfram 数学世界上 zeta article 中的图。我通过 integrand2() 和我在下面发布的配置中实现的标记为 25.5.11 的那个。代码
虽然我有点不愿意发布代码，但由于上述所有因素，这些代码无疑会在某些配置中发现错误的结果 - 我已经对您在下面尝试执行的操作进行了编码，对变量使用了任意精度的浮点对象。
如果您想更改您使用的公式(例如从 25.5.11 到 25.5.12)，您可以更改包装函数
f() 返回的被积函数，或者更好的是，更改 adaptiveQuad 以采用包装在 class 中的任意被积函数方法一个接口(interface)...如果您想使用其他积分公式之一，您还必须更改 findZeta() 中的算术。在开始时使用常量来满足您的需求。我没有测试过很多组合，因为我认为这里的数学问题覆盖了编程问题。
我已经将它设置为在大约 2000 次调用自适应正交方法时执行
2+3i 并匹配 Wolfram 值的前 15 位左右数字。我已经测试过它仍然适用于
PRECISION = 120l 和 EPSILON=1e-15 。对于您提供的三个测试用例，该程序在前 18 个左右的有效数字中匹配 Wolfram alpha。最后一个( 230+30i )即使在一台快速的计算机上也需要很长时间 - 它调用被积函数大约 100,000 次以上。请注意，我在积分中使用 40 作为 INFINITY 的值 - 不是很高 - 但更高的值表现出问题 1)，正如已经讨论过的......注意 这是 不快 (您将在几分钟或几小时内进行测量，而不是几秒钟-但如果您想像大多数人一样接受 10^-15 ~= 10^-70，您只会变得非常快!!) .它会给你一些与 Wolfram Alpha 匹配的数字；) 你可能想把
PRECISION 降到 about 20 ， INFINITY 到 10 和 EPSILON 到 1e-10 以先用小 s 验证一些结果......我已经留下了一些打印品它告诉你每 100 次 adaptiveQuad 被调用以获得舒适感。重复 不管你的精度有多好 - 它不会克服以这种方式计算 zeta 所涉及的函数的数学特征。例如，我 强烈怀疑 这就是 Wolfram alpha 的做法。如果您想要更易于处理的方法，请查找系列求和方法。
import java.io.PrintWriter;
import org.apfloat.ApcomplexMath;
import org.apfloat.Apcomplex;
import org.apfloat.Apfloat;
import org.apfloat.samples.Pi;
public class ZetaFinder
{
//Number of sig figs accuracy. Note that infinite should be reserved
private static long PRECISION = 40l;
// Convergence criterion for integration
static Apfloat EPSILON = new Apfloat("1e-15",PRECISION);
//Value of PI - enhanced using Apfloat library sample calculation of Pi in constructor,
//Fast enough that we don't need to hard code the value in.
//You could code hard value in for perf enhancement
static Apfloat PI = null; //new Apfloat("3.14159");
//Integration limits - I found too high a value for "infinity" causes integration
//to terminate on first iteration. Plot the integrand to see why...
static Apfloat INFINITE_LIMIT = new Apfloat("40",PRECISION);
static Apfloat ZERO_LIMIT = new Apfloat("1e-16",PRECISION); //You can use zero for the 25.5.12
static Apfloat one = new Apfloat("1",PRECISION);
static Apfloat two = new Apfloat("2",PRECISION);
static Apfloat four = new Apfloat("4",PRECISION);
static Apfloat six = new Apfloat("6",PRECISION);
static Apfloat twelve = new Apfloat("12",PRECISION);
static Apfloat fifteen = new Apfloat("15",PRECISION);
static int counter = 0;
Apcomplex s = null;
public ZetaFinder(Apcomplex s)
{
this.s = s;
Pi.setOut(new PrintWriter(System.out, true));
Pi.setErr(new PrintWriter(System.err, true));
PI = (new Pi.RamanujanPiCalculator(PRECISION+10, 10)).execute(); //Get Pi to a higher precision than integer consts
System.out.println("Created a Zeta Finder based on Abel-Plana for s="+s.toString() + " using PI="+PI.toString());
}
public static void main(String[] args)
{
Apfloat re = new Apfloat("2", PRECISION);
Apfloat im = new Apfloat("3", PRECISION);
Apcomplex s = new Apcomplex(re,im);
ZetaFinder finder = new ZetaFinder(s);
System.out.println(finder.findZeta());
}
private Apcomplex findZeta()
{
Apcomplex retval = null;
//Method currently in question (a.k.a. 25.5.12)
//Apcomplex mult = ApcomplexMath.pow(two, this.s);
//Apcomplex firstterm = (ApcomplexMath.pow(two, (this.s.add(one.negate())))).divide(this.s.add(one.negate()));
//Easier integrand method (a.k.a. 25.5.11)
Apcomplex mult = two;
Apcomplex firstterm = (one.divide(two)).add(one.divide(this.s.add(one.negate())));
Apfloat limita = ZERO_LIMIT;//Apfloat.ZERO;
Apfloat limitb = INFINITE_LIMIT;
System.out.println("Trying to integrate between " + limita.toString() + " and " + limitb.toString());
Apcomplex integral = adaptiveQuad(limita, limitb);
retval = firstterm.add((mult.multiply(integral)));
return retval;
}
private Apcomplex adaptiveQuad(Apfloat a, Apfloat b) {
//if (counter % 100 == 0)
{
System.out.println("In here for the " + counter + "th time");
}
counter++;
Apfloat h = b.add(a.negate());
Apfloat c = (a.add(b)).divide(two);
Apfloat d = (a.add(c)).divide(two);
Apfloat e = (b.add(c)).divide(two);
Apcomplex Q1 = (h.divide(six)).multiply(f(a).add(four.multiply(f(c))).add(f(b)));
Apcomplex Q2 = (h.divide(twelve)).multiply(f(a).add(four.multiply(f(d))).add(two.multiply(f(c))).add(four.multiply(f(e))).add(f(b)));
if (ApcomplexMath.abs(Q2.add(Q1.negate())).compareTo(EPSILON) < 0)
{
System.out.println("Returning");
return Q2.add((Q2.add(Q1.negate())).divide(fifteen));
}
else
{
System.out.println("Recursing with intervals "+a+" to " + c + " and " + c + " to " +d);
return adaptiveQuad(a, c).add(adaptiveQuad(c, b));
}
}
private Apcomplex f(Apfloat x)
{
return integrand2(x);
}
/*
* Simple test integrand (z^2)
*
* Can test implementation by asserting that the adaptiveQuad
* with this function evaluates to z^3 / 3
*/
private Apcomplex integrandTest(Apfloat t)
{
return ApcomplexMath.pow(t, two);
}
/*
* Abel-Plana formulation integrand
*/
private Apcomplex integrand1(Apfloat t)
{
Apcomplex numerator = ApcomplexMath.sin(this.s.multiply(ApcomplexMath.atan(t)));
Apcomplex bottomlinefirstbr = one.add(ApcomplexMath.pow(t, two));
Apcomplex D1 = ApcomplexMath.pow(bottomlinefirstbr, this.s.divide(two));
Apcomplex D2 = (ApcomplexMath.exp(PI.multiply(t))).add(one);
Apcomplex denominator = D1.multiply(D2);
Apcomplex retval = numerator.divide(denominator);
//System.out.println("Integrand evaluated at "+t+ " is "+retval);
return retval;
}
/*
* Abel-Plana formulation integrand 25.5.11
*/
private Apcomplex integrand2(Apfloat t)
{
Apcomplex numerator = ApcomplexMath.sin(this.s.multiply(ApcomplexMath.atan(t)));
Apcomplex bottomlinefirstbr = one.add(ApcomplexMath.pow(t, two));
Apcomplex D1 = ApcomplexMath.pow(bottomlinefirstbr, this.s.divide(two));
Apcomplex D2 = ApcomplexMath.exp(two.multiply(PI.multiply(t))).add(one.negate());
Apcomplex denominator = D1.multiply(D2);
Apcomplex retval = numerator.divide(denominator);
//System.out.println("Integrand evaluated at "+t+ " is "+retval);
return retval;
}
}
关于“正确性”的注释
请注意，在您的回答中 - 当
zeta(2+3i) 和 zeta(100) 分别显示 ~ 1e-10 和 ~ 1e-9 的错误时(它们在小数点后第 10 位和第 9 位不同)，您将 zeta(230+30i) 和 10e-14 称为“正确”，但您担心 38e-15 因为它在虚部( 5e-70 与 mpmath 两者都非常接近于零)中表现出 ojit_code 阶错误。因此，在某种意义上，您称为“错误”的值比您称为“正确”的值更接近 Wolfram 值。也许您担心前导数字不同，但这并不是衡量准确性的真正标准。最后一点
除非您这样做是为了了解函数的行为方式以及浮点精度如何与之交互 - 不要这样做 。甚至 Apfloat's own documentation 说:
This package is designed for extreme precision. The result might have a few digits less than you'd expect (about 10) and the last few (about 10) digits in ther result might be inaccurate. If you plan to use numbers with only a few hundred digits, use a program like PARI (it's free and available from ftp://megrez.math.u-bordeaux.fr) or a commercial program like Mathematica or Maple if possible.
我现在将 ojit_code in python 添加到此列表中作为免费替代方案。
关于java - 在 Java 中将 double 转换为 BigDecimal，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/32471327/