4 #include "base/cl_sysdep.h"
7 #include "cln/integer.h"
12 #include "integer/cl_I.h"
13 #include "base/digitseq/cl_DS.h"
14 #include "base/digit/cl_D.h"
15 #include "base/cl_xmacros.h"
19 #define GCD_ALGO 3 // 1: binär, 2: Schulmethode, 3: Lehmer
29 // (abs a) und (abs b) in zwei Buffer packen, als Unsigned Digit Sequences.
30 // [Schreibe oBdA wieder a,b]
35 // (progn (incf j) (setq a (/ a 2)) (setq b (/ b 2)) (go 1))
38 // (while (evenp b) (setq b (/ b 2)))
40 // 2 {a,b >0, beide ungerade}
43 // ((< a b) (rotatef a b))
45 // 3 {a,b >0, beide ungerade, a>b}
47 // 4 {a,b >0, a gerade, b ungerade}
48 // (repeat (setq a (/ a 2)) (until (oddp a)))
53 // weil es oft auftritt (insbesondere bei GCD's mehrerer Zahlen):
56 const cl_I gcd (const cl_I& a, const cl_I& b)
57 { if (eq(a,1)) { return 1; } // a=1 -> 1
58 if (eq(b,1)) { return 1; } // b=1 -> 1
59 if (eq(b,0)) { return abs(a); } // b=0 -> (abs a)
60 if (eq(a,0)) { return abs(b); } // a=0 -> (abs b)
68 // Macro: erzeugt die NUDS zu (abs x), erniedrigt num_stack
69 #define I_abs_to_NUDS(x) \
70 I_to_NDS_1(x, x##_MSDptr = , x##_len = , x##_LSDptr = ); /* (nichtleere) NDS holen */\
71 if ((sintD)mspref(x##_MSDptr,0) < 0) /* falls <0, negieren: */\
72 { neg_loop_lsp(x##_LSDptr,x##_len); } \
73 if (mspref(x##_MSDptr,0) == 0) /* normalisieren (max. 1 Nulldigit entfernen) */\
74 { msshrink(x##_MSDptr); x##_len--; }
75 I_abs_to_NUDS(a); // (abs a) als NUDS erzeugen
76 I_abs_to_NUDS(b); // (abs b) als NUDS erzeugen
77 // Jetzt ist a = a_MSDptr/a_len/a_LSDptr, b = b_MSDptr/b_len/b_LSDptr,
78 // beides NUDS, und a_len>0, b_len>0.
81 { shift1right_loop_msp(x##_MSDptr,x##_len,0); /* um 1 Bit rechts schieben */ \
82 if (mspref(x##_MSDptr,0) == 0) { msshrink(x##_MSDptr); x##_len--; } /* normalisieren */\
84 // Macro: Ob x gerade ist.
86 ((lspref(x##_LSDptr,0) & bit(0)) ==0)
91 { j++; halb(a); halb(b); goto label_1; }
95 while (evenp(b)) { halb(b); }
96 label_2: // a,b >0, beide ungerade
97 // Vergleiche a und b:
98 if (a_len > b_len) goto label_3; // a>b ?
100 { var cl_signean vergleich = compare_loop_msp(a_MSDptr,b_MSDptr,a_len);
101 if (vergleich > 0) goto label_3; // a>b ?
102 if (vergleich == 0) goto label_5; // a=b ?
104 // a<b -> a,b vertauschen:
105 swap(uintD*, a_MSDptr,b_MSDptr);
106 swap(uintC, a_len,b_len);
107 swap(uintD*, a_LSDptr,b_LSDptr);
108 label_3: // a,b >0, beide ungerade, a>b
109 // subtrahiere a := a - b
110 if (!( subfrom_loop_lsp(b_LSDptr,a_LSDptr,b_len) ==0))
111 { dec_loop_lsp(a_LSDptr lspop b_len,a_len-b_len); }
112 while (mspref(a_MSDptr,0) == 0) { msshrink(a_MSDptr); a_len--; } // normalisieren
113 label_4: // a,b >0, a gerade, b ungerade
114 do { halb(a); } while (evenp(a));
117 // a zu einer NDS machen:
118 return ash(NUDS_to_I(a_MSDptr,a_len),j); // ggT der ungeraden Anteile als Integer, mal 2^j
125 #endif /* GCD_ALGO == 1 */
132 // [a:=(abs a), b:=(abs b), while b>0 do (a,b) := (b,(mod a b)), -> a]
138 // a:=(abs a), b:=(abs b)
139 // Falls a=b: return a; falls a<b: vertausche a und b.
141 // Falls b=1, return 1. {spart eine Division durch 1}
142 // Sonst dividieren (divide a b), a:=b, b:=Rest.
143 // Falls b=0, return a, sonst goto (*).
145 const cl_I gcd (const cl_I& a, const cl_I& b)
146 { if (eq(a,1)) { return 1; } // a=1 -> 1
147 if (eq(b,1)) { return 1; } // b=1 -> 1
148 if (eq(b,0)) { return abs(a); } // b=0 -> (abs a)
149 if (eq(a,0)) { return abs(b); } // a=0 -> (abs b)
151 {var cl_I abs_a = abs(a);
152 var cl_I abs_b = abs(b);
155 if (fixnump(a) && fixnump(b)) // ggT zweier Fixnums >0
156 { // bleibt Fixnum, da (gcd a b) <= (min a b)
157 return V_to_FN(gcd(FN_to_UV(a),FN_to_UV(b)));
159 { var cl_signean vergleich = compare(a,b);
160 if (vergleich == 0) { return a; } // a=b -> fertig
161 if (vergleich < 0) { var cl_I tmp = a; a = b; b = a; } // a<b -> a,b vertauschen
163 loop // Hier a > b > 0
164 { if (eq(b,1)) { return 1; } // b=1 -> Ergebnis 1
165 { var cl_I_div_t div = cl_divide(a,b); a = b; b = div.remainder; }
166 if (eq(b,0)) { return a; }
170 #endif /* GCD_ALGO == 2 */
176 // vgl. [ D. E. Knuth: The Art of Computer Programming, Vol. 2: Seminumerical
177 // Algorithms, Sect. 4.5.2., Algorithm L ]
178 // und [ Collins, Loos: SAC-2, Algorithms IGCD, DPCC ].
184 // a:=(abs a), b:=(abs b)
186 // Falls a=b: return a; falls a<b: vertausche a und b.
188 // Falls (- (integer-length a) (integer-length b)) >= intDsize/2,
189 // lohnt sich eine Division: (a,b) := (b , a mod b). Falls b=0: return a.
190 // Falls dagegen 0 <= (- (integer-length a) (integer-length b)) < intDsize/2,
191 // seien a' die führenden intDsize Bits von a
192 // (2^(intDsize-1) <= a' < 2^intDsize) und b' die entsprechenden Bits von b
193 // (2^(intDsize/2) <= b' <= a' < 2^intDsize).
194 // Rechne den Euklid-Algorithmus mit Beifaktoren für ALLE Zahlen (a,b) aus,
195 // die mit a' bzw. b' anfangen; das liefert x1,y1,x2,y2, so daß
196 // ggT(a,b) = ggT(x1*a-y1*b,-x2*a+y2*b) und x1*a-y1*b>=0,-x2*a+y2*b>=0.
197 // Genauer: Mit offensichtlicher Skalierung betrachten wir
198 // a als beliebiges Element des Intervalls [a',a'+1) und
199 // b als beliebiges Element des Intervalls [b',b'+1) und
200 // führen den Euklid-Algorithmus schrittweise durch:
201 // (x1,y1,z1) := (1,0,a'), (x2,y2,z2) := (0,1,b'),
203 // {Hier x1*a'-y1*b'=z1, x1*a-y1*b in [z1-y1,z1+x1), z1-y1>=0, z1>0,
204 // und -x2*a'+y2*b'=z2, -x2*a+y2*b in [z2-x2,z2+y2), z2-x2>=0, z2>0,
205 // x1*y2-x2*y1=1, x1*z2+x2*z1=b', y1*z2+y2*z1=a'.}
206 // Falls z1-y1>=z2+y2:
207 // (x1,y1,z1) := (x1+x2,y1+y2,z1-z2), goto Schleife.
208 // Falls z2-x2>=z1+x1:
209 // (x2,y2,z2) := (x2+x1,y2+y1,z2-z1), goto Schleife.
210 // Sonst muß man abbrechen.
211 // {Zu den Schleifeninvarianten:
212 // 1. Die Gleichungen x1*a'-y1*b'=z1, -x2*a'+y2*b'=z2,
213 // x1*y2-x2*y1=1, x1*z2+x2*z1=b', y1*z2+y2*z1=a' mit Induktion.
214 // 2. Die Ungleichungen x1>0, y1>=0, x2>=0, y2>0 mit Induktion.
215 // 3. Die Ungleichungen z1>=0, z2>=0 nach Fallauswahl.
216 // 4. Die Ungleichung x1+x2>0 aus x1*z2+x2*z1=b'>0,
217 // die Ungleichung y1+y2>0 aus y1*z2+y2*z1=a'>0.
218 // 5. Die Ungleichung z1>0 wegen Fallauswahl und y1+y2>0,
219 // Die Ungleichung z2>0 wegen Fallauswahl und x1+x2>0.
220 // 6. Die Ungleichungen z1-y1>=0, z2-x2>=0 wegen Fallauswahl.
221 // 7. Die Ungleichung max(z1,z2) <= a' mit Induktion.
222 // 8. Die Ungleichung x1+x2 <= x1*z2+x2*z1 = b',
223 // die Ungleichung y1+y2 <= y1*z2+y2*z1 = a'.
224 // Damit bleiben alle Größen im Intervall [0,beta), kein Überlauf.
225 // 9. Die Ungleichungen z1+x1<=beta, z2+y2<=beta mit Induktion.
226 // 10. x1*a-y1*b in (z1-y1,z1+x1) (bzw. [z1,z1+x1) bei y1=0),
227 // -x2*a+y2*b in (z2-x2,z2+y2) (bzw. [z2,z2+y2) bei x2=0),
228 // da a in a'+[0,1) und b in b'+[0,1).
229 // Jedenfalls 0 < x1*a-y1*b < z1+x1 <= x2*z1+x1*z2 = b' falls x2>0,
230 // und 0 < -x2*a+y2*b < z2+y2 <= y1*z2+y2*z1 = a' falls y1>0.}
231 // Man kann natürlich auch mehrere Subtraktionsschritte auf einmal
233 // Falls q := floor((z1-y1)/(z2+y2)) > 0 :
234 // (x1,y1,z1) := (x1+q*x2,y1+q*y2,z1-q*z2), goto Schleife.
235 // Falls q := floor((z2-x2)/(z1+x1)) > 0 :
236 // (x2,y2,z2) := (x2+q*x1,y2+q*y1,z2-q*z1), goto Schleife.
237 // {Am Schluß gilt -(x1+x2) < z1-z2 < y1+y2 und daher
238 // z2-x2 <= b'/(x1+x2) < z1+x1, z1-y1 <= a'/(y1+y2) < z2+y2,
239 // und - unter Berücksichtigung von x1*y2-x2*y1=1 -
240 // z1-y1 <= b'/(x1+x2) < z2+y2, z2-x2 <= a'/(y1+y2) < z1+x1,
241 // also max(z1-y1,z2-x2) <= min(b'/(x1+x2),a'/(y1+y2))
242 // <= max(b'/(x1+x2),a'/(y1+y2)) < min(z1+x1,z2+y2).}
243 // Im Fall y1=x2=0 => x1=y2=1 (der nur bei a'=z1=z2=b' eintreten kann)
244 // ersetze (a,b) := (a-b,b). {Beide >0, da a>b>0 war.}
245 // Der Fall y1=0,x2>0 => x1=y2=1 => a' = z1 < z2+x2*z1 = b'
246 // kann nicht eintreten.
247 // Im Fall x2=0,y1>0 => x1=y2=1 ersetze (a,b) := (a-y1*b,b).
248 // {Das ist OK, da 0 <= z1-y1 = a'-y1*(b'+1) < a-y1*b < a.}
249 // Sonst (y1>0,x2>0) ersetze (a,b) := (x1*a-y1*b,-x2*a+y2*b).
250 // {Das ist OK, da 0 <= z1-y1 = x1*a'-y1*(b'+1) < x1*a-y1*b
251 // und 0 <= z2-x2 = -x2*(a'+1)+y2*b' < -x2*a+y2*b
252 // und x1*a-y1*b < x1*(a'+1)-y1*b' = z1+x1 <= x2*z1+x1*z2 = b' <= b
253 // und -x2*a+y2*b < -x2*a'+y2*(b'+1) = z2+y2 <= y1*z2+y2*z1 = a' <= a.}
256 // Define this to 1 in order to use double-word sized a' and b'.
257 // This gives better x1,y1,x2,y2, because normally the values x1,y1,x2,y2
258 // have only about intDsize/2 bits and so half of the multiplication work
259 // is lost. Actually, this flag multiplies the gcd speed by 1.5, not 2.0.
260 #define DOUBLE_SPEED 1
261 // Speed increases only for large integers. For small ones, computing
262 // with double-word sized a' and b' is too costly. The threshold is
263 // between 12 and 20, around 15.
264 #define cl_gcd_double_threshold 16
266 // gcd of two single-word numbers >0
267 static uintD gcdD (uintD a, uintD b)
268 { var uintD bit_j = (a | b); // endet mit einer 1 und j Nullen
269 bit_j = bit_j ^ (bit_j - 1); // Maske = bit(j) | bit(j-1) | ... | bit(0)
270 if (!((a & bit_j) ==0))
271 { if (!((b & bit_j) ==0)) goto odd_odd; else goto odd_even; }
272 if (!((b & bit_j) ==0)) goto even_odd;
275 { odd_odd: // a,b >0, beide ungerade
276 // Vergleiche a und b:
277 if (a == b) break; // a=b>0 -> fertig
280 even_odd: // a,b >0, a gerade, b ungerade
281 do { a = a>>1; } while ((a & bit_j) ==0);
285 odd_even: // a,b >0, a ungerade, b gerade
286 do { b = b>>1; } while ((b & bit_j) ==0);
293 const cl_I gcd (const cl_I& a, const cl_I& b)
294 { if (eq(a,1)) { return 1; } // a=1 -> 1
295 if (eq(b,1)) { return 1; } // b=1 -> 1
296 if (eq(b,0)) { return abs(a); } // b=0 -> (abs a)
297 if (eq(a,0)) { return abs(b); } // a=0 -> (abs b)
298 if (fixnump(a) && fixnump(b)) // ggT zweier Fixnums /=0
299 { var sintV a_ = FN_to_V(a);
300 if (a_ < 0) { a_ = -a_; }
301 var sintV b_ = FN_to_V(b);
302 if (b_ < 0) { b_ = -b_; }
303 return UV_to_I(gcd((uintV)a_,(uintV)b_));
312 // Macro: erzeugt die NUDS zu (abs x), erniedrigt num_stack
313 #define I_abs_to_NUDS(x) \
314 I_to_NDS_1(x, x##_MSDptr = , x##_len = , x##_LSDptr = ); /* (nichtleere) NDS holen */\
315 if ((sintD)mspref(x##_MSDptr,0) < 0) /* falls <0, negieren: */\
316 { neg_loop_lsp(x##_LSDptr,x##_len); } \
317 if (mspref(x##_MSDptr,0) == 0) /* normalisieren (max. 1 Nulldigit entfernen) */\
318 { msshrink(x##_MSDptr); x##_len--; }
319 I_abs_to_NUDS(a); // (abs a) als NUDS erzeugen
320 I_abs_to_NUDS(b); // (abs b) als NUDS erzeugen
321 // Jetzt ist a = a_MSDptr/a_len/a_LSDptr, b = b_MSDptr/b_len/b_LSDptr,
322 // beides NUDS, und a_len>0, b_len>0.
323 // Platz für zwei Rechenregister besorgen, mit je max(a_len,b_len)+1 Digits:
324 {var uintD* divroomptr; // Platz für Divisionsergebnis
327 {var uintC c_len = (a_len>=b_len ? a_len : b_len) + 1;
328 num_stack_alloc(c_len,divroomptr=,c_LSDptr=);
329 num_stack_alloc(c_len,,d_LSDptr=);
330 // Jetzt ist ../c_len/c_LSDptr, ../c_len/d_LSDptr frei.
333 { // Hier a,b>0, beides NUDS.
334 // Vergleiche a und b:
335 if (a_len > b_len) goto a_greater_b; // a>b ?
337 { var cl_signean vergleich = compare_loop_msp(a_MSDptr,b_MSDptr,a_len);
338 if (vergleich > 0) goto a_greater_b; // a>b ?
339 if (vergleich == 0) break; // a=b ?
341 // a<b -> a,b vertauschen:
342 swap(uintD*, a_MSDptr,b_MSDptr);
343 swap(uintC, a_len,b_len);
344 swap(uintD*, a_LSDptr,b_LSDptr);
346 // Hier a>b>0, beides NUDS.
347 if (b_len==1) // Beschleunigung eines häufigen Falles
348 { var uintD b0 = mspref(b_MSDptr,0);
350 // a>b=1 -> Ergebnis 1.
352 // a>b>1 -> evtl. Division durch b
355 { a0 = mspref(a_MSDptr,0); }
357 { a0 = divu_loop_msp(b0,a_MSDptr,a_len);
359 { return UD_to_I(b0); }
361 return UD_to_I(gcdD(a0,b0));
363 // Entscheidung, ob Division oder Linearkombination:
364 { var uintD a_msd; // führende intDsize Bits von a
365 var uintD b_msd; // entsprechende Bits von b
367 var uintD a_nsd; // nächste intDsize Bits von a
368 var uintD b_nsd; // entsprechende Bits von b
370 { var uintC len_diff = a_len-b_len; // Längendifferenz
371 if (len_diff > 1) goto divide; // >=2 -> Bitlängendifferenz>intDsize -> dividieren
372 #define bitlendiff_limit (intDsize/2) // sollte >0,<intDsize sein
373 {var uintC a_msd_size;
374 a_msd = mspref(a_MSDptr,0); // führendes Digit von a
375 integerlengthD(a_msd,a_msd_size=); // dessen Bit-Länge (>0,<=intDsize) berechnen
376 b_msd = mspref(b_MSDptr,0);
378 {var uintDD b_msdd = // 2 führende Digits von b
380 ? highlowDD(b_msd, mspref(b_MSDptr,1))
383 // a_msd_size+intDsize - b_msdd_size >= bitlendiff_limit -> dividieren:
384 b_msd = lowD(b_msdd >> a_msd_size);
385 if (b_msd < (uintD)bit(intDsize-bitlendiff_limit)) goto divide;
387 b_nsd = lowD(highlowDD(lowD(b_msdd), (b_len<=2-len_diff ? 0 : mspref(b_MSDptr,2-len_diff))) >> a_msd_size);
390 {var uintDD a_msdd = // 2 führende Digits von a
391 highlowDD(a_msd, mspref(a_MSDptr,1));
392 a_msd = lowD(a_msdd >> a_msd_size);
394 a_nsd = lowD(highlowDD(lowD(a_msdd), (a_len<=2 ? 0 : mspref(a_MSDptr,2))) >> a_msd_size);
397 if (a_msd == b_msd) goto subtract;
400 { // a_msd_size - b_msd_size >= bitlendiff_limit -> dividieren:
401 if ((a_msd_size > bitlendiff_limit)
402 && (b_msd < (uintD)bit(a_msd_size-bitlendiff_limit))
405 // Entscheidung für Linearkombination ist gefallen.
406 // a_msd und b_msd so erweitern, daß a_msd die führenden
407 // intDsize Bits von a enthält:
408 {var uintC shiftcount = intDsize-a_msd_size; // Shiftcount nach links (>=0, <intDsize)
410 { a_msd = a_msd << shiftcount;
411 b_msd = b_msd << shiftcount;
412 a_msd |= mspref(a_MSDptr,1) >> a_msd_size;
413 b_msd |= mspref(b_MSDptr,1) >> a_msd_size;
415 if (a_msd == b_msd) goto subtract;
417 a_nsd = mspref(a_MSDptr,1);
418 b_nsd = mspref(b_MSDptr,1);
420 { a_nsd = a_nsd << shiftcount;
421 b_nsd = b_nsd << shiftcount;
423 { a_nsd |= mspref(a_MSDptr,2) >> a_msd_size;
424 b_nsd |= mspref(b_MSDptr,2) >> a_msd_size;
430 { // a_msd_size+intDsize - b_msd_size >= bitlendiff_limit -> dividieren:
431 if ((a_msd_size >= bitlendiff_limit)
432 || (b_msd < (uintD)bit(a_msd_size+intDsize-bitlendiff_limit))
435 // Entscheidung für Linearkombination ist gefallen.
436 // a_msd und b_msd so erweitern, daß a_msd die führenden
437 // intDsize Bits von a enthält:
438 // 0 < a_msd_size < b_msd_size + bitlendiff_limit - intDsize <= bitlendiff_limit < intDsize.
439 a_msd = (a_msd << (intDsize-a_msd_size)) | (mspref(a_MSDptr,1) >> a_msd_size);
441 a_nsd = mspref(a_MSDptr,1) << (intDsize-a_msd_size);
442 b_nsd = b_msd << (intDsize-a_msd_size);
443 a_nsd |= mspref(a_MSDptr,2) >> a_msd_size;
444 b_nsd |= mspref(b_MSDptr,1) >> a_msd_size;
446 b_msd = b_msd >> a_msd_size;
449 #undef bitlendiff_limit
451 // Nun ist a_msd = a' > b' = b_msd.
452 { // Euklid-Algorithmus auf den führenden Digits durchführen:
453 var partial_gcd_result likobi;
455 if (a_len >= cl_gcd_double_threshold)
458 partial_gcd(highlowDD(a_msd,a_nsd),highlowDD(b_msd,b_nsd),&likobi); // liefert x1,y1,x2,y2
460 partial_gcd(a_msd,a_nsd,b_msd,b_nsd,&likobi); // liefert x1,y1,x2,y2
465 { partial_gcd(a_msd,b_msd,&likobi); } // liefert x1,y1,x2,y2, aber nur halb so gut
468 { // Ersetze (a,b) := (a-y1*b,b).
469 if (likobi.y1==1) goto subtract; // einfacherer Fall
470 // Dazu evtl. a um 1 Digit erweitern, so daß a_len=b_len+1:
471 if (a_len == b_len) { lsprefnext(a_MSDptr) = 0; a_len++; }
472 // und y1*b von a subtrahieren:
473 mspref(a_MSDptr,0) -= mulusub_loop_lsp(likobi.y1,b_LSDptr,a_LSDptr,b_len);
476 { // Ersetze (a,b) := (x1*a-y1*b,-x2*a+y2*b).
477 // Dazu evtl. b um 1 Digit erweitern, so daß a_len=b_len:
478 if (!(a_len==b_len)) { lsprefnext(b_MSDptr) = 0; b_len++; }
479 // c := x1*a-y1*b bilden:
480 mulu_loop_lsp(likobi.x1,a_LSDptr,c_LSDptr,a_len);
481 /* lspref(c_LSDptr,a_len) -= */
482 mulusub_loop_lsp(likobi.y1,b_LSDptr,c_LSDptr,a_len);
483 // d := -x2*a+y2*b bilden:
484 mulu_loop_lsp(likobi.y2,b_LSDptr,d_LSDptr,a_len);
485 /* lspref(d_LSDptr,a_len) -= */
486 mulusub_loop_lsp(likobi.x2,a_LSDptr,d_LSDptr,a_len);
487 // Wir wissen, daß 0 < c < b und 0 < d < a. Daher müßten
488 // lspref(c_LSDptr,a_len) und lspref(d_LSDptr,a_len) =0 sein.
489 // a := c und b := d kopieren:
490 copy_loop_lsp(c_LSDptr,a_LSDptr,a_len);
491 copy_loop_lsp(d_LSDptr,b_LSDptr,a_len);
493 while (mspref(b_MSDptr,0)==0) { msshrink(b_MSDptr); b_len--; }
496 { subtract: // Ersetze (a,b) := (a-b,b).
497 if (!( subfrom_loop_lsp(b_LSDptr,a_LSDptr,b_len) ==0))
498 // Übertrag nach b_len Stellen, muß also a_len=b_len+1 sein.
499 { mspref(a_MSDptr,0) -= 1; }
502 while (mspref(a_MSDptr,0)==0) { msshrink(a_MSDptr); a_len--; }
505 { divide: // Ersetze (a,b) := (b , a mod b).
506 {var uintD* old_a_LSDptr = a_LSDptr;
509 cl_UDS_divide(a_MSDptr,a_len,a_LSDptr,b_MSDptr,b_len,b_LSDptr, divroomptr, &q,&r);
510 a_MSDptr = b_MSDptr; a_len = b_len; a_LSDptr = b_LSDptr; // a := b
511 b_len = r.len; if (b_len==0) break; // b=0 -> fertig
512 b_LSDptr = old_a_LSDptr; // b übernimmt den vorherigen Platz von a
513 b_MSDptr = copy_loop_lsp(r.LSDptr,b_LSDptr,b_len); // b := r kopieren
514 goto a_greater_b; // Nun ist a>b>0
517 return NUDS_to_I(a_MSDptr,a_len); // NUDS a als Ergebnis
521 #endif /* GCD_ALGO == 3 */