Changeset 3f8584a2 in sasmodels


Ignore:
Timestamp:
Jul 27, 2016 5:53:07 AM (8 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:
c047acf
Parents:
56547a8
Message:

clean separation between float/double bessel functions

Location:
sasmodels/models/lib
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/lib/sas_J0.c

    r0278e3f r3f8584a2  
    5353 * except YP, YQ which are designed for absolute error. */ 
    5454 
     55#if FLOAT_SIZE>4 
     56//Cephes double precission 
    5557 
    5658 constant double PPJ0[8] = { 
    … …  
    144146  }; 
    145147 
    146  constant double MOJ0[8] = { 
     148double j0(double x); 
     149double j0(double x) 
     150{ 
     151    double w, z, p, q, xn; 
     152 
     153    //const double TWOOPI = 6.36619772367581343075535E-1; 
     154    const double SQ2OPI = 7.9788456080286535587989E-1; 
     155    const double PIO4 = 7.85398163397448309616E-1; 
     156 
     157    const double DR1 = 5.78318596294678452118E0; 
     158    const double DR2 = 3.04712623436620863991E1; 
     159 
     160 
     161    if( x < 0 ) 
     162        x = -x; 
     163 
     164    if( x <= 5.0 ) { 
     165        z = x * x; 
     166        if( x < 1.0e-5 ) 
     167            return( 1.0 - z/4.0 ); 
     168 
     169        p = (z - DR1) * (z - DR2); 
     170        p = p * polevl( z, RPJ0, 3)/p1evl( z, RQJ0, 8 ); 
     171        return( p ); 
     172    } 
     173 
     174    w = 5.0/x; 
     175    q = 25.0/(x*x); 
     176    p = polevl( q, PPJ0, 6)/polevl( q, PQJ0, 6 ); 
     177    q = polevl( q, QPJ0, 7)/p1evl( q, QQJ0, 7 ); 
     178    xn = x - PIO4; 
     179 
     180    double sn, cn; 
     181    SINCOS(xn, sn, cn); 
     182    p = p * cn - w * q * sn; 
     183 
     184    return( p * SQ2OPI / sqrt(x) ); 
     185} 
     186#else 
     187 
     188 constant float MOJ0[8] = { 
    147189        -6.838999669318810E-002, 
    148190        1.864949361379502E-001, 
    … …  
    155197  }; 
    156198 
    157  constant double PHJ0[8] = { 
     199 constant float PHJ0[8] = { 
    158200        3.242077816988247E+001, 
    159201        -3.630592630518434E+001, 
    … …  
    166208  }; 
    167209 
    168  constant double JPJ0[8] = { 
     210 constant float JPJ0[8] = { 
    169211        -6.068350350393235E-008, 
    170212        6.388945720783375E-006, 
    … …  
    177219 }; 
    178220 
    179 double sas_J0(double x); 
    180 double sas_J0(double x) 
     221//Cephes single precission 
     222float j0f(float x); 
     223float j0f(float x) 
    181224{ 
    182  
    183 //Cephes single precission 
    184 #if FLOAT_SIZE>4 
    185     double w, z, p, q, xn; 
    186  
    187     //const double TWOOPI = 6.36619772367581343075535E-1; 
    188     const double SQ2OPI = 7.9788456080286535587989E-1; 
    189     const double PIO4 = 7.85398163397448309616E-1; 
    190  
    191     const double DR1 = 5.78318596294678452118E0; 
    192     const double DR2 = 3.04712623436620863991E1; 
    193  
    194  
    195     if( x < 0 ) 
    196             x = -x; 
    197  
    198     if( x <= 5.0 ) { 
    199             z = x * x; 
    200             if( x < 1.0e-5 ) 
    201                     return( 1.0 - z/4.0 ); 
    202  
    203             p = (z - DR1) * (z - DR2); 
    204             p = p * polevl( z, RPJ0, 3)/p1evl( z, RQJ0, 8 ); 
    205             return( p ); 
    206         } 
    207  
    208     w = 5.0/x; 
    209     q = 25.0/(x*x); 
    210     p = polevl( q, PPJ0, 6)/polevl( q, PQJ0, 6 ); 
    211     q = polevl( q, QPJ0, 7)/p1evl( q, QQJ0, 7 ); 
    212     xn = x - PIO4; 
    213  
    214     double sn, cn; 
    215     SINCOS(xn, sn, cn); 
    216     p = p * cn - w * q * sn; 
    217  
    218     return( p * SQ2OPI / sqrt(x) ); 
    219 //Cephes single precission 
    220 #else 
    221     double xx, w, z, p, q, xn; 
     225    float xx, w, z, p, q, xn; 
    222226 
    223227    //const double YZ1 =  0.43221455686510834878; 
    224228    //const double YZ2 = 22.401876406482861405; 
    225229    //const double YZ3 = 64.130620282338755553; 
    226     const double DR1 =  5.78318596294678452118; 
    227     const double PIO4F = 0.7853981633974483096; 
     230    const float DR1 =  5.78318596294678452118; 
     231    const float PIO4F = 0.7853981633974483096; 
    228232 
    229233    if( x < 0 ) 
    230             xx = -x; 
     234        xx = -x; 
    231235    else 
    232             xx = x; 
     236        xx = x; 
    233237 
    234238    if( x <= 2.0 ) { 
    235             z = xx * xx; 
    236             if( x < 1.0e-3 ) 
    237                     return( 1.0 - 0.25*z ); 
    238  
    239             p = (z-DR1) * polevl( z, JPJ0, 4); 
    240             return( p ); 
    241         } 
     239        z = xx * xx; 
     240        if( x < 1.0e-3 ) 
     241            return( 1.0 - 0.25*z ); 
     242 
     243        p = (z-DR1) * polevl( z, JPJ0, 4); 
     244        return( p ); 
     245    } 
    242246 
    243247    q = 1.0/x; 
    … …  
    249253    p = p * cos(xn + xx); 
    250254    return(p); 
     255} 
    251256#endif 
    252257 
    253 } 
    254  
    255  
     258#if FLOAT_SIZE>4 
     259#define sas_J0 j0 
     260#else 
     261#define sas_J0 j0f 
     262#endif 

  • sasmodels/models/lib/sas_J1.c

    r0278e3f r3f8584a2  
    4040*/ 
    4141 
     42#if FLOAT_SIZE>4 
     43 
     44//Cephes double pression function 
    4245constant double RPJ1[8] = { 
    4346    -8.99971225705559398224E8, 
    … …  
    102105    0.0 }; 
    103106 
    104 constant double JPJ1[8] = { 
     107double j1(double x); 
     108double j1(double x) 
     109{ 
     110 
     111    double w, z, p, q, xn; 
     112 
     113    const double Z1 = 1.46819706421238932572E1; 
     114    const double Z2 = 4.92184563216946036703E1; 
     115    const double THPIO4 =  2.35619449019234492885; 
     116    const double SQ2OPI = 0.79788456080286535588; 
     117 
     118    w = x; 
     119    if( x < 0 ) 
     120        w = -x; 
     121 
     122    if( w <= 5.0 ) { 
     123        z = x * x; 
     124        w = polevl( z, RPJ1, 3 ) / p1evl( z, RQJ1, 8 ); 
     125        w = w * x * (z - Z1) * (z - Z2); 
     126        return( w ); 
     127    } 
     128 
     129    w = 5.0/x; 
     130    z = w * w; 
     131 
     132    p = polevl( z, PPJ1, 6)/polevl( z, PQJ1, 6 ); 
     133    q = polevl( z, QPJ1, 7)/p1evl( z, QQJ1, 7 ); 
     134 
     135    xn = x - THPIO4; 
     136 
     137    double sn, cn; 
     138    SINCOS(xn, sn, cn); 
     139    p = p * cn - w * q * sn; 
     140 
     141    return( p * SQ2OPI / sqrt(x) ); 
     142} 
     143 
     144#else 
     145//Single precission version of cephes 
     146constant float JPJ1[8] = { 
    105147    -4.878788132172128E-009, 
    106148    6.009061827883699E-007, 
    … …  
    113155    }; 
    114156 
    115 constant double MO1J1[8] = { 
     157constant float MO1J1[8] = { 
    116158    6.913942741265801E-002, 
    117159    -2.284801500053359E-001, 
    … …  
    124166    }; 
    125167 
    126 constant double PH1J1[8] = { 
     168constant float PH1J1[8] = { 
    127169    -4.497014141919556E+001, 
    128170    5.073465654089319E+001, 
    … …  
    135177    }; 
    136178 
    137 double sas_J1(double x); 
    138 double sas_J1(double x) 
     179float j1f(float x); 
     180float j1f(float x) 
    139181{ 
    140182 
    141 //Cephes double pression function 
    142 #if FLOAT_SIZE>4 
    143  
    144     double w, z, p, q, xn; 
    145  
    146     const double Z1 = 1.46819706421238932572E1; 
    147     const double Z2 = 4.92184563216946036703E1; 
    148     const double THPIO4 =  2.35619449019234492885; 
    149     const double SQ2OPI = 0.79788456080286535588; 
    150  
    151     w = x; 
    152     if( x < 0 ) 
    153             w = -x; 
    154  
    155     if( w <= 5.0 ) 
    156         { 
    157             z = x * x; 
    158             w = polevl( z, RPJ1, 3 ) / p1evl( z, RQJ1, 8 ); 
    159             w = w * x * (z - Z1) * (z - Z2); 
    160             return( w ); 
    161         } 
    162  
    163     w = 5.0/x; 
    164     z = w * w; 
    165  
    166     p = polevl( z, PPJ1, 6)/polevl( z, PQJ1, 6 ); 
    167     q = polevl( z, QPJ1, 7)/p1evl( z, QQJ1, 7 ); 
    168  
    169     xn = x - THPIO4; 
    170  
    171     double sn, cn; 
    172     SINCOS(xn, sn, cn); 
    173     p = p * cn - w * q * sn; 
    174  
    175     return( p * SQ2OPI / sqrt(x) ); 
    176  
    177  
    178 //Single precission version of cephes 
    179 #else 
    180     double xx, w, z, p, q, xn; 
    181  
    182     const double Z1 = 1.46819706421238932572E1; 
    183     const double THPIO4F =  2.35619449019234492885;    /* 3*pi/4 */ 
     183    float xx, w, z, p, q, xn; 
     184 
     185    const float Z1 = 1.46819706421238932572E1; 
     186    const float THPIO4F =  2.35619449019234492885;    /* 3*pi/4 */ 
    184187 
    185188 
    186189    xx = x; 
    187190    if( xx < 0 ) 
    188             xx = -x; 
    189  
    190     if( xx <= 2.0 ) 
    191         { 
    192             z = xx * xx; 
    193             p = (z-Z1) * xx * polevl( z, JPJ1, 4 ); 
    194             return( p ); 
    195         } 
     191        xx = -x; 
     192 
     193    if( xx <= 2.0 ) { 
     194        z = xx * xx; 
     195        p = (z-Z1) * xx * polevl( z, JPJ1, 4 ); 
     196        return( p ); 
     197    } 
    196198 
    197199    q = 1.0/x; 
    … …  
    204206 
    205207    return(p); 
     208} 
    206209#endif 
    207 } 
     210 
     211#if FLOAT_SIZE>4 
     212#define sas_J1 j1 
     213#else 
     214#define sas_J1 j1f 
     215#endif 
    208216 
    209217//Finally J1c function that equals 2*J1(x)/x 

  • sasmodels/models/lib/sas_JN.c

    r4937980 r3f8584a2  
    4848*/ 
    4949 
    50 double sas_JN( int n, double x ); 
    51  
    52 double sas_JN( int n, double x ) { 
    53  
     50#if FLOAT_SIZE > 4 
     51 
     52double jn( int n, double x ); 
     53double jn( int n, double x ) { 
     54 
     55    // PAK: seems to be machine epsilon/2 
    5456    const double MACHEP = 1.11022302462515654042E-16; 
    5557    double pkm2, pkm1, pk, xk, r, ans; 
    … …  
    5759 
    5860    if( n < 0 ) { 
    59             n = -n; 
    60             if( (n & 1) == 0 )  /* -1**n */ 
    61                     sign = 1; 
    62             else 
    63                     sign = -1; 
    64         } 
    65     else 
    66             sign = 1; 
     61        n = -n; 
     62        if( (n & 1) == 0 )      /* -1**n */ 
     63            sign = 1; 
     64        else 
     65            sign = -1; 
     66    } else { 
     67        sign = 1; 
     68    } 
    6769 
    6870    if( x < 0.0 ) { 
    69             if( n & 1 ) 
    70                     sign = -sign; 
    71             x = -x; 
    72         } 
     71        if( n & 1 ) 
     72            sign = -sign; 
     73    x = -x; 
     74    } 
    7375 
    7476    if( n == 0 ) 
    75             return( sign * sas_J0(x) ); 
     77        return( sign * j0(x) ); 
    7678    if( n == 1 ) 
    77             return( sign * sas_J1(x) ); 
     79        return( sign * j1(x) ); 
    7880    if( n == 2 ) 
    79             return( sign * (2.0 * sas_J1(x) / x  -  sas_J0(x)) ); 
     81        return( sign * (2.0 * j1(x) / x  -  j0(x)) ); 
    8082 
    8183    if( x < MACHEP ) 
    82             return( 0.0 ); 
    83  
    84     #if FLOAT_SIZE > 4 
    85         k = 53; 
    86     #else 
    87         k = 24; 
    88     #endif 
    89  
     84        return( 0.0 ); 
     85 
     86    k = 53; 
    9087    pk = 2 * (n + k); 
    9188    ans = pk; 
    … …  
    9390 
    9491    do { 
    95             pk -= 2.0; 
    96             ans = pk - (xk/ans); 
    97         } while( --k > 0 ); 
     92        pk -= 2.0; 
     93        ans = pk - (xk/ans); 
     94    } while( --k > 0 ); 
    9895 
    9996    /* backward recurrence */ 
    … …  
    10198    pk = 1.0; 
    10299 
    103     #if FLOAT_SIZE > 4 
    104         ans = x/ans; 
    105         pkm1 = 1.0/ans; 
    106  
    107         k = n-1; 
    108         r = 2 * k; 
    109  
    110         do { 
    111             pkm2 = (pkm1 * r  -  pk * x) / x; 
    112                 pk = pkm1; 
    113                 pkm1 = pkm2; 
    114                 r -= 2.0; 
    115             } while( --k > 0 ); 
    116  
    117         if( fabs(pk) > fabs(pkm1) ) 
    118                 ans = sas_J1(x)/pk; 
     100    ans = x/ans; 
     101    pkm1 = 1.0/ans; 
     102 
     103    k = n-1; 
     104    r = 2 * k; 
     105 
     106    do { 
     107        pkm2 = (pkm1 * r  -  pk * x) / x; 
     108        pk = pkm1; 
     109        pkm1 = pkm2; 
     110        r -= 2.0; 
     111    } while( --k > 0 ); 
     112 
     113    if( fabs(pk) > fabs(pkm1) ) 
     114        ans = j1(x)/pk; 
     115    else 
     116        ans = j0(x)/pkm1; 
     117 
     118    return( sign * ans ); 
     119} 
     120 
     121#else 
     122 
     123float jnf(int n, float x) 
     124{ 
     125    // PAK: seems to be machine epsilon/2 
     126    const double MACHEP = 5.9604645e-08; 
     127    float pkm2, pkm1, pk, xk, r, ans; 
     128    int k, sign; 
     129 
     130    if( n < 0 ) { 
     131        n = -n; 
     132        if( (n & 1) == 0 ) /* -1**n */ 
     133            sign = 1; 
    119134        else 
    120                 ans = sas_J0(x)/pkm1; 
    121  
    122             return( sign * ans ); 
    123  
    124     #else 
    125         const double xinv = 1.0/x; 
    126         pkm1 = ans * xinv; 
    127         k = n-1; 
    128         r = (float )(2 * k); 
    129  
    130         do { 
    131                 pkm2 = (pkm1 * r  -  pk * x) * xinv; 
    132                 pk = pkm1; 
    133                 pkm1 = pkm2; 
    134                 r -= 2.0; 
    135             } 
    136         while( --k > 0 ); 
    137  
    138         r = pk; 
    139         if( r < 0 ) 
    140                 r = -r; 
    141         ans = pkm1; 
    142         if( ans < 0 ) 
    143                 ans = -ans; 
    144  
    145         if( r > ans )  /* if( fabs(pk) > fabs(pkm1) ) */ 
    146                 ans = sign * sas_J1(x)/pk; 
    147         else 
    148                 ans = sign * sas_J0(x)/pkm1; 
    149         return( ans ); 
    150     #endif 
     135            sign = -1; 
     136    } else { 
     137        sign = 1; 
     138    } 
     139 
     140    if( x < 0.0 ) { 
     141        if( n & 1 ) 
     142            sign = -sign; 
     143        x = -x; 
     144    } 
     145 
     146    if( n == 0 ) 
     147        return( sign * j0f(x) ); 
     148    if( n == 1 ) 
     149        return( sign * j1f(x) ); 
     150    if( n == 2 ) 
     151        return( sign * (2.0 * j1f(x) / x  -  j0f(x)) ); 
     152 
     153    if( x < MACHEP ) 
     154        return( 0.0 ); 
     155 
     156    k = 24; 
     157    pk = 2 * (n + k); 
     158    ans = pk; 
     159    xk = x * x; 
     160 
     161    do { 
     162        pk -= 2.0; 
     163        ans = pk - (xk/ans); 
     164    } while( --k > 0 ); 
     165 
     166    /* backward recurrence */ 
     167 
     168    pk = 1.0; 
     169 
     170    const float xinv = 1.0/x; 
     171    pkm1 = ans * xinv; 
     172    k = n-1; 
     173    r = (float )(2 * k); 
     174 
     175    do { 
     176        pkm2 = (pkm1 * r  -  pk * x) * xinv; 
     177        pk = pkm1; 
     178        pkm1 = pkm2; 
     179        r -= 2.0; 
     180    } while( --k > 0 ); 
     181 
     182    r = pk; 
     183    if( r < 0 ) 
     184        r = -r; 
     185    ans = pkm1; 
     186    if( ans < 0 ) 
     187        ans = -ans; 
     188 
     189    if( r > ans )  /* if( fabs(pk) > fabs(pkm1) ) */ 
     190        ans = sign * j1f(x)/pk; 
     191    else 
     192        ans = sign * j0f(x)/pkm1; 
     193    return( ans ); 
    151194} 
    152  
     195#endif 
     196 
     197#if FLOAT_SIZE>4 
     198#define sas_JN jn 
     199#else 
     200#define sas_JN jnf 
     201#endif 
Note: See TracChangeset for help on using the changeset viewer.