4 #include "base/cl_sysdep.h"
7 #include "integer/cl_I.h"
12 #include "cln/integer.h"
13 #include "base/digit/cl_D.h"
17 // Dasselbe wie partial_gcd(z1,z2,erg), nur daß z1 und z2 Doppelworte sind.
18 // Bevor im Ergebnis erg ein Überlauf eintritt, wird abgebrochen.
22 static inline uintDD muluDD_unchecked(uintD q, uintDD a)
24 return muluD(q,lowD(a)) + highlowDD_0(muluD_unchecked(q,highD(a)));
27 // Division: liefert min(floor(x / y), 2^intDsize-1).
28 // Vorausgesetzt wird, daß x >= 2 * y > 0.
29 static uintD floorDD (uintDD x, uintDD y)
31 // vgl. Algorithmus für divu_3232_3232().
33 if (y < ((uintDD)1 << intDsize)) {
35 q = ~(uintD)0; // instead of overflow
37 divuD(x, (uintD)y, q =, );
42 integerlengthD(highD(y), shift=);
43 // NB: 0 < shift < intDsize since 2^intDsize <= y < 2^(2*intDsize-1).
44 // Determine q := floor((x>>shift) / ((y>>shift)+1)).
45 var uintD y_shifted = y >> shift;
48 q = highD(x) >> shift;
50 divuD(highD(x) >> shift, y_shifted, q =, );
52 // May need to increment q at most twice.
54 var uintDD p = muluDD_unchecked(q,y);
57 throw runtime_exception();
69 throw runtime_exception();
76 void partial_gcd (uintDD z1, uintDD z2, partial_gcd_result* erg)
83 // Hier ist z1-y1>=z2+y2.
84 // Bestimme q := floor((z1-y1)/(z2+y2)) >= 1 :
86 var uintDD zaehler = z1 - (uintDD)y1;
87 var uintDD nenner = z2 + (uintDD)y2;
88 // z2+y2 <= z1-y1 < beta^2 !
89 if (x2 > (~x1) >> 3) // x1 + 8*x2 >= beta ?
91 if (y2 > (~y1) >> 3) // y1 + 8*y2 >= beta ?
93 // Ist floor(zaehler,8) >= nenner ? Wenn nein -> do_subtract_1.
94 if ((zaehler >> 3) < nenner)
97 var uintD q = floorDD(zaehler,nenner);
99 var uintDD qx = muluD(q,x2);
100 if (qx > (uintDD)(~x1)) {
101 // Choose a smaller value for q, to avoid overflow of x1.
105 var uintDD qy = muluD(q,y2);
106 if (qy > (uintDD)(~y1)) {
107 // Choose a smaller value for q, to avoid overflow of y1.
113 z1 -= muluDD_unchecked(q,z2);
117 // Checks to avoid overflow.
118 if (x2 > ~x1) goto done;
119 if (y2 > ~y1) goto done;
121 if (z1 < z2) throw runtime_exception();
123 // Now really subtract.
127 } while (z1 - (uintDD)y1 >= nenner);
130 if (z2 - (uintDD)x2 <= z1 + (uintDD)(x1-1)) goto done;
131 // Hier ist z2-x2>=z1+x1.
132 // Bestimme q := floor((z2-x2)/(z1+x1)) >= 1 :
134 var uintDD zaehler = z2 - (uintDD)x2;
135 var uintDD nenner = z1 + (uintDD)x1;
136 // z1+x1 <= z2-x2 < beta^2 !
137 if (x1 > (~x2) >> 3) // x2 + 8*x1 >= beta ?
139 if (y1 > (~y2) >> 3) // y2 + 8*y1 >= beta ?
141 // Ist floor(zaehler,8) >= nenner ? Wenn nein -> do_subtract_2.
142 if ((zaehler >> 3) < nenner)
145 var uintD q = floorDD(zaehler,nenner);
147 var uintDD qx = muluD(q,x1);
148 if (qx > (uintDD)(~x2)) {
149 // Choose a smaller value for q, to avoid overflow of x2.
153 var uintDD qy = muluD(q,y1);
154 if (qy > (uintDD)(~y2)) {
155 // Choose a smaller value for q, to avoid overflow of y2.
161 z2 -= muluDD_unchecked(q,z1);
165 // Checks to avoid overflow.
166 if (x1 > ~x2) goto done;
167 if (y1 > ~y2) goto done;
169 if (z2 < z1) throw runtime_exception();
171 // Now really subtract.
175 } while (z2 - (uintDD)x2 >= nenner);
178 if (z1 - (uintDD)y1 <= z2 + (uintDD)(y2-1)) goto done;
181 // Keine Subtraktion (ohne Überlauf) mehr möglich.
182 erg->x1 = x1; erg->y1 = y1; erg->x2 = x2; erg->y2 = y2; // Ergebnis
187 // Division: liefert min(floor(xhi|xlo / yhi|ylo), 2^intDsize-1).
188 // Vorausgesetzt wird, daß xhi|xlo >= 2 * yhi|ylo > 0.
189 static uintD floorDD (uintD xhi, uintD xlo, uintD yhi, uintD ylo)
191 // vgl. Algorithmus für divu_3232_3232().
195 q = ~(uintD)0; // instead of overflow
197 divuD(xhi,xlo, ylo, q =, );
202 integerlengthD(yhi, shift=);
203 // NB: 0 < shift < intDsize since 2^intDsize <= y < 2^(2*intDsize-1).
204 // Determine q := floor((x>>shift) / ((y>>shift)+1)).
205 var uintD y_shifted = (uintD)(yhi << (intDsize-shift)) | (ylo >> shift);
210 divuD(xhi >> shift, (uintD)(xhi << (intDsize-shift)) | (xlo >> shift),
214 // May need to increment q at most twice.
218 muluD(q,ylo, phi =, plo =);
219 muluD(q,yhi, , phi +=);
221 if ((xhi < phi) || ((xhi == phi) && (xlo < plo)))
222 throw runtime_exception();
229 if ((xhi > yhi) || ((xhi == yhi) && (xlo >= ylo))) {
235 if ((xhi > yhi) || ((xhi == yhi) && (xlo >= ylo))) {
242 if ((xhi > yhi) || ((xhi == yhi) && (xlo >= ylo)))
243 throw runtime_exception();
250 void partial_gcd (uintD z1hi, uintD z1lo, uintD z2hi, uintD z2lo, partial_gcd_result* erg)
257 // Hier ist z1-y1>=z2+y2.
258 // Bestimme q := floor((z1-y1)/(z2+y2)) >= 1 :
260 var uintD zaehlerhi = z1hi;
261 var uintD zaehlerlo = z1lo - y1;
262 if (zaehlerlo > z1lo) { zaehlerhi--; }
263 var uintD nennerhi = z2hi;
264 var uintD nennerlo = z2lo + y2;
265 if (nennerlo < z2lo) { nennerhi++; }
266 // z2+y2 <= z1-y1 < beta^2 !
267 if (x2 > (~x1) >> 3) // x1 + 8*x2 >= beta ?
269 if (y2 > (~y1) >> 3) // y1 + 8*y2 >= beta ?
271 // Ist floor(zaehler,8) >= nenner ? Wenn nein -> do_subtract_1.
272 if ((zaehlerhi >> 3) < nennerhi)
274 if ((zaehlerhi >> 3) == nennerhi)
275 if (((uintD)(zaehlerhi << (intDsize-3)) | (zaehlerlo >> 3)) < nennerlo)
278 var uintD q = floorDD(zaehlerhi,zaehlerlo,nennerhi,nennerlo);
284 muluD(q,x2, qxhi =, qx =);
285 if ((qxhi > 0) || (qx > ~x1)) {
286 // Choose a smaller value for q, to avoid overflow of x1.
293 muluD(q,y2, qyhi =, qy =);
294 if ((qyhi > 0) || (qy > ~y1)) {
295 // Choose a smaller value for q, to avoid overflow of y1.
305 muluD(q,z2lo, qzhi =, qzlo =);
306 muluD(q,z2hi, , qzhi +=);
315 // Checks to avoid overflow.
316 if (x2 > ~x1) goto done;
317 if (y2 > ~y1) goto done;
319 if (z1hi < z2hi) throw runtime_exception();
320 if (z1hi == z2hi) if (z1lo < z2lo) throw runtime_exception();
322 // Now really subtract.
329 var uintD z1dec_hi = z1hi;
330 var uintD z1dec_lo = z1lo - y1;
333 if (z1dec_hi < nennerhi)
335 if (z1dec_hi == nennerhi)
336 if (z1dec_lo < nennerlo)
342 var uintD z1inc_hi = z1hi;
343 var uintD z1inc_lo = z1lo + x1-1;
346 var uintD z2dec_hi = z2hi;
347 var uintD z2dec_lo = z2lo - x2;
350 if (z2dec_hi < z1inc_hi) goto done;
351 if (z2dec_hi == z1inc_hi) if (z2dec_lo <= z1inc_lo) goto done;
353 // Hier ist z2-x2>=z1+x1.
354 // Bestimme q := floor((z2-x2)/(z1+x1)) >= 1 :
356 var uintD zaehlerhi = z2hi;
357 var uintD zaehlerlo = z2lo - x2;
358 if (zaehlerlo > z2lo) { zaehlerhi--; }
359 var uintD nennerhi = z1hi;
360 var uintD nennerlo = z1lo + x1;
361 if (nennerlo < z1lo) { nennerhi++; }
362 // z1+x1 <= z2-x2 < beta^2 !
363 if (x1 > (~x2) >> 3) // x2 + 8*x1 >= beta ?
365 if (y1 > (~y2) >> 3) // y2 + 8*y1 >= beta ?
367 // Ist floor(zaehler,8) >= nenner ? Wenn nein -> do_subtract_2.
368 if ((zaehlerhi >> 3) < nennerhi)
370 if ((zaehlerhi >> 3) == nennerhi)
371 if (((uintD)(zaehlerhi << (intDsize-3)) | (zaehlerlo >> 3)) < nennerlo)
374 var uintD q = floorDD(zaehlerhi,zaehlerlo,nennerhi,nennerlo);
380 muluD(q,x1, qxhi =, qx =);
381 if ((qxhi > 0) || (qx > ~x2)) {
382 // Choose a smaller value for q, to avoid overflow of x2.
389 muluD(q,y1, qyhi =, qy =);
390 if ((qyhi > 0) || (qy > ~y2)) {
391 // Choose a smaller value for q, to avoid overflow of y2.
401 muluD(q,z1lo, qzhi =, qzlo =);
402 muluD(q,z1hi, , qzhi +=);
411 // Checks to avoid overflow.
412 if (x1 > ~x2) goto done;
413 if (y1 > ~y2) goto done;
415 if (z2hi < z1hi) throw runtime_exception();
416 if (z2hi == z1hi) if (z2lo < z1lo) throw runtime_exception();
418 // Now really subtract.
425 var uintD z2dec_hi = z2hi;
426 var uintD z2dec_lo = z2lo - x2;
429 if (z2dec_hi < nennerhi)
431 if (z2dec_hi == nennerhi)
432 if (z2dec_lo < nennerlo)
438 var uintD z2inc_hi = z2hi;
439 var uintD z2inc_lo = z2lo + y2-1;
442 var uintD z1dec_hi = z1hi;
443 var uintD z1dec_lo = z1lo - y1;
446 if (z1dec_hi < z2inc_hi) goto done;
447 if (z1dec_hi == z2inc_hi) if (z1dec_lo <= z2inc_lo) goto done;
451 // Keine Subtraktion (ohne Überlauf) mehr möglich.
452 erg->x1 = x1; erg->y1 = y1; erg->x2 = x2; erg->y2 = y2; // Ergebnis