src/base/digitseq/cl_2DS_div.cc

   1 // div2adic().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_2DS.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_2D.h"
  13 #include "cl_DS.h"
  14 #include "cl_abort.h"
  15
  16 // Time for dividing a n word number by a n word number, this is the common
  17 // case and therefore the important one:
  18 // OS: Linux 2.2, intDsize==32,        OS: TRU64/4.0, intDsize==64,
  19 // Machine: P-III/450MHz               Machine: EV5/300MHz:
  20 //      n   standard  Newton             standard  Newton
  21 //      30   0.00002   0.00006            0.00004   0.00020
  22 //     100   0.00009   0.00045            0.00033   0.0015
  23 //     300   0.00069   0.0028             0.0028    0.0085
  24 //    1000   0.018     0.019              0.031     0.065
  25 //    2000   0.028     0.057              0.12      0.20
  26 //    3000   0.078     0.11  <-(~4500)    0.28      0.23  <-(~2700)
  27 //   10000   1.09      0.48               3.14      1.13
  28 //   30000  10.1       1.21              29.7       2.70
  29 // Time for dividing a 2*n word number by a n word number:
  30 // OS: Linux 2.2, intDsize==32,        OS: TRU64/4.0, intDsize==64,
  31 // Machine: P-III/450MHz               Machine: EV5/300MHz:
  32 //      n   standard  Newton             standard  Newton
  33 //      30   0.00004   0.00019            0.00013   0.00067
  34 //     100   0.00032   0.0014             0.0013    0.0046
  35 //     300   0.0027    0.0084             0.011     0.025
  36 //    1000   0.029     0.057              0.12      0.20
  37 //    2000   0.16      0.18  <-(~2400)    0.50      0.46  <-(~1800)
  38 //    3000   0.38      0.22               1.1       0.50
  39 //   10000   4.5       1.05              13.0       2.48
  40 //   30000  51.7       2.67             120.0       6.31
  41 //          Newton faster for:         Newton faster for:
  42 // 1.0*N / N  3300<N<3800, 4400<N        2700<N<3600, N<3800
  43 // 1.1*N / N  3100<N<3700, 4100<N        2450<N<3300, N<3450
  44 // 1.2*N / N  2850<N<3200, 3700<N        2250<N
  45 // 1.3*N / N  2650<N<3000, 3450<N        2050<N
  46 // 1.4*N / N  3400<N                     1850<N
  47 // 1.5*N / N  3100<N                     1750<N
  48 // 1.6*N / N  2850<N                     1650<N
  49 // 1.7*N / N  2650<N                     1600<N
  50 // 1.8*N / N  2550<N                     1500<N
  51 // 1.9*N / N  2450<N                     1400<N
  52 // 2.0*N / N  2400<N                     1350<N
  53 //
  54 // Break-even-point. When in doubt, prefer to choose the standard algorithm.
  55 #if CL_USE_GMP
  56   static inline cl_boolean cl_recip_suitable (uintL m, uintL n) // n <= m
  57     { if (n < 2000)
  58         return cl_false;
  59       else // when n >= 4400/(m/n)^2, i.e. (m/66)^2 > n
  60         { var uintL mq = floor(m,66);
  61           if ((mq >= bit(16)) || ((uintL)(mq*mq) > n))
  62             return cl_true;
  63           else
  64             return cl_false;
  65         }
  66     }
  67 #else
  68 // Use the old default values from CLN version <= 1.0.3 as a crude estimate.
  69 // They came from timings on a i486 33 MHz running Linux:
  70 // Divide N digits by N digits          Divide 2*N digits by N digits
  71 //     N   standard  Newton                 N   standard  Newton
  72 //     10    0.00015 0.00054                10    0.00023 0.00054
  73 //     25    0.00065 0.00256                25    0.00116 0.00256
  74 //     50    0.0024  0.0083                 50    0.0044  0.0082
  75 //    100    0.0089  0.027                 100    0.0172  0.027
  76 //    250    0.054   0.130                 250    0.107   0.130
  77 //    500    0.22    0.42                  500    0.425   0.42  <-(~500)
  78 //   1000    0.86    1.30                 1000    1.72    1.30
  79 //   2500    5.6     4.1  <-(~2070)       2500   11.0     4.1
  80 //   5000   22.3     9.4                  5000   44.7     9.3
  81 //  10000   91.2    20.6                 10000  182      20.5
  82 //
  83 // 1.0*N / N : Newton for N >= 2070 or 1790 >= N >= 1460
  84 // 1.1*N / N : Newton for N >= 1880 or 1790 >= N >= 1320
  85 // 1.2*N / N : Newton for N >= 1250
  86 // 1.3*N / N : Newton for N >= 1010
  87 // 1.4*N / N : Newton for N >=  940
  88 // 1.5*N / N : Newton for N >=  750
  89 // 1.6*N / N : Newton for N >=  625
  90 // 1.7*N / N : Newton for N >=  550
  91 // 1.8*N / N : Newton for N >=  500
  92 // 1.9*N / N : Newton for N >=  500
  93 // 2.0*N / N : Newton for N >=  500
  94   static inline cl_boolean cl_recip_suitable (uintL m, uintL n) // n <= m
  95     { if (n < 500)
  96         return cl_false;
  97       else // when n >= 2100/(m/n)^2, i.e. (m/46)^2 > n
  98         { var uintL mq = floor(m,46);
  99           if ((mq >= bit(16)) || ((uintL)(mq*mq) > n))
 100             return cl_true;
 101           else
 102             return cl_false;
 103         }
 104     }
 105 #endif
 106
 107 void div2adic (uintC a_len, const uintD* a_LSDptr, uintC b_len, const uintD* b_LSDptr, uintD* dest_LSDptr)
 108 {
 109   var uintC lendiff = a_len - b_len;
 110   if (cl_recip_suitable(a_len,b_len))
 111     { // Division using reciprocal (Newton-Hensel algorithm).
 112       CL_ALLOCA_STACK;
 113       // Bestimme Kehrwert c von b mod 2^(intDsize*b_len).
 114       var uintD* c_LSDptr;
 115       num_stack_alloc(b_len,,c_LSDptr=);
 116       recip2adic(b_len,b_LSDptr,c_LSDptr);
 117       // Bestimme q := a * c mod 2^(intDsize*b_len).
 118       var uintD* q_LSDptr;
 119       num_stack_alloc(2*b_len,,q_LSDptr=);
 120       cl_UDS_mul(a_LSDptr,b_len,c_LSDptr,b_len,q_LSDptr);
 121       // Zur Bestimmung des Restes wieder mit b multiplizieren:
 122       var uintD* p_LSDptr;
 123       num_stack_alloc(2*b_len,,p_LSDptr=);
 124       cl_UDS_mul(q_LSDptr,b_len,b_LSDptr,b_len,p_LSDptr);
 125       // Überprüfen, daß p == a mod 2^(intDsize*b_len):
 126       if (compare_loop_msp(a_LSDptr lspop b_len,p_LSDptr lspop b_len,b_len))
 127         cl_abort();
 128       // Quotient q und "Rest" (a-b*q)/2^(intDsize*b_len) ablegen:
 129       copy_loop_lsp(q_LSDptr,dest_LSDptr,b_len);
 130       if (lendiff <= b_len)
 131         { sub_loop_lsp(a_LSDptr lspop b_len,p_LSDptr lspop b_len,dest_LSDptr lspop b_len,lendiff); }
 132         else
 133         { var uintD carry = sub_loop_lsp(a_LSDptr lspop b_len,p_LSDptr lspop b_len,dest_LSDptr lspop b_len,b_len);
 134           copy_loop_lsp(a_LSDptr lspop 2*b_len,dest_LSDptr lspop 2*b_len,lendiff-b_len);
 135           if (carry) { dec_loop_lsp(dest_LSDptr lspop 2*b_len,lendiff-b_len); }
 136         }
 137     }
 138     else
 139     { // Standard division.
 140       var uintD b0inv = div2adic(1,lspref(b_LSDptr,0)); // b'
 141       copy_loop_lsp(a_LSDptr,dest_LSDptr,a_len); // d := a
 142       do { var uintD digit = lspref(dest_LSDptr,0); // nächstes d[j]
 143            digit = mul2adic(b0inv,digit);
 144            // digit = nächstes c[j]
 145            if (a_len <= b_len)
 146              { mulusub_loop_lsp(digit,b_LSDptr,dest_LSDptr,a_len); } // d := d - b * c[j] * beta^j
 147              else
 148              // a_len > b_len, b wird als durch Nullen fortgesetzt gedacht.
 149              { var uintD carry = mulusub_loop_lsp(digit,b_LSDptr,dest_LSDptr,b_len);
 150                if (lspref(dest_LSDptr,b_len) >= carry)
 151                  { lspref(dest_LSDptr,b_len) -= carry; }
 152                else
 153                  { lspref(dest_LSDptr,b_len) -= carry;
 154                    dec_loop_lsp(dest_LSDptr lspop (b_len+1),a_len-(b_len+1));
 155              }   }
 156            // Nun ist lspref(dest_LSDptr,0) = 0.
 157            lspref(dest_LSDptr,0) = digit; // c[j] ablegen
 158            lsshrink(dest_LSDptr); a_len--; // nächstes j
 159          }
 160          until (a_len==lendiff);
 161     }
 162 }
 163 // Bit complexity (N = max(a_len,b_len)): O(M(N)).
 164