examples/atanh_recip.cc

   1 // Computation of artanh(1/m) (m integer) to high precision.
   2
   3 #include "cln/integer.h"
   4 #include "cln/rational.h"
   5 #include "cln/real.h"
   6 #include "cln/complex.h"
   7 #include "cln/lfloat.h"
   8 #include "cl_LF.h"
   9 #include "cl_LF_tran.h"
  10 #include "cl_alloca.h"
  11 #include <cstdlib>
  12 #include <cstring>
  13 #include "cln/timing.h"
  14
  15 #undef floor
  16 #include <cmath>
  17 #define floor cln_floor
  18
  19
  20 // Method 1: atanh(1/m) = sum(n=0..infty, 1/(2n+1) * 1/m^(2n+1))
  21 // Method 2: atanh(1/m) = sum(n=0..infty, (-4)^n*n!^2/(2n+1)! * m/(m^2-1)^(n+1))
  22 // a. Using long floats.                     [N^2]
  23 // b. Simulating long floats using integers. [N^2]
  24 // c. Using integers, no binary splitting.   [N^2]
  25 // d. Using integers, with binary splitting. [FAST]
  26 // Method 3: general built-in algorithm.     [FAST]
  27 // Method 4: atanh(x) = 1/2 ln((1+x)/(1-x)),
  28 // using the general built-in algorithm      [FAST]
  29
  30
  31 // Method 1: atanh(1/m) = sum(n=0..infty, 1/(2n+1) * 1/m^(2n+1))
  32
  33 const cl_LF atanh_recip_1a (cl_I m, uintC len)
  34 {
  35         var uintC actuallen = len + 1;
  36         var cl_LF eps = scale_float(cl_I_to_LF(1,actuallen),-intDsize*(sintL)actuallen);
  37         var cl_I m2 = m*m;
  38         var cl_LF fterm = cl_I_to_LF(1,actuallen)/m;
  39         var cl_LF fsum = fterm;
  40         for (var uintL n = 1; fterm >= eps; n++) {
  41                 fterm = fterm/m2;
  42                 fterm = cl_LF_shortenwith(fterm,eps);
  43                 fsum = fsum + LF_to_LF(fterm/(2*n+1),actuallen);
  44         }
  45         return shorten(fsum,len);
  46 }
  47
  48 const cl_LF atanh_recip_1b (cl_I m, uintC len)
  49 {
  50         var uintC actuallen = len + 1;
  51         var cl_I m2 = m*m;
  52         var cl_I fterm = floor1((cl_I)1 << (intDsize*actuallen), m);
  53         var cl_I fsum = fterm;
  54         for (var uintL n = 1; fterm > 0; n++) {
  55                 fterm = floor1(fterm,m2);
  56                 fsum = fsum + floor1(fterm,2*n+1);
  57         }
  58         return scale_float(cl_I_to_LF(fsum,len),-intDsize*(sintL)actuallen);
  59 }
  60
  61 const cl_LF atanh_recip_1c (cl_I m, uintC len)
  62 {
  63         var uintC actuallen = len + 1;
  64         var cl_I m2 = m*m;
  65         var sintL N = (sintL)(0.69314718*intDsize/2*actuallen/log(double_approx(m))) + 1;
  66         var cl_I num = 0, den = 1; // "lazy rational number"
  67         for (sintL n = N-1; n>=0; n--) {
  68                 // Multiply sum with 1/m^2:
  69                 den = den * m2;
  70                 // Add 1/(2n+1):
  71                 num = num*(2*n+1) + den;
  72                 den = den*(2*n+1);
  73         }
  74         den = den*m;
  75         var cl_LF result = cl_I_to_LF(num,actuallen)/cl_I_to_LF(den,actuallen);
  76         return shorten(result,len);
  77 }
  78
  79 const cl_LF atanh_recip_1d (cl_I m, uintC len)
  80 {
  81         var uintC actuallen = len + 1;
  82         var cl_I m2 = m*m;
  83         var uintL N = (uintL)(0.69314718*intDsize/2*actuallen/log(double_approx(m))) + 1;
  84         CL_ALLOCA_STACK;
  85         var cl_I* bv = (cl_I*) cl_alloca(N*sizeof(cl_I));
  86         var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
  87         var uintL n;
  88         for (n = 0; n < N; n++) {
  89                 new (&bv[n]) cl_I ((cl_I)(2*n+1));
  90                 new (&qv[n]) cl_I (n==0 ? m : m2);
  91         }
  92         var cl_rational_series series;
  93         series.av = NULL; series.bv = bv;
  94         series.pv = NULL; series.qv = qv; series.qsv = NULL;
  95         var cl_LF result = eval_rational_series(N,series,actuallen);
  96         for (n = 0; n < N; n++) {
  97                 bv[n].~cl_I();
  98                 qv[n].~cl_I();
  99         }
 100         return shorten(result,len);
 101 }
 102
 103
 104 // Method 2: atanh(1/m) = sum(n=0..infty, (-4)^n*n!^2/(2n+1)! * m/(m^2-1)^(n+1))
 105
 106 const cl_LF atanh_recip_2a (cl_I m, uintC len)
 107 {
 108         var uintC actuallen = len + 1;
 109         var cl_LF eps = scale_float(cl_I_to_LF(1,actuallen),-intDsize*(sintL)actuallen);
 110         var cl_I m2 = m*m-1;
 111         var cl_LF fterm = cl_I_to_LF(m,actuallen)/m2;
 112         var cl_LF fsum = fterm;
 113         for (var uintL n = 1; fterm >= eps; n++) {
 114                 fterm = The(cl_LF)((2*n)*fterm)/((2*n+1)*m2);
 115                 fterm = cl_LF_shortenwith(fterm,eps);
 116                 if ((n % 2) == 0)
 117                         fsum = fsum + LF_to_LF(fterm,actuallen);
 118                 else
 119                         fsum = fsum - LF_to_LF(fterm,actuallen);
 120         }
 121         return shorten(fsum,len);
 122 }
 123
 124 const cl_LF atanh_recip_2b (cl_I m, uintC len)
 125 {
 126         var uintC actuallen = len + 1;
 127         var cl_I m2 = m*m-1;
 128         var cl_I fterm = floor1((cl_I)m << (intDsize*actuallen), m2);
 129         var cl_I fsum = fterm;
 130         for (var uintL n = 1; fterm > 0; n++) {
 131                 fterm = floor1((2*n)*fterm,(2*n+1)*m2);
 132                 if ((n % 2) == 0)
 133                         fsum = fsum + fterm;
 134                 else
 135                         fsum = fsum - fterm;
 136         }
 137         return scale_float(cl_I_to_LF(fsum,len),-intDsize*(sintL)actuallen);
 138 }
 139
 140 const cl_LF atanh_recip_2c (cl_I m, uintC len)
 141 {
 142         var uintC actuallen = len + 1;
 143         var cl_I m2 = m*m-1;
 144         var uintL N = (uintL)(0.69314718*intDsize*actuallen/log(double_approx(m2))) + 1;
 145         var cl_I num = 0, den = 1; // "lazy rational number"
 146         for (uintL n = N; n>0; n--) {
 147                 // Multiply sum with -(2n)/(2n+1)(m^2+1):
 148                 num = num * (2*n);
 149                 den = - den * ((2*n+1)*m2);
 150                 // Add 1:
 151                 num = num + den;
 152         }
 153         num = num*m;
 154         den = den*m2;
 155         var cl_LF result = cl_I_to_LF(num,actuallen)/cl_I_to_LF(den,actuallen);
 156         return shorten(result,len);
 157 }
 158
 159 const cl_LF atanh_recip_2d (cl_I m, uintC len)
 160 {
 161         var uintC actuallen = len + 1;
 162         var cl_I m2 = m*m-1;
 163         var uintL N = (uintL)(0.69314718*intDsize*actuallen/log(double_approx(m2))) + 1;
 164         CL_ALLOCA_STACK;
 165         var cl_I* pv = (cl_I*) cl_alloca(N*sizeof(cl_I));
 166         var cl_I* qv = (cl_I*) cl_alloca(N*sizeof(cl_I));
 167         var uintL n;
 168         new (&pv[0]) cl_I (m);
 169         new (&qv[0]) cl_I (m2);
 170         for (n = 1; n < N; n++) {
 171                 new (&pv[n]) cl_I (-(cl_I)(2*n));
 172                 new (&qv[n]) cl_I ((2*n+1)*m2);
 173         }
 174         var cl_rational_series series;
 175         series.av = NULL; series.bv = NULL;
 176         series.pv = pv; series.qv = qv; series.qsv = NULL;
 177         var cl_LF result = eval_rational_series(N,series,actuallen);
 178         for (n = 0; n < N; n++) {
 179                 pv[n].~cl_I();
 180                 qv[n].~cl_I();
 181         }
 182         return shorten(result,len);
 183 }
 184
 185
 186 // Main program: Compute and display the timings.
 187
 188 int main (int argc, char * argv[])
 189 {
 190         int repetitions = 1;
 191         if ((argc >= 3) && !strcmp(argv[1],"-r")) {
 192                 repetitions = atoi(argv[2]);
 193                 argc -= 2; argv += 2;
 194         }
 195         if (argc < 2)
 196                 exit(1);
 197         cl_I m = (cl_I)argv[1];
 198         uintL len = atoi(argv[2]);
 199         cl_LF p;
 200         ln(cl_I_to_LF(1000,len+10)); // fill cache
 201         // Method 1.
 202         { CL_TIMING;
 203           for (int rep = repetitions; rep > 0; rep--)
 204             { p = atanh_recip_1a(m,len); }
 205         }
 206         cout << p << endl;
 207         { CL_TIMING;
 208           for (int rep = repetitions; rep > 0; rep--)
 209             { p = atanh_recip_1b(m,len); }
 210         }
 211         cout << p << endl;
 212         { CL_TIMING;
 213           for (int rep = repetitions; rep > 0; rep--)
 214             { p = atanh_recip_1c(m,len); }
 215         }
 216         cout << p << endl;
 217         { CL_TIMING;
 218           for (int rep = repetitions; rep > 0; rep--)
 219             { p = atanh_recip_1d(m,len); }
 220         }
 221         cout << p << endl;
 222         // Method 2.
 223         { CL_TIMING;
 224           for (int rep = repetitions; rep > 0; rep--)
 225             { p = atanh_recip_2a(m,len); }
 226         }
 227         cout << p << endl;
 228         { CL_TIMING;
 229           for (int rep = repetitions; rep > 0; rep--)
 230             { p = atanh_recip_2b(m,len); }
 231         }
 232         cout << p << endl;
 233         { CL_TIMING;
 234           for (int rep = repetitions; rep > 0; rep--)
 235             { p = atanh_recip_2c(m,len); }
 236         }
 237         cout << p << endl;
 238         { CL_TIMING;
 239           for (int rep = repetitions; rep > 0; rep--)
 240             { p = atanh_recip_2d(m,len); }
 241         }
 242         cout << p << endl;
 243         // Method 3.
 244         { CL_TIMING;
 245           for (int rep = repetitions; rep > 0; rep--)
 246             { p = The(cl_LF)(atanh(cl_RA_to_LF(1/(cl_RA)m,len))); }
 247         }
 248         cout << p << endl;
 249         // Method 4.
 250         { CL_TIMING;
 251           for (int rep = repetitions; rep > 0; rep--)
 252             { p = The(cl_LF)(scale_float(ln(cl_RA_to_LF((cl_RA)(m+1)/(cl_RA)(m-1),len)),-1)); }
 253         }
 254         cout << p << endl;
 255 }
 256
 257
 258 // Timings of the above algorithms, on an i486 33 MHz, running Linux.
 259 // m = 3 -> 1/2 ln(2)
 260 //   N    1a     1b     1c     1d     2a     2b     2c     2d      3
 261 //   10   0.021  0.014  0.019  0.012  0.029  0.015  0.023  0.015   0.015
 262 //   25   0.060  0.041  0.073  0.041  0.082  0.046  0.086  0.051   0.066
 263 //   50   0.164  0.110  0.258  0.120  0.203  0.124  0.295  0.142
 264 //  100   0.49   0.35   1.05   0.37   0.60   0.35   1.19   0.42
 265 //  250   2.5    1.9    7.2    1.7    2.9    1.9    8.0    1.8
 266 //  500  10.1    7.2   33.4    5.5   10.7    7.3   36.5    5.9
 267 // 1000  38     30    145     16.1   39     29    158     16.8
 268 // 2500 231    188    976     53    237    186   1081     58
 269 // asymp. N^2    N^2    N^2    FAST   N^2    N^2    N^2    FAST
 270 //
 271 // m = 9 -> 1/2 ln(5/4)
 272 //   N    1a     1b     1c     1d     2a     2b     2c     2d     3      4
 273 //   10   0.0106 0.0072 0.0084 0.0061 0.0139 0.0073 0.0098 0.0073 0.0140 0.0211
 274 //   25   0.031  0.021  0.029  0.019  0.039  0.022  0.031  0.022  0.063  0.081
 275 //   50   0.083  0.057  0.091  0.056  0.098  0.058  0.098  0.060  0.232  0.212
 276 //  100   0.25   0.17   0.32   0.16   0.28   0.17   0.34   0.17   0.60   0.59
 277 //  250   1.28   0.94   2.11   0.77   1.40   0.91   2.18   0.76   2.76   2.76
 278 //  500   5.1    3.6    9.4    2.5    5.2    3.4    9.3    2.4   10.4    9.7
 279 // 1000  19.1   14.7   42      7.8   18.5   13.6   42      7.4   31     30
 280 // 2500 116     93    279     29.6  113     86    278     30.0  129    125
 281 // asymp. N^2    N^2    N^2    FAST   N^2    N^2    N^2    FAST   FAST   FAST