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