Changeset c437dbb in sasmodels

Timestamp:

Mar 1, 2016 1:40:08 PM (9 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

ce0b154

Parents:

cf85329

Message:

Improve accuracy of sph_j1c to 7 digits single, 14 digits double

Location:

sasmodels

Files:

• generate.py (modified) (2 diffs)
• models/lib/sph_j1c.c (modified) (2 diffs)

sasmodels/generate.py

@@ -193 +193 @@
 and adds the parameter table to the top.  The function :func:`model_sources`
 returns a list of files required by the model.
+
+Code follows the C99 standard with the following extensions and conditions::
+
+    M_PI_180 = pi/180
+    M_4PI_3 = 4pi/3
+    square(x) = x*x
+    cube(x) = x*x*x
+    sinc(x) = sin(x)/x, with sin(0)/0 -> 1
+    all double precision constants must include the decimal point
+    all double declarations may be converted to half, float, or long double
+    FLOAT_SIZE is the number of bytes in the converted variables
 """
 from __future__ import print_function
@@ -344 +355 @@
     """
     if dtype == F16:
+        fbytes = 2
         source = _F16_PRAGMA + _convert_type(source, "half", "f")
     elif dtype == F32:
+        fbytes = 4
         source = _convert_type(source, "float", "f")
     elif dtype == F64:
+        fbytes = 8
         source = _F64_PRAGMA + source  # Source is already double
     elif dtype == F128:
+        fbytes = 16
         source = _convert_type(source, "long double", "L")
     else:
         raise ValueError("Unexpected dtype in source conversion: %s"%dtype)
-    return source
+    return ("#define FLOAT_SIZE %d\n"%fbytes)+source
 
 

sasmodels/models/lib/sph_j1c.c

@@ -10 +10 @@
 * Taylor expansion:
 *      1.0 + q2*(-3./30. + q2*(3./840.))+ q2*(-3./45360. + q2*(3./3991680.))))
-* Expression returned from Herbie (herbie.uwpise.org/demo):
-*      const double t = ((1. + 3.*q2*q2/5600.) - q2/20.);
-*      return t*t;
 */
 
@@ -23 +20 @@
     SINCOS(q, sin_q, cos_q);
 
-    const double bessel = (q < 0.384038453352533)
-        ? (1.0 + q2*(-3./30. + q2*(3./840.)))
+#if FLOAT_SIZE>4
+// For double precision, choose a cutoff of 0.18, which is the lower limit
+// for the trig expression to 14 digits.  If only 12 digits are needed, then
+// only 4 terms of the Taylor expansion are needed.
+#define CUTOFF 0.18
+#else
+// For single precision, choose a cutoff halfway between the single precision
+// lower limit for the trig expression of 1.03 and the upper limit of 1.3
+// for the Taylor series expansion with 5 terms (or 9 if you count the zeros).
+// For 5 digits of precision, you can drop two terms of the taylor series
+// and choose a cutoff of 0.3.
+#define CUTOFF 1.15
+#endif
+
+    const double bessel = (q < CUTOFF)
+        ? (1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360. + q2*(3./3991680.)))))
         : 3.0*(sin_q/q - cos_q)/q2;
-
     return bessel;
 
- /*
+/*
     // Code to test various expressions
-    if (sizeof(q2) > 4) {
-        return 3.0*(sin_q/q - cos_q)/q2;
-    } else if (q < 0.384038453352533) {
-        //const double t = ((1. + 3.*q2*q2/5600.) - q2/20.); return t*t;
-        return 1.0 + q2*q2*(3./840.) - q2*(3./30.);
-        //return 1.0 + q2*(-3./30. + q2*(3./840.));
-        //return 1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.)));
-        //return 1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360. + q2*(3./3991680.))));
+    // tested using sympy.mpmath.mp with
+    //    mp.dps = 50 (which is good for q down to 1e-6)
+    //    def j1c(x): return 3*(mp.sin(x)/x - mp.cos(x))/x**2
+    double ret;
+    if (sizeof(q2) > 8) {
+
+        ret = 3.0*(sinl(q)/q - cosl(q))/q2;
+printf("%.15Lf %.15Lg\n", q, ret);
+    //} else if (q < 0.384038453352533) {
+    //} else if (q < 0.18) {
+    } else if (q < 18e-10) {
+        // NB: all are good to 5 digits below 0.18f
+        // good is defined as 14 digits for double and 7 for single
+        //ret = 1.0 + q2*q2*(3./840.) - q2*(3./30.); // good below 0.02d 0.34f
+        ret = square((1. + 3./5600.*q2*q2) - q2/20.); // good below 0.03d 0.08f
+        //ret = 1.0 - q2*(3./30.); // good below 0.02d 0.07f
+        //ret = 1.0 + q2*(-3./30. + q2*(3./840.)); // good below 0.02d 0.34f
+        //ret = 1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))); // good below 0.1d 0.8f, 12 digits below 0.18d
+        //ret = 1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360. + q2*(3./3991680.)))); // good below 0.18d 1.3f
+printf("%.15g %.15g\n", q, ret);
     } else {
-        return 3.0*(sin_q/q - cos_q)/q2;
+        // NB: can use a cutoff of 0.1f if the goal is 5 digits rather than 7
+        ret = 3.0*(sin_q/q - cos_q)/q2; // good above 0.18d, 1.03f
+printf("%.15g %.15g\n", q, ret);
     }
+    return ret;
 */
 }