Changeset 2a602c7 in sasmodels

Timestamp:

Oct 17, 2017 3:52:01 PM (7 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

becded3

Parents:

57eb6a4

Message:

check the single precision accuracy on cephes expm1() replacement function in kernel_header.c

File:

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

explore/precision.py

@@ -309 +309 @@
     np_function=scipy.special.erfc,
     ocl_function=make_ocl("return sas_erfc(q);", "sas_erfc", ["lib/polevl.c", "lib/sas_erf.c"]),
+    limits=(-5., 5.),
+)
+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.),
 )
@@ -426 +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):