7 #include "cl_integer.h"
15 #include "cl_xmacros.h"
17 #define GCD_ALGO 3 // 1: binär, 2: Schulmethode, 3: Lehmer
24 // [a:=(abs A), b:=(abs B), while b>0 do (a,b) := (b,(mod a b)), -> a]
26 // A=1 -> return g=1, (u,v)=(1,0)
27 // B=1 -> return g=1, (u,v)=(0,1)
28 // a:=(abs A), ua:=(signum A), va:=0
29 // b:=(abs B), ub:=0, vb:=(signum B)
30 // A=0 -> return g=b, (u,v) = (ub,vb)
31 // B=0 -> return g=a, (u,v) = (ua,va)
32 // {Stets ua*A+va*B=a, ub*A+vb*B=b, ua*vb-ub*va = +/- 1.}
33 // Falls a=b: return a,ua,va;
34 // falls a<b: vertausche a und b, ua und ub, va und vb.
36 // Falls b=1, return 1,ub,vb. {spart eine Division durch 1}
37 // Sonst dividieren (divide a b) -> q,r.
38 // Falls r=0, return b,ub,vb.
39 // a:=b, b := Rest r = a-q*b, (ua,va,ub,vb) := (ub,vb,ua-q*ub,va-q*vb).
42 const cl_I xgcd (const cl_I& a, const cl_I& b, cl_I* u, cl_I* v)
43 { if (eq(a,1)) // a=1 -> g=1, (u,v)=(1,0)
44 { *u = 1; *v = 0; return a; }
45 if (eq(b,1)) // b=1 -> g=1, (u,v)=(0,1)
46 { *u = 0; *v = 1; return b; }
48 var cl_I ua = (minusp(a) ? cl_I(-1) : cl_I(1)); // ua := +/- 1
51 var cl_I vb = (minusp(b) ? cl_I(-1) : cl_I(1)); // vb := +/- 1
53 {var cl_I abs_a = abs(a);
54 var cl_I abs_b = abs(b);
57 if (eq(b,0)) // b=0 -> g=a, (u,v) = (ua,va)
58 { *u = ua; *v = va; return a; }
59 if (eq(a,0)) // a=0 -> g=b, (u,v) = (ub,vb)
60 { *u = ub; *v = vb; return b; }
61 { var cl_signean vergleich = cl_compare(a,b);
62 if (vergleich == 0) // a=b -> fertig
63 { *u = ua; *v = va; return a; }
64 if (vergleich < 0) // a<b -> a,b vertauschen
65 { swap(cl_I,a,b); swap(cl_I,ua,ub); swap(cl_I,va,vb); }
68 { if (eq(b,1)) // b=1 -> g=b, (u,v) = (ub,vb)
69 { *u = ub; *v = vb; return b; }
70 var cl_I_div_t div = cl_divide(a,b); // Division a / b
71 var cl_I& q = div.quotient;
72 var cl_I& r = div.remainder;
73 if (eq(r,0)) // r=0 -> fertig
74 { *u = ub; *v = vb; return b; }
75 { var cl_I x = ua-q*ub; ua = ub; ub = x; }
76 { var cl_I x = va-q*vb; va = vb; vb = x; }
81 #endif /* GCD_ALGO == 2 */
86 // wie oben bei (gcd A B).
87 // Zusätzlich werden Variablen sA,sB,sk,uAa,uBa,uAb,uBb geführt,
88 // wobei sA,sB,sk Vorzeichen (+/- 1) und uAa,uBa,uAb,uBb Integers >=0 sind mit
89 // uAa * sA*A - uBa * sB*B = a,
90 // - uAb * sA*A + uBb * sB*B = b,
91 // ferner uAa * uBb - uAb * uBa = sk und daher (Cramersche Regel)
92 // uBb * a + uBa * b = sk*sA*A, uAb * a + uAa * b = sk*sB*B.
93 // Zu Beginn (a,b) := (|A|,|B|), (sA,sB) := ((signum A), (signumB)),
94 // (uAa,uBa,uAb,uBb) := (1,0,0,1).
95 // Beim Ersetzen (a,b) := (a-b,b)
96 // ersetzt man (uAa,uBa,uAb,uBb) := (uAa+uAb,uBa+uBb,uAb,uBb).
97 // Beim Ersetzen (a,b) := (a-y1*b,b)
98 // ersetzt man (uAa,uBa,uAb,uBb) := (uAa+y1*uAb,uBa+y1*uBb,uAb,uBb).
99 // Beim Ersetzen (a,b) := (x1*a-y1*b,-x2*a+y2*b) mit x1*y2-x2*y1=1
100 // ersetzt man (uAa,uBa,uAb,uBb) :=
101 // (x1*uAa+y1*uAb,x1*uBa+y1*uBb,x2*uAa+y2*uAb,x2*uBa+y2*uBb).
102 // Beim Ersetzen (a,b) := (b,a)
103 // ersetzt man (uAa,uBa,uAb,uBb) := (uAb,uBb,uAa,uBa),
104 // sk := -sk, (sA,sB) := (-sA,-sB).
105 // Beim Ersetzen (a,b) := (b,a-q*b)
106 // ersetzt man (uAa,uBa,uAb,uBb) := (uAb,uBb,uAa+q*uAb,uBa+q*uBb),
107 // sk := -sk, (sA,sB) := (-sA,-sB).
108 // Zum Schluß ist a der ggT und a = uAa*sA * A + -uBa*sB * B
109 // die gewünschte Linearkombination.
110 // Da stets gilt sk*sA*A = |A|, sk*sB*B = |B|, a>=1, b>=1,
111 // folgt 0 <= uAa <= |B|, 0 <= uAb <= |B|, 0 <= uBa <= |A|, 0 <= uBb <= |A|.
112 // Ferner wird sk nie benutzt, braucht also nicht mitgeführt zu werden.
114 // Define this to 1 in order to use double-word sized a' and b'.
115 // This gives better x1,y1,x2,y2, because normally the values x1,y1,x2,y2
116 // have only about intDsize/2 bits and so half of the multiplication work
117 // is lost. Actually, this flag multiplies the gcd speed by 1.5, not 2.0.
118 #define DOUBLE_SPEED 1
120 // Bildet u := u + v, wobei für u genügend Platz sei:
121 // (Benutzt v.MSDptr nicht.)
122 static void NUDS_likobi0_NUDS (DS* u, DS* v)
123 { var uintC u_len = u->len;
124 var uintC v_len = v->len;
126 { if (!( addto_loop_lsp(v->LSDptr,u->LSDptr,v_len) ==0))
127 { if (!( inc_loop_lsp(u->LSDptr lspop v_len,u_len-v_len) ==0))
128 { lsprefnext(u->MSDptr) = 1; u->len++; }
130 else // u_len <= v_len
131 { u->MSDptr = copy_loop_lsp(v->LSDptr lspop u_len,u->LSDptr lspop u_len,v_len-u_len);
133 if (!( addto_loop_lsp(v->LSDptr,u->LSDptr,u_len) ==0))
134 { if (!( inc_loop_lsp(u->LSDptr lspop u_len,v_len-u_len) ==0))
135 { lsprefnext(u->MSDptr) = 1; u->len++; }
139 // Bildet u := u + q*v, wobei für u genügend Platz sei:
140 // (Dabei sei nachher u>0.)
141 static void NUDS_likobi1_NUDS (DS* u, DS* v, uintD q)
142 { var uintC v_len = v->len;
143 if (v_len>0) // nur nötig, falls v /=0
144 { var uintC u_len = u->len;
146 if (u_len <= v_len) // evtl. u vergrößern
147 { u->MSDptr = clear_loop_lsp(u->MSDptr,v_len-u_len+1);
148 u->len = u_len = v_len+1;
149 } // Nun ist u_len > v_len.
150 carry = muluadd_loop_lsp(q,v->LSDptr,u->LSDptr,v_len);
152 { var uintD* ptr = u->LSDptr lspop v_len;
153 if ((lspref(ptr,0) += carry) < carry)
154 { if (!( inc_loop_lsp(ptr lspop 1,u_len-v_len-1) ==0))
155 { lsprefnext(u->MSDptr) = 1; u->len++; }
157 while (mspref(u->MSDptr,0)==0) { msshrink(u->MSDptr); u->len--; } // normalisieren
160 // Bildet (u,v) := (x1*u+y1*v,x2*u+y2*v), wobei für u,v genügend Platz sei:
161 // (Dabei sei u>0 oder v>0, nachher u>0 und v>0.)
162 static void NUDS_likobi2_NUDS (DS* u, DS* v, partial_gcd_result* q, uintD* c_LSDptr, uintD* d_LSDptr)
163 { var uintC u_len = u->len;
164 var uintC v_len = v->len;
168 { mulu_loop_lsp(q->x1,u->LSDptr,c_LSDptr,u_len); c_len = u_len+1;
169 mulu_loop_lsp(q->x2,u->LSDptr,d_LSDptr,u_len); d_len = u_len+1;
172 muluadd_loop_lsp(q->y1,v->LSDptr,c_LSDptr,v_len);
174 { var uintD* ptr = c_LSDptr lspop v_len;
175 if ((lspref(ptr,0) += carry) < carry)
176 { if (!( inc_loop_lsp(ptr lspop 1,u_len-v_len) ==0))
177 { lspref(c_LSDptr,c_len) = 1; c_len++; }
180 muluadd_loop_lsp(q->y2,v->LSDptr,d_LSDptr,v_len);
182 { var uintD* ptr = d_LSDptr lspop v_len;
183 if ((lspref(ptr,0) += carry) < carry)
184 { if (!(inc_loop_lsp(ptr lspop 1,u_len-v_len) ==0))
185 { lspref(d_LSDptr,d_len) = 1; d_len++; }
189 { mulu_loop_lsp(q->y1,v->LSDptr,c_LSDptr,v_len); c_len = v_len+1;
190 mulu_loop_lsp(q->y2,v->LSDptr,d_LSDptr,v_len); d_len = v_len+1;
193 muluadd_loop_lsp(q->x1,u->LSDptr,c_LSDptr,u_len);
195 { var uintD* ptr = c_LSDptr lspop u_len;
196 if ((lspref(ptr,0) += carry) < carry)
197 { if (!( inc_loop_lsp(ptr lspop 1,v_len-u_len) ==0))
198 { lspref(c_LSDptr,c_len) = 1; c_len++; }
201 muluadd_loop_lsp(q->x2,u->LSDptr,d_LSDptr,u_len);
203 { var uintD* ptr = d_LSDptr lspop u_len;
204 if ((lspref(ptr,0) += carry) < carry)
205 { if (!( inc_loop_lsp(ptr lspop 1,v_len-u_len) ==0))
206 { lspref(d_LSDptr,d_len) = 1; d_len++; }
209 u->MSDptr = copy_loop_lsp(c_LSDptr,u->LSDptr,c_len);
210 while (mspref(u->MSDptr,0)==0) { msshrink(u->MSDptr); c_len--; }
212 v->MSDptr = copy_loop_lsp(d_LSDptr,v->LSDptr,d_len);
213 while (mspref(v->MSDptr,0)==0) { msshrink(v->MSDptr); d_len--; }
218 const cl_I xgcd (const cl_I& a, const cl_I& b, cl_I* u, cl_I* v)
219 { if (eq(a,1)) // a=1 -> g=1, (u,v)=(1,0)
220 { *u = 1; *v = 0; return a; }
221 if (eq(b,1)) // b=1 -> g=1, (u,v)=(0,1)
222 { *u = 0; *v = 1; return b; }
223 var sintL sA = (minusp(a) ? ~0 : 0); // Vorzeichen von A
224 var sintL sB = (minusp(b) ? ~0 : 0); // Vorzeichen von B
232 // Macro: erzeugt die NUDS zu (abs x), erniedrigt num_stack
233 #define I_abs_to_NUDS(x,zero_statement) \
234 I_to_NDS_1(x, x##_MSDptr = , x##_len = , x##_LSDptr = ); /* (nichtleere) NDS holen */\
235 if (x##_len == 0) { zero_statement } /* falls =0, fertig */\
236 if ((sintD)mspref(x##_MSDptr,0) < 0) /* falls <0, negieren: */\
237 { neg_loop_lsp(x##_LSDptr,x##_len); } \
238 if (mspref(x##_MSDptr,0) == 0) /* normalisieren (max. 1 Nulldigit entfernen) */\
239 { msshrink(x##_MSDptr); x##_len--; }
240 I_abs_to_NUDS(a, // (abs A) als NUDS erzeugen
241 // A=0 -> g=|B|, (u,v) = (0,sB)
242 { *u = 0; *v = (sB==0 ? cl_I(1) : cl_I(-1));
245 I_abs_to_NUDS(b, // (abs B) als NUDS erzeugen
246 // B=0 -> g=|A|, (u,v) = (sA,0)
247 { *u = (sA==0 ? cl_I(1) : cl_I(-1)); *v = 0;
250 // Jetzt ist a = a_MSDptr/a_len/a_LSDptr, b = b_MSDptr/b_len/b_LSDptr,
251 // beides NUDS, und a_len>0, b_len>0.
258 var uintD* divroomptr; // Platz für Divisionsergebnis
261 // Platz für uAa,uBa,uAb,uBb besorgen:
262 {var uintC u_len = b_len+1;
263 num_stack_alloc(u_len,,uAa.LSDptr=); uAa.MSDptr = uAa.LSDptr;
264 num_stack_alloc(u_len,,uAb.LSDptr=); uAb.MSDptr = uAb.LSDptr;
266 {var uintC u_len = a_len+1;
267 num_stack_alloc(u_len,,uBa.LSDptr=); uBa.MSDptr = uBa.LSDptr;
268 num_stack_alloc(u_len,,uBb.LSDptr=); uBb.MSDptr = uBb.LSDptr;
270 lsprefnext(uAa.MSDptr) = 1; uAa.len = 1; // uAa := 1
271 uBa.len = 0; // uBa := 0
272 uAb.len = 0; // uAb := 0
273 lsprefnext(uBb.MSDptr) = 1; uBb.len = 1; // uBb := 1
274 // Jetzt ist uAa = uAa.MSDptr/uAa.len/uAa.LSDptr,
275 // uBa = uBa.MSDptr/uBa.len/uBa.LSDptr,
276 // uAb = uAb.MSDptr/uAb.len/uAb.LSDptr,
277 // uBb = uBb.MSDptr/uBb.len/uBb.LSDptr,
279 // Platz für zwei Rechenregister besorgen, mit je max(a_len,b_len)+1 Digits:
280 {var uintL c_len = (uintL)(a_len>=b_len ? a_len : b_len) + 1;
281 num_stack_alloc(c_len,,c_LSDptr=);
282 num_stack_alloc(c_len,divroomptr=,d_LSDptr=);
283 // Jetzt ist ../c_len/c_LSDptr, ../c_len/d_LSDptr frei.
286 { // Hier a,b>0, beides NUDS.
287 // Vergleiche a und b:
288 if (a_len > b_len) goto a_greater_b; // a>b ?
290 { var cl_signean vergleich = compare_loop_msp(a_MSDptr,b_MSDptr,a_len);
291 if (vergleich > 0) goto a_greater_b; // a>b ?
292 if (vergleich == 0) break; // a=b ?
294 // a<b -> a,b vertauschen:
295 swap(uintD*, a_MSDptr,b_MSDptr);
296 swap(uintC, a_len,b_len);
297 swap(uintD*, a_LSDptr,b_LSDptr);
299 swap(DS, uAa,uAb); // und uAa und uAb vertauschen
300 swap(DS, uBa,uBb); // und uBa und uBb vertauschen
301 sA = ~sA; sB = ~sB; // und sA und sB umdrehen
303 // Hier a>b>0, beides NUDS.
304 // Entscheidung, ob Division oder Linearkombination:
305 { var uintD a_msd; // führende intDsize Bits von a
306 var uintD b_msd; // entsprechende Bits von b
308 var uintD a_nsd; // nächste intDsize Bits von a
309 var uintD b_nsd; // entsprechende Bits von b
311 { var uintC len_diff = a_len-b_len; // Längendifferenz
312 if (len_diff > 1) goto divide; // >=2 -> Bitlängendifferenz>intDsize -> dividieren
313 #define bitlendiff_limit (intDsize/2) // sollte >0,<intDsize sein
314 {var uintC a_msd_size;
315 a_msd = mspref(a_MSDptr,0); // führendes Digit von a
316 integerlengthD(a_msd,a_msd_size=); // dessen Bit-Länge (>0,<=intDsize) berechnen
317 b_msd = mspref(b_MSDptr,0);
319 {var uintDD b_msdd = // 2 führende Digits von b
321 ? highlowDD(b_msd, (b_len==1 ? 0 : mspref(b_MSDptr,1)))
324 // a_msd_size+intDsize - b_msdd_size >= bitlendiff_limit -> dividieren:
325 b_msd = lowD(b_msdd >> a_msd_size);
326 if (b_msd < (uintD)bit(intDsize-bitlendiff_limit)) goto divide;
328 b_nsd = lowD(highlowDD(lowD(b_msdd), (b_len<=2-len_diff ? 0 : mspref(b_MSDptr,2-len_diff))) >> a_msd_size);
331 {var uintDD a_msdd = // 2 führende Digits von a
332 highlowDD(a_msd, (a_len==1 ? 0 : mspref(a_MSDptr,1)));
333 a_msd = lowD(a_msdd >> a_msd_size);
335 a_nsd = lowD(highlowDD(lowD(a_msdd), (a_len<=2 ? 0 : mspref(a_MSDptr,2))) >> a_msd_size);
338 if (a_msd == b_msd) goto subtract;
341 { // a_msd_size - b_msd_size >= bitlendiff_limit -> dividieren:
342 if ((a_msd_size > bitlendiff_limit)
343 && (b_msd < (uintD)bit(a_msd_size-bitlendiff_limit))
346 // Entscheidung für Linearkombination ist gefallen.
347 // a_msd und b_msd so erweitern, daß a_msd die führenden
348 // intDsize Bits von a enthält:
349 {var uintC shiftcount = intDsize-a_msd_size; // Shiftcount nach links (>=0, <intDsize)
351 { a_msd = a_msd << shiftcount;
352 b_msd = b_msd << shiftcount;
354 { a_msd |= mspref(a_MSDptr,1) >> a_msd_size;
355 b_msd |= mspref(b_MSDptr,1) >> a_msd_size;
357 if (a_msd == b_msd) goto subtract;
360 { a_nsd = mspref(a_MSDptr,1);
361 b_nsd = mspref(b_MSDptr,1);
363 { a_nsd = a_nsd << shiftcount;
364 b_nsd = b_nsd << shiftcount;
366 { a_nsd |= mspref(a_MSDptr,2) >> a_msd_size;
367 b_nsd |= mspref(b_MSDptr,2) >> a_msd_size;
370 { a_nsd = 0; b_nsd = 0; }
375 { // a_msd_size+intDsize - b_msd_size >= bitlendiff_limit -> dividieren:
376 if ((a_msd_size >= bitlendiff_limit)
377 || (b_msd < (uintD)bit(a_msd_size+intDsize-bitlendiff_limit))
380 // Entscheidung für Linearkombination ist gefallen.
381 // a_msd und b_msd so erweitern, daß a_msd die führenden
382 // intDsize Bits von a enthält:
383 // 0 < a_msd_size < b_msd_size + bitlendiff_limit - intDsize <= bitlendiff_limit < intDsize.
384 a_msd = (a_msd << (intDsize-a_msd_size)) | (mspref(a_MSDptr,1) >> a_msd_size);
386 a_nsd = mspref(a_MSDptr,1) << (intDsize-a_msd_size);
387 b_nsd = b_msd << (intDsize-a_msd_size);
389 { a_nsd |= mspref(a_MSDptr,2) >> a_msd_size;
390 b_nsd |= mspref(b_MSDptr,1) >> a_msd_size;
393 b_msd = b_msd >> a_msd_size;
396 #undef bitlendiff_limit
398 // Nun ist a_msd = a' > b' = b_msd.
399 { // Euklid-Algorithmus auf den führenden Digits durchführen:
400 var partial_gcd_result likobi;
403 partial_gcd(highlowDD(a_msd,a_nsd),highlowDD(b_msd,b_nsd),&likobi); // liefert x1,y1,x2,y2
405 partial_gcd(a_msd,a_nsd,b_msd,b_nsd,&likobi); // liefert x1,y1,x2,y2
408 partial_gcd(a_msd,b_msd,&likobi); // liefert x1,y1,x2,y2, aber nur halb so gut
412 { // Ersetze (a,b) := (a-y1*b,b).
413 if (likobi.y1==1) goto subtract; // einfacherer Fall
414 // Dazu evtl. a um 1 Digit erweitern, so daß a_len=b_len+1:
415 if (a_len == b_len) { lsprefnext(a_MSDptr) = 0; a_len++; }
416 // und y1*b von a subtrahieren:
417 mspref(a_MSDptr,0) -= mulusub_loop_lsp(likobi.y1,b_LSDptr,a_LSDptr,b_len);
418 NUDS_likobi1_NUDS(&uAa,&uAb,likobi.y1); // uAa := uAa + y1 * uAb
419 NUDS_likobi1_NUDS(&uBa,&uBb,likobi.y1); // uBa := uBa + y1 * uBb
422 { // Ersetze (uAa,uAb) := (x1*uAa+y1*uAb,x2*uAa+y2*uAb) :
423 NUDS_likobi2_NUDS(&uAa,&uAb,&likobi,c_LSDptr,d_LSDptr);
424 // Ersetze (uBa,uBb) := (x1*uBa+y1*uBb,x2*uBa+y2*uBb) :
425 NUDS_likobi2_NUDS(&uBa,&uBb,&likobi,c_LSDptr,d_LSDptr);
426 // Ersetze (a,b) := (x1*a-y1*b,-x2*a+y2*b).
427 // Dazu evtl. b um 1 Digit erweitern, so daß a_len=b_len:
428 if (!(a_len==b_len)) { lsprefnext(b_MSDptr) = 0; b_len++; }
429 // c := x1*a-y1*b bilden:
430 mulu_loop_lsp(likobi.x1,a_LSDptr,c_LSDptr,a_len);
431 /* lspref(c_LSDptr,a_len) -= */
432 mulusub_loop_lsp(likobi.y1,b_LSDptr,c_LSDptr,a_len);
433 // d := -x2*a+y2*b bilden:
434 mulu_loop_lsp(likobi.y2,b_LSDptr,d_LSDptr,a_len);
435 /* lspref(d_LSDptr,a_len) -= */
436 mulusub_loop_lsp(likobi.x2,a_LSDptr,d_LSDptr,a_len);
437 // Wir wissen, daß 0 < c < b und 0 < d < a. Daher müßten
438 // lspref(c_LSDptr,a_len) und lspref(d_LSDptr,a_len) =0 sein.
439 // a := c und b := d kopieren:
440 copy_loop_lsp(c_LSDptr,a_LSDptr,a_len);
441 copy_loop_lsp(d_LSDptr,b_LSDptr,a_len);
443 while (mspref(b_MSDptr,0)==0) { msshrink(b_MSDptr); b_len--; }
446 { subtract: // Ersetze (a,b) := (a-b,b).
447 NUDS_likobi0_NUDS(&uAa,&uAb); // uAa := uAa + uAb
448 NUDS_likobi0_NUDS(&uBa,&uBb); // uBa := uBa + uBb
449 if (!( subfrom_loop_lsp(b_LSDptr,a_LSDptr,b_len) ==0))
450 // Übertrag nach b_len Stellen, muß also a_len=b_len+1 sein.
451 { mspref(a_MSDptr,0) -= 1; }
454 while (mspref(a_MSDptr,0)==0) { msshrink(a_MSDptr); a_len--; }
457 { divide: // Ersetze (a,b) := (b , a mod b).
458 {var uintD* old_a_LSDptr = a_LSDptr;
461 cl_UDS_divide(a_MSDptr,a_len,a_LSDptr,b_MSDptr,b_len,b_LSDptr, divroomptr, &q,&r);
462 a_MSDptr = b_MSDptr; a_len = b_len; a_LSDptr = b_LSDptr; // a := b
463 b_len = r.len; if (b_len==0) goto return_a_coeffsb; // b=0 -> fertig
464 b_LSDptr = old_a_LSDptr; // b übernimmt den vorherigen Platz von a
465 b_MSDptr = copy_loop_lsp(r.LSDptr,b_LSDptr,b_len); // b := r kopieren
466 // (uAa,uAb) := (uAb,uAa+q*uAb) :
468 { cl_UDS_mul(q.LSDptr,q.len,uAb.LSDptr,uAb.len,c_LSDptr); // q * uAb
470 c.LSDptr = c_LSDptr; c.len = q.len + uAb.len;
471 if (lspref(c_LSDptr,c.len-1)==0) { c.len--; } // normalisieren
472 NUDS_likobi0_NUDS(&uAa,&c); // zu uAa addieren
473 } // noch uAa,uAb vertauschen (später)
474 // (uBa,uBb) := (uBb,uBa+q*uBb) :
476 { cl_UDS_mul(q.LSDptr,q.len,uBb.LSDptr,uBb.len,c_LSDptr); // q * uBb
478 c.LSDptr = c_LSDptr; c.len = q.len + uBb.len;
479 if (lspref(c_LSDptr,c.len-1)==0) { c.len--; } // normalisieren
480 NUDS_likobi0_NUDS(&uBa,&c); // zu uBa addieren
481 } // noch uBa,uBb vertauschen (später)
482 goto a_greater_b_swap; // Nun ist a>b>0
485 // Nun ist a = b. Wähle diejenige der beiden Linearkombinationen
486 // a = uAa*sA * A + -uBa*sB * B
487 // b = -uAb*sA * A + uBb*sB * B
488 // die die betragsmäßig kleinsten Koeffizienten hat.
489 // Teste auf uBa < uBb. (Das kann auftreten, z.B. bei
490 // A=560014183, B=312839871 wird a=b=1, uAa < uAb, uBa < uBb.)
491 // Falls uBa = uBb, teste auf uAa < uAb. (Das kann auftreten, z.B. bei
492 // A=2, B=3 wird a=b=1, uAa < uAb, uBa = uBb.)
493 if (uBb.len > uBa.len) goto return_a_coeffsa;
494 if (uBb.len < uBa.len) goto return_a_coeffsb;
495 // (uBb.len == uBa.len)
496 { var cl_signean vergleich = compare_loop_msp(uBb.MSDptr,uBa.MSDptr,uBb.len);
497 if (vergleich > 0) goto return_a_coeffsa;
498 if (vergleich < 0) goto return_a_coeffsb;
500 if (uAb.len > uAa.len) goto return_a_coeffsa;
501 if (uAb.len < uAa.len) goto return_a_coeffsb;
502 // (uAb.len == uAa.len)
503 if (compare_loop_msp(uAb.MSDptr,uAa.MSDptr,uAb.len) > 0)
505 { // uAa mit Vorfaktor sA versehen:
506 lsprefnext(uAa.MSDptr) = 0; uAa.len++;
507 if (!(sA==0)) { neg_loop_lsp(uAa.LSDptr,uAa.len); }
508 // uBa mit Vorfaktor -sB versehen:
509 lsprefnext(uBa.MSDptr) = 0; uBa.len++;
510 if (sB==0) { neg_loop_lsp(uBa.LSDptr,uBa.len); }
511 *u = DS_to_I(uAa.MSDptr,uAa.len); // DS uAa als Vorfaktor von A
512 *v = DS_to_I(uBa.MSDptr,uBa.len); // DS uBa als Vorfaktor von B
516 { // uAb mit Vorfaktor -sA versehen:
517 lsprefnext(uAb.MSDptr) = 0; uAb.len++;
518 if (sA==0) { neg_loop_lsp(uAb.LSDptr,uAb.len); }
519 // uBb mit Vorfaktor sB versehen:
520 lsprefnext(uBb.MSDptr) = 0; uBb.len++;
521 if (!(sB==0)) { neg_loop_lsp(uBb.LSDptr,uBb.len); }
522 *u = DS_to_I(uAb.MSDptr,uAb.len); // DS uAb als Vorfaktor von A
523 *v = DS_to_I(uBb.MSDptr,uBb.len); // DS uBb als Vorfaktor von B
526 return NUDS_to_I(a_MSDptr,a_len); // NUDS a als ggT
530 #endif /* GCD_ALGO == 3 */