| 12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061 |
- """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
- """
- from __future__ import division, print_function, absolute_import
- import os
- from scipy._lib._numpy_compat import suppress_warnings
- try:
- # Can remove when sympy #11255 is resolved; see
- # https://github.com/sympy/sympy/issues/11255
- with suppress_warnings() as sup:
- sup.filter(DeprecationWarning, "inspect.getargspec.. is deprecated")
- 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()
|