// gcd().
// General includes.
#include "cl_sysdep.h"
// Specification.
#include "cln/integer.h"
// Implementation.
#include "cl_I.h"
#include "cl_DS.h"
#include "cl_D.h"
#include "cl_xmacros.h"
namespace cln {
#define GCD_ALGO 3 // 1: binär, 2: Schulmethode, 3: Lehmer
#if (GCD_ALGO == 1)
// binäre Methode:
// (gcd a b) :==
// b=0 --> (abs a)
// a=0 --> (abs b)
// sonst:
// (abs a) und (abs b) in zwei Buffer packen, als Unsigned Digit Sequences.
// [Schreibe oBdA wieder a,b]
// (prog ((j 0))
// 1 {a,b >0}
// (if (evenp a)
// (if (evenp b)
// (progn (incf j) (setq a (/ a 2)) (setq b (/ b 2)) (go 1))
// (go 4)
// )
// (while (evenp b) (setq b (/ b 2)))
// )
// 2 {a,b >0, beide ungerade}
// (cond ((> a b))
// ((= a b) (go 5))
// ((< a b) (rotatef a b))
// )
// 3 {a,b >0, beide ungerade, a>b}
// (setq a (- a b))
// 4 {a,b >0, a gerade, b ungerade}
// (repeat (setq a (/ a 2)) (until (oddp a)))
// (go 2)
// 5 {a=b>0}
// (return (ash a j))
// )
// weil es oft auftritt (insbesondere bei GCD's mehrerer Zahlen):
// a=1 oder b=1 --> 1
const cl_I gcd (const cl_I& a, const cl_I& b)
{ if (eq(a,1)) { return 1; } // a=1 -> 1
if (eq(b,1)) { return 1; } // b=1 -> 1
if (eq(b,0)) { return abs(a); } // b=0 -> (abs a)
if (eq(a,0)) { return abs(b); } // a=0 -> (abs b)
CL_ALLOCA_STACK;
var uintD* a_MSDptr;
var uintC a_len;
var uintD* a_LSDptr;
var uintD* b_MSDptr;
var uintC b_len;
var uintD* b_LSDptr;
// Macro: erzeugt die NUDS zu (abs x), erniedrigt num_stack
#define I_abs_to_NUDS(x) \
I_to_NDS_1(x, x##_MSDptr = , x##_len = , x##_LSDptr = ); /* (nichtleere) NDS holen */\
if ((sintD)mspref(x##_MSDptr,0) < 0) /* falls <0, negieren: */\
{ neg_loop_lsp(x##_LSDptr,x##_len); } \
if (mspref(x##_MSDptr,0) == 0) /* normalisieren (max. 1 Nulldigit entfernen) */\
{ msshrink(x##_MSDptr); x##_len--; }
I_abs_to_NUDS(a); // (abs a) als NUDS erzeugen
I_abs_to_NUDS(b); // (abs b) als NUDS erzeugen
// Jetzt ist a = a_MSDptr/a_len/a_LSDptr, b = b_MSDptr/b_len/b_LSDptr,
// beides NUDS, und a_len>0, b_len>0.
// Macro: Halbiere x.
#define halb(x) \
{ shift1right_loop_msp(x##_MSDptr,x##_len,0); /* um 1 Bit rechts schieben */ \
if (mspref(x##_MSDptr,0) == 0) { msshrink(x##_MSDptr); x##_len--; } /* normalisieren */\
}
// Macro: Ob x gerade ist.
#define evenp(x) \
((lspref(x##_LSDptr,0) & bit(0)) ==0)
{ var uintL j = 0;
label_1: // a,b >0
if (evenp(a))
{ if (evenp(b))
{ j++; halb(a); halb(b); goto label_1; }
else
goto label_4;
}
while (evenp(b)) { halb(b); }
label_2: // a,b >0, beide ungerade
// Vergleiche a und b:
if (a_len > b_len) goto label_3; // a>b ?
if (a_len == b_len)
{ var cl_signean vergleich = compare_loop_msp(a_MSDptr,b_MSDptr,a_len);
if (vergleich > 0) goto label_3; // a>b ?
if (vergleich == 0) goto label_5; // a=b ?
}
// a<b -> a,b vertauschen:
swap(uintD*, a_MSDptr,b_MSDptr);
swap(uintC, a_len,b_len);
swap(uintD*, a_LSDptr,b_LSDptr);
label_3: // a,b >0, beide ungerade, a>b
// subtrahiere a := a - b
if (!( subfrom_loop_lsp(b_LSDptr,a_LSDptr,b_len) ==0))
{ dec_loop_lsp(a_LSDptr lspop b_len,a_len-b_len); }
while (mspref(a_MSDptr,0) == 0) { msshrink(a_MSDptr); a_len--; } // normalisieren
label_4: // a,b >0, a gerade, b ungerade
do { halb(a); } while (evenp(a));
goto label_2;
label_5: // a=b>0
// a zu einer NDS machen:
return ash(NUDS_to_I(a_MSDptr,a_len),j); // ggT der ungeraden Anteile als Integer, mal 2^j
}
#undef evenp
#undef halb
#undef I_abs_to_NUDS
}
#endif /* GCD_ALGO == 1 */
#if (GCD_ALGO == 2)
// Schulmethode:
// (gcd a b) :==
// [a:=(abs a), b:=(abs b), while b>0 do (a,b) := (b,(mod a b)), -> a]
// verbessert:
// a=0 -> (abs b)
// b=0 -> (abs a)
// a=1 -> 1
// b=1 -> 1
// a:=(abs a), b:=(abs b)
// Falls a=b: return a; falls a<b: vertausche a und b.
// (*) {Hier a>b>0}
// Falls b=1, return 1. {spart eine Division durch 1}
// Sonst dividieren (divide a b), a:=b, b:=Rest.
// Falls b=0, return a, sonst goto (*).
const cl_I gcd (const cl_I& a, const cl_I& b)
{ if (eq(a,1)) { return 1; } // a=1 -> 1
if (eq(b,1)) { return 1; } // b=1 -> 1
if (eq(b,0)) { return abs(a); } // b=0 -> (abs a)
if (eq(a,0)) { return abs(b); } // a=0 -> (abs b)
// Beträge nehmen:
{var cl_I abs_a = abs(a);
var cl_I abs_b = abs(b);
var cl_I& a = abs_a;
var cl_I& b = abs_b;
if (fixnump(a) && fixnump(b)) // ggT zweier Fixnums >0
{ // bleibt Fixnum, da (gcd a b) <= (min a b)
return L_to_FN(gcd(FN_to_UL(a),FN_to_UL(b)));
}
{ var cl_signean vergleich = compare(a,b);
if (vergleich == 0) { return a; } // a=b -> fertig
if (vergleich < 0) { var cl_I tmp = a; a = b; b = a; } // a<b -> a,b vertauschen
}
loop // Hier a > b > 0
{ if (eq(b,1)) { return 1; } // b=1 -> Ergebnis 1
{ var cl_I_div_t div = cl_divide(a,b); a = b; b = div.remainder; }
if (eq(b,0)) { return a; }
}
}}
#endif /* GCD_ALGO == 2 */
#if (GCD_ALGO == 3)
// Lehmer-Methode:
// vgl. [ D. E. Knuth: The Art of Computer Programming, Vol. 2: Seminumerical
// Algorithms, Sect. 4.5.2., Algorithm L ]
// und [ Collins, Loos: SAC-2, Algorithms IGCD, DPCC ].
// (gcd a b) :==
// a=0 -> (abs b)
// b=0 -> (abs a)
// a=1 -> 1
// b=1 -> 1
// a:=(abs a), b:=(abs b)
// (*) {Hier a,b>0}
// Falls a=b: return a; falls a<b: vertausche a und b.
// {Hier a>b>0}
// Falls (- (integer-length a) (integer-length b)) >= intDsize/2,
// lohnt sich eine Division: (a,b) := (b , a mod b). Falls b=0: return a.
// Falls dagegen 0 <= (- (integer-length a) (integer-length b)) < intDsize/2,
// seien a' die führenden intDsize Bits von a
// (2^(intDsize-1) <= a' < 2^intDsize) und b' die entsprechenden Bits von b
// (2^(intDsize/2) <= b' <= a' < 2^intDsize).
// Rechne den Euklid-Algorithmus mit Beifaktoren für ALLE Zahlen (a,b) aus,
// die mit a' bzw. b' anfangen; das liefert x1,y1,x2,y2, so daß
// ggT(a,b) = ggT(x1*a-y1*b,-x2*a+y2*b) und x1*a-y1*b>=0,-x2*a+y2*b>=0.
// Genauer: Mit offensichtlicher Skalierung betrachten wir
// a als beliebiges Element des Intervalls [a',a'+1) und
// b als beliebiges Element des Intervalls [b',b'+1) und
// führen den Euklid-Algorithmus schrittweise durch:
// (x1,y1,z1) := (1,0,a'), (x2,y2,z2) := (0,1,b'),
// Schleife:
// {Hier x1*a'-y1*b'=z1, x1*a-y1*b in [z1-y1,z1+x1), z1-y1>=0, z1>0,
// und -x2*a'+y2*b'=z2, -x2*a+y2*b in [z2-x2,z2+y2), z2-x2>=0, z2>0,
// x1*y2-x2*y1=1, x1*z2+x2*z1=b', y1*z2+y2*z1=a'.}
// Falls z1-y1>=z2+y2:
// (x1,y1,z1) := (x1+x2,y1+y2,z1-z2), goto Schleife.
// Falls z2-x2>=z1+x1:
// (x2,y2,z2) := (x2+x1,y2+y1,z2-z1), goto Schleife.
// Sonst muß man abbrechen.
// {Zu den Schleifeninvarianten:
// 1. Die Gleichungen x1*a'-y1*b'=z1, -x2*a'+y2*b'=z2,
// x1*y2-x2*y1=1, x1*z2+x2*z1=b', y1*z2+y2*z1=a' mit Induktion.
// 2. Die Ungleichungen x1>0, y1>=0, x2>=0, y2>0 mit Induktion.
// 3. Die Ungleichungen z1>=0, z2>=0 nach Fallauswahl.
// 4. Die Ungleichung x1+x2>0 aus x1*z2+x2*z1=b'>0,
// die Ungleichung y1+y2>0 aus y1*z2+y2*z1=a'>0.
// 5. Die Ungleichung z1>0 wegen Fallauswahl und y1+y2>0,
// Die Ungleichung z2>0 wegen Fallauswahl und x1+x2>0.
// 6. Die Ungleichungen z1-y1>=0, z2-x2>=0 wegen Fallauswahl.
// 7. Die Ungleichung max(z1,z2) <= a' mit Induktion.
// 8. Die Ungleichung x1+x2 <= x1*z2+x2*z1 = b',
// die Ungleichung y1+y2 <= y1*z2+y2*z1 = a'.
// Damit bleiben alle Größen im Intervall [0,beta), kein Überlauf.
// 9. Die Ungleichungen z1+x1<=beta, z2+y2<=beta mit Induktion.
// 10. x1*a-y1*b in (z1-y1,z1+x1) (bzw. [z1,z1+x1) bei y1=0),
// -x2*a+y2*b in (z2-x2,z2+y2) (bzw. [z2,z2+y2) bei x2=0),
// da a in a'+[0,1) und b in b'+[0,1).
// Jedenfalls 0 < x1*a-y1*b < z1+x1 <= x2*z1+x1*z2 = b' falls x2>0,
// und 0 < -x2*a+y2*b < z2+y2 <= y1*z2+y2*z1 = a' falls y1>0.}
// Man kann natürlich auch mehrere Subtraktionsschritte auf einmal
// durchführen:
// Falls q := floor((z1-y1)/(z2+y2)) > 0 :
// (x1,y1,z1) := (x1+q*x2,y1+q*y2,z1-q*z2), goto Schleife.
// Falls q := floor((z2-x2)/(z1+x1)) > 0 :
// (x2,y2,z2) := (x2+q*x1,y2+q*y1,z2-q*z1), goto Schleife.
// {Am Schluß gilt -(x1+x2) < z1-z2 < y1+y2 und daher
// z2-x2 <= b'/(x1+x2) < z1+x1, z1-y1 <= a'/(y1+y2) < z2+y2,
// und - unter Berücksichtigung von x1*y2-x2*y1=1 -
// z1-y1 <= b'/(x1+x2) < z2+y2, z2-x2 <= a'/(y1+y2) < z1+x1,
// also max(z1-y1,z2-x2) <= min(b'/(x1+x2),a'/(y1+y2))
// <= max(b'/(x1+x2),a'/(y1+y2)) < min(z1+x1,z2+y2).}
// Im Fall y1=x2=0 => x1=y2=1 (der nur bei a'=z1=z2=b' eintreten kann)
// ersetze (a,b) := (a-b,b). {Beide >0, da a>b>0 war.}
// Der Fall y1=0,x2>0 => x1=y2=1 => a' = z1 < z2+x2*z1 = b'
// kann nicht eintreten.
// Im Fall x2=0,y1>0 => x1=y2=1 ersetze (a,b) := (a-y1*b,b).
// {Das ist OK, da 0 <= z1-y1 = a'-y1*(b'+1) < a-y1*b < a.}
// Sonst (y1>0,x2>0) ersetze (a,b) := (x1*a-y1*b,-x2*a+y2*b).
// {Das ist OK, da 0 <= z1-y1 = x1*a'-y1*(b'+1) < x1*a-y1*b
// und 0 <= z2-x2 = -x2*(a'+1)+y2*b' < -x2*a+y2*b
// und x1*a-y1*b < x1*(a'+1)-y1*b' = z1+x1 <= x2*z1+x1*z2 = b' <= b
// und -x2*a+y2*b < -x2*a'+y2*(b'+1) = z2+y2 <= y1*z2+y2*z1 = a' <= a.}
// goto (*).
// Define this to 1 in order to use double-word sized a' and b'.
// This gives better x1,y1,x2,y2, because normally the values x1,y1,x2,y2
// have only about intDsize/2 bits and so half of the multiplication work
// is lost. Actually, this flag multiplies the gcd speed by 1.5, not 2.0.
#define DOUBLE_SPEED 1
// Speed increases only for large integers. For small ones, computing
// with double-word sized a' and b' is too costly. The threshold is
// between 12 and 20, around 15.
#define cl_gcd_double_threshold 16
// gcd of two single-word numbers >0
static uintD gcdD (uintD a, uintD b)
{ var uintD bit_j = (a | b); // endet mit einer 1 und j Nullen
bit_j = bit_j ^ (bit_j - 1); // Maske = bit(j) | bit(j-1) | ... | bit(0)
if (!((a & bit_j) ==0))
{ if (!((b & bit_j) ==0)) goto odd_odd; else goto odd_even; }
if (!((b & bit_j) ==0)) goto even_odd;
NOTREACHED;
loop
{ odd_odd: // a,b >0, beide ungerade
// Vergleiche a und b:
if (a == b) break; // a=b>0 -> fertig
if (a > b) // a>b ?
{ a = a-b;
even_odd: // a,b >0, a gerade, b ungerade
do { a = a>>1; } while ((a & bit_j) ==0);
}
else // a<b
{ b = b-a;
odd_even: // a,b >0, a ungerade, b gerade
do { b = b>>1; } while ((b & bit_j) ==0);
}
}
// a=b>0
return a;
}
const cl_I gcd (const cl_I& a, const cl_I& b)
{ if (eq(a,1)) { return 1; } // a=1 -> 1
if (eq(b,1)) { return 1; } // b=1 -> 1
if (eq(b,0)) { return abs(a); } // b=0 -> (abs a)
if (eq(a,0)) { return abs(b); } // a=0 -> (abs b)
if (fixnump(a) && fixnump(b)) // ggT zweier Fixnums /=0
{ var sintL a_ = FN_to_L(a);
if (a_ < 0) { a_ = -a_; }
var sintL b_ = FN_to_L(b);
if (b_ < 0) { b_ = -b_; }
return UL_to_I(gcd((uint32)a_,(uint32)b_));
}
CL_ALLOCA_STACK;
var uintD* a_MSDptr;
var uintC a_len;
var uintD* a_LSDptr;
var uintD* b_MSDptr;
var uintC b_len;
var uintD* b_LSDptr;
// Macro: erzeugt die NUDS zu (abs x), erniedrigt num_stack
#define I_abs_to_NUDS(x) \
I_to_NDS_1(x, x##_MSDptr = , x##_len = , x##_LSDptr = ); /* (nichtleere) NDS holen */\
if ((sintD)mspref(x##_MSDptr,0) < 0) /* falls <0, negieren: */\
{ neg_loop_lsp(x##_LSDptr,x##_len); } \
if (mspref(x##_MSDptr,0) == 0) /* normalisieren (max. 1 Nulldigit entfernen) */\
{ msshrink(x##_MSDptr); x##_len--; }
I_abs_to_NUDS(a); // (abs a) als NUDS erzeugen
I_abs_to_NUDS(b); // (abs b) als NUDS erzeugen
// Jetzt ist a = a_MSDptr/a_len/a_LSDptr, b = b_MSDptr/b_len/b_LSDptr,
// beides NUDS, und a_len>0, b_len>0.
// Platz für zwei Rechenregister besorgen, mit je max(a_len,b_len)+1 Digits:
{var uintD* divroomptr; // Platz für Divisionsergebnis
var uintD* c_LSDptr;
var uintD* d_LSDptr;
{var uintL c_len = (uintL)(a_len>=b_len ? a_len : b_len) + 1;
num_stack_alloc(c_len,divroomptr=,c_LSDptr=);
num_stack_alloc(c_len,,d_LSDptr=);
// Jetzt ist ../c_len/c_LSDptr, ../c_len/d_LSDptr frei.
}
loop
{ // Hier a,b>0, beides NUDS.
// Vergleiche a und b:
if (a_len > b_len) goto a_greater_b; // a>b ?
if (a_len == b_len)
{ var cl_signean vergleich = compare_loop_msp(a_MSDptr,b_MSDptr,a_len);
if (vergleich > 0) goto a_greater_b; // a>b ?
if (vergleich == 0) break; // a=b ?
}
// a<b -> a,b vertauschen:
swap(uintD*, a_MSDptr,b_MSDptr);
swap(uintC, a_len,b_len);
swap(uintD*, a_LSDptr,b_LSDptr);
a_greater_b:
// Hier a>b>0, beides NUDS.
if (b_len==1) // Beschleunigung eines häufigen Falles
{ var uintD b0 = mspref(b_MSDptr,0);
if (b0==1)
// a>b=1 -> Ergebnis 1.
{ return 1; }
// a>b>1 -> evtl. Division durch b
var uintD a0;
if (a_len==1)
{ a0 = mspref(a_MSDptr,0); }
else
{ a0 = divu_loop_msp(b0,a_MSDptr,a_len);
if (a0==0)
{ return UD_to_I(b0); }
}
return UD_to_I(gcdD(a0,b0));
}
// Entscheidung, ob Division oder Linearkombination:
{ var uintD a_msd; // führende intDsize Bits von a
var uintD b_msd; // entsprechende Bits von b
#if DOUBLE_SPEED
var uintD a_nsd; // nächste intDsize Bits von a
var uintD b_nsd; // entsprechende Bits von b
#endif
{ var uintC len_diff = a_len-b_len; // Längendifferenz
if (len_diff > 1) goto divide; // >=2 -> Bitlängendifferenz>intDsize -> dividieren
#define bitlendiff_limit (intDsize/2) // sollte >0,<intDsize sein
{var uintC a_msd_size;
a_msd = mspref(a_MSDptr,0); // führendes Digit von a
integerlengthD(a_msd,a_msd_size=); // dessen Bit-Länge (>0,<=intDsize) berechnen
b_msd = mspref(b_MSDptr,0);
#if HAVE_DD
{var uintDD b_msdd = // 2 führende Digits von b
(len_diff==0
? highlowDD(b_msd, mspref(b_MSDptr,1))
: (uintDD)b_msd
);
// a_msd_size+intDsize - b_msdd_size >= bitlendiff_limit -> dividieren:
b_msd = lowD(b_msdd >> a_msd_size);
if (b_msd < (uintD)bit(intDsize-bitlendiff_limit)) goto divide;
#if DOUBLE_SPEED
b_nsd = lowD(highlowDD(lowD(b_msdd), (b_len<=2-len_diff ? 0 : mspref(b_MSDptr,2-len_diff))) >> a_msd_size);
#endif
}
{var uintDD a_msdd = // 2 führende Digits von a
highlowDD(a_msd, mspref(a_MSDptr,1));
a_msd = lowD(a_msdd >> a_msd_size);
#if DOUBLE_SPEED
a_nsd = lowD(highlowDD(lowD(a_msdd), (a_len<=2 ? 0 : mspref(a_MSDptr,2))) >> a_msd_size);
#endif
}
if (a_msd == b_msd) goto subtract;
#else
if (len_diff==0)
{ // a_msd_size - b_msd_size >= bitlendiff_limit -> dividieren:
if ((a_msd_size > bitlendiff_limit)
&& (b_msd < (uintD)bit(a_msd_size-bitlendiff_limit))
)
goto divide;
// Entscheidung für Linearkombination ist gefallen.
// a_msd und b_msd so erweitern, daß a_msd die führenden
// intDsize Bits von a enthält:
{var uintC shiftcount = intDsize-a_msd_size; // Shiftcount nach links (>=0, <intDsize)
if (shiftcount>0)
{ a_msd = a_msd << shiftcount;
b_msd = b_msd << shiftcount;
a_msd |= mspref(a_MSDptr,1) >> a_msd_size;
b_msd |= mspref(b_MSDptr,1) >> a_msd_size;
}
if (a_msd == b_msd) goto subtract;
#if DOUBLE_SPEED
a_nsd = mspref(a_MSDptr,1);
b_nsd = mspref(b_MSDptr,1);
if (shiftcount>0)
{ a_nsd = a_nsd << shiftcount;
b_nsd = b_nsd << shiftcount;
if (a_len>2)
{ a_nsd |= mspref(a_MSDptr,2) >> a_msd_size;
b_nsd |= mspref(b_MSDptr,2) >> a_msd_size;
} }
#endif
}}
else
// len_diff=1
{ // a_msd_size+intDsize - b_msd_size >= bitlendiff_limit -> dividieren:
if ((a_msd_size >= bitlendiff_limit)
|| (b_msd < (uintD)bit(a_msd_size+intDsize-bitlendiff_limit))
)
goto divide;
// Entscheidung für Linearkombination ist gefallen.
// a_msd und b_msd so erweitern, daß a_msd die führenden
// intDsize Bits von a enthält:
// 0 < a_msd_size < b_msd_size + bitlendiff_limit - intDsize <= bitlendiff_limit < intDsize.
a_msd = (a_msd << (intDsize-a_msd_size)) | (mspref(a_MSDptr,1) >> a_msd_size);
#if DOUBLE_SPEED
a_nsd = mspref(a_MSDptr,1) << (intDsize-a_msd_size);
b_nsd = b_msd << (intDsize-a_msd_size);
a_nsd |= mspref(a_MSDptr,2) >> a_msd_size;
b_nsd |= mspref(b_MSDptr,1) >> a_msd_size;
#endif
b_msd = b_msd >> a_msd_size;
}
#endif
#undef bitlendiff_limit
}}
// Nun ist a_msd = a' > b' = b_msd.
{ // Euklid-Algorithmus auf den führenden Digits durchführen:
var partial_gcd_result likobi;
#if DOUBLE_SPEED
if (a_len >= cl_gcd_double_threshold)
{
#if HAVE_DD
partial_gcd(highlowDD(a_msd,a_nsd),highlowDD(b_msd,b_nsd),&likobi); // liefert x1,y1,x2,y2
#else
partial_gcd(a_msd,a_nsd,b_msd,b_nsd,&likobi); // liefert x1,y1,x2,y2
#endif
}
else
#endif
{ partial_gcd(a_msd,b_msd,&likobi); } // liefert x1,y1,x2,y2, aber nur halb so gut
// Hier y1>0.
if (likobi.x2==0)
{ // Ersetze (a,b) := (a-y1*b,b).
if (likobi.y1==1) goto subtract; // einfacherer Fall
// Dazu evtl. a um 1 Digit erweitern, so daß a_len=b_len+1:
if (a_len == b_len) { lsprefnext(a_MSDptr) = 0; a_len++; }
// und y1*b von a subtrahieren:
mspref(a_MSDptr,0) -= mulusub_loop_lsp(likobi.y1,b_LSDptr,a_LSDptr,b_len);
}
else
{ // Ersetze (a,b) := (x1*a-y1*b,-x2*a+y2*b).
// Dazu evtl. b um 1 Digit erweitern, so daß a_len=b_len:
if (!(a_len==b_len)) { lsprefnext(b_MSDptr) = 0; b_len++; }
// c := x1*a-y1*b bilden:
mulu_loop_lsp(likobi.x1,a_LSDptr,c_LSDptr,a_len);
/* lspref(c_LSDptr,a_len) -= */
mulusub_loop_lsp(likobi.y1,b_LSDptr,c_LSDptr,a_len);
// d := -x2*a+y2*b bilden:
mulu_loop_lsp(likobi.y2,b_LSDptr,d_LSDptr,a_len);
/* lspref(d_LSDptr,a_len) -= */
mulusub_loop_lsp(likobi.x2,a_LSDptr,d_LSDptr,a_len);
// Wir wissen, daß 0 < c < b und 0 < d < a. Daher müßten
// lspref(c_LSDptr,a_len) und lspref(d_LSDptr,a_len) =0 sein.
// a := c und b := d kopieren:
copy_loop_lsp(c_LSDptr,a_LSDptr,a_len);
copy_loop_lsp(d_LSDptr,b_LSDptr,a_len);
// b normalisieren:
while (mspref(b_MSDptr,0)==0) { msshrink(b_MSDptr); b_len--; }
} }
if (cl_false)
{ subtract: // Ersetze (a,b) := (a-b,b).
if (!( subfrom_loop_lsp(b_LSDptr,a_LSDptr,b_len) ==0))
// Übertrag nach b_len Stellen, muß also a_len=b_len+1 sein.
{ mspref(a_MSDptr,0) -= 1; }
}
// a normalisieren:
while (mspref(a_MSDptr,0)==0) { msshrink(a_MSDptr); a_len--; }
}
if (cl_false)
{ divide: // Ersetze (a,b) := (b , a mod b).
{var uintD* old_a_LSDptr = a_LSDptr;
var DS q;
var DS r;
cl_UDS_divide(a_MSDptr,a_len,a_LSDptr,b_MSDptr,b_len,b_LSDptr, divroomptr, &q,&r);
a_MSDptr = b_MSDptr; a_len = b_len; a_LSDptr = b_LSDptr; // a := b
b_len = r.len; if (b_len==0) break; // b=0 -> fertig
b_LSDptr = old_a_LSDptr; // b übernimmt den vorherigen Platz von a
b_MSDptr = copy_loop_lsp(r.LSDptr,b_LSDptr,b_len); // b := r kopieren
goto a_greater_b; // Nun ist a>b>0
}}
} }
return NUDS_to_I(a_MSDptr,a_len); // NUDS a als Ergebnis
#undef I_abs_to_NUDS
}
#endif /* GCD_ALGO == 3 */
} // namespace cln
syntax highlighted by Code2HTML, v. 0.9.1