using for loop to install conda package

.CondaPkg/env/Lib/site-packages/scipy/special/_precompute/cosine_cdf.py

@@ -0,0 +1,18 @@
+
+import mpmath
+
+
+def f(x):
+    return (mpmath.pi + x + mpmath.sin(x)) / (2*mpmath.pi)
+
+
+# Note: 40 digits might be overkill; a few more digits than the default
+# might be sufficient.
+mpmath.mp.dps = 40
+ts = mpmath.taylor(f, -mpmath.pi, 20)
+p, q = mpmath.pade(ts, 9, 10)
+
+p = [float(c) for c in p]
+q = [float(c) for c in q]
+print('p =', p)
+print('q =', q)

.CondaPkg/env/Lib/site-packages/scipy/special/_precompute/expn_asy.py

@@ -0,0 +1,54 @@
+"""Precompute the polynomials for the asymptotic expansion of the
+generalized exponential integral.
+
+Sources
+-------
+[1] NIST, Digital Library of Mathematical Functions,
+    https://dlmf.nist.gov/8.20#ii
+
+"""
+import os
+
+try:
+    import sympy
+    from sympy import Poly
+    x = sympy.symbols('x')
+except ImportError:
+    pass
+
+
+def generate_A(K):
+    A = [Poly(1, x)]
+    for k in range(K):
+        A.append(Poly(1 - 2*k*x, x)*A[k] + Poly(x*(x + 1))*A[k].diff())
+    return A
+
+
+WARNING = """\
+/* This file was automatically generated by _precompute/expn_asy.py.
+ * Do not edit it manually!
+ */
+"""
+
+
+def main():
+    print(__doc__)
+    fn = os.path.join('..', 'cephes', 'expn.h')
+
+    K = 12
+    A = generate_A(K)
+    with open(fn + '.new', 'w') as f:
+        f.write(WARNING)
+        f.write("#define nA {}\n".format(len(A)))
+        for k, Ak in enumerate(A):
+            tmp = ', '.join([str(x.evalf(18)) for x in Ak.coeffs()])
+            f.write("static const double A{}[] = {{{}}};\n".format(k, tmp))
+        tmp = ", ".join(["A{}".format(k) for k in range(K + 1)])
+        f.write("static const double *A[] = {{{}}};\n".format(tmp))
+        tmp = ", ".join([str(Ak.degree()) for Ak in A])
+        f.write("static const int Adegs[] = {{{}}};\n".format(tmp))
+    os.rename(fn + '.new', fn)
+
+
+if __name__ == "__main__":
+    main()

.CondaPkg/env/Lib/site-packages/scipy/special/_precompute/gammainc_asy.py

@@ -0,0 +1,116 @@
+"""
+Precompute coefficients of Temme's asymptotic expansion for gammainc.
+
+This takes about 8 hours to run on a 2.3 GHz Macbook Pro with 4GB ram.
+
+Sources:
+[1] NIST, "Digital Library of Mathematical Functions",
+    https://dlmf.nist.gov/
+
+"""
+import os
+from scipy.special._precompute.utils import lagrange_inversion
+
+try:
+    import mpmath as mp
+except ImportError:
+    pass
+
+
+def compute_a(n):
+    """a_k from DLMF 5.11.6"""
+    a = [mp.sqrt(2)/2]
+    for k in range(1, n):
+        ak = a[-1]/k
+        for j in range(1, len(a)):
+            ak -= a[j]*a[-j]/(j + 1)
+        ak /= a[0]*(1 + mp.mpf(1)/(k + 1))
+        a.append(ak)
+    return a
+
+
+def compute_g(n):
+    """g_k from DLMF 5.11.3/5.11.5"""
+    a = compute_a(2*n)
+    g = [mp.sqrt(2)*mp.rf(0.5, k)*a[2*k] for k in range(n)]
+    return g
+
+
+def eta(lam):
+    """Function from DLMF 8.12.1 shifted to be centered at 0."""
+    if lam > 0:
+        return mp.sqrt(2*(lam - mp.log(lam + 1)))
+    elif lam < 0:
+        return -mp.sqrt(2*(lam - mp.log(lam + 1)))
+    else:
+        return 0
+
+
+def compute_alpha(n):
+    """alpha_n from DLMF 8.12.13"""
+    coeffs = mp.taylor(eta, 0, n - 1)
+    return lagrange_inversion(coeffs)
+
+
+def compute_d(K, N):
+    """d_{k, n} from DLMF 8.12.12"""
+    M = N + 2*K
+    d0 = [-mp.mpf(1)/3]
+    alpha = compute_alpha(M + 2)
+    for n in range(1, M):
+        d0.append((n + 2)*alpha[n+2])
+    d = [d0]
+    g = compute_g(K)
+    for k in range(1, K):
+        dk = []
+        for n in range(M - 2*k):
+            dk.append((-1)**k*g[k]*d[0][n] + (n + 2)*d[k-1][n+2])
+        d.append(dk)
+    for k in range(K):
+        d[k] = d[k][:N]
+    return d
+
+
+header = \
+r"""/* This file was automatically generated by _precomp/gammainc.py.
+ * Do not edit it manually!
+ */
+
+#ifndef IGAM_H
+#define IGAM_H
+
+#define K {}
+#define N {}
+
+static const double d[K][N] =
+{{"""
+
+footer = \
+r"""
+#endif
+"""
+
+
+def main():
+    print(__doc__)
+    K = 25
+    N = 25
+    with mp.workdps(50):
+        d = compute_d(K, N)
+    fn = os.path.join(os.path.dirname(__file__), '..', 'cephes', 'igam.h')
+    with open(fn + '.new', 'w') as f:
+        f.write(header.format(K, N))
+        for k, row in enumerate(d):
+            row = [mp.nstr(x, 17, min_fixed=0, max_fixed=0) for x in row]
+            f.write('{')
+            f.write(", ".join(row))
+            if k < K - 1:
+                f.write('},\n')
+            else:
+                f.write('}};\n')
+        f.write(footer)
+    os.rename(fn + '.new', fn)
+
+
+if __name__ == "__main__":
+    main()

.CondaPkg/env/Lib/site-packages/scipy/special/_precompute/gammainc_data.py

@@ -0,0 +1,124 @@
+"""Compute gammainc and gammaincc for large arguments and parameters
+and save the values to data files for use in tests. We can't just
+compare to mpmath's gammainc in test_mpmath.TestSystematic because it
+would take too long.
+
+Note that mpmath's gammainc is computed using hypercomb, but since it
+doesn't allow the user to increase the maximum number of terms used in
+the series it doesn't converge for many arguments. To get around this
+we copy the mpmath implementation but use more terms.
+
+This takes about 17 minutes to run on a 2.3 GHz Macbook Pro with 4GB
+ram.
+
+Sources:
+[1] Fredrik Johansson and others. mpmath: a Python library for
+    arbitrary-precision floating-point arithmetic (version 0.19),
+    December 2013. http://mpmath.org/.
+
+"""
+import os
+from time import time
+import numpy as np
+from numpy import pi
+
+from scipy.special._mptestutils import mpf2float
+
+try:
+    import mpmath as mp
+except ImportError:
+    pass
+
+
+def gammainc(a, x, dps=50, maxterms=10**8):
+    """Compute gammainc exactly like mpmath does but allow for more
+    summands in hypercomb. See
+
+    mpmath/functions/expintegrals.py#L134
+
+    in the mpmath github repository.
+
+    """
+    with mp.workdps(dps):
+        z, a, b = mp.mpf(a), mp.mpf(x), mp.mpf(x)
+        G = [z]
+        negb = mp.fneg(b, exact=True)
+
+        def h(z):
+            T1 = [mp.exp(negb), b, z], [1, z, -1], [], G, [1], [1+z], b
+            return (T1,)
+
+        res = mp.hypercomb(h, [z], maxterms=maxterms)
+        return mpf2float(res)
+
+
+def gammaincc(a, x, dps=50, maxterms=10**8):
+    """Compute gammaincc exactly like mpmath does but allow for more
+    terms in hypercomb. See
+
+    mpmath/functions/expintegrals.py#L187
+
+    in the mpmath github repository.
+
+    """
+    with mp.workdps(dps):
+        z, a = a, x
+
+        if mp.isint(z):
+            try:
+                # mpmath has a fast integer path
+                return mpf2float(mp.gammainc(z, a=a, regularized=True))
+            except mp.libmp.NoConvergence:
+                pass
+        nega = mp.fneg(a, exact=True)
+        G = [z]
+        # Use 2F0 series when possible; fall back to lower gamma representation
+        try:
+            def h(z):
+                r = z-1
+                return [([mp.exp(nega), a], [1, r], [], G, [1, -r], [], 1/nega)]
+            return mpf2float(mp.hypercomb(h, [z], force_series=True))
+        except mp.libmp.NoConvergence:
+            def h(z):
+                T1 = [], [1, z-1], [z], G, [], [], 0
+                T2 = [-mp.exp(nega), a, z], [1, z, -1], [], G, [1], [1+z], a
+                return T1, T2
+            return mpf2float(mp.hypercomb(h, [z], maxterms=maxterms))
+
+
+def main():
+    t0 = time()
+    # It would be nice to have data for larger values, but either this
+    # requires prohibitively large precision (dps > 800) or mpmath has
+    # a bug. For example, gammainc(1e20, 1e20, dps=800) returns a
+    # value around 0.03, while the true value should be close to 0.5
+    # (DLMF 8.12.15).
+    print(__doc__)
+    pwd = os.path.dirname(__file__)
+    r = np.logspace(4, 14, 30)
+    ltheta = np.logspace(np.log10(pi/4), np.log10(np.arctan(0.6)), 30)
+    utheta = np.logspace(np.log10(pi/4), np.log10(np.arctan(1.4)), 30)
+
+    regimes = [(gammainc, ltheta), (gammaincc, utheta)]
+    for func, theta in regimes:
+        rg, thetag = np.meshgrid(r, theta)
+        a, x = rg*np.cos(thetag), rg*np.sin(thetag)
+        a, x = a.flatten(), x.flatten()
+        dataset = []
+        for i, (a0, x0) in enumerate(zip(a, x)):
+            if func == gammaincc:
+                # Exploit the fast integer path in gammaincc whenever
+                # possible so that the computation doesn't take too
+                # long
+                a0, x0 = np.floor(a0), np.floor(x0)
+            dataset.append((a0, x0, func(a0, x0)))
+        dataset = np.array(dataset)
+        filename = os.path.join(pwd, '..', 'tests', 'data', 'local',
+                                '{}.txt'.format(func.__name__))
+        np.savetxt(filename, dataset)
+
+    print("{} minutes elapsed".format((time() - t0)/60))
+
+
+if __name__ == "__main__":
+    main()

.CondaPkg/env/Lib/site-packages/scipy/special/_precompute/lambertw.py

@@ -0,0 +1,68 @@
+"""Compute a Pade approximation for the principal branch of the
+Lambert W function around 0 and compare it to various other
+approximations.
+
+"""
+import numpy as np
+
+try:
+    import mpmath
+    import matplotlib.pyplot as plt
+except ImportError:
+    pass
+
+
+def lambertw_pade():
+    derivs = [mpmath.diff(mpmath.lambertw, 0, n=n) for n in range(6)]
+    p, q = mpmath.pade(derivs, 3, 2)
+    return p, q
+
+
+def main():
+    print(__doc__)
+    with mpmath.workdps(50):
+        p, q = lambertw_pade()
+        p, q = p[::-1], q[::-1]
+        print("p = {}".format(p))
+        print("q = {}".format(q))
+
+    x, y = np.linspace(-1.5, 1.5, 75), np.linspace(-1.5, 1.5, 75)
+    x, y = np.meshgrid(x, y)
+    z = x + 1j*y
+    lambertw_std = []
+    for z0 in z.flatten():
+        lambertw_std.append(complex(mpmath.lambertw(z0)))
+    lambertw_std = np.array(lambertw_std).reshape(x.shape)
+
+    fig, axes = plt.subplots(nrows=3, ncols=1)
+    # Compare Pade approximation to true result
+    p = np.array([float(p0) for p0 in p])
+    q = np.array([float(q0) for q0 in q])
+    pade_approx = np.polyval(p, z)/np.polyval(q, z)
+    pade_err = abs(pade_approx - lambertw_std)
+    axes[0].pcolormesh(x, y, pade_err)
+    # Compare two terms of asymptotic series to true result
+    asy_approx = np.log(z) - np.log(np.log(z))
+    asy_err = abs(asy_approx - lambertw_std)
+    axes[1].pcolormesh(x, y, asy_err)
+    # Compare two terms of the series around the branch point to the
+    # true result
+    p = np.sqrt(2*(np.exp(1)*z + 1))
+    series_approx = -1 + p - p**2/3
+    series_err = abs(series_approx - lambertw_std)
+    im = axes[2].pcolormesh(x, y, series_err)
+
+    fig.colorbar(im, ax=axes.ravel().tolist())
+    plt.show()
+
+    fig, ax = plt.subplots(nrows=1, ncols=1)
+    pade_better = pade_err < asy_err
+    im = ax.pcolormesh(x, y, pade_better)
+    t = np.linspace(-0.3, 0.3)
+    ax.plot(-2.5*abs(t) - 0.2, t, 'r')
+    fig.colorbar(im, ax=ax)
+    plt.show()
+
+
+if __name__ == '__main__':
+    main()

.CondaPkg/env/Lib/site-packages/scipy/special/_precompute/loggamma.py

@@ -0,0 +1,43 @@
+"""Precompute series coefficients for log-Gamma."""
+
+try:
+    import mpmath
+except ImportError:
+    pass
+
+
+def stirling_series(N):
+    with mpmath.workdps(100):
+        coeffs = [mpmath.bernoulli(2*n)/(2*n*(2*n - 1))
+                  for n in range(1, N + 1)]
+    return coeffs
+
+
+def taylor_series_at_1(N):
+    coeffs = []
+    with mpmath.workdps(100):
+        coeffs.append(-mpmath.euler)
+        for n in range(2, N + 1):
+            coeffs.append((-1)**n*mpmath.zeta(n)/n)
+    return coeffs
+
+
+def main():
+    print(__doc__)
+    print()
+    stirling_coeffs = [mpmath.nstr(x, 20, min_fixed=0, max_fixed=0)
+                       for x in stirling_series(8)[::-1]]
+    taylor_coeffs = [mpmath.nstr(x, 20, min_fixed=0, max_fixed=0)
+                     for x in taylor_series_at_1(23)[::-1]]
+    print("Stirling series coefficients")
+    print("----------------------------")
+    print("\n".join(stirling_coeffs))
+    print()
+    print("Taylor series coefficients")
+    print("--------------------------")
+    print("\n".join(taylor_coeffs))
+    print()
+
+
+if __name__ == '__main__':
+    main()

.CondaPkg/env/Lib/site-packages/scipy/special/_precompute/struve_convergence.py

@@ -0,0 +1,120 @@
+"""
+Convergence regions of the expansions used in ``struve.c``
+
+Note that for v >> z both functions tend rapidly to 0,
+and for v << -z, they tend to infinity.
+
+The floating-point functions over/underflow in the lower left and right
+corners of the figure.
+
+
+Figure legend
+=============
+
+Red region
+    Power series is close (1e-12) to the mpmath result
+
+Blue region
+    Asymptotic series is close to the mpmath result
+
+Green region
+    Bessel series is close to the mpmath result
+
+Dotted colored lines
+    Boundaries of the regions
+
+Solid colored lines
+    Boundaries estimated by the routine itself. These will be used
+    for determining which of the results to use.
+
+Black dashed line
+    The line z = 0.7*|v| + 12
+
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+
+import mpmath
+
+
+def err_metric(a, b, atol=1e-290):
+    m = abs(a - b) / (atol + abs(b))
+    m[np.isinf(b) & (a == b)] = 0
+    return m
+
+
+def do_plot(is_h=True):
+    from scipy.special._ufuncs import (_struve_power_series,
+                                       _struve_asymp_large_z,
+                                       _struve_bessel_series)
+
+    vs = np.linspace(-1000, 1000, 91)
+    zs = np.sort(np.r_[1e-5, 1.0, np.linspace(0, 700, 91)[1:]])
+
+    rp = _struve_power_series(vs[:,None], zs[None,:], is_h)
+    ra = _struve_asymp_large_z(vs[:,None], zs[None,:], is_h)
+    rb = _struve_bessel_series(vs[:,None], zs[None,:], is_h)
+
+    mpmath.mp.dps = 50
+    if is_h:
+        sh = lambda v, z: float(mpmath.struveh(mpmath.mpf(v), mpmath.mpf(z)))
+    else:
+        sh = lambda v, z: float(mpmath.struvel(mpmath.mpf(v), mpmath.mpf(z)))
+    ex = np.vectorize(sh, otypes='d')(vs[:,None], zs[None,:])
+
+    err_a = err_metric(ra[0], ex) + 1e-300
+    err_p = err_metric(rp[0], ex) + 1e-300
+    err_b = err_metric(rb[0], ex) + 1e-300
+
+    err_est_a = abs(ra[1]/ra[0])
+    err_est_p = abs(rp[1]/rp[0])
+    err_est_b = abs(rb[1]/rb[0])
+
+    z_cutoff = 0.7*abs(vs) + 12
+
+    levels = [-1000, -12]
+
+    plt.cla()
+
+    plt.hold(1)
+    plt.contourf(vs, zs, np.log10(err_p).T, levels=levels, colors=['r', 'r'], alpha=0.1)
+    plt.contourf(vs, zs, np.log10(err_a).T, levels=levels, colors=['b', 'b'], alpha=0.1)
+    plt.contourf(vs, zs, np.log10(err_b).T, levels=levels, colors=['g', 'g'], alpha=0.1)
+
+    plt.contour(vs, zs, np.log10(err_p).T, levels=levels, colors=['r', 'r'], linestyles=[':', ':'])
+    plt.contour(vs, zs, np.log10(err_a).T, levels=levels, colors=['b', 'b'], linestyles=[':', ':'])
+    plt.contour(vs, zs, np.log10(err_b).T, levels=levels, colors=['g', 'g'], linestyles=[':', ':'])
+
+    lp = plt.contour(vs, zs, np.log10(err_est_p).T, levels=levels, colors=['r', 'r'], linestyles=['-', '-'])
+    la = plt.contour(vs, zs, np.log10(err_est_a).T, levels=levels, colors=['b', 'b'], linestyles=['-', '-'])
+    lb = plt.contour(vs, zs, np.log10(err_est_b).T, levels=levels, colors=['g', 'g'], linestyles=['-', '-'])
+
+    plt.clabel(lp, fmt={-1000: 'P', -12: 'P'})
+    plt.clabel(la, fmt={-1000: 'A', -12: 'A'})
+    plt.clabel(lb, fmt={-1000: 'B', -12: 'B'})
+
+    plt.plot(vs, z_cutoff, 'k--')
+
+    plt.xlim(vs.min(), vs.max())
+    plt.ylim(zs.min(), zs.max())
+
+    plt.xlabel('v')
+    plt.ylabel('z')
+
+
+def main():
+    plt.clf()
+    plt.subplot(121)
+    do_plot(True)
+    plt.title('Struve H')
+
+    plt.subplot(122)
+    do_plot(False)
+    plt.title('Struve L')
+
+    plt.savefig('struve_convergence.png')
+    plt.show()
+
+
+if __name__ == "__main__":
+    main()

.CondaPkg/env/Lib/site-packages/scipy/special/_precompute/utils.py

@@ -0,0 +1,38 @@
+try:
+    import mpmath as mp
+except ImportError:
+    pass
+
+try:
+    from sympy.abc import x
+except ImportError:
+    pass
+
+
+def lagrange_inversion(a):
+    """Given a series
+
+    f(x) = a[1]*x + a[2]*x**2 + ... + a[n-1]*x**(n - 1),
+
+    use the Lagrange inversion formula to compute a series
+
+    g(x) = b[1]*x + b[2]*x**2 + ... + b[n-1]*x**(n - 1)
+
+    so that f(g(x)) = g(f(x)) = x mod x**n. We must have a[0] = 0, so
+    necessarily b[0] = 0 too.
+
+    The algorithm is naive and could be improved, but speed isn't an
+    issue here and it's easy to read.
+
+    """
+    n = len(a)
+    f = sum(a[i]*x**i for i in range(n))
+    h = (x/f).series(x, 0, n).removeO()
+    hpower = [h**0]
+    for k in range(n):
+        hpower.append((hpower[-1]*h).expand())
+    b = [mp.mpf(0)]
+    for k in range(1, n):
+        b.append(hpower[k].coeff(x, k - 1)/k)
+    b = [mp.mpf(x) for x in b]
+    return b

.CondaPkg/env/Lib/site-packages/scipy/special/_precompute/wright_bessel.py

@@ -0,0 +1,342 @@
+"""Precompute coefficients of several series expansions
+of Wright's generalized Bessel function Phi(a, b, x).
+
+See https://dlmf.nist.gov/10.46.E1 with rho=a, beta=b, z=x.
+"""
+from argparse import ArgumentParser, RawTextHelpFormatter
+import numpy as np
+from scipy.integrate import quad
+from scipy.optimize import minimize_scalar, curve_fit
+from time import time
+
+try:
+    import sympy
+    from sympy import EulerGamma, Rational, S, Sum, \
+        factorial, gamma, gammasimp, pi, polygamma, symbols, zeta
+    from sympy.polys.polyfuncs import horner
+except ImportError:
+    pass
+
+
+def series_small_a():
+    """Tylor series expansion of Phi(a, b, x) in a=0 up to order 5.
+    """
+    order = 5
+    a, b, x, k = symbols("a b x k")
+    A = []  # terms with a
+    X = []  # terms with x
+    B = []  # terms with b (polygammas)
+    # Phi(a, b, x) = exp(x)/gamma(b) * sum(A[i] * X[i] * B[i])
+    expression = Sum(x**k/factorial(k)/gamma(a*k+b), (k, 0, S.Infinity))
+    expression = gamma(b)/sympy.exp(x) * expression
+
+    # nth term of taylor series in a=0: a^n/n! * (d^n Phi(a, b, x)/da^n at a=0)
+    for n in range(0, order+1):
+        term = expression.diff(a, n).subs(a, 0).simplify().doit()
+        # set the whole bracket involving polygammas to 1
+        x_part = (term.subs(polygamma(0, b), 1)
+                  .replace(polygamma, lambda *args: 0))
+        # sign convetion: x part always positive
+        x_part *= (-1)**n
+
+        A.append(a**n/factorial(n))
+        X.append(horner(x_part))
+        B.append(horner((term/x_part).simplify()))
+
+    s = "Tylor series expansion of Phi(a, b, x) in a=0 up to order 5.\n"
+    s += "Phi(a, b, x) = exp(x)/gamma(b) * sum(A[i] * X[i] * B[i], i=0..5)\n"
+    for name, c in zip(['A', 'X', 'B'], [A, X, B]):
+        for i in range(len(c)):
+            s += f"\n{name}[{i}] = " + str(c[i])
+    return s
+
+
+# expansion of digamma
+def dg_series(z, n):
+    """Symbolic expansion of digamma(z) in z=0 to order n.
+
+    See https://dlmf.nist.gov/5.7.E4 and with https://dlmf.nist.gov/5.5.E2
+    """
+    k = symbols("k")
+    return -1/z - EulerGamma + \
+        sympy.summation((-1)**k * zeta(k) * z**(k-1), (k, 2, n+1))
+
+
+def pg_series(k, z, n):
+    """Symbolic expansion of polygamma(k, z) in z=0 to order n."""
+    return sympy.diff(dg_series(z, n+k), z, k)
+
+
+def series_small_a_small_b():
+    """Tylor series expansion of Phi(a, b, x) in a=0 and b=0 up to order 5.
+
+    Be aware of cancellation of poles in b=0 of digamma(b)/Gamma(b) and
+    polygamma functions.
+
+    digamma(b)/Gamma(b) = -1 - 2*M_EG*b + O(b^2)
+    digamma(b)^2/Gamma(b) = 1/b + 3*M_EG + b*(-5/12*PI^2+7/2*M_EG^2) + O(b^2)
+    polygamma(1, b)/Gamma(b) = 1/b + M_EG + b*(1/12*PI^2 + 1/2*M_EG^2) + O(b^2)
+    and so on.
+    """
+    order = 5
+    a, b, x, k = symbols("a b x k")
+    M_PI, M_EG, M_Z3 = symbols("M_PI M_EG M_Z3")
+    c_subs = {pi: M_PI, EulerGamma: M_EG, zeta(3): M_Z3}
+    A = []  # terms with a
+    X = []  # terms with x
+    B = []  # terms with b (polygammas expanded)
+    C = []  # terms that generate B
+    # Phi(a, b, x) = exp(x) * sum(A[i] * X[i] * B[i])
+    # B[0] = 1
+    # B[k] = sum(C[k] * b**k/k!, k=0..)
+    # Note: C[k] can be obtained from a series expansion of 1/gamma(b).
+    expression = gamma(b)/sympy.exp(x) * \
+        Sum(x**k/factorial(k)/gamma(a*k+b), (k, 0, S.Infinity))
+
+    # nth term of taylor series in a=0: a^n/n! * (d^n Phi(a, b, x)/da^n at a=0)
+    for n in range(0, order+1):
+        term = expression.diff(a, n).subs(a, 0).simplify().doit()
+        # set the whole bracket involving polygammas to 1
+        x_part = (term.subs(polygamma(0, b), 1)
+                  .replace(polygamma, lambda *args: 0))
+        # sign convetion: x part always positive
+        x_part *= (-1)**n
+        # expansion of polygamma part with 1/gamma(b)
+        pg_part = term/x_part/gamma(b)
+        if n >= 1:
+            # Note: highest term is digamma^n
+            pg_part = pg_part.replace(polygamma,
+                                      lambda k, x: pg_series(k, x, order+1+n))
+            pg_part = (pg_part.series(b, 0, n=order+1-n)
+                       .removeO()
+                       .subs(polygamma(2, 1), -2*zeta(3))
+                       .simplify()
+                       )
+
+        A.append(a**n/factorial(n))
+        X.append(horner(x_part))
+        B.append(pg_part)
+
+    # Calculate C and put in the k!
+    C = sympy.Poly(B[1].subs(c_subs), b).coeffs()
+    C.reverse()
+    for i in range(len(C)):
+        C[i] = (C[i] * factorial(i)).simplify()
+
+    s = "Tylor series expansion of Phi(a, b, x) in a=0 and b=0 up to order 5."
+    s += "\nPhi(a, b, x) = exp(x) * sum(A[i] * X[i] * B[i], i=0..5)\n"
+    s += "B[0] = 1\n"
+    s += "B[i] = sum(C[k+i-1] * b**k/k!, k=0..)\n"
+    s += "\nM_PI = pi"
+    s += "\nM_EG = EulerGamma"
+    s += "\nM_Z3 = zeta(3)"
+    for name, c in zip(['A', 'X'], [A, X]):
+        for i in range(len(c)):
+            s += f"\n{name}[{i}] = "
+            s += str(c[i])
+    # For C, do also compute the values numerically
+    for i in range(len(C)):
+        s += f"\n# C[{i}] = "
+        s += str(C[i])
+        s += f"\nC[{i}] = "
+        s += str(C[i].subs({M_EG: EulerGamma, M_PI: pi, M_Z3: zeta(3)})
+                 .evalf(17))
+
+    # Does B have the assumed structure?
+    s += "\n\nTest if B[i] does have the assumed structure."
+    s += "\nC[i] are derived from B[1] allone."
+    s += "\nTest B[2] == C[1] + b*C[2] + b^2/2*C[3] + b^3/6*C[4] + .."
+    test = sum([b**k/factorial(k) * C[k+1] for k in range(order-1)])
+    test = (test - B[2].subs(c_subs)).simplify()
+    s += f"\ntest successful = {test==S(0)}"
+    s += "\nTest B[3] == C[2] + b*C[3] + b^2/2*C[4] + .."
+    test = sum([b**k/factorial(k) * C[k+2] for k in range(order-2)])
+    test = (test - B[3].subs(c_subs)).simplify()
+    s += f"\ntest successful = {test==S(0)}"
+    return s
+
+
+def asymptotic_series():
+    """Asymptotic expansion for large x.
+
+    Phi(a, b, x) ~ Z^(1/2-b) * exp((1+a)/a * Z) * sum_k (-1)^k * C_k / Z^k
+    Z = (a*x)^(1/(1+a))
+
+    Wright (1935) lists the coefficients C_0 and C_1 (he calls them a_0 and
+    a_1). With slightly different notation, Paris (2017) lists coefficients
+    c_k up to order k=3.
+    Paris (2017) uses ZP = (1+a)/a * Z  (ZP = Z of Paris) and
+    C_k = C_0 * (-a/(1+a))^k * c_k
+    """
+    order = 8
+
+    class g(sympy.Function):
+        """Helper function g according to Wright (1935)
+
+        g(n, rho, v) = (1 + (rho+2)/3 * v + (rho+2)*(rho+3)/(2*3) * v^2 + ...)
+
+        Note: Wright (1935) uses square root of above definition.
+        """
+        nargs = 3
+
+        @classmethod
+        def eval(cls, n, rho, v):
+            if not n >= 0:
+                raise ValueError("must have n >= 0")
+            elif n == 0:
+                return 1
+            else:
+                return g(n-1, rho, v) \
+                    + gammasimp(gamma(rho+2+n)/gamma(rho+2)) \
+                    / gammasimp(gamma(3+n)/gamma(3))*v**n
+
+    class coef_C(sympy.Function):
+        """Calculate coefficients C_m for integer m.
+
+        C_m is the coefficient of v^(2*m) in the Taylor expansion in v=0 of
+        Gamma(m+1/2)/(2*pi) * (2/(rho+1))^(m+1/2) * (1-v)^(-b)
+            * g(rho, v)^(-m-1/2)
+        """
+        nargs = 3
+
+        @classmethod
+        def eval(cls, m, rho, beta):
+            if not m >= 0:
+                raise ValueError("must have m >= 0")
+
+            v = symbols("v")
+            expression = (1-v)**(-beta) * g(2*m, rho, v)**(-m-Rational(1, 2))
+            res = expression.diff(v, 2*m).subs(v, 0) / factorial(2*m)
+            res = res * (gamma(m + Rational(1, 2)) / (2*pi)
+                         * (2/(rho+1))**(m + Rational(1, 2)))
+            return res
+
+    # in order to have nice ordering/sorting of expressions, we set a = xa.
+    xa, b, xap1 = symbols("xa b xap1")
+    C0 = coef_C(0, xa, b)
+    # a1 = a(1, rho, beta)
+    s = "Asymptotic expansion for large x\n"
+    s += "Phi(a, b, x) = Z**(1/2-b) * exp((1+a)/a * Z) \n"
+    s += "               * sum((-1)**k * C[k]/Z**k, k=0..6)\n\n"
+    s += "Z      = pow(a * x, 1/(1+a))\n"
+    s += "A[k]   = pow(a, k)\n"
+    s += "B[k]   = pow(b, k)\n"
+    s += "Ap1[k] = pow(1+a, k)\n\n"
+    s += "C[0] = 1./sqrt(2. * M_PI * Ap1[1])\n"
+    for i in range(1, order+1):
+        expr = (coef_C(i, xa, b) / (C0/(1+xa)**i)).simplify()
+        factor = [x.denominator() for x in sympy.Poly(expr).coeffs()]
+        factor = sympy.lcm(factor)
+        expr = (expr * factor).simplify().collect(b, sympy.factor)
+        expr = expr.xreplace({xa+1: xap1})
+        s += f"C[{i}] = C[0] / ({factor} * Ap1[{i}])\n"
+        s += f"C[{i}] *= {str(expr)}\n\n"
+    import re
+    re_a = re.compile(r'xa\*\*(\d+)')
+    s = re_a.sub(r'A[\1]', s)
+    re_b = re.compile(r'b\*\*(\d+)')
+    s = re_b.sub(r'B[\1]', s)
+    s = s.replace('xap1', 'Ap1[1]')
+    s = s.replace('xa', 'a')
+    # max integer = 2^31-1 = 2,147,483,647. Solution: Put a point after 10
+    # or more digits.
+    re_digits = re.compile(r'(\d{10,})')
+    s = re_digits.sub(r'\1.', s)
+    return s
+
+
+def optimal_epsilon_integral():
+    """Fit optimal choice of epsilon for integral representation.
+
+    The integrand of
+        int_0^pi P(eps, a, b, x, phi) * dphi
+    can exhibit oscillatory behaviour. It stems from the cosine of P and can be
+    minimized by minimizing the arc length of the argument
+        f(phi) = eps * sin(phi) - x * eps^(-a) * sin(a * phi) + (1 - b) * phi
+    of cos(f(phi)).
+    We minimize the arc length in eps for a grid of values (a, b, x) and fit a
+    parametric function to it.
+    """
+    def fp(eps, a, b, x, phi):
+        """Derivative of f w.r.t. phi."""
+        eps_a = np.power(1. * eps, -a)
+        return eps * np.cos(phi) - a * x * eps_a * np.cos(a * phi) + 1 - b
+
+    def arclength(eps, a, b, x, epsrel=1e-2, limit=100):
+        """Compute Arc length of f.
+
+        Note that the arg length of a function f fro t0 to t1 is given by
+            int_t0^t1 sqrt(1 + f'(t)^2) dt
+        """
+        return quad(lambda phi: np.sqrt(1 + fp(eps, a, b, x, phi)**2),
+                    0, np.pi,
+                    epsrel=epsrel, limit=100)[0]
+
+    # grid of minimal arc length values
+    data_a = [1e-3, 0.1, 0.5, 0.9, 1, 2, 4, 5, 6, 8]
+    data_b = [0, 1, 4, 7, 10]
+    data_x = [1, 1.5, 2, 4, 10, 20, 50, 100, 200, 500, 1e3, 5e3, 1e4]
+    data_a, data_b, data_x = np.meshgrid(data_a, data_b, data_x)
+    data_a, data_b, data_x = (data_a.flatten(), data_b.flatten(),
+                              data_x.flatten())
+    best_eps = []
+    for i in range(data_x.size):
+        best_eps.append(
+            minimize_scalar(lambda eps: arclength(eps, data_a[i], data_b[i],
+                                                  data_x[i]),
+                            bounds=(1e-3, 1000),
+                            method='Bounded', options={'xatol': 1e-3}).x
+        )
+    best_eps = np.array(best_eps)
+    # pandas would be nice, but here a dictionary is enough
+    df = {'a': data_a,
+          'b': data_b,
+          'x': data_x,
+          'eps': best_eps,
+          }
+
+    def func(data, A0, A1, A2, A3, A4, A5):
+        """Compute parametric function to fit."""
+        a = data['a']
+        b = data['b']
+        x = data['x']
+        return (A0 * b * np.exp(-0.5 * a)
+                + np.exp(A1 + 1 / (1 + a) * np.log(x) - A2 * np.exp(-A3 * a)
+                         + A4 / (1 + np.exp(A5 * a))))
+
+    func_params = list(curve_fit(func, df, df['eps'], method='trf')[0])
+
+    s = "Fit optimal eps for integrand P via minimal arc length\n"
+    s += "with parametric function:\n"
+    s += "optimal_eps = (A0 * b * exp(-a/2) + exp(A1 + 1 / (1 + a) * log(x)\n"
+    s += "              - A2 * exp(-A3 * a) + A4 / (1 + exp(A5 * a)))\n\n"
+    s += "Fitted parameters A0 to A5 are:\n"
+    s += ', '.join(['{:.5g}'.format(x) for x in func_params])
+    return s
+
+
+def main():
+    t0 = time()
+    parser = ArgumentParser(description=__doc__,
+                            formatter_class=RawTextHelpFormatter)
+    parser.add_argument('action', type=int, choices=[1, 2, 3, 4],
+                        help='chose what expansion to precompute\n'
+                             '1 : Series for small a\n'
+                             '2 : Series for small a and small b\n'
+                             '3 : Asymptotic series for large x\n'
+                             '    This may take some time (>4h).\n'
+                             '4 : Fit optimal eps for integral representation.'
+                        )
+    args = parser.parse_args()
+
+    switch = {1: lambda: print(series_small_a()),
+              2: lambda: print(series_small_a_small_b()),
+              3: lambda: print(asymptotic_series()),
+              4: lambda: print(optimal_epsilon_integral())
+              }
+    switch.get(args.action, lambda: print("Invalid input."))()
+    print("\n{:.1f} minutes elapsed.\n".format((time() - t0)/60))
+
+
+if __name__ == '__main__':
+    main()

.CondaPkg/env/Lib/site-packages/scipy/special/_precompute/wright_bessel_data.py

@@ -0,0 +1,152 @@
+"""Compute a grid of values for Wright's generalized Bessel function
+and save the values to data files for use in tests. Using mpmath directly in
+tests would take too long.
+
+This takes about 10 minutes to run on a 2.7 GHz i7 Macbook Pro.
+"""
+from functools import lru_cache
+import os
+from time import time
+
+import numpy as np
+from scipy.special._mptestutils import mpf2float
+
+try:
+    import mpmath as mp
+except ImportError:
+    pass
+
+# exp_inf: smallest value x for which exp(x) == inf
+exp_inf = 709.78271289338403
+
+
+# 64 Byte per value
+@lru_cache(maxsize=100_000)
+def rgamma_cached(x, dps):
+    with mp.workdps(dps):
+        return mp.rgamma(x)
+
+
+def mp_wright_bessel(a, b, x, dps=50, maxterms=2000):
+    """Compute Wright's generalized Bessel function as Series with mpmath.
+    """
+    with mp.workdps(dps):
+        a, b, x = mp.mpf(a), mp.mpf(b), mp.mpf(x)
+        res = mp.nsum(lambda k: x**k / mp.fac(k)
+                      * rgamma_cached(a * k + b, dps=dps),
+                      [0, mp.inf],
+                      tol=dps, method='s', steps=[maxterms]
+                      )
+        return mpf2float(res)
+
+
+def main():
+    t0 = time()
+    print(__doc__)
+    pwd = os.path.dirname(__file__)
+    eps = np.finfo(float).eps * 100
+
+    a_range = np.array([eps,
+                        1e-4 * (1 - eps), 1e-4, 1e-4 * (1 + eps),
+                        1e-3 * (1 - eps), 1e-3, 1e-3 * (1 + eps),
+                        0.1, 0.5,
+                        1 * (1 - eps), 1, 1 * (1 + eps),
+                        1.5, 2, 4.999, 5, 10])
+    b_range = np.array([0, eps, 1e-10, 1e-5, 0.1, 1, 2, 10, 20, 100])
+    x_range = np.array([0, eps, 1 - eps, 1, 1 + eps,
+                        1.5,
+                        2 - eps, 2, 2 + eps,
+                        9 - eps, 9, 9 + eps,
+                        10 * (1 - eps), 10, 10 * (1 + eps),
+                        100 * (1 - eps), 100, 100 * (1 + eps),
+                        500, exp_inf, 1e3, 1e5, 1e10, 1e20])
+
+    a_range, b_range, x_range = np.meshgrid(a_range, b_range, x_range,
+                                            indexing='ij')
+    a_range = a_range.flatten()
+    b_range = b_range.flatten()
+    x_range = x_range.flatten()
+
+    # filter out some values, especially too large x
+    bool_filter = ~((a_range < 5e-3) & (x_range >= exp_inf))
+    bool_filter = bool_filter & ~((a_range < 0.2) & (x_range > exp_inf))
+    bool_filter = bool_filter & ~((a_range < 0.5) & (x_range > 1e3))
+    bool_filter = bool_filter & ~((a_range < 0.56) & (x_range > 5e3))
+    bool_filter = bool_filter & ~((a_range < 1) & (x_range > 1e4))
+    bool_filter = bool_filter & ~((a_range < 1.4) & (x_range > 1e5))
+    bool_filter = bool_filter & ~((a_range < 1.8) & (x_range > 1e6))
+    bool_filter = bool_filter & ~((a_range < 2.2) & (x_range > 1e7))
+    bool_filter = bool_filter & ~((a_range < 2.5) & (x_range > 1e8))
+    bool_filter = bool_filter & ~((a_range < 2.9) & (x_range > 1e9))
+    bool_filter = bool_filter & ~((a_range < 3.3) & (x_range > 1e10))
+    bool_filter = bool_filter & ~((a_range < 3.7) & (x_range > 1e11))
+    bool_filter = bool_filter & ~((a_range < 4) & (x_range > 1e12))
+    bool_filter = bool_filter & ~((a_range < 4.4) & (x_range > 1e13))
+    bool_filter = bool_filter & ~((a_range < 4.7) & (x_range > 1e14))
+    bool_filter = bool_filter & ~((a_range < 5.1) & (x_range > 1e15))
+    bool_filter = bool_filter & ~((a_range < 5.4) & (x_range > 1e16))
+    bool_filter = bool_filter & ~((a_range < 5.8) & (x_range > 1e17))
+    bool_filter = bool_filter & ~((a_range < 6.2) & (x_range > 1e18))
+    bool_filter = bool_filter & ~((a_range < 6.2) & (x_range > 1e18))
+    bool_filter = bool_filter & ~((a_range < 6.5) & (x_range > 1e19))
+    bool_filter = bool_filter & ~((a_range < 6.9) & (x_range > 1e20))
+
+    # filter out known values that do not meet the required numerical accuracy
+    # see test test_wright_data_grid_failures
+    failing = np.array([
+        [0.1, 100, 709.7827128933841],
+        [0.5, 10, 709.7827128933841],
+        [0.5, 10, 1000],
+        [0.5, 100, 1000],
+        [1, 20, 100000],
+        [1, 100, 100000],
+        [1.0000000000000222, 20, 100000],
+        [1.0000000000000222, 100, 100000],
+        [1.5, 0, 500],
+        [1.5, 2.220446049250313e-14, 500],
+        [1.5, 1.e-10, 500],
+        [1.5, 1.e-05, 500],
+        [1.5, 0.1, 500],
+        [1.5, 20, 100000],
+        [1.5, 100, 100000],
+        ]).tolist()
+
+    does_fail = np.full_like(a_range, False, dtype=bool)
+    for i in range(x_range.size):
+        if [a_range[i], b_range[i], x_range[i]] in failing:
+            does_fail[i] = True
+
+    # filter and flatten
+    a_range = a_range[bool_filter]
+    b_range = b_range[bool_filter]
+    x_range = x_range[bool_filter]
+    does_fail = does_fail[bool_filter]
+
+    dataset = []
+    print(f"Computing {x_range.size} single points.")
+    print("Tests will fail for the following data points:")
+    for i in range(x_range.size):
+        a = a_range[i]
+        b = b_range[i]
+        x = x_range[i]
+        # take care of difficult corner cases
+        maxterms = 1000
+        if a < 1e-6 and x >= exp_inf/10:
+            maxterms = 2000
+        f = mp_wright_bessel(a, b, x, maxterms=maxterms)
+        if does_fail[i]:
+            print("failing data point a, b, x, value = "
+                  f"[{a}, {b}, {x}, {f}]")
+        else:
+            dataset.append((a, b, x, f))
+    dataset = np.array(dataset)
+
+    filename = os.path.join(pwd, '..', 'tests', 'data', 'local',
+                            'wright_bessel.txt')
+    np.savetxt(filename, dataset)
+
+    print("{:.1f} minutes elapsed".format((time() - t0)/60))
+
+
+if __name__ == "__main__":
+    main()

.CondaPkg/env/Lib/site-packages/scipy/special/_precompute/wrightomega.py

@@ -0,0 +1,41 @@
+import numpy as np
+
+try:
+    import mpmath
+except ImportError:
+    pass
+
+
+def mpmath_wrightomega(x):
+    return mpmath.lambertw(mpmath.exp(x), mpmath.mpf('-0.5'))
+
+
+def wrightomega_series_error(x):
+    series = x
+    desired = mpmath_wrightomega(x)
+    return abs(series - desired) / desired
+
+
+def wrightomega_exp_error(x):
+    exponential_approx = mpmath.exp(x)
+    desired = mpmath_wrightomega(x)
+    return abs(exponential_approx - desired) / desired
+
+
+def main():
+    desired_error = 2 * np.finfo(float).eps
+    print('Series Error')
+    for x in [1e5, 1e10, 1e15, 1e20]:
+        with mpmath.workdps(100):
+            error = wrightomega_series_error(x)
+        print(x, error, error < desired_error)
+
+    print('Exp error')
+    for x in [-10, -25, -50, -100, -200, -400, -700, -740]:
+        with mpmath.workdps(100):
+            error = wrightomega_exp_error(x)
+        print(x, error, error < desired_error)
+
+
+if __name__ == '__main__':
+    main()

.CondaPkg/env/Lib/site-packages/scipy/special/_precompute/zetac.py

@@ -0,0 +1,27 @@
+"""Compute the Taylor series for zeta(x) - 1 around x = 0."""
+try:
+    import mpmath
+except ImportError:
+    pass
+
+
+def zetac_series(N):
+    coeffs = []
+    with mpmath.workdps(100):
+        coeffs.append(-1.5)
+        for n in range(1, N):
+            coeff = mpmath.diff(mpmath.zeta, 0, n)/mpmath.factorial(n)
+            coeffs.append(coeff)
+    return coeffs
+
+
+def main():
+    print(__doc__)
+    coeffs = zetac_series(10)
+    coeffs = [mpmath.nstr(x, 20, min_fixed=0, max_fixed=0)
+              for x in coeffs]
+    print("\n".join(coeffs[::-1]))
+
+
+if __name__ == '__main__':
+    main()