Search code examples
javarecursionnumber-theory

Riemann Zeta Function in Java - Infinite Recursion with Functional Form


Note: Updated on 06/17/2015. Of course this is possible. See the solution below.

Even if anyone copies and pastes this code, you still have a lot of cleanup to do. Also note that you will have problems inside the critical strip from Re(s) = 0 to Re(s) = 1 :). But this is a good start.

import java.util.Scanner;

public class NewTest{

public static void main(String[] args) {
    RiemannZetaMain func = new RiemannZetaMain();
    double s = 0;
    double start, stop, totalTime;
    Scanner scan = new Scanner(System.in);
    System.out.print("Enter the value of s inside the Riemann Zeta Function: ");
    try {
            s = scan.nextDouble();
    } 
    catch (Exception e) {
        System.out.println("You must enter a positive integer greater than 1.");
    }
    start = System.currentTimeMillis();
    if (s <= 0)
        System.out.println("Value for the Zeta Function = " + riemannFuncForm(s));
    else if (s == 1)
        System.out.println("The zeta funxtion is undefined for Re(s) = 1.");
    else if(s >= 2)
        System.out.println("Value for the Zeta Function = " + getStandardSum(s));
    else
        System.out.println("Value for the Zeta Function = " + getNewSum(s));
    stop = System.currentTimeMillis();
    totalTime = (double) (stop-start) / 1000.0;
    System.out.println("Total time taken is " + totalTime + " seconds.");
}

// Standard form the the Zeta function.
public static double standardZeta(double s) {
    int n = 1;
    double currentSum = 0;
    double relativeError = 1;
    double error = 0.000001;
    double remainder;

    while (relativeError > error) {
        currentSum = Math.pow(n, -s) + currentSum;
        remainder = 1 / ((s-1)* Math.pow(n, (s-1)));
        relativeError =  remainder / currentSum;
        n++;
    }
    System.out.println("The number of terms summed was " + n + ".");
    return currentSum;
}

public static double getStandardSum(double s){
    return standardZeta(s);
}

//New Form
// zeta(s) = 2^(-1+2 s)/((-2+2^s) Gamma(1+s)) integral_0^infinity t^s sech^2(t) dt  for Re(s)>-1
public static double Integrate(double start, double end) {
    double currentIntegralValue = 0;
    double dx = 0.0001d; // The size of delta x in the approximation
    double x = start; // A = starting point of integration, B = ending point of integration.

    // Ending conditions for the while loop
    // Condition #1: The value of b - x(i) is less than delta(x).
    // This would throw an out of bounds exception.
    // Condition #2: The value of b - x(i) is greater than 0 (Since you start at A and split the integral 
    // up into "infinitesimally small" chunks up until you reach delta(x)*n.
    while (Math.abs(end - x) >= dx && (end - x) > 0) {
        currentIntegralValue += function(x) * dx; // Use the (Riemann) rectangle sums at xi to compute width * height
        x += dx; // Add these sums together
    }
    return currentIntegralValue;
}

private static double function(double s) {
    double sech = 1 / Math.cosh(s); // Hyperbolic cosecant
    double squared = Math.pow(sech, 2);
    return ((Math.pow(s, 0.5)) * squared);
}

public static double getNewSum(double s){
double constant = Math.pow(2, (2*s)-1) / (((Math.pow(2, s)) -2)*(gamma(1+s)));
    return constant*Integrate(0, 1000);
}

// Gamma Function - Lanczos approximation
public static double gamma(double s){
                double[] p = {0.99999999999980993, 676.5203681218851, -1259.1392167224028,
                                  771.32342877765313, -176.61502916214059, 12.507343278686905,
                                  -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7};
                int g = 7;
                if(s < 0.5) return Math.PI / (Math.sin(Math.PI * s)*gamma(1-s));

                s -= 1;
                double a = p[0];
                double t = s+g+0.5;
                for(int i = 1; i < p.length; i++){
                        a += p[i]/(s+i);
                }

                return Math.sqrt(2*Math.PI)*Math.pow(t, s+0.5)*Math.exp(-t)*a;
        }

//Binomial Co-efficient - NOT CURRENTLY USING
/*
public static double binomial(int n, int k)
{
    if (k>n-k)
        k=n-k;

    long b=1;
    for (int i=1, m=n; i<=k; i++, m--)
        b=b*m/i;
    return b;
} */

// Riemann's Functional Equation
// Tried this initially and utterly failed.
public static double riemannFuncForm(double s) {
double term = Math.pow(2, s)*Math.pow(Math.PI, s-1)*(Math.sin((Math.PI*s)/2))*gamma(1-s);
double nextTerm = Math.pow(2, (1-s))*Math.pow(Math.PI, (1-s)-1)*(Math.sin((Math.PI*(1-s))/2))*gamma(1-(1-s));
double error = Math.abs(term - nextTerm);

if(s == 1.0)
    return 0;
else 
    return Math.pow(2, s)*Math.pow(Math.PI, s-1)*(Math.sin((Math.PI*s)/2))*gamma(1-s)*standardZeta(1-s);
}

}


Solution

  • I think I need to use a different form of the zeta function. When I run the entire program ---

    import java.util.Scanner;
    
    public class Test4{
    
    public static void main(String[] args) {
        RiemannZetaMain func = new RiemannZetaMain();
        double s = 0;
        double start, stop, totalTime;
        Scanner scan = new Scanner(System.in);
        System.out.print("Enter the value of s inside the Riemann Zeta Function: ");
        try {
                s = scan.nextDouble();
        } 
        catch (Exception e) {
            System.out.println("You must enter a positive integer greater than 1.");
        }
        start = System.currentTimeMillis();
        if(s >= 2)
            System.out.println("Value for the Zeta Function = " + getStandardSum(s));
        else
            System.out.println("Value for the Zeta Function = " + getRiemannSum(s));
        stop = System.currentTimeMillis();
        totalTime = (double) (stop-start) / 1000.0;
        System.out.println("Total time taken is " + totalTime + " seconds.");
    }
    
    // Standard form the the Zeta function.
    public static double standardZeta(double s) {
        int n = 1;
        double currentSum = 0;
        double relativeError = 1;
        double error = 0.000001;
        double remainder;
    
        while (relativeError > error) {
            currentSum = Math.pow(n, -s) + currentSum;
            remainder = 1 / ((s-1)* Math.pow(n, (s-1)));
            relativeError =  remainder / currentSum;
            n++;
        }
        System.out.println("The number of terms summed was " + n + ".");
        return currentSum;
    }
    
    public static double getStandardSum(double s){
        return standardZeta(s);
    }
    
    // Riemann's Functional Equation
    public static double riemannFuncForm(double s, double threshold, double currSteps, int maxSteps) {
    double term = Math.pow(2, s)*Math.pow(Math.PI, s-1)*(Math.sin((Math.PI*s)/2))*gamma(1-s);
    //double nextTerm = Math.pow(2, (1-s))*Math.pow(Math.PI, (1-s)-1)*(Math.sin((Math.PI*(1-s))/2))*gamma(1-(1-s));
    //double error = Math.abs(term - nextTerm);
    
    if(s == 1.0)
        return 0;
    else if (s == 0.0)
        return -0.5;
    else if (term < threshold) {//The recursion will stop once the term is less than the threshold
        System.out.println("The number of steps is " + currSteps);
        return term;
    }
    else if (currSteps == maxSteps) {//The recursion will stop if you meet the max steps
        System.out.println("The series did not converge.");
        return term;
    }    
    else //Otherwise just keep recursing
        return term*riemannFuncForm(1-s, threshold, ++currSteps, maxSteps);
    }
    
    public static double getRiemannSum(double s) {
        double threshold = 0.00001;
        double currSteps = 1;
        int maxSteps = 1000;
        return riemannFuncForm(s, threshold, currSteps, maxSteps);
    }
    
    // Gamma Function - Lanczos approximation
    public static double gamma(double s){
                    double[] p = {0.99999999999980993, 676.5203681218851, -1259.1392167224028,
                                      771.32342877765313, -176.61502916214059, 12.507343278686905,
                                      -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7};
                    int g = 7;
                    if(s < 0.5) return Math.PI / (Math.sin(Math.PI * s)*gamma(1-s));
    
                    s -= 1;
                    double a = p[0];
                    double t = s+g+0.5;
                    for(int i = 1; i < p.length; i++){
                            a += p[i]/(s+i);
                    }
    
                    return Math.sqrt(2*Math.PI)*Math.pow(t, s+0.5)*Math.exp(-t)*a;
            }
    
    //Binomial Co-efficient
    public static double binomial(int n, int k)
    {
        if (k>n-k)
            k=n-k;
    
        long b=1;
        for (int i=1, m=n; i<=k; i++, m--)
            b=b*m/i;
        return b;
    }
    

    }

    I notice that plugging in zeta(-1) returns -

    Enter the value of s inside the Riemann Zeta Function: -1
    The number of steps is 1.0
    Value for the Zeta Function = -0.0506605918211689
    Total time taken is 0.0 seconds.
    

    I knew that this value was -1/12. I checked some other values with wolfram alpha and observed that --

    double term = Math.pow(2, s)*Math.pow(Math.PI, s-1)*(Math.sin((Math.PI*s)/2))*gamma(1-s);
    

    Returns the correct value. It is just that I am multiplying this value every time by zeta(1-s). In the case of Zeta(1/2), this will always multiply the result by 0.99999999.

    Enter the value of s inside the Riemann Zeta Function: 0.5
    The series did not converge.
    Value for the Zeta Function = 0.999999999999889
    Total time taken is 0.006 seconds.
    

    I am going to see if I can replace the part for --

        else if (term < threshold) {//The recursion will stop once the term is less than the threshold
        System.out.println("The number of steps is " + currSteps);
        return term;
    }
    

    This difference is the error between two terms in the summation. I may not be thinking about this correctly, it is 1:16am right now. Let me see if I can think better tomorrow ....