Changeset 2a602c7 in sasmodels


Ignore:
Timestamp:
Oct 17, 2017 3:52:01 PM (6 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:
1 edited

Legend:

Unmodified
Added
Removed
  • explore/precision.py

    r57eb6a4 r2a602c7  
    309309    np_function=scipy.special.erfc, 
    310310    ocl_function=make_ocl("return sas_erfc(q);", "sas_erfc", ["lib/polevl.c", "lib/sas_erf.c"]), 
     311    limits=(-5., 5.), 
     312) 
     313add_function( 
     314    name="expm1(x)", 
     315    mp_function=mp.expm1, 
     316    np_function=np.expm1, 
     317    ocl_function=make_ocl("return expm1(q);", "sas_expm1"), 
    311318    limits=(-5., 5.), 
    312319) 
     
    426433) 
    427434 
     435replacement_expm1 = """\ 
     436      double x = (double)q;  // go back to float for single precision kernels 
     437      // Adapted from the cephes math library. 
     438      // Copyright 1984 - 1992 by Stephen L. Moshier 
     439      if (x != x || x == 0.0) { 
     440         return x; // NaN and +/- 0 
     441      } else if (x < -0.5 || x > 0.5) { 
     442         return exp(x) - 1.0; 
     443      } else { 
     444         const double xsq = x*x; 
     445         const double p = ((( 
     446            +1.2617719307481059087798E-4)*xsq 
     447            +3.0299440770744196129956E-2)*xsq 
     448            +9.9999999999999999991025E-1); 
     449         const double q = (((( 
     450            +3.0019850513866445504159E-6)*xsq 
     451            +2.5244834034968410419224E-3)*xsq 
     452            +2.2726554820815502876593E-1)*xsq 
     453            +2.0000000000000000000897E0); 
     454         double r = x * p; 
     455         r =  r / (q - r); 
     456         return r+r; 
     457       } 
     458""" 
     459add_function( 
     460    name="sas_expm1(x)", 
     461    mp_function=mp.expm1, 
     462    np_function=np.expm1, 
     463    ocl_function=make_ocl(replacement_expm1, "sas_expm1"), 
     464) 
     465 
    428466# Alternate versions of 3 j1(x)/x, for posterity 
    429467def taylor_3j1x_x(x): 
Note: See TracChangeset for help on using the changeset viewer.