Changeset 3f8584a2 in sasmodels for sasmodels/models/lib/sas_JN.c

Timestamp:

Jul 27, 2016 3: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

File:

• sasmodels/models/lib/sas_JN.c (modified) (4 diffs)

sasmodels/models/lib/sas_JN.c

@@ -48 +48 @@
 */
 
-double sas_JN( int n, double x );
-
-double sas_JN( int n, double x ) {
-
+#if FLOAT_SIZE > 4
+
+double jn( int n, double x );
+double jn( int n, double x ) {
+
+    // PAK: seems to be machine epsilon/2
     const double MACHEP = 1.11022302462515654042E-16;
     double pkm2, pkm1, pk, xk, r, ans;
@@ -57 +59 @@
 
     if( n < 0 ) {
-            n = -n;
-            if( (n & 1) == 0 )  /* -1**n */
-                    sign = 1;
-            else
-                    sign = -1;
-        }
-    else
-            sign = 1;
+        n = -n;
+        if( (n & 1) == 0 )      /* -1**n */
+            sign = 1;
+        else
+            sign = -1;
+    } else {
+        sign = 1;
+    }
 
     if( x < 0.0 ) {
-            if( n & 1 )
-                    sign = -sign;
-            x = -x;
-        }
+        if( n & 1 )
+            sign = -sign;
+    x = -x;
+    }
 
     if( n == 0 )
-            return( sign * sas_J0(x) );
+        return( sign * j0(x) );
     if( n == 1 )
-            return( sign * sas_J1(x) );
+        return( sign * j1(x) );
     if( n == 2 )
-            return( sign * (2.0 * sas_J1(x) / x  -  sas_J0(x)) );
+        return( sign * (2.0 * j1(x) / x  -  j0(x)) );
 
     if( x < MACHEP )
-            return( 0.0 );
-
-    #if FLOAT_SIZE > 4
-        k = 53;
-    #else
-        k = 24;
-    #endif
-
+        return( 0.0 );
+
+    k = 53;
     pk = 2 * (n + k);
     ans = pk;
@@ -93 +90 @@
 
     do {
-            pk -= 2.0;
-            ans = pk - (xk/ans);
-        } while( --k > 0 );
+        pk -= 2.0;
+        ans = pk - (xk/ans);
+    } while( --k > 0 );
 
     /* backward recurrence */
@@ -101 +98 @@
     pk = 1.0;
 
-    #if FLOAT_SIZE > 4
-        ans = x/ans;
-        pkm1 = 1.0/ans;
-
-        k = n-1;
-        r = 2 * k;
-
-        do {
-            pkm2 = (pkm1 * r  -  pk * x) / x;
-                pk = pkm1;
-                pkm1 = pkm2;
-                r -= 2.0;
-            } while( --k > 0 );
-
-        if( fabs(pk) > fabs(pkm1) )
-                ans = sas_J1(x)/pk;
+    ans = x/ans;
+    pkm1 = 1.0/ans;
+
+    k = n-1;
+    r = 2 * k;
+
+    do {
+        pkm2 = (pkm1 * r  -  pk * x) / x;
+        pk = pkm1;
+        pkm1 = pkm2;
+        r -= 2.0;
+    } while( --k > 0 );
+
+    if( fabs(pk) > fabs(pkm1) )
+        ans = j1(x)/pk;
+    else
+        ans = j0(x)/pkm1;
+
+    return( sign * ans );
+}
+
+#else
+
+float jnf(int n, float x)
+{
+    // PAK: seems to be machine epsilon/2
+    const double MACHEP = 5.9604645e-08;
+    float pkm2, pkm1, pk, xk, r, ans;
+    int k, sign;
+
+    if( n < 0 ) {
+        n = -n;
+        if( (n & 1) == 0 ) /* -1**n */
+            sign = 1;
         else
-                ans = sas_J0(x)/pkm1;
-
-            return( sign * ans );
-
-    #else
-        const double xinv = 1.0/x;
-        pkm1 = ans * xinv;
-        k = n-1;
-        r = (float )(2 * k);
-
-        do {
-                pkm2 = (pkm1 * r  -  pk * x) * xinv;
-                pk = pkm1;
-                pkm1 = pkm2;
-                r -= 2.0;
-            }
-        while( --k > 0 );
-
-        r = pk;
-        if( r < 0 )
-                r = -r;
-        ans = pkm1;
-        if( ans < 0 )
-                ans = -ans;
-
-        if( r > ans )  /* if( fabs(pk) > fabs(pkm1) ) */
-                ans = sign * sas_J1(x)/pk;
-        else
-                ans = sign * sas_J0(x)/pkm1;
-        return( ans );
-    #endif
+            sign = -1;
+    } else {
+        sign = 1;
+    }
+
+    if( x < 0.0 ) {
+        if( n & 1 )
+            sign = -sign;
+        x = -x;
+    }
+
+    if( n == 0 )
+        return( sign * j0f(x) );
+    if( n == 1 )
+        return( sign * j1f(x) );
+    if( n == 2 )
+        return( sign * (2.0 * j1f(x) / x  -  j0f(x)) );
+
+    if( x < MACHEP )
+        return( 0.0 );
+
+    k = 24;
+    pk = 2 * (n + k);
+    ans = pk;
+    xk = x * x;
+
+    do {
+        pk -= 2.0;
+        ans = pk - (xk/ans);
+    } while( --k > 0 );
+
+    /* backward recurrence */
+
+    pk = 1.0;
+
+    const float xinv = 1.0/x;
+    pkm1 = ans * xinv;
+    k = n-1;
+    r = (float )(2 * k);
+
+    do {
+        pkm2 = (pkm1 * r  -  pk * x) * xinv;
+        pk = pkm1;
+        pkm1 = pkm2;
+        r -= 2.0;
+    } while( --k > 0 );
+
+    r = pk;
+    if( r < 0 )
+        r = -r;
+    ans = pkm1;
+    if( ans < 0 )
+        ans = -ans;
+
+    if( r > ans )  /* if( fabs(pk) > fabs(pkm1) ) */
+        ans = sign * j1f(x)/pk;
+    else
+        ans = sign * j0f(x)/pkm1;
+    return( ans );
 }
-
+#endif
+
+#if FLOAT_SIZE>4
+#define sas_JN jn
+#else
+#define sas_JN jnf
+#endif