在Java中将double转换为BigDecimal
Axi*_*004 24 java double bigdecimal complex-numbers
我写了一个Java程序来计算Riemann Zeta函数的值.在程序内部,我创建了一个库来计算必要的复杂函数,如atan,cos等.两个程序中的所有内容都通过double和BigDecimal数据类型访问.这在评估Zeta函数的大值时会产生重大问题.
Zeta函数参考的数值近似
在s具有较大复杂形式的情况下,以高值直接评估此近似会产生问题,例如s = (230+30i).我非常感谢在这里获得有关这方面的信息.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.
这是正确的.
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相比,虚部有点偏.
评估积分的算法称为自适应正交,此处double可以找到Java实现.自适应四元法应用以下内容
// 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的回答中的所有评论,但我得到的印象是你无论如何都希望这样做.尤其要确保你真正理解1)以及这对你意味着什么 - 很容易做大量繁重的计算以产生毫无意义的结果.
Java中的任意精度浮点函数
重申一点,我认为你的问题实际上是你选择的数学和数值方法,但这是使用Apfloat库的实现.我强烈敦促你,因为它避免了任何需要你"滚你自己的"任意精度的数学函数(如使用现成的任意精度库(或类似的一个)pow,exp,sin,atan等).你说
最后,我将开发一种新方法来计算BigDecimal.pow(BigDecimal)
这真的很难得到这一权利.
您还需要观察常量的精度 - 请注意,我使用Apfloat示例实现来计算PI大量(对于某些大型定义!)的sig figs.我在某种程度上相信Apfloat库e在求幂中使用了适当的精确值- 如果你想检查,可以使用源代码.
不同的积分配方来计算zeta
您在其中一个编辑中设置了三种不同的基于集成的方法:
标记为25.5.12的那个是你目前在问题中的那个(虽然可以很容易地计算为零),但由于2)在@ laune的答案中很难使用.我integrand1()在代码中实现了25.5.12 - 我敦促你用t不同的范围绘制它,s = a + 0i并理解它的行为方式.或者看看Wolfram的数学世界中zeta文章中的情节.标记为25.5.11的那个我通过实现integrand2()和我在下面发布的配置中的代码.
码
虽然我有点不愿意发布代码,由于上述所有内容,在某些配置中无疑会发现错误的结果 - 我已经编码了你在下面尝试做的事情,使用任意精度浮点对象作为变量.
如果你想改变你使用的公式(例如从25.5.11到25.5.12),你可以改变包装器函数f()返回哪个被整数,或者更好的是,改变adaptiveQuad以包含一个class包含在接口中的任意被积函数方法. .. findZeta()如果你想使用其他积分配方之一,你还必须改变算术.
在开始时使用常数来玩你心中的内容.我没有测试过很多组合,因为我认为这里的数学问题会覆盖编程问题.
我已经将其设置为2+3i在大约2000次调用自适应四元组方法中进行,并匹配Wolfram值的前15位左右.
我测试过它仍然适用于PRECISION = 120l和EPSILON=1e-15.该程序在您提供的三个测试用例的前18个左右的重要数字中匹配Wolfram alpha.230+30i即使在快速计算机上,最后一个()也需要很长时间 - 它会将被积函数调用大约100,000次.请注意,我使用积分中40的值INFINITY- 不是很高 - 但更高的值表示问题1)已经讨论过......
注意这不是很快(你会在几分钟或几小时内测量,而不是几秒钟 - 但如果你想接受10 ^ -15~ = 10 ^ -70,你就会很快得到,就像大多数人一样!).它会给你一些匹配Wolfram Alpha的是数字;)你可能要考虑PRECISION到约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;
}
}
关于"正确性"的说明
请注意,在您的答案中 - 与Wolfram相比,当您分别显示〜和〜的错误(它们在小数点后第10位和第9位有所不同)时,您正在调用zeta(2+3i)并zeta(100)"正确" ,但您担心的是因为它出现了错误的顺序在虚构的组件中(对比都非常接近于零).因此,在某种意义上,你称之为"错误"的那个比你称之为"正确"的那个更接近Wolfram值.也许你担心领先的数字是不同的,但这并不是衡量准确性的标准.1e-101e-9zeta(230+30i)10e-1438e-155e-70
最后一点
除非你这样做是为了了解函数的行为以及浮点精度如何与它相互作用 - 不要这样做.即使是Apfloat自己的文档也说:
该包装设计用于极高精度.结果可能比您预期的数字少一些(大约10),结果中的最后几个(大约10)数字可能不准确.如果您打算使用只有几百位的数字,请使用像PARI这样的程序(它是免费的,可从ftp://megrez.math.u-bordeaux.fr获得)或者像Mathematica或Maple这样的商业程序.
我现在将mpmathpython作为免费替代添加到此列表中.
(1)积分使用adaptQuad,从区间[0,10]开始.对于z = a + ib,a和b = 0的值越来越大,被积函数是一个越来越频繁的函数,[0,5]中的零个数单独与a成正比,z = 100时上升到43.
因此,以包含一个或多个零的间隔开始近似是有风险的,因为发布的程序非常清楚地显示.对于z = 100,被积函数分别为0,-2.08E-78和7.12E-115,分别为0,5和10.因此,将辛普森公式的结果与1E-20进行比较,结果为真,结果绝对错误.
(2)方法AbelPlana中的计算涉及两个复数,C1和C2.对于z = a + 0i,它们是实数,下表显示了它们的各种值的值:
a C1 C2
10 5.689E1 1.024E3
20 2.759E4 1.048E6
30 1.851E7 1.073E9
40 1.409E10 1.099E12
60 9.770E15 1.152E18
100 6.402E27 1.267E30
现在我们知道ζ(a + 0i)的值随着增加a而减小到1.当1E15以上的两个值显然不可能在彼此相减时在一个附近产生有意义的结果.
该表还表明,通过使用该算法,对于ζ(a + 0i)的良好结果,需要以大约45位有效数字的精度计算C1和C2*I(I是积分).(任意精度数学不能避免(1)中描述的陷阱.)
(3)请注意,当使用具有任意精度的库时,应提供诸如E和PI的值,其精度要比java.lang.Math中提供的double值更好.
编辑 (25.5.11)在[0,10]中具有与(25.5.12)一样多的零.0处的计算很棘手,但它不是奇点.它确实避免了问题(2).
| 归档时间： |
|
| 查看次数： |
2556 次 |
| 最近记录： |