# Changeset e6f1410 in sasmodels

Ignore:
Timestamp:
Mar 6, 2016 8:16:29 PM (7 years ago)
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:
8dca856
Parents:
Message:

examine numerical precision of J1c, sinc, j1c and log gamma

Files:
2 edited
1 moved

Unmodified
Removed
• ## sasmodels/models/lib/lanczos_gamma.c

 r5c1d341 double lanczos_gamma(double q) { // Lanczos approximation to the Gamma function. // Lanczos approximation to the Log Gamma function. double x,y,tmp,ser;
• ## sasmodels/models/lib/sph_j1c.c

 rc437dbb * Note that the values differ from sasview ~ 5e-12 rather than 5e-14, but * in this case it is likely cancellation errors in the original expression * using double precision that are the source.  Single precision only * requires the first 3 terms.  Double precision requires the 4th term. * The fifth term is not needed, and is commented out. * Taylor expansion: *      1.0 + q2*(-3./30. + q2*(3./840.))+ q2*(-3./45360. + q2*(3./3991680.)))) * using double precision that are the source. */ // The choice of the number of terms in the series and the cutoff value for // switching between series and direct calculation depends on the numeric // precision. // // Point where direct calculation reaches machine precision: // //   single machine precision eps 3e-8 at qr=1.1 ** //   double machine precision eps 4e-16 at qr=1.1 // // Point where Taylor series reaches machine precision (eps), where taylor // series matches direct calculation (cross) and the error at that point: // //   prec   n eps  cross  error //   single 3 0.28  0.4   6.2e-7 //   single 4 0.68  0.7   2.3e-7 //   single 5 1.18  1.2   7.5e-8 //   double 3 0.01  0.03  2.3e-13 //   double 4 0.06  0.1   3.1e-14 //   double 5 0.16  0.2   5.0e-15 // // ** Note: relative error on single precision starts increase on the direct // method at qr=1.1, rising from 3e-8 to 5e-5 by qr=1e3.  This should be // safe for the sans range, with objects of 100 nm supported to a q of 0.1 // while maintaining 5 digits of precision.  For usans/sesans, the objects // are larger but the q is smaller, so again it should be fine. // // See explore/sph_j1c.py for code to explore these ranges. // Use 4th order series #if FLOAT_SIZE>4 #define SPH_J1C_CUTOFF 0.1 #else #define SPH_J1C_CUTOFF 0.7 #endif double sph_j1c(double q); double sph_j1c(double q) { const double q2 = q*q; double sin_q, cos_q; SINCOS(q, sin_q, cos_q); #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 // 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); if (q < SPH_J1C_CUTOFF) { const double q2 = q*q; return (1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))));// + q2*(3./3991680.))))); } else { // 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); double sin_q, cos_q; SINCOS(q, sin_q, cos_q); return 3.0*(sin_q/q - cos_q)/(q*q); } return ret; */ }
Note: See TracChangeset for help on using the changeset viewer.