Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Riemann Zeta Function Implemented #2950

Open
wants to merge 12 commits into
base: develop
Choose a base branch
from
2 changes: 2 additions & 0 deletions src/expression/embeddedDocs/embeddedDocs.js
Original file line number Diff line number Diff line change
Expand Up @@ -182,6 +182,7 @@ import { setSizeDocs } from './function/set/setSize.js'
import { setSymDifferenceDocs } from './function/set/setSymDifference.js'
import { setUnionDocs } from './function/set/setUnion.js'
import { erfDocs } from './function/special/erf.js'
import { zetaDocs } from './function/special/riemannZeta.js'
import { madDocs } from './function/statistics/mad.js'
import { maxDocs } from './function/statistics/max.js'
import { meanDocs } from './function/statistics/mean.js'
Expand Down Expand Up @@ -517,6 +518,7 @@ export const embeddedDocs = {

// functions - special
erf: erfDocs,
zeta: zetaDocs,

// functions - statistics
cumsum: cumSumDocs,
Expand Down
14 changes: 14 additions & 0 deletions src/expression/embeddedDocs/function/special/riemannZeta.js
josdejong marked this conversation as resolved.
Show resolved Hide resolved
Original file line number Diff line number Diff line change
@@ -0,0 +1,14 @@
export const zetaDocs = {
name: 'zeta',
category: 'Special',
syntax: [
'zeta(s)'
],
description: 'Compute the Riemann Zeta Function using an infinite series and Riemanns Functional Equation for the entire complex plane',
examples: [
'zeta(0.2)',
'zeta(-0.5)',
'zeta(4)'
],
seealso: []
}
1 change: 1 addition & 0 deletions src/factoriesAny.js
Original file line number Diff line number Diff line change
Expand Up @@ -95,6 +95,7 @@ export { createZeros } from './function/matrix/zeros.js'
export { createFft } from './function/matrix/fft.js'
export { createIfft } from './function/matrix/ifft.js'
export { createErf } from './function/special/erf.js'
export { createRiemannZeta } from './function/special/riemannZeta.js'
export { createMode } from './function/statistics/mode.js'
export { createProd } from './function/statistics/prod.js'
export { createFormat } from './function/string/format.js'
Expand Down
2 changes: 1 addition & 1 deletion src/factoriesNumber.js
Original file line number Diff line number Diff line change
Expand Up @@ -249,7 +249,7 @@ export { createUnequalNumber as createUnequal } from './function/relational/uneq

// special
export { createErf } from './function/special/erf.js'

export { createRiemannZeta } from './function/special/riemannZeta.js'
// statistics
export { createMode } from './function/statistics/mode.js'
export { createProd } from './function/statistics/prod.js'
Expand Down
140 changes: 140 additions & 0 deletions src/function/special/riemannZeta.js
Original file line number Diff line number Diff line change
@@ -0,0 +1,140 @@
import { factory } from '../../utils/factory.js'

const name = 'zeta'
const dependencies = ['typed', 'config', 'multiply', 'pow', 'divide', 'factorial', 'equal', 'gamma', 'sin', 'subtract', 'add', 'Complex', 'BigNumber']

export const createRiemannZeta = /* #__PURE__ */ factory(name, dependencies, ({ typed, config, multiply, pow, divide, factorial, equal, gamma, sin, subtract, add, Complex, BigNumber }) => {
/**
* Compute the Riemann Zeta function of a value using an infinite series for
josdejong marked this conversation as resolved.
Show resolved Hide resolved
* all of the complex plane using Riemann's Functional equation.
*
* Based off the paper by Xavier Gourdon and Pascal Sebah
* ( http://numbers.computation.free.fr/Constants/Miscellaneous/zetaevaluations.pdf )
*
* Implementation and slight modification by Anik Patel
*
*
* Syntax:
*
* math.zeta(n)
*
* Examples:
*
* math.zeta(5) // returns 1.0369277551433699...
* math.zeta(-0.5) // returns -0.2078862249773...
* math.zeta(math.i) // returns 0.0033002..., -0.4181554491...i
*
*
* @param {number | Complex | BigNumber} s A Real, Complex or BigNumber
* @return {number | Complex | BigNumber} The Riemann Zeta of `s`
*/
function zeta (s) {
if (s.re === 0 && s.im === 0) {
return -0.5
}
if (s.re === 1) {
return NaN
}
if (s.re === Infinity && s.im === 0) {
return 1
}
if (s.im === Infinity || s.re === -Infinity) {
return NaN
}
// The "15" is number of decimal places of accuracy (approx)
const n = Math.round(1.3 * 15 + 0.9 * Math.abs(s.im))
if (s.re > -(n - 1) / 2) { return f(s, n) }

// Function Equation for reflection to re(s) < 1
let c = multiply(pow(2, s), pow(Math.PI, subtract(s, 1)))
c = multiply(c, (sin(multiply(Math.PI / 2, s))))
c = multiply(c, gamma(subtract(1, s)))
return multiply(c, zeta(subtract(1, s)))
}
// Big Number alias
function zetaBigNumber (x) {
if (x === 0) {
return -0.5
}
if (equal(x, 1)) {
return NaN
}
if (!x.isFinite() && !x.isNegative()) {
return 1
}
if (!x.isFinite() && x.isNegative()) {
return NaN
}
// 15 decimals of accuracy
const n = 20
if (x > -(n - 1) / 2) { return fBigNumber(x, n) }

// Function Equation for reflection to x < 1
let c = multiply(pow(2, x), pow(BigNumber(Math.PI), subtract(x, 1)))
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This will not result in a true BigNumber version of pi. I think you can simply import pi instead? Or use the function createBigNumberPi?

c = multiply(c, (sin(multiply(BigNumber(Math.PI / 2), x))))
c = multiply(c, gamma(subtract(1, x)))
return multiply(c, zetaBigNumber(subtract(1, x)))
}

/**
* Calculate a portion of the sum
* @param {number} k a positive integer
* @param {number} n a positive integer
* @return {number} the portion of the sum
**/
function d (k, n) {
let S = 0
for (let j = k; j <= n; j++) {
S += factorial(n + j - 1) * (4 ** j) / (factorial(n - j) * factorial(2 * j))
}

return n * S
}
// Big Number alias
function dBigNumber (k, n) {
josdejong marked this conversation as resolved.
Show resolved Hide resolved
let S = BigNumber(0)
const bn = BigNumber(n)
for (let j = k; j <= n; j++) {
const bj = BigNumber(j)
const factor = divide(multiply(factorial(add(bn, subtract(bj, 1))), pow(4, bj)), multiply(factorial(subtract(bn, bj)), factorial(multiply(2, bj))))
S = add(S, factor)
}

return multiply(n, S)
}
/**
* Calculate the positive Riemann Zeta function
* @param {number} s a real or complex number with s.re > 1
* @param {number} n a positive integer
* @return {number} Zeta of s
**/
function f (s, n) {
const c = divide(1, multiply(d(0, n), subtract(1, pow(2, subtract(1, s)))))
let S = new Complex(0, 0)
for (let k = 1; k <= n; k++) {
S = S.add(
divide((-1) ** (k - 1) * d(k, n), pow(k, s))
)
}

return multiply(c, S)
}
// Big Number alias
function fBigNumber (s, n) {
const c = divide(1, multiply(dBigNumber(0, n), subtract(1, pow(2, subtract(1, s)))))
let S = new BigNumber(0)
for (let k = 1; k <= n; k++) {
S = S.add(
divide(multiply((-1) ** (k - 1), dBigNumber(k, n)), pow(k, s))
)
}
return multiply(c, S)
}

// Return the value of the function
return typed(name, {
number: s => zeta(new Complex(s)),
Complex: zeta,
BigNumber: zetaBigNumber
})
})
69 changes: 69 additions & 0 deletions test/unit-tests/function/special/riemannZeta.test.js
Original file line number Diff line number Diff line change
@@ -0,0 +1,69 @@
/* eslint-disable no-loss-of-precision */

import assert from 'assert'
import approx from '../../../../tools/approx.js'
import math from '../../../../src/defaultInstance.js'
const zeta = math.zeta

describe('Riemann Zeta', function () {
it('should calculate the Riemann Zeta Function of a positive integer', function () {
assert.ok(isNaN(zeta(1)))
approx.equal(zeta(2), 1.6449340668482264)
approx.equal(zeta(3), 1.2020569031595942)
approx.equal(zeta(4), 1.0823232337111381)
approx.equal(zeta(5), 1.0369277551433699)
assert.strictEqual(zeta(Infinity), 1) // shouldn't stall
})
it('should calculate the Riemann Zeta Function of a nonpositive integer', function () {
assert.strictEqual(zeta(0), -0.5)
approx.equal(zeta(-1), -1 / 12)
approx.equal(zeta(-2), 0)
approx.equal(zeta(-3), 1 / 120)
approx.equal(zeta(-13), -1 / 12)
assert.ok(isNaN(zeta(-Infinity)))
})
it('should calculate the Riemann Zeta Function of a Big Number', function () {
assert.ok(isNaN(zeta(math.BigNumber(1))))
approx.equal(zeta(math.BigNumber(2)), 1.6449340668482264)
approx.equal(zeta(math.BigNumber(-2)), 0)
approx.equal(zeta(math.BigNumber(20)), 1.000000953962033872)
approx.equal(zeta(math.BigNumber(-21)), -281.4601449275362318)
approx.equal(zeta(math.BigNumber(50)), 1.0000000000000008881)
approx.equal(zeta(math.BigNumber(-211)), 2.727488e231)
approx.equal(zeta(math.BigNumber(100)), 1.000000000000000000)
assert.strictEqual(zeta(Infinity), 1) // shouldn't stall
})

it('should calculate the Riemann Zeta Function of a rational number', function () {
approx.equal(zeta(0.125), -0.6327756234986952552935)
approx.equal(zeta(0.25), -0.81327840526189165652144)
approx.equal(zeta(0.5), -1.460354508809586812889499)
approx.equal(zeta(1.5), 2.61237534868548834334856756)
approx.equal(zeta(2.5), 1.34148725725091717975676969)
approx.equal(zeta(3.5), 1.12673386731705664642781249)
approx.equal(zeta(30.5), 1.00000000065854731257004465)
approx.equal(zeta(144.9), 1.0000000000000000000000000)

approx.equal(zeta(-0.5), -0.2078862249773545660173067)
approx.equal(zeta(-1.5), -0.0254852018898330359495429)
approx.equal(zeta(-2.5), 0.00851692877785033054235856)
})

it('should calculate the Riemann Zeta Function of an irrational number', function () {
approx.equal(zeta(Math.SQRT2), 3.0207376794860326682709)
approx.equal(zeta(Math.PI), 1.17624173838258275887215)
approx.equal(zeta(Math.E), 1.26900960433571711576556)

approx.equal(zeta(-Math.SQRT2), -0.0325059805396893552173896)
approx.equal(zeta(-Math.PI), 0.00744304047846672771406904635)
approx.equal(zeta(-Math.E), 0.00915987755942023170457566822)
})
it('should calculate the Riemann Zeta Function of a Complex number', function () {
approx.equal(zeta(math.complex(0, 1)), math.complex(0.00330022368532410287421711, -0.418155449141321676689274239))
approx.equal(zeta(math.complex(3, 2)), math.complex(0.97304196041894244856408189, -0.1476955930004537946298999860))
approx.equal(zeta(math.complex(-1.5, 3.7)), math.complex(0.244513626137832304395, 0.2077842378226353306923615))
approx.equal(zeta(math.complex(3.9, -5.2)), math.complex(0.952389583517691366229, -0.03276345793831000384775143962))
approx.equal(zeta(math.complex(-1.2, -9.3)), math.complex(2.209608454242663005234, -0.67864225792147162441259999407))
approx.equal(zeta(math.complex(0.5, 14.14)), math.complex(-0.00064921838659084911, 0.004134963322496717323762898714))
})
})
18 changes: 18 additions & 0 deletions types/index.d.ts
Original file line number Diff line number Diff line change
Expand Up @@ -2574,6 +2574,16 @@ declare namespace math {
*/
erf<T extends number | MathCollection>(x: T): NoLiteralType<T>

/**
* Compute the Riemann Zeta function of a value using an infinite series
* and Riemann's Functional equation.
* @param s A real, complex or BigNumber
* @returns The Riemann Zeta of s
*/
riemannZeta<T extends number | Complex | BigNumber | MathCollection>(
x: T
): NoLiteralType<T>

/*************************************************************************
* Statistics functions
************************************************************************/
Expand Down Expand Up @@ -5831,6 +5841,14 @@ declare namespace math {
this: MathJsChain<T>
): MathJsChain<NoLiteralType<T>>

/**
* Compute the Riemann Zeta function of a value using an infinite series
* and Riemann's Functional equation.
*/
riemannZeta<T extends number | Complex | BigNumber | MathCollection>(
this: MathJsChain<T>
): MathJsChain<NoLiteralType<T>>

/*************************************************************************
* Statistics functions
************************************************************************/
Expand Down