home *** CD-ROM | disk | FTP | other *** search
/ Amiga ISO Collection / AmigaUtilCD2.iso / Programming / Misc / CLISP-1.LHA / CLISP960530-sr.lha / src / intgcd.d < prev    next >
Encoding:
Text File  |  1996-04-15  |  50.8 KB  |  1,134 lines
  1. # ggT und kgV von Integers
  2.  
  3. # Liefert den ggT zweier Integers.
  4. # I_I_gcd_I(a,b)
  5. # > a,b: zwei Integers
  6. # < ergebnis: (gcd a b), ein Integer >=0
  7. # kann GC auslösen
  8.   local object I_I_gcd_I (object a, object b);
  9.   #define GCD_ALGO 3  # 1: binär, 2: Schulmethode, 3: Lehmer
  10.  
  11. # Liefert den ggT zweier Integers samt Beifaktoren.
  12. # I_I_xgcd_I_I_I(a,b)
  13. # > a,b: zwei Integers
  14. # < STACK_2=u, STACK_1=v, STACK_0=g: Integers mit u*a+v*b = g >= 0
  15. # erniedrigt STACK um 3
  16. # kann GC auslösen
  17.   local void I_I_xgcd_I_I_I (object a, object b);
  18.   #define XGCD_ALGO 3  # 2: Schulmethode, 3: Lehmer
  19. # Im Fall A/=0, B/=0 genügt das Ergebnis (g,u,v) den Ungleichungen:
  20. #   Falls |A| = |B| : g = |A|, u = (signum A), v = 0.
  21. #   Falls |B| | |A|, |B| < |A| : g = |B|, u = 0, v = (signum B).
  22. #   Falls |A| | |B|, |A| < |B| : g = |A|, u = (signum A), v = 0.
  23. #   Sonst: |u| <= |B| / 2*g, |v| <= |A| / 2*g.
  24. #   In jedem Fall |u| <= |B|/g, |v| < |A|/g.
  25. # (Beweis: Im Prinzip macht man ja mehrere Euklid-Schritte auf einmal. Im
  26. # letzten Fall - oBdA |A| > |B| - braucht man mindestens zwei Euklid-Schritte,
  27. # also gilt im Euklid-Tableau
  28. #                 i         |A|            |B|         Erg.
  29. #                --------------------------------------------
  30. #                 0          1              0          |A|
  31. #                 1          0              1          |B|
  32. #                ...        ...            ...         ...
  33. #                n-1  -(-1)^n*x[n-1]  (-1)^n*y[n-1]   z[n-1]
  34. #                 n    (-1)^n*x[n]    -(-1)^n*y[n]     z[n]
  35. #                n+1  -(-1)^n*x[n+1]  (-1)^n*y[n+1]   z[n+1] = 0
  36. #                --------------------------------------------
  37. #                       g = z[n], |u|=x[n], |v|=y[n]
  38. # n>=2, z[0] > ... > z[n-1] > z[n] = g, g | z[n-1], also z[n-1] >= 2*g.
  39. # Da aber mit  (-1)^i*x[i]*|A| - (-1)^i*y[i]*|B| = z[i]  für i=0..n+1
  40. # und            x[i]*y[i+1] - x[i+1]*y[i] = (-1)^i  für i=0..n,
  41. #                x[i]*z[i+1] - x[i+1]*z[i] = (-1)^i*|B|  für i=0..n,
  42. #                y[i]*z[i+1] - y[i+1]*z[i] = -(-1)^i*|A|  für i=0..n
  43. # auch |A| = y[i+1]*z[i] + y[i]*z[i+1], |B| = x[i+1]*z[i] + x[i]*z[i+1]
  44. # für i=0..n (Cramersche Regel), folgt
  45. # |A| = y[n]*z[n-1] + y[n-1]*z[n] >= y[n]*2*g + 0 = |v|*2*g,
  46. # |B| = x[n]*z[n-1] + x[n-1]*z[n] >= x[n]*2*g + 0 = |u|*2*g.)
  47.  
  48. # Liefert den ggT zweier Langworte.
  49. # UL_UL_gcd_UL(a,b)
  50. # > uintL a,b: zwei Langworte >0
  51. # < ergebnis: (gcd a b), ein Langwort >0
  52. #if GCD_ALGO==2 # nur dann brauchen wir's
  53.   local uintL UL_UL_gcd_UL (uintL a, uintL b);
  54. # binäre Methode:
  55. # (gcd a b) :==
  56. #   (prog ((j 0))
  57. #     1 {a,b >0}
  58. #       (when (oddp a) (if (oddp b) (go 2) (go 4)))
  59. #       (when (oddp b) (go 3))
  60. #       (incf j) (setq a (/ a 2)) (setq b (/ b 2))
  61. #       (go 1)
  62. #     2 {a,b >0, beide ungerade}
  63. #       (cond ((> a b) (setq a (- a b)) (go 3))
  64. #             ((= a b) (go 5))
  65. #             ((< a b) (setq b (- b a)) (go 4))
  66. #       )
  67. #     3 {a,b >0, a gerade, b ungerade}
  68. #       (repeat (setq a (/ a 2)) (until (oddp a)))
  69. #       (go 2)
  70. #     4 {a,b >0, a ungerade, b gerade}
  71. #       (repeat (setq b (/ b 2)) (until (oddp b)))
  72. #       (go 2)
  73. #     5 {a=b>0}
  74. #       (return (ash a j))
  75. #   )
  76. # Statt j zu erhöhen und immer Bit 0 von a und b abfragen,
  77. # fragen wir stattdessen immer Bit j von a und b ab; Bits j-1..0 sind =0.
  78.   local uintL UL_UL_gcd_UL(a,b)
  79.     var reg1 uintL a;
  80.     var reg2 uintL b;
  81.     {
  82.       #ifdef DUMMER_GGT # so macht's ein Mathematiker:
  83.       var reg3 uintL bit_j = bit(0);
  84.       loop
  85.         { # a,b >0
  86.           if (!((a & bit_j) ==0))
  87.             { if (!((b & bit_j) ==0)) goto odd_odd; else goto odd_even; }
  88.           if (!((b & bit_j) ==0)) goto even_odd;
  89.           # a,b >0 gerade
  90.           bit_j = bit_j<<1;
  91.         }
  92.       #else # Trick von B. Degel:
  93.       var reg3 uintL bit_j = (a | b); # endet mit einer 1 und j Nullen
  94.       bit_j = bit_j ^ (bit_j - 1); # Maske = bit(j) | bit(j-1) | ... | bit(0)
  95.       if (!((a & bit_j) ==0))
  96.         { if (!((b & bit_j) ==0)) goto odd_odd; else goto odd_even; }
  97.       if (!((b & bit_j) ==0)) goto even_odd;
  98.       #endif
  99.       loop
  100.         { odd_odd: # a,b >0, beide ungerade
  101.           # Vergleiche a und b:
  102.           if (a == b) break; # a=b>0 -> fertig
  103.           if (a > b) # a>b ?
  104.             { a = a-b;
  105.               even_odd: # a,b >0, a gerade, b ungerade
  106.               do { a = a>>1; } while ((a & bit_j) ==0);
  107.             }
  108.             else # a<b
  109.             { b = b-a;
  110.               odd_even: # a,b >0, a ungerade, b gerade
  111.               do { b = b>>1; } while ((b & bit_j) ==0);
  112.             }
  113.         }
  114.       # a=b>0
  115.       return a;
  116.     }
  117. #endif
  118.  
  119. # binäre Methode:
  120. # (gcd a b) :==
  121. # b=0 --> (abs a)
  122. # a=0 --> (abs b)
  123. # sonst:
  124. #   (abs a) und (abs b) in zwei Buffer packen, als Unsigned Digit Sequences.
  125. #   [Schreibe oBdA wieder a,b]
  126. #   (prog ((j 0))
  127. #     1 {a,b >0}
  128. #       (if (evenp a)
  129. #         (if (evenp b)
  130. #           (progn (incf j) (setq a (/ a 2)) (setq b (/ b 2)) (go 1))
  131. #           (go 4)
  132. #         )
  133. #         (while (evenp b) (setq b (/ b 2)))
  134. #       )
  135. #     2 {a,b >0, beide ungerade}
  136. #       (cond ((> a b))
  137. #             ((= a b) (go 5))
  138. #             ((< a b) (rotatef a b))
  139. #       )
  140. #     3 {a,b >0, beide ungerade, a>b}
  141. #       (setq a (- a b))
  142. #     4 {a,b >0, a gerade, b ungerade}
  143. #       (repeat (setq a (/ a 2)) (until (oddp a)))
  144. #       (go 2)
  145. #     5 {a=b>0}
  146. #       (return (ash a j))
  147. #   )
  148. # weil es oft auftritt (insbesondere bei GCD's mehrerer Zahlen):
  149. # a=1 oder b=1 --> 1
  150. #if GCD_ALGO==1
  151.   local object I_I_gcd_I(a,b)
  152.     var reg8 object a;
  153.     var reg9 object b;
  154.     { if (eq(a,Fixnum_1)) { return a; } # a=1 -> 1
  155.       if (eq(b,Fixnum_1)) { return b; } # b=1 -> 1
  156.       if (eq(b,Fixnum_0)) { return I_abs_I(a); } # b=0 -> (abs a)
  157.       if (eq(a,Fixnum_0)) { return I_abs_I(b); } # a=0 -> (abs b)
  158.      {SAVE_NUM_STACK # num_stack retten
  159.       var reg1 uintD* a_MSDptr;
  160.       var reg3 uintC a_len;
  161.       var reg5 uintD* a_LSDptr;
  162.       var reg2 uintD* b_MSDptr;
  163.       var reg4 uintC b_len;
  164.       var reg6 uintD* b_LSDptr;
  165.       # Macro: erzeugt die NUDS zu (abs x), erniedrigt num_stack
  166.       #define I_abs_to_NUDS(x)  \
  167.         { I_to_NDS_1(x, x##_MSDptr = , x##_len = , x##_LSDptr = ); # (nichtleere) NDS holen \
  168.           if ((sintD)x##_MSDptr[0] < 0) # falls <0, negieren:                               \
  169.             { neg_loop_down(x##_LSDptr,x##_len); }                                          \
  170.           if (x##_MSDptr[0] == 0) # normalisieren (max. 1 Nulldigit entfernen)              \
  171.             { x##_MSDptr++; x##_len--; }                                                    \
  172.         }
  173.       I_abs_to_NUDS(a); # (abs a) als NUDS erzeugen
  174.       I_abs_to_NUDS(b); # (abs b) als NUDS erzeugen
  175.       # Jetzt ist a = a_MSDptr/a_len/a_LSDptr, b = b_MSDptr/b_len/b_LSDptr,
  176.       # beides NUDS, und a_len>0, b_len>0.
  177.       # Macro: Halbiere x.
  178.       #define halb(x)  \
  179.         { shift1right_loop_up(x##_MSDptr,x##_len,0); # um 1 Bit rechts schieben \
  180.           if (x##_MSDptr[0] == 0) { x##_MSDptr++; x##_len--; } # normalisieren  \
  181.         }
  182.       # Macro: Ob x gerade ist.
  183.       #define evenp(x)  \
  184.         ((x##_LSDptr[-1] & bit(0)) ==0)
  185.       { var reg7 uintL j = 0;
  186.         label_1: # a,b >0
  187.           if (evenp(a))
  188.             { if (evenp(b))
  189.                 { j++; halb(a); halb(b); goto label_1; }
  190.                 else
  191.                 goto label_4;
  192.             }
  193.           while (evenp(b)) { halb(b); }
  194.         label_2: # a,b >0, beide ungerade
  195.           # Vergleiche a und b:
  196.           if (a_len > b_len) goto label_3; # a>b ?
  197.           if (a_len == b_len)
  198.             { var reg1 signean vergleich = compare_loop_up(a_MSDptr,b_MSDptr,a_len);
  199.               if (vergleich > 0) goto label_3; # a>b ?
  200.               if (vergleich == 0) goto label_5; # a=b ?
  201.             }
  202.           # a<b -> a,b vertauschen:
  203.           swap(reg1 uintD*, a_MSDptr,b_MSDptr);
  204.           swap(reg1 uintC, a_len,b_len);
  205.           swap(reg1 uintD*, a_LSDptr,b_LSDptr);
  206.         label_3: # a,b >0, beide ungerade, a>b
  207.           # subtrahiere a := a - b
  208.           if (!( subfrom_loop_down(b_LSDptr,a_LSDptr,b_len) ==0))
  209.             { dec_loop_down(&a_LSDptr[-(uintP)b_len],a_len-b_len); }
  210.           while (a_MSDptr[0] == 0) { a_MSDptr++; a_len--; } # normalisieren
  211.         label_4: # a,b >0, a gerade, b ungerade
  212.           do { halb(a); } while (evenp(a));
  213.           goto label_2;
  214.         label_5: # a=b>0
  215.           # a zu einer NDS machen:
  216.           RESTORE_NUM_STACK # num_stack (vorzeitig) zurück
  217.           a = NUDS_to_I(a_MSDptr,a_len); # ggT der ungeraden Anteile als Integer
  218.           return I_I_ash_I(a,fixnum(j)); # (ash a j) als Ergebnis
  219.       }
  220.       #undef evenp
  221.       #undef halb
  222.       #undef I_abs_to_NUDS
  223.     }}
  224. #endif
  225.  
  226. # binäre Methode:
  227. # (xgcd A B) :==
  228. # B=0 --> g = (abs A), u = (signum A), v = 0
  229. # A=0 --> g = (abs B), u = 0, v = (signum B)
  230. # sonst:
  231. # a := (abs A) und b := (abs B) in zwei Buffer packen,
  232. # als Unsigned Digit Sequences.
  233. # (xgcd a b) :==
  234. #   (prog ((j 0) (ua 1) (va 0) (ub 0) (vb 1))
  235. #     {Stets |A|*ua-|B|*va=a*2^j, -|A|*ub+|B|*vb=b*2^j,
  236. #            ua>0, va>=0, ub>=0, vb>0.}
  237. #     1 {a,b >0}
  238. #       (when (oddp a) (setq Aj a Bj b) (if (oddp b) (go 2) (go 4)))
  239. #       (when (oddp b) (setq Aj a Bj b) (go 3))
  240. #       (incf j) (setq a (/ a 2)) (setq b (/ b 2))
  241. #       (go 1)
  242. #     {Ab hier |A| = Aj*2^j, |B| = Bj*2^j, Aj oder Bj ungerade,
  243. #              Aj*ua-Bj*va=a, -Aj*ub+Bj*vb=b, a oder b ungerade,
  244. #              0<ua<=Bj, 0<=va<Aj, 0<=ub<Bj, 0<vb<=Aj. (wieso??)}
  245. #     2 {a,b >0, beide ungerade}
  246. #       (cond ((> a b) (setq a (- a b) ua (+ ua ub) va (+ va vb)) (go 3))
  247. #             ((= a b) (go 5))
  248. #             ((< a b) (setq b (- b a) ub (+ ub ua) vb (+ vb va)) (go 4))
  249. #       )
  250. #     3 {a,b >0, a gerade, b ungerade}
  251. #       (repeat (when (or (oddp ua) (oddp va))
  252. #                 { Falls ua gerade, muß (da Bj*va==0 mod 2) Bj==0, Aj==1 sein.
  253. #                   Falls va gerade, muß (da Aj*ua==0 mod 2) Aj==0, Bj==1 sein.
  254. #                   Falls ua,va beide ungerade, müssen (da Aj*1-Bj*1==0 mod 2)
  255. #                                               Aj und Bj beide ungerade sein.}
  256. #                 (setq ua (+ ua Bj) va (+ va Aj))
  257. #               )
  258. #               {ua,va beide gerade}
  259. #               (setq a (/ a 2) ua (/ ua 2) va (/ va 2))
  260. #               (until (oddp a))
  261. #       )
  262. #       (go 2)
  263. #     4 {a,b >0, a ungerade, b gerade}
  264. #       (repeat (when (or (oddp ub) (oddp vb))
  265. #                 { Falls ub gerade, muß (da Bj*vb==0 mod 2) Bj==0, Aj==1 sein.
  266. #                   Falls vb gerade, muß (da Aj*ub==0 mod 2) Aj==0, Bj==1 sein.
  267. #                   Falls ub,vb beide ungerade, müssen (da -Aj*1+Bj*1==0 mod 2)
  268. #                                               Aj und Bj beide ungerade sein.}
  269. #                 (setq ub (+ ub Bj) vb (+ vb Aj))
  270. #               )
  271. #               {ub,vb beide gerade}
  272. #               (setq b (/ b 2) ub (/ ub 2) vb (/ vb 2))
  273. #               (until (oddp b))
  274. #       )
  275. #       (go 2)
  276. #     5 {a=b>0}
  277. #       (return (values (ash a j) (* (signum A) ua) (- (* (signum B) va))))
  278. #       (return (values (ash a j) (- (* (signum A) ub)) (* (signum B) vb)))
  279. #   )
  280. #if XGCD_ALGO==1
  281. ??
  282. #endif
  283.  
  284. # Schulmethode:
  285. #   (gcd a b) :==
  286. #   [a:=(abs a), b:=(abs b), while b>0 do (a,b) := (b,(mod a b)), -> a]
  287. # verbessert:
  288. # a=0 -> (abs b)
  289. # b=0 -> (abs a)
  290. # a=1 -> 1
  291. # b=1 -> 1
  292. # a:=(abs a), b:=(abs b)
  293. # Falls a=b: return a; falls a<b: vertausche a und b.
  294. # (*) {Hier a>b>0}
  295. # Falls b=1, return 1. {spart eine Division durch 1}
  296. # Sonst dividieren (divide a b), a:=b, b:=Rest.
  297. #       Falls b=0, return a, sonst goto (*).
  298. #if GCD_ALGO==2
  299.   local object I_I_gcd_I(a,b)
  300.     var reg1 object a;
  301.     var reg2 object b;
  302.     { if (eq(a,Fixnum_1)) { return a; } # a=1 -> 1
  303.       if (eq(b,Fixnum_1)) { return b; } # b=1 -> 1
  304.       if (eq(b,Fixnum_0)) { return I_abs_I(a); } # b=0 -> (abs a)
  305.       if (eq(a,Fixnum_0)) { return I_abs_I(b); } # a=0 -> (abs b)
  306.       # Beträge nehmen:
  307.       pushSTACK(b); pushSTACK(I_abs_I(a)); STACK_1 = I_abs_I(STACK_1);
  308.       # Stackaufbau: (abs b), (abs a).
  309.       a = popSTACK(); b = STACK_0;
  310.       # ab jetzt: b in STACK_0.
  311.       if (I_fixnump(a) && I_fixnump(b)) # ggT zweier Fixnums >0
  312.         { # bleibt Fixnum, da (gcd a b) <= (min a b)
  313.           skipSTACK(1);
  314.           return fixnum(UL_UL_gcd_UL(posfixnum_to_L(a),posfixnum_to_L(b)));
  315.         }
  316.       { var reg3 signean vergleich = I_I_comp(a,b);
  317.         if (vergleich == 0) { skipSTACK(1); return a; } # a=b -> fertig
  318.         if (vergleich < 0) { STACK_0 = a; a = b; } # a<b -> a,b vertauschen
  319.       }
  320.       loop # Hier a > STACK_0 > 0
  321.         { if (eq(STACK_0,Fixnum_1)) { return popSTACK(); } # b=1 -> Ergebnis 1
  322.           I_I_divide_I_I(a,STACK_0); b = STACK_0; skipSTACK(2); # b := Rest der Division a / STACK_0
  323.           a = STACK_0; # neues a = altes b
  324.           if (eq(b,Fixnum_0)) { skipSTACK(1); return a; }
  325.           STACK_0 = b;
  326.         }
  327.     }
  328. #endif
  329.  
  330. # Schulmethode:
  331. #   (gcd A B) :==
  332. #   [a:=(abs A), b:=(abs B), while b>0 do (a,b) := (b,(mod a b)), -> a]
  333. # verbessert:
  334. # A=1 -> return g=1, (u,v)=(1,0)
  335. # B=1 -> return g=1, (u,v)=(0,1)
  336. # a:=(abs A), ua:=(signum A), va:=0
  337. # b:=(abs B), ub:=0, vb:=(signum B)
  338. # A=0 -> return g=b, (u,v) = (ub,vb)
  339. # B=0 -> return g=a, (u,v) = (ua,va)
  340. # {Stets ua*A+va*B=a, ub*A+vb*B=b, ua*vb-ub*va = +/- 1.}
  341. # Falls a=b: return a,ua,va;
  342. # falls a<b: vertausche a und b, ua und ub, va und vb.
  343. # (*) {Hier a>b>0}
  344. # Falls b=1, return 1,ub,vb. {spart eine Division durch 1}
  345. # Sonst dividieren (divide a b) -> q,r.
  346. #       Falls r=0, return b,ub,vb.
  347. #       a:=b, b := Rest r = a-q*b, (ua,va,ub,vb) := (ub,vb,ua-q*ub,va-q*vb).
  348. #       goto (*).
  349. #if XGCD_ALGO==2
  350.   local void I_I_xgcd_I_I_I(a,b)
  351.     var reg3 object a;
  352.     var reg2 object b;
  353.     { if (eq(a,Fixnum_1)) # a=1 -> g=1, (u,v)=(1,0)
  354.         { pushSTACK(Fixnum_1); pushSTACK(Fixnum_0);
  355.           pushSTACK(a); return;
  356.         }
  357.       if (eq(b,Fixnum_1)) # b=1 -> g=1, (u,v)=(0,1)
  358.         { pushSTACK(Fixnum_0); pushSTACK(Fixnum_1);
  359.           pushSTACK(b); return;
  360.         }
  361.       # Vorzeichen nehmen:
  362.       pushSTACK(R_minusp(a) ? Fixnum_minus1 : Fixnum_1); # ua := +/- 1
  363.       pushSTACK(Fixnum_0); # va := 0
  364.       pushSTACK(Fixnum_0); # ub := 0
  365.       pushSTACK(R_minusp(b) ? Fixnum_minus1 : Fixnum_1); # vb := +/- 1
  366.       # Beträge nehmen:
  367.       pushSTACK(b); pushSTACK(I_abs_I(a)); STACK_1 = I_abs_I(STACK_1);
  368.       # Stackaufbau: ua, va, ub, vb, b, a.
  369.       a = STACK_0; b = STACK_1;
  370.       if (eq(b,Fixnum_0)) # b=0 -> g=a, (u,v) = (ua,va)
  371.         { skipSTACK(4); pushSTACK(a); return; }
  372.       if (eq(a,Fixnum_0)) # a=0 -> g=b, (u,v) = (ub,vb)
  373.         { STACK_5 = STACK_3; STACK_4 = STACK_2; STACK_3 = STACK_1;
  374.           skipSTACK(3); return;
  375.         }
  376.       { var reg4 signean vergleich = I_I_comp(a,b);
  377.         if (vergleich == 0) # a=b -> fertig
  378.           { skipSTACK(4); pushSTACK(a); return; }
  379.         if (vergleich < 0) # a<b -> a,b vertauschen
  380.           { STACK_0 = b; b = STACK_1 = a; a = STACK_0;
  381.             swap(STACK_5,STACK_3); swap(STACK_4,STACK_2);
  382.       }   }
  383.       loop # Hier a>b>0
  384.         { if (eq(b,Fixnum_1)) # b=1 -> g=b, (u,v) = (ub,vb)
  385.             { STACK_5 = STACK_3; STACK_4 = STACK_2; STACK_3 = b;
  386.               skipSTACK(3); return;
  387.             }
  388.           I_I_divide_I_I(a,b); # Division a / b
  389.           # Stackaufbau: ..., q, r.
  390.           if (eq(STACK_0,Fixnum_0)) # r=0 -> fertig
  391.             { STACK_5 = STACK_3; STACK_4 = STACK_2; STACK_3 = STACK_1;
  392.               skipSTACK(3); return;
  393.             }
  394.           {var reg1 object x = I_I_mal_I(STACK_1,STACK_(3+2)); # q*ub
  395.            x = I_I_minus_I(STACK_(5+2),x); # ua-q*ub
  396.            STACK_(5+2) = STACK_(3+2); STACK_(3+2) = x; # ua := ub, ub := x
  397.           }
  398.           {var reg1 object x = I_I_mal_I(STACK_1,STACK_(2+2)); # q*vb
  399.            x = I_I_minus_I(STACK_(4+2),x); # va-q*vb
  400.            STACK_(4+2) = STACK_(2+2); STACK_(2+2) = x; # va := vb, vb := x
  401.           }
  402.           STACK_(0+2) = a = STACK_(1+2); # a := altes b
  403.           STACK_(1+2) = b = STACK_0; # b := r
  404.           skipSTACK(2);
  405.         }
  406.     }
  407. #endif
  408.  
  409. # Lehmer-Methode:
  410. # vgl. [ D. E. Knuth: The Art of Computer Programming, Vol. 2: Seminumerical
  411. #        Algorithms, Sect. 4.5.2., Algorithm L ]
  412. # und [ Collins, Loos: SAC-2, Algorithms IGCD, DPCC ].
  413. # (gcd a b) :==
  414. # a=0 -> (abs b)
  415. # b=0 -> (abs a)
  416. # a=1 -> 1
  417. # b=1 -> 1
  418. # a:=(abs a), b:=(abs b)
  419. # (*) {Hier a,b>0}
  420. # Falls a=b: return a; falls a<b: vertausche a und b.
  421. # {Hier a>b>0}
  422. # Falls (- (integer-length a) (integer-length b)) >= intDsize/2,
  423. #   lohnt sich eine Division: (a,b) := (b , a mod b). Falls b=0: return a.
  424. # Falls dagegen 0 <= (- (integer-length a) (integer-length b)) < intDsize/2,
  425. #   seien a' die führenden intDsize Bits von a
  426. #   (2^(intDsize-1) <= a' < 2^intDsize) und b' die entsprechenden Bits von b
  427. #   (2^(intDsize/2) <= b' <= a' < 2^intDsize).
  428. #   Rechne den Euklid-Algorithmus mit Beifaktoren für ALLE Zahlen (a,b) aus,
  429. #   die mit a' bzw. b' anfangen; das liefert x1,y1,x2,y2, so daß
  430. #   ggT(a,b) = ggT(x1*a-y1*b,-x2*a+y2*b) und x1*a-y1*b>=0,-x2*a+y2*b>=0.
  431. #   Genauer: Mit offensichtlicher Skalierung betrachten wir
  432. #            a als beliebiges Element des Intervalls [a',a'+1) und
  433. #            b als beliebiges Element des Intervalls [b',b'+1) und
  434. #            führen den Euklid-Algorithmus schrittweise durch:
  435. #            (x1,y1,z1) := (1,0,a'), (x2,y2,z2) := (0,1,b'),
  436. #            Schleife:
  437. #            {Hier x1*a'-y1*b'=z1, x1*a-y1*b in [z1-y1,z1+x1), z1-y1>=0, z1>0,
  438. #             und -x2*a'+y2*b'=z2, -x2*a+y2*b in [z2-x2,z2+y2), z2-x2>=0, z2>0,
  439. #             x1*y2-x2*y1=1, x1*z2+x2*z1=b', y1*z2+y2*z1=a'.}
  440. #            Falls z1-y1>=z2+y2:
  441. #              (x1,y1,z1) := (x1+x2,y1+y2,z1-z2), goto Schleife.
  442. #            Falls z2-x2>=z1+x1:
  443. #              (x2,y2,z2) := (x2+x1,y2+y1,z2-z1), goto Schleife.
  444. #            Sonst muß man abbrechen.
  445. #            {Zu den Schleifeninvarianten:
  446. #             1. Die Gleichungen x1*a'-y1*b'=z1, -x2*a'+y2*b'=z2,
  447. #                x1*y2-x2*y1=1, x1*z2+x2*z1=b', y1*z2+y2*z1=a' mit Induktion.
  448. #             2. Die Ungleichungen x1>0, y1>=0, x2>=0, y2>0 mit Induktion.
  449. #             3. Die Ungleichungen z1>=0, z2>=0 nach Fallauswahl.
  450. #             4. Die Ungleichung x1+x2>0 aus x1*z2+x2*z1=b'>0,
  451. #                die Ungleichung y1+y2>0 aus y1*z2+y2*z1=a'>0.
  452. #             5. Die Ungleichung z1>0 wegen Fallauswahl und y1+y2>0,
  453. #                Die Ungleichung z2>0 wegen Fallauswahl und x1+x2>0.
  454. #             6. Die Ungleichungen z1-y1>=0, z2-x2>=0 wegen Fallauswahl.
  455. #             7. Die Ungleichung max(z1,z2) <= a' mit Induktion.
  456. #             8. Die Ungleichung x1+x2 <= x1*z2+x2*z1 = b',
  457. #                die Ungleichung y1+y2 <= y1*z2+y2*z1 = a'.
  458. #             Damit bleiben alle Größen im Intervall [0,beta), kein Überlauf.
  459. #             9. Die Ungleichungen z1+x1<=beta, z2+y2<=beta mit Induktion.
  460. #             10. x1*a-y1*b in (z1-y1,z1+x1) (bzw. [z1,z1+x1) bei y1=0),
  461. #                -x2*a+y2*b in (z2-x2,z2+y2) (bzw. [z2,z2+y2) bei x2=0),
  462. #                da a in a'+[0,1) und b in b'+[0,1).
  463. #                Jedenfalls 0 < x1*a-y1*b < z1+x1 <= x2*z1+x1*z2 = b' falls x2>0,
  464. #                und        0 < -x2*a+y2*b < z2+y2 <= y1*z2+y2*z1 = a' falls y1>0.}
  465. #            Man kann natürlich auch mehrere Subtraktionsschritte auf einmal
  466. #            durchführen:
  467. #            Falls q := floor((z1-y1)/(z2+y2)) > 0 :
  468. #              (x1,y1,z1) := (x1+q*x2,y1+q*y2,z1-q*z2), goto Schleife.
  469. #            Falls q := floor((z2-x2)/(z1+x1)) > 0 :
  470. #              (x2,y2,z2) := (x2+q*x1,y2+q*y1,z2-q*z1), goto Schleife.
  471. #            {Am Schluß gilt -(x1+x2) < z1-z2 < y1+y2 und daher
  472. #             z2-x2 <= b'/(x1+x2) < z1+x1, z1-y1 <= a'/(y1+y2) < z2+y2,
  473. #             und - unter Berücksichtigung von x1*y2-x2*y1=1 -
  474. #             z1-y1 <= b'/(x1+x2) < z2+y2, z2-x2 <= a'/(y1+y2) < z1+x1,
  475. #             also  max(z1-y1,z2-x2) <= min(b'/(x1+x2),a'/(y1+y2))
  476. #                          <= max(b'/(x1+x2),a'/(y1+y2)) < min(z1+x1,z2+y2).}
  477. #   Im Fall y1=x2=0 => x1=y2=1 (der nur bei a'=z1=z2=b' eintreten kann)
  478. #     ersetze (a,b) := (a-b,b). {Beide >0, da a>b>0 war.}
  479. #   Der Fall y1=0,x2>0 => x1=y2=1 => a' = z1 < z2+x2*z1 = b'
  480. #     kann nicht eintreten.
  481. #   Im Fall x2=0,y1>0 => x1=y2=1 ersetze (a,b) := (a-y1*b,b).
  482. #     {Das ist OK, da 0 <= z1-y1 = a'-y1*(b'+1) < a-y1*b < a.}
  483. #   Sonst (y1>0,x2>0) ersetze (a,b) := (x1*a-y1*b,-x2*a+y2*b).
  484. #     {Das ist OK, da 0 <= z1-y1 = x1*a'-y1*(b'+1) < x1*a-y1*b
  485. #                  und 0 <= z2-x2 = -x2*(a'+1)+y2*b' < -x2*a+y2*b
  486. #      und x1*a-y1*b < x1*(a'+1)-y1*b' = z1+x1 <= x2*z1+x1*z2 = b' <= b
  487. #      und -x2*a+y2*b < -x2*a'+y2*(b'+1) = z2+y2 <= y1*z2+y2*z1 = a' <= a.}
  488. # goto (*).
  489. #if (GCD_ALGO==3) || (XGCD_ALGO==3)
  490.   # Teilfunktion für die Durchführung des Euklid-Algorithmus auf
  491.   # den führenden Ziffern a' und b':
  492.   # partial_gcd(a',b',&erg); mit a'>b'
  493.   # liefert in erg: x1,y1,x2,y2 mit den oben angegebenen Invarianten.
  494.   typedef struct { uintD x1,y1,x2,y2; } partial_gcd_result;
  495.   local void partial_gcd (uintD z1, uintD z2, partial_gcd_result* erg);
  496.   local void partial_gcd(z1,z2,erg) # z1:=a', z2:=b'
  497.     var reg1 uintD z1;
  498.     var reg2 uintD z2;
  499.     var reg7 partial_gcd_result* erg;
  500.     { var reg5 uintD x1 = 1;
  501.       var reg3 uintD y1 = 0;
  502.       var reg6 uintD x2 = 0;
  503.       var reg4 uintD y2 = 1;
  504.       subtract_from_1: # Hier ist z1-y1>=z2+y2.
  505.         # Bestimme q := floor((z1-y1)/(z2+y2)) >= 1 :
  506.         { var reg2 uintD zaehler = z1-y1;
  507.           var reg1 uintD nenner = z2+y2; # z2+y2 <= z1-y1 < beta !
  508.           if (floor(zaehler,8) >= nenner) # zaehler >= 8*nenner ?
  509.             # ja -> Dividieren lohnt sich wohl
  510.             { var reg3 uintD q = floorD(zaehler,nenner);
  511.               x1 += muluD_unchecked(q,x2); # x1 := x1+q*x2
  512.               y1 += muluD_unchecked(q,y2); # y1 := y1+q*y2
  513.               z1 -= muluD_unchecked(q,z2); # z1 := z1-q*z2
  514.             }
  515.             else
  516.             # nein -> ein paarmal subtrahieren ist wohl schneller
  517.             do { x1 += x2; y1 += y2; z1 -= z2; } # (x1,y1,z1) := (x1+x2,y1+y2,z1-z2)
  518.                while (z1-y1 >= nenner);
  519.         }
  520.       if (z2-x2 <= z1+x1-1) goto no_more_subtract;
  521.       subtract_from_2: # Hier ist z2-x2>=z1+x1.
  522.         # Bestimme q := floor((z2-x2)/(z1+x1)) >= 1 :
  523.         { var reg2 uintD zaehler = z2-x2;
  524.           var reg1 uintD nenner = z1+x1; # z1+x1 <= z2-x2 < beta !
  525.           if (floor(zaehler,8) >= nenner) # zaehler >= 8*nenner ?
  526.             # ja -> Dividieren lohnt sich wohl
  527.             { var reg3 uintD q = floorD(zaehler,nenner);
  528.               x2 += muluD_unchecked(q,x1); # x2 := x2+q*x1
  529.               y2 += muluD_unchecked(q,y1); # y2 := y2+q*y1
  530.               z2 -= muluD_unchecked(q,z1); # z2 := z2-q*z1
  531.             }
  532.             else
  533.             # nein -> ein paarmal subtrahieren ist wohl schneller
  534.             do { x2 += x1; y2 += y1; z2 -= z1; } # (x2,y2,z2) := (x2+x1,y2+y1,z2-z1)
  535.                while (z2-x2 >= nenner);
  536.         }
  537.       if (z1-y1 <= z2+y2-1) goto no_more_subtract;
  538.       goto subtract_from_1;
  539.       no_more_subtract: # Keine Subtraktion mehr möglich.
  540.       erg->x1 = x1; erg->y1 = y1; erg->x2 = x2; erg->y2 = y2; # Ergebnis
  541.     }
  542. #endif
  543. #if GCD_ALGO==3
  544.   # Los geht's:
  545.   local object I_I_gcd_I(a,b)
  546.     var reg10 object a;
  547.     var reg10 object b;
  548.     { if (eq(a,Fixnum_1)) { return a; } # a=1 -> 1
  549.       if (eq(b,Fixnum_1)) { return b; } # b=1 -> 1
  550.       if (eq(b,Fixnum_0)) { return I_abs_I(a); } # b=0 -> (abs a)
  551.       if (eq(a,Fixnum_0)) { return I_abs_I(b); } # a=0 -> (abs b)
  552.      {SAVE_NUM_STACK # num_stack retten
  553.       var reg6 uintD* a_MSDptr;
  554.       var reg8 uintC a_len;
  555.       var reg10 uintD* a_LSDptr;
  556.       var reg7 uintD* b_MSDptr;
  557.       var reg9 uintC b_len;
  558.       var reg10 uintD* b_LSDptr;
  559.       # Macro: erzeugt die NUDS zu (abs x), erniedrigt num_stack
  560.       #define I_abs_to_NUDS(x)  \
  561.         { I_to_NDS_1(x, x##_MSDptr = , x##_len = , x##_LSDptr = ); # (nichtleere) NDS holen \
  562.           if ((sintD)x##_MSDptr[0] < 0) # falls <0, negieren:                               \
  563.             { neg_loop_down(x##_LSDptr,x##_len); }                                          \
  564.           if (x##_MSDptr[0] == 0) # normalisieren (max. 1 Nulldigit entfernen)              \
  565.             { x##_MSDptr++; x##_len--; }                                                    \
  566.         }
  567.       I_abs_to_NUDS(a); # (abs a) als NUDS erzeugen
  568.       I_abs_to_NUDS(b); # (abs b) als NUDS erzeugen
  569.       # Jetzt ist a = a_MSDptr/a_len/a_LSDptr, b = b_MSDptr/b_len/b_LSDptr,
  570.       # beides NUDS, und a_len>0, b_len>0.
  571.       # Platz für zwei Rechenregister besorgen, mit je max(a_len,b_len)+1 Digits:
  572.       {var reg10 uintD* divroomptr; # Platz für Divisionsergebnis
  573.        var reg10 uintD* c_LSDptr;
  574.        var reg10 uintD* d_LSDptr;
  575.        {var reg10 uintL c_len = (uintL)(a_len>=b_len ? a_len : b_len) + 1;
  576.         num_stack_need(c_len,divroomptr=,c_LSDptr=);
  577.         num_stack_need(c_len,_EMA_,d_LSDptr=);
  578.         # Jetzt ist ../c_len/c_LSDptr, ../c_len/d_LSDptr frei.
  579.        }
  580.        loop
  581.          { # Hier a,b>0, beides NUDS.
  582.            # Vergleiche a und b:
  583.            if (a_len > b_len) goto a_greater_b; # a>b ?
  584.            if (a_len == b_len)
  585.              { var reg1 signean vergleich = compare_loop_up(a_MSDptr,b_MSDptr,a_len);
  586.                if (vergleich > 0) goto a_greater_b; # a>b ?
  587.                if (vergleich == 0) break; # a=b ?
  588.              }
  589.            # a<b -> a,b vertauschen:
  590.            swap(reg1 uintD*, a_MSDptr,b_MSDptr);
  591.            swap(reg1 uintC, a_len,b_len);
  592.            swap(reg1 uintD*, a_LSDptr,b_LSDptr);
  593.            a_greater_b:
  594.            # Hier a>b>0, beides NUDS.
  595.            # Entscheidung, ob Division oder Linearkombination:
  596.            { var reg1 uintD a_msd; # führende intDsize Bits von a
  597.              var reg2 uintD b_msd; # entsprechende Bits von b
  598.              { var reg4 uintC len_diff = a_len-b_len; # Längendifferenz
  599.                if (len_diff > 1) goto divide; # >=2 -> Bitlängendifferenz>intDsize -> dividieren
  600.                #define bitlendiff_limit  (intDsize/2) # sollte >0,<intDsize sein
  601.               {var reg3 uintC a_msd_size;
  602.                a_msd = a_MSDptr[0]; # führendes Digit von a
  603.                integerlengthD(a_msd,a_msd_size=); # dessen Bit-Länge (>0,<=intDsize) berechnen
  604.                b_msd = b_MSDptr[0];
  605.                #if HAVE_DD
  606.                {var reg5 uintDD b_msdd = # 2 führende Digits von b
  607.                   (len_diff==0
  608.                    ? highlowDD(b_msd, (b_len==1 ? 0 : b_MSDptr[1]))
  609.                    : (uintDD)b_msd
  610.                   );
  611.                 # a_msd_size+intDsize - b_msdd_size >= bitlendiff_limit -> dividieren:
  612.                 b_msdd = b_msdd >> a_msd_size;
  613.                 if (b_msdd < bit(intDsize-bitlendiff_limit)) goto divide;
  614.                 b_msd = lowD(b_msdd);
  615.                }
  616.                {var reg5 uintDD a_msdd = # 2 führende Digits von a
  617.                   highlowDD(a_msd, (a_len==1 ? 0 : a_MSDptr[1]));
  618.                 a_msd = lowD(a_msdd >> a_msd_size);
  619.                }
  620.                if (a_msd == b_msd) goto subtract;
  621.                #else
  622.                if (len_diff==0)
  623.                  { # a_msd_size - b_msd_size >= bitlendiff_limit -> dividieren:
  624.                    if ((a_msd_size > bitlendiff_limit)
  625.                        && (b_msd < bit(a_msd_size-bitlendiff_limit))
  626.                       )
  627.                      goto divide;
  628.                    # Entscheidung für Linearkombination ist gefallen.
  629.                    # a_msd und b_msd so erweitern, daß a_msd die führenden
  630.                    # intDsize Bits von a enthält:
  631.                   {var reg5 uintC shiftcount = intDsize-a_msd_size; # Shiftcount nach links (>=0, <intDsize)
  632.                    if (shiftcount>0)
  633.                      { a_msd = a_msd << shiftcount;
  634.                        b_msd = b_msd << shiftcount;
  635.                        if (a_len>1)
  636.                          { a_msd |= a_MSDptr[1] >> a_msd_size;
  637.                            b_msd |= b_MSDptr[1] >> a_msd_size;
  638.                      }   }
  639.                    if (a_msd == b_msd) goto subtract;
  640.                  }}
  641.                  else
  642.                  # len_diff=1
  643.                  { # a_msd_size+intDsize - b_msd_size >= bitlendiff_limit -> dividieren:
  644.                    if ((a_msd_size >= bitlendiff_limit)
  645.                        || (b_msd < bit(a_msd_size+intDsize-bitlendiff_limit))
  646.                       )
  647.                      goto divide;
  648.                    # Entscheidung für Linearkombination ist gefallen.
  649.                    # a_msd und b_msd so erweitern, daß a_msd die führenden
  650.                    # intDsize Bits von a enthält:
  651.                    # 0 < a_msd_size < b_msd_size + bitlendiff_limit - intDsize <= bitlendiff_limit < intDsize.
  652.                    a_msd = (a_msd << (intDsize-a_msd_size)) | (a_MSDptr[1] >> a_msd_size);
  653.                    b_msd = b_msd >> a_msd_size;
  654.                  }
  655.                #endif
  656.                #undef bitlendiff_limit
  657.              }}
  658.              # Nun ist a_msd = a' > b' = b_msd.
  659.              { # Euklid-Algorithmus auf den führenden Digits durchführen:
  660.                var partial_gcd_result likobi;
  661.                partial_gcd(a_msd,b_msd,&likobi); # liefert x1,y1,x2,y2
  662.                # Hier y1>0.
  663.                if (likobi.x2==0)
  664.                  { # Ersetze (a,b) := (a-y1*b,b).
  665.                    if (likobi.y1==1) goto subtract; # einfacherer Fall
  666.                    # Dazu evtl. a um 1 Digit erweitern, so daß a_len=b_len+1:
  667.                    if (a_len == b_len) { *--a_MSDptr = 0; a_len++; }
  668.                    # und y1*b von a subtrahieren:
  669.                    a_MSDptr[0] -= mulusub_loop_down(likobi.y1,b_LSDptr,a_LSDptr,b_len);
  670.                  }
  671.                  else
  672.                  { # Ersetze (a,b) := (x1*a-y1*b,-x2*a+y2*b).
  673.                    # Dazu evtl. b um 1 Digit erweitern, so daß a_len=b_len:
  674.                    if (!(a_len==b_len)) { *--b_MSDptr = 0; b_len++; }
  675.                    # c := x1*a-y1*b bilden:
  676.                    mulu_loop_down(likobi.x1,a_LSDptr,c_LSDptr,a_len);
  677.                    /* c_LSDptr[-(uintP)a_len-1] -= */
  678.                      mulusub_loop_down(likobi.y1,b_LSDptr,c_LSDptr,a_len);
  679.                    # d := -x2*a+y2*b bilden:
  680.                    mulu_loop_down(likobi.y2,b_LSDptr,d_LSDptr,a_len);
  681.                    /* d_LSDptr[-(uintP)a_len-1] -= */
  682.                      mulusub_loop_down(likobi.x2,a_LSDptr,d_LSDptr,a_len);
  683.                    # Wir wissen, daß 0 < c < b und 0 < d < a. Daher müßten
  684.                    # c_LSDptr[-a_len-1] und d_LSDptr[-a_len-1] =0 sein.
  685.                    # a := c und b := d kopieren:
  686.                    copy_loop_down(c_LSDptr,a_LSDptr,a_len);
  687.                    copy_loop_down(d_LSDptr,b_LSDptr,a_len);
  688.                    # b normalisieren:
  689.                    while (b_MSDptr[0]==0) { b_MSDptr++; b_len--; }
  690.              }   }
  691.              if (FALSE)
  692.                { subtract: # Ersetze (a,b) := (a-b,b).
  693.                  if (!( subfrom_loop_down(b_LSDptr,a_LSDptr,b_len) ==0))
  694.                    # Übertrag nach b_len Stellen, muß also a_len=b_len+1 sein.
  695.                    { a_MSDptr[0] -= 1; }
  696.                }
  697.              # a normalisieren:
  698.              while (a_MSDptr[0]==0) { a_MSDptr++; a_len--; }
  699.            }
  700.            if (FALSE)
  701.              { divide: # Ersetze (a,b) := (b , a mod b).
  702.               {var reg10 uintD* old_a_LSDptr = a_LSDptr;
  703.                var DS q;
  704.                var DS r;
  705.                UDS_divide_(a_MSDptr,a_len,a_LSDptr,b_MSDptr,b_len,b_LSDptr, divroomptr, &q,&r);
  706.                a_MSDptr = b_MSDptr; a_len = b_len; a_LSDptr = b_LSDptr; # a := b
  707.                b_len = r.len; if (b_len==0) break; # b=0 -> fertig
  708.                b_LSDptr = old_a_LSDptr; # b übernimmt den vorherigen Platz von a
  709.                b_MSDptr = copy_loop_down(r.LSDptr,b_LSDptr,b_len); # b := r kopieren
  710.                goto a_greater_b; # Nun ist a>b>0
  711.              }}
  712.       }  }
  713.       RESTORE_NUM_STACK # num_stack (vorzeitig) zurück
  714.       return NUDS_to_I(a_MSDptr,a_len); # NUDS a als Ergebnis
  715.       #undef I_abs_to_NUDS
  716.     }}
  717. #endif
  718. #if XGCD_ALGO==3
  719. # (xgcd A B) :==
  720. # wie oben bei (gcd A B).
  721. # Zusätzlich werden Variablen sA,sB,sk,uAa,uBa,uAb,uBb geführt,
  722. # wobei sA,sB,sk Vorzeichen (+/- 1) und uAa,uBa,uAb,uBb Integers >=0 sind mit
  723. #     uAa * sA*A - uBa * sB*B = a,
  724. #   - uAb * sA*A + uBb * sB*B = b,
  725. # ferner  uAa * uBb - uAb * uBa = sk  und daher (Cramersche Regel)
  726. #   uBb * a + uBa * b = sk*sA*A, uAb * a + uAa * b = sk*sB*B.
  727. # Zu Beginn (a,b) := (|A|,|B|), (sA,sB) := ((signum A), (signumB)),
  728. #           (uAa,uBa,uAb,uBb) := (1,0,0,1).
  729. # Beim Ersetzen (a,b) := (a-b,b)
  730. #   ersetzt man (uAa,uBa,uAb,uBb) := (uAa+uAb,uBa+uBb,uAb,uBb).
  731. # Beim Ersetzen (a,b) := (a-y1*b,b)
  732. #   ersetzt man (uAa,uBa,uAb,uBb) := (uAa+y1*uAb,uBa+y1*uBb,uAb,uBb).
  733. # Beim Ersetzen (a,b) := (x1*a-y1*b,-x2*a+y2*b) mit x1*y2-x2*y1=1
  734. #   ersetzt man (uAa,uBa,uAb,uBb) :=
  735. #               (x1*uAa+y1*uAb,x1*uBa+y1*uBb,x2*uAa+y2*uAb,x2*uBa+y2*uBb).
  736. # Beim Ersetzen (a,b) := (b,a)
  737. #   ersetzt man (uAa,uBa,uAb,uBb) := (uAb,uBb,uAa,uBa),
  738. #               sk := -sk, (sA,sB) := (-sA,-sB).
  739. # Beim Ersetzen (a,b) := (b,a-q*b)
  740. #   ersetzt man (uAa,uBa,uAb,uBb) := (uAb,uBb,uAa+q*uAb,uBa+q*uBb),
  741. #               sk := -sk, (sA,sB) := (-sA,-sB).
  742. # Zum Schluß ist a der ggT und a = uAa*sA * A + -uBa*sB * B
  743. # die gewünschte Linearkombination.
  744. # Da stets gilt sk*sA*A = |A|, sk*sB*B = |B|, a>=1, b>=1,
  745. # folgt 0 <= uAa <= |B|, 0 <= uAb <= |B|, 0 <= uBa <= |A|, 0 <= uBb <= |A|.
  746. # Ferner wird sk nie benutzt, braucht also nicht mitgeführt zu werden.
  747.   # Bildet u := u + v, wobei für u genügend Platz sei:
  748.   # (Benutzt v.LSDptr nicht.)
  749.   local void NUDS_likobi0_NUDS (DS* u, DS* v);
  750.   local void NUDS_likobi0_NUDS(u,v)
  751.     var reg1 DS* u;
  752.     var reg2 DS* v;
  753.     { var reg3 uintC u_len = u->len;
  754.       var reg4 uintC v_len = v->len;
  755.       if (u_len >= v_len)
  756.         { if (!( addto_loop_down(v->LSDptr,u->LSDptr,v_len) ==0))
  757.             { if (!( inc_loop_down(&(u->LSDptr)[-(uintP)v_len],u_len-v_len) ==0))
  758.                 { *--(u->MSDptr) = 1; u->len++; }
  759.         }   }
  760.         else # u_len <= v_len
  761.         { u->MSDptr = copy_loop_down(&(v->LSDptr)[-(uintP)u_len],&(u->LSDptr)[-(uintP)u_len],v_len-u_len);
  762.           u->len = v_len;
  763.           if (!( addto_loop_down(v->LSDptr,u->LSDptr,u_len) ==0))
  764.             { if (!( inc_loop_down(&(u->LSDptr)[-(uintP)u_len],v_len-u_len) ==0))
  765.                 { *--(u->MSDptr) = 1; u->len++; }
  766.         }   }
  767.     }
  768.   # Bildet u := u + q*v, wobei für u genügend Platz sei:
  769.   # (Dabei sei nachher u>0.)
  770.   local void NUDS_likobi1_NUDS (DS* u, DS* v, uintD q);
  771.   local void NUDS_likobi1_NUDS(u,v,q)
  772.     var reg2 DS* u;
  773.     var reg1 DS* v;
  774.     var reg6 uintD q;
  775.     { var reg3 uintC v_len = v->len;
  776.       if (v_len>0) # nur nötig, falls v /=0
  777.         { var reg4 uintC u_len = u->len;
  778.           var reg5 uintD carry;
  779.           if (u_len <= v_len) # evtl. u vergrößern
  780.             { u->MSDptr = clear_loop_down(u->MSDptr,v_len-u_len+1);
  781.               u->len = u_len = v_len+1;
  782.             } # Nun ist u_len > v_len.
  783.           carry = muluadd_loop_down(q,v->LSDptr,u->LSDptr,v_len);
  784.           if (!(carry==0))
  785.             { var reg5 uintD* ptr = &(u->LSDptr)[-(uintP)v_len-1];
  786.               if ((ptr[0] += carry) < carry)
  787.                 { if (!( inc_loop_down(ptr,u_len-v_len-1) ==0))
  788.                     { *--(u->MSDptr) = 1; u->len++; }
  789.             }   }
  790.           while ((u->MSDptr)[0]==0) { (u->MSDptr)++; u->len--; } # normalisieren
  791.     }   }
  792.   # Bildet (u,v) := (x1*u+y1*v,x2*u+y2*v), wobei für u,v genügend Platz sei:
  793.   # (Dabei sei u>0 oder v>0, nachher u>0 und v>0.)
  794.   local void NUDS_likobi2_NUDS (DS* u, DS* v, partial_gcd_result* q, uintD* c_LSDptr, uintD* d_LSDptr);
  795.   local void NUDS_likobi2_NUDS(u,v,q,c_LSDptr,d_LSDptr)
  796.     var reg1 DS* u;
  797.     var reg2 DS* v;
  798.     var reg3 partial_gcd_result* q;
  799.     var reg9 uintD* c_LSDptr;
  800.     var reg10 uintD* d_LSDptr;
  801.     { var reg4 uintC u_len = u->len;
  802.       var reg5 uintC v_len = v->len;
  803.       var reg6 uintC c_len;
  804.       var reg7 uintC d_len;
  805.       if (u_len >= v_len)
  806.         { mulu_loop_down(q->x1,u->LSDptr,c_LSDptr,u_len); c_len = u_len+1;
  807.           mulu_loop_down(q->x2,u->LSDptr,d_LSDptr,u_len); d_len = u_len+1;
  808.           if (!(v_len==0))
  809.             {{var reg6 uintD carry =
  810.                 muluadd_loop_down(q->y1,v->LSDptr,c_LSDptr,v_len);
  811.               if (!(carry==0))
  812.                 { var reg1 uintD* ptr = &c_LSDptr[-(uintP)v_len-1];
  813.                   if ((ptr[0] += carry) < carry)
  814.                     { if (!( inc_loop_down(ptr,u_len-v_len) ==0))
  815.                         { c_len++; c_LSDptr[-(uintP)c_len] = 1; }
  816.              }  }   }
  817.              {var reg6 uintD carry =
  818.                 muluadd_loop_down(q->y2,v->LSDptr,d_LSDptr,v_len);
  819.               if (!(carry==0))
  820.                 { var reg1 uintD* ptr = &d_LSDptr[-(uintP)v_len-1];
  821.                   if ((ptr[0] += carry) < carry)
  822.                     { if (!(inc_loop_down(ptr,u_len-v_len) ==0))
  823.                         { d_len++; d_LSDptr[-(uintP)d_len] = 1; }
  824.             }}  }   }
  825.         }
  826.         else
  827.         { mulu_loop_down(q->y1,v->LSDptr,c_LSDptr,v_len); c_len = v_len+1;
  828.           mulu_loop_down(q->y2,v->LSDptr,d_LSDptr,v_len); d_len = v_len+1;
  829.           if (!(u_len==0))
  830.             {{var reg6 uintD carry =
  831.                 muluadd_loop_down(q->x1,u->LSDptr,c_LSDptr,u_len);
  832.               if (!(carry==0))
  833.                 { var reg1 uintD* ptr = &c_LSDptr[-(uintP)u_len-1];
  834.                   if ((ptr[0] += carry) < carry)
  835.                     { if (!( inc_loop_down(ptr,v_len-u_len) ==0))
  836.                         { c_len++; c_LSDptr[-(uintP)c_len] = 1; }
  837.              }  }   }
  838.              {var reg6 uintD carry =
  839.                 muluadd_loop_down(q->x2,u->LSDptr,d_LSDptr,u_len);
  840.               if (!(carry==0))
  841.                 { var reg1 uintD* ptr = &d_LSDptr[-(uintP)u_len-1];
  842.                   if ((ptr[0] += carry) < carry)
  843.                     { if (!( inc_loop_down(ptr,v_len-u_len) ==0))
  844.                         { d_len++; d_LSDptr[-(uintP)d_len] = 1; }
  845.             }}  }   }
  846.         }
  847.       u->MSDptr = copy_loop_down(c_LSDptr,u->LSDptr,c_len);
  848.       while ((u->MSDptr)[0]==0) { u->MSDptr++; c_len--; }
  849.       u->len = c_len;
  850.       v->MSDptr = copy_loop_down(d_LSDptr,v->LSDptr,d_len);
  851.       while ((v->MSDptr)[0]==0) { v->MSDptr++; d_len--; }
  852.       v->len = d_len;
  853.     }
  854.   # Los geht's:
  855.   local void I_I_xgcd_I_I_I(a,b)
  856.     var reg10 object a;
  857.     var reg10 object b;
  858.     { if (eq(a,Fixnum_1)) # a=1 -> g=1, (u,v)=(1,0)
  859.         { pushSTACK(Fixnum_1); pushSTACK(Fixnum_0);
  860.           pushSTACK(a); return;
  861.         }
  862.       if (eq(b,Fixnum_1)) # b=1 -> g=1, (u,v)=(0,1)
  863.         { pushSTACK(Fixnum_0); pushSTACK(Fixnum_1);
  864.           pushSTACK(b); return;
  865.         }
  866.      {var reg10 boolean sA = (R_minusp(a) ? ~0 : 0); # Vorzeichen von A
  867.       var reg10 boolean sB = (R_minusp(b) ? ~0 : 0); # Vorzeichen von B
  868.       SAVE_NUM_STACK # num_stack retten
  869.       var reg6 uintD* a_MSDptr;
  870.       var reg8 uintC a_len;
  871.       var reg10 uintD* a_LSDptr;
  872.       var reg7 uintD* b_MSDptr;
  873.       var reg9 uintC b_len;
  874.       var reg10 uintD* b_LSDptr;
  875.       # Macro: erzeugt die NUDS zu (abs x), erniedrigt num_stack
  876.       #define I_abs_to_NUDS(x,zero_statement)  \
  877.         { I_to_NDS_1(x, x##_MSDptr = , x##_len = , x##_LSDptr = ); # (nichtleere) NDS holen \
  878.           if (x##_len == 0) { zero_statement } # falls =0, fertig                           \
  879.           if ((sintD)x##_MSDptr[0] < 0) # falls <0, negieren:                               \
  880.             { neg_loop_down(x##_LSDptr,x##_len); }                                          \
  881.           if (x##_MSDptr[0] == 0) # normalisieren (max. 1 Nulldigit entfernen)              \
  882.             { x##_MSDptr++; x##_len--; }                                                    \
  883.         }
  884.       I_abs_to_NUDS(a, # (abs A) als NUDS erzeugen
  885.                     # A=0 -> g=|B|, (u,v) = (0,sB)
  886.                     { pushSTACK(Fixnum_0); # u
  887.                       pushSTACK(sB==0 ? Fixnum_1 : Fixnum_minus1); # v
  888.                       pushSTACK(I_abs_I(b)); # g
  889.                       return;
  890.                     });
  891.       I_abs_to_NUDS(b, # (abs B) als NUDS erzeugen
  892.                     # B=0 -> g=|A|, (u,v) = (sA,0)
  893.                     { pushSTACK(sA==0 ? Fixnum_1 : Fixnum_minus1); # u
  894.                       pushSTACK(Fixnum_0); # v
  895.                       pushSTACK(I_abs_I(a)); # g
  896.                       return;
  897.                     });
  898.       # Jetzt ist a = a_MSDptr/a_len/a_LSDptr, b = b_MSDptr/b_len/b_LSDptr,
  899.       # beides NUDS, und a_len>0, b_len>0.
  900.       {# Beifaktoren:
  901.        var DS uAa;
  902.        var DS uBa;
  903.        var DS uAb;
  904.        var DS uBb;
  905.        # Rechenregister:
  906.        var reg10 uintD* divroomptr; # Platz für Divisionsergebnis
  907.        var reg10 uintD* c_LSDptr;
  908.        var reg10 uintD* d_LSDptr;
  909.        # Platz für uAa,uBa,uAb,uBb besorgen:
  910.        {var reg1 uintC u_len = b_len+1;
  911.         num_stack_need(u_len,_EMA_,uAa.LSDptr=); uAa.MSDptr = uAa.LSDptr;
  912.         num_stack_need(u_len,_EMA_,uAb.LSDptr=); uAb.MSDptr = uAb.LSDptr;
  913.        }
  914.        {var reg1 uintC u_len = a_len+1;
  915.         num_stack_need(u_len,_EMA_,uBa.LSDptr=); uBa.MSDptr = uBa.LSDptr;
  916.         num_stack_need(u_len,_EMA_,uBb.LSDptr=); uBb.MSDptr = uBb.LSDptr;
  917.        }
  918.        *--uAa.MSDptr = 1; uAa.len = 1; # uAa := 1
  919.        uBa.len = 0; # uBa := 0
  920.        uAb.len = 0; # uAb := 0
  921.        *--uBb.MSDptr = 1; uBb.len = 1; # uBb := 1
  922.        # Jetzt ist uAa = uAa.MSDptr/uAa.len/uAa.LSDptr,
  923.        #           uBa = uBa.MSDptr/uBa.len/uBa.LSDptr,
  924.        #           uAb = uAb.MSDptr/uAb.len/uAb.LSDptr,
  925.        #           uBb = uBb.MSDptr/uBb.len/uBb.LSDptr,
  926.        # alles NUDS.
  927.        # Platz für zwei Rechenregister besorgen, mit je max(a_len,b_len)+1 Digits:
  928.        {var reg10 uintL c_len = (uintL)(a_len>=b_len ? a_len : b_len) + 1;
  929.         num_stack_need(c_len,_EMA_,c_LSDptr=);
  930.         num_stack_need(c_len,divroomptr=,d_LSDptr=);
  931.         # Jetzt ist ../c_len/c_LSDptr, ../c_len/d_LSDptr frei.
  932.        }
  933.        loop
  934.          { # Hier a,b>0, beides NUDS.
  935.            # Vergleiche a und b:
  936.            if (a_len > b_len) goto a_greater_b; # a>b ?
  937.            if (a_len == b_len)
  938.              { var reg1 signean vergleich = compare_loop_up(a_MSDptr,b_MSDptr,a_len);
  939.                if (vergleich > 0) goto a_greater_b; # a>b ?
  940.                if (vergleich == 0) break; # a=b ?
  941.              }
  942.            # a<b -> a,b vertauschen:
  943.            swap(reg1 uintD*, a_MSDptr,b_MSDptr);
  944.            swap(reg1 uintC, a_len,b_len);
  945.            swap(reg1 uintD*, a_LSDptr,b_LSDptr);
  946.            a_greater_b_swap:
  947.            swap(DS, uAa,uAb); # und uAa und uAb vertauschen
  948.            swap(DS, uBa,uBb); # und uBa und uBb vertauschen
  949.            sA = ~sA; sB = ~sB; # und sA und sB umdrehen
  950.            a_greater_b:
  951.            # Hier a>b>0, beides NUDS.
  952.            # Entscheidung, ob Division oder Linearkombination:
  953.            { var reg1 uintD a_msd; # führende intDsize Bits von a
  954.              var reg2 uintD b_msd; # entsprechende Bits von b
  955.              { var reg4 uintC len_diff = a_len-b_len; # Längendifferenz
  956.                if (len_diff > 1) goto divide; # >=2 -> Bitlängendifferenz>intDsize -> dividieren
  957.                #define bitlendiff_limit  (intDsize/2) # sollte >0,<intDsize sein
  958.               {var reg3 uintC a_msd_size;
  959.                a_msd = a_MSDptr[0]; # führendes Digit von a
  960.                integerlengthD(a_msd,a_msd_size=); # dessen Bit-Länge (>0,<=intDsize) berechnen
  961.                b_msd = b_MSDptr[0];
  962.                #if HAVE_DD
  963.                {var reg5 uintDD b_msdd = # 2 führende Digits von b
  964.                   (len_diff==0
  965.                    ? highlowDD(b_msd, (b_len==1 ? 0 : b_MSDptr[1]))
  966.                    : (uintDD)b_msd
  967.                   );
  968.                 # a_msd_size+intDsize - b_msdd_size >= bitlendiff_limit -> dividieren:
  969.                 b_msdd = b_msdd >> a_msd_size;
  970.                 if (b_msdd < bit(intDsize-bitlendiff_limit)) goto divide;
  971.                 b_msd = lowD(b_msdd);
  972.                }
  973.                {var reg5 uintDD a_msdd = # 2 führende Digits von a
  974.                   highlowDD(a_msd, (a_len==1 ? 0 : a_MSDptr[1]));
  975.                 a_msd = lowD(a_msdd >> a_msd_size);
  976.                }
  977.                if (a_msd == b_msd) goto subtract;
  978.                #else
  979.                if (len_diff==0)
  980.                  { # a_msd_size - b_msd_size >= bitlendiff_limit -> dividieren:
  981.                    if ((a_msd_size > bitlendiff_limit)
  982.                        && (b_msd < bit(a_msd_size-bitlendiff_limit))
  983.                       )
  984.                      goto divide;
  985.                    # Entscheidung für Linearkombination ist gefallen.
  986.                    # a_msd und b_msd so erweitern, daß a_msd die führenden
  987.                    # intDsize Bits von a enthält:
  988.                   {var reg5 uintC shiftcount = intDsize-a_msd_size; # Shiftcount nach links (>=0, <intDsize)
  989.                    if (shiftcount>0)
  990.                      { a_msd = a_msd << shiftcount;
  991.                        b_msd = b_msd << shiftcount;
  992.                        if (a_len>1)
  993.                          { a_msd |= a_MSDptr[1] >> a_msd_size;
  994.                            b_msd |= b_MSDptr[1] >> a_msd_size;
  995.                      }   }
  996.                    if (a_msd == b_msd) goto subtract;
  997.                  }}
  998.                  else
  999.                  # len_diff=1
  1000.                  { # a_msd_size+intDsize - b_msd_size >= bitlendiff_limit -> dividieren:
  1001.                    if ((a_msd_size >= bitlendiff_limit)
  1002.                        || (b_msd < bit(a_msd_size+intDsize-bitlendiff_limit))
  1003.                       )
  1004.                      goto divide;
  1005.                    # Entscheidung für Linearkombination ist gefallen.
  1006.                    # a_msd und b_msd so erweitern, daß a_msd die führenden
  1007.                    # intDsize Bits von a enthält:
  1008.                    # 0 < a_msd_size < b_msd_size + bitlendiff_limit - intDsize <= bitlendiff_limit < intDsize.
  1009.                    a_msd = (a_msd << (intDsize-a_msd_size)) | (a_MSDptr[1] >> a_msd_size);
  1010.                    b_msd = b_msd >> a_msd_size;
  1011.                  }
  1012.                #endif
  1013.                #undef bitlendiff_limit
  1014.              }}
  1015.              # Nun ist a_msd = a' > b' = b_msd.
  1016.              { # Euklid-Algorithmus auf den führenden Digits durchführen:
  1017.                var partial_gcd_result likobi;
  1018.                partial_gcd(a_msd,b_msd,&likobi); # liefert x1,y1,x2,y2
  1019.                # Hier y1>0.
  1020.                if (likobi.x2==0)
  1021.                  { # Ersetze (a,b) := (a-y1*b,b).
  1022.                    if (likobi.y1==1) goto subtract; # einfacherer Fall
  1023.                    # Dazu evtl. a um 1 Digit erweitern, so daß a_len=b_len+1:
  1024.                    if (a_len == b_len) { *--a_MSDptr = 0; a_len++; }
  1025.                    # und y1*b von a subtrahieren:
  1026.                    a_MSDptr[0] -= mulusub_loop_down(likobi.y1,b_LSDptr,a_LSDptr,b_len);
  1027.                    NUDS_likobi1_NUDS(&uAa,&uAb,likobi.y1); # uAa := uAa + y1 * uAb
  1028.                    NUDS_likobi1_NUDS(&uBa,&uBb,likobi.y1); # uBa := uBa + y1 * uBb
  1029.                  }
  1030.                  else
  1031.                  { # Ersetze (uAa,uAb) := (x1*uAa+y1*uAb,x2*uAa+y2*uAb) :
  1032.                    NUDS_likobi2_NUDS(&uAa,&uAb,&likobi,c_LSDptr,d_LSDptr);
  1033.                    # Ersetze (uBa,uBb) := (x1*uBa+y1*uBb,x2*uBa+y2*uBb) :
  1034.                    NUDS_likobi2_NUDS(&uBa,&uBb,&likobi,c_LSDptr,d_LSDptr);
  1035.                    # Ersetze (a,b) := (x1*a-y1*b,-x2*a+y2*b).
  1036.                    # Dazu evtl. b um 1 Digit erweitern, so daß a_len=b_len:
  1037.                    if (!(a_len==b_len)) { *--b_MSDptr = 0; b_len++; }
  1038.                    # c := x1*a-y1*b bilden:
  1039.                    mulu_loop_down(likobi.x1,a_LSDptr,c_LSDptr,a_len);
  1040.                    /* c_LSDptr[-(uintP)a_len-1] -= */
  1041.                      mulusub_loop_down(likobi.y1,b_LSDptr,c_LSDptr,a_len);
  1042.                    # d := -x2*a+y2*b bilden:
  1043.                    mulu_loop_down(likobi.y2,b_LSDptr,d_LSDptr,a_len);
  1044.                    /* d_LSDptr[-(uintP)a_len-1] -= */
  1045.                      mulusub_loop_down(likobi.x2,a_LSDptr,d_LSDptr,a_len);
  1046.                    # Wir wissen, daß 0 < c < b und 0 < d < a. Daher müßten
  1047.                    # c_LSDptr[-a_len-1] und d_LSDptr[-a_len-1] =0 sein.
  1048.                    # a := c und b := d kopieren:
  1049.                    copy_loop_down(c_LSDptr,a_LSDptr,a_len);
  1050.                    copy_loop_down(d_LSDptr,b_LSDptr,a_len);
  1051.                    # b normalisieren:
  1052.                    while (b_MSDptr[0]==0) { b_MSDptr++; b_len--; }
  1053.              }   }
  1054.              if (FALSE)
  1055.                { subtract: # Ersetze (a,b) := (a-b,b).
  1056.                  NUDS_likobi0_NUDS(&uAa,&uAb); # uAa := uAa + uAb
  1057.                  NUDS_likobi0_NUDS(&uBa,&uBb); # uBa := uBa + uBb
  1058.                  if (!( subfrom_loop_down(b_LSDptr,a_LSDptr,b_len) ==0))
  1059.                    # Übertrag nach b_len Stellen, muß also a_len=b_len+1 sein.
  1060.                    { a_MSDptr[0] -= 1; }
  1061.                }
  1062.              # a normalisieren:
  1063.              while (a_MSDptr[0]==0) { a_MSDptr++; a_len--; }
  1064.            }
  1065.            if (FALSE)
  1066.              { divide: # Ersetze (a,b) := (b , a mod b).
  1067.               {var reg10 uintD* old_a_LSDptr = a_LSDptr;
  1068.                var DS q;
  1069.                var DS r;
  1070.                UDS_divide_(a_MSDptr,a_len,a_LSDptr,b_MSDptr,b_len,b_LSDptr, divroomptr, &q,&r);
  1071.                a_MSDptr = b_MSDptr; a_len = b_len; a_LSDptr = b_LSDptr; # a := b
  1072.                b_len = r.len; if (b_len==0) break; # b=0 -> fertig
  1073.                b_LSDptr = old_a_LSDptr; # b übernimmt den vorherigen Platz von a
  1074.                b_MSDptr = copy_loop_down(r.LSDptr,b_LSDptr,b_len); # b := r kopieren
  1075.                # (uAa,uAb) := (uAb,uAa+q*uAb) :
  1076.                if (!(uAb.len==0))
  1077.                  { mulu_2loop_down(q.LSDptr,q.len,uAb.LSDptr,uAb.len,c_LSDptr); # q * uAb
  1078.                   {var DS c;
  1079.                    c.LSDptr = c_LSDptr; c.len = q.len + uAb.len;
  1080.                    if (c_LSDptr[-(uintP)c.len]==0) { c.len--; } # normalisieren
  1081.                    NUDS_likobi0_NUDS(&uAa,&c); # zu uAa addieren
  1082.                  }} # noch uAa,uAb vertauschen (später)
  1083.                # (uBa,uBb) := (uBb,uBa+q*uBb) :
  1084.                if (!(uBb.len==0))
  1085.                  { mulu_2loop_down(q.LSDptr,q.len,uBb.LSDptr,uBb.len,c_LSDptr); # q * uBb
  1086.                   {var DS c;
  1087.                    c.LSDptr = c_LSDptr; c.len = q.len + uBb.len;
  1088.                    if (c_LSDptr[-(uintP)c.len]==0) { c.len--; } # normalisieren
  1089.                    NUDS_likobi0_NUDS(&uBa,&c); # zu uBa addieren
  1090.                  }} # noch uBa,uBb vertauschen (später)
  1091.                goto a_greater_b_swap; # Nun ist a>b>0
  1092.              }}
  1093.          }
  1094.        # uAb mit Vorfaktor -sA versehen:
  1095.        *--uAb.MSDptr = 0; uAb.len++;
  1096.        if (sA==0) { neg_loop_down(uAb.LSDptr,uAb.len); }
  1097.        # uBb mit Vorfaktor sB versehen:
  1098.        *--uBb.MSDptr = 0; uBb.len++;
  1099.        if (!(sB==0)) { neg_loop_down(uBb.LSDptr,uBb.len); }
  1100.        pushSTACK(DS_to_I(uAb.MSDptr,uAb.len)); # DS uAb als Vorfaktor von A
  1101.        pushSTACK(DS_to_I(uBb.MSDptr,uBb.len)); # DS uBb als Vorfaktor von B
  1102.       }
  1103.       pushSTACK(NUDS_to_I(a_MSDptr,a_len)); # NUDS a als ggT
  1104.       RESTORE_NUM_STACK # num_stack zurück
  1105.       #undef I_abs_to_NUDS
  1106.     }}
  1107. #endif
  1108.  
  1109. # Liefert das kgV zweier Integers.
  1110. # I_I_lcm_I(a,b)
  1111. # > a,b: zwei Integers
  1112. # < ergebnis: (lcm a b), ein Integer >=0
  1113. # kann GC auslösen
  1114.   local object I_I_lcm_I (object a, object b);
  1115. # Methode:
  1116. # a=0 oder b=0 -> Ergebnis 0.
  1117. # a:=(abs a), b:=(abs b).
  1118. # g:=ggT(a,b)>0.
  1119. # Falls g=1, Ergebnis a*b, sonst Ergebnis (a/g)*b.
  1120.   local object I_I_lcm_I(a,b)
  1121.     var reg1 object a;
  1122.     var reg3 object b;
  1123.     { if (eq(a,Fixnum_0)) { return a; }
  1124.       if (eq(b,Fixnum_0)) { return b; }
  1125.       # Beträge nehmen:
  1126.       pushSTACK(b); pushSTACK(I_abs_I(a)); STACK_1 = b = I_abs_I(STACK_1);
  1127.       # Stackaufbau: b := (abs b), a := (abs a).
  1128.      {var reg2 object g = I_I_gcd_I(STACK_0,b); # g = (gcd a b)
  1129.       a = popSTACK();
  1130.       if (!eq(g,Fixnum_1)) { a = I_I_exquopos_I(a,g); } # a durch g (beide >0) dividieren
  1131.       return I_I_mal_I(a,popSTACK()); # mit b multiplizieren
  1132.     }}
  1133.  
  1134.