Annotation of OpenXM_contrib/pari-2.2/src/kernel/none/mp.c, Revision 1.1.1.1
1.1 noro 1: /* $Id: mp.c,v 1.47 2001/09/28 17:16:45 xavier Exp $
2:
3: Copyright (C) 2000 The PARI group.
4:
5: This file is part of the PARI/GP package.
6:
7: PARI/GP is free software; you can redistribute it and/or modify it under the
8: terms of the GNU General Public License as published by the Free Software
9: Foundation. It is distributed in the hope that it will be useful, but WITHOUT
10: ANY WARRANTY WHATSOEVER.
11:
12: Check the License for details. You should have received a copy of it, along
13: with the package; see the file 'COPYING'. If not, write to the Free Software
14: Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15:
16: /***********************************************************************/
17: /** **/
18: /** MULTIPRECISION KERNEL **/
19: /** **/
20: /***********************************************************************/
21: /* most of the routines in this file are commented out in 68k */
22: /* version (#ifdef __M68K__) since they are defined in mp.s */
23: #include "pari.h"
24:
25: /* NOTE: arguments of "spec" routines (muliispec, addiispec, etc.) aren't
26: * GENs but pairs (long *a, long na) representing a list of digits (in basis
27: * BITS_IN_LONG) : a[0], ..., a[na-1]. [ In ordre to facilitate splitting: no
28: * need to reintroduce codewords ]
29: * Use speci(a,na) to visualize the coresponding GEN.
30: */
31:
32: /* z2 := z1[imin..imax].f shifted left sh bits (feeding f from the right) */
33: #define shift_left2(z2,z1,imin,imax,f, sh,m) {\
34: register ulong _l,_k = ((ulong)f)>>m;\
35: GEN t1 = z1 + imax, t2 = z2 + imax, T = z1 + imin;\
36: while (t1 > T) {\
37: _l = *t1--;\
38: *t2-- = (_l<<(ulong)sh) | _k;\
39: _k = _l>>(ulong)m;\
40: }\
41: *t2 = (*t1<<(ulong)sh) | _k;\
42: }
43: #define shift_left(z2,z1,imin,imax,f, sh) {\
44: register const ulong _m = BITS_IN_LONG - sh;\
45: shift_left2((z2),(z1),(imin),(imax),(f),(sh),(_m));\
46: }
47:
48: /* z2 := f.z1[imin..imax-1] shifted right sh bits (feeding f from the left) */
49: #define shift_right2(z2,z1,imin,imax,f, sh,m) {\
50: register GEN t1 = z1 + imin, t2 = z2 + imin, T = z1 + imax;\
51: register ulong _k,_l = *t1++;\
52: *t2++ = (_l>>(ulong)sh) | ((ulong)f<<(ulong)m);\
53: while (t1 < T) {\
54: _k = _l<<(ulong)m; _l = *t1++;\
55: *t2++ = (_l>>(ulong)sh) | _k;\
56: }\
57: }
58: #define shift_right(z2,z1,imin,imax,f, sh) {\
59: register const ulong _m = BITS_IN_LONG - sh;\
60: shift_right2((z2),(z1),(imin),(imax),(f),(sh),(_m));\
61: }
62:
63: #define MASK(x) (((ulong)(x)) & (LGEFINTBITS | SIGNBITS))
64: int
65: egalii(GEN x, GEN y)
66: {
67: long i;
68: if (MASK(x[1]) != MASK(y[1])) return 0;
69: i = lgefint(x)-1; while (i>1 && x[i]==y[i]) i--;
70: return i==1;
71: }
72: #undef MASK
73:
74: #ifdef INLINE
75: INLINE
76: #endif
77: GEN
78: addsispec(long s, GEN x, long nx)
79: {
80: GEN xd, zd = (GEN)avma;
81: long lz;
82: LOCAL_OVERFLOW;
83:
84: lz = nx+3; (void)new_chunk(lz);
85: xd = x + nx;
86: *--zd = addll(*--xd, s);
87: if (overflow)
88: for(;;)
89: {
90: if (xd == x) { *--zd = 1; break; } /* enlarge z */
91: *--zd = ((ulong)*--xd) + 1;
92: if (*zd) { lz--; break; }
93: }
94: else lz--;
95: while (xd > x) *--zd = *--xd;
96: *--zd = evalsigne(1) | evallgefint(lz);
97: *--zd = evaltyp(t_INT) | evallg(lz);
98: avma=(long)zd; return zd;
99: }
100:
101: #define swapspec(x,y, nx,ny) {long _a=nx;GEN _z=x; nx=ny; ny=_a; x=y; y=_z;}
102:
103: #ifdef INLINE
104: INLINE
105: #endif
106: GEN
107: addiispec(GEN x, GEN y, long nx, long ny)
108: {
109: GEN xd,yd,zd;
110: long lz;
111: LOCAL_OVERFLOW;
112:
113: if (nx < ny) swapspec(x,y, nx,ny);
114: if (ny == 1) return addsispec(*y,x,nx);
115: zd = (GEN)avma;
116: lz = nx+3; (void)new_chunk(lz);
117: xd = x + nx;
118: yd = y + ny;
119: *--zd = addll(*--xd, *--yd);
120: while (yd > y) *--zd = addllx(*--xd, *--yd);
121: if (overflow)
122: for(;;)
123: {
124: if (xd == x) { *--zd = 1; break; } /* enlarge z */
125: *--zd = ((ulong)*--xd) + 1;
126: if (*zd) { lz--; break; }
127: }
128: else lz--;
129: while (xd > x) *--zd = *--xd;
130: *--zd = evalsigne(1) | evallgefint(lz);
131: *--zd = evaltyp(t_INT) | evallg(lz);
132: avma=(long)zd; return zd;
133: }
134:
135: /* assume x >= y */
136: #ifdef INLINE
137: INLINE
138: #endif
139: GEN
140: subisspec(GEN x, long s, long nx)
141: {
142: GEN xd, zd = (GEN)avma;
143: long lz;
144: LOCAL_OVERFLOW;
145:
146: lz = nx+2; (void)new_chunk(lz);
147: xd = x + nx;
148: *--zd = subll(*--xd, s);
149: if (overflow)
150: for(;;)
151: {
152: *--zd = ((ulong)*--xd) - 1;
153: if (*xd) break;
154: }
155: if (xd == x)
156: while (*zd == 0) { zd++; lz--; } /* shorten z */
157: else
158: do *--zd = *--xd; while (xd > x);
159: *--zd = evalsigne(1) | evallgefint(lz);
160: *--zd = evaltyp(t_INT) | evallg(lz);
161: avma=(long)zd; return zd;
162: }
163:
164: /* assume x > y */
165: #ifdef INLINE
166: INLINE
167: #endif
168: GEN
169: subiispec(GEN x, GEN y, long nx, long ny)
170: {
171: GEN xd,yd,zd;
172: long lz;
173: LOCAL_OVERFLOW;
174:
175: if (ny==1) return subisspec(x,*y,nx);
176: zd = (GEN)avma;
177: lz = nx+2; (void)new_chunk(lz);
178: xd = x + nx;
179: yd = y + ny;
180: *--zd = subll(*--xd, *--yd);
181: while (yd > y) *--zd = subllx(*--xd, *--yd);
182: if (overflow)
183: for(;;)
184: {
185: *--zd = ((ulong)*--xd) - 1;
186: if (*xd) break;
187: }
188: if (xd == x)
189: while (*zd == 0) { zd++; lz--; } /* shorten z */
190: else
191: do *--zd = *--xd; while (xd > x);
192: *--zd = evalsigne(1) | evallgefint(lz);
193: *--zd = evaltyp(t_INT) | evallg(lz);
194: avma=(long)zd; return zd;
195: }
196:
197: #ifndef __M68K__
198:
199: /* prototype of positive small ints */
200: static long pos_s[] = {
201: evaltyp(t_INT) | m_evallg(3), evalsigne(1) | evallgefint(3), 0 };
202:
203: /* prototype of negative small ints */
204: static long neg_s[] = {
205: evaltyp(t_INT) | m_evallg(3), evalsigne(-1) | evallgefint(3), 0 };
206:
207: void
208: affir(GEN x, GEN y)
209: {
210: const long s=signe(x),ly=lg(y);
211: long lx,sh,i;
212:
213: if (!s)
214: {
215: y[1] = evalexpo(-bit_accuracy(ly));
216: y[2] = 0; return;
217: }
218:
219: lx=lgefint(x); sh=bfffo(x[2]);
220: y[1] = evalsigne(s) | evalexpo(bit_accuracy(lx)-sh-1);
221: if (sh)
222: {
223: if (lx>ly)
224: { shift_left(y,x,2,ly-1, x[ly],sh); }
225: else
226: {
227: for (i=lx; i<ly; i++) y[i]=0;
228: shift_left(y,x,2,lx-1, 0,sh);
229: }
230: return;
231: }
232:
233: if (lx>=ly)
234: for (i=2; i<ly; i++) y[i]=x[i];
235: else
236: {
237: for (i=2; i<lx; i++) y[i]=x[i];
238: for ( ; i<ly; i++) y[i]=0;
239: }
240: }
241:
242: void
243: affrr(GEN x, GEN y)
244: {
245: long lx,ly,i;
246:
247: y[1] = x[1]; if (!signe(x)) { y[2]=0; return; }
248:
249: lx=lg(x); ly=lg(y);
250: if (lx>=ly)
251: for (i=2; i<ly; i++) y[i]=x[i];
252: else
253: {
254: for (i=2; i<lx; i++) y[i]=x[i];
255: for ( ; i<ly; i++) y[i]=0;
256: }
257: }
258:
259: GEN
260: shifti(GEN x, long n)
261: {
262: const long s=signe(x);
263: long lx,ly,i,m;
264: GEN y;
265:
266: if (!s) return gzero;
267: if (!n) return icopy(x);
268: lx = lgefint(x);
269: if (n > 0)
270: {
271: GEN z = (GEN)avma;
272: long d = n>>TWOPOTBITS_IN_LONG;
273:
274: ly = lx+d; y = new_chunk(ly);
275: for ( ; d; d--) *--z = 0;
276: m = n & (BITS_IN_LONG-1);
277: if (!m) for (i=2; i<lx; i++) y[i]=x[i];
278: else
279: {
280: register const ulong sh = BITS_IN_LONG - m;
281: shift_left2(y,x, 2,lx-1, 0,m,sh);
282: i = ((ulong)x[2]) >> sh;
283: if (i) { ly++; y = new_chunk(1); y[2] = i; }
284: }
285: }
286: else
287: {
288: n = -n;
289: ly = lx - (n>>TWOPOTBITS_IN_LONG);
290: if (ly<3) return gzero;
291: y = new_chunk(ly);
292: m = n & (BITS_IN_LONG-1);
293: if (!m) for (i=2; i<ly; i++) y[i]=x[i];
294: else
295: {
296: shift_right(y,x, 2,ly, 0,m);
297: if (y[2] == 0)
298: {
299: if (ly==3) { avma = (long)(y+3); return gzero; }
300: ly--; avma = (long)(++y);
301: }
302: }
303: }
304: y[1] = evalsigne(s)|evallgefint(ly);
305: y[0] = evaltyp(t_INT)|evallg(ly); return y;
306: }
307:
308: GEN
309: mptrunc(GEN x)
310: {
311: long d,e,i,s,m;
312: GEN y;
313:
314: if (typ(x)==t_INT) return icopy(x);
315: if ((s=signe(x)) == 0 || (e=expo(x)) < 0) return gzero;
316: d = (e>>TWOPOTBITS_IN_LONG) + 3;
317: m = e & (BITS_IN_LONG-1);
318: if (d > lg(x)) err(precer, "mptrunc (precision loss in truncation)");
319:
320: y=cgeti(d); y[1] = evalsigne(s) | evallgefint(d);
321: if (++m == BITS_IN_LONG)
322: for (i=2; i<d; i++) y[i]=x[i];
323: else
324: {
325: register const ulong sh = BITS_IN_LONG - m;
326: shift_right2(y,x, 2,d,0, sh,m);
327: }
328: return y;
329: }
330:
331: /* integral part */
332: GEN
333: mpent(GEN x)
334: {
335: long d,e,i,lx,m;
336: GEN y;
337:
338: if (typ(x)==t_INT) return icopy(x);
339: if (signe(x) >= 0) return mptrunc(x);
340: if ((e=expo(x)) < 0) return stoi(-1);
341: d = (e>>TWOPOTBITS_IN_LONG) + 3;
342: m = e & (BITS_IN_LONG-1);
343: lx=lg(x); if (d>lx) err(precer, "mpent (precision loss in trucation)");
344: y = new_chunk(d);
345: if (++m == BITS_IN_LONG)
346: {
347: for (i=2; i<d; i++) y[i]=x[i];
348: i=d; while (i<lx && !x[i]) i++;
349: if (i==lx) goto END;
350: }
351: else
352: {
353: register const ulong sh = BITS_IN_LONG - m;
354: shift_right2(y,x, 2,d,0, sh,m);
355: if (x[d-1]<<m == 0)
356: {
357: i=d; while (i<lx && !x[i]) i++;
358: if (i==lx) goto END;
359: }
360: }
361: /* set y:=y+1 */
362: for (i=d-1; i>=2; i--) { y[i]++; if (y[i]) goto END; }
363: y=new_chunk(1); y[2]=1; d++;
364: END:
365: y[1] = evalsigne(-1) | evallgefint(d);
366: y[0] = evaltyp(t_INT) | evallg(d); return y;
367: }
368:
369: int
370: cmpsi(long x, GEN y)
371: {
372: ulong p;
373:
374: if (!x) return -signe(y);
375:
376: if (x>0)
377: {
378: if (signe(y)<=0) return 1;
379: if (lgefint(y)>3) return -1;
380: p=y[2]; if (p == (ulong)x) return 0;
381: return p < (ulong)x ? 1 : -1;
382: }
383:
384: if (signe(y)>=0) return -1;
385: if (lgefint(y)>3) return 1;
386: p=y[2]; if (p == (ulong)-x) return 0;
387: return p < (ulong)(-x) ? -1 : 1;
388: }
389:
390: int
391: cmpii(GEN x, GEN y)
392: {
393: const long sx = signe(x), sy = signe(y);
394: long lx,ly,i;
395:
396: if (sx<sy) return -1;
397: if (sx>sy) return 1;
398: if (!sx) return 0;
399:
400: lx=lgefint(x); ly=lgefint(y);
401: if (lx>ly) return sx;
402: if (lx<ly) return -sx;
403: i=2; while (i<lx && x[i]==y[i]) i++;
404: if (i==lx) return 0;
405: return ((ulong)x[i] > (ulong)y[i]) ? sx : -sx;
406: }
407:
408: int
409: cmprr(GEN x, GEN y)
410: {
411: const long sx = signe(x), sy = signe(y);
412: long ex,ey,lx,ly,lz,i;
413:
414: if (sx<sy) return -1;
415: if (sx>sy) return 1;
416: if (!sx) return 0;
417:
418: ex=expo(x); ey=expo(y);
419: if (ex>ey) return sx;
420: if (ex<ey) return -sx;
421:
422: lx=lg(x); ly=lg(y); lz = (lx<ly)?lx:ly;
423: i=2; while (i<lz && x[i]==y[i]) i++;
424: if (i<lz) return ((ulong)x[i] > (ulong)y[i]) ? sx : -sx;
425: if (lx>=ly)
426: {
427: while (i<lx && !x[i]) i++;
428: return (i==lx) ? 0 : sx;
429: }
430: while (i<ly && !y[i]) i++;
431: return (i==ly) ? 0 : -sx;
432: }
433:
434: GEN
435: addss(long x, long y)
436: {
437: if (!x) return stoi(y);
438: if (x>0) { pos_s[2] = x; return addsi(y,pos_s); }
439: neg_s[2] = -x; return addsi(y,neg_s);
440: }
441:
442: GEN
443: addsi(long x, GEN y)
444: {
445: long sx,sy,ly;
446: GEN z;
447:
448: if (!x) return icopy(y);
449: sy=signe(y); if (!sy) return stoi(x);
450: if (x<0) { sx=-1; x=-x; } else sx=1;
451: if (sx==sy)
452: {
453: z = addsispec(x,y+2, lgefint(y)-2);
454: setsigne(z,sy); return z;
455: }
456: ly=lgefint(y);
457: if (ly==3)
458: {
459: const long d = y[2] - x;
460: if (!d) return gzero;
461: z=cgeti(3);
462: if (y[2] < 0 || d > 0) {
463: z[1] = evalsigne(sy) | evallgefint(3);
464: z[2] = d;
465: }
466: else {
467: z[1] = evalsigne(-sy) | evallgefint(3);
468: z[2] =-d;
469: }
470: return z;
471: }
472: z = subisspec(y+2,x, ly-2);
473: setsigne(z,sy); return z;
474: }
475:
476: GEN
477: addii(GEN x, GEN y)
478: {
479: long sx = signe(x);
480: long sy = signe(y);
481: long lx,ly;
482: GEN z;
483:
484: if (!sx) return sy? icopy(y): gzero;
485: if (!sy) return icopy(x);
486: lx=lgefint(x); ly=lgefint(y);
487:
488: if (sx==sy)
489: z = addiispec(x+2,y+2,lx-2,ly-2);
490: else
491: { /* sx != sy */
492: long i = lx - ly;
493: if (i==0)
494: {
495: i = absi_cmp(x,y);
496: if (!i) return gzero;
497: }
498: if (i<0) { sx=sy; swapspec(x,y, lx,ly); } /* ensure |x| >= |y| */
499: z = subiispec(x+2,y+2,lx-2,ly-2);
500: }
501: setsigne(z,sx); return z;
502: }
503:
504: GEN
505: addsr(long x, GEN y)
506: {
507: if (!x) return rcopy(y);
508: if (x>0) { pos_s[2]=x; return addir(pos_s,y); }
509: neg_s[2] = -x; return addir(neg_s,y);
510: }
511:
512: GEN
513: addir(GEN x, GEN y)
514: {
515: long e,l,ly;
516: GEN z;
517:
518: if (!signe(x)) return rcopy(y);
519: e = expo(y)-expi(x);
520: if (!signe(y))
521: {
522: #if 0
523: if (e>0) err(adder3);
524: #else /* design issue: make 0.0 "absorbing" */
525: if (e>0) return rcopy(y);
526: #endif
527: z = cgetr(3 + ((-e)>>TWOPOTBITS_IN_LONG));
528: affir(x,z); return z;
529: }
530:
531: ly=lg(y);
532: if (e>0)
533: {
534: l = ly - (e>>TWOPOTBITS_IN_LONG);
535: if (l<3) return rcopy(y);
536: }
537: else l = ly + ((-e)>>TWOPOTBITS_IN_LONG)+1;
538: z=cgetr(l); affir(x,z); y=addrr(z,y);
539: z = y+l; ly = lg(y); while (ly--) z[ly] = y[ly];
540: avma=(long)z; return z;
541: }
542:
543: GEN
544: addrr(GEN x, GEN y)
545: {
546: long sx=signe(x),sy=signe(y),ex=expo(x),ey=expo(y);
547: long e,m,flag,i,j,f2,lx,ly,lz;
548: GEN z;
549: LOCAL_OVERFLOW;
550:
551: e = ey-ex;
552: if (!sy)
553: {
554: if (!sx)
555: {
556: if (e > 0) ex=ey;
557: z=cgetr(3); z[1]=evalexpo(ex); z[2]=0; return z;
558: }
559: if (e > 0) { z=cgetr(3); z[1]=evalexpo(ey); z[2]=0; return z; }
560: lz = 3 + ((-e)>>TWOPOTBITS_IN_LONG);
561: lx = lg(x); if (lz>lx) lz=lx;
562: z = cgetr(lz); while(--lz) z[lz]=x[lz];
563: return z;
564: }
565: if (!sx)
566: {
567: if (e < 0) { z=cgetr(3); z[1]=evalexpo(ex); z[2]=0; return z; }
568: lz = 3 + (e>>TWOPOTBITS_IN_LONG);
569: ly = lg(y); if (lz>ly) lz=ly;
570: z = cgetr(lz); while (--lz) z[lz]=y[lz];
571: return z;
572: }
573:
574: if (e < 0) { z=x; x=y; y=z; ey=ex; i=sx; sx=sy; sy=i; e=-e; }
575: /* now ey >= ex */
576: lx=lg(x); ly=lg(y);
577: if (e)
578: {
579: long d = e >> TWOPOTBITS_IN_LONG, l = ly-d;
580: if (l > lx) { flag=0; lz=lx+d+1; }
581: else if (l > 2) { flag=1; lz=ly; lx=l; }
582: else return rcopy(y);
583: m = e & (BITS_IN_LONG-1);
584: }
585: else
586: {
587: if (lx > ly) lx = ly;
588: flag=2; lz=lx; m=0;
589: }
590: z = cgetr(lz);
591: if (m)
592: { /* shift right x m bits */
593: const ulong sh = BITS_IN_LONG-m;
594: GEN p1 = x; x = new_chunk(lx+1);
595: shift_right2(x,p1,2,lx, 0,m,sh);
596: if (flag==0)
597: {
598: x[lx] = p1[lx-1] << sh;
599: if (x[lx]) { flag = 3; lx++; }
600: }
601: }
602:
603: if (sx==sy)
604: { /* addition */
605: i=lz-1; avma = (long)z;
606: if (flag==0) { z[i] = y[i]; i--; }
607: overflow=0;
608: for (j=lx-1; j>=2; i--,j--) z[i] = addllx(x[j],y[i]);
609:
610: if (overflow)
611: for (;;)
612: {
613: if (i==1)
614: {
615: shift_right(z,z, 2,lz, 1,1);
616: ey=evalexpo(ey+1);
617: z[1] = evalsigne(sx) | ey; return z;
618: }
619: z[i] = y[i]+1; if (z[i--]) break;
620: }
621: for (; i>=2; i--) z[i]=y[i];
622: z[1] = evalsigne(sx) | evalexpo(ey); return z;
623: }
624:
625: /* subtraction */
626: if (e) f2 = 1;
627: else
628: {
629: i=2; while (i<lx && x[i]==y[i]) i++;
630: if (i==lx)
631: {
632: avma = (long)(z+lz);
633: e = evalexpo(ey - bit_accuracy(lx));
634: z=cgetr(3); z[1]=e; z[2]=0; return z;
635: }
636: f2 = ((ulong)y[i] > (ulong)x[i]);
637: }
638:
639: /* result is non-zero. f2 = (y > x) */
640: i=lz-1;
641: if (f2)
642: {
643: if (flag==0) { z[i] = y[i]; i--; }
644: j=lx-1; z[i] = subll(y[i],x[j]); i--; j--;
645: for (; j>=2; i--,j--) z[i] = subllx(y[i],x[j]);
646: if (overflow)
647: for (;;) { z[i] = y[i]-1; if (y[i--]) break; }
648: for (; i>=2; i--) z[i]=y[i];
649: sx = sy;
650: }
651: else
652: {
653: if (flag==0) { z[i] = -y[i]; i--; overflow=1; } else overflow=0;
654: for (; i>=2; i--) z[i]=subllx(x[i],y[i]);
655: }
656:
657: x = z+2; i=0; while (!x[i]) i++;
658: lz -= i; z += i; m = bfffo(z[2]);
659: e = evalexpo(ey - (m | (i<<TWOPOTBITS_IN_LONG)));
660: if (m) shift_left(z,z,2,lz-1, 0,m);
661: z[1] = evalsigne(sx) | e;
662: z[0] = evaltyp(t_REAL) | evallg(lz);
663: avma = (long)z; return z;
664: }
665:
666: GEN
667: mulss(long x, long y)
668: {
669: long s,p1;
670: GEN z;
671: LOCAL_HIREMAINDER;
672:
673: if (!x || !y) return gzero;
674: if (x<0) { s = -1; x = -x; } else s=1;
675: if (y<0) { s = -s; y = -y; }
676: p1 = mulll(x,y);
677: if (hiremainder)
678: {
679: z=cgeti(4); z[1] = evalsigne(s) | evallgefint(4);
680: z[2]=hiremainder; z[3]=p1; return z;
681: }
682: z=cgeti(3); z[1] = evalsigne(s) | evallgefint(3);
683: z[2]=p1; return z;
684: }
685: #endif
686:
687: GEN
688: muluu(ulong x, ulong y)
689: {
690: long p1;
691: GEN z;
692: LOCAL_HIREMAINDER;
693:
694: if (!x || !y) return gzero;
695: p1 = mulll(x,y);
696: if (hiremainder)
697: {
698: z=cgeti(4); z[1] = evalsigne(1) | evallgefint(4);
699: z[2]=hiremainder; z[3]=p1; return z;
700: }
701: z=cgeti(3); z[1] = evalsigne(1) | evallgefint(3);
702: z[2]=p1; return z;
703: }
704:
705: /* assume ny > 0 */
706: #ifdef INLINE
707: INLINE
708: #endif
709: GEN
710: mulsispec(long x, GEN y, long ny)
711: {
712: GEN yd, z = (GEN)avma;
713: long lz = ny+3;
714: LOCAL_HIREMAINDER;
715:
716: (void)new_chunk(lz);
717: yd = y + ny; *--z = mulll(x, *--yd);
718: while (yd > y) *--z = addmul(x,*--yd);
719: if (hiremainder) *--z = hiremainder; else lz--;
720: *--z = evalsigne(1) | evallgefint(lz);
721: *--z = evaltyp(t_INT) | evallg(lz);
722: avma=(long)z; return z;
723: }
724:
725: GEN
726: mului(ulong x, GEN y)
727: {
728: long s = signe(y);
729: GEN z;
730:
731: if (!s || !x) return gzero;
732: z = mulsispec(x, y+2, lgefint(y)-2);
733: setsigne(z,s); return z;
734: }
735:
736: /* a + b*Y, assume Y >= 0, 0 <= a,b <= VERYBIGINT */
737: GEN
738: addsmulsi(long a, long b, GEN Y)
739: {
740: GEN yd,y,z;
741: long ny,lz;
742: LOCAL_HIREMAINDER;
743: LOCAL_OVERFLOW;
744:
745: if (!signe(Y)) return stoi(a);
746:
747: y = Y+2; z = (GEN)avma;
748: ny = lgefint(Y)-2;
749: lz = ny+3;
750:
751: (void)new_chunk(lz);
752: yd = y + ny; *--z = addll(a, mulll(b, *--yd));
753: if (overflow) hiremainder++; /* can't overflow */
754: while (yd > y) *--z = addmul(b,*--yd);
755: if (hiremainder) *--z = hiremainder; else lz--;
756: *--z = evalsigne(1) | evallgefint(lz);
757: *--z = evaltyp(t_INT) | evallg(lz);
758: avma=(long)z; return z;
759: }
760:
761: #ifndef __M68K__
762:
763: GEN
764: mulsi(long x, GEN y)
765: {
766: long s = signe(y);
767: GEN z;
768:
769: if (!s || !x) return gzero;
770: if (x<0) { s = -s; x = -x; }
771: z = mulsispec(x, y+2, lgefint(y)-2);
772: setsigne(z,s); return z;
773: }
774:
775: GEN
776: mulsr(long x, GEN y)
777: {
778: long lx,i,s,garde,e,sh,m;
779: GEN z;
780: LOCAL_HIREMAINDER;
781:
782: if (!x) return gzero;
783: s = signe(y);
784: if (!s)
785: {
786: if (x<0) x = -x;
787: e = y[1] + (BITS_IN_LONG-1)-bfffo(x);
788: if (e & ~EXPOBITS) err(muler2);
789: z=cgetr(3); z[1]=e; z[2]=0; return z;
790: }
791: if (x<0) { s = -s; x = -x; }
792: if (x==1) { z=rcopy(y); setsigne(z,s); return z; }
793:
794: lx = lg(y); e = expo(y);
795: z=cgetr(lx); y--; garde=mulll(x,y[lx]);
796: for (i=lx-1; i>=3; i--) z[i]=addmul(x,y[i]);
797: z[2]=hiremainder;
798:
799: sh = bfffo(hiremainder); m = BITS_IN_LONG-sh;
800: if (sh) shift_left2(z,z, 2,lx-1, garde,sh,m);
801: e = evalexpo(m+e);
802: z[1] = evalsigne(s) | e; return z;
803: }
804:
805: GEN
806: mulrr(GEN x, GEN y)
807: {
808: long sx = signe(x), sy = signe(y);
809: long i,j,ly,lz,lzz,e,flag,p1;
810: ulong garde;
811: GEN z, y1;
812: LOCAL_HIREMAINDER;
813: LOCAL_OVERFLOW;
814:
815: e = evalexpo(expo(x)+expo(y));
816: if (!sx || !sy) { z=cgetr(3); z[1]=e; z[2]=0; return z; }
817: if (sy<0) sx = -sx;
818: lz=lg(x); ly=lg(y);
819: if (lz>ly) { lz=ly; z=x; x=y; y=z; flag=1; } else flag = (lz!=ly);
820: z=cgetr(lz); z[1] = evalsigne(sx) | e;
821: if (lz==3)
822: {
823: if (flag)
824: {
825: (void)mulll(x[2],y[3]);
826: garde = addmul(x[2],y[2]);
827: }
828: else
829: garde = mulll(x[2],y[2]);
830: if ((long)hiremainder<0) { z[2]=hiremainder; z[1]++; }
831: else z[2]=(hiremainder<<1) | (garde>>(BITS_IN_LONG-1));
832: return z;
833: }
834:
835: if (flag) { (void)mulll(x[2],y[lz]); garde = hiremainder; } else garde=0;
836: lzz=lz-1; p1=x[lzz];
837: if (p1)
838: {
839: (void)mulll(p1,y[3]);
840: garde = addll(addmul(p1,y[2]), garde);
841: z[lzz] = overflow+hiremainder;
842: }
843: else z[lzz]=0;
844: for (j=lz-2, y1=y-j; j>=3; j--)
845: {
846: p1 = x[j]; y1++;
847: if (p1)
848: {
849: (void)mulll(p1,y1[lz+1]);
850: garde = addll(addmul(p1,y1[lz]), garde);
851: for (i=lzz; i>j; i--)
852: {
853: hiremainder += overflow;
854: z[i] = addll(addmul(p1,y1[i]), z[i]);
855: }
856: z[j] = hiremainder+overflow;
857: }
858: else z[j]=0;
859: }
860: p1=x[2]; y1++;
861: garde = addll(mulll(p1,y1[lz]), garde);
862: for (i=lzz; i>2; i--)
863: {
864: hiremainder += overflow;
865: z[i] = addll(addmul(p1,y1[i]), z[i]);
866: }
867: z[2] = hiremainder+overflow;
868: if (z[2] < 0) z[1]++;
869: else
870: shift_left(z,z,2,lzz,garde, 1);
871: return z;
872: }
873:
874: GEN
875: mulir(GEN x, GEN y)
876: {
877: long sx=signe(x),sy,lz,lzz,ey,e,p1,i,j;
878: ulong garde;
879: GEN z,y1;
880: LOCAL_HIREMAINDER;
881: LOCAL_OVERFLOW;
882:
883: if (!sx) return gzero;
884: if (!is_bigint(x)) return mulsr(itos(x),y);
885: sy=signe(y); ey=expo(y);
886: if (!sy)
887: {
888: e = evalexpo(expi(x)+ey);
889: z=cgetr(3); z[1]=e; z[2]=0; return z;
890: }
891:
892: if (sy<0) sx = -sx;
893: lz=lg(y); z=cgetr(lz);
894: y1=cgetr(lz+1);
895: affir(x,y1); x=y; y=y1;
896: e = evalexpo(expo(y)+ey);
897: z[1] = evalsigne(sx) | e;
898: if (lz==3)
899: {
900: (void)mulll(x[2],y[3]);
901: garde=addmul(x[2],y[2]);
902: if ((long)hiremainder < 0) { z[2]=hiremainder; z[1]++; }
903: else z[2]=(hiremainder<<1) | (garde>>(BITS_IN_LONG-1));
904: avma=(long)z; return z;
905: }
906:
907: (void)mulll(x[2],y[lz]); garde=hiremainder;
908: lzz=lz-1; p1=x[lzz];
909: if (p1)
910: {
911: (void)mulll(p1,y[3]);
912: garde=addll(addmul(p1,y[2]),garde);
913: z[lzz] = overflow+hiremainder;
914: }
915: else z[lzz]=0;
916: for (j=lz-2, y1=y-j; j>=3; j--)
917: {
918: p1=x[j]; y1++;
919: if (p1)
920: {
921: (void)mulll(p1,y1[lz+1]);
922: garde = addll(addmul(p1,y1[lz]), garde);
923: for (i=lzz; i>j; i--)
924: {
925: hiremainder += overflow;
926: z[i] = addll(addmul(p1,y1[i]), z[i]);
927: }
928: z[j] = hiremainder+overflow;
929: }
930: else z[j]=0;
931: }
932: p1=x[2]; y1++;
933: garde = addll(mulll(p1,y1[lz]), garde);
934: for (i=lzz; i>2; i--)
935: {
936: hiremainder += overflow;
937: z[i] = addll(addmul(p1,y1[i]), z[i]);
938: }
939: z[2] = hiremainder+overflow;
940: if (z[2] < 0) z[1]++;
941: else
942: shift_left(z,z,2,lzz,garde, 1);
943: avma=(long)z; return z;
944: }
945:
946: /* written by Bruno Haible following an idea of Robert Harley */
947: long
948: vals(ulong z)
949: {
950: static char tab[64]={-1,0,1,12,2,6,-1,13,3,-1,7,-1,-1,-1,-1,14,10,4,-1,-1,8,-1,-1,25,-1,-1,-1,-1,-1,21,27,15,31,11,5,-1,-1,-1,-1,-1,9,-1,-1,24,-1,-1,20,26,30,-1,-1,-1,-1,23,-1,19,29,-1,22,18,28,17,16,-1};
951: #ifdef LONG_IS_64BIT
952: long s;
953: #endif
954:
955: if (!z) return -1;
956: #ifdef LONG_IS_64BIT
957: if (! (z&0xffffffff)) { s = 32; z >>=32; } else s = 0;
958: #endif
959: z |= ~z + 1;
960: z += z << 4;
961: z += z << 6;
962: z ^= z << 16; /* or z -= z<<16 */
963: #ifdef LONG_IS_64BIT
964: return s + tab[(z&0xffffffff)>>26];
965: #else
966: return tab[z>>26];
967: #endif
968: }
969:
970: GEN
971: modss(long x, long y)
972: {
973: LOCAL_HIREMAINDER;
974:
975: if (!y) err(moder1);
976: if (y<0) y=-y;
977: hiremainder=0; divll(labs(x),y);
978: if (!hiremainder) return gzero;
979: return (x < 0) ? stoi(y-hiremainder) : stoi(hiremainder);
980: }
981:
982: GEN
983: resss(long x, long y)
984: {
985: LOCAL_HIREMAINDER;
986:
987: if (!y) err(reser1);
988: hiremainder=0; divll(labs(x),labs(y));
989: return (x < 0) ? stoi(-((long)hiremainder)) : stoi(hiremainder);
990: }
991:
992: GEN
993: divsi(long x, GEN y)
994: {
995: long p1, s = signe(y);
996: LOCAL_HIREMAINDER;
997:
998: if (!s) err(diver2);
999: if (!x || lgefint(y)>3 || ((long)y[2])<0)
1000: {
1001: hiremainder=x; SAVE_HIREMAINDER; return gzero;
1002: }
1003: hiremainder=0; p1=divll(labs(x),y[2]);
1004: if (x<0) { hiremainder = -((long)hiremainder); p1 = -p1; }
1005: if (s<0) p1 = -p1;
1006: SAVE_HIREMAINDER; return stoi(p1);
1007: }
1008: #endif
1009:
1010: GEN
1011: modui(ulong x, GEN y)
1012: {
1013: LOCAL_HIREMAINDER;
1014:
1015: if (!signe(y)) err(diver2);
1016: if (!x || lgefint(y)>3) hiremainder=x;
1017: else
1018: {
1019: hiremainder=0; (void)divll(x,y[2]);
1020: }
1021: return utoi(hiremainder);
1022: }
1023:
1024: ulong
1025: umodiu(GEN y, ulong x)
1026: {
1027: long sy=signe(y),ly,i;
1028: LOCAL_HIREMAINDER;
1029:
1030: if (!x) err(diver4);
1031: if (!sy) return 0;
1032: ly = lgefint(y);
1033: if (x <= (ulong)y[2]) hiremainder=0;
1034: else
1035: {
1036: if (ly==3) return (sy > 0)? (ulong)y[2]: x - (ulong)y[2];
1037: hiremainder=y[2]; ly--; y++;
1038: }
1039: for (i=2; i<ly; i++) (void)divll(y[i],x);
1040: if (!hiremainder) return 0;
1041: return (sy > 0)? hiremainder: x - hiremainder;
1042: }
1043:
1044: GEN
1045: modiu(GEN y, ulong x) { return utoi(umodiu(y,x)); }
1046:
1047: #ifndef __M68K__
1048:
1049: GEN
1050: modsi(long x, GEN y)
1051: {
1052: long s = signe(y);
1053: GEN p1;
1054: LOCAL_HIREMAINDER;
1055:
1056: if (!s) err(diver2);
1057: if (!x || lgefint(y)>3 || ((long)y[2])<0) hiremainder=x;
1058: else
1059: {
1060: hiremainder=0; (void)divll(labs(x),y[2]);
1061: if (x<0) hiremainder = -((long)hiremainder);
1062: }
1063: if (!hiremainder) return gzero;
1064: if (x>0) return stoi(hiremainder);
1065: if (s<0)
1066: { setsigne(y,1); p1=addsi(hiremainder,y); setsigne(y,-1); }
1067: else
1068: p1=addsi(hiremainder,y);
1069: return p1;
1070: }
1071:
1072: GEN
1073: divis(GEN y, long x)
1074: {
1075: long sy=signe(y),ly,s,i;
1076: GEN z;
1077: LOCAL_HIREMAINDER;
1078:
1079: if (!x) err(diver4);
1080: if (!sy) { hiremainder=0; SAVE_HIREMAINDER; return gzero; }
1081: if (x<0) { s = -sy; x = -x; } else s=sy;
1082:
1083: ly = lgefint(y);
1084: if ((ulong)x <= (ulong)y[2]) hiremainder=0;
1085: else
1086: {
1087: if (ly==3) { hiremainder=itos(y); SAVE_HIREMAINDER; return gzero; }
1088: hiremainder=y[2]; ly--; y++;
1089: }
1090: z = cgeti(ly); z[1] = evallgefint(ly) | evalsigne(s);
1091: for (i=2; i<ly; i++) z[i]=divll(y[i],x);
1092: if (sy<0) hiremainder = - ((long)hiremainder);
1093: SAVE_HIREMAINDER; return z;
1094: }
1095:
1096: GEN
1097: divir(GEN x, GEN y)
1098: {
1099: GEN xr,z;
1100: long av,ly;
1101:
1102: if (!signe(y)) err(diver5);
1103: if (!signe(x)) return gzero;
1104: ly=lg(y); z=cgetr(ly); av=avma; xr=cgetr(ly+1); affir(x,xr);
1105: affrr(divrr(xr,y),z); avma=av; return z;
1106: }
1107:
1108: GEN
1109: divri(GEN x, GEN y)
1110: {
1111: GEN yr,z;
1112: long av,lx,s=signe(y);
1113:
1114: if (!s) err(diver8);
1115: if (!signe(x))
1116: {
1117: const long ex = evalexpo(expo(x) - expi(y));
1118:
1119: if (ex<0) err(diver12);
1120: z=cgetr(3); z[1]=ex; z[2]=0; return z;
1121: }
1122: if (!is_bigint(y)) return divrs(x, s>0? y[2]: -y[2]);
1123:
1124: lx=lg(x); z=cgetr(lx);
1125: av=avma; yr=cgetr(lx+1); affir(y,yr);
1126: affrr(divrr(x,yr),z); avma=av; return z;
1127: }
1128:
1129: void
1130: diviiz(GEN x, GEN y, GEN z)
1131: {
1132: long av=avma,lz;
1133: GEN xr,yr;
1134:
1135: if (typ(z)==t_INT) { affii(divii(x,y),z); avma=av; return; }
1136: lz=lg(z); xr=cgetr(lz); affir(x,xr); yr=cgetr(lz); affir(y,yr);
1137: affrr(divrr(xr,yr),z); avma=av;
1138: }
1139:
1140: void
1141: mpdivz(GEN x, GEN y, GEN z)
1142: {
1143: long av=avma;
1144:
1145: if (typ(z)==t_INT)
1146: {
1147: if (typ(x)==t_REAL || typ(y)==t_REAL) err(divzer1);
1148: affii(divii(x,y),z); avma=av; return;
1149: }
1150: if (typ(x)==t_INT)
1151: {
1152: GEN xr,yr;
1153: long lz;
1154:
1155: if (typ(y)==t_REAL) { affrr(divir(x,y),z); avma=av; return; }
1156: lz=lg(z); xr=cgetr(lz); affir(x,xr); yr=cgetr(lz); affir(y,yr);
1157: affrr(divrr(xr,yr),z); avma=av; return;
1158: }
1159: if (typ(y)==t_REAL) { affrr(divrr(x,y),z); avma=av; return; }
1160: affrr(divri(x,y),z); avma=av;
1161: }
1162:
1163: GEN
1164: divsr(long x, GEN y)
1165: {
1166: long av,ly;
1167: GEN p1,z;
1168:
1169: if (!signe(y)) err(diver3);
1170: if (!x) return gzero;
1171: ly=lg(y); z=cgetr(ly); av=avma;
1172: p1=cgetr(ly+1); affsr(x,p1); affrr(divrr(p1,y),z);
1173: avma=av; return z;
1174: }
1175:
1176: GEN
1177: modii(GEN x, GEN y)
1178: {
1179: switch(signe(x))
1180: {
1181: case 0: return gzero;
1182: case 1: return resii(x,y);
1183: default:
1184: {
1185: long av = avma;
1186: (void)new_chunk(lgefint(y));
1187: x = resii(x,y); avma=av;
1188: if (x==gzero) return x;
1189: return subiispec(y+2,x+2,lgefint(y)-2,lgefint(x)-2);
1190: }
1191: }
1192: }
1193:
1194: void
1195: modiiz(GEN x, GEN y, GEN z)
1196: {
1197: const long av = avma;
1198: affii(modii(x,y),z); avma=av;
1199: }
1200:
1201: GEN
1202: divrs(GEN x, long y)
1203: {
1204: long i,lx,ex,garde,sh,s=signe(x);
1205: GEN z;
1206: LOCAL_HIREMAINDER;
1207:
1208: if (!y) err(diver6);
1209: if (!s)
1210: {
1211: z=cgetr(3); z[1] = x[1] - (BITS_IN_LONG-1) + bfffo(y);
1212: if (z[1]<0) err(diver7);
1213: z[2]=0; return z;
1214: }
1215: if (y<0) { s = -s; y = -y; }
1216: if (y==1) { z=rcopy(x); setsigne(z,s); return z; }
1217:
1218: z=cgetr(lx=lg(x)); hiremainder=0;
1219: for (i=2; i<lx; i++) z[i]=divll(x[i],y);
1220:
1221: /* we may have hiremainder != 0 ==> garde */
1222: garde=divll(0,y); sh=bfffo(z[2]); ex=evalexpo(expo(x)-sh);
1223: if (sh) shift_left(z,z, 2,lx-1, garde,sh);
1224: z[1] = evalsigne(s) | ex;
1225: return z;
1226: }
1227:
1228: GEN
1229: divrr(GEN x, GEN y)
1230: {
1231: long sx=signe(x), sy=signe(y), lx,ly,lz,e,i,j;
1232: ulong si,saux;
1233: GEN z,x1;
1234:
1235: if (!sy) err(diver9);
1236: e = evalexpo(expo(x) - expo(y));
1237: if (!sx) { z=cgetr(3); z[1]=e; z[2]=0; return z; }
1238: if (sy<0) sx = -sx;
1239: e = evalsigne(sx) | e;
1240: lx=lg(x); ly=lg(y);
1241: if (ly==3)
1242: {
1243: ulong k = x[2], l = (lx>3)? x[3]: 0;
1244: LOCAL_HIREMAINDER;
1245: if (k < (ulong)y[2]) e--;
1246: else
1247: {
1248: l >>= 1; if (k&1) l |= HIGHBIT;
1249: k >>= 1;
1250: }
1251: z=cgetr(3); z[1]=e;
1252: hiremainder=k; z[2]=divll(l,y[2]); return z;
1253: }
1254:
1255: lz = (lx<=ly)? lx: ly;
1256: x1 = (z=new_chunk(lz))-1;
1257: x1[1]=0; for (i=2; i<lz; i++) x1[i]=x[i];
1258: x1[lz] = (lx>lz)? x[lz]: 0;
1259: x=z; si=y[2]; saux=y[3];
1260: for (i=0; i<lz-1; i++)
1261: { /* x1 = x + (i-1) */
1262: ulong k,qp;
1263: LOCAL_HIREMAINDER;
1264: LOCAL_OVERFLOW;
1265: if ((ulong)x1[1] == si)
1266: {
1267: qp = MAXULONG; k=addll(si,x1[2]);
1268: }
1269: else
1270: {
1271: if ((ulong)x1[1] > si) /* can't happen if i=0 */
1272: {
1273: GEN y1 = y+1;
1274: overflow=0;
1275: for (j=lz-i; j>0; j--) x1[j] = subllx(x1[j],y1[j]);
1276: j=i; do x[--j]++; while (j && !x[j]);
1277: }
1278: hiremainder=x1[1]; overflow=0;
1279: qp=divll(x1[2],si); k=hiremainder;
1280: }
1281: if (!overflow)
1282: {
1283: long k3 = subll(mulll(qp,saux), x1[3]);
1284: long k4 = subllx(hiremainder,k);
1285: while (!overflow && k4) { qp--; k3=subll(k3,saux); k4=subllx(k4,si); }
1286: }
1287: j = lz-i+1;
1288: if (j<ly) (void)mulll(qp,y[j]); else { hiremainder=0 ; j=ly; }
1289: for (j--; j>1; j--)
1290: {
1291: x1[j] = subll(x1[j], addmul(qp,y[j]));
1292: hiremainder += overflow;
1293: }
1294: if ((ulong)x1[1] != hiremainder)
1295: {
1296: if ((ulong)x1[1] < hiremainder)
1297: {
1298: qp--;
1299: overflow=0; for (j=lz-i; j>1; j--) x1[j]=addllx(x1[j], y[j]);
1300: }
1301: else
1302: {
1303: x1[1] -= hiremainder;
1304: while (x1[1])
1305: {
1306: qp++; if (!qp) { j=i; do x[--j]++; while (j && !x[j]); }
1307: overflow=0; for (j=lz-i; j>1; j--) x1[j]=subllx(x1[j],y[j]);
1308: x1[1] -= overflow;
1309: }
1310: }
1311: }
1312: x1[1]=qp; x1++;
1313: }
1314: x1 = x-1; for (j=lz-1; j>=2; j--) x[j]=x1[j];
1315: if (*x) { shift_right(x,x, 2,lz, 1,1); } else e--;
1316: x[0]=evaltyp(t_REAL)|evallg(lz);
1317: x[1]=e; return x;
1318: }
1319: #endif /* !defined(__M68K__) */
1320: /* The following ones are not in mp.s (mulii is, with a different algorithm) */
1321:
1322: /* Integer division x / y:
1323: * if z = ONLY_REM return remainder, otherwise return quotient
1324: * if z != NULL set *z to remainder
1325: * *z is the last object on stack (and thus can be disposed of with cgiv
1326: * instead of gerepile)
1327: * If *z is zero, we put gzero here and no copy.
1328: * space needed: lx + ly
1329: */
1330: GEN
1331: dvmdii(GEN x, GEN y, GEN *z)
1332: {
1333: long sx=signe(x),sy=signe(y);
1334: long av,lx,ly,lz,i,j,sh,k,lq,lr;
1335: ulong si,qp,saux, *xd,*rd,*qd;
1336: GEN r,q,x1;
1337:
1338: if (!sy) err(dvmer1);
1339: if (!sx)
1340: {
1341: if (!z || z == ONLY_REM) return gzero;
1342: *z=gzero; return gzero;
1343: }
1344: lx=lgefint(x);
1345: ly=lgefint(y); lz=lx-ly;
1346: if (lz <= 0)
1347: {
1348: if (lz == 0)
1349: {
1350: for (i=2; i<lx; i++)
1351: if (x[i] != y[i])
1352: {
1353: if ((ulong)x[i] > (ulong)y[i]) goto DIVIDE;
1354: goto TRIVIAL;
1355: }
1356: if (z == ONLY_REM) return gzero;
1357: if (z) *z = gzero;
1358: if (sx < 0) sy = -sy;
1359: return stoi(sy);
1360: }
1361: TRIVIAL:
1362: if (z == ONLY_REM) return icopy(x);
1363: if (z) *z = icopy(x);
1364: return gzero;
1365: }
1366: DIVIDE: /* quotient is non-zero */
1367: av=avma; if (sx<0) sy = -sy;
1368: if (ly==3)
1369: {
1370: LOCAL_HIREMAINDER;
1371: si = y[2];
1372: if (si <= (ulong)x[2]) hiremainder=0;
1373: else
1374: {
1375: hiremainder = x[2]; lx--; x++;
1376: }
1377: q = new_chunk(lx); for (i=2; i<lx; i++) q[i]=divll(x[i],si);
1378: if (z == ONLY_REM)
1379: {
1380: avma=av; if (!hiremainder) return gzero;
1381: r=cgeti(3);
1382: r[1] = evalsigne(sx) | evallgefint(3);
1383: r[2]=hiremainder; return r;
1384: }
1385: q[1] = evalsigne(sy) | evallgefint(lx);
1386: q[0] = evaltyp(t_INT) | evallg(lx);
1387: if (!z) return q;
1388: if (!hiremainder) { *z=gzero; return q; }
1389: r=cgeti(3);
1390: r[1] = evalsigne(sx) | evallgefint(3);
1391: r[2] = hiremainder; *z=r; return q;
1392: }
1393:
1394: x1 = new_chunk(lx); sh = bfffo(y[2]);
1395: if (sh)
1396: { /* normalize so that highbit(y) = 1 (shift left x and y by sh bits)*/
1397: register const ulong m = BITS_IN_LONG - sh;
1398: r = new_chunk(ly);
1399: shift_left2(r, y,2,ly-1, 0,sh,m); y = r;
1400: shift_left2(x1,x,2,lx-1, 0,sh,m);
1401: x1[1] = ((ulong)x[2]) >> m;
1402: }
1403: else
1404: {
1405: x1[1]=0; for (j=2; j<lx; j++) x1[j]=x[j];
1406: }
1407: x=x1; si=y[2]; saux=y[3];
1408: for (i=0; i<=lz; i++)
1409: { /* x1 = x + i */
1410: LOCAL_HIREMAINDER;
1411: LOCAL_OVERFLOW;
1412: if ((ulong)x1[1] == si)
1413: {
1414: qp = MAXULONG; k=addll(si,x1[2]);
1415: }
1416: else
1417: {
1418: hiremainder=x1[1]; overflow=0;
1419: qp=divll(x1[2],si); k=hiremainder;
1420: }
1421: if (!overflow)
1422: {
1423: long k3 = subll(mulll(qp,saux), x1[3]);
1424: long k4 = subllx(hiremainder,k);
1425: while (!overflow && k4) { qp--; k3=subll(k3,saux); k4=subllx(k4,si); }
1426: }
1427: hiremainder=0;
1428: for (j=ly-1; j>1; j--)
1429: {
1430: x1[j] = subll(x1[j], addmul(qp,y[j]));
1431: hiremainder+=overflow;
1432: }
1433: if ((ulong)x1[1] < hiremainder)
1434: {
1435: overflow=0; qp--;
1436: for (j=ly-1; j>1; j--) x1[j] = addllx(x1[j],y[j]);
1437: }
1438: x1[1]=qp; x1++;
1439: }
1440:
1441: lq = lz+2;
1442: if (!z)
1443: {
1444: qd = (ulong*)av;
1445: xd = (ulong*)(x + lq);
1446: if (x[1]) { lz++; lq++; }
1447: while (lz--) *--qd = *--xd;
1448: *--qd = evalsigne(sy) | evallgefint(lq);
1449: *--qd = evaltyp(t_INT) | evallg(lq);
1450: avma = (long)qd; return (GEN)qd;
1451: }
1452:
1453: j=lq; while (j<lx && !x[j]) j++;
1454: lz = lx-j;
1455: if (z == ONLY_REM)
1456: {
1457: if (lz==0) { avma = av; return gzero; }
1458: rd = (ulong*)av; lr = lz+2;
1459: xd = (ulong*)(x + lx);
1460: if (!sh) while (lz--) *--rd = *--xd;
1461: else
1462: { /* shift remainder right by sh bits */
1463: const ulong shl = BITS_IN_LONG - sh;
1464: ulong l;
1465: xd--;
1466: while (--lz) /* fill r[3..] */
1467: {
1468: l = *xd >> sh;
1469: *--rd = l | (*--xd << shl);
1470: }
1471: l = *xd >> sh;
1472: if (l) *--rd = l; else lr--;
1473: }
1474: *--rd = evalsigne(sx) | evallgefint(lr);
1475: *--rd = evaltyp(t_INT) | evallg(lr);
1476: avma = (long)rd; return (GEN)rd;
1477: }
1478:
1479: lr = lz+2;
1480: rd = NULL; /* gcc -Wall */
1481: if (lz)
1482: { /* non zero remainder: initialize rd */
1483: xd = (ulong*)(x + lx);
1484: if (!sh)
1485: {
1486: rd = (ulong*)avma; (void)new_chunk(lr);
1487: while (lz--) *--rd = *--xd;
1488: }
1489: else
1490: { /* shift remainder right by sh bits */
1491: const ulong shl = BITS_IN_LONG - sh;
1492: ulong l;
1493: rd = (ulong*)x; /* overwrite shifted y */
1494: xd--;
1495: while (--lz)
1496: {
1497: l = *xd >> sh;
1498: *--rd = l | (*--xd << shl);
1499: }
1500: l = *xd >> sh;
1501: if (l) *--rd = l; else lr--;
1502: }
1503: *--rd = evalsigne(sx) | evallgefint(lr);
1504: *--rd = evaltyp(t_INT) | evallg(lr);
1505: rd += lr;
1506: }
1507: qd = (ulong*)av;
1508: xd = (ulong*)(x + lq);
1509: if (x[1]) lq++;
1510: j = lq-2; while (j--) *--qd = *--xd;
1511: *--qd = evalsigne(sy) | evallgefint(lq);
1512: *--qd = evaltyp(t_INT) | evallg(lq);
1513: q = (GEN)qd;
1514: if (lr==2) *z = gzero;
1515: else
1516: { /* rd has been properly initialized: we had lz > 0 */
1517: while (lr--) *--qd = *--rd;
1518: *z = (GEN)qd;
1519: }
1520: avma = (long)qd; return q;
1521: }
1522:
1523: /* assume y > x > 0. return y mod x */
1524: ulong
1525: mppgcd_resiu(GEN y, ulong x)
1526: {
1527: long i, ly = lgefint(y);
1528: LOCAL_HIREMAINDER;
1529:
1530: hiremainder = 0;
1531: for (i=2; i<ly; i++) (void)divll(y[i],x);
1532: return hiremainder;
1533: }
1534:
1535: /* Assume x>y>0, both of them odd. return x-y if x=y mod 4, x+y otherwise */
1536: void
1537: mppgcd_plus_minus(GEN x, GEN y, GEN res)
1538: {
1539: long av = avma;
1540: long lx = lgefint(x)-1;
1541: long ly = lgefint(y)-1, lt,m,i;
1542: GEN t;
1543:
1544: if ((x[lx]^y[ly]) & 3) /* x != y mod 4*/
1545: t = addiispec(x+2,y+2,lx-1,ly-1);
1546: else
1547: t = subiispec(x+2,y+2,lx-1,ly-1);
1548:
1549: lt = lgefint(t)-1; while (!t[lt]) lt--;
1550: m = vals(t[lt]); lt++;
1551: if (m == 0) /* 2^32 | t */
1552: {
1553: for (i = 2; i < lt; i++) res[i] = t[i];
1554: }
1555: else if (t[2] >> m)
1556: {
1557: shift_right(res,t, 2,lt, 0,m);
1558: }
1559: else
1560: {
1561: lt--; t++;
1562: shift_right(res,t, 2,lt, t[1],m);
1563: }
1564: res[1] = evalsigne(1)|evallgefint(lt);
1565: avma = av;
1566: }
1567:
1568: /* private version of mulss */
1569: ulong
1570: smulss(ulong x, ulong y, ulong *rem)
1571: {
1572: LOCAL_HIREMAINDER;
1573: x=mulll(x,y); *rem=hiremainder; return x;
1574: }
1575:
1576: #ifdef LONG_IS_64BIT
1577: # define DIVCONVI 7
1578: #else
1579: # define DIVCONVI 14
1580: #endif
1581:
1582: /* convert integer --> base 10^9 [assume x > 0] */
1583: GEN
1584: convi(GEN x)
1585: {
1586: ulong av=avma, lim;
1587: long lz;
1588: GEN z,p1;
1589:
1590: lz = 3 + ((lgefint(x)-2)*15)/DIVCONVI;
1591: z=new_chunk(lz); z[1] = -1; p1 = z+2;
1592: lim = stack_lim(av,1);
1593: for(;;)
1594: {
1595: x = divis(x,1000000000); *p1++ = hiremainder;
1596: if (!signe(x)) { avma=av; return p1; }
1597: if (low_stack(lim, stack_lim(av,1))) x = gerepileuptoint((long)z,x);
1598: }
1599: }
1600: #undef DIVCONVI
1601:
1602: /* convert fractional part --> base 10^9 [assume expo(x) < 0] */
1603: GEN
1604: confrac(GEN x)
1605: {
1606: long lx=lg(x), ex = -expo(x)-1, d,m,ly,ey,lr,i,j;
1607: GEN y,y1;
1608:
1609: ey = ex + bit_accuracy(lx);
1610: ly = 1 + ((ey-1) >> TWOPOTBITS_IN_LONG);
1611: d = ex >> TWOPOTBITS_IN_LONG;
1612: m = ex & (BITS_IN_LONG-1);
1613: y = new_chunk(ly); y1 = y + (d-2);
1614: while (d--) y[d]=0;
1615: if (!m)
1616: for (i=2; i<lx; i++) y1[i] = x[i];
1617: else
1618: { /* shift x left sh bits */
1619: const ulong sh=BITS_IN_LONG-m;
1620: ulong k=0, l;
1621: for (i=2; i<lx; i++)
1622: {
1623: l = x[i];
1624: y1[i] = (l>>m)|k;
1625: k = l<<sh;
1626: }
1627: y1[i] = k;
1628: }
1629: lr = 1 + ((long) (ey*L2SL10))/9;
1630: y1 = new_chunk(lr);
1631: for (j=0; j<lr; j++)
1632: {
1633: LOCAL_HIREMAINDER;
1634: hiremainder=0;
1635: for (i=ly-1; i>=0; i--) y[i]=addmul(y[i],1000000000);
1636: y1[j]=hiremainder;
1637: }
1638: return y1;
1639: }
1640:
1641: GEN
1642: truedvmdii(GEN x, GEN y, GEN *z)
1643: {
1644: long av = avma;
1645: GEN r, q = dvmdii(x,y,&r); /* assume that r is last on stack */
1646: GEN *gptr[2];
1647:
1648: if (signe(r)>=0)
1649: {
1650: if (z == ONLY_REM) return gerepileuptoint(av,r);
1651: if (z) *z = r; else cgiv(r);
1652: return q;
1653: }
1654:
1655: if (z == ONLY_REM)
1656: {
1657: r = subiispec(y+2,r+2, lgefint(y)-2,lgefint(r)-2);
1658: return gerepileuptoint(av, r);
1659: }
1660: q = addsi(-signe(y),q);
1661: if (!z) return gerepileuptoint(av, q);
1662:
1663: *z = subiispec(y+2,r+2, lgefint(y)-2,lgefint(r)-2);
1664: gptr[0]=&q; gptr[1]=z; gerepilemanysp(av,(long)r,gptr,2);
1665: return q;
1666: }
1667:
1668: /* EXACT INTEGER DIVISION */
1669:
1670: /* Find c such that 1=c*b mod 2^BITS_IN_LONG, assuming b odd (unchecked) */
1671: static ulong
1672: invrev(ulong b)
1673: {
1674: ulong x;
1675: switch(b & 7)
1676: {
1677: case 3:
1678: case 5: x = b+8; break;
1679: default: x = b; break;
1680: } /* x = b^(-1) mod 2^4 */
1681: x = x*(2-b*x);
1682: x = x*(2-b*x);
1683: x = x*(2-b*x); /* b^(-1) mod 2^32 (Newton applied to 1/x - b = 0) */
1684: #ifdef LONG_IS_64BIT
1685: x = x*(2-b*x); /* b^(-1) mod 2^64 */
1686: #endif
1687: return x;
1688: }
1689:
1690: /* assume xy>0, y odd */
1691: GEN
1692: diviuexact(GEN x, ulong y)
1693: {
1694: long i,lz,lx;
1695: ulong q, yinv;
1696: GEN z, z0, x0, x0min;
1697:
1698: if (y == 1) return icopy(x);
1699: lx = lgefint(x);
1700: if (lx == 3) return stoi((ulong)x[2] / y);
1701: yinv = invrev(y);
1702: lz = (y <= (ulong)x[2]) ? lx : lx-1;
1703: z = new_chunk(lz);
1704: z0 = z + lz;
1705: x0 = x + lx; x0min = x + lx-lz+2;
1706:
1707: while (x0 > x0min)
1708: {
1709: *--z0 = q = yinv*((ulong)*--x0); /* i-th quotient */
1710: if (!q) continue;
1711: /* x := x - q * y */
1712: { /* update neither lowest word (could set it to 0) nor highest ones */
1713: register GEN x1 = x0 - 1;
1714: LOCAL_HIREMAINDER;
1715: (void)mulll(q,y);
1716: if (hiremainder)
1717: {
1718: if ((ulong)*x1 < hiremainder)
1719: {
1720: *x1 -= hiremainder;
1721: do (*--x1)--; while ((ulong)*x1 == MAXULONG);
1722: }
1723: else
1724: *x1 -= hiremainder;
1725: }
1726: }
1727: }
1728: i=2; while(!z[i]) i++;
1729: z += i-2; lz -= i-2;
1730: z[0] = evaltyp(t_INT)|evallg(lz);
1731: z[1] = evalsigne(1)|evallg(lz);
1732: avma = (ulong)z; return z;
1733: }
1734:
1735: /* Find z such that x=y*z, knowing that y | x (unchecked)
1736: * Method: y0 z0 = x0 mod B = 2^BITS_IN_LONG ==> z0 = 1/y0 mod B.
1737: * Set x := (x - z0 y) / B, updating only relevant words, and repeat */
1738: GEN
1739: diviiexact(GEN x, GEN y)
1740: {
1741: long av,lx,ly,lz,vy,i,ii,sx = signe(x), sy = signe(y);
1742: ulong y0inv,q;
1743: GEN z;
1744:
1745: if (!sy) err(dvmer1);
1746: if (!sx) return gzero;
1747: vy = vali(y); av = avma;
1748: (void)new_chunk(lgefint(x)); /* enough room for z */
1749: if (vy)
1750: { /* make y odd */
1751: #if 0
1752: if (vali(x) < vy) err(talker,"impossible division in diviiexact");
1753: #endif
1754: y = shifti(y,-vy);
1755: x = shifti(x,-vy);
1756: }
1757: else x = icopy(x); /* necessary because we destroy x */
1758: avma = av; /* will erase our x,y when exiting */
1759: /* now y is odd */
1760: ly = lgefint(y);
1761: if (ly == 3)
1762: {
1763: x = diviuexact(x,(ulong)y[2]);
1764: if (signe(x)) setsigne(x,sx*sy); /* should have x != 0 at this point */
1765: return x;
1766: }
1767: lx = lgefint(x); if (ly>lx) err(talker,"impossible division in diviiexact");
1768: y0inv = invrev(y[ly-1]);
1769: i=2; while (i<ly && y[i]==x[i]) i++;
1770: lz = (i==ly || (ulong)y[i] < (ulong)x[i]) ? lx-ly+3 : lx-ly+2;
1771: z = new_chunk(lz);
1772:
1773: y += ly - 1; /* now y[-i] = i-th word of y */
1774: for (ii=lx-1,i=lz-1; i>=2; i--,ii--)
1775: {
1776: long limj;
1777: LOCAL_HIREMAINDER;
1778: LOCAL_OVERFLOW;
1779:
1780: z[i] = q = y0inv*((ulong)x[ii]); /* i-th quotient */
1781: if (!q) continue;
1782:
1783: /* x := x - q * y */
1784: (void)mulll(q,y[0]); limj = max(lx - lz, ii+3-ly);
1785: { /* update neither lowest word (could set it to 0) nor highest ones */
1786: register GEN x0 = x + (ii - 1), y0 = y - 1, xlim = x + limj;
1787: for (; x0 >= xlim; x0--, y0--)
1788: {
1789: *x0 = subll(*x0, addmul(q,*y0));
1790: hiremainder += overflow;
1791: }
1792: if (hiremainder && limj != lx - lz)
1793: {
1794: if ((ulong)*x0 < hiremainder)
1795: {
1796: *x0 -= hiremainder;
1797: do (*--x0)--; while ((ulong)*x0 == MAXULONG);
1798: }
1799: else
1800: *x0 -= hiremainder;
1801: }
1802: }
1803: }
1804: #if 0
1805: i=2; while(i<lz && !z[i]) i++;
1806: if (i==lz) err(bugparier,"diviiexact");
1807: #else
1808: i=2; while(!z[i]) i++;
1809: #endif
1810: z += i-2; lz -= (i-2);
1811: z[0] = evaltyp(t_INT)|evallg(lz);
1812: z[1] = evalsigne(sx*sy)|evallg(lz);
1813: avma = (ulong)z; return z;
1814: }
1815:
1816: long
1817: smodsi(long x, GEN y)
1818: {
1819: if (x<0) err(talker,"negative small integer in smodsi");
1820: (void)divsi(x,y); return hiremainder;
1821: }
1822:
1823: /* x and y are integers. Return 1 if |x| == |y|, 0 otherwise */
1824: int
1825: absi_equal(GEN x, GEN y)
1826: {
1827: long lx,i;
1828:
1829: if (!signe(x)) return !signe(y);
1830: if (!signe(y)) return 0;
1831:
1832: lx=lgefint(x); if (lx != lgefint(y)) return 0;
1833: i=2; while (i<lx && x[i]==y[i]) i++;
1834: return (i==lx);
1835: }
1836:
1837: /* x and y are integers. Return sign(|x| - |y|) */
1838: int
1839: absi_cmp(GEN x, GEN y)
1840: {
1841: long lx,ly,i;
1842:
1843: if (!signe(x)) return signe(y)? -1: 0;
1844: if (!signe(y)) return 1;
1845:
1846: lx=lgefint(x); ly=lgefint(y);
1847: if (lx>ly) return 1;
1848: if (lx<ly) return -1;
1849: i=2; while (i<lx && x[i]==y[i]) i++;
1850: if (i==lx) return 0;
1851: return ((ulong)x[i] > (ulong)y[i])? 1: -1;
1852: }
1853:
1854: /* x and y are reals. Return sign(|x| - |y|) */
1855: int
1856: absr_cmp(GEN x, GEN y)
1857: {
1858: long ex,ey,lx,ly,lz,i;
1859:
1860: if (!signe(x)) return signe(y)? -1: 0;
1861: if (!signe(y)) return 1;
1862:
1863: ex=expo(x); ey=expo(y);
1864: if (ex>ey) return 1;
1865: if (ex<ey) return -1;
1866:
1867: lx=lg(x); ly=lg(y); lz = (lx<ly)?lx:ly;
1868: i=2; while (i<lz && x[i]==y[i]) i++;
1869: if (i<lz) return ((ulong)x[i] > (ulong)y[i])? 1: -1;
1870: if (lx>=ly)
1871: {
1872: while (i<lx && !x[i]) i++;
1873: return (i==lx)? 0: 1;
1874: }
1875: while (i<ly && !y[i]) i++;
1876: return (i==ly)? 0: -1;
1877: }
1878:
1879: /********************************************************************/
1880: /** **/
1881: /** INTEGER MULTIPLICATION (KARATSUBA) **/
1882: /** **/
1883: /********************************************************************/
1884: #define _sqri_l 47
1885: #define _muli_l 25 /* optimal on PII 350MHz + gcc 2.7.2.1 (gp-dyn) */
1886:
1887: #if 0 /* for tunings */
1888: long KARATSUBA_SQRI_LIMIT = _sqri_l;
1889: long KARATSUBA_MULI_LIMIT = _muli_l;
1890:
1891: void setsqri(long a) { KARATSUBA_SQRI_LIMIT = a; }
1892: void setmuli(long a) { KARATSUBA_MULI_LIMIT = a; }
1893:
1894: GEN
1895: speci(GEN x, long nx)
1896: {
1897: GEN z;
1898: long i;
1899: if (!nx) return gzero;
1900: z = cgeti(nx+2); z[1] = evalsigne(1)|evallgefint(nx+2);
1901: for (i=0; i<nx; i++) z[i+2] = x[i];
1902: return z;
1903: }
1904: #else
1905: # define KARATSUBA_SQRI_LIMIT _sqri_l
1906: # define KARATSUBA_MULI_LIMIT _muli_l
1907: #endif
1908:
1909: /* nx >= ny = num. of digits of x, y (not GEN, see mulii) */
1910: #ifdef INLINE
1911: INLINE
1912: #endif
1913: GEN
1914: muliispec(GEN x, GEN y, long nx, long ny)
1915: {
1916: GEN z2e,z2d,yd,xd,ye,zd;
1917: long p1,lz;
1918: LOCAL_HIREMAINDER;
1919:
1920: if (!ny) return gzero;
1921: zd = (GEN)avma; lz = nx+ny+2;
1922: (void)new_chunk(lz);
1923: xd = x + nx;
1924: yd = y + ny;
1925: ye = yd; p1 = *--xd;
1926:
1927: *--zd = mulll(p1, *--yd); z2e = zd;
1928: while (yd > y) *--zd = addmul(p1, *--yd);
1929: *--zd = hiremainder;
1930:
1931: while (xd > x)
1932: {
1933: LOCAL_OVERFLOW;
1934: yd = ye; p1 = *--xd;
1935:
1936: z2d = --z2e;
1937: *z2d = addll(mulll(p1, *--yd), *z2d); z2d--;
1938: while (yd > y)
1939: {
1940: hiremainder += overflow;
1941: *z2d = addll(addmul(p1, *--yd), *z2d); z2d--;
1942: }
1943: *--zd = hiremainder + overflow;
1944: }
1945: if (*zd == 0) { zd++; lz--; } /* normalize */
1946: *--zd = evalsigne(1) | evallgefint(lz);
1947: *--zd = evaltyp(t_INT) | evallg(lz);
1948: avma=(long)zd; return zd;
1949: }
1950:
1951: #ifdef INLINE
1952: INLINE
1953: #endif
1954: GEN
1955: sqrispec(GEN x, long nx)
1956: {
1957: GEN z2e,z2d,yd,xd,zd,x0,z0;
1958: long p1,lz;
1959: LOCAL_HIREMAINDER;
1960:
1961: if (!nx) return gzero;
1962: zd = (GEN)avma; lz = (nx+1) << 1;
1963: z0 = new_chunk(lz);
1964: if (nx == 1)
1965: {
1966: *--zd = mulll(*x, *x);
1967: *--zd = hiremainder; goto END;
1968: }
1969: xd = x + nx;
1970:
1971: /* compute double products --> zd */
1972: p1 = *--xd; yd = xd; --zd;
1973: *--zd = mulll(p1, *--yd); z2e = zd;
1974: while (yd > x) *--zd = addmul(p1, *--yd);
1975: *--zd = hiremainder;
1976:
1977: x0 = x+1;
1978: while (xd > x0)
1979: {
1980: LOCAL_OVERFLOW;
1981: p1 = *--xd; yd = xd;
1982:
1983: z2e -= 2; z2d = z2e;
1984: *z2d = addll(mulll(p1, *--yd), *z2d); z2d--;
1985: while (yd > x)
1986: {
1987: hiremainder += overflow;
1988: *z2d = addll(addmul(p1, *--yd), *z2d); z2d--;
1989: }
1990: *--zd = hiremainder + overflow;
1991: }
1992: /* multiply zd by 2 (put result in zd - 1) */
1993: zd[-1] = ((*zd & HIGHBIT) != 0);
1994: shift_left(zd, zd, 0, (nx<<1)-3, 0, 1);
1995:
1996: /* add the squares */
1997: xd = x + nx; zd = z0 + lz;
1998: p1 = *--xd;
1999: zd--; *zd = mulll(p1,p1);
2000: zd--; *zd = addll(hiremainder, *zd);
2001: while (xd > x)
2002: {
2003: p1 = *--xd;
2004: zd--; *zd = addll(mulll(p1,p1)+ overflow, *zd);
2005: zd--; *zd = addll(hiremainder + overflow, *zd);
2006: }
2007:
2008: END:
2009: if (*zd == 0) { zd++; lz--; } /* normalize */
2010: *--zd = evalsigne(1) | evallgefint(lz);
2011: *--zd = evaltyp(t_INT) | evallg(lz);
2012: avma=(long)zd; return zd;
2013: }
2014:
2015: /* return (x shifted left d words) + y. Assume d > 0, x > 0 and y >= 0 */
2016: static GEN
2017: addshiftw(GEN x, GEN y, long d)
2018: {
2019: GEN z,z0,y0,yd, zd = (GEN)avma;
2020: long a,lz,ly = lgefint(y);
2021:
2022: z0 = new_chunk(d);
2023: a = ly-2; yd = y+ly;
2024: if (a >= d)
2025: {
2026: y0 = yd-d; while (yd > y0) *--zd = *--yd; /* copy last d words of y */
2027: a -= d;
2028: if (a)
2029: z = addiispec(x+2, y+2, lgefint(x)-2, a);
2030: else
2031: z = icopy(x);
2032: }
2033: else
2034: {
2035: y0 = yd-a; while (yd > y0) *--zd = *--yd; /* copy last a words of y */
2036: while (zd >= z0) *--zd = 0; /* complete with 0s */
2037: z = icopy(x);
2038: }
2039: lz = lgefint(z)+d;
2040: z[1] = evalsigne(1) | evallgefint(lz);
2041: z[0] = evaltyp(t_INT) | evallg(lz); return z;
2042: }
2043:
2044: /* Fast product (Karatsuba) of integers. a and b are "special" GENs
2045: * c,c0,c1,c2 are genuine GENs.
2046: */
2047: static GEN
2048: quickmulii(GEN a, GEN b, long na, long nb)
2049: {
2050: GEN a0,c,c0;
2051: long av,n0,n0a,i;
2052:
2053: if (na < nb) swapspec(a,b, na,nb);
2054: if (nb == 1) return mulsispec(*b, a, na);
2055: if (nb == 0) return gzero;
2056: if (nb < KARATSUBA_MULI_LIMIT) return muliispec(a,b,na,nb);
2057: i=(na>>1); n0=na-i; na=i;
2058: av=avma; a0=a+na; n0a=n0;
2059: while (!*a0 && n0a) { a0++; n0a--; }
2060:
2061: if (n0a && nb > n0)
2062: { /* nb <= na <= n0 */
2063: GEN b0,c1,c2;
2064: long n0b;
2065:
2066: nb -= n0;
2067: c = quickmulii(a,b,na,nb);
2068: b0 = b+nb; n0b = n0;
2069: while (!*b0 && n0b) { b0++; n0b--; }
2070: if (n0b)
2071: {
2072: c0 = quickmulii(a0,b0, n0a,n0b);
2073:
2074: c2 = addiispec(a0,a, n0a,na);
2075: c1 = addiispec(b0,b, n0b,nb);
2076: c1 = quickmulii(c1+2,c2+2, lgefint(c1)-2,lgefint(c2)-2);
2077: c2 = addiispec(c0+2, c+2, lgefint(c0)-2,lgefint(c) -2);
2078:
2079: c1 = subiispec(c1+2,c2+2, lgefint(c1)-2,lgefint(c2)-2);
2080: }
2081: else
2082: {
2083: c0 = gzero;
2084: c1 = quickmulii(a0,b, n0a,nb);
2085: }
2086: c = addshiftw(c,c1, n0);
2087: }
2088: else
2089: {
2090: c = quickmulii(a,b,na,nb);
2091: c0 = quickmulii(a0,b,n0a,nb);
2092: }
2093: return gerepileuptoint(av, addshiftw(c,c0, n0));
2094: }
2095:
2096: /* actual operations will take place on a+2 and b+2: we strip the codewords */
2097: GEN
2098: mulii(GEN a,GEN b)
2099: {
2100: long sa,sb;
2101: GEN z;
2102:
2103: sa=signe(a); if (!sa) return gzero;
2104: sb=signe(b); if (!sb) return gzero;
2105: if (sb<0) sa = -sa;
2106: z = quickmulii(a+2,b+2, lgefint(a)-2,lgefint(b)-2);
2107: setsigne(z,sa); return z;
2108: }
2109:
2110: GEN
2111: resiimul(GEN x, GEN sy)
2112: {
2113: GEN r, q, y = (GEN)sy[1], invy;
2114: long av = avma, k;
2115:
2116: k = cmpii(x, y);
2117: if (k <= 0) return k? icopy(x): gzero;
2118: invy = (GEN)sy[2];
2119: q = mulir(x,invy);
2120: q = mptrunc(q); /* <= divii(x, y) (at most 1 less) */
2121: r = subii(x, mulii(y,q));
2122: /* resii(x,y) + y >= r >= resii(x,y) */
2123: k = cmpii(r, y);
2124: if (k >= 0)
2125: {
2126: if (k == 0) { avma = av; return gzero; }
2127: r = subiispec(r+2, y+2, lgefint(r)-2, lgefint(y)-2);
2128: }
2129: #if 0
2130: q = subii(r,resii(x,y));
2131: if (signe(q))
2132: err(talker,"bug in resiimul: x = %Z\ny = %Z\ndif = %Z", x,y,q);
2133: #endif
2134: return gerepileuptoint(av, r); /* = resii(x, y) */
2135: }
2136:
2137: /* x % (2^n), assuming x, n >= 0 */
2138: GEN
2139: resmod2n(GEN x, long n)
2140: {
2141: long hi,l,k,lx,ly;
2142: GEN z, xd, zd;
2143:
2144: if (!signe(x) || !n) return gzero;
2145:
2146: l = n & (BITS_IN_LONG-1); /* n % BITS_IN_LONG */
2147: k = n >> TWOPOTBITS_IN_LONG; /* n / BITS_IN_LONG */
2148: lx = lgefint(x);
2149: if (lx < k+3) return icopy(x);
2150:
2151: xd = x + (lx-k-1);
2152: /* x = |_|...|#|1|...|k| : copy the last l bits of # and the last k words
2153: * ^--- initial xd */
2154: hi = *xd & ((1<<l) - 1); /* last l bits of # = top bits of result */
2155: if (!hi)
2156: { /* strip leading zeroes from result */
2157: xd++; while (k && !*xd) { k--; xd++; }
2158: if (!k) return gzero;
2159: ly = k+2; xd--;
2160: }
2161: else
2162: ly = k+3;
2163:
2164: zd = z = cgeti(ly);
2165: *++zd = evalsigne(1) | evallgefint(ly);
2166: if (hi) *++zd = hi;
2167: for ( ;k; k--) *++zd = *++xd;
2168: return z;
2169: }
2170:
2171: static GEN
2172: quicksqri(GEN a, long na)
2173: {
2174: GEN a0,c;
2175: long av,n0,n0a,i;
2176:
2177: if (na < KARATSUBA_SQRI_LIMIT) return sqrispec(a,na);
2178: i=(na>>1); n0=na-i; na=i;
2179: av=avma; a0=a+na; n0a=n0;
2180: while (!*a0 && n0a) { a0++; n0a--; }
2181: c = quicksqri(a,na);
2182: if (n0a)
2183: {
2184: GEN t, c1, c0 = quicksqri(a0,n0a);
2185: #if 0
2186: c1 = shifti(quickmulii(a0,a, n0a,na),1);
2187: #else /* slower !!! */
2188: t = addiispec(a0,a,n0a,na);
2189: t = quicksqri(t+2,lgefint(t)-2);
2190: c1= addiispec(c0+2,c+2, lgefint(c0)-2, lgefint(c)-2);
2191: c1= subiispec(t+2, c1+2, lgefint(t)-2, lgefint(c1)-2);
2192: #endif
2193: c = addshiftw(c,c1, n0);
2194: c = addshiftw(c,c0, n0);
2195: }
2196: else
2197: c = addshiftw(c,gzero,n0<<1);
2198: return gerepileuptoint(av, c);
2199: }
2200:
2201: GEN
2202: sqri(GEN a) { return quicksqri(a+2, lgefint(a)-2); }
2203:
2204: #define MULRR_LIMIT 32
2205: #define MULRR2_LIMIT 32
2206:
2207: #if 0
2208: GEN
2209: karamulrr1(GEN x, GEN y)
2210: {
2211: long sx,sy;
2212: long i,i1,i2,lx=lg(x),ly=lg(y),e,flag,garde;
2213: long lz2,lz3,lz4;
2214: GEN z,lo1,lo2,hi;
2215:
2216: /* ensure that lg(y) >= lg(x) */
2217: if (lx>ly) { lx=ly; z=x; x=y; y=z; flag=1; } else flag = (lx!=ly);
2218: if (lx < MULRR_LIMIT) return mulrr(x,y);
2219: sx=signe(x); sy=signe(y);
2220: e = evalexpo(expo(x)+expo(y));
2221: if (!sx || !sy) { z=cgetr(3); z[2]=0; z[1]=e; return z; }
2222: if (sy<0) sx = -sx;
2223: ly=lx+flag; z=cgetr(lx);
2224: lz2 = (lx>>1); lz3 = lx-lz2;
2225: x += 2; lx -= 2;
2226: y += 2; ly -= 2;
2227: hi = quickmulii(x,y,lz2,lz2);
2228: i1=lz2; while (i1<lx && !x[i1]) i1++;
2229: lo1 = quickmulii(y,x+i1,lz2,lx-i1);
2230: i2=lz2; while (i2<ly && !y[i2]) i2++;
2231: lo2 = quickmulii(x,y+i2,lz2,ly-i2);
2232:
2233: if (signe(lo1))
2234: {
2235: if (flag) { lo2 = addshiftw(lo1,lo2,1); lz3++; } else lo2=addii(lo1,lo2);
2236: }
2237: lz4=lgefint(lo2)-lz3;
2238: if (lz4>0) hi = addiispec(hi+2,lo2+2, lgefint(hi)-2,lz4);
2239: if (hi[2] < 0)
2240: {
2241: e++; garde=hi[lx+2];
2242: for (i=lx+1; i>=2 ; i--) z[i]=hi[i];
2243: }
2244: else
2245: {
2246: garde = (hi[lx+2] << 1);
2247: shift_left(z,hi,2,lx+1, garde, 1);
2248: }
2249: z[1]=evalsigne(sx) | e;
2250: if (garde < 0)
2251: { /* round to nearest */
2252: i=lx+2; do z[--i]++; while (z[i]==0);
2253: if (i==1) z[2]=HIGHBIT;
2254: }
2255: avma=(long)z; return z;
2256: }
2257:
2258: GEN
2259: karamulrr2(GEN x, GEN y)
2260: {
2261: long sx,sy,i,lx=lg(x),ly=lg(y),e,flag,garde;
2262: GEN z,hi;
2263:
2264: if (lx>ly) { lx=ly; z=x; x=y; y=z; flag=1; } else flag = (lx!=ly);
2265: if (lx < MULRR2_LIMIT) return mulrr(x,y);
2266: ly=lx+flag; sx=signe(x); sy=signe(y);
2267: e = evalexpo(expo(x)+expo(y));
2268: if (!sx || !sy) { z=cgetr(3); z[2]=0; z[1]=e; return z; }
2269: if (sy<0) sx = -sx;
2270: z=cgetr(lx);
2271: hi=quickmulii(y+2,x+2,ly-2,lx-2);
2272: if (hi[2] < 0)
2273: {
2274: e++; garde=hi[lx];
2275: for (i=2; i<lx ; i++) z[i]=hi[i];
2276: }
2277: else
2278: {
2279: garde = (hi[lx] << 1);
2280: shift_left(z,hi,2,lx-1, garde, 1);
2281: }
2282: z[1]=evalsigne(sx) | e;
2283: if (garde < 0)
2284: { /* round to nearest */
2285: i=lx; do z[--i]++; while (z[i]==0);
2286: if (i==1) z[2]=HIGHBIT;
2287: }
2288: avma=(long)z; return z;
2289: }
2290:
2291: GEN
2292: karamulrr(GEN x, GEN y, long flag)
2293: {
2294: switch(flag)
2295: {
2296: case 1: return karamulrr1(x,y);
2297: case 2: return karamulrr2(x,y);
2298: default: err(flagerr,"karamulrr");
2299: }
2300: return NULL; /* not reached */
2301: }
2302:
2303: GEN
2304: karamulir(GEN x, GEN y, long flag)
2305: {
2306: long sx=signe(x),lz,i;
2307: GEN z,temp,z1;
2308:
2309: if (!sx) return gzero;
2310: lz=lg(y); z=cgetr(lz);
2311: temp=cgetr(lz+1); affir(x,temp);
2312: z1=karamulrr(temp,y,flag);
2313: for (i=1; i<lz; i++) z[i]=z1[i];
2314: avma=(long)z; return z;
2315: }
2316: #endif
2317:
2318: #ifdef LONG_IS_64BIT
2319:
2320: #if PARI_BYTE_ORDER == LITTLE_ENDIAN_64 || PARI_BYTE_ORDER == BIG_ENDIAN_64
2321: #else
2322: error... unknown machine
2323: #endif
2324:
2325: GEN
2326: dbltor(double x)
2327: {
2328: GEN z = cgetr(3);
2329: long e;
2330: union { double f; ulong i; } fi;
2331: const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
2332: const int exp_mid = 0x3ff;/* exponent bias */
2333: const int expo_len = 11; /* number of bits of exponent */
2334: LOCAL_HIREMAINDER;
2335:
2336: if (x==0) { z[1]=evalexpo(-308); z[2]=0; return z; }
2337: fi.f = x;
2338: e = ((fi.i & (HIGHBIT-1)) >> mant_len) - exp_mid;
2339: z[1] = evalexpo(e) | evalsigne(x<0? -1: 1);
2340: z[2] = (fi.i << expo_len) | HIGHBIT;
2341: return z;
2342: }
2343:
2344: double
2345: rtodbl(GEN x)
2346: {
2347: long ex,s=signe(x);
2348: ulong a;
2349: union { double f; ulong i; } fi;
2350: const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
2351: const int exp_mid = 0x3ff;/* exponent bias */
2352: const int expo_len = 11; /* number of bits of exponent */
2353: LOCAL_HIREMAINDER;
2354:
2355: if (typ(x)==t_INT && !s) return 0.0;
2356: if (typ(x)!=t_REAL) err(typeer,"rtodbl");
2357: if (!s || (ex=expo(x)) < - exp_mid) return 0.0;
2358:
2359: /* start by rounding to closest */
2360: a = (x[2] & (HIGHBIT-1)) + 0x400;
2361: if (a & HIGHBIT) { ex++; a=0; }
2362: if (ex >= exp_mid) err(rtodber);
2363: fi.i = ((ex + exp_mid) << mant_len) | (a >> expo_len);
2364: if (s<0) fi.i |= HIGHBIT;
2365: return fi.f;
2366: }
2367:
2368: #else
2369:
2370: #if PARI_BYTE_ORDER == LITTLE_ENDIAN
2371: # define INDEX0 1
2372: # define INDEX1 0
2373: #elif PARI_BYTE_ORDER == BIG_ENDIAN
2374: # define INDEX0 0
2375: # define INDEX1 1
2376: #else
2377: error... unknown machine
2378: #endif
2379:
2380: GEN
2381: dbltor(double x)
2382: {
2383: GEN z;
2384: long e;
2385: union { double f; ulong i[2]; } fi;
2386: const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
2387: const int exp_mid = 0x3ff;/* exponent bias */
2388: const int shift = mant_len-32;
2389: const int expo_len = 11; /* number of bits of exponent */
2390:
2391: if (x==0) { z=cgetr(3); z[1]=evalexpo(-308); z[2]=0; return z; }
2392: fi.f = x; z=cgetr(4);
2393: {
2394: const ulong a = fi.i[INDEX0];
2395: const ulong b = fi.i[INDEX1];
2396: e = ((a & (HIGHBIT-1)) >> shift) - exp_mid;
2397: z[1] = evalexpo(e) | evalsigne(x<0? -1: 1);
2398: z[3] = b << expo_len;
2399: z[2] = HIGHBIT | b >> (BITS_IN_LONG-expo_len) | (a << expo_len);
2400: }
2401: return z;
2402: }
2403:
2404: double
2405: rtodbl(GEN x)
2406: {
2407: long ex,s=signe(x),lx=lg(x);
2408: ulong a,b,k;
2409: union { double f; ulong i[2]; } fi;
2410: const int mant_len = 52; /* mantissa bits (excl. hidden bit) */
2411: const int exp_mid = 0x3ff;/* exponent bias */
2412: const int shift = mant_len-32;
2413: const int expo_len = 11; /* number of bits of exponent */
2414:
2415: if (typ(x)==t_INT && !s) return 0.0;
2416: if (typ(x)!=t_REAL) err(typeer,"rtodbl");
2417: if (!s || (ex=expo(x)) < - exp_mid) return 0.0;
2418:
2419: /* start by rounding to closest */
2420: a = x[2] & (HIGHBIT-1);
2421: if (lx > 3)
2422: {
2423: b = x[3] + 0x400UL; if (b < 0x400UL) a++;
2424: if (a & HIGHBIT) { ex++; a=0; }
2425: }
2426: else b = 0;
2427: if (ex > exp_mid) err(rtodber);
2428: ex += exp_mid;
2429: k = (a >> expo_len) | (ex << shift);
2430: if (s<0) k |= HIGHBIT;
2431: fi.i[INDEX0] = k;
2432: fi.i[INDEX1] = (a << (BITS_IN_LONG-expo_len)) | (b >> expo_len);
2433: return fi.f;
2434: }
2435: #endif
2436:
2437: /* Old cgiv without reference count (which was not used anyway)
2438: * Should be a macro.
2439: */
2440: void
2441: cgiv(GEN x)
2442: {
2443: if (x == (GEN) avma)
2444: avma = (long) (x+lg(x));
2445: }
2446:
2447: /********************************************************************/
2448: /** **/
2449: /** INTEGER EXTENDED GCD (AND INVMOD) **/
2450: /** **/
2451: /********************************************************************/
2452:
2453: /* GN 1998Oct25, originally developed in January 1998 under 2.0.4.alpha,
2454: * in the context of trying to improve elliptic curve cryptosystem attacking
2455: * algorithms. 2001Jan02 -- added bezout() functionality.
2456: *
2457: * Two basic ideas - (1) avoid many integer divisions, especially when the
2458: * quotient is 1 (which happens more than 40% of the time). (2) Use Lehmer's
2459: * trick as modified by Jebelean of extracting a couple of words' worth of
2460: * leading bits from both operands, and compute partial quotients from them
2461: * as long as we can be sure of their values. The Jebelean modifications
2462: * consist in reliable inequalities from which we can decide fast whether
2463: * to carry on or to return to the outer loop, and in re-shifting after the
2464: * first word's worth of bits has been used up. All of this is described
2465: * in R. Lercier's these [pp148-153 & 163f.], except his outer loop isn't
2466: * quite right (the catch-up divisions needed when one partial quotient is
2467: * larger than a word are missing).
2468: *
2469: * The API consists of invmod() and bezout() below; the single-word routines
2470: * xgcduu and xxgcduu may be called directly if desired; lgcdii() probably
2471: * doesn't make much sense out of context.
2472: *
2473: * The whole lot is a factor 6 .. 8 faster on word-sized operands, and asym-
2474: * ptotically about a factor 2.5 .. 3, depending on processor architecture,
2475: * than the naive continued-division code. Unfortunately, thanks to the
2476: * unrolled loops and all, the code is a bit lengthy.
2477: */
2478:
2479: /*==================================
2480: * xgcduu(d,d1,f,v,v1,s)
2481: * xxgcduu(d,d1,f,u,u1,v,v1,s)
2482: * rgcduu(d,d1,vmax,u,u1,v,v1,s)
2483: *==================================*/
2484: /*
2485: * Fast `final' extended gcd algorithm, acting on two ulongs. Ideally this
2486: * should be replaced with assembler versions wherever possible. The present
2487: * code essentially does `subtract, compare, and possibly divide' at each step,
2488: * which is reasonable when hardware division (a) exists, (b) is a bit slowish
2489: * and (c) does not depend a lot on the operand values (as on i486). When
2490: * wordsize division is in fact an assembler routine based on subtraction,
2491: * this strategy may not be the most efficient one.
2492: *
2493: * xxgcduu() should be called with d > d1 > 0, returns gcd(d,d1), and assigns
2494: * the usual signless cont.frac. recurrence matrix to [u, u1; v, v1] (i.e.,
2495: * the product of all the [0, 1; 1 q_j] where the leftmost factor arises from
2496: * the quotient of the first division step), and the information about the
2497: * implied signs to s (-1 when an odd number of divisions has been done,
2498: * 1 otherwise). xgcduu() is exactly the same except that u,u1 are not com-
2499: * puted (and not returned, of course).
2500: *
2501: * The input flag f should be set to 1 if we know in advance that gcd(d,d1)==1
2502: * (so we can stop the chain division one step early: as soon as the remainder
2503: * equals 1). Use this when you intend to use only what would be v, or only
2504: * what would be u and v, after that final division step, but not u1 and v1.
2505: * With the flag in force and thus without that final step, the interesting
2506: * quantity/ies will still sit in [u1 and] v1, of course.
2507: *
2508: * For computing the inverse of a single-word INTMOD known to exist, pass f=1
2509: * to xgcduu(), and obtain the result from s and v1. (The routine does the
2510: * right thing when d1==1 already.) For finishing a multiword modinv known
2511: * to exist, pass f=1 to xxgcduu(), and multiply the returned matrix (with
2512: * properly adjusted signs) onto the values v' and v1' previously obtained
2513: * from the multiword division steps. Actually, just take the scalar product
2514: * of [v',v1'] with [u1,-v1], and change the sign if s==-1. (If the final
2515: * step had been carried out, it would be [-u,v], and s would also change.)
2516: * For reducing a rational number to lowest terms, pass f=0 to xgcduu().
2517: * Finally, f=0 with xxgcduu() is useful for Bezout computations.
2518: * [Harrumph. In the above prescription, the sign turns out to be precisely
2519: * wrong.]
2520: * (It is safe for invmod() to call xgcduu() with f=1, because f&1 doesn't
2521: * make a difference when gcd(d,d1)>1. The speedup is negligible.)
2522: *
2523: * In principle, when gcd(d,d1) is known to be 1, it is straightforward to
2524: * recover the final u,u1 given only v,v1 and s. However, it probably isn't
2525: * worthwhile, as it trades a few multiplications for a division.
2526: *
2527: * Note that these routines do not know and do not need to know about the
2528: * PARI stack.
2529: *
2530: * Added 2001Jan15:
2531: * rgcduu() is a variant of xxgcduu() which does not have f (the effect is
2532: * that of f=0), but instead has a ulong vmax parameter, for use in rational
2533: * reconstruction below. It returns when v1 exceeds vmax; v will never
2534: * exceed vmax. (vmax=0 is taken as a synonym of MAXULONG i.e. unlimited,
2535: * in which case rgcduu behaves exactly like xxgcduu with f=0.) The return
2536: * value of rgcduu() is typically meaningless; the interesting part is the
2537: * matrix.
2538: */
2539:
2540: ulong
2541: xgcduu(ulong d, ulong d1, int f, ulong* v, ulong* v1, long *s)
2542: {
2543: ulong xv,xv1, xs, q,res;
2544: LOCAL_HIREMAINDER;
2545:
2546: /* The above blurb contained a lie. The main loop always stops when d1
2547: * has become equal to 1. If (d1 == 1 && !(f&1)) after the loop, we do
2548: * the final `division' of d by 1 `by hand' as it were.
2549: *
2550: * The loop has already been unrolled once. Aggressive optimization could
2551: * well lead to a totally unrolled assembler version...
2552: *
2553: * On modern x86 architectures, this loop is a pig anyway. The division
2554: * instruction always puts its result into the same pair of registers, and
2555: * we always want to use one of them straight away, so pipeline performance
2556: * will suck big time. An assembler version should probably do a first loop
2557: * computing and storing all the quotients -- their number is bounded in
2558: * advance -- and then assembling the matrix in a second pass. On other
2559: * architectures where we can cycle through four or so groups of registers
2560: * and exploit a fast ALU result-to-operand feedback path, this is much less
2561: * of an issue. (Intel sucks. See http://www.x86.org/ ...)
2562: */
2563: xs = res = 0;
2564: xv = 0UL; xv1 = 1UL;
2565: while (d1 > 1UL)
2566: {
2567: d -= d1; /* no need to use subll */
2568: if (d >= d1)
2569: {
2570: hiremainder = 0; q = 1 + divll(d,d1); d = hiremainder;
2571: xv += q * xv1;
2572: }
2573: else
2574: xv += xv1;
2575: /* possible loop exit */
2576: if (d <= 1UL) { xs=1; break; }
2577: /* repeat with inverted roles */
2578: d1 -= d;
2579: if (d1 >= d)
2580: {
2581: hiremainder = 0; q = 1 + divll(d1,d); d1 = hiremainder;
2582: xv1 += q * xv;
2583: }
2584: else
2585: xv1 += xv;
2586: } /* while */
2587:
2588: if (!(f&1)) /* division by 1 postprocessing if needed */
2589: {
2590: if (xs && d==1)
2591: { xv1 += d1 * xv; xs = 0; res = 1UL; }
2592: else if (!xs && d1==1)
2593: { xv += d * xv1; xs = 1; res = 1UL; }
2594: }
2595:
2596: if (xs)
2597: {
2598: *s = -1; *v = xv1; *v1 = xv;
2599: return (res ? res : (d==1 ? 1UL : d1));
2600: }
2601: else
2602: {
2603: *s = 1; *v = xv; *v1 = xv1;
2604: return (res ? res : (d1==1 ? 1UL : d));
2605: }
2606: }
2607:
2608:
2609: ulong
2610: xxgcduu(ulong d, ulong d1, int f,
2611: ulong* u, ulong* u1, ulong* v, ulong* v1, long *s)
2612: {
2613: ulong xu,xu1, xv,xv1, xs, q,res;
2614: LOCAL_HIREMAINDER;
2615:
2616: xs = res = 0;
2617: xu = xv1 = 1UL;
2618: xu1 = xv = 0UL;
2619: while (d1 > 1UL)
2620: {
2621: d -= d1; /* no need to use subll */
2622: if (d >= d1)
2623: {
2624: hiremainder = 0; q = 1 + divll(d,d1); d = hiremainder;
2625: xv += q * xv1;
2626: xu += q * xu1;
2627: }
2628: else
2629: { xv += xv1; xu += xu1; }
2630: /* possible loop exit */
2631: if (d <= 1UL) { xs=1; break; }
2632: /* repeat with inverted roles */
2633: d1 -= d;
2634: if (d1 >= d)
2635: {
2636: hiremainder = 0; q = 1 + divll(d1,d); d1 = hiremainder;
2637: xv1 += q * xv;
2638: xu1 += q * xu;
2639: }
2640: else
2641: { xv1 += xv; xu1 += xu; }
2642: } /* while */
2643:
2644: if (!(f&1)) /* division by 1 postprocessing if needed */
2645: {
2646: if (xs && d==1)
2647: {
2648: xv1 += d1 * xv;
2649: xu1 += d1 * xu;
2650: xs = 0; res = 1UL;
2651: }
2652: else if (!xs && d1==1)
2653: {
2654: xv += d * xv1;
2655: xu += d * xu1;
2656: xs = 1; res = 1UL;
2657: }
2658: }
2659:
2660: if (xs)
2661: {
2662: *s = -1; *u = xu1; *u1 = xu; *v = xv1; *v1 = xv;
2663: return (res ? res : (d==1 ? 1UL : d1));
2664: }
2665: else
2666: {
2667: *s = 1; *u = xu; *u1 = xu1; *v = xv; *v1 = xv1;
2668: return (res ? res : (d1==1 ? 1UL : d));
2669: }
2670: }
2671:
2672: ulong
2673: rgcduu(ulong d, ulong d1, ulong vmax,
2674: ulong* u, ulong* u1, ulong* v, ulong* v1, long *s)
2675: {
2676: ulong xu,xu1, xv,xv1, xs, q, res=0;
2677: int f = 0;
2678: LOCAL_HIREMAINDER;
2679:
2680: if (vmax == 0) vmax = MAXULONG;
2681: xs = res = 0;
2682: xu = xv1 = 1UL;
2683: xu1 = xv = 0UL;
2684: while (d1 > 1UL)
2685: {
2686: d -= d1; /* no need to use subll */
2687: if (d >= d1)
2688: {
2689: hiremainder = 0; q = 1 + divll(d,d1); d = hiremainder;
2690: xv += q * xv1;
2691: xu += q * xu1;
2692: }
2693: else
2694: { xv += xv1; xu += xu1; }
2695: /* possible loop exit */
2696: if (xv > vmax) { f=xs=1; break; }
2697: if (d <= 1UL) { xs=1; break; }
2698: /* repeat with inverted roles */
2699: d1 -= d;
2700: if (d1 >= d)
2701: {
2702: hiremainder = 0; q = 1 + divll(d1,d); d1 = hiremainder;
2703: xv1 += q * xv;
2704: xu1 += q * xu;
2705: }
2706: else
2707: { xv1 += xv; xu1 += xu; }
2708: /* possible loop exit */
2709: if (xv1 > vmax) { f=1; break; }
2710: } /* while */
2711:
2712: if (!(f&1)) /* division by 1 postprocessing if needed */
2713: {
2714: if (xs && d==1)
2715: {
2716: xv1 += d1 * xv;
2717: xu1 += d1 * xu;
2718: xs = 0; res = 1UL;
2719: }
2720: else if (!xs && d1==1)
2721: {
2722: xv += d * xv1;
2723: xu += d * xu1;
2724: xs = 1; res = 1UL;
2725: }
2726: }
2727:
2728: if (xs)
2729: {
2730: *s = -1; *u = xu1; *u1 = xu; *v = xv1; *v1 = xv;
2731: return (res ? res : (d==1 ? 1UL : d1));
2732: }
2733: else
2734: {
2735: *s = 1; *u = xu; *u1 = xu1; *v = xv; *v1 = xv1;
2736: return (res ? res : (d1==1 ? 1UL : d));
2737: }
2738: }
2739:
2740: /*==================================
2741: * lgcdii(d,d1,u,u1,v,v1,vmax)
2742: *==================================*/
2743: /* Lehmer's partial extended gcd algorithm, acting on two t_INT GENs.
2744: *
2745: * Tries to determine, using the leading 2*BITS_IN_LONG significant bits of d
2746: * and a quantity of bits from d1 obtained by a shift of the same displacement,
2747: * as many partial quotients of d/d1 as possible, and assigns to [u,u1;v,v1]
2748: * the product of all the [0, 1; 1, q_j] thus obtained, where the leftmost
2749: * factor arises from the quotient of the first division step.
2750: *
2751: * For use in rational reconstruction, input param vmax can be given a
2752: * nonzero value. In this case, we will return early as soon as v1 > vmax
2753: * (i.e. it is guaranteed that v <= vmax). --2001Jan15
2754: *
2755: * MUST be called with d > d1 > 0, and with d occupying more than one
2756: * significant word (if it doesn't, the caller has no business with us;
2757: * he/she/it should use xgcduu() instead). Returns the number of reduction/
2758: * swap steps carried out, possibly zero, or under certain conditions minus
2759: * that number. When the return value is nonzero, the caller should use the
2760: * returned recurrence matrix to update its own copies of d,d1. When the
2761: * return value is non-positive, and the latest remainder after updating
2762: * turns out to be nonzero, the caller should at once attempt a full division,
2763: * rather than first trying lgcdii() again -- this typically happens when we
2764: * are about to encounter a quotient larger than half a word. (This is not
2765: * detected infallibly -- after a positive return value, it is perfectly
2766: * possible that the next stage will end up needing a full division. After
2767: * a negative return value, however, this is certain, and should be acted
2768: * upon.)
2769: *
2770: * (The sign information, for which xgcduu() has its return argument s, is now
2771: * implicit in the LSB of our return value, and the caller may take advantage
2772: * of the fact that a return value of +-1 implies u==0,u1==v==1 [only v1 pro-
2773: * vides interesting information in this case]. One might also use the fact
2774: * that if the return value is +-2, then u==1, but this is rather marginal.)
2775: *
2776: * If it was not possible to determine even the first quotient, either because
2777: * we're too close to an integer quotient or because the quotient would be
2778: * larger than one word (if the `leading digit' of d1 after shifting is all
2779: * zeros), we return 0 and do not bother to assign anything to the last four
2780: * args.
2781: *
2782: * The division chain might (sometimes) even run to completion. It will be
2783: * up to the caller to detect this case.
2784: *
2785: * This routine does _not_ change d or d1; this will also be up to the caller.
2786: *
2787: * Note that this routine does not know and does not need to know about the
2788: * PARI stack.
2789: */
2790: /*#define DEBUG_LEHMER 1 */
2791: int
2792: lgcdii(ulong* d, ulong* d1,
2793: ulong* u, ulong* u1, ulong* v, ulong* v1,
2794: ulong vmax)
2795: {
2796: /* Strategy: (1) Extract/shift most significant bits. We assume that d
2797: * has at least two significant words, but we can cope with a one-word d1.
2798: * Let dd,dd1 be the most significant dividend word and matching part of the
2799: * divisor.
2800: * (2) Check for overflow on the first division. For our purposes, this
2801: * happens when the upper half of dd1 is zero. (Actually this is detected
2802: * during extraction.)
2803: * (3) Get a fix on the first quotient. We compute q = floor(dd/dd1), which
2804: * is an upper bound for floor(d/d1), and which gives the true value of the
2805: * latter if (and-almost-only-if) the remainder dd' = dd-q*dd1 is >= q.
2806: * (If it isn't, we give up. This is annoying because the subsequent full
2807: * division will repeat some work already done, but it happens very infre-
2808: * quently. Doing the extra-bit-fetch in this case would be awkward.)
2809: * (4) Finish initializations.
2810: *
2811: * The remainder of the action is comparatively boring... The main loop has
2812: * been unrolled once (so we don't swap things and we can apply Jebelean's
2813: * termination conditions which alternatingly take two different forms during
2814: * successive iterations). When we first run out of sufficient bits to form
2815: * a quotient, and have an extra word of each operand, we pull out two whole
2816: * word's worth of dividend bits, and divisor bits of matching significance;
2817: * to these we apply our partial matrix (disregarding overflow because the
2818: * result mod 2^(2*BITS_IN_LONG) will in fact give the correct values), and
2819: * re-extract one word's worth of the current dividend and a matching amount
2820: * of divisor bits. The affair will normally terminate with matrix entries
2821: * just short of a whole word. (We terminate the inner loop before these can
2822: * possibly overflow.)
2823: */
2824: ulong dd,dd1,ddlo,dd1lo, sh,shc; /* `digits', shift count */
2825: ulong xu,xu1, xv,xv1, q,res; /* recurrences, partial quotient, count */
2826: ulong tmp0,tmp1,tmp2,tmpd,tmpu,tmpv; /* temps */
2827: long ld, ld1, lz; /* t_INT effective lengths */
2828: int skip = 0; /* a boolean flag */
2829: LOCAL_OVERFLOW;
2830: LOCAL_HIREMAINDER;
2831:
2832: #ifdef DEBUG_LEHMER
2833: voir(d, -1); voir(d1, -1);
2834: #endif
2835: /* following is just for convenience: vmax==0 means no bound */
2836: if (vmax == 0) vmax = MAXULONG;
2837: ld = lgefint(d); ld1 = lgefint(d1); lz = ld - ld1; /* >= 0 */
2838: if (lz > 1) return 0; /* rare, quick and desperate exit */
2839:
2840: d += 2; d1 += 2; /* point at the leading `digits' */
2841: dd1lo = 0; /* unless we find something better */
2842: sh = bfffo(*d); /* obtain dividend left shift count */
2843:
2844: if (sh)
2845: { /* do the shifting */
2846: shc = BITS_IN_LONG - sh;
2847: if (lz)
2848: { /* dividend longer than divisor */
2849: dd1 = (*d1 >> shc);
2850: if (!(HIGHMASK & dd1)) return 0; /* overflow detected */
2851: if (ld1 > 3)
2852: dd1lo = (*d1 << sh) + (d1[1] >> shc);
2853: else
2854: dd1lo = (*d1 << sh);
2855: }
2856: else
2857: { /* dividend and divisor have the same length */
2858: /* dd1 = shiftl(d1,sh) would have left hiremainder==0, and dd1 != 0. */
2859: dd1 = (*d1 << sh);
2860: if (!(HIGHMASK & dd1)) return 0;
2861: if (ld1 > 3)
2862: {
2863: dd1 += (d1[1] >> shc);
2864: if (ld1 > 4)
2865: dd1lo = (d1[1] << sh) + (d1[2] >> shc);
2866: else
2867: dd1lo = (d1[1] << sh);
2868: }
2869: }
2870: /* following lines assume d to have 2 or more significant words */
2871: dd = (*d << sh) + (d[1] >> shc);
2872: if (ld > 4)
2873: ddlo = (d[1] << sh) + (d[2] >> shc);
2874: else
2875: ddlo = (d[1] << sh);
2876: }
2877: else
2878: { /* no shift needed */
2879: if (lz) return 0; /* div'd longer than div'r: o'flow automatic */
2880: dd1 = *d1;
2881: if (!(HIGHMASK & dd1)) return 0;
2882: if(ld1 > 3) dd1lo = d1[1];
2883: /* assume again that d has another significant word */
2884: dd = *d; ddlo = d[1];
2885: }
2886: #ifdef DEBUG_LEHMER
2887: fprintferr(" %lx:%lx, %lx:%lx\n", dd, ddlo, dd1, dd1lo);
2888: #endif
2889:
2890: /* First subtraction/division stage. (If a subtraction initially suffices,
2891: * we don't divide at all.) If a Jebelean condition is violated, and we
2892: * can't fix it even by looking at the low-order bits in ddlo,dd1lo, we
2893: * give up and ask for a full division. Otherwise we commit the result,
2894: * possibly deciding to re-shift immediately afterwards.
2895: */
2896: dd -= dd1;
2897: if (dd < dd1)
2898: { /* first quotient known to be == 1 */
2899: xv1 = 1UL;
2900: if (!dd) /* !(Jebelean condition), extraspecial case */
2901: { /* note this can actually happen... Now
2902: * q==1 is known, but we underflow already.
2903: * OTOH we've just shortened d by a whole word.
2904: * Thus we feel pleased with ourselves and
2905: * return. (The re-shift code below would
2906: * do so anyway.) */
2907: *u = 0; *v = *u1 = *v1 = 1UL;
2908: return -1; /* Next step will be a full division. */
2909: } /* if !(Jebelean) then */
2910: }
2911: else
2912: { /* division indicated */
2913: hiremainder = 0;
2914: xv1 = 1 + divll(dd, dd1); /* xv1: alternative spelling of `q', here ;) */
2915: dd = hiremainder;
2916: if (dd < xv1) /* !(Jebelean cond'), non-extra special case */
2917: { /* Attempt to complete the division using the
2918: * less significant bits, before skipping right
2919: * past the 1st loop to the reshift stage. */
2920: ddlo = subll(ddlo, mulll(xv1, dd1lo));
2921: dd = subllx(dd, hiremainder);
2922:
2923: /* If we now have an overflow, q was _certainly_ too large. Thanks to
2924: * our decision not to get here unless the original dd1 had bits set in
2925: * the upper half of the word, however, we now do know that the correct
2926: * quotient is in fact q-1. Adjust our data accordingly. */
2927: if (overflow)
2928: {
2929: xv1--;
2930: ddlo = addll(ddlo,dd1lo);
2931: dd = addllx(dd,dd1); /* overflows again which cancels the borrow */
2932: /* ...and fall through to skip=1 below */
2933: }
2934: else
2935: /* Test Jebelean condition anew, at this point using _all_ the extracted
2936: * bits we have. This is clutching at straws; we have a more or less
2937: * even chance of succeeding this time. Note that if we fail, we really
2938: * do not know whether the correct quotient would have been q or some
2939: * smaller value. */
2940: if (!dd && ddlo < xv1) return 0;
2941:
2942: /* Otherwise, we now know that q is correct, but we cannot go into the
2943: * 1st loop. Raise a flag so we'll remember to skip past the loop.
2944: * Get here also after the q-1 adjustment case. */
2945: skip = 1;
2946: } /* if !(Jebelean) then */
2947: }
2948: res = 1;
2949: #ifdef DEBUG_LEHMER
2950: fprintferr(" q = %ld, %lx, %lx\n", xv1, dd1, dd);
2951: #endif
2952: if (xv1 > vmax)
2953: { /* gone past the bound already */
2954: *u = 0UL; *u1 = 1UL; *v = 1UL; *v1 = xv1;
2955: return res;
2956: }
2957: xu = 0UL; xv = xu1 = 1UL;
2958:
2959: /* Some invariants from here across the first loop:
2960: *
2961: * At this point, and again after we are finished with the first loop and
2962: * subsequent conditional, a division and the associated update of the
2963: * recurrence matrix have just been carried out completely. The matrix
2964: * xu,xu1;xv,xv1 has been initialized (or updated, possibly with permuted
2965: * columns), and the current remainder == next divisor (dd at the moment)
2966: * is nonzero (it might be zero here, but then skip will have been set).
2967: *
2968: * After the first loop, or when skip is set already, it will also be the
2969: * case that there aren't sufficiently many bits to continue without re-
2970: * shifting. If the divisor after reshifting is zero, or indeed if it
2971: * doesn't have more than half a word's worth of bits, we will have to
2972: * return at that point. Otherwise, we proceed into the second loop.
2973: *
2974: * Furthermore, when we reach the re-shift stage, dd:ddlo and dd1:dd1lo will
2975: * already reflect the result of applying the current matrix to the old
2976: * ddorig:ddlo and dd1orig:dd1lo. (For the first iteration above, this
2977: * was easy to achieve, and we didn't even need to peek into the (now
2978: * no longer existent!) saved words. After the loop, we'll stop for a
2979: * moment to merge in the ddlo,dd1lo contributions.)
2980: *
2981: * Note that after the first division, even an a priori quotient of 1 cannot
2982: * be trusted until we've checked Jebelean's condition -- it cannot be too
2983: * large, of course, but it might be too small.
2984: */
2985:
2986: if (!skip)
2987: {
2988: for(;;)
2989: {
2990: /* First half of loop divides dd into dd1, and leaves the recurrence
2991: * matrix xu,...,xv1 groomed the wrong way round (xu,xv will be the newer
2992: * entries) when successful. */
2993: tmpd = dd1 - dd;
2994: if (tmpd < dd)
2995: { /* quotient suspected to be 1 */
2996: #ifdef DEBUG_LEHMER
2997: q = 1;
2998: #endif
2999: tmpu = xu + xu1; /* cannot overflow -- everything bounded by
3000: * the original dd during first loop */
3001: tmpv = xv + xv1;
3002: }
3003: else
3004: { /* division indicated */
3005: hiremainder = 0;
3006: q = 1 + divll(tmpd, dd);
3007: tmpd = hiremainder;
3008: tmpu = xu + q*xu1; /* can't overflow, but may need to be undone */
3009: tmpv = xv + q*xv1;
3010: }
3011:
3012: tmp0 = addll(tmpv, xv1);
3013: if ((tmpd < tmpu) || overflow ||
3014: (dd - tmpd < tmp0)) /* !(Jebelean cond.) */
3015: break; /* skip ahead to reshift stage */
3016: else
3017: { /* commit dd1, xu, xv */
3018: res++;
3019: dd1 = tmpd; xu = tmpu; xv = tmpv;
3020: #ifdef DEBUG_LEHMER
3021: fprintferr(" q = %ld, %lx, %lx [%lu,%lu;%lu,%lu]\n",
3022: q, dd, dd1, xu1, xu, xv1, xv);
3023: #endif
3024: if (xv > vmax)
3025: { /* time to return */
3026: *u = xu1; *u1 = xu; *v = xv1; *v1 = xv;
3027: return res;
3028: }
3029: }
3030:
3031: /* Second half of loop divides dd1 into dd, and the matrix returns to its
3032: * normal arrangement. */
3033: tmpd = dd - dd1;
3034: if (tmpd < dd1)
3035: { /* quotient suspected to be 1 */
3036: #ifdef DEBUG_LEHMER
3037: q = 1;
3038: #endif
3039: tmpu = xu1 + xu; /* cannot overflow */
3040: tmpv = xv1 + xv;
3041: }
3042: else
3043: { /* division indicated */
3044: hiremainder = 0;
3045: q = 1 + divll(tmpd, dd1);
3046: tmpd = hiremainder;
3047: tmpu = xu1 + q*xu;
3048: tmpv = xv1 + q*xv;
3049: }
3050:
3051: tmp0 = addll(tmpu, xu);
3052: if ((tmpd < tmpv) || overflow ||
3053: (dd1 - tmpd < tmp0)) /* !(Jebelean cond.) */
3054: break; /* skip ahead to reshift stage */
3055: else
3056: { /* commit dd, xu1, xv1 */
3057: res++;
3058: dd = tmpd; xu1 = tmpu; xv1 = tmpv;
3059: #ifdef DEBUG_LEHMER
3060: fprintferr(" q = %ld, %lx, %lx [%lu,%lu;%lu,%lu]\n",
3061: q, dd1, dd, xu, xu1, xv, xv1);
3062: #endif
3063: if (xv1 > vmax)
3064: { /* time to return */
3065: *u = xu; *u1 = xu1; *v = xv; *v1 = xv1;
3066: return res;
3067: }
3068: }
3069:
3070: } /* end of first loop */
3071:
3072: /* Intermezzo: update dd:ddlo, dd1:dd1lo. (But not if skip is set.) */
3073:
3074: if (res&1)
3075: {
3076: /* after failed division in 1st half of loop:
3077: * [dd1:dd1lo,dd:ddlo] = [ddorig:ddlo,dd1orig:dd1lo]
3078: * * [ -xu, xu1 ; xv, -xv1 ]
3079: * (Actually, we only multiply [ddlo,dd1lo] onto the matrix and
3080: * add the high-order remainders + overflows onto [dd1,dd].)
3081: */
3082: tmp1 = mulll(ddlo, xu); tmp0 = hiremainder;
3083: tmp1 = subll(mulll(dd1lo,xv), tmp1);
3084: dd1 += subllx(hiremainder, tmp0);
3085: tmp2 = mulll(ddlo, xu1); tmp0 = hiremainder;
3086: ddlo = subll(tmp2, mulll(dd1lo,xv1));
3087: dd += subllx(tmp0, hiremainder);
3088: dd1lo = tmp1;
3089: }
3090: else
3091: {
3092: /* after failed division in 2nd half of loop:
3093: * [dd:ddlo,dd1:dd1lo] = [ddorig:ddlo,dd1orig:dd1lo]
3094: * * [ xu1, -xu ; -xv1, xv ]
3095: * (Actually, we only multiply [ddlo,dd1lo] onto the matrix and
3096: * add the high-order remainders + overflows onto [dd,dd1].)
3097: */
3098: tmp1 = mulll(ddlo, xu1); tmp0 = hiremainder;
3099: tmp1 = subll(tmp1, mulll(dd1lo,xv1));
3100: dd += subllx(tmp0, hiremainder);
3101: tmp2 = mulll(ddlo, xu); tmp0 = hiremainder;
3102: dd1lo = subll(mulll(dd1lo,xv), tmp2);
3103: dd1 += subllx(hiremainder, tmp0);
3104: ddlo = tmp1;
3105: }
3106: #ifdef DEBUG_LEHMER
3107: fprintferr(" %lx:%lx, %lx:%lx\n", dd, ddlo, dd1, dd1lo);
3108: #endif
3109: } /* end of skip-pable section: get here also, with res==1, when there
3110: * was a problem immediately after the very first division. */
3111:
3112: /* Re-shift. Note: the shift count _can_ be zero, viz. under the following
3113: * precise conditions: The original dd1 had its topmost bit set, so the 1st
3114: * q was 1, and after subtraction, dd had its topmost bit unset. If now
3115: * dd==0, we'd have taken the return exit already, so we couldn't have got
3116: * here. If not, then it must have been the second division which has gone
3117: * amiss (because dd1 was very close to an exact multiple of the remainder
3118: * dd value, so this will be very rare). At this point, we'd have a fairly
3119: * slim chance of fixing things by re-examining dd1:dd1lo vs. dd:ddlo, but
3120: * this is not guaranteed to work. Instead of trying, we return at once.
3121: * The caller will see to it that the initial subtraction is re-done using
3122: * _all_ the bits of both operands, which already helps, and the next round
3123: * will either be a full division (if dd occupied a halfword or less), or
3124: * another llgcdii() first step. In the latter case, since we try a little
3125: * harder during our first step, we may actually be able to fix the problem,
3126: * and get here again with improved low-order bits and with another step
3127: * under our belt. Otherwise we'll have given up above and forced a full-
3128: * blown division.
3129: *
3130: * If res is even, the shift count _cannot_ be zero. (The first step forces
3131: * a zero into the remainder's MSB, and all subsequent remainders will have
3132: * inherited it.)
3133: *
3134: * The re-shift stage exits if the next divisor has at most half a word's
3135: * worth of bits.
3136: *
3137: * For didactic reasons, the second loop will be arranged in the same way
3138: * as the first -- beginning with the division of dd into dd1, as if res
3139: * was odd. To cater for this, if res is actually even, we swap things
3140: * around during reshifting. (During the second loop, the parity of res
3141: * does not matter; we know in which half of the loop we are when we decide
3142: * to return.)
3143: */
3144: #ifdef DEBUG_LEHMER
3145: fprintferr("(sh)");
3146: #endif
3147:
3148: if (res&1)
3149: { /* after odd number of division(s) */
3150: if (dd1 && (sh = bfffo(dd1)))
3151: {
3152: shc = BITS_IN_LONG - sh;
3153: dd = (ddlo >> shc) + (dd << sh);
3154: if (!(HIGHMASK & dd))
3155: {
3156: *u = xu; *u1 = xu1; *v = xv; *v1 = xv1;
3157: return -res; /* full division asked for */
3158: }
3159: dd1 = (dd1lo >> shc) + (dd1 << sh);
3160: }
3161: else
3162: { /* time to return: <= 1 word left, or sh==0 */
3163: *u = xu; *u1 = xu1; *v = xv; *v1 = xv1;
3164: return res;
3165: }
3166: }
3167: else
3168: { /* after even number of divisions */
3169: if (dd)
3170: {
3171: sh = bfffo(dd); /* Known to be positive. */
3172: shc = BITS_IN_LONG - sh;
3173: /* dd:ddlo will become the new dd1, and v.v. */
3174: tmpd = (ddlo >> shc) + (dd << sh);
3175: dd = (dd1lo >> shc) + (dd1 << sh);
3176: dd1 = tmpd;
3177: /* This has completed the swap; now dd is again the current divisor.
3178: * The following test originally inspected dd1 -- a most subtle and
3179: * most annoying bug. The Management. */
3180: if (HIGHMASK & dd)
3181: {
3182: /* recurrence matrix is the wrong way round; swap it. */
3183: tmp0 = xu; xu = xu1; xu1 = tmp0;
3184: tmp0 = xv; xv = xv1; xv1 = tmp0;
3185: }
3186: else
3187: {
3188: /* recurrence matrix is the wrong way round; fix this. */
3189: *u = xu1; *u1 = xu; *v = xv1; *v1 = xv;
3190: return -res; /* full division asked for */
3191: }
3192: }
3193: else
3194: { /* time to return: <= 1 word left */
3195: *u = xu1; *u1 = xu; *v = xv1; *v1 = xv;
3196: return res;
3197: }
3198: } /* end reshift */
3199:
3200: #ifdef DEBUG_LEHMER
3201: fprintferr(" %lx:%lx, %lx:%lx\n", dd, ddlo, dd1, dd1lo);
3202: #endif
3203:
3204: /* The Second Loop. Rip-off of the first, but we now check for overflow
3205: * in the recurrences. Returns instead of breaking when we cannot fix the
3206: * quotient any longer.
3207: */
3208:
3209: for(;;)
3210: {
3211: /* First half of loop divides dd into dd1, and leaves the recurrence
3212: * matrix xu,...,xv1 groomed the wrong way round (xu,xv will be the newer
3213: * entries) when successful. */
3214: tmpd = dd1 - dd;
3215: if (tmpd < dd)
3216: { /* quotient suspected to be 1 */
3217: #ifdef DEBUG_LEHMER
3218: q = 1;
3219: #endif
3220: tmpu = xu + xu1;
3221: tmpv = addll(xv, xv1); /* xv,xv1 will overflow first */
3222: tmp1 = overflow;
3223: }
3224: else
3225: { /* division indicated */
3226: hiremainder = 0;
3227: q = 1 + divll(tmpd, dd);
3228: tmpd = hiremainder;
3229: tmpu = xu + q*xu1;
3230: tmpv = addll(xv, mulll(q,xv1));
3231: tmp1 = overflow | hiremainder;
3232: }
3233:
3234: tmp0 = addll(tmpv, xv1);
3235: if ((tmpd < tmpu) || overflow || tmp1 ||
3236: (dd - tmpd < tmp0)) /* !(Jebelean cond.) */
3237: {
3238: /* The recurrence matrix has not yet been warped... */
3239: *u = xu; *u1 = xu1; *v = xv; *v1 = xv1;
3240: break;
3241: }
3242: /* commit dd1, xu, xv */
3243: res++;
3244: dd1 = tmpd; xu = tmpu; xv = tmpv;
3245: #ifdef DEBUG_LEHMER
3246: fprintferr(" q = %ld, %lx, %lx\n", q, dd, dd1);
3247: #endif
3248: if (xv > vmax)
3249: { /* time to return */
3250: *u = xu1; *u1 = xu; *v = xv1; *v1 = xv;
3251: return res;
3252: }
3253:
3254: /* Second half of loop divides dd1 into dd, and the matrix returns to its
3255: * normal arrangement. */
3256: tmpd = dd - dd1;
3257: if (tmpd < dd1)
3258: { /* quotient suspected to be 1 */
3259: #ifdef DEBUG_LEHMER
3260: q = 1;
3261: #endif
3262: tmpu = xu1 + xu;
3263: tmpv = addll(xv1, xv);
3264: tmp1 = overflow;
3265: }
3266: else
3267: { /* division indicated */
3268: hiremainder = 0;
3269: q = 1 + divll(tmpd, dd1);
3270: tmpd = hiremainder;
3271: tmpu = xu1 + q*xu;
3272: tmpv = addll(xv1, mulll(q, xv));
3273: tmp1 = overflow | hiremainder;
3274: }
3275:
3276: tmp0 = addll(tmpu, xu);
3277: if ((tmpd < tmpv) || overflow || tmp1 ||
3278: (dd1 - tmpd < tmp0)) /* !(Jebelean cond.) */
3279: {
3280: /* The recurrence matrix has not yet been unwarped, so it is
3281: * the wrong way round; fix this. */
3282: *u = xu1; *u1 = xu; *v = xv1; *v1 = xv;
3283: break;
3284: }
3285:
3286: res++; /* commit dd, xu1, xv1 */
3287: dd = tmpd; xu1 = tmpu; xv1 = tmpv;
3288: #ifdef DEBUG_LEHMER
3289: fprintferr(" q = %ld, %lx, %lx\n", q, dd1, dd);
3290: #endif
3291: if (xv1 > vmax)
3292: { /* time to return */
3293: *u = xu; *u1 = xu1; *v = xv; *v1 = xv1;
3294: return res;
3295: }
3296: } /* end of second loop */
3297:
3298: return res;
3299: }
3300:
3301: /*==================================
3302: * invmod(a,b,res)
3303: *==================================
3304: * If a is invertible, return 1, and set res = a^{ -1 }
3305: * Otherwise, return 0, and set res = gcd(a,b)
3306: *
3307: * This is sufficiently different from bezout() to be implemented separately
3308: * instead of having a bunch of extra conditionals in a single function body
3309: * to meet both purposes.
3310: */
3311:
3312: int
3313: invmod(GEN a, GEN b, GEN *res)
3314: {
3315: GEN v,v1,d,d1,q,r;
3316: long av,av1,lim;
3317: long s;
3318: ulong g;
3319: ulong xu,xu1,xv,xv1; /* Lehmer stage recurrence matrix */
3320: int lhmres; /* Lehmer stage return value */
3321:
3322: if (typ(a) != t_INT || typ(b) != t_INT) err(arither1);
3323: if (!signe(b)) { *res=absi(a); return 0; }
3324: av = avma;
3325: if (lgefint(b) == 3) /* single-word affair */
3326: {
3327: d1 = modiu(a, (ulong)(b[2]));
3328: if (d1 == gzero)
3329: {
3330: if (b[2] == 1L)
3331: { *res = gzero; return 1; }
3332: else
3333: { *res = absi(b); return 0; }
3334: }
3335: g = xgcduu((ulong)(b[2]), (ulong)(d1[2]), 1, &xv, &xv1, &s);
3336: #ifdef DEBUG_LEHMER
3337: fprintferr(" <- %lu,%lu\n", (ulong)(b[2]), (ulong)(d1[2]));
3338: fprintferr(" -> %lu,%ld,%lu; %lx\n", g,s,xv1,avma);
3339: #endif
3340: avma = av;
3341: if (g != 1UL) { *res = utoi(g); return 0; }
3342: xv = xv1 % (ulong)(b[2]); if (s < 0) xv = ((ulong)(b[2])) - xv;
3343: *res = utoi(xv); return 1;
3344: }
3345:
3346: (void)new_chunk(lgefint(b));
3347: d = absi(b); d1 = modii(a,d);
3348:
3349: v=gzero; v1=gun; /* general case */
3350: #ifdef DEBUG_LEHMER
3351: fprintferr("INVERT: -------------------------\n");
3352: output(d1);
3353: #endif
3354: av1 = avma; lim = stack_lim(av,1);
3355:
3356: while (lgefint(d) > 3 && signe(d1))
3357: {
3358: #ifdef DEBUG_LEHMER
3359: fprintferr("Calling Lehmer:\n");
3360: #endif
3361: lhmres = lgcdii((ulong*)d, (ulong*)d1, &xu, &xu1, &xv, &xv1, MAXULONG);
3362: if (lhmres != 0) /* check progress */
3363: { /* apply matrix */
3364: #ifdef DEBUG_LEHMER
3365: fprintferr("Lehmer returned %d [%lu,%lu;%lu,%lu].\n",
3366: lhmres, xu, xu1, xv, xv1);
3367: #endif
3368: if ((lhmres == 1) || (lhmres == -1))
3369: {
3370: if (xv1 == 1)
3371: {
3372: r = subii(d,d1); d=d1; d1=r;
3373: a = subii(v,v1); v=v1; v1=a;
3374: }
3375: else
3376: {
3377: r = subii(d, mului(xv1,d1)); d=d1; d1=r;
3378: a = subii(v, mului(xv1,v1)); v=v1; v1=a;
3379: }
3380: }
3381: else
3382: {
3383: r = subii(muliu(d,xu), muliu(d1,xv));
3384: a = subii(muliu(v,xu), muliu(v1,xv));
3385: d1 = subii(muliu(d,xu1), muliu(d1,xv1)); d = r;
3386: v1 = subii(muliu(v,xu1), muliu(v1,xv1)); v = a;
3387: if (lhmres&1)
3388: {
3389: setsigne(d,-signe(d));
3390: setsigne(v,-signe(v));
3391: }
3392: else
3393: {
3394: if (signe(d1)) { setsigne(d1,-signe(d1)); }
3395: setsigne(v1,-signe(v1));
3396: }
3397: }
3398: }
3399: #ifdef DEBUG_LEHMER
3400: else
3401: fprintferr("Lehmer returned 0.\n");
3402: output(d); output(d1); output(v); output(v1);
3403: sleep(1);
3404: #endif
3405:
3406: if (lhmres <= 0 && signe(d1))
3407: {
3408: q = dvmdii(d,d1,&r);
3409: #ifdef DEBUG_LEHMER
3410: fprintferr("Full division:\n");
3411: printf(" q = "); output(q); sleep (1);
3412: #endif
3413: a = subii(v,mulii(q,v1));
3414: v=v1; v1=a;
3415: d=d1; d1=r;
3416: }
3417: if (low_stack(lim, stack_lim(av,1)))
3418: {
3419: GEN *gptr[4]; gptr[0]=&d; gptr[1]=&d1; gptr[2]=&v; gptr[3]=&v1;
3420: if(DEBUGMEM>1) err(warnmem,"invmod");
3421: gerepilemany(av1,gptr,4);
3422: }
3423: } /* end while */
3424:
3425: /* Postprocessing - final sprint */
3426: if (signe(d1))
3427: {
3428: /* Assertions: lgefint(d)==lgefint(d1)==3, and
3429: * gcd(d,d1) is nonzero and fits into one word
3430: */
3431: g = xxgcduu((ulong)(d[2]), (ulong)(d1[2]), 1, &xu, &xu1, &xv, &xv1, &s);
3432: #ifdef DEBUG_LEHMER
3433: output(d);output(d1);output(v);output(v1);
3434: fprintferr(" <- %lu,%lu\n", (ulong)(d[2]), (ulong)(d1[2]));
3435: fprintferr(" -> %lu,%ld,%lu; %lx\n", g,s,xv1,avma);
3436: #endif
3437: if (g != 1UL) { avma = av; *res = utoi(g); return 0; }
3438: /* (From the xgcduu() blurb:)
3439: * For finishing the multiword modinv, we now have to multiply the
3440: * returned matrix (with properly adjusted signs) onto the values
3441: * v' and v1' previously obtained from the multiword division steps.
3442: * Actually, it is sufficient to take the scalar product of [v',v1']
3443: * with [u1,-v1], and change the sign if s==1.
3444: */
3445: v = subii(muliu(v,xu1),muliu(v1,xv1));
3446: if (s > 0) setsigne(v,-signe(v));
3447: avma = av; *res = modii(v,b);
3448: #ifdef DEBUG_LEHMER
3449: output(*res); fprintfderr("============================Done.\n");
3450: sleep(1);
3451: #endif
3452: return 1;
3453: }
3454: /* get here when the final sprint was skipped (d1 was zero already) */
3455: avma = av;
3456: if (!egalii(d,gun)) { *res = icopy(d); return 0; }
3457: *res = modii(v,b);
3458: #ifdef DEBUG_LEHMER
3459: output(*res); fprintferr("============================Done.\n");
3460: sleep(1);
3461: #endif
3462: return 1;
3463: }
3464:
3465: /*==================================
3466: * bezout(a,b,pu,pv)
3467: *==================================
3468: * Return g = gcd(a,b) >= 0, and assign GENs u,v through pointers pu,pv
3469: * such that g = u*a + v*b.
3470: * Special cases:
3471: * a == b == 0 ==> pick u=1, v=0
3472: * a != 0 == b ==> keep v=0
3473: * a == 0 != b ==> keep u=0
3474: * |a| == |b| != 0 ==> keep u=0, set v=+-1
3475: * Assignments through pu,pv will be suppressed when the corresponding
3476: * pointer is NULL (but the computations will happen nonetheless).
3477: */
3478:
3479: GEN
3480: bezout(GEN a, GEN b, GEN *pu, GEN *pv)
3481: {
3482: GEN t,u,u1,v,v1,d,d1,q,r;
3483: GEN *pt;
3484: long av,av1,lim;
3485: long s;
3486: int sa, sb;
3487: ulong g;
3488: ulong xu,xu1,xv,xv1; /* Lehmer stage recurrence matrix */
3489: int lhmres; /* Lehmer stage return value */
3490:
3491: if (typ(a) != t_INT || typ(b) != t_INT) err(arither1);
3492: s = absi_cmp(a,b);
3493: if (s < 0)
3494: {
3495: t=b; b=a; a=t;
3496: pt=pu; pu=pv; pv=pt;
3497: }
3498: /* now |a| >= |b| */
3499:
3500: sa = signe(a); sb = signe(b);
3501: if (!sb)
3502: {
3503: if (pv != NULL) *pv = gzero;
3504: switch(sa)
3505: {
3506: case 0:
3507: if (pu != NULL) *pu = gun; return gzero;
3508: case 1:
3509: if (pu != NULL) *pu = gun; return icopy(a);
3510: case -1:
3511: if (pu != NULL) *pu = negi(gun); return(negi(a));
3512: }
3513: }
3514: if (s == 0) /* |a| == |b| != 0 */
3515: {
3516: if (pu != NULL) *pu = gzero;
3517: if (sb > 0)
3518: {
3519: if (pv != NULL) *pv = gun; return icopy(b);
3520: }
3521: else
3522: {
3523: if (pv != NULL) *pv = negi(gun); return(negi(b));
3524: }
3525: }
3526: /* now |a| > |b| > 0 */
3527:
3528: if (lgefint(a) == 3) /* single-word affair */
3529: {
3530: g = xxgcduu((ulong)(a[2]), (ulong)(b[2]), 0, &xu, &xu1, &xv, &xv1, &s);
3531: sa = s > 0 ? sa : -sa;
3532: sb = s > 0 ? -sb : sb;
3533: if (pu != NULL)
3534: {
3535: if (xu == 0) *pu = gzero; /* can happen when b divides a */
3536: else if (xu == 1) *pu = sa < 0 ? negi(gun) : gun;
3537: else if (xu == 2) *pu = sa < 0 ? negi(gdeux) : gdeux;
3538: else
3539: {
3540: *pu = cgeti(3);
3541: (*pu)[1] = evalsigne(sa)|evallgefint(3);
3542: (*pu)[2] = xu;
3543: }
3544: }
3545: if (pv != NULL)
3546: {
3547: if (xv == 1) *pv = sb < 0 ? negi(gun) : gun;
3548: else if (xv == 2) *pv = sb < 0 ? negi(gdeux) : gdeux;
3549: else
3550: {
3551: *pv = cgeti(3);
3552: (*pv)[1] = evalsigne(sb)|evallgefint(3);
3553: (*pv)[2] = xv;
3554: }
3555: }
3556: if (g == 1) return gun;
3557: else if (g == 2) return gdeux;
3558: else
3559: {
3560: r = cgeti(3);
3561: r[1] = evalsigne(1)|evallgefint(3);
3562: r[2] = g;
3563: return r;
3564: }
3565: }
3566:
3567: /* general case */
3568: av = avma;
3569: (void)new_chunk(lgefint(b) + (lgefint(a)<<1)); /* room for u,v,gcd */
3570: /* if a is significantly larger than b, calling lgcdii() is not the best
3571: * way to start -- reduce a mod b first
3572: */
3573: if (lgefint(a) > lgefint(b))
3574: {
3575: d = absi(b), q = dvmdii(absi(a), d, &d1);
3576: if (!signe(d1)) /* a == qb */
3577: {
3578: avma = av;
3579: if (pu != NULL) *pu = gzero;
3580: if (pv != NULL) *pv = sb < 0 ? negi(gun) : gun;
3581: return (icopy(d));
3582: }
3583: else
3584: {
3585: u = gzero;
3586: u1 = v = gun;
3587: v1 = negi(q);
3588: }
3589: /* if this results in lgefint(d) == 3, will fall past main loop */
3590: }
3591: else
3592: {
3593: d = absi(a); d1 = absi(b);
3594: u = v1 = gun; u1 = v = gzero;
3595: }
3596: av1 = avma; lim = stack_lim(av, 1);
3597:
3598: /* main loop is almost identical to that of invmod() */
3599: while (lgefint(d) > 3 && signe(d1))
3600: {
3601: lhmres = lgcdii((ulong *)d, (ulong *)d1, &xu, &xu1, &xv, &xv1, MAXULONG);
3602: if (lhmres != 0) /* check progress */
3603: { /* apply matrix */
3604: if ((lhmres == 1) || (lhmres == -1))
3605: {
3606: if (xv1 == 1)
3607: {
3608: r = subii(d,d1); d=d1; d1=r;
3609: a = subii(u,u1); u=u1; u1=a;
3610: a = subii(v,v1); v=v1; v1=a;
3611: }
3612: else
3613: {
3614: r = subii(d, mului(xv1,d1)); d=d1; d1=r;
3615: a = subii(u, mului(xv1,u1)); u=u1; u1=a;
3616: a = subii(v, mului(xv1,v1)); v=v1; v1=a;
3617: }
3618: }
3619: else
3620: {
3621: r = subii(muliu(d,xu), muliu(d1,xv));
3622: d1 = subii(muliu(d,xu1), muliu(d1,xv1)); d = r;
3623: a = subii(muliu(u,xu), muliu(u1,xv));
3624: u1 = subii(muliu(u,xu1), muliu(u1,xv1)); u = a;
3625: a = subii(muliu(v,xu), muliu(v1,xv));
3626: v1 = subii(muliu(v,xu1), muliu(v1,xv1)); v = a;
3627: if (lhmres&1)
3628: {
3629: setsigne(d,-signe(d));
3630: setsigne(u,-signe(u));
3631: setsigne(v,-signe(v));
3632: }
3633: else
3634: {
3635: if (signe(d1)) { setsigne(d1,-signe(d1)); }
3636: setsigne(u1,-signe(u1));
3637: setsigne(v1,-signe(v1));
3638: }
3639: }
3640: }
3641: if (lhmres <= 0 && signe(d1))
3642: {
3643: q = dvmdii(d,d1,&r);
3644: #ifdef DEBUG_LEHMER
3645: fprintferr("Full division:\n");
3646: printf(" q = "); output(q); sleep (1);
3647: #endif
3648: a = subii(u,mulii(q,u1));
3649: u=u1; u1=a;
3650: a = subii(v,mulii(q,v1));
3651: v=v1; v1=a;
3652: d=d1; d1=r;
3653: }
3654: if (low_stack(lim, stack_lim(av,1)))
3655: {
3656: GEN *gptr[6]; gptr[0]=&d; gptr[1]=&d1; gptr[2]=&u; gptr[3]=&u1;
3657: gptr[4]=&v; gptr[5]=&v1;
3658: if(DEBUGMEM>1) err(warnmem,"bezout");
3659: gerepilemany(av1,gptr,6);
3660: }
3661: } /* end while */
3662:
3663: /* Postprocessing - final sprint */
3664: if (signe(d1))
3665: {
3666: /* Assertions: lgefint(d)==lgefint(d1)==3, and
3667: * gcd(d,d1) is nonzero and fits into one word
3668: */
3669: g = xxgcduu((ulong)(d[2]), (ulong)(d1[2]), 0, &xu, &xu1, &xv, &xv1, &s);
3670: u = subii(muliu(u,xu), muliu(u1, xv));
3671: v = subii(muliu(v,xu), muliu(v1, xv));
3672: if (s < 0) { sa = -sa; sb = -sb; }
3673: avma = av;
3674: if (pu != NULL) *pu = sa < 0 ? negi(u) : icopy(u);
3675: if (pv != NULL) *pv = sb < 0 ? negi(v) : icopy(v);
3676: if (g == 1) return gun;
3677: else if (g == 2) return gdeux;
3678: else
3679: {
3680: r = cgeti(3);
3681: r[1] = evalsigne(1)|evallgefint(3);
3682: r[2] = g;
3683: return r;
3684: }
3685: }
3686: /* get here when the final sprint was skipped (d1 was zero already).
3687: * Now the matrix is final, and d contains the gcd.
3688: */
3689: avma = av;
3690: if (pu != NULL) *pu = sa < 0 ? negi(u) : icopy(u);
3691: if (pv != NULL) *pv = sb < 0 ? negi(v) : icopy(v);
3692: return icopy(d);
3693: }
3694:
3695: /*==================================
3696: * cbezout(a,b,uu,vv)
3697: *==================================
3698: * Same as bezout() but for C longs.
3699: * Return g = gcd(a,b) >= 0, and assign longs u,v through pointers uu,vv
3700: * such that g = u*a + v*b.
3701: * Special cases:
3702: * a == b == 0 ==> pick u=1, v=0 (and return 1, surprisingly)
3703: * a != 0 == b ==> keep v=0
3704: * a == 0 != b ==> keep u=0
3705: * |a| == |b| != 0 ==> keep u=0, set v=+-1
3706: * Assignments through uu,vv happen unconditionally; non-NULL pointers
3707: * _must_ be used.
3708: */
3709: long
3710: cbezout(long a,long b,long *uu,long *vv)
3711: {
3712: long s,*t;
3713: ulong d = labs(a), d1 = labs(b);
3714: ulong r,u,u1,v,v1;
3715:
3716: #ifdef DEBUG_CBEZOUT
3717: fprintferr("> cbezout(%ld,%ld,%p,%p)\n", a, b, (void *)uu, (void *)vv);
3718: #endif
3719: if (!b)
3720: {
3721: *vv=0L;
3722: if (!a)
3723: {
3724: *uu=1L;
3725: #ifdef DEBUG_CBEZOUT
3726: fprintferr("< %ld (%ld, %ld)\n", 1L, *uu, *vv);
3727: #endif
3728: return 0L;
3729: }
3730: *uu = a < 0 ? -1L : 1L;
3731: #ifdef DEBUG_CBEZOUT
3732: fprintferr("< %ld (%ld, %ld)\n", (long)d, *uu, *vv);
3733: #endif
3734: return (long)d;
3735: }
3736: else if (!a || (d == d1))
3737: {
3738: *uu = 0L; *vv = b < 0 ? -1L : 1L;
3739: #ifdef DEBUG_CBEZOUT
3740: fprintferr("< %ld (%ld, %ld)\n", (long)d1, *uu, *vv);
3741: #endif
3742: return (long)d1;
3743: }
3744: else if (d == 1) /* frequently used by nfinit */
3745: {
3746: *uu = a; *vv = 0L;
3747: #ifdef DEBUG_CBEZOUT
3748: fprintferr("< %ld (%ld, %ld)\n", 1L, *uu, *vv);
3749: #endif
3750: return 1L;
3751: }
3752: else if (d < d1)
3753: {
3754: /* bug in gcc-2.95.3:
3755: * s = a; a = b; b = s; produces wrong result a = b. This is OK: */
3756: { long _x = a; a = b; b = _x; } /* in order to keep the right signs */
3757: r = d; d = d1; d1 = r;
3758: t = uu; uu = vv; vv = t;
3759: #ifdef DEBUG_CBEZOUT
3760: fprintferr(" swapping\n");
3761: #endif
3762: }
3763: /* d > d1 > 0 */
3764: r = xxgcduu(d, d1, 0, &u, &u1, &v, &v1, &s);
3765: if (s < 0)
3766: {
3767: *uu = a < 0 ? u : -(long)u;
3768: *vv = b < 0 ? -(long)v : v;
3769: }
3770: else
3771: {
3772: *uu = a < 0 ? -(long)u : u;
3773: *vv = b < 0 ? v : -(long)v;
3774: }
3775: #ifdef DEBUG_CBEZOUT
3776: fprintferr("< %ld (%ld, %ld)\n", (long)r, *uu, *vv);
3777: #endif
3778: return (long)r;
3779: }
3780:
3781: /*==========================================================
3782: * ratlift(GEN x, GEN m, GEN *a, GEN *b, GEN amax, GEN bmax)
3783: *==========================================================
3784: * Reconstruct rational number from its residue x mod m
3785: * Given t_INT x, m, amax>=0, bmax>0 such that
3786: * 0 <= x < m; 2*amax*bmax < m
3787: * attempts to find t_INT a, b such that
3788: * (1) a = b*x (mod m)
3789: * (2) |a| <= amax, 0 < b <= bmax
3790: * (3) gcd(m, b) = gcd(a, b)
3791: * If unsuccessful, it will return 0 and leave a,b unchanged (and
3792: * caller may deduce no such a,b exist). If successful, sets a,b
3793: * and returns 1. If there exist a,b satisfying (1), (2), and
3794: * (3') gcd(m, b) = 1
3795: * then they are uniquely determined subject to (1),(2) and
3796: * (3'') gcd(a, b) = 1,
3797: * and will be returned by the routine. (The caller may wish to
3798: * check gcd(a,b)==1, either directly or based on known prime
3799: * divisors of m, depending on the application.)
3800: * Reference:
3801: @article {MR97c:11116,
3802: AUTHOR = {Collins, George E. and Encarnaci{\'o}n, Mark J.},
3803: TITLE = {Efficient rational number reconstruction},
3804: JOURNAL = {J. Symbolic Comput.},
3805: FJOURNAL = {Journal of Symbolic Computation},
3806: VOLUME = {20},
3807: YEAR = {1995},
3808: NUMBER = {3},
3809: PAGES = {287--297},
3810: ISSN = {0747-7171},
3811: MRCLASS = {11Y16 (68Q40)},
3812: MRNUMBER = {97c:11116},
3813: MRREVIEWER = {Maurice Mignotte},
3814: }
3815: * Preprint available from:
3816: * ftp://ftp.risc.uni-linz.ac.at/pub/techreports/1994/94-64.ps.gz
3817: */
3818:
3819: /* #define DEBUG_RATLIFT */
3820:
3821: int
3822: ratlift(GEN x, GEN m, GEN *a, GEN *b, GEN amax, GEN bmax)
3823: {
3824: GEN d,d1,v,v1,q,r;
3825: ulong av = avma, av1, lim;
3826: long lb,lr,lbb,lbr,s,s0;
3827: ulong vmax;
3828: ulong xu,xu1,xv,xv1; /* Lehmer stage recurrence matrix */
3829: int lhmres; /* Lehmer stage return value */
3830:
3831: if ((typ(x) | typ(m) | typ(amax) | typ(bmax)) != t_INT) err(arither1);
3832: if (signe(bmax) <= 0)
3833: err(talker, "ratlift: bmax must be > 0, found\n\tbmax=%Z\n", bmax);
3834: if (signe(amax) < 0)
3835: err(talker, "ratilft: amax must be >= 0, found\n\tamax=%Z\n", amax);
3836: /* check 2*amax*bmax < m */
3837: if (cmpii(shifti(mulii(amax, bmax), 1), m) >= 0)
3838: err(talker, "ratlift: must have 2*amax*bmax < m, found\n\tamax=%Z\n\tbmax=%Z\n\tm=%Z\n", amax,bmax,m);
3839: /* we _could_ silently replace x with modii(x,m) instead of the following,
3840: * but let's leave this up to the caller
3841: */
3842: avma = av; s = signe(x);
3843: if (s < 0 || cmpii(x,m) >= 0)
3844: err(talker, "ratlift: must have 0 <= x < m, found\n\tx=%Z\n\tm=%Z\n", x,m);
3845:
3846: /* special cases x=0 and/or amax=0 */
3847: if (s == 0)
3848: {
3849: if (a != NULL) *a = gzero;
3850: if (b != NULL) *b = gun;
3851: return 1;
3852: }
3853: else if (signe(amax)==0)
3854: return 0;
3855: /* assert: m > x > 0, amax > 0 */
3856:
3857: /* check here whether a=x, b=1 is a solution */
3858: if (cmpii(x,amax) <= 0)
3859: {
3860: if (a != NULL) *a = icopy(x);
3861: if (b != NULL) *b = gun;
3862: return 1;
3863: }
3864:
3865: /* There is no special case for single-word numbers since this is
3866: * mainly meant to be used with large moduli.
3867: */
3868: (void)new_chunk(lgefint(bmax) + lgefint(amax)); /* room for a,b */
3869: d = m; d1 = x;
3870: v = gzero; v1 = gun;
3871: /* assert d1 > amax, v1 <= bmax here */
3872: lb = lgefint(bmax);
3873: lbb = bfffo(bmax[2]);
3874: s = 1;
3875: av1 = avma; lim = stack_lim(av, 1);
3876:
3877: /* general case: Euclidean division chain starting with m div x, and
3878: * with bounds on the sequence of convergents' denoms v_j.
3879: * Just to be different from what invmod and bezout are doing, we work
3880: * here with the all-nonnegative matrices [u,u1;v,v1]=prod_j([0,1;1,q_j]).
3881: * Loop invariants:
3882: * (a) (sign)*[-v,v1]*x = [d,d1] (mod m) (componentwise)
3883: * (sign initially +1, changes with each Euclidean step)
3884: * so [a,b] will be obtained in the form [-+d,v] or [+-d1,v1];
3885: * this congruence is a consequence of
3886: * (b) [x,m]~ = [u,u1;v,v1]*[d1,d]~,
3887: * where u,u1 is the usual numerator sequence starting with 1,0
3888: * instead of 0,1 (just multiply the eqn on the left by the inverse
3889: * matrix, which is det*[v1,-u1;-v,u], where "det" is the same as the
3890: * "(sign)" in (a)). From m = v*d1 + v1*d and
3891: * (c) d > d1 >= 0, 0 <= v < v1,
3892: * we have d >= m/(2*v1), so while v1 remains smaller than m/(2*amax),
3893: * the pair [-(sign)*d,v] satisfies (1) but violates (2) (d > amax).
3894: * Conversely, v1 > bmax indicates that no further solutions will be
3895: * forthcoming; [-(sign)*d,v] will be the last, and first, candidate.
3896: * Thus there's at most one point in the chain division where a solution
3897: * can live: v < bmax, v1 >= m/(2*amax) > bmax, and this is acceptable
3898: * iff in fact d <= amax (e.g. m=221, x=34 or 35, amax=bmax=10 fail on
3899: * this count while x=32,33,36,37 succeed). However, a division may leave
3900: * a zero residue before we ever reach this point (consider m=210, x=35,
3901: * amax=bmax=10), and our caller may find that gcd(d,v) > 1 (numerous
3902: * examples -- keep m=210 and consider any of x=29,31,32,33,34,36,37,38,
3903: * 39,40,41).
3904: * Furthermore, at the start of the loop body we have in fact
3905: * (c') 0 <= v < v1 <= bmax, d > d1 > amax >= 0,
3906: * (and are never done already).
3907: *
3908: * Main loop is similar to those of invmod() and bezout(), except for
3909: * having to determine appropriate vmax bounds, and checking termination
3910: * conditions. The signe(d1) condition is only for paranoia
3911: */
3912: while (lgefint(d) > 3 && signe(d1))
3913: {
3914: /* determine vmax for lgcdii so as to ensure v won't overshoot.
3915: * If v+v1 > bmax, the next step would take v1 beyond the limit, so
3916: * since [+-d1,v1] is not a solution, we give up. Otherwise if v+v1
3917: * is way shorter than bmax, use vmax=MAXULUNG. Otherwise, set vmax
3918: * to a crude lower approximation of bmax/(v+v1), or to 1, which will
3919: * allow the inner loop to do one step
3920: */
3921: r = addii(v,v1);
3922: lr = lb - lgefint(r);
3923: lbr = bfffo(r[2]);
3924: if (cmpii(r,bmax) > 0) /* done, not found */
3925: {
3926: avma = av;
3927: return 0;
3928: }
3929: else if (lr > 1) /* still more than a word's worth to go */
3930: {
3931: vmax = MAXULONG;
3932: }
3933: else /* take difference of bit lengths */
3934: {
3935: lr = (lr << TWOPOTBITS_IN_LONG) - lbb + lbr;
3936: if (lr > BITS_IN_LONG)
3937: vmax = MAXULONG;
3938: else if (lr == 0)
3939: vmax = 1UL;
3940: else
3941: vmax = 1UL << (lr-1);
3942: /* the latter is pessimistic but faster than a division */
3943: }
3944: /* do a Lehmer-Jebelean round */
3945: lhmres = lgcdii((ulong *)d, (ulong *)d1, &xu, &xu1, &xv, &xv1, vmax);
3946: if (lhmres != 0) /* check progress */
3947: { /* apply matrix */
3948: if ((lhmres == 1) || (lhmres == -1))
3949: {
3950: s = -s;
3951: if (xv1 == 1)
3952: {
3953: /* re-use v+v1 computed above */
3954: v=v1; v1=r;
3955: r = subii(d,d1); d=d1; d1=r;
3956: }
3957: else
3958: {
3959: r = subii(d, mului(xv1,d1)); d=d1; d1=r;
3960: r = addii(v, mului(xv1,v1)); v=v1; v1=r;
3961: }
3962: }
3963: else
3964: {
3965: r = subii(muliu(d,xu), muliu(d1,xv));
3966: d1 = subii(muliu(d,xu1), muliu(d1,xv1)); d = r;
3967: r = addii(muliu(v,xu), muliu(v1,xv));
3968: v1 = addii(muliu(v,xu1), muliu(v1,xv1)); v = r;
3969: if (lhmres&1)
3970: {
3971: setsigne(d,-signe(d));
3972: s = -s;
3973: }
3974: else if (signe(d1))
3975: {
3976: setsigne(d1,-signe(d1));
3977: }
3978: }
3979: /* check whether we're done. Assert v <= bmax here. Examine v1:
3980: * if v1 > bmax, check d and return 0 or 1 depending on the outcome;
3981: * if v1 <= bmax, check d1 and return 1 if d1 <= amax, otherwise
3982: * proceed.
3983: */
3984: if (cmpii(v1,bmax) > 0) /* certainly done */
3985: {
3986: avma = av;
3987: if (cmpii(d,amax) <= 0) /* done, found */
3988: {
3989: if (a != NULL)
3990: {
3991: *a = icopy(d);
3992: setsigne(*a,-s); /* sign opposite to s */
3993: }
3994: if (b != NULL) *b = icopy(v);
3995: return 1;
3996: }
3997: else /* done, not found */
3998: return 0;
3999: }
4000: else if (cmpii(d1,amax) <= 0) /* also done, found */
4001: {
4002: avma = av;
4003: if (a != NULL)
4004: {
4005: if (signe(d1))
4006: {
4007: *a = icopy(d1);
4008: setsigne(*a,s); /* same sign as s */
4009: }
4010: else
4011: *a = gzero;
4012: }
4013: if (b != NULL) *b = icopy(v1);
4014: return 1;
4015: }
4016: } /* lhmres != 0 */
4017:
4018: if (lhmres <= 0 && signe(d1))
4019: {
4020: q = dvmdii(d,d1,&r);
4021: #ifdef DEBUG_LEHMER
4022: fprintferr("Full division:\n");
4023: printf(" q = "); output(q); sleep (1);
4024: #endif
4025: d=d1; d1=r;
4026: r = addii(v,mulii(q,v1));
4027: v=v1; v1=r;
4028: s = -s;
4029: /* check whether we are done now. Since we weren't before the div, it
4030: * suffices to examine v1 and d1 -- the new d (former d1) cannot cut it
4031: */
4032: if (cmpii(v1,bmax) > 0) /* done, not found */
4033: {
4034: avma = av;
4035: return 0;
4036: }
4037: else if (cmpii(d1,amax) <= 0) /* done, found */
4038: {
4039: avma = av;
4040: if (a != NULL)
4041: {
4042: if (signe(d1))
4043: {
4044: *a = icopy(d1);
4045: setsigne(*a,s); /* same sign as s */
4046: }
4047: else
4048: *a = gzero;
4049: }
4050: if (b != NULL) *b = icopy(v1);
4051: return 1;
4052: }
4053: }
4054:
4055: if (low_stack(lim, stack_lim(av,1)))
4056: {
4057: GEN *gptr[4]; gptr[0]=&d; gptr[1]=&d1; gptr[2]=&v; gptr[3]=&v1;
4058: if(DEBUGMEM>1) err(warnmem,"ratlift");
4059: gerepilemany(av1,gptr,4);
4060: }
4061: } /* end while */
4062:
4063: /* Postprocessing - final sprint. Since we usually underestimate vmax,
4064: * this function needs a loop here instead of a simple conditional.
4065: * Note we can only get here when amax fits into one word (which will
4066: * typically not be the case!). The condition is bogus -- d1 is never
4067: * zero at the start of the loop. There will be at most a few iterations,
4068: * so we don't bother collecting garbage
4069: */
4070: while (signe(d1))
4071: {
4072: /* Assertions: lgefint(d)==lgefint(d1)==3.
4073: * Moreover, we aren't done already, or we would have returned by now.
4074: * Recompute vmax...
4075: */
4076: #ifdef DEBUG_RATLIFT
4077: fprintferr("rl-fs: d,d1=%Z,%Z\n", d, d1);
4078: fprintferr("rl-fs: v,v1=%Z,%Z\n", v, v1);
4079: #endif
4080: r = addii(v,v1);
4081: lr = lb - lgefint(r);
4082: lbr = bfffo(r[2]);
4083: if (cmpii(r,bmax) > 0) /* done, not found */
4084: {
4085: avma = av;
4086: return 0;
4087: }
4088: else if (lr > 1) /* still more than a word's worth to go */
4089: {
4090: vmax = MAXULONG; /* (cannot in fact happen) */
4091: }
4092: else /* take difference of bit lengths */
4093: {
4094: lr = (lr << TWOPOTBITS_IN_LONG) - lbb + lbr;
4095: if (lr > BITS_IN_LONG)
4096: vmax = MAXULONG;
4097: else if (lr == 0)
4098: vmax = 1UL;
4099: else
4100: vmax = 1UL << (lr-1); /* as above */
4101: }
4102: #ifdef DEBUG_RATLIFT
4103: fprintferr("rl-fs: vmax=%lu\n", vmax);
4104: #endif
4105: /* single-word "Lehmer", discarding the gcd or whatever it returns */
4106: (void)rgcduu((ulong)d[2], (ulong)d1[2], vmax, &xu, &xu1, &xv, &xv1, &s0);
4107: #ifdef DEBUG_RATLIFT
4108: fprintferr("rl-fs: [%lu,%lu; %lu,%lu] %s\n",
4109: xu, xu1, xv, xv1,
4110: s0 < 0 ? "-" : "+");
4111: #endif
4112: if (xv1 == 1) /* avoid multiplications */
4113: {
4114: /* re-use v+v1 computed above */
4115: v=v1; v1=r;
4116: r = subii(d,d1); d=d1; d1=r;
4117: s = -s;
4118: }
4119: else if (xu == 0) /* and xv==1, xu1==1, xv1 > 1 */
4120: {
4121: r = subii(d, mului(xv1,d1)); d=d1; d1=r;
4122: r = addii(v, mului(xv1,v1)); v=v1; v1=r;
4123: s = -s;
4124: }
4125: else
4126: {
4127: r = subii(muliu(d,xu), muliu(d1,xv));
4128: d1 = subii(muliu(d,xu1), muliu(d1,xv1)); d = r;
4129: r = addii(muliu(v,xu), muliu(v1,xv));
4130: v1 = addii(muliu(v,xu1), muliu(v1,xv1)); v = r;
4131: if (s0 < 0)
4132: {
4133: setsigne(d,-signe(d));
4134: s = -s;
4135: }
4136: else if (signe(d1)) /* sic: might vanish now */
4137: {
4138: setsigne(d1,-signe(d1));
4139: }
4140: }
4141: /* check whether we're done, as above. Assert v <= bmax. Examine v1:
4142: * if v1 > bmax, check d and return 0 or 1 depending on the outcome;
4143: * if v1 <= bmax, check d1 and return 1 if d1 <= amax, otherwise proceed.
4144: */
4145: if (cmpii(v1,bmax) > 0) /* certainly done */
4146: {
4147: avma = av;
4148: if (cmpii(d,amax) <= 0) /* done, found */
4149: {
4150: if (a != NULL)
4151: {
4152: *a = icopy(d);
4153: setsigne(*a,-s); /* sign opposite to s */
4154: }
4155: if (b != NULL) *b = icopy(v);
4156: return 1;
4157: }
4158: else /* done, not found */
4159: return 0;
4160: }
4161: else if (cmpii(d1,amax) <= 0) /* also done, found */
4162: {
4163: avma = av;
4164: if (a != NULL)
4165: {
4166: if (signe(d1))
4167: {
4168: *a = icopy(d1);
4169: setsigne(*a,s); /* same sign as s */
4170: }
4171: else
4172: *a = gzero;
4173: }
4174: if (b != NULL) *b = icopy(v1);
4175: return 1;
4176: }
4177: } /* while */
4178:
4179: /* get here when we have run into d1 == 0 before returning... in fact,
4180: * this cannot happen.
4181: */
4182: err(talker, "ratlift failed to catch d1 == 0\n");
4183: /* NOTREACHED */
4184: return 0;
4185: }
4186:
4187: /********************************************************************/
4188: /** **/
4189: /** INTEGER SQUARE ROOT **/
4190: /** **/
4191: /********************************************************************/
4192:
4193: /* sqrt()'s result may be off by 1 when a is not representable exactly as a
4194: * double [64bit machine] */
4195: ulong
4196: usqrtsafe(ulong a)
4197: {
4198: ulong x = (ulong)sqrt((double)a);
4199: #ifdef LONG_IS_64BIT
4200: ulong y = x+1; if (y <= MAXHALFULONG && y*y <= a) x = y;
4201: #endif
4202: return x;
4203: }
4204:
4205: /* Assume a >= 0 has <= 4 words, return trunc[sqrt(a)] */
4206: ulong
4207: mpsqrtl(GEN a)
4208: {
4209: long l = lgefint(a);
4210: ulong x,y,z,k,m;
4211: int L, s;
4212:
4213: if (l <= 3) return l==2? 0: usqrtsafe((ulong)a[2]);
4214: s = bfffo(a[2]);
4215: if (s > 1)
4216: {
4217: s &= ~1UL; /* make it even */
4218: z = BITS_IN_LONG - s;
4219: m = (ulong)(a[2] << s) | (a[3] >> z);
4220: L = z>>1;
4221: }
4222: else
4223: {
4224: m = (ulong)a[2];
4225: L = BITS_IN_LONG/2;
4226: }
4227: /* m = BIL-1 (bffo odd) or BIL first bits of a */
4228: k = (long) sqrt((double)m);
4229: if (L == BITS_IN_LONG/2 && k == MAXHALFULONG)
4230: x = MAXULONG;
4231: else
4232: x = (k+1) << L;
4233: do
4234: {
4235: LOCAL_HIREMAINDER;
4236: LOCAL_OVERFLOW;
4237: hiremainder = a[2];
4238: if (hiremainder >= x) return x;
4239: z = addll(x, divll(a[3],x));
4240: z >>= 1; if (overflow) z |= HIGHBIT;
4241: y = x; x = z;
4242: }
4243: while (x < y);
4244: return y;
4245: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to mp.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari-2.2 / src / kernel / none