# 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

mpmath.fibonacci(x, **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 for `fibonacci()`.

  Basic examples

  Some small Fibonacci numbers are:

      >>> from mpmath import mp, fibonacci, fib, pi, findroot, nsum, sqrt, inf
      >>> 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

  1. [Weisstein] http://mathworld.wolfram.com/FibonacciNumber.html

## Bernoulli numbers and polynomials

mpmath.bernoulli(n, plus=False)

  Computes the nth Bernoulli number, \(B_n\), for any integer \(n \ge 0\).

  The Bernoulli numbers are rational numbers, but this function returns a floating-point approximation. To obtain an exact fraction, use `bernfrac()` instead.

  Optional `plus` flag (default: False) control the sign choice of the \(B_1\) value (default: \(-0.5\)).

  Examples

  Numerical values of the first few Bernoulli numbers:

      >>> from mpmath import mp, bernoulli, pi, fac, zeta
      >>> mp.pretty = True
      >>> for n in range(15):
      ...     print("%s %s" % (n, bernoulli(n)))
      ...
      0 1.0
      1 -0.5
      2 0.166666666666667
      3 0.0
      4 -0.0333333333333333
      5 0.0
      6 0.0238095238095238
      7 0.0
      8 -0.0333333333333333
      9 0.0
      10 0.0757575757575758
      11 0.0
      12 -0.253113553113553
      13 0.0
      14 1.16666666666667

  Bernoulli numbers can be approximated with arbitrary precision:

      >>> mp.dps = 50
      >>> bernoulli(100)
      -2.8382249570693706959264156336481764738284680928013e+78

  Arbitrarily large \(n\) are supported:

      >>> mp.dps = 15
      >>> bernoulli(10**20 + 2)
      3.09136296657021e+1876752564973863312327

  The Bernoulli numbers are related to the Riemann zeta function at integer arguments:

      >>> -bernoulli(8) * (2*pi)**8 / (2*fac(8))
      1.00407735619794
      >>> zeta(8)
      1.00407735619794

  Algorithm

  For small \(n\) (\(n < 3000\)) `bernoulli()` uses a recurrence formula due to Ramanujan. All results in this range are cached, so sequential computation of small Bernoulli numbers is guaranteed to be fast.

  For larger \(n\), \(B_n\) is evaluated in terms of the Riemann zeta function.

  References

  1. [Wikipedia] https://en.wikipedia.org/wiki/Bernoulli_number

mpmath.bernfrac(n, plus=False)

  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.

  Optional `plus` flag (default: False) control the sign choice of the \(B_1\) value (default: \((-1, 2)\)).

  Examples

  The first few Bernoulli numbers are exactly:

      >>> from mpmath import bernfrac, mp, mpf, bernoulli
      >>> 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\):

      >>> import sys
      >>> if hasattr(sys, 'set_int_max_str_digits'):
      ...     sys.set_int_max_str_digits(30000)
      >>> del sys
      >>> 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

  1. [Weisstein] http://mathworld.wolfram.com/vonStaudt-ClausenTheorem.html

  2. [Bernoulli]

  3. [Wikipedia] https://en.wikipedia.org/wiki/Bernoulli_number

mpmath.bernpoly(n, z)

  Evaluates the Bernoulli polynomial \(B_n(z)\).

  The first few Bernoulli polynomials are:

      >>> from mpmath import mp, nprint, chop, taylor, bernpoly, bernoulli
      >>> 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

mpmath.eulernum(n, exact=False)

  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, eulernum, chop, diffs, sech, sqrt, pi, e
      >>> 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=True` to 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

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, chop, taylor, eulerpoly, inf, eulernum
      >>> 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

mpmath.bell(n, x=1)

  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, bell, nprint, taylor, det, superfac
      >>> 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

  1. [Weisstein] http://mathworld.wolfram.com/DobinskisFormula.html

## Stirling numbers

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, taylor, ff, stirling1, matrix, stirling2
      >>> 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=True` to obtain exact values of Stirling numbers as integers:

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

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, stirling2, taylor, ff
      >>> 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

## Prime counting functions

mpmath.primepi(x)

  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 function:

      >>> from mpmath import primepi
      >>> [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.primepi2(x)

  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.

  Examples

  Exact values of the prime counting function for small \(x\):

      >>> from mpmath import mp, iv, primepi2, primepi, mpf, riemannr
      >>> mp.pretty = True
      >>> 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 by `primepi2()` essentially amounts to a small perturbation of the value returned by `riemannr()`:

      >>> primepi2(10**100)
      [4.3619719871407024816e+97, 4.3619719871407032404e+97]
      >>> riemannr(10**100)
      4.3619719871407e+97

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, primepi, riemannr, nstr, li, diff, mpf
      >>> 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

mpmath.cyclotomic(n, z)

  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, chop, taylor, cyclotomic, nstr, fprod,
      ...                     unitroots, polyroots)
      >>> 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, asc=True):
      ...     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

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, mangoldt, fsum
      >>> 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