Changeset 0a9d219 in sasmodels


Ignore:
Timestamp:
Mar 17, 2016 8:45:24 AM (8 years ago)
Author:
wojciech
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:
16afd49
Parents:
3936ad3
Message:

Single precisson cephes version has been added

Location:
sasmodels/models/lib
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/lib/j0d.c

    r3936ad3 r0a9d219  
    9696double j0( double ); 
    9797 
    98 double j0(double x) 
    99 { 
     98double j0(double x) { 
     99 
     100//Cephes single precission 
     101if ( FLOAT_SIZE>4 ) { 
    100102    double w, z, p, q, xn; 
    101103 
    … …  
    103105    const double SQ2OPI = 7.9788456080286535587989E-1; 
    104106    const double PIO4 = 7.85398163397448309616E-1; 
     107 
     108    const double DR1 = 5.78318596294678452118E0; 
     109    const double DR2 = 3.04712623436620863991E1; 
    105110 
    106111    const double PP[8] = { 
    … …  
    193198    }; 
    194199 
    195  
    196     /*  5.783185962946784521175995758455807035071 */ 
    197     const double DR1 = 5.78318596294678452118E0; 
    198     /* 30.47126234366208639907816317502275584842 */ 
    199     const double DR2 = 3.04712623436620863991E1; 
    200  
    201  
    202200    if( x < 0 ) 
    203201            x = -x; 
    … …  
    225223    return( p * SQ2OPI / sqrt(x) ); 
    226224} 
    227  
    228  
     225//Cephes single precission 
     226else { 
     227    double xx, w, z, p, q, xn; 
     228 
     229    const double YZ1 =  0.43221455686510834878; 
     230    const double YZ2 = 22.401876406482861405; 
     231    const double YZ3 = 64.130620282338755553; 
     232    const double DR1 =  5.78318596294678452118; 
     233    const double PIO4F = 0.7853981633974483096; 
     234 
     235    double MO[8] = { 
     236        -6.838999669318810E-002, 
     237        1.864949361379502E-001, 
     238        -2.145007480346739E-001, 
     239        1.197549369473540E-001, 
     240        -3.560281861530129E-003, 
     241        -4.969382655296620E-002, 
     242        -3.355424622293709E-006, 
     243        7.978845717621440E-001 
     244    }; 
     245 
     246    double PH[8] = { 
     247        3.242077816988247E+001, 
     248        -3.630592630518434E+001, 
     249        1.756221482109099E+001, 
     250        -4.974978466280903E+000, 
     251        1.001973420681837E+000, 
     252        -1.939906941791308E-001, 
     253        6.490598792654666E-002, 
     254        -1.249992184872738E-001 
     255    }; 
     256 
     257    double YP[8] = { 
     258        9.454583683980369E-008, 
     259        -9.413212653797057E-006, 
     260        5.344486707214273E-004, 
     261        -1.584289289821316E-002, 
     262        1.707584643733568E-001, 
     263        0.0, 
     264        0.0, 
     265        0.0 
     266    }; 
     267 
     268    double JP[8] = { 
     269        -6.068350350393235E-008, 
     270        6.388945720783375E-006, 
     271        -3.969646342510940E-004, 
     272        1.332913422519003E-002, 
     273        -1.729150680240724E-001, 
     274        0.0, 
     275        0.0, 
     276        0.0 
     277    }; 
     278 
     279    if( x < 0 ) 
     280            xx = -x; 
     281    else 
     282            xx = x; 
     283 
     284    if( x <= 2.0 ) { 
     285            z = xx * xx; 
     286            if( x < 1.0e-3 ) 
     287                    return( 1.0 - 0.25*z ); 
     288 
     289            p = (z-DR1) * polevl( z, JP, 4); 
     290            return( p ); 
     291        } 
     292 
     293    q = 1.0/x; 
     294    w = sqrtf(q); 
     295 
     296    p = w * polevl( q, MO, 7); 
     297    w = q*q; 
     298    xn = q * polevl( w, PH, 7) - PIO4F; 
     299    p = p * cosf(xn + xx); 
     300    return(p); 
     301} 
     302} 
     303 
     304 

  • sasmodels/models/lib/j1d.c

    r3936ad3 r0a9d219  
    7979 
    8080double j1(double x) { 
     81 
     82//Cephes double pression function 
     83if (FLOAT_SIZE>4) { 
    8184 
    8285    double w, z, p, q, xn; 
    … …  
    156159    }; 
    157160 
    158     //const double Z1 = 1.46819706421238932572E1; 
    159     //const double Z2 = 4.92184563216946036703E1; 
    160  
    161161    w = x; 
    162162    if( x < 0 ) 
    … …  
    182182 
    183183    return( p * SQ2OPI / sqrt(x) ); 
     184 
    184185} 
     186//Single precission version of cephes 
     187else { 
     188    double xx, w, z, p, q, xn; 
     189 
     190    const double Z1 = 1.46819706421238932572E1; 
     191    const double THPIO4F =  2.35619449019234492885;    /* 3*pi/4 */ 
     192 
     193 
     194    double JP[8] = { 
     195        -4.878788132172128E-009, 
     196        6.009061827883699E-007, 
     197        -4.541343896997497E-005, 
     198        1.937383947804541E-003, 
     199        -3.405537384615824E-002, 
     200        0.0, 
     201        0.0, 
     202        0.0 
     203    }; 
     204 
     205    double MO1[8] = { 
     206        6.913942741265801E-002, 
     207        -2.284801500053359E-001, 
     208        3.138238455499697E-001, 
     209        -2.102302420403875E-001, 
     210        5.435364690523026E-003, 
     211        1.493389585089498E-001, 
     212        4.976029650847191E-006, 
     213        7.978845453073848E-001 
     214    }; 
     215 
     216    double PH1[8] = { 
     217        -4.497014141919556E+001, 
     218        5.073465654089319E+001, 
     219        -2.485774108720340E+001, 
     220        7.222973196770240E+000, 
     221        -1.544842782180211E+000, 
     222        3.503787691653334E-001, 
     223        -1.637986776941202E-001, 
     224        3.749989509080821E-001 
     225    }; 
     226 
     227    xx = x; 
     228    if( xx < 0 ) 
     229            xx = -x; 
     230 
     231    if( xx <= 2.0 ) 
     232        { 
     233            z = xx * xx; 
     234            p = (z-Z1) * xx * polevl( z, JP, 4 ); 
     235            return( p ); 
     236        } 
     237 
     238    q = 1.0/x; 
     239    w = sqrt(q); 
     240 
     241    p = w * polevl( q, MO1, 7); 
     242    w = q*q; 
     243    xn = q * polevl( w, PH1, 7) - THPIO4F; 
     244    p = p * cos(xn + xx); 
     245 
     246    return(p); 
     247} 
     248} 
     249 

  • sasmodels/models/lib/jnd.c

    r3936ad3 r0a9d219  
    4949 
    5050double jn( int n, double x ); 
    51 #define MACHEP 1.11022302462515654042E-16 
    5251 
    53 double jn( int n, double x ) 
    54 { 
    55 double pkm2, pkm1, pk, xk, r, ans; 
    56 int k, sign; 
     52double jn( int n, double x ) { 
    5753 
    58 if( n < 0 ) 
    59         { 
    60         n = -n; 
    61         if( (n & 1) == 0 )      /* -1**n */ 
    62                 sign = 1; 
    63         else 
    64                 sign = -1; 
     54    const double MACHEP = 1.11022302462515654042E-16; 
     55    double pkm2, pkm1, pk, xk, r, ans, xinv; 
     56    int k, sign; 
     57 
     58    if( n < 0 ) { 
     59            n = -n; 
     60            if( (n & 1) == 0 )  /* -1**n */ 
     61                    sign = 1; 
     62            else 
     63                    sign = -1; 
    6564        } 
    66 else 
    67         sign = 1; 
     65    else 
     66            sign = 1; 
    6867 
    69 if( x < 0.0 ) 
    70         { 
    71         if( n & 1 ) 
    72                 sign = -sign; 
    73         x = -x; 
     68    if( x < 0.0 ) { 
     69            if( n & 1 ) 
     70                    sign = -sign; 
     71            x = -x; 
    7472        } 
    7573 
    76 if( n == 0 ) 
    77         return( sign * j0(x) ); 
    78 if( n == 1 ) 
    79         return( sign * j1(x) ); 
    80 if( n == 2 ) 
    81         return( sign * (2.0 * j1(x) / x  -  j0(x)) ); 
     74    if( n == 0 ) 
     75            return( sign * j0(x) ); 
     76    if( n == 1 ) 
     77            return( sign * j1(x) ); 
     78    if( n == 2 ) 
     79            return( sign * (2.0 * j1(x) / x  -  j0(x)) ); 
    8280 
    83 if( x < MACHEP ) 
    84         return( 0.0 ); 
     81    if( x < MACHEP ) 
     82            return( 0.0 ); 
    8583 
     84    if (FLOAT_SIZE > 4) 
     85        k = 53; 
     86    else 
     87        k = 24; 
    8688 
    87 pk = 2 * (n + k); 
    88 ans = pk; 
    89 xk = x * x; 
     89    pk = 2 * (n + k); 
     90    ans = pk; 
     91    xk = x * x; 
    9092 
    91 do 
    92         { 
    93         pk -= 2.0; 
    94         ans = pk - (xk/ans); 
     93    do { 
     94            pk -= 2.0; 
     95            ans = pk - (xk/ans); 
     96        } while( --k > 0 ); 
     97 
     98    /* backward recurrence */ 
     99 
     100    pk = 1.0; 
     101 
     102    if (FLOAT_SIZE > 4) { 
     103        ans = x/ans; 
     104        pkm1 = 1.0/ans; 
     105 
     106        k = n-1; 
     107        r = 2 * k; 
     108 
     109        do { 
     110            pkm2 = (pkm1 * r  -  pk * x) / x; 
     111                pk = pkm1; 
     112                pkm1 = pkm2; 
     113                r -= 2.0; 
     114            } while( --k > 0 ); 
     115 
     116        if( fabs(pk) > fabs(pkm1) ) 
     117                ans = j1(x)/pk; 
     118        else 
     119                ans = j0(x)/pkm1; 
     120 
     121            return( sign * ans ); 
    95122        } 
    96 while( --k > 0 ); 
    97 ans = x/ans; 
     123    else { 
     124        xinv = 1.0/x; 
     125        pkm1 = ans * xinv; 
     126        k = n-1; 
     127        r = (float )(2 * k); 
    98128 
    99 /* backward recurrence */ 
     129        do { 
     130                pkm2 = (pkm1 * r  -  pk * x) * xinv; 
     131                pk = pkm1; 
     132                pkm1 = pkm2; 
     133                r -= 2.0; 
     134            } 
     135        while( --k > 0 ); 
    100136 
    101 pk = 1.0; 
    102 pkm1 = 1.0/ans; 
    103 k = n-1; 
    104 r = 2 * k; 
     137        r = pk; 
     138        if( r < 0 ) 
     139                r = -r; 
     140        ans = pkm1; 
     141        if( ans < 0 ) 
     142                ans = -ans; 
    105143 
    106 do 
    107         { 
    108         pkm2 = (pkm1 * r  -  pk * x) / x; 
    109         pk = pkm1; 
    110         pkm1 = pkm2; 
    111         r -= 2.0; 
    112         } 
    113 while( --k > 0 ); 
    114  
    115 if( fabs(pk) > fabs(pkm1) ) 
    116         ans = j1(x)/pk; 
    117 else 
    118         ans = j0(x)/pkm1; 
    119 return( sign * ans ); 
     144        if( r > ans )  /* if( fabs(pk) > fabs(pkm1) ) */ 
     145                ans = sign * j1(x)/pk; 
     146        else 
     147                ans = sign * j0(x)/pkm1; 
     148        return( ans ); 
     149    } 
    120150} 
    121151 
Note: See TracChangeset for help on using the changeset viewer.