src/base/digitseq/cl_DS_mul_kara.h

   1
   2   static void mulu_karatsuba (const uintD* sourceptr1, uintC len1,
   3                               const uintD* sourceptr2, uintC len2,
   4                               uintD* destptr)
   5     // Karatsuba-Multiplikation
   6     // Prinzip: (x1*b^k+x0) * (y1*b^k+y0)
   7     //        = x1*y1 * b^2k + ((x1+x0)*(y1+y0)-x1*y1-x0*y0) * b^k + x0*y0
   8     // Methode 1 (Collins/Loos, Degel):
   9     // source2 wird in floor(len2/len1) einzelne UDS mit je einer
  10     // Länge len3 (len1 <= len3 < 2*len1) unterteilt,
  11     // jeweils k=floor(len3/2).
  12     // Methode 2 (Haible):
  13     // source2 wird in ceiling(len2/len1) einzelne UDS mit je einer
  14     // Länge len3 (0 < len3 <= len1) unterteilt, jeweils k=floor(len1/2).
  15     // Aufwand für die hinteren Einzelteile:
  16     // bei beiden Methoden jeweils 3*len1^2.
  17     // Aufwand für das vorderste Teil (alles, falls len1 <= len2 < 2*len1)
  18     // mit r = len1, s = (len2 mod len1) + len1 (>= len1, < 2*len1):
  19     // bei Methode 1:
  20     //                       |   :       |  r
  21     //                    |      :       |  s
  22     //      (r-s/2)*s/2 + s/2*s/2 + s/2*s/2 = r*s/2 + s^2/4 .
  23     // bei Methode 2:
  24     //                       |     :     |  r
  25     //                    |  |     :     |  s
  26     //      (s-r)*r + r/2*r/2 + r/2*r/2 + r/2*r/2 = r*s - r^2/4 .
  27     // Wegen (r*s/2 + s^2/4) - (r*s - r^2/4) = (r-s)^2/4 >= 0
  28     // ist Methode 2 günstiger.
  29     // Denkfehler! Dies gilt - wenn überhaupt - nur knapp oberhalb des
  30     // Break-Even-Points.
  31     // Im allgemeinen ist der Multiplikationsaufwand für zwei Zahlen der
  32     // Längen u bzw. v nämlich gegeben durch  min(u,v)^c * max(u,v),
  33     // wobei  c = log3/log2 - 1 = 0.585...
  34     // Dadurch wird der Aufwand in Abhängigkeit des Parameters t = k,
  35     // r/2 <= t <= s/2 (der einzig sinnvolle Bereich), zu
  36     // (r-t)^c*(s-t) + t^c*(s-t) + t^(1+c).
  37     // Dessen Optimum liegt (im Bereich r <= s <= 2*r)
  38     // - im klassischen Fall c=1 tatsächlich stets bei t=r/2 [Methode 2],
  39     // - im Karatsuba-Fall c=0.6 aber offenbar bei t=s/2 [Methode 1]
  40     //   oder ganz knapp darunter.
  41     // Auch erweist sich Methode 1 im Experiment als effizienter.
  42     // Daher implementieren wir Methode 1 :
  43     { // Es ist 2 <= len1 <= len2.
  44       // Spezialfall Quadrieren abfangen (häufig genug, daß sich das lohnt):
  45       if (sourceptr1 == sourceptr2)
  46         if (len1 == len2)
  47           { mulu_karatsuba_square(sourceptr1,len1,destptr); return; }
  48       var cl_boolean first_part = cl_true; // Flag, ob jetzt das erste Teilprodukt berechnet wird
  49       if (len2 >= 2*len1)
  50         { CL_SMALL_ALLOCA_STACK;
  51           // Teilprodukte von jeweils len1 mal len1 Digits bilden:
  52           var uintC k_lo = floor(len1,2); // Länge der Low-Teile: floor(len1/2) >0
  53           var uintC k_hi = len1 - k_lo; // Länge der High-Teile: ceiling(len1/2) >0
  54           // Es gilt k_lo <= k_hi <= len1, k_lo + k_hi = len1.
  55           // Summe x1+x0 berechnen:
  56           var uintD* sum1_MSDptr;
  57           var uintC sum1_len = k_hi; // = max(k_lo,k_hi)
  58           var uintD* sum1_LSDptr;
  59           num_stack_small_alloc_1(sum1_len,sum1_MSDptr=,sum1_LSDptr=);
  60           {var uintD carry = // Hauptteile von x1 und x0 addieren:
  61              add_loop_lsp(sourceptr1 lspop k_lo,sourceptr1,sum1_LSDptr,k_lo);
  62            if (!(k_lo==k_hi))
  63              // noch k_hi-k_lo = 1 Digits abzulegen
  64              { mspref(sum1_MSDptr,0) = lspref(sourceptr1,len1-1); // = lspref(sourceptr1,2*k_lo)
  65                if (!(carry==0)) { if (++(mspref(sum1_MSDptr,0)) == 0) carry=1; else carry=0; }
  66              }
  67            if (carry) { lsprefnext(sum1_MSDptr) = 1; sum1_len++; }
  68           }
  69          {  // Platz für Summe y1+y0 belegen:
  70             var uintC sum2_maxlen = k_hi+1;
  71             var uintD* sum2_LSDptr;
  72             num_stack_small_alloc(sum2_maxlen,,sum2_LSDptr=);
  73             // Platz für Produkte x0*y0, x1*y1 belegen:
  74           { var uintD* prod_MSDptr;
  75             var uintD* prod_LSDptr;
  76             var uintD* prodhi_LSDptr;
  77             num_stack_small_alloc(2*(uintL)len1,prod_MSDptr=,prod_LSDptr=);
  78             prodhi_LSDptr = prod_LSDptr lspop 2*k_lo;
  79             // prod_MSDptr/2*len1/prod_LSDptr wird zuerst die beiden
  80             // Produkte x1*y1 in prod_MSDptr/2*k_hi/prodhi_LSDptr
  81             //      und x0*y0 in prodhi_LSDptr/2*k_lo/prod_LSDptr,
  82             // dann das Produkt (b^k*x1+x0)*(b^k*y1+y0) enthalten.
  83             // Platz fürs Produkt (x1+x0)*(y1+y0) belegen:
  84            {var uintD* prodmid_MSDptr;
  85             var uintD* prodmid_LSDptr;
  86             num_stack_small_alloc(sum1_len+sum2_maxlen,prodmid_MSDptr=,prodmid_LSDptr=);
  87             // Schleife über die hinteren Einzelteile:
  88             do { // Produkt x0*y0 berechnen:
  89                  cl_UDS_mul(sourceptr1,k_lo,sourceptr2,k_lo,prod_LSDptr);
  90                  // Produkt x1*y1 berechnen:
  91                  cl_UDS_mul(sourceptr1 lspop k_lo,k_hi,sourceptr2 lspop k_lo,k_hi,prodhi_LSDptr);
  92                  // Summe y1+y0 berechnen:
  93                 {var uintC sum2_len = k_hi; // = max(k_lo,k_hi)
  94                  var uintD* sum2_MSDptr = sum2_LSDptr lspop sum2_len;
  95                  {var uintD carry = // Hauptteile von y1 und y0 addieren:
  96                     add_loop_lsp(sourceptr2 lspop k_lo,sourceptr2,sum2_LSDptr,k_lo);
  97                   if (!(k_lo==k_hi))
  98                     // noch k_hi-k_lo = 1 Digits abzulegen
  99                     { mspref(sum2_MSDptr,0) = lspref(sourceptr2,len1-1); // = lspref(sourceptr2,2*k_lo)
 100                       if (!(carry==0)) { if (++(mspref(sum2_MSDptr,0)) == 0) carry=1; else carry=0; }
 101                     }
 102                   if (carry) { lsprefnext(sum2_MSDptr) = 1; sum2_len++; }
 103                  }
 104                  // Produkt (x1+x0)*(y1+y0) berechnen:
 105                  cl_UDS_mul(sum1_LSDptr,sum1_len,sum2_LSDptr,sum2_len,prodmid_LSDptr);
 106                  // Das Produkt beansprucht  2*k_hi + (0 oder 1) <= sum1_len + sum2_len  Digits.
 107                  {var uintC prodmid_len = sum1_len+sum2_len;
 108                   // Davon x1*y1 abziehen:
 109                   {var uintD carry =
 110                      subfrom_loop_lsp(prodhi_LSDptr,prodmid_LSDptr,2*k_hi);
 111                    // Falls Carry: Produkt beansprucht 2*k_hi+1 Digits.
 112                    // Carry um maximal 1 Digit weitertragen:
 113                    if (!(carry==0)) { lspref(prodmid_LSDptr,2*k_hi) -= 1; }
 114                   }
 115                   // Und x0*y0 abziehen:
 116                   {var uintD carry =
 117                      subfrom_loop_lsp(prod_LSDptr,prodmid_LSDptr,2*k_lo);
 118                    // Carry um maximal prodmid_len-2*k_lo Digits weitertragen:
 119                    if (!(carry==0))
 120                      { dec_loop_lsp(prodmid_LSDptr lspop 2*k_lo,prodmid_len-2*k_lo); }
 121                   }
 122                   // prodmid_LSDptr[-prodmid_len..-1] enthält nun x0*y1+x1*y0.
 123                   // Dies wird zu prod = x1*y1*b^(2*k) + x0*y0 addiert:
 124                   {var uintD carry =
 125                      addto_loop_lsp(prodmid_LSDptr,prod_LSDptr lspop k_lo,prodmid_len);
 126                      // (Benutze dabei k_lo+prodmid_len <= k_lo+2*(k_hi+1) = 2*len1-k_lo+2 <= 2*len1 .)
 127                    if (!(carry==0))
 128                      { inc_loop_lsp(prod_LSDptr lspop (k_lo+prodmid_len),2*len1-(k_lo+prodmid_len)); }
 129                 }}}
 130                  // Das Teilprodukt zum Gesamtprodukt addieren:
 131                  if (first_part)
 132                    { copy_loop_lsp(prod_LSDptr,destptr,2*len1);
 133                      destptr = destptr lspop len1;
 134                      first_part = cl_false;
 135                    }
 136                    else
 137                    { var uintD carry =
 138                        addto_loop_lsp(prod_LSDptr,destptr,len1);
 139                      destptr = destptr lspop len1;
 140                      copy_loop_lsp(prod_LSDptr lspop len1,destptr,len1);
 141                      if (!(carry==0)) { inc_loop_lsp(destptr,len1); }
 142                    }
 143                  sourceptr2 = sourceptr2 lspop len1; len2 -= len1;
 144                }
 145                while (len2 >= 2*len1);
 146          }}}
 147         }
 148       // Nun ist len1 <= len2 < 2*len1.
 149       // letztes Teilprodukt von len1 mal len2 Digits bilden:
 150      {CL_SMALL_ALLOCA_STACK;
 151       var uintD* prod_MSDptr;
 152       var uintC prod_len = len1+len2;
 153       var uintD* prod_LSDptr;
 154       num_stack_small_alloc((uintL)prod_len,prod_MSDptr=,prod_LSDptr=);
 155       { var uintC k_hi = floor(len2,2); // Länge der High-Teile: floor(len2/2) >0
 156         var uintC k_lo = len2 - k_hi; // Länge der Low-Teile: ceiling(len2/2) >0
 157         // Es gilt k_hi <= k_lo <= len1 <= len2, k_lo + k_hi = len2.
 158         var uintC x1_len = len1-k_lo; // <= len2-k_lo = k_hi <= k_lo
 159         // Summe x1+x0 berechnen:
 160         var uintD* sum1_MSDptr;
 161         var uintC sum1_len = k_lo; // = max(k_lo,k_hi)
 162         var uintD* sum1_LSDptr;
 163         num_stack_small_alloc_1(sum1_len,sum1_MSDptr=,sum1_LSDptr=);
 164         {var uintD carry = // x1 und unteren Teil von x0 addieren:
 165            add_loop_lsp(sourceptr1 lspop k_lo,sourceptr1,sum1_LSDptr,x1_len);
 166          // und den oberen Teil von x0 dazu:
 167          copy_loop_lsp(sourceptr1 lspop x1_len,sum1_LSDptr lspop x1_len,k_lo-x1_len);
 168          if (!(carry==0))
 169            { carry = inc_loop_lsp(sum1_LSDptr lspop x1_len,k_lo-x1_len);
 170              if (carry) { lsprefnext(sum1_MSDptr) = 1; sum1_len++; }
 171            }
 172         }
 173        {// Summe y1+y0 berechnen:
 174         var uintD* sum2_MSDptr;
 175         var uintC sum2_len = k_lo; // = max(k_lo,k_hi)
 176         var uintD* sum2_LSDptr;
 177         num_stack_small_alloc_1(sum2_len,sum2_MSDptr=,sum2_LSDptr=);
 178         {var uintD carry = // Hauptteile von y1 und y0 addieren:
 179            add_loop_lsp(sourceptr2 lspop k_lo,sourceptr2,sum2_LSDptr,k_hi);
 180          if (!(k_lo==k_hi))
 181            // noch k_lo-k_hi = 1 Digits abzulegen
 182            { mspref(sum2_MSDptr,0) = lspref(sourceptr2,k_lo-1); // = lspref(sourceptr2,k_hi)
 183              if (!(carry==0)) { if (++(mspref(sum2_MSDptr,0)) == 0) carry=1; else carry=0; }
 184            }
 185          if (carry) { lsprefnext(sum2_MSDptr) = 1; sum2_len++; }
 186         }
 187         // Platz für Produkte x0*y0, x1*y1:
 188         { var uintC prodhi_len = x1_len+k_hi;
 189           var uintD* prodhi_LSDptr = prod_LSDptr lspop 2*k_lo;
 190           // prod_MSDptr/len1+len2/prod_LSDptr wird zuerst die beiden
 191           // Produkte x1*y1 in prod_MSDptr/x1_len+k_hi/prodhi_LSDptr
 192           //      und x0*y0 in prodhi_LSDptr/2*k_lo/prod_LSDptr,
 193           // dann das Produkt (b^k*x1+x0)*(b^k*y1+y0) enthalten.
 194           // Platz fürs Produkt (x1+x0)*(y1+y0) belegen:
 195          {var uintD* prodmid_MSDptr;
 196           var uintC prodmid_len = sum1_len+sum2_len;
 197           var uintD* prodmid_LSDptr;
 198           num_stack_small_alloc(prodmid_len,prodmid_MSDptr=,prodmid_LSDptr=);
 199           // Produkt (x1+x0)*(y1+y0) berechnen:
 200           cl_UDS_mul(sum1_LSDptr,sum1_len,sum2_LSDptr,sum2_len,prodmid_LSDptr);
 201           // Das Produkt beansprucht  2*k_lo + (0 oder 1) <= sum1_len + sum2_len = prodmid_len  Digits.
 202           // Produkt x0*y0 berechnen:
 203           cl_UDS_mul(sourceptr1,k_lo,sourceptr2,k_lo,prod_LSDptr);
 204           // Produkt x1*y1 berechnen:
 205           if (!(x1_len==0))
 206             { cl_UDS_mul(sourceptr1 lspop k_lo,x1_len,sourceptr2 lspop k_lo,k_hi,prodhi_LSDptr);
 207              // Und x1*y1 abziehen:
 208              {var uintD carry =
 209                 subfrom_loop_lsp(prodhi_LSDptr,prodmid_LSDptr,prodhi_len);
 210               // Carry um maximal prodmid_len-prodhi_len Digits weitertragen:
 211               if (!(carry==0))
 212                 { dec_loop_lsp(prodmid_LSDptr lspop prodhi_len,prodmid_len-prodhi_len); }
 213             }}
 214             else
 215             // Produkt x1*y1=0, nichts abzuziehen
 216             { clear_loop_lsp(prodhi_LSDptr,prodhi_len); }
 217           // Und x0*y0 abziehen:
 218           {var uintD carry =
 219              subfrom_loop_lsp(prod_LSDptr,prodmid_LSDptr,2*k_lo);
 220            // Falls Carry: Produkt beansprucht 2*k_lo+1 Digits.
 221            // Carry um maximal 1 Digit weitertragen:
 222            if (!(carry==0)) { lspref(prodmid_LSDptr,2*k_lo) -= 1; }
 223           }
 224           // prodmid_LSDptr[-prodmid_len..-1] enthält nun x0*y1+x1*y0.
 225           // Dies ist < b^k_lo * b^k_hi + b^x1_len * b^k_lo
 226           //          = b^len2 + b^len1 <= 2 * b^len2,
 227           // paßt also in len2+1 Digits.
 228           // Im Fall x1_len=0 ist es sogar < b^k_lo * b^k_hi = b^len2,
 229           // es paßt also in len2 Digits.
 230           // prodmid_len, wenn möglich, um maximal 2 verkleinern:
 231           // (benutzt prodmid_len >= 2*k_lo >= len2 >= 2)
 232           if (mspref(prodmid_MSDptr,0)==0)
 233             { prodmid_len--;
 234               if (mspref(prodmid_MSDptr,1)==0) { prodmid_len--; }
 235             }
 236           // Nun ist k_lo+prodmid_len <= len1+len2 .
 237           // (Denn es war prodmid_len = sum1_len+sum2_len <= 2*(k_lo+1)
 238           //  <= len2+3, und nach 2-maliger Verkleinerung jedenfalls
 239           //  prodmid_len <= len2+1. Im Falle k_lo < len1 also
 240           //  k_lo + prodmid_len <= (len1-1)+(len2+1) = len1+len2.
 241           //  Im Falle k_lo = len1 aber ist x1_len=0, sum1_len = k_lo, also
 242           //  war prodmid_len = sum1_len+sum2_len <= 2*k_lo+1 <= len2+2,
 243           //  nach 2-maliger Verkleinerung jedenfalls prodmid_len <= len2.)
 244           // prodmid*b^k = (x0*y1+x1*y0)*b^k zu prod = x1*y1*b^(2*k) + x0*y0 addieren:
 245           {var uintD carry =
 246              addto_loop_lsp(prodmid_LSDptr,prod_LSDptr lspop k_lo,prodmid_len);
 247            if (!(carry==0))
 248              { inc_loop_lsp(prod_LSDptr lspop (k_lo+prodmid_len),prod_len-(k_lo+prodmid_len)); }
 249       }}}}}
 250       // Das Teilprodukt zum Gesamtprodukt addieren:
 251       if (first_part)
 252         { copy_loop_lsp(prod_LSDptr,destptr,prod_len); }
 253         else
 254         { var uintD carry =
 255             addto_loop_lsp(prod_LSDptr,destptr,len1);
 256           destptr = destptr lspop len1;
 257           copy_loop_lsp(prod_LSDptr lspop len1,destptr,len2);
 258           if (!(carry==0)) { inc_loop_lsp(destptr,len2); }
 259         }
 260     }}