/cln-1.3.2/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    } }}}