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
  14 namespace cln {
  15
  16 // Tabelle: enthält zu jeder Basis b (2 <= b <= 36)
  17 // - eine Kettenbruchapproximation num/den von intDsize*log(2)/log(b)
  18 //   (num/den >= intDsize*log(2)/log(b), mit num <= 2^10)
  19 // - k-1 und b^k mit b^k < 2^intDsize, k maximal.
  20   typedef struct { /* uintW num,den; */ uintC k_1; uintD b_hoch_k; } power_table_entry;
  21   static power_table_entry table [36-2+1] = {
  22     #if (intDsize==8)
  23       { /*    8,  1, */ 7-1, 2*2*2*2*2*2*2},
  24       { /*  106, 21, */ 5-1, 3*3*3*3*3},
  25       { /*    4,  1, */ 3-1, 4*4*4},
  26       { /*  789,229, */ 3-1, 5*5*5},
  27       { /*  359,116, */ 3-1, 6*6*6},
  28       { /*  436,153, */ 2-1, 7*7},
  29       { /*    8,  3, */ 2-1, 8*8},
  30       { /*   53, 21, */ 2-1, 9*9},
  31       { /*  525,218, */ 2-1, 10*10},
  32       { /* 1006,435, */ 2-1, 11*11},
  33       { /*  665,298, */ 2-1, 12*12},
  34       { /*  988,457, */ 2-1, 13*13},
  35       { /*  872,415, */ 2-1, 14*14},
  36       { /*  987,482, */ 2-1, 15*15},
  37       { /*    2,  1, */ 1-1, 16},
  38       { /*  869,444, */ 1-1, 17},
  39       { /*  871,454, */ 1-1, 18},
  40       { /*  597,317, */ 1-1, 19},
  41       { /*   87, 47, */ 1-1, 20},
  42       { /*  989,543, */ 1-1, 21},
  43       { /*  949,529, */ 1-1, 22},
  44       { /*  191,108, */ 1-1, 23},
  45       { /*  930,533, */ 1-1, 24},
  46       { /*  789,458, */ 1-1, 25},
  47       { /*  691,406, */ 1-1, 26},
  48       { /*  461,274, */ 1-1, 27},
  49       { /*  218,131, */ 1-1, 28},
  50       { /*  690,419, */ 1-1, 29},
  51       { /*  494,303, */ 1-1, 30},
  52       { /*  633,392, */ 1-1, 31},
  53       { /*    8,  5, */ 1-1, 32},
  54       { /*  766,483, */ 1-1, 33},
  55       { /*  629,400, */ 1-1, 34},
  56       { /*  967,620, */ 1-1, 35},
  57       { /*  359,232, */ 1-1, 36},
  58     #endif
  59     #if (intDsize==16)
  60       { /*   16,  1, */ 15-1, 2*2*2*2*2*2*2*2*2*2*2*2*2*2*2},
  61       { /*  212, 21, */ 10-1, 3*3*3*3*3*3*3*3*3*3},
  62       { /*    8,  1, */  7-1, 4*4*4*4*4*4*4},
  63       { /*  379, 55, */  6-1, 5*5*5*5*5*5},
  64       { /*  359, 58, */  6-1, 6*6*6*6*6*6},
  65       { /*  872,153, */  5-1, 7*7*7*7*7},
  66       { /*   16,  3, */  5-1, 8*8*8*8*8},
  67       { /*  106, 21, */  5-1, 9*9*9*9*9},
  68       { /*  525,109, */  4-1, 10*10*10*10},
  69       { /* 1013,219, */  4-1, 11*11*11*11},
  70       { /*  665,149, */  4-1, 12*12*12*12},
  71       { /*  761,176, */  4-1, 13*13*13*13},
  72       { /*  685,163, */  4-1, 14*14*14*14},
  73       { /*  987,241, */  4-1, 15*15*15*15},
  74       { /*    4,  1, */  3-1, 16*16*16},
  75       { /*  869,222, */  3-1, 17*17*17},
  76       { /*  871,227, */  3-1, 18*18*18},
  77       { /*  113, 30, */  3-1, 19*19*19},
  78       { /*  174, 47, */  3-1, 20*20*20},
  79       { /*   51, 14, */  3-1, 21*21*21},
  80       { /*  653,182, */  3-1, 22*22*22},
  81       { /*  191, 54, */  3-1, 23*23*23},
  82       { /*  677,194, */  3-1, 24*24*24},
  83       { /*  789,229, */  3-1, 25*25*25},
  84       { /*  691,203, */  3-1, 26*26*26},
  85       { /*  461,137, */  3-1, 27*27*27},
  86       { /*  436,131, */  3-1, 28*28*28},
  87       { /*  359,109, */  3-1, 29*29*29},
  88       { /*  988,303, */  3-1, 30*30*30},
  89       { /*  633,196, */  3-1, 31*31*31},
  90       { /*   16,  5, */  3-1, 32*32*32},
  91       { /*  203, 64, */  3-1, 33*33*33},
  92       { /*  629,200, */  3-1, 34*34*34},
  93       { /*  967,310, */  3-1, 35*35*35},
  94       { /*  359,116, */  3-1, 36*36*36},
  95     #endif
  96     #if (intDsize==32)
  97       { /*   32,  1, */ 31-1, 2UL*2UL*2UL*2UL*2UL*2UL*2UL*2UL*2UL*2UL*2UL*2UL*2UL*2UL*2UL*2UL*2UL*2UL*2UL*2UL*2UL*2UL*2UL*2UL*2UL*2UL*2UL*2UL*2UL*2UL*2UL},
  98       { /*  424, 21, */ 20-1, 3UL*3UL*3UL*3UL*3UL*3UL*3UL*3UL*3UL*3UL*3UL*3UL*3UL*3UL*3UL*3UL*3UL*3UL*3UL*3UL},
  99       { /*   16,  1, */ 15-1, 4UL*4UL*4UL*4UL*4UL*4UL*4UL*4UL*4UL*4UL*4UL*4UL*4UL*4UL*4UL},
 100       { /*  758, 55, */ 13-1, 5UL*5UL*5UL*5UL*5UL*5UL*5UL*5UL*5UL*5UL*5UL*5UL*5UL},
 101       { /*  359, 29, */ 12-1, 6UL*6UL*6UL*6UL*6UL*6UL*6UL*6UL*6UL*6UL*6UL*6UL},
 102       { /*   57,  5, */ 11-1, 7UL*7UL*7UL*7UL*7UL*7UL*7UL*7UL*7UL*7UL*7UL},
 103       { /*   32,  3, */ 10-1, 8UL*8UL*8UL*8UL*8UL*8UL*8UL*8UL*8UL*8UL},
 104       { /*  212, 21, */ 10-1, 9UL*9UL*9UL*9UL*9UL*9UL*9UL*9UL*9UL*9UL},
 105       { /*  289, 30, */  9-1, 10UL*10UL*10UL*10UL*10UL*10UL*10UL*10UL*10UL},
 106       { /*  990,107, */  9-1, 11UL*11UL*11UL*11UL*11UL*11UL*11UL*11UL*11UL},
 107       { /*  848, 95, */  8-1, 12UL*12UL*12UL*12UL*12UL*12UL*12UL*12UL},
 108       { /*  761, 88, */  8-1, 13UL*13UL*13UL*13UL*13UL*13UL*13UL*13UL},
 109       { /* 1017,121, */  8-1, 14UL*14UL*14UL*14UL*14UL*14UL*14UL*14UL},
 110       { /*  901,110, */  8-1, 15UL*15UL*15UL*15UL*15UL*15UL*15UL*15UL},
 111       { /*    8,  1, */  7-1, 16UL*16UL*16UL*16UL*16UL*16UL*16UL},
 112       { /*  869,111, */  7-1, 17UL*17UL*17UL*17UL*17UL*17UL*17UL},
 113       { /*  683, 89, */  7-1, 18UL*18UL*18UL*18UL*18UL*18UL*18UL},
 114       { /*  113, 15, */  7-1, 19UL*19UL*19UL*19UL*19UL*19UL*19UL},
 115       { /*  348, 47, */  7-1, 20UL*20UL*20UL*20UL*20UL*20UL*20UL},
 116       { /*   51,  7, */  7-1, 21UL*21UL*21UL*21UL*21UL*21UL*21UL},
 117       { /*  653, 91, */  7-1, 22UL*22UL*22UL*22UL*22UL*22UL*22UL},
 118       { /*  191, 27, */  7-1, 23UL*23UL*23UL*23UL*23UL*23UL*23UL},
 119       { /*  677, 97, */  6-1, 24UL*24UL*24UL*24UL*24UL*24UL},
 120       { /*  379, 55, */  6-1, 25UL*25UL*25UL*25UL*25UL*25UL},
 121       { /*  851,125, */  6-1, 26UL*26UL*26UL*26UL*26UL*26UL},
 122       { /*  922,137, */  6-1, 27UL*27UL*27UL*27UL*27UL*27UL},
 123       { /*  872,131, */  6-1, 28UL*28UL*28UL*28UL*28UL*28UL},
 124       { /*  718,109, */  6-1, 29UL*29UL*29UL*29UL*29UL*29UL},
 125       { /*  150, 23, */  6-1, 30UL*30UL*30UL*30UL*30UL*30UL},
 126       { /*  633, 98, */  6-1, 31UL*31UL*31UL*31UL*31UL*31UL},
 127       { /*   32,  5, */  6-1, 32UL*32UL*32UL*32UL*32UL*32UL},
 128       { /*  203, 32, */  6-1, 33UL*33UL*33UL*33UL*33UL*33UL},
 129       { /*  629,100, */  6-1, 34UL*34UL*34UL*34UL*34UL*34UL},
 130       { /*  967,155, */  6-1, 35UL*35UL*35UL*35UL*35UL*35UL},
 131       { /*  359, 58, */  6-1, 36UL*36UL*36UL*36UL*36UL*36UL},
 132     #endif
 133     #if (intDsize==64)
 134       { /*   64,  1, */ 63-1, 2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL*2ULL},
 135       { /*  848, 21, */ 40-1, 3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL*3ULL},
 136       { /*   32,  1, */ 31-1, 4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL*4ULL},
 137       { /*  634, 23, */ 27-1, 5ULL*5ULL*5ULL*5ULL*5ULL*5ULL*5ULL*5ULL*5ULL*5ULL*5ULL*5ULL*5ULL*5ULL*5ULL*5ULL*5ULL*5ULL*5ULL*5ULL*5ULL*5ULL*5ULL*5ULL*5ULL*5ULL*5ULL},
 138       { /*  718, 29, */ 24-1, 6ULL*6ULL*6ULL*6ULL*6ULL*6ULL*6ULL*6ULL*6ULL*6ULL*6ULL*6ULL*6ULL*6ULL*6ULL*6ULL*6ULL*6ULL*6ULL*6ULL*6ULL*6ULL*6ULL*6ULL},
 139       { /*  114,  5, */ 22-1, 7ULL*7ULL*7ULL*7ULL*7ULL*7ULL*7ULL*7ULL*7ULL*7ULL*7ULL*7ULL*7ULL*7ULL*7ULL*7ULL*7ULL*7ULL*7ULL*7ULL*7ULL*7ULL},
 140       { /*   64,  3, */ 21-1, 8ULL*8ULL*8ULL*8ULL*8ULL*8ULL*8ULL*8ULL*8ULL*8ULL*8ULL*8ULL*8ULL*8ULL*8ULL*8ULL*8ULL*8ULL*8ULL*8ULL*8ULL},
 141       { /*  424, 21, */ 20-1, 9ULL*9ULL*9ULL*9ULL*9ULL*9ULL*9ULL*9ULL*9ULL*9ULL*9ULL*9ULL*9ULL*9ULL*9ULL*9ULL*9ULL*9ULL*9ULL*9ULL},
 142       { /*  289, 15, */ 19-1, 10ULL*10ULL*10ULL*10ULL*10ULL*10ULL*10ULL*10ULL*10ULL*10ULL*10ULL*10ULL*10ULL*10ULL*10ULL*10ULL*10ULL*10ULL*10ULL},
 143       { /* 1018, 55, */ 18-1, 11ULL*11ULL*11ULL*11ULL*11ULL*11ULL*11ULL*11ULL*11ULL*11ULL*11ULL*11ULL*11ULL*11ULL*11ULL*11ULL*11ULL*11ULL},
 144       { /*  607, 34, */ 17-1, 12ULL*12ULL*12ULL*12ULL*12ULL*12ULL*12ULL*12ULL*12ULL*12ULL*12ULL*12ULL*12ULL*12ULL*12ULL*12ULL*12ULL},
 145       { /*  761, 44, */ 17-1, 13ULL*13ULL*13ULL*13ULL*13ULL*13ULL*13ULL*13ULL*13ULL*13ULL*13ULL*13ULL*13ULL*13ULL*13ULL*13ULL*13ULL},
 146       { /*  975, 58, */ 16-1, 14ULL*14ULL*14ULL*14ULL*14ULL*14ULL*14ULL*14ULL*14ULL*14ULL*14ULL*14ULL*14ULL*14ULL*14ULL*14ULL},
 147       { /*  901, 55, */ 16-1, 15ULL*15ULL*15ULL*15ULL*15ULL*15ULL*15ULL*15ULL*15ULL*15ULL*15ULL*15ULL*15ULL*15ULL*15ULL*15ULL},
 148       { /*   16,  1, */ 15-1, 16ULL*16ULL*16ULL*16ULL*16ULL*16ULL*16ULL*16ULL*16ULL*16ULL*16ULL*16ULL*16ULL*16ULL*16ULL},
 149       { /*  595, 38, */ 15-1, 17ULL*17ULL*17ULL*17ULL*17ULL*17ULL*17ULL*17ULL*17ULL*17ULL*17ULL*17ULL*17ULL*17ULL*17ULL},
 150       { /* 1013, 66, */ 15-1, 18ULL*18ULL*18ULL*18ULL*18ULL*18ULL*18ULL*18ULL*18ULL*18ULL*18ULL*18ULL*18ULL*18ULL*18ULL},
 151       { /*  226, 15, */ 15-1, 19ULL*19ULL*19ULL*19ULL*19ULL*19ULL*19ULL*19ULL*19ULL*19ULL*19ULL*19ULL*19ULL*19ULL*19ULL},
 152       { /*  696, 47, */ 14-1, 20ULL*20ULL*20ULL*20ULL*20ULL*20ULL*20ULL*20ULL*20ULL*20ULL*20ULL*20ULL*20ULL*20ULL},
 153       { /*  102,  7, */ 14-1, 21ULL*21ULL*21ULL*21ULL*21ULL*21ULL*21ULL*21ULL*21ULL*21ULL*21ULL*21ULL*21ULL*21ULL},
 154       { /*  775, 54, */ 14-1, 22ULL*22ULL*22ULL*22ULL*22ULL*22ULL*22ULL*22ULL*22ULL*22ULL*22ULL*22ULL*22ULL*22ULL},
 155       { /*  382, 27, */ 14-1, 23ULL*23ULL*23ULL*23ULL*23ULL*23ULL*23ULL*23ULL*23ULL*23ULL*23ULL*23ULL*23ULL*23ULL},
 156       { /* 1019, 73, */ 13-1, 24ULL*24ULL*24ULL*24ULL*24ULL*24ULL*24ULL*24ULL*24ULL*24ULL*24ULL*24ULL*24ULL},
 157       { /*  758, 55, */ 13-1, 25ULL*25ULL*25ULL*25ULL*25ULL*25ULL*25ULL*25ULL*25ULL*25ULL*25ULL*25ULL*25ULL},
 158       { /*  994, 73, */ 13-1, 26ULL*26ULL*26ULL*26ULL*26ULL*26ULL*26ULL*26ULL*26ULL*26ULL*26ULL*26ULL*26ULL},
 159       { /*  673, 50, */ 13-1, 27ULL*27ULL*27ULL*27ULL*27ULL*27ULL*27ULL*27ULL*27ULL*27ULL*27ULL*27ULL*27ULL},
 160       { /*  892, 67, */ 13-1, 28ULL*28ULL*28ULL*28ULL*28ULL*28ULL*28ULL*28ULL*28ULL*28ULL*28ULL*28ULL*28ULL},
 161       { /*  830, 63, */ 13-1, 29ULL*29ULL*29ULL*29ULL*29ULL*29ULL*29ULL*29ULL*29ULL*29ULL*29ULL*29ULL*29ULL},
 162       { /*  300, 23, */ 13-1, 30ULL*30ULL*30ULL*30ULL*30ULL*30ULL*30ULL*30ULL*30ULL*30ULL*30ULL*30ULL*30ULL},
 163       { /*  633, 49, */ 12-1, 31ULL*31ULL*31ULL*31ULL*31ULL*31ULL*31ULL*31ULL*31ULL*31ULL*31ULL*31ULL},
 164       { /*   64,  5, */ 12-1, 32ULL*32ULL*32ULL*32ULL*32ULL*32ULL*32ULL*32ULL*32ULL*32ULL*32ULL*32ULL},
 165       { /*  203, 16, */ 12-1, 33ULL*33ULL*33ULL*33ULL*33ULL*33ULL*33ULL*33ULL*33ULL*33ULL*33ULL*33ULL},
 166       { /*  629, 50, */ 12-1, 34ULL*34ULL*34ULL*34ULL*34ULL*34ULL*34ULL*34ULL*34ULL*34ULL*34ULL*34ULL},
 167       { /*  836, 67, */ 12-1, 35ULL*35ULL*35ULL*35ULL*35ULL*35ULL*35ULL*35ULL*35ULL*35ULL*35ULL*35ULL},
 168       { /*  359, 29, */ 12-1, 36ULL*36ULL*36ULL*36ULL*36ULL*36ULL*36ULL*36ULL*36ULL*36ULL*36ULL*36ULL},
 169     #endif
 170     };
 171
 172 // Timing für Dezimal-Umwandlung einer Zahl mit N Digits = (N*32) Bits,
 173 // auf einem i486 33 MHz unter Linux:
 174 //    N      standard  dnq(div)  dnq(mul)  combined
 175 //     10     0.00031   0.00043   0.00059   0.00031
 176 //     25     0.00103   0.00125   0.00178   0.00103
 177 //     50     0.0030    0.0034    0.0051    0.0030
 178 //    100     0.0100    0.0108    0.0155    0.0100
 179 //    250     0.054     0.055     0.064     0.054
 180 //    500     0.207     0.209     0.229     0.207
 181 //    750     0.47      0.48      0.47      0.47
 182 //   1000     0.81      0.81      0.86      0.81
 183 //   1250     1.25      1.12      1.20      1.12
 184 //   1500     1.81      1.60      1.64      1.61
 185 //   1750     2.45      2.24      2.15      2.25
 186 //   1940     3.01      3.03      3.12      2.80
 187 //   2000     3.20      3.11      3.30      2.89
 188 //   2500     5.00      4.11      4.38      3.91
 189 //   3000     7.3       5.8       5.7       5.5
 190 //   4000    13.0      12.4      12.9       9.7
 191 //   5000    20.3      15.3      15.1      12.4
 192 //  10000    81.4      57.8      56.4      32.5
 193 //  25000                                 112
 194 //  50000                                 265
 195 // dnq(div) means divide-and-conquer using division by B at the topmost call,
 196 //                threshold = 1015.
 197 // dnq(mul) means divide-and-conquer using multiplication by 1/B at the topmost
 198 //                call, threshold = 2050.
 199 // combined means divide-and-conquer as long as length >= threshold.
 200   const unsigned int cl_digits_div_threshold = 1015;
 201   //#define MUL_REPLACES_DIV
 202   const int cl_digits_algo = 1;
 203
 204 // Tabelle: enthält zu jeder Basis b (2 <= b <= 36)
 205 // NULL oder einen Vektor von lazy berechneten b^(k*2^i) und 1/b^(k*2^i).
 206   typedef struct cached_power_table_entry {
 207     ALLOCATE_ANYWHERE(cached_power_table_entry)
 208     cl_I base_pow; // 0 or b^(k*2^i)
 209     #ifdef MUL_REPLACES_DIV
 210     cl_I inv_base_pow; // if base_pow: floor(2^(2*integer_length(base_pow))/base_pow)
 211     #endif
 212   } cached_power_table_entry;
 213   struct cached_power_table {
 214         cached_power_table_entry element[30];
 215         // Constructor and destructor - nothing special.
 216         cached_power_table () {}
 217         ~cached_power_table () {}
 218         // Allocation and deallocation.
 219         void* operator new (size_t size) { return malloc_hook(size); }
 220         void operator delete (void* ptr) { free_hook(ptr); }
 221   };
 222   static cached_power_table* ctable [36-2+1] =
 223     { NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
 224       NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
 225       NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
 226       NULL, NULL, NULL, NULL, NULL
 227     };
 228   static const cached_power_table_entry * cached_power (uintD base, uintL i)
 229     { var cached_power_table* ptr;
 230       if (!(ptr = ctable[base-2]))
 231         { ctable[base-2] = ptr = new cached_power_table (); }
 232       var uintL j;
 233       for (j = 0; j <= i; j++)
 234         if (zerop(ptr->element[j].base_pow))
 235           { // Compute b^(k*2^j) and its inverse.
 236             cl_I x =
 237               (j==0 ? (cl_I)(unsigned long)(table[base-2].b_hoch_k)
 238                     : ptr->element[j-1].base_pow * ptr->element[j-1].base_pow
 239               );
 240             ptr->element[j].base_pow = x;
 241             #ifdef MUL_REPLACES_DIV
 242             ptr->element[j].inv_base_pow = floor1(ash(1,2*integer_length(x)),x);
 243             #endif
 244           }
 245       return &ptr->element[i];
 246     }
 247   AT_DESTRUCTION(cached_power)
 248     { for (var uintD base = 2; base <= 36; base++)
 249         { var cached_power_table* ptr = ctable[base-2];
 250           if (ptr)
 251             { delete ptr; ctable[base-2] = NULL; }
 252         }
 253     }
 254
 255 // like I_to_digits, except that the result has exactly erg_len characters.
 256 static inline void I_to_digits_noshrink (const cl_I& X, uintD base, uintL erg_len, cl_digits* erg)
 257 {
 258   I_to_digits(X,base,erg);
 259   if (erg->len > erg_len) cl_abort();
 260   var uintL count = erg_len - erg->len;
 261   if (count > 0)
 262     { var uintB* ptr = erg->MSBptr;
 263       do { *--ptr = '0'; } while (--count > 0);
 264       erg->MSBptr = ptr; erg->len = erg->len;
 265     }
 266 }
 267
 268 void I_to_digits (const cl_I& X, uintD base, cl_digits* erg)
 269 {
 270 // Methode:
 271 // Umwandlung ins Stellensystem der Basis b geht durch Umwandlung ins Stellen-
 272 // system der Basis b^k (k>=1, b^k<2^intDsize, k maximal) vor sich.
 273 // Aufsuchen von k und b^k aus einer Tabelle.
 274 // Reduktion der UDS zu einer NUDS X.
 275 // Falls X=0: die eine Ziffer 0.
 276 // Falls X>0:
 277 //   Dividiere X durch das Wort b^k,
 278 //   (Single-Precision-Division, vgl. UDS_DIVIDE mit n=1:
 279 //     r:=0, j:=m=Länge(X),
 280 //     while j>0 do
 281 //       j:=j-1, r:=r*beta+X[j], X[j]:=floor(r/b^k), r:=r-b^k*q[j].
 282 //     r=Rest.)
 283 //   zerlege den Rest (mit k-1 Divisionen durch b) in k Ziffern, wandle diese
 284 //   Ziffern einzeln in Ascii um und lege sie an die DIGITS an.
 285 //   Teste auf Speicherüberlauf.
 286 //   X := Quotient.
 287 //   Mache aus X wieder eine NUDS (maximal 1 Nulldigit streichen).
 288 //   Dies solange bis X=0.
 289 //   Streiche die führenden Nullen.
 290       // Aufsuchen von k-1 und b^k aus der Tabelle:
 291       var power_table_entry* tableptr = &table[base-2];
 292       var uintC k_1 = tableptr->k_1; // k-1
 293       var uintD b_hoch_k = tableptr->b_hoch_k; // b^k
 294       var uintB* erg_ptr = erg->LSBptr;
 295       #define next_digit(d)  { *--erg_ptr = (d<10 ? '0'+d : 'A'-10+d); }
 296       // Spezialfälle:
 297       if (zerop(X))
 298         { next_digit(0); goto fertig; } // 0 -> eine Ziffer '0'
 299       else if ((base & (base-1)) == 0)
 300         { // Schneller Algorithmus für Zweierpotenzen
 301           var const uintD* MSDptr;
 302           var uintC len;
 303           var const uintD* LSDptr;
 304           I_to_NDS_nocopy(X, MSDptr=,len=,LSDptr=,cl_false,);
 305           var int b = (base==2 ? 1 : base==4 ? 2 : base==8 ? 3 : base==16 ? 4 : /*base==32*/ 5);
 306           var uintD carry = 0;
 307           var int carrybits = 0;
 308           loop
 309             { if (carrybits >= b)
 310                 { var uintD d = carry & (base-1);
 311                   next_digit(d);
 312                   carry = carry >> b; carrybits -= b;
 313                 }
 314                 else
 315                 { var uintD d = carry;
 316                   if (LSDptr != MSDptr)
 317                     { carry = lsprefnext(LSDptr);
 318                       d |= (carry << carrybits) & (base-1);
 319                       next_digit(d);
 320                       carry = carry >> (b-carrybits); carrybits = intDsize - (b-carrybits);
 321                     }
 322                     else
 323                     { next_digit(d); break; }
 324                 }
 325             }
 326         }
 327       else if (fixnump(X) || TheBignum(X)->length < cl_digits_div_threshold
 328                || !cl_digits_algo)
 329         { // Standard-Algorithmus
 330           CL_ALLOCA_STACK;
 331           var uintD* MSDptr;
 332           var uintC len;
 333           var uintD* LSDptr;
 334           I_to_NDS(X, MSDptr=,len=,LSDptr=);
 335           // normalisiere zu einer NUDS:
 336           if (mspref(MSDptr,0)==0) { msshrink(MSDptr); len--; }
 337           loop
 338             { // Noch die NUDS MSDptr/len/.. mit len>0 abzuarbeiten.
 339               // Single-Precision-Division durch b^k:
 340               var uintD rest = divu_loop_msp(b_hoch_k,MSDptr,len);
 341               // Zerlegen des Restes in seine k Ziffern:
 342               var uintC count = k_1;
 343               if ((intDsize>=11) || (count>0))
 344                 // (Bei intDsize>=11 ist wegen b<=36 zwangsläufig
 345                 // k = ceiling(intDsize*log(2)/log(b))-1 >= 2, also count = k_1 > 0.)
 346                 do { var uintD d;
 347                      #if HAVE_DD
 348                        divuD((uintDD)rest,base,rest=,d=);
 349                      #else
 350                        divuD(0,rest,base,rest=,d=);
 351                      #endif
 352                      next_digit(d);
 353                    }
 354                    until (--count == 0);
 355               next_digit(rest); // letzte der k Ziffern ablegen
 356               // Quotienten normalisieren (max. 1 Digit streichen):
 357               if (mspref(MSDptr,0)==0) { msshrink(MSDptr); len--; if (len==0) break; }
 358         }   }
 359       else
 360         { // Divide-and-conquer:
 361           // Find largest i such that B = base^(k*2^i) satisfies B <= X.
 362           // Divide by B: X = X1*B + X0. Convert X0 to string, fill up
 363           // for k*2^i characters, convert X1 to string. (Have to convert
 364           // X0 first because the conversion may temporarily prepend some
 365           // zero characters.)
 366           var uintL ilen_X = integer_length(X);
 367           var const cached_power_table_entry * p;
 368           var uintL ilen_B;
 369           var uintL i;
 370           for (i = 0; ; i++)
 371             { p = cached_power(base,i);
 372               ilen_B = integer_length(p->base_pow);
 373               if (2*ilen_B >= ilen_X) break;
 374               // 2*ilen_B < ilen_X, so certainly B^2 < X, let's continue with i+1.
 375             }
 376           // 2*ilen_B >= ilen_X, implies X < 2*B^2.
 377           // Of course also X >= B, implies ilen_X >= ilen_B.
 378           #ifdef MUL_REPLACES_DIV
 379           // Divide by B by computing
 380           //   q := floor((X * floor(2^ilen_X/B)) / 2^ilen_X).
 381           // We have q <= floor(X/B) <= q+1, so we may have to increment q.
 382           // Note also that
 383           // floor(2^ilen_X/B) = floor(floor(2^(2*ilen_B)/B)/2^(2*ilen_B-ilen_X))
 384           var cl_I q = (X * (p->inv_base_pow >> (2*ilen_B-ilen_X))) >> ilen_X;
 385           var cl_I r = X - q * p->base_pow;
 386           if (r < 0) cl_abort();
 387           if (r >= p->base_pow)
 388             { q = q+1; r = r - p->base_pow;
 389               if (r >= p->base_pow) cl_abort();
 390             }
 391           #else
 392           var cl_I_div_t q_r = floor2(X,p->base_pow);
 393           var const cl_I& q = q_r.quotient;
 394           var const cl_I& r = q_r.remainder;
 395           #endif
 396           var const cl_I& X1 = q;
 397           var const cl_I& X0 = r;
 398           var uintL B_baselen = (uintL)(k_1+1)<<i;
 399           I_to_digits_noshrink(X0,base,B_baselen,erg);
 400           erg->LSBptr -= B_baselen;
 401           I_to_digits(X1,base,erg);
 402           erg->LSBptr += B_baselen;
 403           erg_ptr = erg->MSBptr;
 404         }
 405       #undef next_digit
 406       // Streiche führende Nullen:
 407       while (*erg_ptr == '0') { erg_ptr++; }
 408       fertig:
 409       erg->MSBptr = erg_ptr;
 410       erg->len = erg->LSBptr - erg_ptr;
 411 }
 412 // Bit complexity (N := length(X)): O(log(N)*M(N)).
 413
 414 }  // namespace cln