Changeset 3a220e6 in sasmodels

Timestamp:

Oct 3, 2017 1:56:55 PM (7 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

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

Children:

237c9cf

Parents:

6e72989

Message:

tweak cutoff point for debye taylor series expansion

Files:

• explore/precision.py (modified) (3 diffs)
• sasmodels/models/lib/wrc_cyl.c (modified) (2 diffs)

explore/precision.py

@@ -81 +81 @@
         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
@@ -328 +329 @@
     ocl_function=make_ocl("""
     const double qsq = q*q;
-    if (qsq < 0.1) {
+    if (qsq < -0.75*0.75) { // Pade approximation around 0
+        const double x = qsq;
+        if (0) {
+            // 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) {
+            // padeapproximant[2*exp[-x^2] + x^2-1)/x^4, {x, 0, 6}]
+            // works for single precision, with qsq < 1
+            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 (0) {
+            // 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 {
+            // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}]
+            // works for double precision, with qsq < 9/16
+            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 < 0.1) { // Taylor series around 0.
         const double x = qsq;
         const double C0 = +1.;
@@ -339 +366 @@
         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 if (qsq < 1.0) {
-        // taylor series around q^2 = 0.5
-        const double x = qsq - 0.5;
-        const double sqrt_e = sqrt(M_E);
-        const double C0 = 8./sqrt_e - 4.;
-        const double C1 = 24. - 40./sqrt_e;
-        const double C2 = 132./sqrt_e - 80.;
-        const double C3 = 224. - 1108./(3.*sqrt_e);
-        const double C4 = 2849./(3.*sqrt_e) - 576.;
-        const double C5 = 1408 - 11607./(5.*sqrt_e);
-        return ((((C5*x + C4)*x + C3)*x + C2)*x + C1)*x + C0;
-    } */ else {
+        //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.*(exp(-qsq) + qsq - 1.)/(qsq*qsq);
     }

sasmodels/models/lib/wrc_cyl.c

@@ -121 +121 @@
 {
 #if FLOAT_SIZE>4
-#define DEBYE_CUTOFF 0.01  // 1e-13 error
+#define DEBYE_CUTOFF 0.1  // 1e-14 error
 #else
-#define DEBYE_CUTOFF 1.0  // 1e-7 error
+#define DEBYE_CUTOFF 0.9  // 4e-7 error
 #endif
+
+/* For double precision, the following gets 1e-15 error rather than 1e-14
+    if (qsq < 9./16.) {
+        // 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.;
+        const double x = qsq;
+        return (((((A5*x + A4)*x + A3)*x + A2)*x + A1)*x + 1.)
+                /(((((B5*x + B4)*x + B3)*x + B2)*x + B1)*x + 1.);
+    }
+*/
 
     if (qsq < DEBYE_CUTOFF) {
@@ -137 +148 @@
         const double C7 = -1./181440.;
         return ((((((C7*x + C6)*x + C5)*x + C4)*x + C3)*x + C2)*x + C1)*x + C0;
-    /* failed attempt to improve double precision better than 1e-13
-    } else if (qsq < 1.0) {
-        // taylor series around q^2 = 0.5
-        const double x = qsq - 0.5;
-        const double sqrt_e = sqrt(M_E);
-        const double C0 = 8./sqrt_e - 4.;
-        const double C1 = 24. - 40./sqrt_e;
-        const double C2 = 132./sqrt_e - 80.;
-        const double C3 = 224. - 1108./(3.*sqrt_e);
-        const double C4 = 2849./(3.*sqrt_e) - 576.;
-        const double C5 = 1408 - 11607./(5.*sqrt_e);
-        return ((((C5*x + C4)*x + C3)*x + C2)*x + C1)*x + C0;
-    */
     } else {
         return 2.*(exp(-qsq) + qsq - 1.)/(qsq*qsq);