src/base/digitseq/cl_DS_mul_kara_sqr.h

   1   // Eine vereinfachte Version von mulu_karatsuba für den Fall
   2   // sourceptr1 == sourceptr2 && len1 == len2.
   3   // Weniger Variablen, eine Additionsschleife weniger, eine Kopierschleife
   4   // weniger, und bei rekursiven Aufrufen ist wieder
   5   // sourceptr1 == sourceptr2 && len1 == len2.
   6   static void mulu_karatsuba_square (const uintD* sourceptr, uintC len,
   7                                      uintD* destptr)
   8     { // Es ist 2 <= len.
   9       CL_SMALL_ALLOCA_STACK;
  10       var uintC prod_len = 2*len;
  11       var uintD* prod_LSDptr = destptr;
  12       var uintC k_hi = floor(len,2); // Länge der High-Teile: floor(len/2) >0
  13       var uintC k_lo = len - k_hi; // Länge der Low-Teile: ceiling(len/2) >0
  14       // Es gilt k_hi <= k_lo <= len, k_lo + k_hi = len.
  15       // Summe x1+x0 berechnen:
  16       var uintD* sum_MSDptr;
  17       var uintC sum_len = k_lo; // = max(k_lo,k_hi)
  18       var uintD* sum_LSDptr;
  19       num_stack_small_alloc_1(sum_len,sum_MSDptr=,sum_LSDptr=);
  20       {var uintD carry = // Hauptteile von x1 und x0 addieren:
  21          add_loop_lsp(sourceptr lspop k_lo,sourceptr,sum_LSDptr,k_hi);
  22        if (!(k_lo==k_hi))
  23          // noch k_lo-k_hi = 1 Digits abzulegen
  24          { mspref(sum_MSDptr,0) = lspref(sourceptr,k_lo-1); // = lspref(sourceptr,k_hi)
  25            if (!(carry==0)) { if (++(mspref(sum_MSDptr,0)) == 0) carry=1; else carry=0; }
  26          }
  27        if (carry) { lsprefnext(sum_MSDptr) = 1; sum_len++; }
  28       }
  29       // Platz für Produkte x0*x0, x1*x1:
  30       { var uintC prodhi_len = 2*k_hi;
  31         var uintD* prodhi_LSDptr = prod_LSDptr lspop 2*k_lo;
  32         // prod_MSDptr/2*len/prod_LSDptr wird zuerst die beiden
  33         // Produkte x1*x1 in prod_MSDptr/2*k_hi/prodhi_LSDptr
  34         //      und x0*x0 in prodhi_LSDptr/2*k_lo/prod_LSDptr,
  35         // dann das Produkt (b^k*x1+x0)*(b^k*x1+x0) enthalten.
  36         // Platz fürs Produkt (x1+x0)*(x1+x0) belegen:
  37        {var uintD* prodmid_MSDptr;
  38         var uintC prodmid_len = 2*sum_len;
  39         var uintD* prodmid_LSDptr;
  40         num_stack_small_alloc(prodmid_len,prodmid_MSDptr=,prodmid_LSDptr=);
  41         // Produkt (x1+x0)*(x1+x0) berechnen:
  42         cl_UDS_mul_square(sum_LSDptr,sum_len,prodmid_LSDptr);
  43         // Das Produkt beansprucht  2*k_lo + (0 oder 1) <= 2*sum_len = prodmid_len  Digits.
  44         // Produkt x0*x0 berechnen:
  45         cl_UDS_mul_square(sourceptr,k_lo,prod_LSDptr);
  46         // Produkt x1*x1 berechnen:
  47         cl_UDS_mul_square(sourceptr lspop k_lo,k_hi,prodhi_LSDptr);
  48         // Und x1*x1 abziehen:
  49         {var uintD carry =
  50            subfrom_loop_lsp(prodhi_LSDptr,prodmid_LSDptr,prodhi_len);
  51          // Carry um maximal prodmid_len-prodhi_len Digits weitertragen:
  52          if (!(carry==0))
  53            { dec_loop_lsp(prodmid_LSDptr lspop prodhi_len,prodmid_len-prodhi_len); }
  54         }
  55         // Und x0*x0 abziehen:
  56         {var uintD carry =
  57            subfrom_loop_lsp(prod_LSDptr,prodmid_LSDptr,2*k_lo);
  58          // Falls Carry: Produkt beansprucht 2*k_lo+1 Digits.
  59          // Carry um maximal 1 Digit weitertragen:
  60          if (!(carry==0)) { lspref(prodmid_LSDptr,2*k_lo) -= 1; }
  61         }
  62         // prodmid_LSDptr[-prodmid_len..-1] enthält nun 2*x0*x1.
  63         // Dies ist < 2 * b^k_lo * b^k_hi = 2 * b^len,
  64         // paßt also in len+1 Digits.
  65         // prodmid_len, wenn möglich, um maximal 2 verkleinern:
  66         // (benutzt prodmid_len >= 2*k_lo >= len >= 2)
  67         if (mspref(prodmid_MSDptr,0)==0)
  68           { prodmid_len--;
  69             if (mspref(prodmid_MSDptr,1)==0) { prodmid_len--; }
  70           }
  71         // Nun ist k_lo+prodmid_len <= 2*len .
  72         // (Denn es war prodmid_len = 2*sum_len <= 2*(k_lo+1)
  73         //  <= len+3, und nach 2-maliger Verkleinerung jedenfalls
  74         //  prodmid_len <= len+1. Wegen k_lo < len also
  75         //  k_lo + prodmid_len <= (len-1)+(len+1) = 2*len.)
  76         // prodmid*b^k = 2*x0*x1*b^k zu prod = x1*x1*b^(2*k) + x0*x0 addieren:
  77         {var uintD carry =
  78            addto_loop_lsp(prodmid_LSDptr,prod_LSDptr lspop k_lo,prodmid_len);
  79          if (!(carry==0))
  80            { inc_loop_lsp(prod_LSDptr lspop (k_lo+prodmid_len),prod_len-(k_lo+prodmid_len)); }
  81     } }}}