src/base/digitseq/cl_DS_mul.cc

   1 // cl_UDS_mul().
   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_low.h"
  13 #include "cl_malloc.h"
  14 #include "cl_abort.h"
  15
  16 // Multiplikations-Doppelschleife:
  17 // Multipliziert zwei UDS und legt das Ergebnis in einer dritten UDS ab.
  18 // cl_UDS_mul(sourceptr1,len1,sourceptr2,len2,destptr);
  19 // multipliziert die UDS  sourceptr1[-len1..-1]  (len1>0)
  20 //           mit der UDS  sourceptr2[-len1..-1]  (len2>0)
  21 // und legt das Ergebnis in der UDS  destptr[-len..-1]  (len=len1+len2) ab.
  22 // Unterhalb von destptr werden len Digits Platz benötigt.
  23   void cl_UDS_mul (const uintD* sourceptr1, uintC len1,
  24                    const uintD* sourceptr2, uintC len2,
  25                    uintD* destptr);
  26 // Spezialfall sourceptr1 == sourceptr2 && len1 == len2.
  27   void cl_UDS_mul_square (const uintD* sourceptr, uintC len,
  28                           uintD* destptr);
  29
  30 // Multiplikation nach Schulmethode:
  31   static inline void mulu_2loop (const uintD* sourceptr1, uintC len1,
  32                                  const uintD* sourceptr2, uintC len2,
  33                                  uintD* destptr)
  34   { // Es ist 2 <= len1 <= len2.
  35     // Erster Schleifendurchlauf:
  36     mulu_loop_lsp(lsprefnext(sourceptr1),sourceptr2,destptr,len2);
  37     lsshrink(destptr);
  38     var uintD* destptr2 = destptr lspop len2;
  39     // äußere Schleife läuft über source1 :
  40     dotimespC(len1,len1-1,
  41       { // innere Schleife läuft über source2 :
  42         var uintD carry =
  43           muluadd_loop_lsp(lsprefnext(sourceptr1),sourceptr2,destptr,len2);
  44         lsprefnext(destptr2) = carry; // UDS um das Carry-Digit verlängern
  45         lsshrink(destptr);
  46       });
  47   }
  48   static inline void mulu_2loop_square (const uintD* sourceptr, uintC len,
  49                                         uintD* destptr)
  50   { // Es ist 2 <= len.
  51     // Gemischte Produkte:
  52     #if 0
  53     // 2*(  x[1] * x[0..0] * b^1
  54     //    + x[2] * x[1..0] * b^2
  55     //    + ...
  56     //    + x[n-1] * x[n-2..0]*b^(n-1))
  57     { var const uintD* sourceptr1 = sourceptr lspop 1;
  58       var uintD* destptr1 = destptr;
  59       lsprefnext(destptr1) = 0;
  60       var uintD* destptr2 = destptr1;
  61       var uintC count;
  62       for (count = 1; count < len; count++)
  63         { // sourceptr1 = sourceptr lspop count, destptr1 = destptr lspop count,
  64           // destptr2 = destptr lspop (2*count-1).
  65           lsprefnext(destptr2) = 0;
  66           var uintD carry =
  67             muluadd_loop_lsp(lsprefnext(sourceptr1),sourceptr,destptr1,count);
  68           lsprefnext(destptr2) = carry;
  69           destptr1 = destptr1 lspop 1;
  70         }
  71       { var uintD carry = shift1left_loop_lsp(destptr lspop 1,2*len-2);
  72         lspref(destptr2,0) = (carry==0 ? 0 : 1);
  73     } }
  74     #else
  75     // 2*(  x[n-1..1] * x[0] * b^1
  76     //    + x[n-1..2] * x[1] * b^3
  77     //    + ...
  78     //    + x[n-1..n-1] * x[n-2] * b^(2*n-3))
  79     { var const uintD* sourceptr1 = sourceptr;
  80       var uintD* destptr2 = destptr;
  81       lsprefnext(destptr2) = 0;
  82       var uintC count = len-1;
  83       { var uintD digit = lsprefnext(sourceptr1);
  84         mulu_loop_lsp(digit,sourceptr1,destptr2,count);
  85       }
  86       var uintD* destptr1 = destptr lspop (len+1);
  87       while (--count > 0)
  88         { destptr2 = destptr2 lspop 2;
  89           var uintD digit = lsprefnext(sourceptr1);
  90           var uintD carry = muluadd_loop_lsp(digit,sourceptr1,destptr2,count);
  91           lsprefnext(destptr1) = carry;
  92         }
  93       { var uintD carry = shift1left_loop_lsp(destptr lspop 1,2*len-2);
  94         lspref(destptr1,0) = (carry==0 ? 0 : 1);
  95     } }
  96     #endif
  97     // Quadrate:
  98     len = 2*len;
  99     do { len -= 2;
 100          var uintD digit = lsprefnext(sourceptr);
 101          #if HAVE_DD
 102          var uintDD prod = muluD(digit,digit);
 103          var uintDD accu = highlowDD(lspref(destptr,1),lspref(destptr,0));
 104          accu += prod;
 105          lspref(destptr,0) = lowD(accu); lspref(destptr,1) = highD(accu);
 106          destptr = destptr lspop 2;
 107          if (accu < prod) { inc_loop_lsp(destptr,len); }
 108          #else
 109          var uintD hi;
 110          var uintD lo;
 111          muluD(digit,digit, hi=,lo=);
 112          var uintD tmp;
 113          tmp = lspref(destptr,0) + lo; lspref(destptr,0) = tmp;
 114          if (tmp < lo) hi++;
 115          tmp = lspref(destptr,1) + hi; lspref(destptr,1) = tmp;
 116          destptr = destptr lspop 2;
 117          if (tmp < hi) { inc_loop_lsp(destptr,len); }
 118          #endif
 119        } while (len > 0);
 120   }
 121
 122 // Karatsuba-multiplication: O(n^(log 3 / log 2))
 123   static void mulu_karatsuba (const uintD* sourceptr1, uintC len1,
 124                               const uintD* sourceptr2, uintC len2,
 125                               uintD* destptr);
 126   static void mulu_karatsuba_square (const uintD* sourceptr, uintC len,
 127                                      uintD* destptr);
 128 #include "cl_DS_mul_kara.h"
 129   // karatsuba_threshold = length, from which on Karatsuba-multiplication is a
 130   // gain and will be preferred.  The break-even point is determined from
 131   // timings.  The test is (progn (time (! 5000)) nil), which does many small
 132   // and some very large multiplications.  The measured runtimes are:
 133   // OS: Linux 2.2, intDsize==32,        OS: TRU64/4.0, intDsize==64,
 134   // Machine: P-III/450MHz               Machine: EV5/300MHz:
 135   // threshold  time in 0.01 sec.        time in 0.01 sec.
 136   //      5          3.55                     2.29
 137   //     10          2.01                     1.71
 138   //     15          1.61                     1.61
 139   //     20          1.51                     1.60  <-
 140   //     25          1.45                     1.63
 141   //     30          1.39                     1.66
 142   //     35          1.39  <-                 1.67
 143   //     40          1.39                     1.71
 144   //     45          1.40                     1.75
 145   //     50          1.41                     1.78
 146   //     55          1.41                     1.79
 147   //     60          1.44                     1.84
 148   //     65          1.44                     1.85
 149   //     70          1.43                     1.85
 150   //     75          1.45                     1.89
 151   //     80          1.47                     1.91
 152   //     90          1.51                     1.96
 153   //    100          1.53                     1.97
 154   //    150          1.62                     2.13
 155   //    250          1.75                     2.19
 156   //    500          1.87                     2.17
 157   //   1000          1.87                     2.18
 158   //   2000          1.88                     2.17
 159   // The optimum appears to be between 20 and 40.  But since that optimum
 160   // depends on the ratio time(uintD-mul)/time(uintD-add) and the measured
 161   // times are more sensitive to a shift towards lower thresholds we are
 162   // careful and choose a value at the upper end:
 163 #if CL_USE_GMP
 164   const unsigned int cl_karatsuba_threshold = 35;
 165 #else
 166   const unsigned int cl_karatsuba_threshold = 16;
 167   // (In CLN version <= 1.0.3 cl_karatsuba_threshold was always 16)
 168 #endif
 169
 170 #if 0 // Doesn't seem to be worth the effort
 171
 172 // FFT-Multiplikation nach Nussbaumer: O(n log n log log n)
 173 #include "cl_DS_mul_nuss.h"
 174   // nuss_threshold = Länge, ab der die Nussbaumer-Multiplikation bevorzugt
 175   // wird. Der Break-Even-Point bestimmt sich aus Zeitmessungen.
 176   // Multiplikation zweier N-Wort-Zahlen unter Linux mit einem 80486:
 177   //    N     kara   nuss  nuss-asm  (time in sec.)
 178   //   1000    0.36   1.05   0.70
 179   //   5000    4.69  10.0    6.71
 180   //  25000   61.6   62.7   40.2
 181   //  32500   91.8   62.7   40.3
 182   //  35000  102.7  124.7   80.4
 183   //  50000  185    132     85.2
 184   int cl_nuss_threshold = 1000000;
 185
 186 // FFT-Multiplikation in Z/pZ: O(n log n log log n)
 187 #include "cl_DS_mul_fftp.h"
 188   // fftp_threshold = Länge, ab der die FFT-Multiplikation mod p bevorzugt
 189   // wird. Der Break-Even-Point bestimmt sich aus Zeitmessungen.
 190   // Multiplikation zweier N-Wort-Zahlen unter Linux mit einem 80486:
 191   //    N     kara   fftp  (time in sec.)
 192   //   1000    0.36   1.57
 193   //   5000    4.66  14.86
 194   //  25000   61.1   75.0
 195   //  32500   90.8   75.5
 196   //  35000  101.6  150.1
 197   //  50000  183    160
 198   int cl_fftp_threshold = 1000000;
 199
 200 // FFT-Multiplikation in Z/pZ: O(n log n log log n)
 201 // für drei verschiedene Primzahlen p1,p2,p3 < 2^32.
 202 #include "cl_DS_mul_fftp3.h"
 203   // fftp3_threshold = Länge, ab der die FFT-Multiplikation mod p_i bevorzugt
 204   // wird. Der Break-Even-Point bestimmt sich aus Zeitmessungen.
 205   // Multiplikation zweier N-Wort-Zahlen unter Linux mit einem 80486:
 206   //    N     kara   fftp3  fftp  (time in sec.)
 207   //   1000    0.36   0.59   1.57
 208   //   5000    4.66   5.44  14.89
 209   //  10000   13.98  11.91  32.43
 210   //  25000   61.1   27.4   75.4
 211   //  32500   90.5   28.1   75.5
 212   //  35000  101.4   54.8  150.4
 213   //  50000  183     58.9  161.6
 214   int cl_fftp3_threshold = 1000000;
 215
 216 // FFT-Multiplikation in Z/pZ: O(n log n log log n)
 217 // für drei verschiedene Primzahlen p1,p2,p3 < 2^32,
 218 // mit Montgomery-Multiplikation.
 219 #include "cl_DS_mul_fftp3m.h"
 220   // fftp3_threshold = Länge, ab der die FFT-Multiplikation mod p_i bevorzugt
 221   // wird. Der Break-Even-Point bestimmt sich aus Zeitmessungen.
 222   // Multiplikation zweier N-Wort-Zahlen unter
 223   // Linux mit einem 80486, 33 MHz, mit Benutzung der GMP-Low-Level-Funktionen:
 224   //    N     kara    fftm  fftp3m  fftp3  fftp  (time in sec.)
 225   //   1000    0.35   0.49   0.54   0.59   1.58
 226   //   2500    1.48   0.97   2.34   2.52   6.99
 227   //   5000    4.43   2.19   5.08   5.48  15.16
 228   //  10000   13.33   4.68  10.93  11.82  32.94
 229   //  25000   58.5   12.0   25.3   27.4   77.0
 230   //  32500   86.0   25.0   26.1   28.0   77.3
 231   //  35000   96.5   25.0   50.8   54.9  152.8
 232   //  50000  176     25.2   54.2   58.5  163.4
 233   // und auf einer SPARC 20 mit 75 MHz, ohne GMP-Low-Level-Funktionen:
 234   //    N     kara    fftm  fftp3m  fftp3  fftp  (time in sec.)
 235   //   1000    0.076  0.096  0.113  0.233  0.415
 236   //   2500    0.32   0.21   0.48   1.03   1.82
 237   //   5000    0.97   0.51   1.03   2.22   3.96
 238   //  10000    2.99   1.03   2.23   4.72   8.59
 239   //  25000   13.22   2.73   4.99  10.78  19.73
 240   //  32500   19.3    5.7    5.2   10.9   19.7
 241   //  35000   21.5    5.9   10.0   21.7   39.4
 242   //  50000   39.5    6.0   11.3   23.1   42.7
 243   int cl_fftp3m_threshold = 1000000;
 244
 245 #endif
 246
 247 // FFT-Multiplikation in Z/pZ: O(n^1.29)
 248 #include "cl_DS_mul_fftm.h"
 249   // fftm_threshold = length, from which on FFT multiplication mod m is a gain
 250   // and will be preferred.  The break-even point is determined from timings.
 251   // The times to multiply two N-limb numbers are:
 252   // OS: Linux 2.2, intDsize==32,        OS: TRU64/4.0, intDsize==64,
 253   // Machine: P-III/450MHz               Machine: EV5/300MHz:
 254   //    N     kara   fftm  (time in sec.)    kara   fftm
 255   //   1000    0.005  0.016                   0.018  0.028
 256   //   1500    0.009  0.012                   0.032  0.028
 257   //   2000    0.015  0.025                   0.053  0.052  <-
 258   //   2500    0.022  0.026                   0.067  0.052
 259   //   3000    0.029  0.027  <-               0.093  0.053
 260   //   3500    0.035  0.037                   0.12   0.031
 261   //   4000    0.045  0.028                   0.16   0.12
 262   //   5000    0.064  0.050                   0.20   0.11
 263   //   7000    0.110  0.051                   0.37   0.20
 264   //  10000    0.19   0.11                    0.61   0.26
 265   //  20000    0.59   0.23                    1.85   0.55
 266   //  30000    1.10   0.25                    3.79   0.56
 267   //  50000    2.52   1.76                    8.15   1.37
 268   //  70000    4.41   2.30                   14.09   2.94
 269   // 100000    7.55   1.53                   24.48   2.96
 270   // More playing around with timings reveals that there are some values where
 271   // FFT multiplication is somewhat slower than Karatsuba, both for len1==len2
 272   // and also if len1<len2.
 273   // Here are the timigs from CLN version <= 1.0.3:
 274   // //  Linux mit einem 80486:               Solaris, Sparc 10/20:
 275   // //    N     kara   fftm  (time in sec.)    kara   fftm
 276   // //   1000    0.36   0.54                    0.08   0.10
 277   // //   5000    4.66   2.48                    1.01   0.51
 278   // //  25000   61.1   13.22                   13.23   2.73
 279   // //  32500   91.0   27.5                    20.0    5.8
 280   // //  35000  102.1   27.5                    21.5    5.6
 281   // //  50000  183     27.6                    40.7    5.6
 282   // // Multiplikation zweier N-Wort-Zahlen unter
 283   // //  Linux mit einem 80486:               Solaris, Sparc 10/20:
 284   // //    N     kara   fftm  (time in sec.)    kara   fftm
 285   // //   1000    0.36   0.54                    0.08   0.10
 286   // //   1260    0.52   0.50                    0.11   0.10
 287   // //   1590    0.79   0.51                    0.16   0.10
 288   // //   2000    1.09   1.07                    0.23   0.21
 289   // //   2520    1.57   1.08                    0.33   0.21
 290   // //   3180    2.32   1.08                    0.50   0.21
 291   // //   4000    3.29   2.22                    0.70   0.41
 292   // //   5040    4.74   2.44                    0.99   0.50
 293   // //    N1    N2    kara   fftm  (time in sec.)    kara   fftm
 294   // //   1250  1250    0.51   0.50                    0.11   0.10
 295   // //   1250  1580    0.70   0.50                    0.15   0.10
 296   // //   1250  2000    0.89   0.51                    0.18   0.10
 297   // //   1250  2250    0.99   0.51                    0.21   0.10
 298   // //   1250  2500    1.08   1.03     <---           0.22   0.21
 299   // //   1250  2800    1.20   1.07                    0.26   0.21
 300   // //   1250  3100    1.35   1.07                    0.28   0.21
 301   // // Es gibt also noch Werte von (len1,len2) mit 1250 <= len1 <= len2, bei
 302   // // denen "kara" schneller ist als "fftm", aber nicht viele und dort auch
 303   // // nur um 5%. Darum wählen wir ab hier die FFT-Multiplikation.
 304   // // 140000: 4.15s  12.53  23.7
 305   // // 14000:  4.16s
 306   // // 11000:  4.16s
 307   // // 9000:   1.47s
 308   // // 7000:   1.48s
 309   // // 1400:   1.42s   2.80   6.5
 310 #if CL_USE_GMP
 311   const unsigned int cl_fftm_threshold = 2500;
 312   // must be >= 6 (else infinite recursion)
 313 #else
 314   // Use the old default value from CLN version <= 1.0.3 as a crude estimate.
 315   const unsigned int cl_fftm_threshold = 1250;
 316 #endif
 317   // This is the threshold for multiplication of equally sized factors.
 318   // When the lengths differ much, the threshold varies:
 319   //                OS: Linux 2.2, intDsize==32,  OS: TRU64/4.0, intDsize==64,
 320   //                Machine: P-III/450MHz         Machine: EV5/300MHz:
 321   // len2 =  3000   len1 >= 2600                  len1 >= 800
 322   // len2 =  4000   len1 >= 1500                  len1 >= 700
 323   // len2 =  5000   len1 >= 1100                  len1 >= 600
 324   // len2 =  6000   len1 >= 1300                  len1 >= 700
 325   // len2 =  7000   len1 >= 1100                  len1 >= 600
 326   // len2 =  8000   len1 >= 900                   len1 >= 500
 327   // len2 =  9000   len1 >= 1300                  len1 >= 600
 328   // len2 = 10000   len1 >= 1100                  len1 >= 500
 329   // len2 = 11000   len1 >= 1000                  len1 >= 500
 330   // len2 = 12000   len1 >= 900                   len1 >= 700
 331   // len2 = 13000   len1 >= 900                   len1 >= 500
 332   // len2 = 14000   len1 >= 900                   len1 >= 600
 333   // Here are the timigs from CLN version <= 1.0.3:
 334   // //   len2 = 3000   len1 >= 800
 335   // //   len2 = 3500   len1 >= 700
 336   // //   len2 = 4000   len1 >= 580
 337   // //   len2 = 4500   len1 >= 430
 338   // //   len2 = 5000   len1 >= 370
 339   // //   len2 = 5500   len1 >= 320
 340   // //   len2 = 6000   len1 >= 500
 341   // //   len2 = 7000   len1 >= 370
 342   // //   len2 = 8000   len1 >= 330
 343   // //   len2 = 9000   len1 >= 420
 344   // //   len2 =10000   len1 >= 370
 345   // //   len2 =11000   len1 >= 330
 346   // //   len2 =12000   len1 >= 330
 347   // //   len2 =13000   len1 >= 350
 348   // Let's choose the following condition:
 349 #if CL_USE_GMP
 350   const unsigned int cl_fftm_threshold1 = 600;
 351 #else
 352   // Use the old default values from CLN version <= 1.0.3 as a crude estimate.
 353   const unsigned int cl_fftm_threshold1 = 330;
 354 #endif
 355   const unsigned int cl_fftm_threshold2 = 2*cl_fftm_threshold;
 356   //   len1 > cl_fftm_threshold1 && len2 > cl_fftm_threshold2
 357   //   && len1 >= cl_fftm_threshold1 + cl_fftm_threshold/(len2-cl_fftm_threshold1)*(cl_fftm_threshold-cl_fftm_threshold1).
 358   static inline cl_boolean cl_fftm_suitable (uintL len1, uintL len2)
 359     { if (len1 >= cl_fftm_threshold)
 360         return cl_true;
 361       if (len1 > cl_fftm_threshold1)
 362         if (len2 > cl_fftm_threshold2)
 363           { var uint32 hi;
 364             var uint32 lo;
 365             mulu32(len1-cl_fftm_threshold1,len2-cl_fftm_threshold1, hi=,lo=);
 366             if (hi > 0 || lo >= cl_fftm_threshold*(cl_fftm_threshold-cl_fftm_threshold1))
 367               return cl_true;
 368           }
 369       return cl_false;
 370     }
 371
 372 #if 0 // Doesn't seem to be worth the effort
 373
 374 // FFT-Multiplikation über den komplexen Zahlen.
 375 #include "cl_DS_mul_fftc.h"
 376   // Multiplikation zweier N-Wort-Zahlen unter
 377   //  Linux mit einem i486 33 MHz
 378   //    N     kara/fftm  fftc   fftclong
 379   //   1000      0.35     1.52     0.94
 380   //   2500      0.98     7.6/8.4  4.7
 381   //   5000      2.2     18.2     10.2
 382   //  10000      4.7     34       22
 383   // Multiplikation zweier N-Wort-Zahlen unter
 384   //  Linux mit einem i586 90/100 MHz
 385   //    N     kara/fftm  fftc   fftclong
 386   //   1000      0.03     0.20   0.16
 387   //   2500      0.16     1.6    0.92
 388   //   5000      0.3      2.6    2.2
 389   //  10000      0.7      7.1    4.8
 390   //  25000      1.6    (50MB)  20.7(22MB)
 391   // Multiplikation zweier N-Wort-Zahlen unter
 392   //  Solaris, Sparc 20, 75 MHz
 393   //    N     kara/fftm  fftc
 394   //   1000      0.07     0.14
 395   //   2500      0.21     0.76
 396   //   5000      0.44     1.75
 397   //  10000      0.88     4.95
 398   //  25000      2.3     (15MB)
 399
 400 // FFT-Multiplikation über den komplexen Zahlen, Symmetrie ausnutzend.
 401 #include "cl_DS_mul_fftcs.h"
 402   // Multiplikation zweier N-Wort-Zahlen unter
 403   //  Linux mit einem i486 33 MHz
 404   //    N     kara/fftm  fftcs  fftcslong
 405   //   1000      0.34     0.71     0.43
 406   //   2500      0.98     3.4      2.1
 407   //   5000      2.2      8.0      4.7
 408   //  10000      4.7     16.1     10.4
 409   // Multiplikation zweier N-Wort-Zahlen unter
 410   //  Solaris, Sparc 20, 75 MHz
 411   //    N     kara/fftm  fftcs
 412   //    300      0.010    0.012
 413   //    400      0.018    0.027
 414   //    500      0.023    0.027
 415   //    600      0.031    0.027
 416   //    700      0.031    0.027
 417   //    800      0.051    0.058
 418   //    900      0.064    0.059
 419   //   1000      0.069    0.059
 420   //   1250      0.088    0.13
 421   //   1500      0.088    0.13
 422   //   1750      0.088    0.13
 423   //   2000      0.19     0.13
 424   //   2500      0.19     0.29
 425   //   3000      0.19     0.33
 426   //   3500      0.20     0.31
 427   //   4000      0.37     0.70
 428   //   4500      0.38     0.70
 429   //   5000      0.43     0.69
 430   //   6000      0.43     0.69
 431   //   7000      0.43     1.62
 432   //   8000      0.88     1.60
 433   //   9000      0.88     1.6
 434   //  10000      0.90     1.55
 435   //  12000      0.89     4.7
 436   //  14000      0.90     5.2
 437   //  16000      1.43     5.2
 438
 439 #endif
 440
 441 #if 0 // Keine gute Fehlerabschätzung
 442
 443 // FFT-Multiplikation über den komplexen Zahlen, mit reellen Zahlen rechnend.
 444 #include "cl_DS_mul_fftr.h"
 445   // Multiplikation zweier N-Wort-Zahlen unter
 446   //  Linux mit einem i486 33 MHz
 447   //    N     kara/fftm  fftr   fftrlong
 448   //   1000      0.34     0.64     0.40
 449   //   2500      0.98     3.5      2.0
 450   //   5000      2.2      7.2/7.7  4.6
 451   //  10000      4.7     16.6     10.0
 452
 453 #endif
 454
 455 //  int cl_mul_algo = 0;
 456   void cl_UDS_mul (const uintD* sourceptr1, uintC len1,
 457                    const uintD* sourceptr2, uintC len2,
 458                    uintD* destptr)
 459     { // len1<=len2 erzwingen:
 460       if (len1>len2)
 461         {{var const uintD* temp;
 462           temp = sourceptr1; sourceptr1 = sourceptr2; sourceptr2 = temp;
 463          }
 464          {var uintC temp;
 465           temp = len1; len1 = len2; len2 = temp;
 466         }}
 467       if (len1==1)
 468         // nur eine Einfachschleife
 469         { mulu_loop_lsp(lsprefnext(sourceptr1),sourceptr2,destptr,len2); }
 470       else
 471         {
 472 //          if (cl_mul_algo > 0)
 473 //            mulu_fftcs(sourceptr1,len1,sourceptr2,len2,destptr);
 474 //          else
 475 //          if (cl_mul_algo > 0)
 476 //            mulu_nussbaumer(sourceptr1,len1,sourceptr2,len2,destptr);
 477 //          else
 478           if (len1 < cl_karatsuba_threshold)
 479             // Multiplikation nach Schulmethode
 480             mulu_2loop(sourceptr1,len1,sourceptr2,len2,destptr);
 481           else // len1 groß
 482           if (!cl_fftm_suitable(len1,len2))
 483             // Karatsuba-Multiplikation
 484             // (ausgelagert, um die eigentliche Multiplikationsfunktion nicht
 485             // durch zu viele Registervariablen zu belasten):
 486             mulu_karatsuba(sourceptr1,len1,sourceptr2,len2,destptr);
 487           else
 488             //mulu_fft_modp(sourceptr1,len1,sourceptr2,len2,destptr);
 489             //mulu_nussbaumer(sourceptr1,len1,sourceptr2,len2,destptr);
 490             //mulu_fft_modp3(sourceptr1,len1,sourceptr2,len2,destptr);
 491             mulu_fft_modm(sourceptr1,len1,sourceptr2,len2,destptr);
 492           #ifdef DEBUG_MUL_XXX
 493           { // Check the correctness of an other multiplication algorithm:
 494             CL_ALLOCA_STACK;
 495             var uintD tmpprod_xxx = cl_alloc_array(uintD,len1+len2);
 496             mulu_xxx(sourceptr1,len1,sourceptr2,len2,arrayLSDptr(tmpprod_xxx,len1+len2));
 497             if (compare_loop_msp(destptr lspop (len1+len2),arrayMSDptr(tmpprod_xxx,len1+len2),len1+len2))
 498               cl_abort();
 499           }
 500           #endif
 501         }
 502     }
 503
 504 // Special support for squaring.
 505 // Squaring takes approximately 69% of the time of a normal multiplication.
 506   #include "cl_DS_mul_kara_sqr.h" // defines mulu_karatsuba_square()
 507   void cl_UDS_mul_square (const uintD* sourceptr, uintC len,
 508                           uintD* destptr)
 509   { if (len==1)
 510       { var uintD digit = lspref(sourceptr,0);
 511         #if HAVE_DD
 512         var uintDD prod = muluD(digit,digit);
 513         lspref(destptr,0) = lowD(prod); lspref(destptr,1) = highD(prod);
 514         #else
 515         muluD(digit,digit, lspref(destptr,1)=,lspref(destptr,0)=);
 516         #endif
 517       }
 518     else
 519       { if (len < cl_karatsuba_threshold)
 520           mulu_2loop_square(sourceptr,len,destptr);
 521         else
 522           if (!(len >= cl_fftm_threshold))
 523             mulu_karatsuba_square(sourceptr,len,destptr);
 524           else
 525             mulu_fft_modm(sourceptr,len,sourceptr,len,destptr);
 526       }
 527   }