Changes in explore/precision.py [237c9cf:2a602c7] in sasmodels

File:

• explore/precision.py (modified) (5 diffs)

explore/precision.py

@@ -1 +1 @@
 #!/usr/bin/env python
 r"""
-Show numerical precision of $2 J_1(x)/x$.
+Show numerical precision of various expressions.
+
+Evaluates the same function(s) in single and double precision and compares
+the results to 500 digit mpmath evaluation of the same function.
+
+Note: a quick way to generation C and python code for taylor series
+expansions from sympy:
+
+    import sympy as sp
+    x = sp.var("x")
+    f = sp.sin(x)/x
+    t = sp.series(f, n=12).removeO()  # taylor series with no O(x^n) term
+    p = sp.horner(t)   # Horner representation
+    p = p.replace(x**2, sp.var("xsq")  # simplify if alternate terms are zero
+    p = p.n(15)  # evaluate coefficients to 15 digits (optional)
+    c_code = sp.ccode(p, assign_to=sp.var("p"))  # convert to c code
+    py_code = c[:-1]  # strip semicolon to convert c to python
+
+    # mpmath has pade() rational function approximation, which might work
+    # better than the taylor series for some functions:
+    P, Q = mp.pade(sp.Poly(t.n(15),x).coeffs(), L, M)
+    P = sum(a*x**n for n,a in enumerate(reversed(P)))
+    Q = sum(a*x**n for n,a in enumerate(reversed(Q)))
+    c_code = sp.ccode(sp.horner(P)/sp.horner(Q), assign_to=sp.var("p"))
+
+    # There are richardson and shanks series accelerators in both sympy
+    # and mpmath that may be helpful.
 """
 from __future__ import division, print_function
@@ -81 +107 @@
         elif xrange == "linear":
             lin_min, lin_max, lin_steps = 1, 1000, 2000
-            lin_min, lin_max, lin_steps = 0.001, 2, 2000
         elif xrange == "log":
             log_min, log_max, log_steps = -3, 5, 400
@@ -287 +312 @@
 )
 add_function(
+    name="expm1(x)",
+    mp_function=mp.expm1,
+    np_function=np.expm1,
+    ocl_function=make_ocl("return expm1(q);", "sas_expm1"),
+    limits=(-5., 5.),
+)
+add_function(
     name="arctan(x)",
     mp_function=mp.atan,
@@ -322 +354 @@
     np_function=lambda x: np.fmod(x, 2*np.pi),
     ocl_function=make_ocl("return fmod(q, 2*M_PI);", "sas_fmod"),
-)
-add_function(
-    name="debye",
-    mp_function=lambda x: 2*(mp.exp(-x**2) + x**2 - 1)/x**4,
-    np_function=lambda x: 2*(np.expm1(-x**2) + x**2)/x**4,
-    ocl_function=make_ocl("""
-    const double qsq = q*q;
-    if (qsq < 1.0) { // Pade approximation
-        const double x = qsq;
-        if (0) { // 0.36 single
-            // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 4}]
-            return (x*x/180. + 1.)/((1./30.*x + 1./3.)*x + 1);
-        } else if (0) { // 1.0 for single
-            // padeapproximant[2*exp[-x^2] + x^2-1)/x^4, {x, 0, 6}]
-            const double A1=1./24., A2=1./84, A3=-1./3360;
-            const double B1=3./8., B2=3./56., B3=1./336.;
-            return (((A3*x + A2)*x + A1)*x + 1.)/(((B3*x + B2)*x + B1)*x + 1.);
-        } else if (1) { // 1.0 for single, 0.25 for double
-            // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}]
-            const double A1=1./15., A2=1./60, A3=0., A4=1./75600.;
-            const double B1=2./5., B2=1./15., B3=1./180., B4=1./5040.;
-            return ((((A4*x + A3)*x + A2)*x + A1)*x + 1.)
-                  /((((B4*x + B3)*x + B2)*x + B1)*x + 1.);
-        } else { // 1.0 for single, 0.5 for double
-            // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}]
-            const double A1=1./12., A2=2./99., A3=1./2640., A4=1./23760., A5=-1./1995840.;
-            const double B1=5./12., B2=5./66., B3=1./132., B4=1./2376., B5=1./95040.;
-            return (((((A5*x + A4)*x + A3)*x + A2)*x + A1)*x + 1.)
-                  /(((((B5*x + B4)*x + B3)*x + B2)*x + B1)*x + 1.);
-        }
-    } else if (qsq < 1.) { // Taylor series; 0.9 for single, 0.25 for double
-        const double x = qsq;
-        const double C0 = +1.;
-        const double C1 = -1./3.;
-        const double C2 = +1./12.;
-        const double C3 = -1./60.;
-        const double C4 = +1./360.;
-        const double C5 = -1./2520.;
-        const double C6 = +1./20160.;
-        const double C7 = -1./181440.;
-        //return ((((C5*x + C4)*x + C3)*x + C2)*x + C1)*x + C0;
-        //return (((((C6*x + C5)*x + C4)*x + C3)*x + C2)*x + C1)*x + C0;
-        return ((((((C7*x + C6)*x + C5)*x + C4)*x + C3)*x + C2)*x + C1)*x + C0;
-    } else {
-        return 2.*(expm1(-qsq) + qsq)/(qsq*qsq);
-    }
-    """, "sas_debye"),
 )
 
@@ -448 +433 @@
 )
 
+replacement_expm1 = """\
+      double x = (double)q;  // go back to float for single precision kernels
+      // Adapted from the cephes math library.
+      // Copyright 1984 - 1992 by Stephen L. Moshier
+      if (x != x || x == 0.0) {
+         return x; // NaN and +/- 0
+      } else if (x < -0.5 || x > 0.5) {
+         return exp(x) - 1.0;
+      } else {
+         const double xsq = x*x;
+         const double p = (((
+            +1.2617719307481059087798E-4)*xsq
+            +3.0299440770744196129956E-2)*xsq
+            +9.9999999999999999991025E-1);
+         const double q = ((((
+            +3.0019850513866445504159E-6)*xsq
+            +2.5244834034968410419224E-3)*xsq
+            +2.2726554820815502876593E-1)*xsq
+            +2.0000000000000000000897E0);
+         double r = x * p;
+         r =  r / (q - r);
+         return r+r;
+       }
+"""
+add_function(
+    name="sas_expm1(x)",
+    mp_function=mp.expm1,
+    np_function=np.expm1,
+    ocl_function=make_ocl(replacement_expm1, "sas_expm1"),
+)
+
 # Alternate versions of 3 j1(x)/x, for posterity
 def taylor_3j1x_x(x):