Riemann Zeta Function Implemented

src/expression/embeddedDocs/embeddedDocs.js

@@ -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'
@@ -517,6 +518,7 @@ export const embeddedDocs = {
 
   // functions - special
   erf: erfDocs,
+  zeta: zetaDocs,
 
   // functions - statistics
   cumsum: cumSumDocs,

src/expression/embeddedDocs/function/special/riemannZeta.js

@@ -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: []
+}

src/factoriesAny.js

@@ -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'

src/factoriesNumber.js

@@ -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'

src/function/special/riemannZeta.js

@@ -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
+   * 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 or complex number
+   * @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)))
+    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) {
+    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
+  })
+})

test/unit-tests/function/special/riemannZeta.test.js

@@ -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))
+  })
+})

types/index.d.ts

@@ -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 x A real or complex number
+     * @returns The Riemann Zeta of x
+     */
+    riemannZeta<T extends number | Complex | MathCollection>(
+      x: T
+    ): NoLiteralType<T>
+
     /*************************************************************************
      * Statistics functions
      ************************************************************************/
@@ -5833,6 +5843,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 | MathCollection>(
+      this: MathJsChain<T>
+    ): MathJsChain<NoLiteralType<T>>
+
     /*************************************************************************
      * Statistics functions
      ************************************************************************/