Riemann Zeta Function Implemented

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,87 @@
+import { factory } from '../../utils/factory.js'
+
+const name = 'riemannZeta'
+const dependencies = ['typed', 'config', 'multiply', 'pow', 'divide', 'factorial', 'gamma', 'sin', 'subtract', 'add', 'Complex']
+
+export const createRiemannZeta = /* #__PURE__ */ factory(name, dependencies, ({ typed, config, multiply, pow, divide, factorial, gamma, sin, subtract, add, Complex }) => {
+  /**
+   * Compute the Riemann Zeta function of a value using an infinite series for
+   * all of the complex plane using Riemann's Functional equation.
+   *
+   *
+   * Syntax:
+   *
+   *    math.riemannZeta(n)
+   *
+   * Examples:
+   *  FIX::
+   *    math.riemannZeta(5)       // returns 1.0369277551433699...
+   *    math.riemannZeta(-0.5)    // returns -0.2078862249773...
+   *    math.riemannZeta(math.i)  // returns 0.0033002..., -0.4181554491...i
+   *
+   *
+   * @param {number | Complex} z   A real or complex number
+   * @return {number | Complex}    The zeta of `z`
+   */
+  function zeta (s) {
+    s = Complex(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)))
+  }
+
+  /**
+   * 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
+  }
+  /**
+   * 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)
+  }
+
+  return typed(name, {
+    number: zeta,
+    Complex: zeta
+  })
+})

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

@@ -0,0 +1,59 @@
+/* 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.riemannZeta
+
+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 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
      ************************************************************************/