expn_asy.py 1.6 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061
  1. """Precompute the polynomials for the asymptotic expansion of the
  2. generalized exponential integral.
  3. Sources
  4. -------
  5. [1] NIST, Digital Library of Mathematical Functions,
  6. https://dlmf.nist.gov/8.20#ii
  7. """
  8. from __future__ import division, print_function, absolute_import
  9. import os
  10. from scipy._lib._numpy_compat import suppress_warnings
  11. try:
  12. # Can remove when sympy #11255 is resolved; see
  13. # https://github.com/sympy/sympy/issues/11255
  14. with suppress_warnings() as sup:
  15. sup.filter(DeprecationWarning, "inspect.getargspec.. is deprecated")
  16. import sympy
  17. from sympy import Poly
  18. x = sympy.symbols('x')
  19. except ImportError:
  20. pass
  21. def generate_A(K):
  22. A = [Poly(1, x)]
  23. for k in range(K):
  24. A.append(Poly(1 - 2*k*x, x)*A[k] + Poly(x*(x + 1))*A[k].diff())
  25. return A
  26. WARNING = """\
  27. /* This file was automatically generated by _precompute/expn_asy.py.
  28. * Do not edit it manually!
  29. */
  30. """
  31. def main():
  32. print(__doc__)
  33. fn = os.path.join('..', 'cephes', 'expn.h')
  34. K = 12
  35. A = generate_A(K)
  36. with open(fn + '.new', 'w') as f:
  37. f.write(WARNING)
  38. f.write("#define nA {}\n".format(len(A)))
  39. for k, Ak in enumerate(A):
  40. tmp = ', '.join([str(x.evalf(18)) for x in Ak.coeffs()])
  41. f.write("static const double A{}[] = {{{}}};\n".format(k, tmp))
  42. tmp = ", ".join(["A{}".format(k) for k in range(K + 1)])
  43. f.write("static const double *A[] = {{{}}};\n".format(tmp))
  44. tmp = ", ".join([str(Ak.degree()) for Ak in A])
  45. f.write("static const int Adegs[] = {{{}}};\n".format(tmp))
  46. os.rename(fn + '.new', fn)
  47. if __name__ == "__main__":
  48. main()