Conversion double en BigDecimal en Java
j'ai écrit un programme Java qui calcule des valeurs pour la fonction Zeta de Riemann . À l'intérieur du programme, j'ai fait une bibliothèque pour calculer les fonctions complexes nécessaires telles que atan, cos, etc. Tout ce qui se trouve dans les deux programmes est accessible par les types de données double et BigDecimal . Cela crée des problèmes majeurs lors de l'évaluation des valeurs de la fonction Zêta.
l'approximation numérique pour les références de la fonction Zeta
l'évaluation directe de cette approximation à des valeurs élevées crée des problèmes lorsque s a une grande forme complexe, comme s = (230+30i) . Je suis très reconnaissant d'obtenir des informations sur ce ici . L'évaluation de S2.minus(S1) crée des erreurs parce que j'ai écrit quelque chose de mal dans la méthode adaptiveQuad .
à titre d'exemple, Zeta(2+3i) à travers ce programme génère
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.
qui est correct .
Zeta(100+0i) génère
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.
, ce qui est également correct par rapport à Wolfram . Le problème est dû à quelque chose à l'intérieur de la méthode appelée adaptiveQuad .
Zeta(230+30i) génère
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.
l'imaginaire partie est un peu hors par rapport à Wolfram .
L'algorithme pour évaluer l'intégrale est connu comme "1519560920 Adaptatives" Quadrature et un double Java application ici . La méthode du quad adaptatif s'applique comme suit:
// 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);
}
voici ma quatrième tentative d'écrire le programme
/**************************************************************************
**
** 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);
}
}
nombres complexes ( 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);
}
}
j'examine la réponse ci-dessous.
un problème peut être dû à
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);
Je n'applique pas BigDecimal.pow(BigDecimal) . Bien que, je ne pense pas que ce soit la question directe qui provoque l'arithmétique flottante de créer des différences.
Edit : j'ai essayé une nouvelle approximation intégrale de la fonction Zeta. En fin de compte, je développerai une nouvelle méthode pour calculer BigDecimal.pow(BigDecimal) .
3 réponses
"Caveat je suis d'accord avec tous les commentaires dans la réponse de @laune , mais j'ai l'impression que vous pourriez vouloir poursuivre ceci de toute façon. Assurez - vous surtout que vous comprenez vraiment 1) et ce que cela signifie pour vous-il est très facile de faire beaucoup de calculs lourds pour produire des résultats sans signification.
à virgule flottante en précision Arbitraire fonctions en Java
pour réitérer un peu, je pense que votre problème est vraiment avec les mathématiques et la méthode numérique que vous avez choisi, mais voici une implémentation en utilisant la Bibliothèque Apfloat . Je vous conseille vivement d'utiliser la bibliothèque de précision arbitraire (ou similaire), car elle vous évite d'avoir à "lancer vos propres" fonctions de précision arbitraire (telles que pow , exp , sin , atan etc). Vous dites
finalement, je développerai une nouvelle méthode pour calculer BigDecimal.pow (BigDecimal)
c'est vraiment difficile d'obtenir ce droit.
vous devez surveiller la précision de vos constantes, aussi - notez que j'utilise une implémentation D'échantillon Apfloat pour calculer PI à un grand nombre (pour une certaine définition de grand!) de sig figues. Je suis dans une certaine mesure confiant que la bibliothèque Apfloat utilise des valeurs suffisamment précises pour e en exponentiation - la source est disponible si vous voulez vérifier.
différentes formulations intégrales pour calculer zeta
vous mettez en place trois différentes méthodes basées sur l'intégration dans l'une de vos éditions:
Celui étiqueté 25.5.12 est celui que vous avez actuellement dans la question et (bien que cela puisse être calculé à zéro facilement), il est difficile de travailler avec en raison de 2) dans la réponse de @laune . J'ai mis en œuvre 25.5.12 comme integrand1() dans le code - je vous exhorte à tracer avec la gamme de t pour différents s = a + 0i et de comprendre comment il se comporte. Ou regardez les tracés dans le article zeta sur le MathWorld de Wolfram. Celui étiqueté 25.5.11 I implémenté via integrand2() et le code dans la configuration Je publie ci-dessous.
Code
alors que je suis un peu réticent à poster du code qui trouvera sans doute des résultats erronés dans certaines configurations en raison de toutes les choses ci - dessus-j'ai encodé ce que vous essayez de faire ci-dessous, en utilisant des objets de précision flottants arbitraires pour les variables.
si vous voulez changer la formulation que vous utilisez (par exemple de 25.5.11 à 25.5.12), vous pouvez changer quelle fonction d'enrubannage f() retourne ou, mieux encore, change adaptiveQuad pour prendre dans un arbitraire méthode d'integrand enveloppé dans un class avec une interface... Vous devrez également modifier l'arithmétique dans findZeta() si vous voulez utiliser une des autres formulations intégrales.
Jouez avec les constantes au début du contenu de votre cœur. Je n'ai pas testé beaucoup de combinaisons, je pense que les problèmes mathématiques ici remplacer la programmation.
Je l'ai laissé configuré pour faire 2+3i dans environ 2000 appels à la méthode de la quadrature adaptative et correspondent aux 15 premiers chiffres de la valeur de Wolfram.
j'ai testé qu'il fonctionne toujours avec PRECISION = 120l et EPSILON=1e-15 . Le programme fait correspondre Wolfram alpha dans les 18 premiers chiffres significatifs pour les trois cas de test que vous fournissez. Le dernier ( 230+30i ) prend beaucoup de temps, même sur un ordinateur rapide - il appelle l'integrand fucntion certains de plus de 100 000 fois. Notez que j'utilise 40 pour la valeur de INFINITY dans l'intégrale - pas très haut, mais des valeurs plus élevées présenter le problème 1) comme déjà discuté...
N.B. C'est pas rapide (vous mesurerez en minutes ou en heures, pas en secondes - mais vous obtenez vraiment rapide si vous voulez accepter que 10^-15 ~= 10^-70 comme la plupart des gens serait!!). Il vous donnera quelques chiffres qui correspondent à Wolfram Alpha ;) vous pourriez vouloir prendre PRECISION vers le bas à environ 20 , INFINITY à 10 et EPSILON à 1e-10 pour vérifier quelques résultats avec petit s d'abord... Je l'ai laissé dans une certaine impression ainsi il vous dit chaque 100e fois adaptiveQuad est appelé pour le confort.
Reiteration quelle que soit la qualité de votre précision - il est ne va pas surmonter les caractéristiques mathématiques des fonctions impliquées dans cette façon de calculer zeta. I doute fort C'est ainsi que Wolfram alpha le fait, par exemple. Recherchez les méthodes de sommation série si vous voulez plus de méthodes tractable.
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;
}
}
Une note sur la "correction"
notez que dans votre réponse - vous êtes appelant zeta(2+3i) et zeta(100) "correct" par rapport à Wolfram quand ils présentent des erreurs de ~ 1e-10 et ~ 1e-9 respectivement (ils diffèrent dans la 10ème et 9ème décimale), mais vous êtes inquiet au sujet de zeta(230+30i) parce qu'il présente une erreur d'ordre 10e-14 dans la composante imaginaire ( 38e-15 vs 5e-70 qui sont tous deux très près de zéro). Donc, dans certains sens, celui que vous appelez "mauvais" est plus proche de la valeur du Wolfram que ceux que vous appelez "correct." Peut-être vous êtes inquiet que les chiffres sont différents, mais ce n'est pas vraiment une mesure de précision.
note finale
sauf si vous le faites pour apprendre comment les fonctions se comportent et comment la précision de la virgule flottante interagit avec elle - ne faites pas les choses de cette façon . Même Apfloat la documentation propre dit:
ce paquet est conçu pour une précision extrême. Le résultat aurait quelques chiffres de moins que prévu (environ 10) et les derniers (environ 10) chiffres dans le résultat peut être inexact. Si vous prévoyez d'utiliser nombre avec seulement quelques centaines de chiffres, utilisez un programme comme PARI (il est gratuit et disponible à partir de ftp://megrez.math.u-bordeaux.fr ) ou programme commercial comme Mathematica ou Maple si possible.
j'ajouterais mpmath en python à cette liste comme alternative libre maintenant.
(1) l'intégration utilise adaptQuad, en commençant par un intervalle [0,10]. Pour z=a + ib avec des valeurs de plus en plus grandes de a et b=0, l'integrand est une fonction de plus en plus oscillante, avec le nombre de zéros dans [0,5] seul étant proportionnel à a et augmentant à 43 pour z=100.
par conséquent, commencer l'approximation avec un intervalle qui inclut un ou plusieurs zéros est risqué, car le programme tel que affiché montre très clairement. Pour z=100, l'integrand est 0, -2,08 E-78 et 7,12 E-115 à 0, 5 et 10, respectivement. Par conséquent, en comparant le résultat de la formule de Simpson à 1E-20 retourne vrai, et le résultat est absolument faux.
(2) le calcul dans la méthode AbelPlana implique deux nombres complexes, C1 et C2. Pour z=a+0i, ils sont réels, et le tableau ci-dessous montre leurs valeurs pour diverses valeurs de A:
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
maintenant nous savons que les valeurs de ζ(a+0i) diminuent vers 1 pour augmenter A. Il est clair il est impossible pour deux valeurs supérieures à 1E15 de produire un résultat significatif près de l'une lorsqu'elles sont soustraites l'une de l'autre.
le tableau suggère également que pour un bon résultat de ζ(a+0i) en utilisant cet algorithme, C1 et C2*I (I est l'intégrale) doivent être calculés avec une précision d'environ 45 chiffres significatifs. (Les calculs de précision arbitraires n'évitent pas l'écueil décrit dans (1).)
(3) noter que lors de l'utilisation d'une bibliothèque avec une précision arbitraire, des valeurs telles que E et PI doit être fourni avec une meilleure précision que les valeurs doubles en java.lang.Les maths peuvent offrir.
Modifier (25.5.11) a autant de zéros dans [0,10] (25.5.12). Le calcul à 0 est difficile, mais ce n'est pas une singularité. Il élude la question (2).
pour une réponse concernant l'utilisation de l'arithmétique de précision arbitraire avec la méthode intégrale décrite dans L'OP-voir mon autre réponse
cependant , j'ai été intrigué par cela et j'ai pensé qu'une méthode de somme de série devrait être plus numériquement stable. J'ai trouvé la représentation série de Dirichlet sur Wikipedia et l'ai mis en œuvre (code entièrement exécutable ci-dessous).
Cela m'a donné un aperçu intéressant. Si je mets la convergence EPSILON à 1e-30 j'obtiens exactement les mêmes chiffres et l'exposant (i.e. 1e-70 dans la partie imaginaire) que Wolfram pour zeta(100) et zeta(230+ 30i) et l'algorithme se termine après seulement 1 ou 2 Termes s'ajoutant à la somme. Cela me suggère deux choses:
- Wolfram alpha utilise cette méthode de Somme ou quelque chose de similaire pour calculer les valeurs qu'il renvoie.
- la"justesse" de ces valeurs est difficile à évaluer. Par exemple-zeta(100) a une valeur exacte en termes de PI, donc peut être jugé. Je ne sais pas si cette estimation de
zeta(230+30i)est meilleure ou pire que celle trouvée par la méthode intégrale - cette méthode est vraiment assez lente à converger à
zeta(2+3i)et peut avoir besoinEPSILONprendre plus bas pour être utilisable.
j'ai aussi trouvé un papier académique qui est un compendium de méthodes numériques pour calculer zeta . Cela m'indique que le problème sous-jacent ici est certainement "non trivial"!!
quoi qu'il en soit - je laisse l'implémentation de la série sum ici comme une alternative pour tous ceux qui pourraient LA croiser à l'avenir.
import java.io.PrintWriter;
import org.apfloat.ApcomplexMath;
import org.apfloat.Apcomplex;
import org.apfloat.Apfloat;
import org.apfloat.ApfloatMath;
import org.apfloat.samples.Pi;
public class ZetaSeries {
//Number of sig figs accuracy. Note that infinite should be reserved
private static long PRECISION = 100l;
// Convergence criterion for integration
static Apfloat EPSILON = new Apfloat("1e-30",PRECISION);
static Apfloat one = new Apfloat("1",PRECISION);
static Apfloat two = new Apfloat("2",PRECISION);
static Apfloat minus_one = one.negate();
static Apfloat three = new Apfloat("3",PRECISION);
private Apcomplex s = null;
private Apcomplex s_plus_two = null;
public ZetaSeries(Apcomplex s) {
this.s = s;
this.s_plus_two = two.add(s);
}
public static void main(String[] args) {
Apfloat re = new Apfloat("230", PRECISION);
Apfloat im = new Apfloat("30", PRECISION);
Apcomplex s = new Apcomplex(re,im);
ZetaSeries z = new ZetaSeries(s);
System.out.println(z.findZeta());
}
private Apcomplex findZeta() {
Apcomplex series_sum = Apcomplex.ZERO;
Apcomplex multiplier = (one.divide(this.s.add(minus_one)));
int stop_condition = 1;
long n = 1;
while (stop_condition > 0)
{
Apcomplex term_to_add = sum_term(n);
stop_condition = ApcomplexMath.abs(term_to_add).compareTo(EPSILON);
series_sum = series_sum.add(term_to_add);
//if(n%50 == 0)
{
System.out.println("At iteration " + n + " : " + multiplier.multiply(series_sum));
}
n+=1;
}
return multiplier.multiply(series_sum);
}
private Apcomplex sum_term(long n_long) {
Apfloat n = new Apfloat(n_long, PRECISION);
Apfloat n_plus_one = n.add(one);
Apfloat two_n = two.multiply(n);
Apfloat t1 = (n.multiply(n_plus_one)).divide(two);
Apcomplex t2 = (two_n.add(three).add(this.s)).divide(ApcomplexMath.pow(n_plus_one,s_plus_two));
Apcomplex t3 = (two_n.add(minus_one).add(this.s.negate())).divide(ApcomplexMath.pow(n,this.s_plus_two));
return t1.multiply(t2.add(t3.negate()));
}
}