Powers and logarithms

Nth roots

sqrt()

mpmath.sqrt(x, **kwargs)

sqrt(x) gives the principal square root of \(x\), \(\sqrt x\). For positive real numbers, the principal root is simply the positive square root. For arbitrary complex numbers, the principal square root is defined to satisfy \(\sqrt x = \exp(\log(x)/2)\). The function thus has a branch cut along the negative half real axis.

For all mpmath numbers x, calling sqrt(x) is equivalent to performing x**0.5.

Examples

Basic examples and limits:

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> sqrt(10)
3.16227766016838
>>> sqrt(100)
10.0
>>> sqrt(-4)
(0.0 + 2.0j)
>>> sqrt(1+1j)
(1.09868411346781 + 0.455089860562227j)
>>> sqrt(inf)
+inf

Square root evaluation is fast at huge precision:

>>> mp.dps = 50000
>>> a = sqrt(3)
>>> str(a)[-10:]
'9329332815'

mpmath.iv.sqrt() supports interval arguments:

>>> iv.dps = 15; iv.pretty = True
>>> iv.sqrt([16,100])
[4.0, 10.0]
>>> iv.sqrt(2)
[1.4142135623730949234, 1.4142135623730951455]
>>> iv.sqrt(2) ** 2
[1.9999999999999995559, 2.0000000000000004441]

hypot()

mpmath.hypot(x, y)

Computes the Euclidean norm of the vector \((x, y)\), equal to \(\sqrt{x^2 + y^2}\). Both \(x\) and \(y\) must be real.

cbrt()

mpmath.cbrt(x, **kwargs)

cbrt(x) computes the cube root of \(x\), \(x^{1/3}\). This function is faster and more accurate than raising to a floating-point fraction:

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = False
>>> 125**(mpf(1)/3)
mpf('4.9999999999999991')
>>> cbrt(125)
mpf('5.0')

Every nonzero complex number has three cube roots. This function returns the cube root defined by \(\exp(\log(x)/3)\) where the principal branch of the natural logarithm is used. Note that this does not give a real cube root for negative real numbers:

>>> mp.pretty = True
>>> cbrt(-1)
(0.5 + 0.866025403784439j)

root()

mpmath.root(z, n, k=0)

root(z, n, k=0) computes an \(n\)-th root of \(z\), i.e. returns a number \(r\) that (up to possible approximation error) satisfies \(r^n = z\). (nthroot is available as an alias for root.)

Every complex number \(z \ne 0\) has \(n\) distinct \(n\)-th roots, which are equidistant points on a circle with radius \(|z|^{1/n}\), centered around the origin. A specific root may be selected using the optional index \(k\). The roots are indexed counterclockwise, starting with \(k = 0\) for the root closest to the positive real half-axis.

The \(k = 0\) root is the so-called principal \(n\)-th root, often denoted by \(\sqrt[n]{z}\) or \(z^{1/n}\), and also given by \(\exp(\log(z) / n)\). If \(z\) is a positive real number, the principal root is just the unique positive \(n\)-th root of \(z\). Under some circumstances, non-principal real roots exist: for positive real \(z\), \(n\) even, there is a negative root given by \(k = n/2\); for negative real \(z\), \(n\) odd, there is a negative root given by \(k = (n-1)/2\).

To obtain all roots with a simple expression, use [root(z,n,k) for k in range(n)].

An important special case, root(1, n, k) returns the \(k\)-th \(n\)-th root of unity, \(\zeta_k = e^{2 \pi i k / n}\). Alternatively, unitroots() provides a slightly more convenient way to obtain the roots of unity, including the option to compute only the primitive roots of unity.

Both \(k\) and \(n\) should be integers; \(k\) outside of range(n) will be reduced modulo \(n\). If \(n\) is negative, \(x^{-1/n} = 1/x^{1/n}\) (or the equivalent reciprocal for a non-principal root with \(k \ne 0\)) is computed.

root() is implemented to use Newton’s method for small \(n\). At high precision, this makes \(x^{1/n}\) not much more expensive than the regular exponentiation, \(x^n\). For very large \(n\), nthroot() falls back to use the exponential function.

Examples

nthroot()/root() is faster and more accurate than raising to a floating-point fraction:

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = False
>>> 16807 ** (mpf(1)/5)
mpf('7.0000000000000009')
>>> root(16807, 5)
mpf('7.0')
>>> nthroot(16807, 5)    # Alias
mpf('7.0')

A high-precision root:

>>> mp.dps = 50; mp.pretty = True
>>> nthroot(10, 5)
1.584893192461113485202101373391507013269442133825
>>> nthroot(10, 5) ** 5
10.0

Computing principal and non-principal square and cube roots:

>>> mp.dps = 15
>>> root(10, 2)
3.16227766016838
>>> root(10, 2, 1)
-3.16227766016838
>>> root(-10, 3)
(1.07721734501594 + 1.86579517236206j)
>>> root(-10, 3, 1)
-2.15443469003188
>>> root(-10, 3, 2)
(1.07721734501594 - 1.86579517236206j)

All the 7th roots of a complex number:

>>> for r in [root(3+4j, 7, k) for k in range(7)]:
...     print("%s %s" % (r, r**7))
...
(1.24747270589553 + 0.166227124177353j) (3.0 + 4.0j)
(0.647824911301003 + 1.07895435170559j) (3.0 + 4.0j)
(-0.439648254723098 + 1.17920694574172j) (3.0 + 4.0j)
(-1.19605731775069 + 0.391492658196305j) (3.0 + 4.0j)
(-1.05181082538903 - 0.691023585965793j) (3.0 + 4.0j)
(-0.115529328478668 - 1.25318497558335j) (3.0 + 4.0j)
(0.907748109144957 - 0.871672518271819j) (3.0 + 4.0j)

Cube roots of unity:

>>> for k in range(3): print(root(1, 3, k))
...
1.0
(-0.5 + 0.866025403784439j)
(-0.5 - 0.866025403784439j)

Some exact high order roots:

>>> root(75**210, 105)
5625.0
>>> root(1, 128, 96)
(0.0 - 1.0j)
>>> root(4**128, 128, 96)
(0.0 - 4.0j)

unitroots()

mpmath.unitroots(n, primitive=False)

unitroots(n) returns \(\zeta_0, \zeta_1, \ldots, \zeta_{n-1}\), all the distinct \(n\)-th roots of unity, as a list. If the option primitive=True is passed, only the primitive roots are returned.

Every \(n\)-th root of unity satisfies \((\zeta_k)^n = 1\). There are \(n\) distinct roots for each \(n\) (\(\zeta_k\) and \(\zeta_j\) are the same when \(k = j \pmod n\)), which form a regular polygon with vertices on the unit circle. They are ordered counterclockwise with increasing \(k\), starting with \(\zeta_0 = 1\).

Examples

The roots of unity up to \(n = 4\):

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> nprint(unitroots(1))
[1.0]
>>> nprint(unitroots(2))
[1.0, -1.0]
>>> nprint(unitroots(3))
[1.0, (-0.5 + 0.866025j), (-0.5 - 0.866025j)]
>>> nprint(unitroots(4))
[1.0, (0.0 + 1.0j), -1.0, (0.0 - 1.0j)]

Roots of unity form a geometric series that sums to 0:

>>> mp.dps = 50
>>> chop(fsum(unitroots(25)))
0.0

Primitive roots up to \(n = 4\):

>>> mp.dps = 15
>>> nprint(unitroots(1, primitive=True))
[1.0]
>>> nprint(unitroots(2, primitive=True))
[-1.0]
>>> nprint(unitroots(3, primitive=True))
[(-0.5 + 0.866025j), (-0.5 - 0.866025j)]
>>> nprint(unitroots(4, primitive=True))
[(0.0 + 1.0j), (0.0 - 1.0j)]

There are only four primitive 12th roots:

>>> nprint(unitroots(12, primitive=True))
[(0.866025 + 0.5j), (-0.866025 + 0.5j), (-0.866025 - 0.5j), (0.866025 - 0.5j)]

The \(n\)-th roots of unity form a group, the cyclic group of order \(n\). Any primitive root \(r\) is a generator for this group, meaning that \(r^0, r^1, \ldots, r^{n-1}\) gives the whole set of unit roots (in some permuted order):

>>> for r in unitroots(6): print(r)
...
1.0
(0.5 + 0.866025403784439j)
(-0.5 + 0.866025403784439j)
-1.0
(-0.5 - 0.866025403784439j)
(0.5 - 0.866025403784439j)
>>> r = unitroots(6, primitive=True)[1]
>>> for k in range(6): print(chop(r**k))
...
1.0
(0.5 - 0.866025403784439j)
(-0.5 - 0.866025403784439j)
-1.0
(-0.5 + 0.866025403784438j)
(0.5 + 0.866025403784438j)

The number of primitive roots equals the Euler totient function \(\phi(n)\):

>>> [len(unitroots(n, primitive=True)) for n in range(1,20)]
[1, 1, 2, 2, 4, 2, 6, 4, 6, 4, 10, 4, 12, 6, 8, 8, 16, 6, 18]

Exponentiation

exp()

mpmath.exp(x, **kwargs)

Computes the exponential function,

\[\exp(x) = e^x = \sum_{k=0}^{\infty} \frac{x^k}{k!}.\]

For complex numbers, the exponential function also satisfies

\[\exp(x+yi) = e^x (\cos y + i \sin y).\]

Basic examples

Some values of the exponential function:

>>> from mpmath import *
>>> mp.dps = 25; mp.pretty = True
>>> exp(0)
1.0
>>> exp(1)
2.718281828459045235360287
>>> exp(-1)
0.3678794411714423215955238
>>> exp(inf)
+inf
>>> exp(-inf)
0.0

Arguments can be arbitrarily large:

>>> exp(10000)
8.806818225662921587261496e+4342
>>> exp(-10000)
1.135483865314736098540939e-4343

Evaluation is supported for interval arguments via mpmath.iv.exp():

>>> iv.dps = 25; iv.pretty = True
>>> iv.exp([-inf,0])
[0.0, 1.0]
>>> iv.exp([0,1])
[1.0, 2.71828182845904523536028749558]

The exponential function can be evaluated efficiently to arbitrary precision:

>>> mp.dps = 10000
>>> exp(pi)  
23.140692632779269005729...8984304016040616

Functional properties

Numerical verification of Euler’s identity for the complex exponential function:

>>> mp.dps = 15
>>> exp(j*pi)+1
(0.0 + 1.22464679914735e-16j)
>>> chop(exp(j*pi)+1)
0.0

This recovers the coefficients (reciprocal factorials) in the Maclaurin series expansion of exp:

>>> nprint(taylor(exp, 0, 5))
[1.0, 1.0, 0.5, 0.166667, 0.0416667, 0.00833333]

The exponential function is its own derivative and antiderivative:

>>> exp(pi)
23.1406926327793
>>> diff(exp, pi)
23.1406926327793
>>> quad(exp, [-inf, pi])
23.1406926327793

The exponential function can be evaluated using various methods, including direct summation of the series, limits, and solving the defining differential equation:

>>> nsum(lambda k: pi**k/fac(k), [0,inf])
23.1406926327793
>>> limit(lambda k: (1+pi/k)**k, inf)
23.1406926327793
>>> odefun(lambda t, x: x, 0, 1)(pi)
23.1406926327793

power()

mpmath.power(x, y)

Converts \(x\) and \(y\) to mpmath numbers and evaluates \(x^y = \exp(y \log(x))\):

>>> from mpmath import *
>>> mp.dps = 30; mp.pretty = True
>>> power(2, 0.5)
1.41421356237309504880168872421

This shows the leading few digits of a large Mersenne prime (performing the exact calculation 2**43112609-1 and displaying the result in Python would be very slow):

>>> power(2, 43112609)-1
3.16470269330255923143453723949e+12978188

expj()

mpmath.expj(x, **kwargs)

Convenience function for computing \(e^{ix}\):

>>> from mpmath import *
>>> mp.dps = 25; mp.pretty = True
>>> expj(0)
(1.0 + 0.0j)
>>> expj(-1)
(0.5403023058681397174009366 - 0.8414709848078965066525023j)
>>> expj(j)
(0.3678794411714423215955238 + 0.0j)
>>> expj(1+j)
(0.1987661103464129406288032 + 0.3095598756531121984439128j)

expjpi()

mpmath.expjpi(x, **kwargs)

Convenience function for computing \(e^{i \pi x}\). Evaluation is accurate near zeros (see also cospi(), sinpi()):

>>> from mpmath import *
>>> mp.dps = 25; mp.pretty = True
>>> expjpi(0)
(1.0 + 0.0j)
>>> expjpi(1)
(-1.0 + 0.0j)
>>> expjpi(0.5)
(0.0 + 1.0j)
>>> expjpi(-1)
(-1.0 + 0.0j)
>>> expjpi(j)
(0.04321391826377224977441774 + 0.0j)
>>> expjpi(1+j)
(-0.04321391826377224977441774 + 0.0j)

expm1()

mpmath.expm1(x)

Computes \(e^x - 1\), accurately for small \(x\).

Unlike the expression exp(x) - 1, expm1(x) does not suffer from potentially catastrophic cancellation:

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> exp(1e-10)-1; print(expm1(1e-10))
1.00000008274037e-10
1.00000000005e-10
>>> exp(1e-20)-1; print(expm1(1e-20))
0.0
1.0e-20
>>> 1/(exp(1e-20)-1)
Traceback (most recent call last):
  ...
ZeroDivisionError
>>> 1/expm1(1e-20)
1.0e+20

Evaluation works for extremely tiny values:

>>> expm1(0)
0.0
>>> expm1('1e-10000000')
1.0e-10000000

powm1()

mpmath.powm1(x, y)

Computes \(x^y - 1\), accurately when \(x^y\) is very close to 1.

This avoids potentially catastrophic cancellation:

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> power(0.99999995, 1e-10) - 1
0.0
>>> powm1(0.99999995, 1e-10)
-5.00000012791934e-18

Powers exactly equal to 1, and only those powers, yield 0 exactly:

>>> powm1(-j, 4)
(0.0 + 0.0j)
>>> powm1(3, 0)
0.0
>>> powm1(fadd(-1, 1e-100, exact=True), 4)
-4.0e-100

Evaluation works for extremely tiny \(y\):

>>> powm1(2, '1e-100000')
6.93147180559945e-100001
>>> powm1(j, '1e-1000')
(-1.23370055013617e-2000 + 1.5707963267949e-1000j)

Logarithms

log()

mpmath.log(x, b=None)

Computes the base-\(b\) logarithm of \(x\), \(\log_b(x)\). If \(b\) is unspecified, log() computes the natural (base \(e\)) logarithm and is equivalent to ln(). In general, the base \(b\) logarithm is defined in terms of the natural logarithm as \(\log_b(x) = \ln(x)/\ln(b)\).

By convention, we take \(\log(0) = -\infty\).

The natural logarithm is real if \(x > 0\) and complex if \(x < 0\) or if \(x\) is complex. The principal branch of the complex logarithm is used, meaning that \(\Im(\ln(x)) = -\pi < \arg(x) \le \pi\).

Examples

Some basic values and limits:

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> log(1)
0.0
>>> log(2)
0.693147180559945
>>> log(1000,10)
3.0
>>> log(4, 16)
0.5
>>> log(j)
(0.0 + 1.5707963267949j)
>>> log(-1)
(0.0 + 3.14159265358979j)
>>> log(0)
-inf
>>> log(inf)
+inf

The natural logarithm is the antiderivative of \(1/x\):

>>> quad(lambda x: 1/x, [1, 5])
1.6094379124341
>>> log(5)
1.6094379124341
>>> diff(log, 10)
0.1

The Taylor series expansion of the natural logarithm around \(x = 1\) has coefficients \((-1)^{n+1}/n\):

>>> nprint(taylor(log, 1, 7))
[0.0, 1.0, -0.5, 0.333333, -0.25, 0.2, -0.166667, 0.142857]

log() supports arbitrary precision evaluation:

>>> mp.dps = 50
>>> log(pi)
1.1447298858494001741434273513530587116472948129153
>>> log(pi, pi**3)
0.33333333333333333333333333333333333333333333333333
>>> mp.dps = 25
>>> log(3+4j)
(1.609437912434100374600759 + 0.9272952180016122324285125j)

ln()

mpmath.ln(x, **kwargs)

Computes the base-\(b\) logarithm of \(x\), \(\log_b(x)\). If \(b\) is unspecified, log() computes the natural (base \(e\)) logarithm and is equivalent to ln(). In general, the base \(b\) logarithm is defined in terms of the natural logarithm as \(\log_b(x) = \ln(x)/\ln(b)\).

By convention, we take \(\log(0) = -\infty\).

The natural logarithm is real if \(x > 0\) and complex if \(x < 0\) or if \(x\) is complex. The principal branch of the complex logarithm is used, meaning that \(\Im(\ln(x)) = -\pi < \arg(x) \le \pi\).

Examples

Some basic values and limits:

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> log(1)
0.0
>>> log(2)
0.693147180559945
>>> log(1000,10)
3.0
>>> log(4, 16)
0.5
>>> log(j)
(0.0 + 1.5707963267949j)
>>> log(-1)
(0.0 + 3.14159265358979j)
>>> log(0)
-inf
>>> log(inf)
+inf

The natural logarithm is the antiderivative of \(1/x\):

>>> quad(lambda x: 1/x, [1, 5])
1.6094379124341
>>> log(5)
1.6094379124341
>>> diff(log, 10)
0.1

The Taylor series expansion of the natural logarithm around \(x = 1\) has coefficients \((-1)^{n+1}/n\):

>>> nprint(taylor(log, 1, 7))
[0.0, 1.0, -0.5, 0.333333, -0.25, 0.2, -0.166667, 0.142857]

log() supports arbitrary precision evaluation:

>>> mp.dps = 50
>>> log(pi)
1.1447298858494001741434273513530587116472948129153
>>> log(pi, pi**3)
0.33333333333333333333333333333333333333333333333333
>>> mp.dps = 25
>>> log(3+4j)
(1.609437912434100374600759 + 0.9272952180016122324285125j)

log10()

mpmath.log10(x)

Computes the base-10 logarithm of \(x\), \(\log_{10}(x)\). log10(x) is equivalent to log(x, 10).

log1p()

mpmath.log1p(x)

Computes \(\log(1+x)\), accurately for small \(x\).

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> log(1+1e-10); print(mp.log1p(1e-10))
1.00000008269037e-10
9.9999999995e-11
>>> mp.log1p(1e-100j)
(5.0e-201 + 1.0e-100j)
>>> mp.log1p(0)
0.0

Lambert W function

lambertw()

mpmath.lambertw(z, k=0)

The Lambert W function \(W(z)\) is defined as the inverse function of \(w \exp(w)\). In other words, the value of \(W(z)\) is such that \(z = W(z) \exp(W(z))\) for any complex number \(z\).

The Lambert W function is a multivalued function with infinitely many branches \(W_k(z)\), indexed by \(k \in \mathbb{Z}\). Each branch gives a different solution \(w\) of the equation \(z = w \exp(w)\). All branches are supported by lambertw():

• lambertw(z) gives the principal solution (branch 0)

• lambertw(z, k) gives the solution on branch \(k\)

The Lambert W function has two partially real branches: the principal branch (\(k = 0\)) is real for real \(z > -1/e\), and the \(k = -1\) branch is real for \(-1/e < z < 0\). All branches except \(k = 0\) have a logarithmic singularity at \(z = 0\).

The definition, implementation and choice of branches is based on [Corless].

Plots

# Branches 0 and -1 of the Lambert W function
plot([lambertw, lambda x: lambertw(x,-1)], [-2,2], [-5,2], points=2000)

# Principal branch of the Lambert W function W(z)
cplot(lambertw, [-1,1], [-1,1], points=50000)

Basic examples

The Lambert W function is the inverse of \(w \exp(w)\):

>>> from mpmath import *
>>> mp.dps = 25; mp.pretty = True
>>> w = lambertw(1)
>>> w
0.5671432904097838729999687
>>> w*exp(w)
1.0

Any branch gives a valid inverse:

>>> w = lambertw(1, k=3)
>>> w
(-2.853581755409037807206819 + 17.11353553941214591260783j)
>>> w = lambertw(1, k=25)
>>> w
(-5.047020464221569709378686 + 155.4763860949415867162066j)
>>> chop(w*exp(w))
1.0

Applications to equation-solving

The Lambert W function may be used to solve various kinds of equations, such as finding the value of the infinite power tower \(z^{z^{z^{\ldots}}}\):

>>> def tower(z, n):
...     if n == 0:
...         return z
...     return z ** tower(z, n-1)
...
>>> tower(mpf(0.5), 100)
0.6411857445049859844862005
>>> -lambertw(-log(0.5))/log(0.5)
0.6411857445049859844862005

Properties

The Lambert W function grows roughly like the natural logarithm for large arguments:

>>> lambertw(1000); log(1000)
5.249602852401596227126056
6.907755278982137052053974
>>> lambertw(10**100); log(10**100)
224.8431064451185015393731
230.2585092994045684017991

The principal branch of the Lambert W function has a rational Taylor series expansion around \(z = 0\):

>>> nprint(taylor(lambertw, 0, 6), 10)
[0.0, 1.0, -1.0, 1.5, -2.666666667, 5.208333333, -10.8]

Some special values and limits are:

>>> lambertw(0)
0.0
>>> lambertw(1)
0.5671432904097838729999687
>>> lambertw(e)
1.0
>>> lambertw(inf)
+inf
>>> lambertw(0, k=-1)
-inf
>>> lambertw(0, k=3)
-inf
>>> lambertw(inf, k=2)
(+inf + 12.56637061435917295385057j)
>>> lambertw(inf, k=3)
(+inf + 18.84955592153875943077586j)
>>> lambertw(-inf, k=3)
(+inf + 21.9911485751285526692385j)

The \(k = 0\) and \(k = -1\) branches join at \(z = -1/e\) where \(W(z) = -1\) for both branches. Since \(-1/e\) can only be represented approximately with binary floating-point numbers, evaluating the Lambert W function at this point only gives \(-1\) approximately:

>>> lambertw(-1/e, 0)
-0.9999999999998371330228251
>>> lambertw(-1/e, -1)
-1.000000000000162866977175

If \(-1/e\) happens to round in the negative direction, there might be a small imaginary part:

>>> mp.dps = 15
>>> lambertw(-1/e)
(-1.0 + 8.22007971483662e-9j)
>>> lambertw(-1/e+eps)
-0.999999966242188

References

1. [Corless]

Arithmetic-geometric mean

agm()

mpmath.agm(a, b=1)

agm(a, b) computes the arithmetic-geometric mean of \(a\) and \(b\), defined as the limit of the following iteration:

\[ \begin{align}\begin{aligned}a_0 = a\\b_0 = b\\a_{n+1} = \frac{a_n+b_n}{2}\\b_{n+1} = \sqrt{a_n b_n}\end{aligned}\end{align} \]

This function can be called with a single argument, computing \(\mathrm{agm}(a,1) = \mathrm{agm}(1,a)\).

Examples

It is a well-known theorem that the geometric mean of two distinct positive numbers is less than the arithmetic mean. It follows that the arithmetic-geometric mean lies between the two means:

>>> from mpmath import *
>>> mp.dps = 15; mp.pretty = True
>>> a = mpf(3)
>>> b = mpf(4)
>>> sqrt(a*b)
3.46410161513775
>>> agm(a,b)
3.48202767635957
>>> (a+b)/2
3.5

The arithmetic-geometric mean is scale-invariant:

>>> agm(10*e, 10*pi)
29.261085515723
>>> 10*agm(e, pi)
29.261085515723

As an order-of-magnitude estimate, \(\mathrm{agm}(1,x) \approx x\) for large \(x\):

>>> agm(10**10)
643448704.760133
>>> agm(10**50)
1.34814309345871e+48

For tiny \(x\), \(\mathrm{agm}(1,x) \approx -\pi/(2 \log(x/4))\):

>>> agm('0.01')
0.262166887202249
>>> -pi/2/log('0.0025')
0.262172347753122

The arithmetic-geometric mean can also be computed for complex numbers:

>>> agm(3, 2+j)
(2.51055133276184 + 0.547394054060638j)

The AGM iteration converges very quickly (each step doubles the number of correct digits), so agm() supports efficient high-precision evaluation:

>>> mp.dps = 10000
>>> a = agm(1,2)
>>> str(a)[-10:]
'1679581912'

Mathematical relations

The arithmetic-geometric mean may be used to evaluate the following two parametric definite integrals:

\[ \begin{align}\begin{aligned}I_1 = \int_0^{\infty} \frac{1}{\sqrt{(x^2+a^2)(x^2+b^2)}} \,dx\\I_2 = \int_0^{\pi/2} \frac{1}{\sqrt{a^2 \cos^2(x) + b^2 \sin^2(x)}} \,dx\end{aligned}\end{align} \]

We have:

>>> mp.dps = 15
>>> a = 3
>>> b = 4
>>> f1 = lambda x: ((x**2+a**2)*(x**2+b**2))**-0.5
>>> f2 = lambda x: ((a*cos(x))**2 + (b*sin(x))**2)**-0.5
>>> quad(f1, [0, inf])
0.451115405388492
>>> quad(f2, [0, pi/2])
0.451115405388492
>>> pi/(2*agm(a,b))
0.451115405388492

A formula for \(\Gamma(1/4)\):

>>> gamma(0.25)
3.62560990822191
>>> sqrt(2*sqrt(2*pi**3)/agm(1,sqrt(2)))
3.62560990822191

Possible issues

The branch cut chosen for complex \(a\) and \(b\) is somewhat arbitrary.