I wrote a program to calculate values for the Zeta function through the Abel-Plana formula.
Depending on the value of the error term referenced, Java returns an error message for Exception in thread "main" java.lang.StackOverflowError.
Here is a working copy of the program with test values. The Complex.Java program is a helper program.
/**************************************************************************
**
** 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 AbelMain extends Complex {
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(re, 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^(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, 0).pow(z.divide(2.0));
double D2 = Math.pow(Math.E, Math.PI * t) + 1.0;
Complex den = D1.multiply(D2);
return num.divide(den);
}
/**
* 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) {
BigDecimal EPSILON = new BigDecimal(1E-20);
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(4)).add(f(b, z))).
multiply(step / 6.0);
Complex S2 = (f(a, z).add(f(d, z).multiply(4)).add(f(c, z).multiply(2))
.add(f(e, z).multiply(4)).add(f(b, z))).multiply(step / 12.0);
Complex result = (S2.minus(S1)).divide(15.0);
if(BigDecimal.valueOf(S2.minus(S1).mod()).compareTo(EPSILON) == -1
|| BigDecimal.valueOf(S2.minus(S1).mod()).compareTo(EPSILON)
== 0)
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> C1 = 2^(z-1) / (z-1)
* <br> C2 = 2^(z)
* @param z - The complex value of z = a + i*b
* @return the value of Zeta(z) through C1, C2, and the
* quadrature approximation.
*/
public static Complex AbelPlana(Complex z) {
Complex two = new Complex(2.0, 0.0);
Complex C1 = two.pow(z.minus(1.0)).divide(z.minus(1.0));
Complex C2 = two.pow(z);
Complex mult = C2.multiply(adaptiveQuad(0, 10, z));
if ( z.re < 0 && z.re % 2 == 0 && z.im == 0)
return new Complex(0.0, 0.0);
else
return C1.minus(mult);
}
// Needed to reference the super class
public AbelMain(double re, double im) {
super(re, im);
}
}
The BigDecimal EPSILON = new BigDecimal(1E-20); variable works fine and returns values.
Calculation of the Riemann Zeta Function in the form Zeta(s) = a + ib.
Enter the value of [a] inside the Riemann Zeta Function: 54
Enter the value of [b] inside the Riemann Zeta Function: 2
The value for Zeta(s) is 3.742894786684807E13 + 1.6565035537902216E14*i
Total time taken is 0.01 seconds.
Helper program
/**************************************************************************
**
** 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.
**************************************************************************/
public class Complex extends Object{
public double re;
public double im;
/**
Constructor for the complex number z = a + i*b
@param re Real part
@param im Imaginary part
*/
public Complex (double re, double im) {
this.re = re;
this.im = im;
}
/**
Real part of the Complex number
@return Re[z] where z = a + i*b.
*/
public double real() {
return re;
}
/**
Imaginary part of the Complex number
@return Im[z] where z = a + i*b.
*/
public double 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);
}
/**
Addition of two Complex numbers (z is unchanged).
<br>(a+i*b) + (c+i*d) = (a+c) + i*(b+d)
@param t is the number to add.
@return z+t where z = a+i*b and t = c+i*d
*/
public Complex add(Complex t) {
return new Complex(re + t.real(), im + t.imag());
}
/**
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){
return new Complex(this.re + d, this.im);
}
/**
Subtraction of two Complex numbers (z is unchanged).
<br>(a+i*b) - (c+i*d) = (a-c) + i*(b-d)
@param t is the number to subtract.
@return z-t where z = a+i*b and t = c+i*d
*/
public Complex minus(Complex t) {
return new Complex(re - t.real(), im - t.imag());
}
/**
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 minus(double d){
return new Complex(this.re - d, this.im);
}
/**
Complex multiplication (z is unchanged).
<br> (a+i*b) * (c+i*d) = (a*c)+ i(b*c) + i(a*d) - (b*d)
@param t is the number to multiply by.
@return z*t where z = a+i*b and t = c+i*d
*/
public Complex multiply(Complex t) {
return new Complex(re * t.real() - im * t.imag(), re *
t.imag() + im * t.real());
}
/**
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){
return new Complex(this.re * d,this.im * d);}
/**
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 != 0.0 || im != 0.0)
return Math.sqrt(re*re + im*im);
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) {
return Math.sqrt(Math.pow(re*re, 2) + Math.pow(im*im, 2));
}
/**
Division of Complex numbers (doesn't change this Complex number).
<br>(a+i*b) / (c+i*d) = (a+i*b)*(c-i*d) / (c+i*d)*(c-i*d) =
* (a*c+b*d) + i*(b*c-a*d) / (c^2-d^2)
@param t is the number to divide by
@return new Complex number z/t where z = a+i*b
*/
public Complex divide (Complex t) {
double denom = Math.pow(t.mod(), 2); // Square the modulus
return new Complex((re * t.real() + im * t.imag()) / denom,
(im * t.real() - re * t.imag()) / denom);
}
/**
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){
return new Complex(this.re / d, this.im / d);
}
/**
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) * Math.cos(im), Math.exp(re) *
Math.sin(im));
}
/**
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, re);
}
/**
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) * Math.sin(re), sinh(im)*
Math.cos(re));
}
/**
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) * Math.cos(re), -sinh(im) *
Math.sin(re));
}
/**
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) * Math.cos(im), cosh(re) *
Math.sin(im));
}
/**
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) *Math.cos(im), sinh(re) *
Math.sin(im));
}
/**
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,this.im-1.0);
Complex den = new Complex(-this.re,-this.im-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());
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() {
if (re != 0.0 && im > 0.0) {
return re + " + " + im +"*i";
}
if (re !=0.0 && im < 0.0) {
return re + " - "+ (-im) + "*i";
}
if (im == 0.0) {
return String.valueOf(re);
}
if (re == 0.0) {
return im + "*i";
}
return re + " + i*" + im;
}
}
Increasing the error term to BigDecimal EPSILON = new BigDecimal(1E-50); causes the program to return java.lang.StackOverflowError.
/**************************************************************************
**
** 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 AbelMain extends Complex {
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(re, 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^(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, 0).pow(z.divide(2.0));
double D2 = Math.pow(Math.E, Math.PI * t) + 1.0;
Complex den = D1.multiply(D2);
return num.divide(den);
}
/**
* 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) {
BigDecimal EPSILON = new BigDecimal(1E-50);
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(4)).add(f(b, z))).
multiply(step / 6.0);
Complex S2 = (f(a, z).add(f(d, z).multiply(4)).add(f(c, z).multiply(2))
.add(f(e, z).multiply(4)).add(f(b, z))).multiply(step / 12.0);
Complex result = (S2.minus(S1)).divide(15.0);
if(BigDecimal.valueOf(S2.minus(S1).mod()).compareTo(EPSILON) == -1
|| BigDecimal.valueOf(S2.minus(S1).mod()).compareTo(EPSILON)
== 0)
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> C1 = 2^(z-1) / (z-1)
* <br> C2 = 2^(z)
* @param z - The complex value of z = a + i*b
* @return the value of Zeta(z) through C1, C2, and the
* quadrature approximation.
*/
public static Complex AbelPlana(Complex z) {
Complex two = new Complex(2.0, 0.0);
Complex C1 = two.pow(z.minus(1.0)).divide(z.minus(1.0));
Complex C2 = two.pow(z);
Complex mult = C2.multiply(adaptiveQuad(0, 10, z));
if ( z.re < 0 && z.re % 2 == 0 && z.im == 0)
return new Complex(0.0, 0.0);
else
return C1.minus(mult);
}
// Needed to reference the super class
public AbelMain(double re, double im) {
super(re, im);
}
}
Test Run
Calculation of the Riemann Zeta Function in the form Zeta(s) = a + ib.
Enter the value of [a] inside the Riemann Zeta Function: 54
Enter the value of [b] inside the Riemann Zeta Function: 2
Exception in thread "main" java.lang.StackOverflowError
at sun.misc.FpUtils.getExponent(FpUtils.java:147)
at java.lang.Math.getExponent(Math.java:1310)
at java.lang.StrictMath.floorOrCeil(StrictMath.java:355)
at java.lang.StrictMath.floor(StrictMath.java:340)
at java.lang.Math.floor(Math.java:424)
at sun.misc.FloatingDecimal.dtoa(FloatingDecimal.java:620)
at sun.misc.FloatingDecimal.<init>(FloatingDecimal.java:459)
at java.lang.Double.toString(Double.java:196)
at java.math.BigDecimal.valueOf(BigDecimal.java:1069)
at AbelMain.adaptiveQuad(AbelMain.java:92)
at AbelMain.adaptiveQuad(AbelMain.java:97)
at AbelMain.adaptiveQuad(AbelMain.java:97)
at AbelMain.adaptiveQuad(AbelMain.java:97)
at AbelMain.adaptiveQuad(AbelMain.java:97)
at AbelMain.adaptiveQuad(AbelMain.java:97)
at AbelMain.adaptiveQuad(AbelMain.java:97)
at AbelMain.adaptiveQuad(AbelMain.java:97)
at AbelMain.adaptiveQuad(AbelMain.java:97)
at AbelMain.adaptiveQuad(AbelMain.java:97)
at AbelMain.adaptiveQuad(AbelMain.java:97)
at AbelMain.adaptiveQuad(AbelMain.java:97)
at AbelMain.adaptiveQuad(AbelMain.java:97)
at AbelMain.adaptiveQuad(AbelMain.java:97)
at AbelMain.adaptiveQuad(AbelMain.java:97)
at AbelMain.adaptiveQuad(AbelMain.java:97)
at AbelMain.adaptiveQuad(AbelMain.java:97)
at AbelMain.adaptiveQuad(AbelMain.java:97)
at AbelMain.adaptiveQuad(AbelMain.java:97)
These two lines point to
if(BigDecimal.valueOf(S2.minus(S1).mod()).compareTo(EPSILON) == -1
|| BigDecimal.valueOf(S2.minus(S1).mod()).compareTo(EPSILON)
== 0)
and
return adaptiveQuad(a, c, z).add(adaptiveQuad(c, b, z));
Is this error due to infinite recursion? That value for epsilon might have to be above a certain threshold?
You mix "float" and "BigDecimal", with terrific side effects... in particular (I tried with input data 54 and 2), when you use "step" float variable, as result of difference between a an b, after a while it reaches a "fixed point" of 1.0587911840678754E-22, due to the internal representation of Complex data, where you use floats.
Well, at this point, this "step" never change, S1 and S2 neither, and difference remains the same ("result" evaluates 0.0, -2.3915493791143544E-44) so it never goes under EPSILON (as 1E-50) and algorithm recurses infinitely (or better, until explodes). In other words, you avoided approximations from a side, using BigDecimal, but you reintroduced it via float usage :)
Try to remove float usage everywhere (in Complex too), and everything should run fine.