src/integer/conv/cl_I_to_digits.cc

   1 // UDS_to_digits().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_I.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_DS.h"
  13 #include "cl_I_cached_power.h"
  14
  15 namespace cln {
  16
  17 // Timing für Dezimal-Umwandlung einer Zahl mit N Digits = (N*32) Bits,
  18 // auf einem i486 33 MHz unter Linux:
  19 //    N      standard  dnq(div)  dnq(mul)  combined
  20 //     10     0.00031   0.00043   0.00059   0.00031
  21 //     25     0.00103   0.00125   0.00178   0.00103
  22 //     50     0.0030    0.0034    0.0051    0.0030
  23 //    100     0.0100    0.0108    0.0155    0.0100
  24 //    250     0.054     0.055     0.064     0.054
  25 //    500     0.207     0.209     0.229     0.207
  26 //    750     0.47      0.48      0.47      0.47
  27 //   1000     0.81      0.81      0.86      0.81
  28 //   1250     1.25      1.12      1.20      1.12
  29 //   1500     1.81      1.60      1.64      1.61
  30 //   1750     2.45      2.24      2.15      2.25
  31 //   1940     3.01      3.03      3.12      2.80
  32 //   2000     3.20      3.11      3.30      2.89
  33 //   2500     5.00      4.11      4.38      3.91
  34 //   3000     7.3       5.8       5.7       5.5
  35 //   4000    13.0      12.4      12.9       9.7
  36 //   5000    20.3      15.3      15.1      12.4
  37 //  10000    81.4      57.8      56.4      32.5
  38 //  25000                                 112
  39 //  50000                                 265
  40 // dnq(div) means divide-and-conquer using division by B at the topmost call,
  41 //                threshold = 1015.
  42 // dnq(mul) means divide-and-conquer using multiplication by 1/B at the topmost
  43 //                call, threshold = 2050.
  44 // combined means divide-and-conquer as long as length >= threshold.
  45   const unsigned int cl_digits_div_threshold = 1015;
  46   const int cl_digits_algo = 1;
  47
  48 // like I_to_digits, except that the result has exactly erg_len characters.
  49 static inline void I_to_digits_noshrink (const cl_I& X, uintD base, uintC erg_len, cl_digits* erg)
  50 {
  51   I_to_digits(X,base,erg);
  52   if (erg->len > erg_len) throw runtime_exception();
  53   var uintC count = erg_len - erg->len;
  54   if (count > 0)
  55     { var uintB* ptr = erg->MSBptr;
  56       do { *--ptr = '0'; } while (--count > 0);
  57       erg->MSBptr = ptr; erg->len = erg_len;
  58     }
  59 }
  60
  61 void I_to_digits (const cl_I& X, uintD base, cl_digits* erg)
  62 {
  63 // Methode:
  64 // Umwandlung ins Stellensystem der Basis b geht durch Umwandlung ins Stellen-
  65 // system der Basis b^k (k>=1, b^k<2^intDsize, k maximal) vor sich.
  66 // Aufsuchen von k und b^k aus einer Tabelle.
  67 // Reduktion der UDS zu einer NUDS X.
  68 // Falls X=0: die eine Ziffer 0.
  69 // Falls X>0:
  70 //   Dividiere X durch das Wort b^k,
  71 //   (Single-Precision-Division, vgl. UDS_DIVIDE mit n=1:
  72 //     r:=0, j:=m=Länge(X),
  73 //     while j>0 do
  74 //       j:=j-1, r:=r*beta+X[j], X[j]:=floor(r/b^k), r:=r-b^k*q[j].
  75 //     r=Rest.)
  76 //   zerlege den Rest (mit k-1 Divisionen durch b) in k Ziffern, wandle diese
  77 //   Ziffern einzeln in Ascii um und lege sie an die DIGITS an.
  78 //   Teste auf Speicherüberlauf.
  79 //   X := Quotient.
  80 //   Mache aus X wieder eine NUDS (maximal 1 Nulldigit streichen).
  81 //   Dies solange bis X=0.
  82 //   Streiche die führenden Nullen.
  83       // Aufsuchen von k-1 und b^k aus der Tabelle:
  84       var const power_table_entry* tableptr = &power_table[base-2];
  85       var uintC k = tableptr->k;
  86       var uintD b_hoch_k = tableptr->b_to_the_k; // b^k
  87       var uintB* erg_ptr = erg->LSBptr;
  88       #define next_digit(d)  { *--erg_ptr = (d<10 ? '0'+d : 'A'-10+d); }
  89       // Spezialfälle:
  90       if (zerop(X))
  91         { next_digit(0); goto fertig; } // 0 -> eine Ziffer '0'
  92       else if ((base & (base-1)) == 0)
  93         { // Schneller Algorithmus für Zweierpotenzen
  94           var const uintD* MSDptr;
  95           var uintC len;
  96           var const uintD* LSDptr;
  97           I_to_NDS_nocopy(X, MSDptr=,len=,LSDptr=,false,);
  98           var int b = (base==2 ? 1 : base==4 ? 2 : base==8 ? 3 : base==16 ? 4 : /*base==32*/ 5);
  99           var uintD carry = 0;
 100           var int carrybits = 0;
 101           loop
 102             { if (fixnump(X) && erg->LSBptr-erg_ptr>=cl_value_len)
 103                 break;
 104               if (carrybits >= b)
 105                 { var uintD d = carry & (base-1);
 106                   next_digit(d);
 107                   carry = carry >> b; carrybits -= b;
 108                 }
 109                 else
 110                 { var uintD d = carry;
 111                   if (LSDptr != MSDptr)
 112                     { carry = lsprefnext(LSDptr);
 113                       d |= (carry << carrybits) & (base-1);
 114                       next_digit(d);
 115                       carry = carry >> (b-carrybits); carrybits = intDsize - (b-carrybits);
 116                     }
 117                     else
 118                     { next_digit(d); break; }
 119                 }
 120             }
 121         }
 122       else if (fixnump(X) || TheBignum(X)->length < cl_digits_div_threshold
 123                || !cl_digits_algo)
 124         { // Standard-Algorithmus
 125           CL_ALLOCA_STACK;
 126           var uintD* MSDptr;
 127           var uintC len;
 128           var uintD* LSDptr;
 129           I_to_NDS(X, MSDptr=,len=,LSDptr=);
 130           // normalisiere zu einer NUDS:
 131           if (mspref(MSDptr,0)==0) { msshrink(MSDptr); len--; }
 132           loop
 133             { // Noch die NUDS MSDptr/len/.. mit len>0 abzuarbeiten.
 134               // Single-Precision-Division durch b^k:
 135               var uintD rest = divu_loop_msp(b_hoch_k,MSDptr,len);
 136               // Zerlegen des Restes in seine k Ziffern:
 137               var uintC count = k-1;
 138               if (fixnump(X) && count>cl_value_len-1)
 139                   count = cl_value_len-1;
 140               if ((intDsize>=11) || (count>0))
 141                 // (Bei intDsize>=11 ist wegen b<=36 zwangsläufig
 142                 // k = ceiling(intDsize*log(2)/log(b))-1 >= 2, also count = k-1 > 0.)
 143                 do { var uintD d;
 144                      #if HAVE_DD
 145                        divuD((uintDD)rest,base,rest=,d=);
 146                      #else
 147                        divuD(0,rest,base,rest=,d=);
 148                      #endif
 149                      next_digit(d);
 150                 } until (--count == 0);
 151               next_digit(rest); // letzte der k Ziffern ablegen
 152               // Quotienten normalisieren (max. 1 Digit streichen):
 153               if (mspref(MSDptr,0)==0) { msshrink(MSDptr); len--; if (len==0) break; }
 154         }   }
 155       else
 156         { // Divide-and-conquer:
 157           // Find largest i such that B = base^(k*2^i) satisfies B <= X.
 158           // Divide by B: X = X1*B + X0. Convert X0 to string, fill up
 159           // for k*2^i characters, convert X1 to string. (Have to convert
 160           // X0 first because the conversion may temporarily prepend some
 161           // zero characters.)
 162           var uintC ilen_X = integer_length(X);
 163           var const cached_power_table_entry * p;
 164           var uintC ilen_B;
 165           var uintL i;
 166           for (i = 0; ; i++)
 167             { p = cached_power(base,i);
 168               ilen_B = integer_length(p->base_pow);
 169               if (2*ilen_B >= ilen_X) break;
 170               // 2*ilen_B < ilen_X, so certainly B^2 < X, let's continue with i+1.
 171             }
 172           // 2*ilen_B >= ilen_X, implies X < 2*B^2.
 173           // Of course also X >= B, implies ilen_X >= ilen_B.
 174           #ifdef MUL_REPLACES_DIV
 175           // Divide by B by computing
 176           //   q := floor((X * floor(2^ilen_X/B)) / 2^ilen_X).
 177           // We have q <= floor(X/B) <= q+1, so we may have to increment q.
 178           // Note also that
 179           // floor(2^ilen_X/B) = floor(floor(2^(2*ilen_B)/B)/2^(2*ilen_B-ilen_X))
 180           var cl_I q = (X * (p->inv_base_pow >> (2*ilen_B-ilen_X))) >> ilen_X;
 181           var cl_I r = X - q * p->base_pow;
 182           if (r < 0) throw runtime_exception();
 183           if (r >= p->base_pow)
 184             { q = q+1; r = r - p->base_pow;
 185               if (r >= p->base_pow) throw runtime_exception();
 186             }
 187           #else
 188           var cl_I_div_t q_r = floor2(X,p->base_pow);
 189           var const cl_I& q = q_r.quotient;
 190           var const cl_I& r = q_r.remainder;
 191           #endif
 192           var const cl_I& X1 = q;
 193           var const cl_I& X0 = r;
 194           var uintL B_baselen = (uintL)(k)<<i;
 195           I_to_digits_noshrink(X0,base,B_baselen,erg);
 196           erg->LSBptr -= B_baselen;
 197           I_to_digits(X1,base,erg);
 198           erg->LSBptr += B_baselen;
 199           erg_ptr = erg->MSBptr;
 200         }
 201       #undef next_digit
 202       // Streiche führende Nullen:
 203       while (*erg_ptr == '0') { erg_ptr++; }
 204       fertig:
 205       erg->MSBptr = erg_ptr;
 206       erg->len = erg->LSBptr - erg_ptr;
 207 }
 208 // Bit complexity (N := length(X)): O(log(N)*M(N)).
 209
 210 }  // namespace cln