-
-
Notifications
You must be signed in to change notification settings - Fork 1.2k
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
Bobingstern
wants to merge
12
commits into
josdejong:develop
Choose a base branch
from
Bobingstern:develop
base: develop
Could not load branches
Branch not found: {{ refName }}
Could not load tags
Nothing to show
Are you sure you want to change the base?
Some commits from the old base branch may be removed from the timeline,
and old review comments may become outdated.
Open
Changes from 10 commits
Commits
Show all changes
12 commits
Select commit
Hold shift + click to select a range
23b4819
Riemann Zeta Function
Bobingstern b67e2c3
Big Number zeta and added docs
Bobingstern 0c4aea5
Merge branch 'develop' into develop
Bobingstern c507716
Original algorithm paper credited
Bobingstern 1db74c1
Merge branch 'develop' of https://github.com/Bobingstern/mathjs into …
Bobingstern 3ddb7dc
Update index.d.ts
Bobingstern 19cbf65
Update riemannZeta.js
Bobingstern dc7db48
Update index.d.ts
Bobingstern f1fa319
Merge branch 'develop' into develop
Bobingstern 59290c5
Renamed files to reflect zeta
Bobingstern 19aae26
Merge branch 'develop' into develop
josdejong 2e67f57
Merge branch 'develop' into develop
josdejong File filter
Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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: [] | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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', 'pi'] | ||
|
||
export const createZeta = /* #__PURE__ */ factory(name, dependencies, ({ typed, config, multiply, pow, divide, factorial, equal, gamma, sin, subtract, add, Complex, BigNumber, pi }) => { | ||
/** | ||
* Compute the Riemann Zeta function of a value using an infinite series for | ||
* 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 parameter to to the Riemann Zeta Function | ||
* @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) } | ||
|
||
// Functional Equation for reflection to re(s) < -( n - 1) / 2 | ||
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) } | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. There are still a few places where you implicitly convert a BigNumber into anumber, like |
||
|
||
// Function Equation for reflection to x < 1 | ||
let c = multiply(pow(2, x), pow(new BigNumber(pi), subtract(x, 1))) | ||
c = multiply(c, (sin(multiply(new BigNumber(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) { | ||
let S = new BigNumber(0) | ||
const bn = new BigNumber(n) | ||
for (let j = k; j <= n; j++) { | ||
const bj = new 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} Riemann 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 | ||
}) | ||
}) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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)) | ||
}) | ||
}) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
by default,
BigNumbers
are configured to have a precision of 64 digits. Can we utilizeconfig.precision
here?