Number-theoretical, combinatorial and integer functions¶
For factorial-type functions, including binomial coefficients, double factorials, etc., see the separate section Factorials and gamma functions.
Fibonacci numbers¶
fibonacci()/fib()¶
-
mpmath.fibonacci(n, **kwargs)¶ fibonacci(n)computes the \(n\)-th Fibonacci number, \(F(n)\). The Fibonacci numbers are defined by the recurrence \(F(n) = F(n-1) + F(n-2)\) with the initial values \(F(0) = 0\), \(F(1) = 1\).fibonacci()extends this definition to arbitrary real and complex arguments using the formula\[F(z) = \frac{\phi^z - \cos(\pi z) \phi^{-z}}{\sqrt 5}\]where \(\phi\) is the golden ratio.
fibonacci()also uses this continuous formula to compute \(F(n)\) for extremely large \(n\), where calculating the exact integer would be wasteful.For convenience,
fib()is available as an alias forfibonacci().Basic examples
Some small Fibonacci numbers are:
>>> from mpmath import * >>> mp.dps = 15; mp.pretty = True >>> for i in range(10): ... print(fibonacci(i)) ... 0.0 1.0 1.0 2.0 3.0 5.0 8.0 13.0 21.0 34.0 >>> fibonacci(50) 12586269025.0
The recurrence for \(F(n)\) extends backwards to negative \(n\):
>>> for i in range(10): ... print(fibonacci(-i)) ... 0.0 1.0 -1.0 2.0 -3.0 5.0 -8.0 13.0 -21.0 34.0
Large Fibonacci numbers will be computed approximately unless the precision is set high enough:
>>> fib(200) 2.8057117299251e+41 >>> mp.dps = 45 >>> fib(200) 280571172992510140037611932413038677189525.0
fibonacci()can compute approximate Fibonacci numbers of stupendous size:>>> mp.dps = 15 >>> fibonacci(10**25) 3.49052338550226e+2089876402499787337692720
Real and complex arguments
The extended Fibonacci function is an analytic function. The property \(F(z) = F(z-1) + F(z-2)\) holds for arbitrary \(z\):
>>> mp.dps = 15 >>> fib(pi) 2.1170270579161 >>> fib(pi-1) + fib(pi-2) 2.1170270579161 >>> fib(3+4j) (-5248.51130728372 - 14195.962288353j) >>> fib(2+4j) + fib(1+4j) (-5248.51130728372 - 14195.962288353j)
The Fibonacci function has infinitely many roots on the negative half-real axis. The first root is at 0, the second is close to -0.18, and then there are infinitely many roots that asymptotically approach \(-n+1/2\):
>>> findroot(fib, -0.2) -0.183802359692956 >>> findroot(fib, -2) -1.57077646820395 >>> findroot(fib, -17) -16.4999999596115 >>> findroot(fib, -24) -23.5000000000479
Mathematical relationships
For large \(n\), \(F(n+1)/F(n)\) approaches the golden ratio:
>>> mp.dps = 50 >>> fibonacci(101)/fibonacci(100) 1.6180339887498948482045868343656381177203127439638 >>> +phi 1.6180339887498948482045868343656381177203091798058
The sum of reciprocal Fibonacci numbers converges to an irrational number for which no closed form expression is known:
>>> mp.dps = 15 >>> nsum(lambda n: 1/fib(n), [1, inf]) 3.35988566624318
Amazingly, however, the sum of odd-index reciprocal Fibonacci numbers can be expressed in terms of a Jacobi theta function:
>>> nsum(lambda n: 1/fib(2*n+1), [0, inf]) 1.82451515740692 >>> sqrt(5)*jtheta(2,0,(3-sqrt(5))/2)**2/4 1.82451515740692
Some related sums can be done in closed form:
>>> nsum(lambda k: 1/(1+fib(2*k+1)), [0, inf]) 1.11803398874989 >>> phi - 0.5 1.11803398874989 >>> f = lambda k:(-1)**(k+1) / sum(fib(n)**2 for n in range(1,int(k+1))) >>> nsum(f, [1, inf]) 0.618033988749895 >>> phi-1 0.618033988749895
References
Bernoulli numbers and polynomials¶
bernfrac()¶
-
mpmath.bernfrac(n)¶ Returns a tuple of integers \((p, q)\) such that \(p/q = B_n\) exactly, where \(B_n\) denotes the \(n\)-th Bernoulli number. The fraction is always reduced to lowest terms. Note that for \(n > 1\) and \(n\) odd, \(B_n = 0\), and \((0, 1)\) is returned.
Examples
The first few Bernoulli numbers are exactly:
>>> from mpmath import * >>> for n in range(15): ... p, q = bernfrac(n) ... print("%s %s/%s" % (n, p, q)) ... 0 1/1 1 -1/2 2 1/6 3 0/1 4 -1/30 5 0/1 6 1/42 7 0/1 8 -1/30 9 0/1 10 5/66 11 0/1 12 -691/2730 13 0/1 14 7/6
This function works for arbitrarily large \(n\):
>>> p, q = bernfrac(10**4) >>> print(q) 2338224387510 >>> print(len(str(p))) 27692 >>> mp.dps = 15 >>> print(mpf(p) / q) -9.04942396360948e+27677 >>> print(bernoulli(10**4)) -9.04942396360948e+27677
Note
bernoulli()computes a floating-point approximation directly, without computing the exact fraction first. This is much faster for large \(n\).Algorithm
bernfrac()works by computing the value of \(B_n\) numerically and then using the von Staudt-Clausen theorem [1] to reconstruct the exact fraction. For large \(n\), this is significantly faster than computing \(B_1, B_2, \ldots, B_2\) recursively with exact arithmetic. The implementation has been tested for \(n = 10^m\) up to \(m = 6\).In practice,
bernfrac()appears to be about three times slower than the specialized program calcbn.exe [2]References
MathWorld, von Staudt-Clausen Theorem: http://mathworld.wolfram.com/vonStaudt-ClausenTheorem.html
The Bernoulli Number Page: http://www.bernoulli.org/
bernpoly()¶
-
mpmath.bernpoly(n, z)¶ Evaluates the Bernoulli polynomial \(B_n(z)\).
The first few Bernoulli polynomials are:
>>> from mpmath import * >>> mp.dps = 15; mp.pretty = True >>> for n in range(6): ... nprint(chop(taylor(lambda x: bernpoly(n,x), 0, n))) ... [1.0] [-0.5, 1.0] [0.166667, -1.0, 1.0] [0.0, 0.5, -1.5, 1.0] [-0.0333333, 0.0, 1.0, -2.0, 1.0] [0.0, -0.166667, 0.0, 1.66667, -2.5, 1.0]
At \(z = 0\), the Bernoulli polynomial evaluates to a Bernoulli number (see
bernoulli()):>>> bernpoly(12, 0), bernoulli(12) (-0.253113553113553, -0.253113553113553) >>> bernpoly(13, 0), bernoulli(13) (0.0, 0.0)
Evaluation is accurate for large \(n\) and small \(z\):
>>> mp.dps = 25 >>> bernpoly(100, 0.5) 2.838224957069370695926416e+78 >>> bernpoly(1000, 10.5) 5.318704469415522036482914e+1769
Euler numbers and polynomials¶
eulernum()¶
-
mpmath.eulernum(n)¶ Gives the \(n\)-th Euler number, defined as the \(n\)-th derivative of \(\mathrm{sech}(t) = 1/\cosh(t)\) evaluated at \(t = 0\). Equivalently, the Euler numbers give the coefficients of the Taylor series
\[\mathrm{sech}(t) = \sum_{n=0}^{\infty} \frac{E_n}{n!} t^n.\]The Euler numbers are closely related to Bernoulli numbers and Bernoulli polynomials. They can also be evaluated in terms of Euler polynomials (see
eulerpoly()) as \(E_n = 2^n E_n(1/2)\).Examples
Computing the first few Euler numbers and verifying that they agree with the Taylor series:
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> [eulernum(n) for n in range(11)] [1.0, 0.0, -1.0, 0.0, 5.0, 0.0, -61.0, 0.0, 1385.0, 0.0, -50521.0] >>> chop(diffs(sech, 0, 10)) [1.0, 0.0, -1.0, 0.0, 5.0, 0.0, -61.0, 0.0, 1385.0, 0.0, -50521.0]
Euler numbers grow very rapidly.
eulernum()efficiently computes numerical approximations for large indices:>>> eulernum(50) -6.053285248188621896314384e+54 >>> eulernum(1000) 3.887561841253070615257336e+2371 >>> eulernum(10**20) 4.346791453661149089338186e+1936958564106659551331
Comparing with an asymptotic formula for the Euler numbers:
>>> n = 10**5 >>> (-1)**(n//2) * 8 * sqrt(n/(2*pi)) * (2*n/(pi*e))**n 3.69919063017432362805663e+436961 >>> eulernum(n) 3.699193712834466537941283e+436961
Pass
exact=Trueto obtain exact values of Euler numbers as integers:>>> print(eulernum(50, exact=True)) -6053285248188621896314383785111649088103498225146815121 >>> print(eulernum(200, exact=True) % 10**10) 1925859625 >>> eulernum(1001, exact=True) 0
eulerpoly()¶
-
mpmath.eulerpoly(n, z)¶ Evaluates the Euler polynomial \(E_n(z)\), defined by the generating function representation
\[\frac{2e^{zt}}{e^t+1} = \sum_{n=0}^\infty E_n(z) \frac{t^n}{n!}.\]The Euler polynomials may also be represented in terms of Bernoulli polynomials (see
bernpoly()) using various formulas, for example\[E_n(z) = \frac{2}{n+1} \left( B_n(z)-2^{n+1}B_n\left(\frac{z}{2}\right) \right).\]Special values include the Euler numbers \(E_n = 2^n E_n(1/2)\) (see
eulernum()).Examples
Computing the coefficients of the first few Euler polynomials:
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> for n in range(6): ... chop(taylor(lambda z: eulerpoly(n,z), 0, n)) ... [1.0] [-0.5, 1.0] [0.0, -1.0, 1.0] [0.25, 0.0, -1.5, 1.0] [0.0, 1.0, 0.0, -2.0, 1.0] [-0.5, 0.0, 2.5, 0.0, -2.5, 1.0]
Evaluation for arbitrary \(z\):
>>> eulerpoly(2,3) 6.0 >>> eulerpoly(5,4) 423.5 >>> eulerpoly(35, 11111111112) 3.994957561486776072734601e+351 >>> eulerpoly(4, 10+20j) (-47990.0 - 235980.0j) >>> eulerpoly(2, '-3.5e-5') 0.000035001225 >>> eulerpoly(3, 0.5) 0.0 >>> eulerpoly(55, -10**80) -1.0e+4400 >>> eulerpoly(5, -inf) -inf >>> eulerpoly(6, -inf) +inf
Computing Euler numbers:
>>> 2**26 * eulerpoly(26,0.5) -4087072509293123892361.0 >>> eulernum(26) -4087072509293123892361.0
Evaluation is accurate for large \(n\) and small \(z\):
>>> eulerpoly(100, 0.5) 2.29047999988194114177943e+108 >>> eulerpoly(1000, 10.5) 3.628120031122876847764566e+2070 >>> eulerpoly(10000, 10.5) 1.149364285543783412210773e+30688
Bell numbers and polynomials¶
bell()¶
-
mpmath.bell(n, x)¶ For \(n\) a nonnegative integer,
bell(n,x)evaluates the Bell polynomial \(B_n(x)\), the first few of which are\[ \begin{align}\begin{aligned}B_0(x) = 1\\B_1(x) = x\\B_2(x) = x^2+x\\B_3(x) = x^3+3x^2+x\end{aligned}\end{align} \]If \(x = 1\) or
bell()is called with only one argument, it gives the \(n\)-th Bell number \(B_n\), which is the number of partitions of a set with \(n\) elements. By setting the precision to at least \(\log_{10} B_n\) digits,bell()provides fast calculation of exact Bell numbers.In general,
bell()computes\[B_n(x) = e^{-x} \left(\mathrm{sinc}(\pi n) + E_n(x)\right)\]where \(E_n(x)\) is the generalized exponential function implemented by
polyexp(). This is an extension of Dobinski’s formula [1], where the modification is the sinc term ensuring that \(B_n(x)\) is continuous in \(n\);bell()can thus be evaluated, differentiated, etc for arbitrary complex arguments.Examples
Simple evaluations:
>>> from mpmath import * >>> mp.dps = 25; mp.pretty = True >>> bell(0, 2.5) 1.0 >>> bell(1, 2.5) 2.5 >>> bell(2, 2.5) 8.75
Evaluation for arbitrary complex arguments:
>>> bell(5.75+1j, 2-3j) (-10767.71345136587098445143 - 15449.55065599872579097221j)
The first few Bell polynomials:
>>> for k in range(7): ... nprint(taylor(lambda x: bell(k,x), 0, k)) ... [1.0] [0.0, 1.0] [0.0, 1.0, 1.0] [0.0, 1.0, 3.0, 1.0] [0.0, 1.0, 7.0, 6.0, 1.0] [0.0, 1.0, 15.0, 25.0, 10.0, 1.0] [0.0, 1.0, 31.0, 90.0, 65.0, 15.0, 1.0]
The first few Bell numbers and complementary Bell numbers:
>>> [int(bell(k)) for k in range(10)] [1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147] >>> [int(bell(k,-1)) for k in range(10)] [1, -1, 0, 1, 1, -2, -9, -9, 50, 267]
Large Bell numbers:
>>> mp.dps = 50 >>> bell(50) 185724268771078270438257767181908917499221852770.0 >>> bell(50,-1) -29113173035759403920216141265491160286912.0
Some even larger values:
>>> mp.dps = 25 >>> bell(1000,-1) -1.237132026969293954162816e+1869 >>> bell(1000) 2.989901335682408421480422e+1927 >>> bell(1000,2) 6.591553486811969380442171e+1987 >>> bell(1000,100.5) 9.101014101401543575679639e+2529
A determinant identity satisfied by Bell numbers:
>>> mp.dps = 15 >>> N = 8 >>> det([[bell(k+j) for j in range(N)] for k in range(N)]) 125411328000.0 >>> superfac(N-1) 125411328000.0
References
Stirling numbers¶
stirling1()¶
-
mpmath.stirling1(n, k, exact=False)¶ 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 function:
>>> from mpmath 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
The matrices of Stirling numbers of first and second kind are inverses of each other:
>>> A = matrix(5, 5); B = matrix(5, 5) >>> for n in range(5): ... for k in range(5): ... A[n,k] = stirling1(n,k) ... B[n,k] = stirling2(n,k) ... >>> A * B [1.0 0.0 0.0 0.0 0.0] [0.0 1.0 0.0 0.0 0.0] [0.0 0.0 1.0 0.0 0.0] [0.0 0.0 0.0 1.0 0.0] [0.0 0.0 0.0 0.0 1.0]
Pass
exact=Trueto obtain exact values of Stirling numbers as integers:>>> stirling1(42, 5) -2.864498971768501633736628e+50 >>> print(stirling1(42, 5, exact=True)) -286449897176850163373662803014001546235808317440000
stirling2()¶
-
mpmath.stirling2(n, k, exact=False)¶ 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 function:
>>> from mpmath 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=Trueto obtain exact values of Stirling numbers as integers:>>> stirling2(52, 10) 2.641822121003543906807485e+45 >>> print(stirling2(52, 10, exact=True)) 2641822121003543906807485307053638921722527655
Prime counting functions¶
primepi2()¶
-
mpmath.primepi2(x)¶ Returns an interval (as an
mpiinstance) providing bounds for the value of the prime counting function \(\pi(x)\). For small \(x\),primepi2()returns an exact interval based on the output ofprimepi(). 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.
Examples
Exact values of the prime counting function for small \(x\):
>>> from mpmath import * >>> mp.dps = 15; mp.pretty = True >>> iv.dps = 15; iv.pretty = True >>> primepi2(10) [4.0, 4.0] >>> primepi2(100) [25.0, 25.0] >>> primepi2(1000) [168.0, 168.0]
Loose intervals are generated for moderately large \(x\):
>>> primepi2(10000), primepi(10000) ([1209.0, 1283.0], 1229) >>> primepi2(50000), primepi(50000) ([5070.0, 5263.0], 5133)
As \(x\) increases, the absolute error gets worse while the relative error improves. The exact value of \(\pi(10^{23})\) is 1925320391606803968923, and
primepi2()gives 9 significant digits:>>> p = primepi2(10**23) >>> p [1.9253203909477020467e+21, 1.925320392280406229e+21] >>> mpf(p.delta) / mpf(p.a) 6.9219865355293e-10
A more precise, nonrigorous estimate for \(\pi(x)\) can be obtained using the Riemann R function (
riemannr()). For large enough \(x\), the value returned byprimepi2()essentially amounts to a small perturbation of the value returned byriemannr():>>> primepi2(10**100) [4.3619719871407024816e+97, 4.3619719871407032404e+97] >>> riemannr(10**100) 4.3619719871407e+97
riemannr()¶
-
mpmath.riemannr(x)¶ Evaluates the Riemann R function, a smooth approximation of the prime counting function \(\pi(x)\) (see
primepi()). 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.
The Riemann R function gives a very accurate approximation of the prime counting function. For example, it is wrong by at most 2 for \(x < 1000\), and for \(x = 10^9\) differs from the exact value of \(\pi(x)\) by 79, or less than two parts in a million. It is about 10 times more accurate than the logarithmic integral estimate (see
li()), which however is even faster to evaluate. It is orders of magnitude more accurate than the extremely fast \(x/\log x\) estimate.Examples
For small arguments, the Riemann R function almost exactly gives the prime counting function if rounded to the nearest integer:
>>> from mpmath import * >>> mp.dps = 15; mp.pretty = True >>> primepi(50), riemannr(50) (15, 14.9757023241462) >>> max(abs(primepi(n)-int(round(riemannr(n)))) for n in range(100)) 1 >>> max(abs(primepi(n)-int(round(riemannr(n)))) for n in range(300)) 2
The Riemann R function can be evaluated for arguments far too large for exact determination of \(\pi(x)\) to be computationally feasible with any presently known algorithm:
>>> riemannr(10**30) 1.46923988977204e+28 >>> riemannr(10**100) 4.3619719871407e+97 >>> riemannr(10**1000) 4.3448325764012e+996
A comparison of the Riemann R function and logarithmic integral estimates for \(\pi(x)\) using exact values of \(\pi(10^n)\) up to \(n = 9\). The fractional error is shown in parentheses:
>>> exact = [4,25,168,1229,9592,78498,664579,5761455,50847534] >>> for n, p in enumerate(exact): ... n += 1 ... r, l = riemannr(10**n), li(10**n) ... rerr, lerr = nstr((r-p)/p,3), nstr((l-p)/p,3) ... print("%i %i %s(%s) %s(%s)" % (n, p, r, rerr, l, lerr)) ... 1 4 4.56458314100509(0.141) 6.1655995047873(0.541) 2 25 25.6616332669242(0.0265) 30.1261415840796(0.205) 3 168 168.359446281167(0.00214) 177.609657990152(0.0572) 4 1229 1226.93121834343(-0.00168) 1246.13721589939(0.0139) 5 9592 9587.43173884197(-0.000476) 9629.8090010508(0.00394) 6 78498 78527.3994291277(0.000375) 78627.5491594622(0.00165) 7 664579 664667.447564748(0.000133) 664918.405048569(0.000511) 8 5761455 5761551.86732017(1.68e-5) 5762209.37544803(0.000131) 9 50847534 50847455.4277214(-1.55e-6) 50849234.9570018(3.35e-5)
The derivative of the Riemann R function gives the approximate probability for a number of magnitude \(x\) to be prime:
>>> diff(riemannr, 1000) 0.141903028110784 >>> mpf(primepi(1050) - primepi(950)) / 100 0.15
Evaluation is supported for arbitrary arguments and at arbitrary precision:
>>> mp.dps = 30 >>> riemannr(7.5) 3.72934743264966261918857135136 >>> riemannr(-4+2j) (-0.551002208155486427591793957644 + 2.16966398138119450043195899746j)
Cyclotomic polynomials¶
cyclotomic()¶
-
mpmath.cyclotomic(n, x)¶ 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 cyclotomic polynomials are integer polynomials, the first of which can be written explicitly as
\[ \begin{align}\begin{aligned}\Phi_0(x) = 1\\\Phi_1(x) = x - 1\\\Phi_2(x) = x + 1\\\Phi_3(x) = x^3 + x^2 + 1\\\Phi_4(x) = x^2 + 1\\\Phi_5(x) = x^4 + x^3 + x^2 + x + 1\\\Phi_6(x) = x^2 - x + 1\end{aligned}\end{align} \]Examples
The coefficients of low-order cyclotomic polynomials can be recovered using Taylor expansion:
>>> from mpmath import * >>> mp.dps = 15; mp.pretty = True >>> for n in range(9): ... p = chop(taylor(lambda x: cyclotomic(n,x), 0, 10)) ... print("%s %s" % (n, nstr(p[:10+1-p[::-1].index(1)]))) ... 0 [1.0] 1 [-1.0, 1.0] 2 [1.0, 1.0] 3 [1.0, 1.0, 1.0] 4 [1.0, 0.0, 1.0] 5 [1.0, 1.0, 1.0, 1.0, 1.0] 6 [1.0, -1.0, 1.0] 7 [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] 8 [1.0, 0.0, 0.0, 0.0, 1.0]
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)
Up to permutation, the roots of a given cyclotomic polynomial can be checked to agree with the list of primitive roots:
>>> p = taylor(lambda x: cyclotomic(6,x), 0, 6)[:3] >>> for r in polyroots(p[::-1]): ... print(r) ... (0.5 - 0.8660254037844386467637232j) (0.5 + 0.8660254037844386467637232j) >>> >>> for r in unitroots(6, primitive=True): ... print(r) ... (0.5 + 0.8660254037844386467637232j) (0.5 - 0.8660254037844386467637232j)
Arithmetic functions¶
mangoldt()¶
-
mpmath.mangoldt(n)¶ Evaluates the von Mangoldt function \(\Lambda(n) = \log p\) if \(n = p^k\) a power of a prime, and \(\Lambda(n) = 0\) otherwise.
Examples
>>> from mpmath 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