]> www.ginac.de Git - cln.git/blobdiff - examples/atan_recip.cc
Revert "Convert complex numbers to real numbers if imaginary part is floating-point...
[cln.git] / examples / atan_recip.cc
index 1540bd5761fbe6aa9d6ce73705f60e9969c4e306..f061fae2fff916c01b72573672cc0cef69da876f 100644 (file)
@@ -31,11 +31,11 @@ using namespace cln;
 const cl_LF atan_recip_1a (cl_I m, uintC len)
 {
        var uintC actuallen = len + 1;
-       var cl_LF eps = scale_float(cl_I_to_LF(1,actuallen),-intDsize*(sintL)actuallen);
+       var cl_LF eps = scale_float(cl_I_to_LF(1,actuallen),-intDsize*(sintC)actuallen);
        var cl_I m2 = m*m;
        var cl_LF fterm = cl_I_to_LF(1,actuallen)/m;
        var cl_LF fsum = fterm;
-       for (var uintL n = 1; fterm >= eps; n++) {
+       for (var uintC n = 1; fterm >= eps; n++) {
                fterm = fterm/m2;
                fterm = cl_LF_shortenwith(fterm,eps);
                if ((n % 2) == 0)
@@ -52,23 +52,23 @@ const cl_LF atan_recip_1b (cl_I m, uintC len)
        var cl_I m2 = m*m;
        var cl_I fterm = floor1((cl_I)1 << (intDsize*actuallen), m);
        var cl_I fsum = fterm;
-       for (var uintL n = 1; fterm > 0; n++) {
+       for (var uintC n = 1; fterm > 0; n++) {
                fterm = floor1(fterm,m2);
                if ((n % 2) == 0)
                        fsum = fsum + floor1(fterm,2*n+1);
                else
                        fsum = fsum - floor1(fterm,2*n+1);
        }
-       return scale_float(cl_I_to_LF(fsum,len),-intDsize*(sintL)actuallen);
+       return scale_float(cl_I_to_LF(fsum,len),-intDsize*(sintC)actuallen);
 }
 
 const cl_LF atan_recip_1c (cl_I m, uintC len)
 {
        var uintC actuallen = len + 1;
        var cl_I m2 = m*m;
-       var sintL N = (sintL)(0.69314718*intDsize/2*actuallen/log(double_approx(m))) + 1;
+       var sintC N = (sintC)(0.69314718*intDsize/2*actuallen/log(double_approx(m))) + 1;
        var cl_I num = 0, den = 1; // "lazy rational number"
-       for (sintL n = N-1; n>=0; n--) {
+       for (sintC n = N-1; n>=0; n--) {
                // Multiply sum with 1/m^2:
                den = den * m2;
                // Add (-1)^n/(2n+1):
@@ -87,11 +87,11 @@ const cl_LF atan_recip_1d (cl_I m, uintC len)
 {
        var uintC actuallen = len + 1;
        var cl_I m2 = m*m;
-       var uintL N = (uintL)(0.69314718*intDsize/2*actuallen/log(double_approx(m))) + 1;
+       var uintC N = (uintC)(0.69314718*intDsize/2*actuallen/log(double_approx(m))) + 1;
        CL_ALLOCA_STACK;
        var cl_I* bv = (cl_I*) cl_alloca(N*sizeof(cl_I));
        var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
-       var uintL n;
+       var uintC n;
        for (n = 0; n < N; n++) {
                new (&bv[n]) cl_I ((n % 2) == 0 ? (cl_I)(2*n+1) : -(cl_I)(2*n+1));
                new (&qv[n]) cl_I (n==0 ? m : m2);
@@ -113,11 +113,11 @@ const cl_LF atan_recip_1d (cl_I m, uintC len)
 const cl_LF atan_recip_2a (cl_I m, uintC len)
 {
        var uintC actuallen = len + 1;
-       var cl_LF eps = scale_float(cl_I_to_LF(1,actuallen),-intDsize*(sintL)actuallen);
+       var cl_LF eps = scale_float(cl_I_to_LF(1,actuallen),-intDsize*(sintC)actuallen);
        var cl_I m2 = m*m+1;
        var cl_LF fterm = cl_I_to_LF(m,actuallen)/m2;
        var cl_LF fsum = fterm;
-       for (var uintL n = 1; fterm >= eps; n++) {
+       for (var uintC n = 1; fterm >= eps; n++) {
                fterm = The(cl_LF)((2*n)*fterm)/((2*n+1)*m2);
                fterm = cl_LF_shortenwith(fterm,eps);
                fsum = fsum + LF_to_LF(fterm,actuallen);
@@ -131,20 +131,20 @@ const cl_LF atan_recip_2b (cl_I m, uintC len)
        var cl_I m2 = m*m+1;
        var cl_I fterm = floor1((cl_I)m << (intDsize*actuallen), m2);
        var cl_I fsum = fterm;
-       for (var uintL n = 1; fterm > 0; n++) {
+       for (var uintC n = 1; fterm > 0; n++) {
                fterm = floor1((2*n)*fterm,(2*n+1)*m2);
                fsum = fsum + fterm;
        }
-       return scale_float(cl_I_to_LF(fsum,len),-intDsize*(sintL)actuallen);
+       return scale_float(cl_I_to_LF(fsum,len),-intDsize*(sintC)actuallen);
 }
 
 const cl_LF atan_recip_2c (cl_I m, uintC len)
 {
        var uintC actuallen = len + 1;
        var cl_I m2 = m*m+1;
-       var uintL N = (uintL)(0.69314718*intDsize*actuallen/log(double_approx(m2))) + 1;
+       var uintC N = (uintC)(0.69314718*intDsize*actuallen/log(double_approx(m2))) + 1;
        var cl_I num = 0, den = 1; // "lazy rational number"
-       for (uintL n = N; n>0; n--) {
+       for (uintC n = N; n>0; n--) {
                // Multiply sum with (2n)/(2n+1)(m^2+1):
                num = num * (2*n);
                den = den * ((2*n+1)*m2);
@@ -161,11 +161,11 @@ const cl_LF atan_recip_2d (cl_I m, uintC len)
 {
        var uintC actuallen = len + 1;
        var cl_I m2 = m*m+1;
-       var uintL N = (uintL)(0.69314718*intDsize*actuallen/log(double_approx(m2))) + 1;
+       var uintC N = (uintC)(0.69314718*intDsize*actuallen/log(double_approx(m2))) + 1;
        CL_ALLOCA_STACK;
        var cl_I* pv = (cl_I*) cl_alloca(N*sizeof(cl_I));
        var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
-       var uintL n;
+       var uintC n;
        new (&pv[0]) cl_I (m);
        new (&qv[0]) cl_I (m2);
        for (n = 1; n < N; n++) {
@@ -196,7 +196,7 @@ int main (int argc, char * argv[])
        if (argc < 2)
                exit(1);
        cl_I m = (cl_I)argv[1];
-       uintL len = atoi(argv[2]);
+       uintC len = atol(argv[2]);
        cl_LF p;
        ln(cl_I_to_LF(1000,len+10)); // fill cache
        // Method 1.
@@ -251,7 +251,7 @@ int main (int argc, char * argv[])
 
 
 // Timings of the above algorithms, on an i486 33 MHz, running Linux.
-// m = 390112. (For Jörg Arndt's formula (4.15).)
+// m = 390112. (For Jörg Arndt's formula (4.15).)
 //    N      1a     1b     1c     1d      2a     2b     2c     2d      3
 //    10     0.0027 0.0018 0.0019 0.0019  0.0032 0.0022 0.0019 0.0019  0.0042
 //    25     0.0085 0.0061 0.0058 0.0061  0.0095 0.0069 0.0056 0.0061  0.028
@@ -265,10 +265,10 @@ int main (int argc, char * argv[])
 // 10000
 // asymp.    N^2    N^2    N^2    FAST    N^2    N^2    N^2    FAST    FAST
 //
-// m = 319. (For Jörg Arndt's formula (4.7).)
+// m = 319. (For Jörg Arndt's formula (4.7).)
 //    N      1a     1b     1c     1d      2a     2b     2c     2d      3
 //  1000     6.06   4.40   9.17   3.82    5.29   3.90   7.50   3.53   50.3
 //
-// m = 18. (For Jörg Arndt's formula (4.4).)
+// m = 18. (For Jörg Arndt's formula (4.4).)
 //    N      1a     1b     1c     1d      2a     2b     2c     2d      3
 //  1000    11.8    9.0   22.3    6.0    10.2    7.7   17.1    5.7    54.3