home *** CD-ROM | disk | FTP | other *** search
/ Il CD di internet / CD.iso / SOURCE / D / CLISP / CLISPSRC.TAR / clisp-1995-01-01 / src / intmal.d < prev    next >
Encoding:
Text File  |  1994-08-12  |  28.3 KB  |  625 lines
  1. # Multiplikation ganzer Zahlen
  2.  
  3. # meldet ▄berlauf bei der Multiplikation:
  4.   nonreturning_function(local, mal_ueberlauf, (void));
  5.   local void mal_ueberlauf()
  6.     { fehler(error,
  7.              DEUTSCH ? "▄berlauf bei Multiplikation langer Zahlen" :
  8.              ENGLISH ? "overflow during multiplication of large numbers" :
  9.              FRANCAIS ? "DΘbordement de capacitΘ lors d'une multiplication de grands nombres." :
  10.              ""
  11.             );
  12.     }
  13.  
  14. # Multipliziert zwei Unsigned-Digit-sequences.
  15. # UDS_UDS_mal_UDS(len1,LSDptr1, len2,LSDptr2, MSDptr=,len=,LSDptr=);
  16. # multipliziert die UDS ../len1/LSDptr1 und ../len2/LSDptr2.
  17. # Dabei sollte len1>0 und len2>0 sein.
  18. # Ergebnis ist die UDS MSDptr/len/LSDptr, mit len=len1+len2, im Stack.
  19. # Dabei wird num_stack erniedrigt.
  20.   #define UDS_UDS_mal_UDS(len1,LSDptr1,len2,LSDptr2, MSDptr_zuweisung,len_zuweisung,LSDptr_zuweisung)  \
  21.     { var reg1 uintL len_from_UDSmal = (uintL)(len1) + (uintL)(len2);                            \
  22.       var reg2 uintD* LSDptr_from_UDSmal;                                                        \
  23.       if ((intCsize < 32) && (len_from_UDSmal > (uintL)(bitc(intCsize)-1))) { mal_ueberlauf(); } \
  24.       nowarn len_zuweisung len_from_UDSmal;                                                      \
  25.       num_stack_need(len_from_UDSmal,MSDptr_zuweisung,LSDptr_zuweisung LSDptr_from_UDSmal =);    \
  26.       mulu_2loop_down((LSDptr1),(len1),(LSDptr2),(len2),LSDptr_from_UDSmal);                     \
  27.     }
  28.  
  29. # Multiplikations-Doppelschleife:
  30. # Multipliziert zwei UDS und legt das Ergebnis in einer dritten UDS ab.
  31. # mulu_2loop_down(sourceptr1,len1,sourceptr2,len2,destptr);
  32. # multipliziert die UDS  sourceptr1[-len1..-1]  (len1>0)
  33. #           mit der UDS  sourceptr2[-len1..-1]  (len2>0)
  34. # und legt das Ergebnis in der UDS  destptr[-len..-1]  (len=len1+len2) ab.
  35. # Unterhalb von destptr werden len Digits Platz ben÷tigt.
  36.   local void mulu_2loop_down (uintD* sourceptr1, uintC len1,
  37.                               uintD* sourceptr2, uintC len2,
  38.                               uintD* destptr);
  39.   local void mulu_2bigloop_down (uintD* sourceptr1, uintC len1,
  40.                                  uintD* sourceptr2, uintC len2,
  41.                                  uintD* destptr);
  42.   # karatsuba_threshold = LΣnge, ab der die Karatsuba-Multiplikation bevorzugt
  43.   # wird. Der Break-Even-Point bestimmt sich aus Zeitmessungen.
  44.   # Als Test dient (progn (time (! 5000)) nil), das viele kleine und einige
  45.   # ganz gro▀e Multiplikationen durchfⁿhrt. Mi▀ die Runtime.
  46.   # Unter Linux mit einem 80486:      Auf einer Sparc 2:
  47.   # threshold  time in 0.01 sec.
  48.   #      5        125
  49.   #      6        116
  50.   #      7        107
  51.   #      8        101
  52.   #      9         99
  53.   #     10         98
  54.   #     11         97
  55.   #     12         96
  56.   #     13         97
  57.   #     14         97
  58.   #     15         97
  59.   #     16         98
  60.   #     17         98
  61.   #     18         98
  62.   #     19         98
  63.   #     20         99
  64.   #     25        103
  65.   #     30        109
  66.   #     40        115
  67.   #     50        122
  68.   #     70        132
  69.   #    100        151
  70.   #    150        164
  71.   #    250        183
  72.   #    500        203
  73.   #   1000        203
  74.   # Das Optimum scheint bei karatsuba_threshold = 12 zu liegen.
  75.   # Da das Optimum aber vom VerhΣltnis
  76.   #         Zeit fⁿr uintD-Multiplikation / Zeit fⁿr uintD-Addition
  77.   # abhΣngt und die gemessenen Zeiten auf eine Unterschreitung des Optimums
  78.   # empfindlicher reagieren als auf eine ▄berschreitung des Optimums,
  79.   # sind wir vorsichtig und wΣhlen einen Wert etwas ⁿber dem Optimum:
  80.   #define karatsuba_threshold  16
  81.   local void mulu_2loop_down(sourceptr1,len1,sourceptr2,len2,destptr)
  82.     var reg3 uintD* sourceptr1;
  83.     var reg4 uintC len1;
  84.     var reg6 uintD* sourceptr2;
  85.     var reg7 uintC len2;
  86.     var reg5 uintD* destptr;
  87.     { # len1<=len2 erzwingen:
  88.       if (len1>len2)
  89.         {{var reg1 uintD* temp;
  90.           temp = sourceptr1; sourceptr1 = sourceptr2; sourceptr2 = temp;
  91.          }
  92.          {var reg1 uintC temp;
  93.           temp = len1; len1 = len2; len2 = temp;
  94.         }}
  95.       if (len1==1)
  96.         # nur eine Einfachschleife
  97.         { mulu_loop_down(*--sourceptr1,sourceptr2,destptr,len2); }
  98.       elif (len1 < karatsuba_threshold)
  99.         # Multiplikation nach Schulmethode
  100.         { # len2 Digits auf 0 setzen:
  101.           var reg2 uintD* destptr2 = clear_loop_down(destptr,len2);
  102.           # Σu▀ere Schleife lΣuft ⁿber source1 :
  103.           dotimespC(len1,len1,
  104.             { # innere Schleife lΣuft ⁿber source2 :
  105.               var reg1 uintD carry =
  106.                 muluadd_loop_down(*--sourceptr1,sourceptr2,destptr,len2);
  107.               *--destptr2 = carry; # UDS um das Carry-Digit verlΣngern
  108.               destptr--;
  109.             });
  110.         }
  111.       else # len1 gro▀
  112.         # Karatsuba-Multiplikation
  113.         # (ausgelagert, um die eigentliche Multiplikationsfunktion nicht
  114.         # durch zu viele Registervariablen zu belasten):
  115.         mulu_2bigloop_down(sourceptr1,len1,sourceptr2,len2,destptr);
  116.     }
  117.   local void mulu_2bigloop_down(sourceptr1,len1,sourceptr2,len2,destptr)
  118.     var reg4 uintD* sourceptr1;
  119.     var reg5 uintC len1;
  120.     var reg6 uintD* sourceptr2;
  121.     var reg7 uintC len2;
  122.     var reg8 uintD* destptr;
  123.     # Karatsuba-Multiplikation
  124.     # Prinzip: (x1*b^k+x0) * (y1*b^k+y0)
  125.     #        = x1*y1 * b^2k + ((x1+x0)*(y1+y0)-x1*y1-x0*y0) * b^k + x0*y0
  126.     # Methode 1 (Collins/Loos, Degel):
  127.     # source2 wird in floor(len2/len1) einzelne UDS mit je einer
  128.     # LΣnge len3 (len1 <= len3 < 2*len1) unterteilt,
  129.     # jeweils k=floor(len3/2).
  130.     # Methode 2 (Haible):
  131.     # source2 wird in ceiling(len2/len1) einzelne UDS mit je einer
  132.     # LΣnge len3 (0 < len3 <= len1) unterteilt, jeweils k=floor(len1/2).
  133.     # Aufwand fⁿr die hinteren Einzelteile:
  134.     # bei beiden Methoden jeweils 3*len1^2.
  135.     # Aufwand fⁿr das vorderste Teil (alles, falls len1 <= len2 < 2*len1)
  136.     # mit r = len1, s = (len2 mod len1) + len1 (>= len1, < 2*len1):
  137.     # bei Methode 1:
  138.     #                       |   :       |  r
  139.     #                    |      :       |  s
  140.     #      (r-s/2)*s/2 + s/2*s/2 + s/2*s/2 = r*s/2 + s^2/4 .
  141.     # bei Methode 2:
  142.     #                       |     :     |  r
  143.     #                    |  |     :     |  s
  144.     #      (s-r)*r + r/2*r/2 + r/2*r/2 + r/2*r/2 = r*s - r^2/4 .
  145.     # Wegen (r*s/2 + s^2/4) - (r*s - r^2/4) = (r-s)^2/4 >= 0
  146.     # ist Methode 2 gⁿnstiger.
  147.     # Denkfehler! Dies gilt - wenn ⁿberhaupt - nur knapp oberhalb des
  148.     # Break-Even-Points.
  149.     # Im allgemeinen ist der Multiplikationsaufwand fⁿr zwei Zahlen der
  150.     # LΣngen u bzw. v nΣmlich gegeben durch  min(u,v)^c * max(u,v),
  151.     # wobei  c = log3/log2 - 1 = 0.585...
  152.     # Dadurch wird der Aufwand in AbhΣngigkeit des Parameters t = k,
  153.     # r/2 <= t <= s/2 (der einzig sinnvolle Bereich), zu
  154.     # (r-t)^c*(s-t) + t^c*(s-t) + t^(1+c).
  155.     # Dessen Optimum liegt (im Bereich r <= s <= 2*r)
  156.     # - im klassischen Fall c=1 tatsΣchlich stets bei t=r/2 [Methode 2],
  157.     # - im Karatsuba-Fall c=0.6 aber offenbar bei t=s/2 [Methode 1]
  158.     #   oder ganz knapp darunter.
  159.     # Auch erweist sich Methode 1 im Experiment als effizienter.
  160.     # Daher implementieren wir Methode 1 :
  161.     { # Es ist 2 <= len2 <= len1.
  162.       var reg9 boolean first_part = TRUE; # Flag, ob jetzt das erste Teilprodukt berechnet wird
  163.       if (len2 >= 2*len1)
  164.         { SAVE_NUM_STACK
  165.           # Teilprodukte von jeweils len1 mal len1 Digits bilden:
  166.           var reg8 uintC k_lo = floor(len1,2); # LΣnge der Low-Teile: floor(len1/2) >0
  167.           var reg8 uintC k_hi = len1 - k_lo; # LΣnge der High-Teile: ceiling(len1/2) >0
  168.           # Es gilt k_lo <= k_hi <= len1, k_lo + k_hi = len1.
  169.           # Summe x1+x0 berechnen:
  170.           var reg6 uintD* sum1_MSDptr;
  171.           var reg6 uintC sum1_len = k_hi; # = max(k_lo,k_hi)
  172.           var reg6 uintD* sum1_LSDptr;
  173.           num_stack_need(sum1_len+1,sum1_MSDptr=,sum1_LSDptr=);
  174.           sum1_MSDptr++; # 1 Digit vorne als Reserve
  175.           {var reg1 uintD carry = # Hauptteile von x1 und x0 addieren:
  176.              add_loop_down(&sourceptr1[-(uintP)k_lo],sourceptr1,sum1_LSDptr,k_lo);
  177.            if (!(k_lo==k_hi))
  178.              # noch k_hi-k_lo = 1 Digits abzulegen
  179.              { sum1_MSDptr[0] = sourceptr1[-(uintP)len1]; # = sourceptr1[-2*(uintP)k_lo-1]
  180.                if (!(carry==0)) { if (++(sum1_MSDptr[0]) == 0) carry=1; else carry=0; }
  181.              }
  182.            if (carry) { *--sum1_MSDptr = 1; sum1_len++; }
  183.           }
  184.          {  # Platz fⁿr Summe y1+y0 belegen:
  185.             var reg6 uintC sum2_maxlen = k_hi+1;
  186.             var reg6 uintD* sum2_LSDptr;
  187.             num_stack_need(sum2_maxlen,,sum2_LSDptr=);
  188.             # Platz fⁿr Produkte x0*y0, x1*y1 belegen:
  189.           { var reg6 uintD* prod_MSDptr;
  190.             var reg6 uintD* prod_LSDptr;
  191.             var reg6 uintD* prodhi_LSDptr;
  192.             num_stack_need(2*(uintL)len1,prod_MSDptr=,prod_LSDptr=);
  193.             prodhi_LSDptr = &prod_LSDptr[-2*(uintP)k_lo];
  194.             # prod_MSDptr/2*len1/prod_LSDptr wird zuerst die beiden
  195.             # Produkte x1*y1 in prod_MSDptr/2*k_hi/prodhi_LSDptr
  196.             #      und x0*y0 in prodhi_LSDptr/2*k_lo/prod_LSDptr,
  197.             # dann das Produkt (b^k*x1+x0)*(b^k*y1+y0) enthalten.
  198.             # Platz fⁿrs Produkt (x1+x0)*(y1+y0) belegen:
  199.            {var reg6 uintD* prodmid_MSDptr;
  200.             var reg6 uintD* prodmid_LSDptr;
  201.             num_stack_need(sum1_len+sum2_maxlen,prodmid_MSDptr=,prodmid_LSDptr=);
  202.             # Schleife ⁿber die hinteren Einzelteile:
  203.             do { # Produkt x0*y0 berechnen:
  204.                  mulu_2loop_down(sourceptr1,k_lo,sourceptr2,k_lo,prod_LSDptr);
  205.                  # Produkt x1*y1 berechnen:
  206.                  mulu_2loop_down(&sourceptr1[-(uintP)k_lo],k_hi,&sourceptr2[-(uintP)k_lo],k_hi,prodhi_LSDptr);
  207.                  # Summe y1+y0 berechnen:
  208.                 {var reg6 uintC sum2_len = k_hi; # = max(k_lo,k_hi)
  209.                  var reg6 uintD* sum2_MSDptr = &sum2_LSDptr[-(uintP)sum2_len];
  210.                  {var reg1 uintD carry = # Hauptteile von y1 und y0 addieren:
  211.                     add_loop_down(&sourceptr2[-(uintP)k_lo],sourceptr2,sum2_LSDptr,k_lo);
  212.                   if (!(k_lo==k_hi))
  213.                     # noch k_hi-k_lo = 1 Digits abzulegen
  214.                     { sum2_MSDptr[0] = sourceptr2[-(uintP)len1]; # = sourceptr2[-2*(uintP)k_lo-1]
  215.                       if (!(carry==0)) { if (++(sum2_MSDptr[0]) == 0) carry=1; else carry=0; }
  216.                     }
  217.                   if (carry) { *--sum2_MSDptr = 1; sum2_len++; }
  218.                  }
  219.                  # Produkt (x1+x0)*(y1+y0) berechnen:
  220.                  mulu_2loop_down(sum1_LSDptr,sum1_len,sum2_LSDptr,sum2_len,prodmid_LSDptr);
  221.                  # Das Produkt beansprucht  2*k_hi + (0 oder 1) <= sum1_len + sum2_len  Digits.
  222.                  {var reg4 uintC prodmid_len = sum1_len+sum2_len;
  223.                   # Davon x1*y1 abziehen:
  224.                   {var reg1 uintD carry =
  225.                      subfrom_loop_down(prodhi_LSDptr,prodmid_LSDptr,2*k_hi);
  226.                    # Falls Carry: Produkt beansprucht 2*k_hi+1 Digits.
  227.                    # Carry um maximal 1 Digit weitertragen:
  228.                    if (!(carry==0)) { prodmid_LSDptr[-2*(uintP)k_hi-1] -= 1; }
  229.                   }
  230.                   # Und x0*y0 abziehen:
  231.                   {var reg1 uintD carry =
  232.                      subfrom_loop_down(prod_LSDptr,prodmid_LSDptr,2*k_lo);
  233.                    # Carry um maximal prodmid_len-2*k_lo Digits weitertragen:
  234.                    if (!(carry==0))
  235.                      { dec_loop_down(&prodmid_LSDptr[-2*(uintP)k_lo],prodmid_len-2*k_lo); }
  236.                   }
  237.                   # prodmid_LSDptr[-prodmid_len..-1] enthΣlt nun x0*y1+x1*y0.
  238.                   # Dies wird zu prod = x1*y1*b^(2*k) + x0*y0 addiert:
  239.                   {var reg1 uintD carry =
  240.                      addto_loop_down(prodmid_LSDptr,&prod_LSDptr[-(uintP)k_lo],prodmid_len);
  241.                      # (Benutze dabei k_lo+prodmid_len <= k_lo+2*(k_hi+1) = 2*len1-k_lo+2 <= 2*len1 .)
  242.                    if (!(carry==0))
  243.                      { inc_loop_down(&prod_LSDptr[-(uintP)(k_lo+prodmid_len)],2*len1-(k_lo+prodmid_len)); }
  244.                 }}}
  245.                  # Das Teilprodukt zum Gesamtprodukt addieren:
  246.                  if (first_part)
  247.                    { copy_loop_down(prod_LSDptr,destptr,2*len1);
  248.                      destptr -= len1;
  249.                      first_part = FALSE;
  250.                    }
  251.                    else
  252.                    { var reg1 uintD carry =
  253.                        addto_loop_down(prod_LSDptr,destptr,len1);
  254.                      destptr -= len1;
  255.                      copy_loop_down(&prod_LSDptr[-(uintP)len1],destptr,len1);
  256.                      if (!(carry==0)) { inc_loop_down(destptr,len1); }
  257.                    }
  258.                  sourceptr2 -= len1; len2 -= len1;
  259.                }
  260.                while (len2 >= 2*len1);
  261.          }}}
  262.          RESTORE_NUM_STACK
  263.         }
  264.       # Nun ist len1 <= len2 < 2*len1.
  265.       # letztes Teilprodukt von len1 mal len2 Digits bilden:
  266.      {SAVE_NUM_STACK
  267.       var reg6 uintD* prod_MSDptr;
  268.       var reg6 uintC prod_len = len1+len2;
  269.       var reg6 uintD* prod_LSDptr;
  270.       num_stack_need((uintL)prod_len,prod_MSDptr=,prod_LSDptr=);
  271.       { var reg8 uintC k_hi = floor(len2,2); # LΣnge der High-Teile: floor(len2/2) >0
  272.         var reg8 uintC k_lo = len2 - k_hi; # LΣnge der Low-Teile: ceiling(len2/2) >0
  273.         # Es gilt k_hi <= k_lo <= len1 <= len2, k_lo + k_hi = len2.
  274.         var reg1 uintC x1_len = len1-k_lo; # <= len2-k_lo = k_hi <= k_lo
  275.         # Summe x1+x0 berechnen:
  276.         var reg6 uintD* sum1_MSDptr;
  277.         var reg6 uintC sum1_len = k_lo; # = max(k_lo,k_hi)
  278.         var reg6 uintD* sum1_LSDptr;
  279.         num_stack_need(sum1_len+1,sum1_MSDptr=,sum1_LSDptr=);
  280.         sum1_MSDptr++; # 1 Digit vorne als Reserve
  281.         {var reg1 uintD carry = # x1 und unteren Teil von x0 addieren:
  282.            add_loop_down(&sourceptr1[-(uintP)k_lo],sourceptr1,sum1_LSDptr,x1_len);
  283.          # und den oberen Teil von x0 dazu:
  284.          copy_loop_down(&sourceptr1[-(uintP)x1_len],&sum1_LSDptr[-(uintP)x1_len],k_lo-x1_len);
  285.          if (!(carry==0))
  286.            { carry = inc_loop_down(&sum1_LSDptr[-(uintP)x1_len],k_lo-x1_len);
  287.              if (carry) { *--sum1_MSDptr = 1; sum1_len++; }
  288.            }
  289.         }
  290.        {# Summe y1+y0 berechnen:
  291.         var reg6 uintD* sum2_MSDptr;
  292.         var reg6 uintC sum2_len = k_lo; # = max(k_lo,k_hi)
  293.         var reg6 uintD* sum2_LSDptr;
  294.         num_stack_need(sum2_len+1,sum2_MSDptr=,sum2_LSDptr=);
  295.         sum2_MSDptr++; # 1 Digit vorne als Reserve
  296.         {var reg1 uintD carry = # Hauptteile von y1 und y0 addieren:
  297.            add_loop_down(&sourceptr2[-(uintP)k_lo],sourceptr2,sum2_LSDptr,k_hi);
  298.          if (!(k_lo==k_hi))
  299.            # noch k_lo-k_hi = 1 Digits abzulegen
  300.            { sum2_MSDptr[0] = sourceptr2[-(uintP)k_lo]; # = sourceptr2[-(uintP)k_hi-1]
  301.              if (!(carry==0)) { if (++(sum2_MSDptr[0]) == 0) carry=1; else carry=0; }
  302.            }
  303.          if (carry) { *--sum2_MSDptr = 1; sum2_len++; }
  304.         }
  305.         # Platz fⁿr Produkte x0*y0, x1*y1:
  306.         { var reg6 uintC prodhi_len = x1_len+k_hi;
  307.           var reg6 uintD* prodhi_LSDptr = &prod_LSDptr[-2*(uintP)k_lo];
  308.           # prod_MSDptr/len1+len2/prod_LSDptr wird zuerst die beiden
  309.           # Produkte x1*y1 in prod_MSDptr/x1_len+k_hi/prodhi_LSDptr
  310.           #      und x0*y0 in prodhi_LSDptr/2*k_lo/prod_LSDptr,
  311.           # dann das Produkt (b^k*x1+x0)*(b^k*y1+y0) enthalten.
  312.           # Platz fⁿrs Produkt (x1+x0)*(y1+y0) belegen:
  313.          {var reg6 uintD* prodmid_MSDptr;
  314.           var reg6 uintC prodmid_len = sum1_len+sum2_len;
  315.           var reg6 uintD* prodmid_LSDptr;
  316.           num_stack_need(prodmid_len,prodmid_MSDptr=,prodmid_LSDptr=);
  317.           # Produkt (x1+x0)*(y1+y0) berechnen:
  318.           mulu_2loop_down(sum1_LSDptr,sum1_len,sum2_LSDptr,sum2_len,prodmid_LSDptr);
  319.           # Das Produkt beansprucht  2*k_lo + (0 oder 1) <= sum1_len + sum2_len = prodmid_len  Digits.
  320.           # Produkt x0*y0 berechnen:
  321.           mulu_2loop_down(sourceptr1,k_lo,sourceptr2,k_lo,prod_LSDptr);
  322.           # Produkt x1*y1 berechnen:
  323.           if (!(x1_len==0))
  324.             { mulu_2loop_down(&sourceptr1[-(uintP)k_lo],x1_len,&sourceptr2[-(uintP)k_lo],k_hi,prodhi_LSDptr);
  325.              # Und x1*y1 abziehen:
  326.              {var reg1 uintD carry =
  327.                 subfrom_loop_down(prodhi_LSDptr,prodmid_LSDptr,prodhi_len);
  328.               # Carry um maximal prodmid_len-prodhi_len Digits weitertragen:
  329.               if (!(carry==0))
  330.                 { dec_loop_down(&prodmid_LSDptr[-(uintP)prodhi_len],prodmid_len-prodhi_len); }
  331.             }}
  332.             else
  333.             # Produkt x1*y1=0, nichts abzuziehen
  334.             { clear_loop_down(prodhi_LSDptr,prodhi_len); }
  335.           # Und x0*y0 abziehen:
  336.           {var reg1 uintD carry =
  337.              subfrom_loop_down(prod_LSDptr,prodmid_LSDptr,2*k_lo);
  338.            # Falls Carry: Produkt beansprucht 2*k_lo+1 Digits.
  339.            # Carry um maximal 1 Digit weitertragen:
  340.            if (!(carry==0)) { prodmid_LSDptr[-2*(uintP)k_lo-1] -= 1; }
  341.           }
  342.           # prodmid_LSDptr[-prodmid_len..-1] enthΣlt nun x0*y1+x1*y0.
  343.           # Dies ist < b^k_lo * b^k_hi + b^x1_len * b^k_lo
  344.           #          = b^len2 + b^len1 <= 2 * b^len2,
  345.           # pa▀t also in len2+1 Digits.
  346.           # Im Fall x1_len=0 ist es sogar < b^k_lo * b^k_hi = b^len2,
  347.           # es pa▀t also in len2 Digits.
  348.           # prodmid_len, wenn m÷glich, um maximal 2 verkleinern:
  349.           # (benutzt prodmid_len >= 2*k_lo >= len2 >= 2)
  350.           if (prodmid_MSDptr[0]==0)
  351.             { prodmid_len--;
  352.               if (prodmid_MSDptr[1]==0) { prodmid_len--; }
  353.             }
  354.           # Nun ist k_lo+prodmid_len <= len1+len2 .
  355.           # (Denn es war prodmid_len = sum1_len+sum2_len <= 2*(k_lo+1)
  356.           #  <= len2+3, und nach 2-maliger Verkleinerung jedenfalls
  357.           #  prodmid_len <= len2+1. Im Falle k_lo < len1 also
  358.           #  k_lo + prodmid_len <= (len1-1)+(len2+1) = len1+len2.
  359.           #  Im Falle k_lo = len1 aber ist x1_len=0, sum1_len = k_lo, also
  360.           #  war prodmid_len = sum1_len+sum2_len <= 2*k_lo+1 <= len2+2,
  361.           #  nach 2-maliger Verkleinerung jedenfalls prodmid_len <= len2.)
  362.           # prodmid*b^k = (x0*y1+x1*y0)*b^k zu prod = x1*y1*b^(2*k) + x0*y0 addieren:
  363.           {var reg1 uintD carry =
  364.              addto_loop_down(prodmid_LSDptr,&prod_LSDptr[-(uintP)k_lo],prodmid_len);
  365.            if (!(carry==0))
  366.              { inc_loop_down(&prod_LSDptr[-(uintP)(k_lo+prodmid_len)],prod_len-(k_lo+prodmid_len)); }
  367.       }}}}}
  368.       # Das Teilprodukt zum Gesamtprodukt addieren:
  369.       if (first_part)
  370.         { copy_loop_down(prod_LSDptr,destptr,prod_len); }
  371.         else
  372.         { var reg1 uintD carry =
  373.             addto_loop_down(prod_LSDptr,destptr,len1);
  374.           destptr -= len1;
  375.           copy_loop_down(&prod_LSDptr[-(uintP)len1],destptr,len2);
  376.           if (!(carry==0)) { inc_loop_down(destptr,len2); }
  377.         }
  378.       RESTORE_NUM_STACK
  379.     }}
  380.  
  381. # Multipliziert zwei Digit-sequences.
  382. # DS_DS_mal_DS(MSDptr1,len1,LSDptr1, MSDptr2,len2,LSDptr2, MSDptr=,len=,LSDptr=);
  383. # multipliziert die DS MSDptr1/len1/LSDptr1 und MSDptr2/len2/LSDptr2.
  384. # Dabei sollte len1>0 und len2>0 sein. Alles sollten Variablen sein!
  385. # Ergebnis ist die DS MSDptr/len/LSDptr, mit len=len1+len2, im Stack.
  386. # Dabei wird num_stack erniedrigt.
  387.   # Methode:
  388.   # Erst unsigned multiplizieren. Dann bis zu zwei Subtraktionen.
  389.   # Sei b=2^intDsize, k=len1, l=len2, n=DS1, m=DS2.
  390.   # Gesucht ist n * m.
  391.   # Wir errechnen erst das unsigned-product p (mod b^(k+l)).
  392.   # n>0, m>0: p = n*m,             n*m = p
  393.   # n<0, m>0: p = (n+b^k)*m,       n*m + b^(k+l) = p - b^k * m (mod b^(k+l)).
  394.   # n>0, m<0: p = n*(m+b^l),       n*m + b^(k+l) = p - b^l * n (mod b^(k+l)).
  395.   # n<0, m<0: p = (n+b^k)*(m+b^l),
  396.   #           n*m = p - b^k * (m+b^l) - b^l * (n+b^k) (mod b^(k+l)).
  397.   #define DS_DS_mal_DS(MSDptr1,len1,LSDptr1,MSDptr2,len2,LSDptr2, MSDptr_zuweisung,len_zuweisung,LSDptr_zuweisung)  \
  398.     { var reg3 uintD* LSDptr0;                                      \
  399.       UDS_UDS_mal_UDS(len1,LSDptr1,len2,LSDptr2, MSDptr_zuweisung,len_zuweisung,LSDptr_zuweisung LSDptr0 = ); \
  400.       if ((sintD)(MSDptr1[0]) < 0) # n<0 ?                          \
  401.         # mu▀ m bzw. m+b^l subtrahieren, um k Digits verschoben:    \
  402.         { subfrom_loop_down(LSDptr2,&LSDptr0[-(uintP)len1],len2); } \
  403.       if ((sintD)(MSDptr2[0]) < 0) # m<0 ?                          \
  404.         # mu▀ n bzw. n+b^k subtrahieren, um l Digits verschoben:    \
  405.         { subfrom_loop_down(LSDptr1,&LSDptr0[-(uintP)len2],len1); } \
  406.     }
  407.  
  408. # (* x y), wo x und y Integers sind. Ergebnis Integer.
  409. # kann GC ausl÷sen
  410.   # Methode:
  411.   # x=0 oder y=0 -> Ergebnis 0
  412.   # x und y beide Fixnums -> direkt multiplizieren
  413.   # sonst: zu WS machen, multiplizieren.
  414.   local object I_I_mal_I (object x, object y);
  415.   local object I_I_mal_I(x,y)
  416.     var reg1 object x;
  417.     var reg2 object y;
  418.     { if (eq(x,Fixnum_0) || eq(y,Fixnum_0))
  419.         { return Fixnum_0; }
  420.       if (I_fixnump(x) && I_fixnump(y))
  421.         { var reg4 sint32 x_ = FN_to_L(x);
  422.           var reg5 sint32 y_ = FN_to_L(y);
  423.          #if (oint_data_len+1 > intLsize)
  424.           # nur falls x und y Integers mit h÷chstens 32 Bit sind:
  425.           if ((((sint32)R_sign(x) ^ x_) >= 0) && (((sint32)R_sign(y) ^ y_) >= 0))
  426.          #endif
  427.          {# Werte direkt multiplizieren:
  428.           var reg3 uint32 hi;
  429.           var reg6 uint32 lo;
  430.           mulu32((uint32)x_,(uint32)y_,hi=,lo=); # erst unsigned multiplizieren
  431.           if (x_ < 0) { hi -= (uint32)y_; } # dann Korrektur fⁿr Vorzeichen
  432.           if (y_ < 0) { hi -= (uint32)x_; } # (vgl. DS_DS_mal_DS)
  433.           return L2_to_I(hi,lo);
  434.         }}
  435.      {SAVE_NUM_STACK # num_stack retten
  436.       var reg6 uintD* xMSDptr;
  437.       var reg6 uintC xlen;
  438.       var reg6 uintD* xLSDptr;
  439.       var reg6 uintD* yMSDptr;
  440.       var reg6 uintC ylen;
  441.       var reg6 uintD* yLSDptr;
  442.       var reg5 uintD* ergMSDptr;
  443.       var reg5 uintC erglen;
  444.       I_to_NDS_nocopy(x, xMSDptr = , xlen = , xLSDptr = );
  445.       I_to_NDS_nocopy(y, yMSDptr = , ylen = , yLSDptr = );
  446.       DS_DS_mal_DS(xMSDptr,xlen,xLSDptr,yMSDptr,ylen,yLSDptr, ergMSDptr=,erglen=,);
  447.       RESTORE_NUM_STACK # num_stack (vorzeitig) zurⁿck
  448.       return DS_to_I(ergMSDptr,erglen);
  449.     }}
  450.  
  451. # (EXPT x y), wo x Integer, y Integer >0 ist.
  452. # kann GC ausl÷sen
  453.   # Methode:
  454.   #   a:=x, b:=y, c:=1. [a^b*c bleibt invariant, = x^y.]
  455.   #   Solange b>1,
  456.   #     falls b ungerade, setze c:=a*c,
  457.   #     setze b:=floor(b/2),
  458.   #     setze a:=a*a.
  459.   #   Wenn b=1, setze c:=a*c.
  460.   #   Liefere c.
  461.   # Oder optimiert:
  462.   #   a:=x, b:=y.
  463.   #   Solange b gerade, setze a:=a*a, b:=b/2. [a^b bleibt invariant, = x^y.]
  464.   #   c:=a.
  465.   #   Solange b:=floor(b/2) >0 ist,
  466.   #     setze a:=a*a, und falls b ungerade, setze c:=a*c.
  467.   #   Liefere c.
  468.   local object I_I_expt_I (object x, object y);
  469.   #if 0 # unoptimiert
  470.   local object I_I_expt_I(x,y)
  471.     var reg1 object x;
  472.     var reg2 object y;
  473.     { pushSTACK(x); pushSTACK(Fixnum_1); pushSTACK(y);
  474.       # Stackaufbau: a, c, b.
  475.       until (eq(STACK_0,Fixnum_1)) # solange bis b=1
  476.         { if (I_oddp(STACK_0)) # b ungerade?
  477.             { STACK_1 = I_I_mal_I(STACK_2,STACK_1); } # c:=a*c
  478.           STACK_0 = I_I_ash_I(STACK_0,Fixnum_minus1); # b := (ash b -1) = (floor b 2)
  479.           STACK_2 = I_I_mal_I(STACK_2,STACK_2); # a:=a*a
  480.         }
  481.       skipSTACK(1);
  482.      {var reg1 object c = popSTACK();
  483.       var reg2 object a = popSTACK();
  484.       return I_I_mal_I(a,c); # a*c als Ergebnis
  485.     }}
  486.   #else # optimiert
  487.   local object I_I_expt_I(x,y)
  488.     var reg3 object x;
  489.     var reg1 object y;
  490.     { pushSTACK(x); pushSTACK(y);
  491.       # Stackaufbau: a, b.
  492.       while (!I_oddp(y))
  493.         { var reg2 object a = STACK_1; STACK_1 = I_I_mal_I(a,a); # a:=a*a
  494.           STACK_0 = y = I_I_ash_I(STACK_0,Fixnum_minus1); # b := (ash b -1)
  495.         }
  496.       pushSTACK(STACK_1); # c:=a
  497.       # Stackaufbau: a, b, c.
  498.       until (eq(y=STACK_1,Fixnum_1)) # Solange b/=1
  499.         { STACK_1 = I_I_ash_I(y,Fixnum_minus1); # b := (ash b -1)
  500.          {var reg2 object a = STACK_2; STACK_2 = a = I_I_mal_I(a,a); # a:=a*a
  501.           if (I_oddp(STACK_1)) { STACK_0 = I_I_mal_I(a,STACK_0); } # evtl. c:=a*c
  502.         }}
  503.       {var reg2 object erg = STACK_0; skipSTACK(3); return erg; }
  504.     }
  505.   #endif
  506.  
  507. # FakultΣt (! n), wo n Fixnum >=0 ist. Ergebnis Integer.
  508. # kann GC ausl÷sen
  509.   # Methode:
  510.   # n <= 10 -> Ergebnis (Fixnum) aus Tabelle
  511.   # Sonst:
  512.   #   Zweierpotenzen extra am Schlu▀ durch einen Shift um
  513.   #   ord2(n!) = sum(k>=1, floor(n/2^k) ) = n - logcount(n)  Bits.
  514.   #   Fⁿr k>=1 wird jede ungerade Zahl m im Intervall n/2^k < m <= n/2^(k-1)
  515.   #   genau k mal gebraucht (als ungerader Anteil von m*2^0,...,m*2^(k-1) ).
  516.   #   Zur Bestimmung des Produkts aller ungeraden Zahlen in einem Intervall
  517.   #   a < m <= b verwenden wir eine rekursive Funktion, die nach Divide-and-
  518.   #   Conquer das Produkt ⁿber die Intervalle a < m <= c und c < m <= b
  519.   #   (c := floor((a+b)/2)) bestimmt und beide zusammenmultipliziert. Dies
  520.   #   vermeidet, da▀ oft gro▀e Zahlen mit ganz kleinen Zahlen multipliziert
  521.   #   werden.
  522.   local object FN_fak_I (object n);
  523.   # UP fⁿr FakultΣt:
  524.   # Bilde das Produkt prod(a < i <= b, 2*i+1), wobei 0 <= a < b klein.
  525.     local object prod_ungerade (uintL a, uintL b);
  526.     local object prod_ungerade(a,b)
  527.       var reg6 uintL a;
  528.       var reg5 uintL b;
  529.       { var reg4 uintL diff = b-a; # Anzahl der Faktoren
  530.         if (diff <= 4)
  531.           # Produkt iterativ bilden
  532.           { var reg1 object faktor = fixnum(2*b+1); # 2*b+1 als letzter Faktor
  533.             var reg2 object produkt = faktor;
  534.             var reg3 uintC count;
  535.             dotimesC(count,diff-1,
  536.               { faktor = fixnum_inc(faktor,-2); # nΣchster Faktor
  537.                 produkt = I_I_mal_I(faktor,produkt); # mit bisherigem Produkt multiplizieren
  538.               });
  539.             return produkt;
  540.           }
  541.           else
  542.           # Produkt rekursiv bilden
  543.           { var reg2 uintL c = floor(a+b,2); # c:=floor((a+b)/2)
  544.             var reg1 object teil = prod_ungerade(a,c); # erstes Teilprodukt
  545.             pushSTACK(teil);
  546.             teil = prod_ungerade(c,b); # zweites Teilprodukt
  547.             return I_I_mal_I(popSTACK(),teil); # und beide multiplizieren
  548.           }
  549.       }
  550.   local object FN_fak_I(n)
  551.     var reg2 object n;
  552.     { local var object fakul_table [] = {
  553.         Fixnum_1,
  554.         fixnum(1UL),
  555.         fixnum(1UL*2),
  556.         #if (oint_data_len>=3)
  557.         fixnum(1UL*2*3),
  558.         #if (oint_data_len>=5)
  559.         fixnum(1UL*2*3*4),
  560.         #if (oint_data_len>=7)
  561.         fixnum(1UL*2*3*4*5),
  562.         #if (oint_data_len>=10)
  563.         fixnum(1UL*2*3*4*5*6),
  564.         #if (oint_data_len>=13)
  565.         fixnum(1UL*2*3*4*5*6*7),
  566.         #if (oint_data_len>=16)
  567.         fixnum(1UL*2*3*4*5*6*7*8),
  568.         #if (oint_data_len>=19)
  569.         fixnum(1UL*2*3*4*5*6*7*8*9),
  570.         #if (oint_data_len>=22)
  571.         fixnum(1UL*2*3*4*5*6*7*8*9*10),
  572.         #if (oint_data_len>=26)
  573.         fixnum(1UL*2*3*4*5*6*7*8*9*10*11),
  574.         #if (oint_data_len>=29)
  575.         fixnum(1UL*2*3*4*5*6*7*8*9*10*11*12),
  576.         #if (oint_data_len>=33)
  577.         ...
  578.         #endif
  579.         #endif
  580.         #endif
  581.         #endif
  582.         #endif
  583.         #endif
  584.         #endif
  585.         #endif
  586.         #endif
  587.         #endif
  588.         #endif
  589.         };
  590.       var reg1 uintL n_ = posfixnum_to_L(n);
  591.       if (n_ < sizeof(fakul_table)/sizeof(object))
  592.         { return fakul_table[n_]; }
  593.         else
  594.         { pushSTACK(Fixnum_1); # bisheriges Produkt := 1
  595.           pushSTACK(n);        # n
  596.           pushSTACK(Fixnum_1); # k := 1
  597.           pushSTACK(n);        # obere Intervallgrenze floor(n/2^(k-1))
  598.           loop
  599.             { # Stackaufbau: prod, n, k, floor(n/2^(k-1)).
  600.               # 'n' enthΣlt floor(n/2^(k-1)).
  601.               n = I_I_ash_I(n,Fixnum_minus1); # untere Grenze floor(n/2^k)
  602.               # 'n' enthΣlt floor(n/2^k).
  603.               # Bilde Teilprodukt prod(A < i <= B & oddp(i), i)
  604.               #       = prod(floor((A-1)/2) < i <= floor((B-1)/2), 2*i+1)
  605.               # wobei B = floor(n/2^(k-1)), A = floor(n/2^k) = floor(B/2).
  606.               { var reg1 uintL b = floor(posfixnum_to_L(STACK_0)-1,2);
  607.                 if (b==0) break; # B=2 oder B=1 -> Produkt fertig
  608.                {var reg1 uintL a = floor(posfixnum_to_L(n)-1,2);
  609.                 pushSTACK(n);
  610.                 {var reg1 object temp = prod_ungerade(a,b);
  611.                  temp = I_I_expt_I(temp,STACK_2); # hoch k nehmen
  612.                  STACK_4 = I_I_mal_I(temp,STACK_4); # und aufmultiplizieren
  613.               }}}
  614.               STACK_2 = fixnum_inc(STACK_2,1); # k:=k+1
  615.               n = popSTACK(); STACK_0 = n;
  616.             }
  617.           skipSTACK(2);
  618.           # Stackaufbau: prod, n.
  619.          {var reg1 object temp = I_logcount_I(STACK_0); # (logcount n)
  620.           temp = I_I_minus_I(popSTACK(),temp); # (- n (logcount n))
  621.           return I_I_ash_I(popSTACK(),temp); # (ash prod (- n (logcount n)))
  622.         }}
  623.     }
  624.  
  625.