Riemann Zeta function implementation template #2904

Bobingstern · 2023-02-23T02:28:28Z

I've coded up an implementation to calculate the riemann zeta function over the entire complex plane using math.js. It works very quickly for real numbers and uses a number of terms proportional to $\approx (1.3 * \ln (d * 0.1) + 0.9 * | t |)$ where $d$ is the number of decimal places of accuracy and $t$ is the imaginary part of the input. I've implemented this in my small js library called special.js. Here is a description of the algorithm used:

Computation makes use of the function $d_{k} = n \sum_{j = k}^{n} \frac{(n + j - 1)! 4^{j}}{(n - j)! (2 j)!}$ where $n$ is chosen arbitrarily making use of the precision and uses $ζ (s) = \frac{1}{d_{0} (1 - 2^{1 - s})} \sum_{k = 1}^{n} \frac{(- 1)^{k - 1} d_{k}}{k^{s}}$ for $ℜ (s) > 1$ and uses the functional equation $ζ (s) = 2^{s} π^{s - 1} \sin (\frac{π s}{2}) Γ (1 - s) ζ (1 - s)$ where $Γ (s)$ is the gamma function to extend it to the entire complex plane.

The js code in math.js is as follows:

const math = require("mathjs")

function zeta(s, prec, only){
	s = math.complex(s)
    if (s.re == 0 && s.im == 0){
        return -0.5
    }
    prec = prec == undefined ? 1e-3 : prec
    let dec = -Math.round(Math.log(prec*0.1)/Math.log(10))
    let n = Math.round(1.3*dec + 0.9*Math.abs(s.im))
    n = Math.min(n, 60)
    
    function d(k){
        let S = 0
        for (let j=k;j<=n;j++){
            S += (math.factorial(n+j-1) * (4**j) / (math.factorial(n-j)*math.factorial(2*j)  ))
        }
        return n*S
    }
    function f(s){
        let c = math.divide(1, 
            math.multiply(d(0), math.subtract(1, math.pow(2, math.subtract(1, s)))))
        let S = math.complex(0, 0)
        for (let k=1;k<=n;k++){
            S = S.add(
                math.divide((-1)**(k-1) * d(k), math.pow(k, s))
            )
        }
        return math.multiply(c, S)
    }
    if (s.re > 1 || s.re == 0 || only){
        //console.log("Running "+ n +" iterations")
        return f(s)
    }
    else{
        let c = math.multiply(math.pow(2, s), math.pow(Math.PI, math.subtract(s, 1)))
        c = math.multiply(c, (math.sin(math.multiply(Math.PI/2, s))))
        c = math.multiply(c, math.gamma(math.subtract(1, s)))
        //Recursion
        return math.multiply(c, zeta(math.subtract(1, s), prec, true))
    }
}

function benchmark(s, n){
	let avg = 0
	for (let i=0;i<n;i++){
		let now = performance.now()
		zeta(s, 0.0000001)
		avg += performance.now() - now
	}
	return avg/n
}
console.log(benchmark(math.complex(1/2, 14.134725), 100)) //First zero of the zeta function

I'm not very experienced in node.js libraries and such so I hope that someone more experienced can add it to the library.

The text was updated successfully, but these errors were encountered:

josdejong · 2023-02-23T15:20:28Z

Thanks for sharing Anik! Would be nice indeed to implement the Riemann Zeta function in mathjs itself.

Anyone able to work out the function in a PR and add tests and documentation? Help would be welcome.

Bobingstern · 2023-02-23T20:44:50Z

I've added a new enhancement using the Riemann Siegel formula to evaluate on the critical line much faster. The code itself is very ugly since the remainder term calculation is horrible and so I used a taylor series for the derivatives and such. The graph is that of the zeta function on the critical line with varying imaginary inputs with blue and orange representing the real and imaginary output respectively. You can read about the Riemann Siegel Z(t) function here. The new code is slightly changed and is as follows. I'll try my best to figure out how to incorporate this into math.js and try to send a PR with tests and docs.

const math = require("mathjs")

function zeta(s, prec, only){
	s = math.complex(s)
    if (s.re == 0 && s.im == 0){
        return -0.5
    }
    prec = prec == undefined ? 1e-3 : prec
    let dec = -Math.round(Math.log(prec*0.1)/Math.log(10))
    let n = Math.round(1.3*dec + 0.9*Math.abs(s.im))
    //n = Math.min(n, 60)
    if (s.re == 0.5 && s.im >= 1){
        return math.divide(RiemannZ(s.im, 5), math.exp(math.complex(0, theta(s.im))))
    }
    function d(k){
        let S = 0
        for (let j=k;j<=n;j++){
            S += (math.factorial(n+j-1) * (4**j) / (math.factorial(n-j)*math.factorial(2*j)  ))
        }
        return n*S
    }
    function f(s){
        let c = math.divide(1, 
            math.multiply(d(0), math.subtract(1, math.pow(2, math.subtract(1, s)))))
        let S = math.complex(0, 0)
        for (let k=1;k<=n;k++){
            S = S.add(
                math.divide((-1)**(k-1) * d(k), math.pow(k, s))
            )
        }
        return math.multiply(c, S)
    }
    if (s.re > 1 || s.re == 0 || only){
        //console.log("Running "+ n +" iterations")
        return f(s)
    }
    else{
        let c = math.multiply(math.pow(2, s), math.pow(Math.PI, math.subtract(s, 1)))
        c = math.multiply(c, (math.sin(math.multiply(Math.PI/2, s))))
        c = math.multiply(c, math.gamma(math.subtract(1, s)))
        //Recursion
        return math.multiply(c, zeta(math.subtract(1, s), prec, true))
    }
}

//Hellish coefficents function. It is minifed since its just a bunch of return statements
function C($,_){return 0==$?.3826834323650898*_**0+.43724046807752043*_**2+.1323765754803435*_**4-.013605026047674188*_**6-.013567621970103581*_**8-.0016237253231444653*_**10+2970535373337969e-19*_**12+794330087952147e-19*_**14+46556124614505e-20*_**16-143272516309551e-20*_**18-10354847112313e-20*_**20+1235792708386e-20*_**22+17881083858e-19*_**24-339141439e-19*_**26-163266339e-19*_**28-37851093e-20*_**30+9327423e-20*_**32+522184e-20*_**34-33507e-20*_**36-3412e-20*_**38+58e-20*_**40+15e-20*_**42:1==$?-.026825102628375348*_**1+.013784773426351853*_**3+.03849125048223508*_**5+.009871066299062077*_**7-.0033107597608584044*_**9-.0014647808577954152*_**11-1320794062487696e-20*_**13+5922748701847141e-20*_**15+598024258537345e-20*_**17-96413224561698e-20*_**19-18334733722714e-20*_**21+446708756272e-20*_**23+270963508218e-20*_**25+7785288654e-20*_**27-2343762601e-20*_**29-158301728e-20*_**31+12119942e-20*_**33+1458378e-20*_**35-28786e-20*_**37-8663e-20*_**39-84e-20*_**41+36e-20*_**43+1e-20*_**45:2==$?.005188542830293168*_**0+30946583880634744e-20*_**2-.011335941078229373*_**4+.0022330457419581446*_**6+.00519663740886233*_**8+3439914407620834e-19*_**10-5910648427470583e-19*_**12-10229972547935857e-20*_**14+2088839221699276e-20*_**16+592766549309654e-20*_**18-16423838362436e-20*_**20-15161199700941e-20*_**22-590780369821e-20*_**24+209115148595e-20*_**26+17815649583e-20*_**28-1616407246e-20*_**30-238069625e-20*_**32+5398265e-20*_**34+1975014e-20*_**36+23333e-20*_**38-11188e-20*_**40-416e-20*_**42+44e-20*_**44+3e-20*_**46:3==$?-.0013397160907194568*_**1+.003744215136379394*_**3-.0013303178919321468*_**5-.0022654660765471786*_**7+9548499998506731e-19*_**9+6010038458963604e-19*_**11-10128858286776622e-20*_**13-6865733449299826e-20*_**15+59853667915386e-20*_**17+333165985123995e-20*_**19+21919289102435e-20*_**21-7890884245681e-20*_**23-94146850813e-19*_**25+95701162109e-20*_**27+18763137453e-20*_**29-443783768e-20*_**31-224267385e-20*_**33-3627687e-20*_**35+1763981e-20*_**37+79608e-20*_**39-942e-19*_**41-713e-20*_**43+33e-20*_**45+4e-20*_**47:46483389361763383e-20*_**0-.001005660736534047*_**2+24044856573725794e-20*_**4+.0010283086149702322*_**6-7657861071755644e-19*_**8-20365286803084818e-20*_**10+2321229049106873e-19*_**12+326021442438652e-19*_**14-2557906251794953e-20*_**16-410746443891574e-20*_**18+117811136403713e-20*_**20+24456561422485e-20*_**22-2391582476734e-20*_**24-750521420704e-20*_**26+13312279416e-20*_**28+13440626754e-20*_**30+351377004e-20*_**32-151915445e-20*_**34-8915418e-20*_**36+1119589e-20*_**38+10516e-19*_**40-5179e-20*_**42-807e-20*_**44+11e-20*_**46+4e-20*_**48}

function theta(t){
    twopi = 2 * math.pi
    return (t/2.0 * math.log(t/twopi) - t/2.0 - math.pi/8.0
            + 1.0/(48.0*t) + 7.0/(5760*t**3))
}
//Riemann Z(t) formula
function RiemannZ(t, r){
    tDivTwoPi = t / (2 * math.pi)
    val = math.sqrt(tDivTwoPi)
    N = math.floor(math.abs(val))
    SUM = 0

    i = 1
    while (i<=N){
        SUM += (math.cos(theta(t) - t * math.log(i))) / math.sqrt(i)
        i += 1
    }
    SUM *= 2
    frac = val - N
    k = 0
    R = 0

    //Remainder Calculation
    while (k <= r){
        R += C(k, 2.0*frac-1.0) * tDivTwoPi**(k * (-0.5))
        k += 1
    }
    remainder = (-1)**(N-1) * tDivTwoPi**(-0.25) * R
    return SUM + remainder
}
function benchmark(s, n){
	let avg = 0
	for (let i=0;i<n;i++){
		let now = performance.now()
		zeta(s, 0.0000001)
		avg += performance.now() - now
	}
	return avg/n
}

console.log(benchmark(math.complex(1, 167.235364), 100)) //Normal zeta function testing

console.log(benchmark(math.complex(1/2, 167.235364), 100)) //Riemann Siegel Critical zeta function testing (Triggered if Re(s) is 1/2). Much Much faster

HanchaiN · 2023-03-23T10:12:14Z

Maybe we can make the hellish coefficients function more readable by storing the whole coefficient like in erf.js. This may cost more space and time since it will have to either store unused coefficient (0s) or figures out if you will use odd or even coefficients. The latter might be better in terms of performance though.
Here's an example (the coefficient is not included but will be referred to as c):

function C(n, x){
    const x2 = x * x;
    const m = n < 4 ? n : 4;
    const y = c[m].reduceRight((acc, v) => acc * x2 + v, 0);
    if (m % 2 === 1) return x * y;
    return y;
}

It should be noted that RiemannZ(t, r) seems to be called only with RiemannZ(s.im, 5) so the while loop in the remainder calculation (while (k <= r) {...}) might as well be replaced as a for-loop. Also, the function C only provides cases for n from 1 to 4 or more while the function will receive n from 0 to 5.

Bobingstern · 2023-03-23T16:13:05Z

I'm thinking about replacing the coefficients function with the actual exact function since currently it is a Taylor series and will not be accurate everywhere. The actual expressions is, unsurprisingly, extremely complicated and requires high order derivatives of Psi(x) defined by cos(2pi(x^2-x-1/16))/cos(2pi x). It can be calculated beforehand but the expression is massive. I'll try to figure a better way to calculate them

Bobingstern · 2023-03-24T12:34:25Z

Actually, scratch that. The remainder terms coefficients are only computed from values between 0 and 1 so the taylor polynomial will always be accurate

Bobingstern · 2023-03-25T20:14:50Z

Here is a new version that uses an array of Chebyshev coefficients that I recalculated since they are easier and more efficient to calculate as they dont need derivatives of the function. It should be slightly more accurate. I also increased the size of the theta function asymptotic series. If this ever becomes a part of the library make the documentation says that on the critical strip Re(s)=1/2 it is more accurate for larger imaginary parts of s.

const math = require("mathjs")

function zeta(s, prec, only){
	s = math.complex(s)
    if (s.re == 0 && s.im == 0){
        return -0.5
    }
    prec = prec == undefined ? 1e-3 : prec
    let dec = -Math.round(Math.log(prec*0.1)/Math.log(10))
    let n = Math.round(1.3*dec + 0.9*Math.abs(s.im))
    //n = Math.min(n, 60)
    if (s.re == 0.5 && s.im >= 1){
        return math.divide(RiemannZ(s.im, 5), math.exp(math.complex(0, theta(s.im))))
    }
    function d(k){
        let S = 0
        for (let j=k;j<=n;j++){
            S += (math.factorial(n+j-1) * (4**j) / (math.factorial(n-j)*math.factorial(2*j)  ))
        }
        return n*S
    }
    function f(s){
        let c = math.divide(1, 
            math.multiply(d(0), math.subtract(1, math.pow(2, math.subtract(1, s)))))
        let S = math.complex(0, 0)
        for (let k=1;k<=n;k++){
            S = S.add(
                math.divide((-1)**(k-1) * d(k), math.pow(k, s))
            )
        }
        return math.multiply(c, S)
    }
    if (s.re > 1 || s.re == 0 || only){
        //console.log("Running "+ n +" iterations")
        return f(s)
    }
    else{
        let c = math.multiply(math.pow(2, s), math.pow(Math.PI, math.subtract(s, 1)))
        c = math.multiply(c, (math.sin(math.multiply(Math.PI/2, s))))
        c = math.multiply(c, math.gamma(math.subtract(1, s)))
        //Recursion
        return math.multiply(c, zeta(math.subtract(1, s), prec, true))
    }
}

//Recalculated Chebyshev Coefficients of the C_n remainder terms

let COEFFS = [
    [0.923879532511287, -2.40447091953738, 2.40447091953731, 4.83173297425403, -18.2366510002884, 27.7681391978478, -15.8207845139374, -20.3367348174362, 59.9961770203178, -68.9551740675586, 31.2289762240246, 32.5922407636694, -78.9517947267673, 76.1241393720219, -29.3541472364891, -26.4278439633590, 55.6592159924519, -47.1547592849097, 16.0946292814341, 12.9266552804525, -24.4211331642284, 18.6477145049252, -5.73684112759751, -4.23325603983216, 7.21979739242415, -4.99805538450325, 1.44686340611579, 0.799581563264023, -1.31257324199806, 0.911170990271288, -0.412224988085833, 0.128698140636259, -0.0270029163232092, 0.00344714227880782, -0.000202773075223989], 
    [-0.0305973064996957, 0.461939766255362, -1.75843795691946, 2.00372576637813, 4.50743571933152, -21.2760928439984, 36.6797542376395, -23.7311762621744, -34.0548079959763, 109.992977770541, -137.869658800580, 67.6629804489852, 76.1471060202849, -197.381017573592, 203.055481067454, -83.1670713430562, -79.3213718207383, 176.299440913682, -157.057082049896, 55.9466150913115, 47.4759178652798, -92.5370798084723, 72.7962775735004, -23.8221608424208, -14.8104724965542, 27.2300044031592, -21.0837908437584, 10.5984076439006, -3.66338228115894, 0.848149907535393, -0.119100570732473, 0.00768390778919178, 0, 0, 0], 
    [0.00126887416428184, -0.0111913805347592, 0.327207334757912, -1.44355359345758, 1.87849284960183, 4.68609783021407, -23.9356008766487, 44.5242414810130, -30.9000728632431, -47.3669886936200, 162.241028400620, -215.226813669303, 111.357648674754, 131.860162876894, -358.437979444606, 385.481292937290, -163.854692271498, -164.368569789547, 374.996985792918, -342.561372714489, 129.386633986255, 91.5227366906853, -191.992073538253, 168.821481501502, -96.6670956058399, 38.8879995294781, -11.3837276978800, 2.70237445867694, -0.677132748432874, 0.189485290404773, -0.0424074953767696, 0.00576293084189383, -0.000360183177618364, 0, 0], 
    [-0.000198685215762818, 0.000634437558714868, -0.0133989039223216, 0.272672627933153, -1.30262308602626, 1.80420030669118, 4.79983592017692, -25.9749113933860, 51.0616883030611, -37.2822241555971, -60.0160104336138, 215.160016643964, -297.689882486609, 159.317039674423, 198.384065753389, -551.165445390966, 605.577082784030, -273.442623542416, -221.155939212526, 546.948246506788, -564.287770461232, 389.025949390999, -205.040274582806, 94.9650486292110, -43.9691141455737, 19.5426422565593, -7.04944653893387, 1.80267548660601, -0.304299103221043, 0.0376794554422968, -0.00631617634682575, 0.00136798372183128, -0.000180091588809182, 1.09146417460110e-5, 0], 
    [-5.07834970801379e-6, 0.000511211112458644, 0.000462415235288165, -0.00893129701989592, 0.237505483503315, -1.19643463310244, 1.74132246219081, 4.86185383715175, -27.5544279242915, 56.3935328904713, -42.3593283545426, -72.1759566071907, 262.119631415809, -368.201606160634, 211.927733277448, 191.409048782631, -596.068295951527, 769.456754995061, -714.178102108556, 581.496539107120, -455.489248129327, 324.152279512017, -187.349666282090, 82.0787824396755, -27.3285528783793, 7.86174336264113, -2.43869230997746, 0.771987683027644, -0.186001845301658, 0.0275901015701257, -0.00211756194455121, 0.000101873812045577, -2.13747456536137e-5, 2.72866043650276e-6, -1.60509437441339e-7], 
    [-7.39673470098183e-5, -2.64032272925124e-6, -0.000187509769015465, 0.000392835654663165, -0.00817869022817752, 0.214165277226668, -1.09885504050019, 1.59128848777676, 4.93593360797150, -27.1871636716867, 55.3861362588388, -47.4485864374421, -35.9730875561453, 181.701798013587, -356.942344773391, 579.622529836452, -855.333847937984, 1055.12686998679, -1009.05886184206, 732.568881108251, -417.683070958224, 208.229753023005, -103.326483363813, 49.6379792638476, -19.5926031763341, 5.62999428113937, -1.17615289381254, 0.230185499695570, -0.0579965982595923, 0.0134274300840040, -0.00188553293735814, 0.000121647286281170, 0, 0, 0], 
    [1.86387963190633e-6, 4.26943314699186e-5, -4.41281106403069e-5, -0.000294472792644385, 0.00382003543773591, -0.00991822657860930, 0.0527134501127473, -0.338480513590947, 1.65830840099570, -6.12839795343795, 10.0862751985791, 26.6537161100196, -201.593649430621, 580.364809726527, -1036.58913268528, 1297.24311797875, -1241.26229478129, 1018.07530395798, -798.346157613283, 592.123797543887, -371.487874584203, 181.605499101827, -70.0979821459174, 24.7438297827556, -9.33859794630445, 3.32579133488958, -0.880073238034347, 0.151978087934585, -0.0185196703765880, 0.00299982404790995, -0.000671371504200198, 9.12354647108777e-5, -5.70221654442985e-6, 0, 0]
]


function C(n, x){
    const m = n <= 5 ? n : 5;
    const y = COEFFS[m].reduceRight((acc, v) => acc * x + v, 0);
    return y;
}

//function C($,_){return 0==$?.3826834323650898*_**0+.43724046807752043*_**2+.1323765754803435*_**4-.013605026047674188*_**6-.013567621970103581*_**8-.0016237253231444653*_**10+2970535373337969e-19*_**12+794330087952147e-19*_**14+46556124614505e-20*_**16-143272516309551e-20*_**18-10354847112313e-20*_**20+1235792708386e-20*_**22+17881083858e-19*_**24-339141439e-19*_**26-163266339e-19*_**28-37851093e-20*_**30+9327423e-20*_**32+522184e-20*_**34-33507e-20*_**36-3412e-20*_**38+58e-20*_**40+15e-20*_**42:1==$?-.026825102628375348*_**1+.013784773426351853*_**3+.03849125048223508*_**5+.009871066299062077*_**7-.0033107597608584044*_**9-.0014647808577954152*_**11-1320794062487696e-20*_**13+5922748701847141e-20*_**15+598024258537345e-20*_**17-96413224561698e-20*_**19-18334733722714e-20*_**21+446708756272e-20*_**23+270963508218e-20*_**25+7785288654e-20*_**27-2343762601e-20*_**29-158301728e-20*_**31+12119942e-20*_**33+1458378e-20*_**35-28786e-20*_**37-8663e-20*_**39-84e-20*_**41+36e-20*_**43+1e-20*_**45:2==$?.005188542830293168*_**0+30946583880634744e-20*_**2-.011335941078229373*_**4+.0022330457419581446*_**6+.00519663740886233*_**8+3439914407620834e-19*_**10-5910648427470583e-19*_**12-10229972547935857e-20*_**14+2088839221699276e-20*_**16+592766549309654e-20*_**18-16423838362436e-20*_**20-15161199700941e-20*_**22-590780369821e-20*_**24+209115148595e-20*_**26+17815649583e-20*_**28-1616407246e-20*_**30-238069625e-20*_**32+5398265e-20*_**34+1975014e-20*_**36+23333e-20*_**38-11188e-20*_**40-416e-20*_**42+44e-20*_**44+3e-20*_**46:3==$?-.0013397160907194568*_**1+.003744215136379394*_**3-.0013303178919321468*_**5-.0022654660765471786*_**7+9548499998506731e-19*_**9+6010038458963604e-19*_**11-10128858286776622e-20*_**13-6865733449299826e-20*_**15+59853667915386e-20*_**17+333165985123995e-20*_**19+21919289102435e-20*_**21-7890884245681e-20*_**23-94146850813e-19*_**25+95701162109e-20*_**27+18763137453e-20*_**29-443783768e-20*_**31-224267385e-20*_**33-3627687e-20*_**35+1763981e-20*_**37+79608e-20*_**39-942e-19*_**41-713e-20*_**43+33e-20*_**45+4e-20*_**47:46483389361763383e-20*_**0-.001005660736534047*_**2+24044856573725794e-20*_**4+.0010283086149702322*_**6-7657861071755644e-19*_**8-20365286803084818e-20*_**10+2321229049106873e-19*_**12+326021442438652e-19*_**14-2557906251794953e-20*_**16-410746443891574e-20*_**18+117811136403713e-20*_**20+24456561422485e-20*_**22-2391582476734e-20*_**24-750521420704e-20*_**26+13312279416e-20*_**28+13440626754e-20*_**30+351377004e-20*_**32-151915445e-20*_**34-8915418e-20*_**36+1119589e-20*_**38+10516e-19*_**40-5179e-20*_**42-807e-20*_**44+11e-20*_**46+4e-20*_**48}

function theta(t){
    twopi = 2 * math.pi
    return (t/2.0 * math.log(t/twopi) - t/2.0 - math.pi/8.0
            + 1.0/(48.0*t) + 7.0/(5760*t**3) + 31/(80640*t**5) + 
            127/(430080*t**7) + 511/(1216512*t**9) + 1414477/(1476034560*t**11) + 
            8191/(2555904*t**13) + 118518239/(8021606400*t**15) + 5749691557/(64012419072*t**17) + 
            91546277357/(131491430400*t**19) + 23273283019/(3472883712*t**21))
}
//Riemann Z(t) formula
function RiemannZ(t, r){
    let tDivTwoPi = t / (2 * math.pi)
    let val = math.sqrt(tDivTwoPi)
    let N = math.floor(math.abs(val))
    let SUM = 0

    let i = 1
    while (i<=N){
        SUM += (math.cos(theta(t) - t * math.log(i))) / math.sqrt(i)
        i += 1
    }
    SUM *= 2
    let frac = val - N
    let k = 0
    let R = 0

    //Remainder Calculation
    while (k <= r){
        R += C(k, frac) * tDivTwoPi**(k * (-0.5))
        k += 1
    }
    let remainder = (-1)**(N-1) * tDivTwoPi**(-0.25) * R
    return SUM + remainder
}
function benchmark(s, n){
	let avg = 0
	for (let i=0;i<n;i++){
		let now = performance.now()
		zeta(s, 0.0000001)
		avg += performance.now() - now
	}
	return avg/n
}

//console.log(zeta(math.complex(1/2, 874)))
console.log(benchmark(math.complex(1, 167.235364), 100)) //Normal zeta function testing

console.log(benchmark(math.complex(1/2, 167.235364), 100)) //Riemann Siegel Critical zeta function testing (Triggered if Re(s) is 1/2). Much Much faster

josdejong · 2023-03-27T07:47:18Z

Thanks for sharing @Bobingstern . The large array with coefficients takes a couple of kilobytes, which can be an downside when using the library in the browser. Is this really worth it in your opinion? I understand this will trade kilobytes for runtime performance.

Bobingstern · 2023-03-27T12:17:41Z

Theres no good way I can think of to truncate the array but I think theres a slightly better method for storing it. I used sympy to compute the expressions symoblically. What I could do instead is store only the coefficient of the Chebyshev POLYNOMIAL rather than have sympy expand the entire thing into one fat polynomial. Unforutnately, I will need to figure out a fast and efficient way to compute high derivatives of the chebyshev polynomail buts I could just use the trig definitions of them to make it easier

josdejong · 2023-03-29T10:48:02Z

Ok. I'm not sure what would be "best" here, just trying to understand the need and importance of the large array with pre-computed coefficients.

Bobingstern · 2023-04-03T17:42:02Z

They are important in the formulation of remainder terms in the Riemann-Siegel formula. The remainders are like this
$C_{0} = Ψ (p)$
$C_{1} = - \frac{1}{96 π^{2}} Ψ^{(3)} (p)$
$C_{2} = \frac{1}{18432 π^{4}} Ψ^{(6)} (p) + \frac{1}{64 π^{2}} Ψ^{(2)} (p)$
where
$Ψ (p) = \frac{\cos (2 π (p^{2} - p - 1 / 16))}{\cos (2 π p)}$
Its clear that the high derivatives of $Ψ (p)$ are very problematic. I used a Chebyshev approximation of this function for each remainder term $C_{n}$ . This is why the array is so big. The remainder terms are very important as they added around 5-6 digits of accuracy and even more for larger values of t.

josdejong · 2023-04-04T06:34:11Z

👍 thanks for the clarification.

josdejong · 2023-08-23T14:09:52Z

Function zeta is published now in v11.10.0, thanks again Anik.

josdejong added feature help wanted labels Feb 23, 2023

Bobingstern mentioned this issue May 11, 2023

Riemann Zeta Function Implemented #2950

Closed

josdejong closed this as completed Aug 23, 2023

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Riemann Zeta function implementation template #2904

Riemann Zeta function implementation template #2904

Bobingstern commented Feb 23, 2023

josdejong commented Feb 23, 2023

Bobingstern commented Feb 23, 2023

HanchaiN commented Mar 23, 2023

Bobingstern commented Mar 23, 2023

Bobingstern commented Mar 24, 2023

Bobingstern commented Mar 25, 2023 •

edited

Loading

josdejong commented Mar 27, 2023

Bobingstern commented Mar 27, 2023

josdejong commented Mar 29, 2023

Bobingstern commented Apr 3, 2023

josdejong commented Apr 4, 2023

josdejong commented Aug 23, 2023

Riemann Zeta function implementation template #2904

Riemann Zeta function implementation template #2904

Comments

Bobingstern commented Feb 23, 2023

josdejong commented Feb 23, 2023

Bobingstern commented Feb 23, 2023

HanchaiN commented Mar 23, 2023

Bobingstern commented Mar 23, 2023

Bobingstern commented Mar 24, 2023

Bobingstern commented Mar 25, 2023 • edited Loading

josdejong commented Mar 27, 2023

Bobingstern commented Mar 27, 2023

josdejong commented Mar 29, 2023

Bobingstern commented Apr 3, 2023

josdejong commented Apr 4, 2023

josdejong commented Aug 23, 2023

Bobingstern commented Mar 25, 2023 •

edited

Loading