cl_DS_div.cc

   1 // cl_UDS_divide().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cl_DS.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cl_N.h"
  13 #include "cln/abort.h"
  14
  15 namespace cln {
  16
  17 // We observe the following timings in seconds:
  18 // Time for dividing a 2*n word number by a n word number:
  19 // OS: Linux 2.2, intDsize==32,        OS: TRU64/4.0, intDsize==64,
  20 // Machine: P-III/450MHz               Machine: EV5/300MHz:
  21 //      n   standard  Newton               standard  Newton
  22 //      10   0.000010  0.000024             0.000036  0.000058
  23 //      30   0.000026  0.000080             0.00012   0.00027
  24 //     100   0.00018   0.00048              0.00084   0.0016
  25 //     300   0.0013    0.0028               0.0062    0.0090
  26 //    1000   0.014     0.019                0.064     0.066  <-(~2200)
  27 //    2000   0.058     0.058  <-(~2000)     0.26      0.20
  28 //    3000   0.20      0.11                 0.57      0.24
  29 //   10000   2.3       0.50                 6.7       1.2
  30 //   30000  24.4       1.2                 62.0       2.8
  31 // Time for dividing a 3*n word number by a n word number:
  32 // OS: Linux 2.2, intDsize==32,        OS: TRU64/4.0, intDsize==64,
  33 // Machine: P-III/450MHz               Machine: EV5/300MHz:
  34 //      n   standard  Newton               standard  Newton
  35 //      10   0.000013  0.000040             0.000063   0.00011
  36 //      30   0.000046  0.00018              0.00024    0.00062
  37 //     100   0.00035   0.0012               0.0016     0.0040
  38 //     300   0.0027    0.0071               0.012      0.021
  39 //    1000   0.029     0.047                0.13       0.16
  40 //    2000   0.12      0.14  <-(~2200)      0.51       0.45  <-(~1600)
  41 //    3000   0.40      0.22                 1.1        0.52
  42 //   10000   4.5       0.76                13.2        2.0
  43 //   30000  42.0       2.8                123.0        6.0
  44 // Time for dividing m digits by n digits:
  45 // OS: Linux 2.2, intDsize==32,        OS: TRU64/4.0, intDsize==64,
  46 // Machine: P-III/450MHz               Machine: EV5/300MHz:
  47 //      n   Newton faster for:         Newton faster for:
  48 //    2-400 never                      never
  49 //    600   never                       670<m<900 (definitly negligible)
  50 //    800   never                       850<m<1500
  51 //   1000   never                      1030<m<2000
  52 //   1200   1400<m<1700                1230<m<2700
  53 //   1500   1590<m<2500                1530<m<5100, 6600<m (ridge negligible)
  54 //   2000   2060<m<3600                2030<m
  55 //   3000   3040<m                     3030<m
  56 //   4000   4030<m                     4030<m
  57 //   5000   5030<m                     5030<m
  58 //   8000   8030<m                     8030<m
  59 // Break-even-point, should be acceptable for both architectures.
  60 // When in doubt, prefer to choose the standard algorithm.
  61 #if CL_USE_GMP
  62   static inline cl_boolean cl_recip_suitable (uintL m, uintL n) // m > n
  63     { if (n < 900)
  64         return cl_false;
  65       else
  66         if (n < 2200)
  67           return (cl_boolean)((m >= n+50) && (m < 2*n-600));
  68         else
  69           return (cl_boolean)(m >= n+30);
  70     }
  71 #else
  72 // Use the old default values from CLN version <= 1.0.3 as a crude estimate.
  73 // They came from the timings for dividing m digits by n digits on an i486/33:
  74 // Dividing 2*N digits by N digits:    Dividing 3*N digits by N digits:
  75 //      N   standard  Newton           standard  Newton
  76 //      10    0.0003  0.0012             0.0006  0.0025
  77 //      25    0.0013  0.0044             0.0026  0.0103
  78 //      50    0.0047  0.0125             0.0092  0.030
  79 //     100    0.017   0.037              0.035   0.089
  80 //     250    0.108   0.146              0.215   0.362
  81 //     500    0.43    0.44  <-(~550)     0.85    1.10
  82 //    1000    1.72    1.32               3.44    3.21  <-(~740)
  83 //    2500   11.2     4.1               23.3     7.9
  84 //    5000   44.3     9.5               89.0    15.6
  85 //   10000  187      20.6              362      33.1
  86 // Time for dividing m digits by n digits:
  87 //   n = 2,3,5,10,25,50,100,250: Newton never faster.
  88 //   n = 400: Newton faster for m >= 440, m < 600
  89 //   n = 500: Newton faster for m >= 530, m < 900
  90 //   n = 600: Newton faster for m >= 630, m < 1250
  91 //   n = 700: Newton faster for m >= 730, m < 1530
  92 //   n = 800: Newton faster for m >= 825, m < 2600 or m >= 5300
  93 //   n = 900: Newton faster for m >= 925, m < 2700 or m >= 3400
  94 //   n = 1000: Newton faster for m >= 1020
  95 //   n = 1500: Newton faster for m >= 1520
  96 //   n = 2000: Newton faster for m >= 2020
  97 //   n = 2500: Newton faster for m >= 2520
  98 //   n = 5000: Newton faster for m >= 5020
  99   static inline cl_boolean cl_recip_suitable (uintL m, uintL n) // m > n
 100     { if (n < 500)
 101         return cl_false;
 102       else
 103         if (n < 1000)
 104           return (cl_boolean)((m >= n+30) && (m < 3*n-600));
 105         else
 106           return (cl_boolean)(m >= n+20);
 107     }
 108 #endif
 109
 110 // Dividiert zwei Unsigned Digit sequences durcheinander.
 111 // UDS_divide(a_MSDptr,a_len,a_LSDptr, b_MSDptr,b_len,b_LSDptr, &q,&r);
 112 // Die UDS a = a_MSDptr/a_len/a_LSDptr (a>=0) wird durch
 113 // die UDS b = b_MSDptr/b_len/b_LSDptr (b>=0) dividiert:
 114 // a = q * b + r mit 0 <= r < b. Bei b=0 Error.
 115 // q der Quotient, r der Rest.
 116 // q = q_MSDptr/q_len/q_LSDptr, r = r_MSDptr/r_len/r_LSDptr beides
 117 // Normalized Unsigned Digit sequences.
 118 // Vorsicht: q_LSDptr <= r_MSDptr,
 119 //           Vorzeichenerweiterung von r kann q zerstören!
 120 //           Vorzeichenerweiterung von q ist erlaubt.
 121 // a und b werden nicht modifiziert.
 122 //
 123 // Methode:
 124 // erst a und b normalisieren: a=[a[m-1],...,a[0]], b=[b[n-1],...,b[0]]
 125 // mit m>=0 und n>0 (Stellensystem der Basis beta=2^intDsize).
 126 // Falls m<n, ist q:=0 und r:=a.
 127 // Falls m>=n=1, Single-Precision-Division:
 128 //   r:=0, j:=m,
 129 //   while j>0 do
 130 //     {Hier (q[m-1]*beta^(m-1)+...+q[j]*beta^j) * b[0] + r*beta^j =
 131 //           = a[m-1]*beta^(m-1)+...+a[j]*beta^j und 0<=r<b[0]<beta}
 132 //     j:=j-1, r:=r*beta+a[j], q[j]:=floor(r/b[0]), r:=r-b[0]*q[j].
 133 //   Normalisiere [q[m-1],...,q[0]], liefert q.
 134 // Falls m>=n>1, Multiple-Precision-Division:
 135 //   Es gilt a/b < beta^(m-n+1).
 136 //   s:=intDsize-1-(Nummer des höchsten Bits in b[n-1]), 0<=s<intDsize.
 137 //   Schiebe a und b um s Bits nach links und kopiere sie dabei. r:=a.
 138 //   r=[r[m],...,r[0]], b=[b[n-1],...,b[0]] mit b[n-1]>=beta/2.
 139 //   Für j=m-n,...,0: {Hier 0 <= r < b*beta^(j+1).}
 140 //     Berechne q* :
 141 //       q* := floor((r[j+n]*beta+r[j+n-1])/b[n-1]).
 142 //       Bei Überlauf (q* >= beta) setze q* := beta-1.
 143 //       Berechne c2 := ((r[j+n]*beta+r[j+n-1]) - q* * b[n-1])*beta + r[j+n-2]
 144 //       und c3 := b[n-2] * q*.
 145 //       {Es ist 0 <= c2 < 2*beta^2, sogar 0 <= c2 < beta^2 falls kein
 146 //        Überlauf aufgetreten war. Ferner 0 <= c3 < beta^2.
 147 //        Bei Überlauf und r[j+n]*beta+r[j+n-1] - q* * b[n-1] >= beta,
 148 //        das heißt c2 >= beta^2, kann man die nächste Abfrage überspringen.}
 149 //       Solange c3 > c2, {hier 0 <= c2 < c3 < beta^2} setze
 150 //         q* := q* - 1, c2 := c2 + b[n-1]*beta, c3 := c3 - b[n-2].
 151 //       Falls q* > 0:
 152 //         Setze r := r - b * q* * beta^j, im einzelnen:
 153 //           [r[n+j],...,r[j]] := [r[n+j],...,r[j]] - q* * [b[n-1],...,b[0]].
 154 //           also: u:=0, for i:=0 to n-1 do
 155 //                         u := u + q* * b[i],
 156 //                         r[j+i]:=r[j+i]-(u mod beta) (+ beta, falls Carry),
 157 //                         u:=u div beta (+ 1, falls bei der Subtraktion Carry)
 158 //                 r[n+j]:=r[n+j]-u.
 159 //           {Da stets u = (q* * [b[i-1],...,b[0]] div beta^i) + 1
 160 //                       < q* + 1 <= beta, läuft der Übertrag u nicht über.}
 161 //         Tritt dabei ein negativer Übertrag auf, so setze q* := q* - 1
 162 //           und [r[n+j],...,r[j]] := [r[n+j],...,r[j]] + [0,b[n-1],...,b[0]].
 163 //     Setze q[j] := q*.
 164 //   Normalisiere [q[m-n],..,q[0]] und erhalte den Quotienten q,
 165 //   Schiebe [r[n-1],...,r[0]] um s Bits nach rechts, normalisiere und
 166 //   erhalte den Rest r.
 167 //   Dabei kann q[j] auf dem Platz von r[n+j] liegen.
 168   void cl_UDS_divide (const uintD* a_MSDptr, uintC a_len, const uintD* a_LSDptr,
 169                       const uintD* b_MSDptr, uintC b_len, const uintD* b_LSDptr,
 170                       uintD* roomptr, // ab roomptr kommen a_len+1 freie Digits
 171                       DS* q_, DS* r_)
 172     { // a normalisieren (a_MSDptr erhöhen, a_len erniedrigen):
 173       while ((a_len>0) && (mspref(a_MSDptr,0)==0)) { msshrink(a_MSDptr); a_len--; }
 174       // b normalisieren (b_MSDptr erhöhen, b_len erniedrigen):
 175       loop
 176         { if (b_len==0) { cl_error_division_by_0(); }
 177           if (mspref(b_MSDptr,0)==0) { msshrink(b_MSDptr); b_len--; }
 178           else break;
 179         }
 180       // jetzt m=a_len >=0 und n=b_len >0.
 181       if (a_len < b_len)
 182         // m<n: Trivialfall, q=0, r:= Kopie von a.
 183         { var uintD* r_MSDptr = roomptr;
 184           var uintD* r_LSDptr = roomptr mspop a_len;
 185           // Speicheraufbau: r_MSDptr/0/r_MSDptr/a_len/r_LSDptr
 186           //                    |     q    |       r      |
 187           copy_loop_lsp(a_LSDptr,r_LSDptr,a_len);
 188           q_->MSDptr = r_MSDptr; q_->len = 0; q_->LSDptr = r_MSDptr; // q = 0, eine NUDS
 189           r_->MSDptr = r_MSDptr; r_->len = a_len; r_->LSDptr = r_LSDptr; // r = Kopie von a, eine NUDS
 190           return;
 191         }
 192       elif (b_len==1)
 193         // n=1: Single-Precision-Division
 194         { // beta^(m-1) <= a < beta^m  ==>  beta^(m-2) <= a/b < beta^m
 195           var uintD* q_MSDptr = roomptr;
 196           var uintD* q_LSDptr = q_MSDptr mspop a_len;
 197           var uintD* r_MSDptr = q_LSDptr;
 198           var uintD* r_LSDptr = r_MSDptr mspop 1;
 199           // Speicheraufbau: q_MSDptr/a_len/q_LSDptr    r_MSDptr/1/r_LSDptr
 200           //                     |      q      |           |     r    |
 201          {var uintD rest = divucopy_loop_msp(mspref(b_MSDptr,0),a_MSDptr,q_MSDptr,a_len); // Division durch b[0]
 202           var uintC r_len;
 203           if (!(rest==0))
 204             { mspref(r_MSDptr,0) = rest; r_len=1; } // Rest als r ablegen
 205             else
 206             { r_MSDptr = r_LSDptr; r_len=0; } // Rest auf 0 normalisieren
 207           if (mspref(q_MSDptr,0)==0)
 208             { msshrink(q_MSDptr); a_len--; } // q normalisieren
 209           q_->MSDptr = q_MSDptr; q_->len = a_len; q_->LSDptr = q_LSDptr; // q ablegen
 210           r_->MSDptr = r_MSDptr; r_->len = r_len; r_->LSDptr = r_LSDptr; // r ablegen
 211           return;
 212         }}
 213       else
 214         // n>1: Multiple-Precision-Division
 215         { // beta^(m-1) <= a < beta^m, beta^(n-1) <= b < beta^n  ==>
 216           // beta^(m-n-1) <= a/b < beta^(m-n+1).
 217           var uintL s;
 218           CL_ALLOCA_STACK;
 219           // s bestimmen:
 220           { var uintD msd = mspref(b_MSDptr,0); // b[n-1]
 221             #if 0
 222             s = 0;
 223             while ((sintD)msd >= 0) { msd = msd<<1; s++; }
 224             #else // ein wenig effizienter, Abfrage auf s=0 vorwegnehmen
 225             if ((sintD)msd < 0)
 226               { s = 0; goto shift_ok; }
 227               else
 228               { integerlengthD(msd, s = intDsize - ); goto shift; }
 229             #endif
 230           }
 231           // 0 <= s < intDsize.
 232           // Kopiere b und schiebe es dabei um s Bits nach links:
 233           if (!(s==0))
 234             shift:
 235             { var uintD* new_b_MSDptr;
 236               var uintD* new_b_LSDptr;
 237               num_stack_alloc(b_len,new_b_MSDptr=,new_b_LSDptr=);
 238               shiftleftcopy_loop_lsp(b_LSDptr,new_b_LSDptr,b_len,s);
 239               b_MSDptr = new_b_MSDptr; b_LSDptr = new_b_LSDptr;
 240             }
 241           shift_ok:
 242           // Wieder b = b_MSDptr/b_len/b_LSDptr.
 243           // Kopiere a und schiebe es dabei um s Bits nach links, erhalte r:
 244          {var uintD* r_MSDptr = roomptr;
 245           var uintD* r_LSDptr = roomptr mspop (a_len+1);
 246           // Speicheraufbau:  r_MSDptr/          a_len+1           /r_LSDptr
 247           //                     |                  r                  |
 248           // später:          q_MSDptr/a_len-b_len+1/r_MSDptr/b_len/r_LSDptr
 249           //                     |           q          |       r      |
 250           if (s==0)
 251             { copy_loop_lsp(a_LSDptr,r_LSDptr,a_len); mspref(r_MSDptr,0) = 0; }
 252             else
 253             { mspref(r_MSDptr,0) = shiftleftcopy_loop_lsp(a_LSDptr,r_LSDptr,a_len,s); }
 254           // Nun r = r_MSDptr/a_len+1/r_LSDptr.
 255           var uintC j = a_len-b_len; // m-n
 256           var uintD* q_MSDptr = r_MSDptr;
 257           var uintC q_len = j+1; // q wird m-n+1 Digits haben
 258           if (cl_recip_suitable(a_len,b_len))
 259             { // Bestimme Kehrwert c von b.
 260               var uintD* c_MSDptr;
 261               var uintD* c_LSDptr;
 262               num_stack_alloc(j+3,c_MSDptr=,c_LSDptr=);
 263               cl_UDS_recip(b_MSDptr,b_len,c_MSDptr,j+1);
 264               // c hat j+3 Digits, | beta^(m+2)/b - c | < beta.
 265               // Mit a' = floor(a/beta^n) multiplizieren, liefert d':
 266               var uintD* d_MSDptr;
 267               UDS_UDS_mul_UDS(j+1,r_MSDptr mspop (j+1), j+3,c_MSDptr mspop (j+3),
 268                               d_MSDptr=,,);
 269               // d' has 2*j+4 digits, d := floor(d'/beta^(j+2)) has j+2 digits.
 270               // | beta^(m+2)/b - c | < beta  ==> (since a < beta^(m+1))
 271               // | beta^(m+2)*a/b - a*c | < beta^(m+2),
 272               // 0 <= a - a'*beta^n < beta^n ==> (since c <= 2*beta^(j+2))
 273               // 0 <= a*c - a'*c*beta^n < 2*beta^(m+2) ==>
 274               // -beta^(m+2) < beta^(m+2)*a/b - a'*c*beta^n < 3*beta^(m+2) ==>
 275               // -1 < a/b - a'*c*beta^(-j-2) < 3 ==>
 276               // -1 < a/b - d'*beta^(-j-2) < 3,
 277               // -1 < d'*beta^(-j-2) - d <= 0 ==>
 278               // -2 < a/b - d < 3 ==>
 279               // -2 <= q - d < 3 ==> |q-d| <= 2.
 280               var uintD* d_LSDptr = d_MSDptr mspop (j+2);
 281               // Zur Bestimmung des Restes wieder mit b multiplizieren:
 282               var uintD* p_MSDptr;
 283               var uintD* p_LSDptr;
 284               UDS_UDS_mul_UDS(j+2,d_LSDptr, b_len,b_LSDptr, p_MSDptr=,,p_LSDptr=);
 285               // d ist um höchstens 2 zu groß, muß also evtl. zweimal um 1
 286               // decrementieren, bis das Produkt <= a wird.
 287               if ((mspref(p_MSDptr,0) > 0) || (compare_loop_msp(p_MSDptr mspop 1,r_MSDptr,a_len+1) > 0))
 288                 { dec_loop_lsp(d_LSDptr,j+2);
 289                   if (subfrom_loop_lsp(b_LSDptr,p_LSDptr,b_len))
 290                     dec_loop_lsp(p_LSDptr lspop b_len,j+2);
 291                   if ((mspref(p_MSDptr,0) > 0) || (compare_loop_msp(p_MSDptr mspop 1,r_MSDptr,a_len+1) > 0))
 292                     { dec_loop_lsp(d_LSDptr,j+2);
 293                       if (subfrom_loop_lsp(b_LSDptr,p_LSDptr,b_len))
 294                         dec_loop_lsp(p_LSDptr lspop b_len,j+2);
 295                       if ((mspref(p_MSDptr,0) > 0) || (compare_loop_msp(p_MSDptr mspop 1,r_MSDptr,a_len+1) > 0))
 296                         cl_abort();
 297                 }   }
 298               // Rest bestimmen:
 299               subfrom_loop_lsp(p_LSDptr,r_LSDptr,a_len+1);
 300               if (test_loop_msp(r_MSDptr,j)) cl_abort();
 301               r_MSDptr = r_LSDptr lspop b_len; // = r_MSDptr mspop (j+1);
 302               // d ist um höchstens 2 zu klein, muß also evtl. zweimal um 1
 303               // incrementieren, bis der Rest < b wird.
 304               if ((lspref(r_MSDptr,0) > 0) || (compare_loop_msp(r_MSDptr,b_MSDptr,b_len) >= 0))
 305                 { inc_loop_lsp(d_LSDptr,j+2);
 306                   if (subfrom_loop_lsp(b_LSDptr,r_LSDptr,b_len))
 307                     lspref(r_LSDptr,b_len) -= 1;
 308                   if ((lspref(r_MSDptr,0) > 0) || (compare_loop_msp(r_MSDptr,b_MSDptr,b_len) >= 0))
 309                     { inc_loop_lsp(d_LSDptr,j+2);
 310                       if (subfrom_loop_lsp(b_LSDptr,r_LSDptr,b_len))
 311                         lspref(r_LSDptr,b_len) -= 1;
 312                       if ((lspref(r_MSDptr,0) > 0) || (compare_loop_msp(r_MSDptr,b_MSDptr,b_len) >= 0))
 313                         cl_abort();
 314                 }   }
 315               // r ist fertig, q := d.
 316               if (mspref(d_MSDptr,0) > 0) cl_abort();
 317               q_len = j+1; copy_loop_msp(d_MSDptr mspop 1,q_MSDptr,q_len);
 318             }
 319             else
 320             { var uintD* r_ptr = r_LSDptr lspop j; // Pointer oberhalb von r[j]
 321               j = j+1;
 322               var uintD b_msd = mspref(b_MSDptr,0); // b[n-1]
 323               var uintD b_2msd = mspref(b_MSDptr,1); // b[n-2]
 324               #if HAVE_DD
 325               var uintDD b_msdd = highlowDD(b_msd,b_2msd); // b[n-1]*beta+b[n-2]
 326               #endif
 327               // Divisions-Schleife: (wird m-n+1 mal durchlaufen)
 328               // j = Herabzähler, b_MSDptr/b_len/b_LSDptr = [b[n-1],...,b[0]], b_len=n,
 329               // r_MSDptr = Pointer auf r[n+j] = Pointer auf q[j],
 330               // r_ptr = Pointer oberhalb von r[j].
 331               do { var uintD q_stern;
 332                    var uintD c1;
 333                    if (mspref(r_MSDptr,0) < b_msd) // r[j+n] < b[n-1] ?
 334                      { // Dividiere r[j+n]*beta+r[j+n-1] durch b[n-1], ohne Überlauf:
 335                        #if HAVE_DD
 336                          divuD(highlowDD(mspref(r_MSDptr,0),mspref(r_MSDptr,1)),b_msd, q_stern=,c1=);
 337                        #else
 338                          divuD(mspref(r_MSDptr,0),mspref(r_MSDptr,1),b_msd, q_stern=,c1=);
 339                        #endif
 340                      }
 341                      else
 342                      { // Überlauf, also r[j+n]*beta+r[j+n-1] >= beta*b[n-1]
 343                        q_stern = bitm(intDsize)-1; // q* = beta-1
 344                        // Teste ob r[j+n]*beta+r[j+n-1] - (beta-1)*b[n-1] >= beta
 345                        // <==> r[j+n]*beta+r[j+n-1] + b[n-1] >= beta*b[n-1]+beta
 346                        // <==> b[n-1] < floor((r[j+n]*beta+r[j+n-1]+b[n-1])/beta) {<= beta !} ist.
 347                        // Wenn ja, direkt zur Subtraktionschleife.
 348                        // (Andernfalls ist r[j+n]*beta+r[j+n-1] - (beta-1)*b[n-1] < beta
 349                        //  <==> floor((r[j+n]*beta+r[j+n-1]+b[n-1])/beta) = b[n-1] ).
 350                        if ((mspref(r_MSDptr,0) > b_msd) || ((c1 = mspref(r_MSDptr,1)+b_msd) < b_msd))
 351                          // r[j+n] >= b[n-1]+1 oder
 352                          // r[j+n] = b[n-1] und Addition r[j+n-1]+b[n-1] gibt Carry ?
 353                          { goto subtract; } // ja -> direkt in die Subtraktion
 354                      }
 355                    // q_stern = q*,
 356                    // c1 = (r[j+n]*beta+r[j+n-1]) - q* * b[n-1] (>=0, <beta).
 357                    #if HAVE_DD
 358                      { var uintDD c2 = highlowDD(c1,mspref(r_MSDptr,2)); // c1*beta+r[j+n-2]
 359                        var uintDD c3 = muluD(b_2msd,q_stern); // b[n-2] * q*
 360                        // Solange c2 < c3, c2 erhöhen, c3 erniedrigen:
 361                        // Rechne dabei mit c3-c2:
 362                        // solange >0, um b[n-1]*beta+b[n-2] erniedrigen.
 363                        // Dies kann wegen b[n-1]*beta+b[n-2] >= beta^2/2
 364                        // höchstens zwei mal auftreten.
 365                        if (c3 > c2)
 366                          { q_stern = q_stern-1; // q* := q* - 1
 367                            if (c3-c2 > b_msdd)
 368                              { q_stern = q_stern-1; } // q* := q* - 1
 369                      }   }
 370                    #else
 371                      // Wie oben, nur mit zweigeteilten c2=[c2hi|c2lo] und c3=[c3hi|c3lo]:
 372                      #define c2hi c1
 373                      { var uintD c2lo = mspref(r_MSDptr,2); // c2hi*beta+c2lo = c1*beta+r[j+n-2]
 374                        var uintD c3hi;
 375                        var uintD c3lo;
 376                        muluD(b_2msd,q_stern, c3hi=,c3lo=); // c3hi*beta+c3lo = b[n-2] * q*
 377                        if ((c3hi > c2hi) || ((c3hi == c2hi) && (c3lo > c2lo)))
 378                          { q_stern = q_stern-1; // q* := q* - 1
 379                            c3hi -= c2hi; if (c3lo < c2lo) { c3hi--; }; c3lo -= c2lo; // c3 := c3-c2
 380                            if ((c3hi > b_msd) || ((c3hi == b_msd) && (c3lo > b_2msd)))
 381                              { q_stern = q_stern-1; } // q* := q* - 1
 382                      }   }
 383                      #undef c2hi
 384                    #endif
 385                    if (!(q_stern==0))
 386                      subtract:
 387                      { // Subtraktionsschleife: r := r - b * q* * beta^j
 388                        var uintD carry = mulusub_loop_lsp(q_stern,b_LSDptr,r_ptr,b_len);
 389                        // Noch r_ptr[-b_len-1] -= carry, d.h. r_MSDptr[0] -= carry
 390                        // durchführen und danach r_MSDptr[0] vergessen:
 391                        if (carry > mspref(r_MSDptr,0))
 392                          // Subtraktion ergab Übertrag
 393                          { q_stern = q_stern-1; // q* := q* - 1
 394                            addto_loop_lsp(b_LSDptr,r_ptr,b_len); // Additionsschleife
 395                            // r[n+j] samt Carry kann vergessen werden...
 396                      }   }
 397                    // Berechnung von q* ist fertig.
 398                    msprefnext(r_MSDptr) = q_stern; // als q[j] ablegen
 399                    r_ptr = r_ptr mspop 1;
 400                  }
 401                  until (--j == 0);
 402             }
 403           // Nun ist q = [q[m-n],..,q[0]] = q_MSDptr/q_len/r_MSDptr
 404           // und r = [r[n-1],...,r[0]] = r_MSDptr/b_len/r_LSDptr.
 405           // q normalisieren und ablegen:
 406           if (mspref(q_MSDptr,0)==0)
 407             { msshrink(q_MSDptr); q_len--; }
 408           q_->MSDptr = q_MSDptr; q_->len = q_len; q_->LSDptr = r_MSDptr;
 409           // Schiebe [r[n-1],...,r[0]] um s Bits nach rechts:
 410           if (!(s==0))
 411             { shiftright_loop_msp(r_MSDptr,b_len,s); }
 412           // r normalisieren und ablegen:
 413           while ((b_len>0) && (mspref(r_MSDptr,0)==0))
 414             { msshrink(r_MSDptr); b_len--; }
 415           r_->MSDptr = r_MSDptr; r_->len = b_len; r_->LSDptr = r_LSDptr;
 416           return;
 417         }}
 418     }
 419 // Bit complexity (N = a_len): O(M(N)).
 420
 421 }  // namespace cln