Additional numbertheoretic functions#

Mpmath: Prime counting function#

ctx.primepi(x)#

where ctx is dec, mpm, fpm, gmp or apm.

Returns the prime counting function.

See also Wikipedia [1395], MathWorld [1004], Mpmath [761].

Evaluates the prime counting function, \(\pi(x)\), which gives the number of primes less than or equal to \(x\). The argument \(x\) may be fractional.

The prime counting function is very expensive to evaluate precisely for large \(x\), and the present implementation is not optimized in any way. For numerical approximation of the prime counting function, it is better to use primepi2() or riemannr().

Some values of the prime counting method:

>>> from mpfunlab import *
>>> [primepi(k) for k in range(20)]
[0, 0, 1, 2, 2, 3, 3, 4, 4, 4, 4, 5, 5, 6, 6, 6, 6, 7, 7, 8]
>>> primepi(3.5)
2
>>> primepi(100000)
9592

Mpmath: Mangoldt function#

ctx.mangoldt(n)#

where ctx is dec, mpm, fpm, gmp or apm.

Returns the Mangoldt function. See also Wikipedia [1479], MathWorld [1101], Mpmath [760].

Evaluates the von Mangoldt function \(\Lambda(n) = \log p\) if \(n = p^k\) a power of a prime, and \(\Lambda(n) = 0\) otherwise:

>>> from mpfunlab import *
>>> mp.dps = 25; mp.pretty = True
>>> [mangoldt(n) for n in range(-2,3)]
[0.0, 0.0, 0.0, 0.0, 0.6931471805599453094172321]
>>> mangoldt(6)
0.0
>>> mangoldt(7)
1.945910149055313305105353
>>> mangoldt(8)
0.6931471805599453094172321
>>> fsum(mangoldt(n) for n in range(101))
94.04531122935739224600493
>>> fsum(mangoldt(n) for n in range(10001))
10013.39669326311478372032

Mpmath: Upper bound for the value of the prime counting function#

ctx.primepi2_upper(x)#

where ctx is dec, mpm, fpm, gmp or apm.

Returns bounds for the value of the prime counting function. See also Wikipedia [1394], Wikipedia [1395], MathWorld [1004], MathWorld [1005].

Returns an interval (as an mpi instance) providing bounds for the value of the prime counting function \(\pi(x)\). For small \(x\), primepi2() returns an exact interval based on the output of primepi(). For \(x > 2656\), a loose interval based on Schoenfeld’s inequality

\[|\pi(x) - \mathrm{li}(x)| < \frac{\sqrt x \log x}{8 \pi}\]

is returned. This estimate is rigorous assuming the truth of the Riemann hypothesis, and can be computed very quickly.

An example:

>>> from mpfunlab import dec, mpm, gmp, fpm, apm
>>> mpm.dps = 40; x = '1E+10'
>>> dx = dec.primepi2_upper(x); mx = mpm.primepi2_upper(x); gx = gmp.primepi2_upper(x)
>>> fx = fpm.primepi2_upper(x); ax = apm.primepi2_upper(x)
>>> mpm.show([dx, mx, gx, fx, ax])
dec:  4.551472320000000000000000000000000000000E+8
mpm:  4.551472320000000000000000000000000000000e+8
gmp:  4.551472320000000000000000000000000000000E+08
fpm:  4.55147232000000E+08
apm:  4.551472320000000000000000000000000000000e+8 (0.0%)

Mpmath: Lower dound for the value of the prime counting function#

ctx.primepi2_lower(x)#

where ctx is dec, mpm, fpm, gmp or apm.

Returns bounds for the value of the prime counting function. See also Wikipedia [1394], Wikipedia [1395], MathWorld [1004], MathWorld [1005].

\[|\pi(x) - \mathrm{li}(x)| < \frac{\sqrt x \log x}{8 \pi}\]

is returned. This estimate is rigorous assuming the truth of the Riemann hypothesis, and can be computed very quickly.

An example:

>>> from mpfunlab import dec, mpm, gmp, fpm, apm
>>> mpm.dps = 40; x = '1E+10'
>>> dx = dec.primepi2_lower(x); mx = mpm.primepi2_lower(x); gx = gmp.primepi2_lower(x)
>>> fx = fpm.primepi2_lower(x); ax = apm.primepi2_lower(x)
>>> mpm.show([dx, mx, gx, fx, ax])
dec:  4.549639970000000000000000000000000000000E+8
mpm:  4.549639970000000000000000000000000000000e+8
gmp:  4.549639970000000000000000000000000000000E+08
fpm:  4.54963997000000E+08
apm:  4.549639970000000000000000000000000000000e+8 (0.0%)

Mpmath, DAMath: Riemann R function#

ctx.riemann_r(x)#

where ctx is dec, mpm, fpm, gmp or apm.

Returns the Riemann R function. See also MathWorld [1005], Mpmath [648]., Mpmath [649].

Evaluates the Riemann R function, a smooth approximation of the prime counting function. The Riemann R function gives a fast numerical approximation useful e.g. to roughly estimate the number of primes in a given interval.

The Riemann R function is computed using the rapidly convergent Gram series,

\[R(x) = 1 + \sum_{k=1}^{\infty} \frac{\log^k x}{k k! \zeta(k+1)}.\]

From the Gram series, one sees that the Riemann R function is a well-defined analytic function (except for a branch cut along the negative real half-axis); it can be evaluated for arbitrary real or complex arguments.

An example with real input:

>>> from xlcalcnet import dec, mpm, gmp, fpm, apm
>>> mpm.dps = 40; x = '7.5'
>>> \mathrm{d}x = dec.riemannr(x); mx = mpm.riemannr(x); gx = gmp.riemannr(x)
>>> fx = fpm.riemannr(x); ax = apm.riemannr(x)
>>> mpm.show([\mathrm{d}x, mx, gx, fx, ax])
dec:  3.729347432649662619188571351358357974070E+0
mpm:  3.729347432649662619188571351358357974070e+0
gmp:  3.729347432649662619188571351358357974070E+00
fpm:  3.72934743264966E+00
mpm:  3.729347432649662619188571351358357974070e+0

An example with complex input:

>>> from xlcalcnet import dec, mpm, gmp, fpm, apm
>>> mpm.dps = 20; z = '-4+2j'
>>> \mathrm{d}z = dec.riemannr(z); mz = mpm.riemannr(z); gz = gmp.riemannr(z)
>>> fz = fpm.riemannr(z); az = apm.riemannr(z)
>>> mpm.show([\mathrm{d}z, mz, gz, fz, az],  aligned=True)
dec: -5.5100220815548642759E-1  + 2.1696639813811945004E+0j
mpm: -5.5100220815548642759e-1  + 2.1696639813811945004e+0j
gmp: -5.5100220815548642759E-01 + 2.1696639813811945004E+00j
fpm: -5.51002208155486E-01      + 2.16966398138119E+00j
mpm: -5.5100220815548642759e-1  + 2.1696639813811945004e+0j

Mpmath, DAMath: Prime zeta function#

ctx.primezeta(s)#

where ctx is math53, dec, mpm, fpm, gmp or apm.

Returns for \(x > 0.2\) the prime zeta function \(\displaystyle P(x) = \sum_{p \text{ prime}} p^{-x}\), or its real part for \(x<1\).

Mpmath: Mertens constant#

ctx.const_mertens()#

where ctx is dec, mpm, fpm, gmp or apm.

Returns the Mertens constant. See also Wikipedia [1470], MathWorld [1102], Mpmath [739].

Represents the Mertens or Meissel-Mertens constant, which is the prime number analog of Euler’s constant:

\[B_1 = \lim_{N\to\infty} \left(\sum_{p_k \le N} \frac{1}{p_k} - \log \log N \right)\]

Here \(p_k\) denotes the \(k\)-th prime number. Other names for this constant include the Hadamard-de la Vallee-Poussin constant or the prime reciprocal constant.

Beispielsweise:

\[M=\gamma +\sum _{k=2}^{\infty }{\frac {\mu (k)}{k}}\log {\bigg (}\zeta (k){\bigg )}\]

where \(\mu (n)\) denotes the Möbius function and \(\zeta (n)\) the Riemann zeta function.

The following gives the Mertens constant to 50 digits:

>>> from mpfunlab import *
>>> mp.dps = 50; mp.pretty = True
>>> +mertens
0.2614972128476427837554268386086958590515666482612

Mpmath: Twin prime constant#

ctx.const_twinprime()#

where ctx is dec, mpm, fpm, gmp or apm.

Returns the Twin prime constant. See also Wikipedia [1478], MathWorld [1108], Mpmath [748].

Represents the twin prime constant, which is the factor \(C_2\) featuring in the Hardy-Littlewood conjecture for the growth of the twin prime counting function,

\[\pi_2(n) \sim 2 C_2 \frac{n}{\log^2 n}.\]

It is given by the product over primes

\[C_2 = \prod_{p\ge3} \frac{p(p-2)}{(p-1)^2} \approx 0.66016\]

See also code in mpmath: mpmath/libmp/gammazeta.py, def twinprime_fixed(prec)

See also code in mpmath: mpmath/libmp/libintmath.py, def moebius(n)

Flajolet and Vardi (1996) give series with accelerated convergence (see mathworld)

\[C_2 = \prod_{n=2} ^ {\infty} \left[ \zeta(n)(1-2^{-n}) \right]^{-I_n}, \quad \text{with } I_n = \frac{1}{n} \sum_{d|n} \mu(d) 2^{n/d}.\]

Computing \(C_2\) to 50 digits:

>>> from mpfunlab import *
>>> mp.dps = 50; mp.pretty = True
>>> +twinprime
0.66016181584686957392781211001455577843262336028473

Mpmath: Cyclotomic polynomial#

ctx.cyclotomic(n, x)#

where ctx is fpm, mpm, ipm, dec, gmp or apm.

Returns the cyclotomic polynomial. See also Wikipedia [1466], MathWorld [1096], Mpmath [750].

Evaluates the cyclotomic polynomial \(\Phi_n(x)\), defined by

\[\Phi_n(x) = \prod_{\zeta} (x - \zeta)\]

where \(\zeta\) ranges over all primitive \(n\)-th roots of unity (see unitroots()). An equivalent representation, used for computation, is

\[\Phi_n(x) = \prod_{d\mid n}(x^d-1)^{\mu(n/d)} = \Phi_n(x)\]

where \(\mu(m)\) denotes the Moebius function.

The definition as a product over primitive roots may be checked by computing the product explicitly (for a real argument, this method will generally introduce numerical noise in the imaginary part):

>>> mp.dps = 25
>>> z = 3+4j
>>> cyclotomic(10, z)
(-419.0 - 360.0j)
>>> fprod(z-r for r in unitroots(10, primitive=True))
(-419.0 - 360.0j)
>>> z = 3
>>> cyclotomic(10, z)
61.0
>>> fprod(z-r for r in unitroots(10, primitive=True))
(61.0 - 3.146045605088568607055454e-25j)

Mpmath: Stirling number of the first kind#

ctx.stirling1(n, k, exact=False)#

where ctx is fpm, mpm, ipm, dec, gmp or apm.

Returns the Stirling number of the first kind. See also Wikipedia [1476], MathWorld [1105], NIST [43], Mpmath [746].

Gives the Stirling number of the first kind \(s(n,k)\), defined by

\[x(x-1)(x-2)\cdots(x-n+1) = \sum_{k=0}^n s(n,k) x^k.\]

The value is computed using an integer recurrence. The implementation is not optimized for approximating large values quickly.

Examples

Comparing with the generating method:

>>> from xlcalcnet import *
>>> mp.dps = 25; mp.pretty = True
>>> taylor(lambda x: ff(x, 5), 0, 5)
[0.0, 24.0, -50.0, 35.0, -10.0, 1.0]
>>> [stirling1(5, k) for k in range(6)]
[0.0, 24.0, -50.0, 35.0, -10.0, 1.0]

Recurrence relation:

>>> n, k = 5, 3
>>> stirling1(n+1,k) + n*stirling1(n,k) - stirling1(n,k-1)
0.0

Pass exact=True to obtain exact values of Stirling numbers as integers:

>>> stirling1(42, 5)
-2.864498971768501633736628e+50
>>> print(stirling1(42, 5, exact=True))
-286449897176850163373662803014001546235808317440000

Mpmath: Stirling number of the second kind#

ctx.stirling2(n, k, exact=False)#

where ctx is fpm, mpm, ipm, dec, gmp or apm.

Returns the Stirling number of the second kind. See also Wikipedia [1477], MathWorld [1106], NIST [43], Mpmath [747].

Gives the Stirling number of the second kind \(S(n,k)\), defined by

\[x^n = \sum_{k=0}^n S(n,k) x(x-1)(x-2)\cdots(x-k+1)\]

The value is computed using integer arithmetic to evaluate a power sum. The implementation is not optimized for approximating large values quickly.

Examples

Comparing with the generating method:

>>> from xlcalcnet import *
>>> mp.dps = 25; mp.pretty = True
>>> taylor(lambda x: sum(stirling2(5,k) * ff(x,k) for k in range(6)), 0, 5)
[0.0, 0.0, 0.0, 0.0, 0.0, 1.0]

Recurrence relation:

>>> n, k = 5, 3
>>> stirling2(n+1,k) - k*stirling2(n,k) - stirling2(n,k-1)
0.0

Pass exact=True to obtain exact values of Stirling numbers as integers:

>>> stirling2(52, 10)
2.641822121003543906807485e+45
>>> print(stirling2(52, 10, exact=True))
2641822121003543906807485307053638921722527655

Mpmath: CODE!! Polyexponential function#

ctx.polyexp(s, z)#

where ctx is fpm, mpm, ipm, dec, gmp or apm.

Returns the polyexponential function. See Mpmath [741].

Evaluates the polyexponential function, defined for arbitrary complex \(s\), \(z\) by the series

\[E_s(z) = \sum_{k=1}^{\infty} \frac{k^s}{k!} z^k.\]

\(E_s(z)\) is constructed from the exponential function analogously to how the polylogarithm is constructed from the ordinary logarithm; as a function of \(s\) (with \(z\) fixed), \(E_s\) is an L-series It is an entire function of both \(s\) and \(z\).

The polyexponential function provides a generalization of the Touchard polynomials to noninteger orders \(n\). In terms of the Bell polynomials,

\[E_s(z) = e^z B_s(z) - \mathrm{sinc}(\pi s).\]

Note that \(B_n(x)\) and \(e^{-x} E_n(x)\) are identical if \(n\) is a nonzero integer, but not otherwise. In particular, they differ at \(n = 0\).

Examples

Evaluating a series:

>>> from mpfunlab import *
>>> mp.dps = 25; mp.pretty = True
>>> nsum(lambda k: sqrt(k)/fac(k), [1,inf])
2.101755547733791780315904
>>> polyexp(0.5,1)
2.101755547733791780315904

Evaluation for arbitrary arguments:

>>> polyexp(-3-4j, 2.5+2j)
(2.351660261190434618268706 + 1.202966666673054671364215j)

Evaluation is accurate for tiny function values:

>>> polyexp(4, -100)
3.499471750566824369520223e-36

If \(n\) is a nonpositive integer, \(E_n\) reduces to a special instance of the hypergeometric function \(\,_pF_q\):

>>> n = 3
>>> x = pi
>>> polyexp(-n,x)
4.042192318847986561771779
>>> x*hyper([1]*(n+1), [2]*(n+1), x)
4.042192318847986561771779

Mpmath, DAMath: Moebius function, \(\mu(n)\)#

math53.moebius(n)#

Returns \(\mu(n)\), the Möbius function.

For integer \(n\), the Moebius \(\mu\) function \(\mu(n)\) equals 0 if \(n\) has repeated integer factors. Otherwise, if \(n\) is the product of \(k\) distinct primes, the Moebius \(\mu\) function \(\mu(n)\) equals \((-1)^k\).

https://en.wikipedia.org/wiki/M%C3%B6bius_function

An example in Python

>>> from xlcalcnet import xreal
>>> xreal.Moebius(3)
xreal('5.2359877559829887307E-1')
>>> xreal.Moebius)
xreal('5.3518479027559984754E-1')

An example in Visual Basic

>>> from xlcalcnet import Gpr
>>> Gpr.Moebius(3)
Gpr('5.2359877559829887307E-1')
>>> Gpr.Moebius)
Gpr('5.3518479027559984754E-1')

Additional numbertheoretic functions

Contents

Additional numbertheoretic functions#

Mpmath: Prime counting function#

Mpmath: Mangoldt function#

Mpmath: Upper bound for the value of the prime counting function#

Mpmath: Lower dound for the value of the prime counting function#

Mpmath, DAMath: Riemann R function#

Mpmath, DAMath: Prime zeta function#

Mpmath: Mertens constant#

Mpmath: Twin prime constant#

Mpmath: Cyclotomic polynomial#

Mpmath: Stirling number of the first kind#

Mpmath: Stirling number of the second kind#

Mpmath: CODE!! Polyexponential function#

Mpmath, DAMath: Moebius function, \(\mu(n)\)#