Annotation of OpenXM_contrib/pari-2.2/src/basemath/arith1.c, Revision 1.1.1.1
1.1 noro 1: /* $Id: arith1.c,v 1.54 2001/10/01 12:11:28 karim 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: /** ARITHMETIC FUNCTIONS **/
19: /** (first part) **/
20: /** **/
21: /*********************************************************************/
22: #include "pari.h"
23:
24: /*********************************************************************/
25: /** **/
26: /** ARITHMETIC FUNCTION PROTOTYPES **/
27: /** **/
28: /*********************************************************************/
29: GEN
30: garith_proto(GEN f(GEN), GEN x, int do_error)
31: {
32: long tx=typ(x),lx,i;
33: GEN y;
34:
35: if (is_matvec_t(tx))
36: {
37: lx=lg(x); y=cgetg(lx,tx);
38: for (i=1; i<lx; i++) y[i] = (long) garith_proto(f,(GEN) x[i], do_error);
39: return y;
40: }
41: if (tx != t_INT && do_error) err(arither1);
42: return f(x);
43: }
44:
45: GEN
46: arith_proto(long f(GEN), GEN x, int do_error)
47: {
48: long tx=typ(x),lx,i;
49: GEN y;
50:
51: if (is_matvec_t(tx))
52: {
53: lx=lg(x); y=cgetg(lx,tx);
54: for (i=1; i<lx; i++) y[i] = (long) arith_proto(f,(GEN) x[i], do_error);
55: return y;
56: }
57: if (tx != t_INT && do_error) err(arither1);
58: return stoi(f(x));
59: }
60:
61: GEN
62: arith_proto2(long f(GEN,GEN), GEN x, GEN n)
63: {
64: long l,i,tx = typ(x);
65: GEN y;
66:
67: if (is_matvec_t(tx))
68: {
69: l=lg(x); y=cgetg(l,tx);
70: for (i=1; i<l; i++) y[i] = (long) arith_proto2(f,(GEN) x[i],n);
71: return y;
72: }
73: if (tx != t_INT) err(arither1);
74: tx=typ(n);
75: if (is_matvec_t(tx))
76: {
77: l=lg(n); y=cgetg(l,tx);
78: for (i=1; i<l; i++) y[i] = (long) arith_proto2(f,x,(GEN) n[i]);
79: return y;
80: }
81: if (tx != t_INT) err(arither1);
82: return stoi(f(x,n));
83: }
84:
85: static GEN
86: arith_proto2gs(long f(GEN,long), GEN x, long y)
87: {
88: long l,i,tx = typ(x);
89: GEN t;
90:
91: if (is_matvec_t(tx))
92: {
93: l=lg(x); t=cgetg(l,tx);
94: for (i=1; i<l; i++) t[i]= (long) arith_proto2gs(f,(GEN) x[i],y);
95: return t;
96: }
97: if (tx != t_INT) err(arither1);
98: return stoi(f(x,y));
99: }
100:
101: GEN
102: garith_proto2gs(GEN f(GEN,long), GEN x, long y)
103: {
104: long l,i,tx = typ(x);
105: GEN t;
106:
107: if (is_matvec_t(tx))
108: {
109: l=lg(x); t=cgetg(l,tx);
110: for (i=1; i<l; i++) t[i]= (long) garith_proto2gs(f,(GEN) x[i],y);
111: return t;
112: }
113: if (tx != t_INT) err(arither1);
114: return f(x,y);
115: }
116:
117: GEN
118: gassoc_proto(GEN f(GEN,GEN), GEN x, GEN y)
119: {
120: int tx=typ(x);
121: if (!y)
122: {
123: ulong av=avma;
124: if (tx!=t_VEC && tx!=t_COL)
125: err(typeer,"association");
126: return gerepileupto(av,divide_conquer_prod(x,f));
127: }
128: return f(x,y);
129: }
130:
131: /*********************************************************************/
132: /** **/
133: /** ORDER of INTEGERMOD x in (Z/nZ)* **/
134: /** **/
135: /*********************************************************************/
136:
137: GEN
138: order(GEN x)
139: {
140: long av=avma,av1,i,e;
141: GEN o,m,p;
142:
143: if (typ(x) != t_INTMOD || !gcmp1(mppgcd((GEN) x[1],(GEN) x[2])))
144: err(talker,"not an element of (Z/nZ)* in order");
145: o=phi((GEN) x[1]); m=decomp(o);
146: for (i = lg(m[1])-1; i; i--)
147: {
148: p=gcoeff(m,i,1); e=itos(gcoeff(m,i,2));
149: do
150: {
151: GEN o1=divii(o,p), y=powgi(x,o1);
152: if (!gcmp1((GEN)y[2])) break;
153: e--; o = o1;
154: }
155: while (e);
156: }
157: av1=avma; return gerepile(av,av1,icopy(o));
158: }
159:
160: /******************************************************************/
161: /* */
162: /* GENERATOR of (Z/mZ)* */
163: /* */
164: /******************************************************************/
165:
166: GEN
167: ggener(GEN m)
168: {
169: return garith_proto(gener,m,1);
170: }
171:
172: GEN
173: gener(GEN m)
174: {
175: long av=avma,av1,k,i,e;
176: GEN x,t,q,p;
177:
178: if (typ(m) != t_INT) err(arither1);
179: e = signe(m);
180: if (!e) err(talker,"zero modulus in znprimroot");
181: if (is_pm1(m)) { avma=av; return gmodulss(0,1); }
182: if (e<0) m = absi(m);
183:
184: e = mod4(m);
185: if (e == 0) /* m = 0 mod 4 */
186: {
187: if (cmpis(m,4)) err(generer); /* m != 4, non cyclic */
188: return gmodulsg(3,m);
189: }
190: if (e == 2) /* m = 0 mod 2 */
191: {
192: q=shifti(m,-1); x = (GEN) gener(q)[2];
193: if (!mod2(x)) x = addii(x,q);
194: av1=avma; return gerepile(av,av1,gmodulcp(x,m));
195: }
196:
197: t=decomp(m); if (lg(t[1]) != 2) err(generer);
198: p=gcoeff(t,1,1); e=itos(gcoeff(t,1,2)); q=subis(p,1);
199: if (e>=2)
200: {
201: x = (GEN) gener(p)[2];
202: if (gcmp1(powmodulo(x,q, sqri(p)))) x = addii(x,p);
203: av1=avma; return gerepile(av,av1,gmodulcp(x,m));
204: }
205:
206: p=(GEN)decomp(q)[1]; k=lg(p)-1; x=stoi(1);
207: for(;;)
208: {
209: x[2]++;
210: if (gcmp1(mppgcd(m,x)))
211: {
212: for (i=k; i; i--)
213: if (gcmp1(powmodulo(x, divii(q,(GEN)p[i]), m))) break;
214: if (!i) break;
215: }
216: }
217: av1=avma; return gerepile(av,av1,gmodulcp(x,m));
218: }
219:
220: GEN
221: znstar(GEN n)
222: {
223: GEN p1,z,q,u,v,d,list,ep,h,gen,moduli,p,a;
224: long i,j,nbp,sizeh;
225: ulong av;
226:
227: if (typ(n) != t_INT) err(arither1);
228: if (!signe(n))
229: {
230: z=cgetg(4,t_VEC);
231: z[1]=deux; p1=cgetg(2,t_VEC);
232: z[2]=(long)p1; p1[1]=deux; p1=cgetg(2,t_VEC);
233: z[3]=(long)p1; p1[1]=lneg(gun);
234: return z;
235: }
236: av=avma; n=absi(n);
237: if (cmpis(n,2)<=0)
238: {
239: avma=av; z=cgetg(4,t_VEC);
240: z[1]=un;
241: z[2]=lgetg(1,t_VEC);
242: z[3]=lgetg(1,t_VEC);
243: return z;
244: }
245: list=factor(n); ep=(GEN)list[2]; list=(GEN)list[1]; nbp=lg(list)-1;
246: h = cgetg(nbp+2,t_VEC);
247: gen = cgetg(nbp+2,t_VEC);
248: moduli = cgetg(nbp+2,t_VEC);
249: switch(mod8(n))
250: {
251: case 0:
252: h[1] = lmul2n(gun,itos((GEN)ep[1])-2); h[2] = deux;
253: gen[1] = lstoi(5); gen[2] = laddis(gmul2n((GEN)h[1],1),-1);
254: moduli[1] = moduli[2] = lmul2n(gun,itos((GEN)ep[1]));
255: sizeh = nbp+1; i=3; j=2; break;
256: case 4:
257: h[1] = deux;
258: gen[1] = lstoi(3);
259: moduli[1] = lmul2n(gun,itos((GEN)ep[1]));
260: sizeh = nbp; i=j=2; break;
261: case 2: case 6:
262: sizeh = nbp-1; i=1; j=2; break;
263: default: /* 1, 3, 5, 7 */
264: sizeh = nbp; i=j=1;
265: }
266: for ( ; j<=nbp; i++,j++)
267: {
268: p = (GEN)list[j]; q = gpuigs(p, itos((GEN)ep[j])-1);
269: h[i] = lmulii(addis(p,-1),q); p1 = mulii(p,q);
270: gen[i] = gener(p1)[2];
271: moduli[i] = (long)p1;
272: }
273: #if 0
274: if (nbp == 1 && is_pm1(ep[1]))
275: gen[1] = lmodulcp((GEN)gen[1],n);
276: else
277: #endif
278: for (i=1; i<=sizeh; i++)
279: {
280: q = (GEN)moduli[i]; a = (GEN)gen[i];
281: u = mpinvmod(q,divii(n,q));
282: gen[i] = lmodulcp(addii(a,mulii(mulii(subsi(1,a),u),q)), n);
283: }
284:
285: for (i=sizeh; i>=2; i--)
286: for (j=i-1; j>=1; j--)
287: if (!divise((GEN)h[j],(GEN)h[i]))
288: {
289: d=bezout((GEN)h[i],(GEN)h[j],&u,&v);
290: q=divii((GEN)h[j],d);
291: h[j]=(long)mulii((GEN)h[i],q); h[i]=(long)d;
292: gen[j]=ldiv((GEN)gen[j], (GEN)gen[i]);
293: gen[i]=lmul((GEN)gen[i], powgi((GEN)gen[j], mulii(v,q)));
294: }
295: q=gun; for (i=1; i<=sizeh && !gcmp1((GEN)h[i]); i++) q=mulii(q,(GEN)h[i]);
296: setlg(h,i); setlg(gen,i); z=cgetg(4,t_VEC);
297: z[1]=(long)q; z[2]=(long)h; z[3]=(long)gen;
298: return gerepilecopy(av,z);
299: }
300:
301: /*********************************************************************/
302: /** **/
303: /** INTEGRAL SQUARE ROOT **/
304: /** **/
305: /*********************************************************************/
306: extern ulong mpsqrtl(GEN a);
307:
308: GEN
309: gracine(GEN a)
310: {
311: return garith_proto(racine,a,1); /* hm. --GN */
312: }
313:
314: /* Use l as lgefint(a) [a may have more digits] */
315: static GEN
316: racine_r(GEN a, long l)
317: {
318: long av,s;
319: GEN x,y,z;
320:
321: if (l <= 4) return utoi(mpsqrtl(a));
322: av = avma;
323: s = 2 + ((l-1) >> 1);
324: setlgefint(a, s);
325: x = addis(racine_r(a, s), 1); setlgefint(a, l);
326: /* x = good approx (from above) of sqrt(a): about l/2 correct bits */
327: x = shifti(x, (l - s)*(BITS_IN_LONG/2));
328: do
329: { /* one or two iterations should do the trick */
330: z = shifti(addii(x,divii(a,x)), -1);
331: y = x; x = z;
332: }
333: while (cmpii(x,y) < 0);
334: avma = (long)y;
335: return gerepileuptoint(av,y);
336: }
337:
338: GEN
339: racine_i(GEN a, int roundup)
340: {
341: long k,m,l = lgefint(a), av = avma;
342: GEN x = racine_r(a, l);
343: if (roundup && signe(x))
344: {
345: m = modBIL(x);
346: k = (m * m != a[l-1] || !egalii(sqri(x),a));
347: avma = (long)x;
348: if (k) x = gerepileuptoint(av, addis(x,1));
349: }
350: return x;
351: }
352:
353: GEN
354: racine(GEN a)
355: {
356: if (typ(a) != t_INT) err(arither1);
357: switch (signe(a))
358: {
359: case 1: return racine_i(a,0);
360: case 0: return gzero;
361: case -1:
362: {
363: GEN x = cgetg(3,t_COMPLEX);
364: x[1] = zero;
365: x[2] = (long)racine_i(a,0); return x;
366: }
367: }
368: return NULL; /* not reached */
369: }
370:
371: /*********************************************************************/
372: /** **/
373: /** PERFECT SQUARE **/
374: /** **/
375: /*********************************************************************/
376: extern ulong usqrtsafe(ulong a);
377:
378: static int
379: carremod(ulong A)
380: {
381: static int carresmod64[]={
382: 1,1,0,0,1,0,0,0,0,1, 0,0,0,0,0,0,1,1,0,0, 0,0,0,0,0,1,0,0,0,0,
383: 0,0,0,1,0,0,1,0,0,0, 0,1,0,0,0,0,0,0,0,1, 0,0,0,0,0,0,0,1,0,0, 0,0,0,0};
384: static int carresmod63[]={
385: 1,1,0,0,1,0,0,1,0,1, 0,0,0,0,0,0,1,0,1,0, 0,0,1,0,0,1,0,0,1,0,
386: 0,0,0,0,0,0,1,1,0,0, 0,0,0,1,0,0,1,0,0,1, 0,0,0,0,0,0,0,0,1,0, 0,0,0};
387: static int carresmod65[]={
388: 1,1,0,0,1,0,0,0,0,1, 1,0,0,0,1,0,1,0,0,0, 0,0,0,0,0,1,1,0,0,1,
389: 1,0,0,0,0,1,1,0,0,1, 1,0,0,0,0,0,0,0,0,1, 0,1,0,0,0,1,1,0,0,0, 0,1,0,0,1};
390: static int carresmod11[]={1,1,0,1,1,1,0,0,0,1, 0};
391: return (carresmod64[A & 0x3fUL]
392: && carresmod63[A % 63UL]
393: && carresmod65[A % 65UL]
394: && carresmod11[A % 11UL]);
395: }
396:
397: /* emulate carrecomplet on single-word positive integers */
398: ulong
399: ucarrecomplet(ulong A)
400: {
401: if (carremod(A))
402: {
403: ulong a = usqrtsafe(A);
404: if (a * a == A) return a;
405: }
406: return 0;
407: }
408:
409: long
410: carrecomplet(GEN x, GEN *pt)
411: {
412: long av;
413: GEN y;
414:
415: switch(signe(x))
416: {
417: case -1: return 0;
418: case 0: if (pt) *pt=gzero; return 1;
419: }
420: if (lgefint(x) == 3)
421: {
422: long a = ucarrecomplet((ulong)x[2]);
423: if (!a) return 0;
424: if (pt) *pt = stoi(a);
425: return 1;
426: }
427: if (!carremod((ulong)smodis(x, 64*63*65*11))) return 0;
428: av=avma; y = racine(x);
429: if (!egalii(sqri(y),x)) { avma=av; return 0; }
430: if (pt) { *pt = y; avma=(long)y; } else avma=av;
431: return 1;
432: }
433:
434: static GEN
435: polcarrecomplet(GEN x, GEN *pt)
436: {
437: long i,l,av,av2;
438: GEN y,a,b;
439:
440: if (!signe(x)) return gun;
441: l=lgef(x); if ((l&1) == 0) return gzero; /* odd degree */
442: i=2; while (isexactzero((GEN)x[i])) i++;
443: if (i&1) return gzero;
444: av2 = avma; a = (GEN)x[i];
445: switch (typ(a))
446: {
447: case t_POL: case t_INT:
448: y = gcarrecomplet(a,&b); break;
449: default:
450: y = gcarreparfait(a); b = NULL; break;
451: }
452: if (y == gzero) { avma = av2; return gzero; }
453: av = avma; x = gdiv(x,a);
454: y = gtrunc(gsqrt(greffe(x,l,1),0)); av2 = avma;
455: if (!gegal(gsqr(y), x)) { avma=av; return gzero; }
456: if (pt)
457: {
458: avma = av2;
459: if (!gcmp1(a))
460: {
461: if (!b) b = gsqrt(a,DEFAULTPREC);
462: y = gmul(b,y);
463: }
464: *pt = gerepileupto(av,y);
465: }
466: else avma = av;
467: return gun;
468: }
469:
470: GEN
471: gcarrecomplet(GEN x, GEN *pt)
472: {
473: long tx = typ(x),l,i;
474: GEN z,y,p,t;
475:
476: if (!pt) return gcarreparfait(x);
477: if (is_matvec_t(tx))
478: {
479: l=lg(x); y=cgetg(l,tx); z=cgetg(l,tx);
480: for (i=1; i<l; i++)
481: {
482: t=gcarrecomplet((GEN)x[i],&p);
483: y[i] = (long)t;
484: z[i] = gcmp0(t)? zero: (long)p;
485: }
486: *pt=z; return y;
487: }
488: if (tx == t_POL) return polcarrecomplet(x,pt);
489: if (tx != t_INT) err(arither1);
490: return stoi(carrecomplet(x,pt));
491: }
492:
493: GEN
494: gcarreparfait(GEN x)
495: {
496: GEN p1,a,p;
497: long tx=typ(x),l,i,av,v;
498:
499: switch(tx)
500: {
501: case t_INT:
502: return carreparfait(x)? gun: gzero;
503:
504: case t_REAL:
505: return (signe(x)>=0)? gun: gzero;
506:
507: case t_INTMOD:
508: {
509: GEN b, q, e;
510: long w;
511: a = (GEN)x[2]; if (!signe(a)) return gun;
512: av = avma;
513: q = absi((GEN)x[1]); v = vali(q);
514: if (v) /* > 0 */
515: {
516: long dv;
517: w = vali(a); dv = v - w;
518: if (dv > 0)
519: {
520: if (w & 1) { avma = av; return gzero; }
521: if (dv >= 2)
522: {
523: b = w? shifti(a,-w): a;
524: if ((dv>=3 && mod8(b) != 1) ||
525: (dv==2 && mod4(b) != 1)) { avma = av; return gzero; }
526: }
527: }
528: q = shifti(q, -v);
529: }
530: /* q is now odd */
531: i = kronecker(a,q);
532: if (i < 0) { avma = av; return gzero; }
533: if (i==0)
534: {
535: GEN d = mppgcd(a,q);
536: p1 = factor(d);
537: p = (GEN)p1[1];
538: e = (GEN)p1[2]; l = lg(e);
539: for (i=1; i<l; i++)
540: {
541: v = pvaluation(a,(GEN)p[i],&p1);
542: w = itos((GEN)e[i]);
543: if (v < w && (v&1 || kronecker(p1,(GEN)p[i]) == -1))
544: { avma = av; return gzero; }
545: q = diviiexact(q, gpowgs((GEN)p[i], w));
546: }
547: if (kronecker(a,q) == -1) { avma = av; return gzero; }
548: }
549: /* kro(a,q) = 1, q odd: need to factor q */
550: p1 = factor(q);
551: p = (GEN)p1[1];
552: e = (GEN)p1[2]; l = lg(e) - 1;
553: /* kro(a,q) = 1, check all p|q but the last (product formula) */
554: for (i=1; i<l; i++)
555: if (kronecker(a,(GEN)p[i]) == -1) { avma = av; return gzero; }
556: return gun;
557: }
558:
559: case t_FRAC: case t_FRACN:
560: av=avma; l=carreparfait(mulii((GEN)x[1],(GEN)x[2]));
561: avma=av; return l? gun: gzero;
562:
563: case t_COMPLEX:
564: return gun;
565:
566: case t_PADIC:
567: a = (GEN)x[4]; if (!signe(a)) return gun;
568: if (valp(x)&1) return gzero;
569: p = (GEN)x[2];
570: if (!egalii(p, gdeux))
571: return (kronecker(a,p)== -1)? gzero: gun;
572:
573: v = precp(x); /* here p=2, a is odd */
574: if ((v>=3 && mod8(a) != 1 ) ||
575: (v==2 && mod4(a) != 1)) return gzero;
576: return gun;
577:
578: case t_POL:
579: return polcarrecomplet(x,NULL);
580:
581: case t_SER:
582: if (!signe(x)) return gun;
583: if (valp(x)&1) return gzero;
584: return gcarreparfait((GEN)x[2]);
585:
586: case t_RFRAC: case t_RFRACN:
587: av=avma;
588: l=itos(gcarreparfait(gmul((GEN)x[1],(GEN)x[2])));
589: avma=av; return stoi(l);
590:
591: case t_QFR: case t_QFI:
592: return gcarreparfait((GEN)x[1]);
593:
594: case t_VEC: case t_COL: case t_MAT:
595: l=lg(x); p1=cgetg(l,tx);
596: for (i=1; i<l; i++) p1[i]=(long)gcarreparfait((GEN)x[i]);
597: return p1;
598: }
599: err(typeer,"issquare");
600: return NULL; /* not reached */
601: }
602:
603: /*********************************************************************/
604: /** **/
605: /** KRONECKER SYMBOL **/
606: /** **/
607: /*********************************************************************/
608: /* u = 3,5 mod 8 ? (= 2 not a square mod u) */
609: #define ome(t) (labs(((t)&7)-4) == 1)
610: #define gome(t) (ome(modBIL(t)))
611:
612: GEN
613: gkronecker(GEN x, GEN y)
614: {
615: return arith_proto2(kronecker,x,y);
616: }
617:
618: long
619: kronecker(GEN x, GEN y)
620: {
621: GEN z;
622: long av,r,s=1;
623:
624: av=avma;
625: switch (signe(y))
626: {
627: case -1: y=negi(y); if (signe(x)<0) s= -1; break;
628: case 0: return is_pm1(x);
629: }
630: r=vali(y);
631: if (r)
632: {
633: if (mpodd(x))
634: {
635: if (odd(r) && gome(x)) s= -s;
636: y=shifti(y,-r);
637: }
638: else { avma=av; return 0; }
639: }
640: x=modii(x,y);
641: while (signe(x))
642: {
643: r=vali(x);
644: if (r)
645: {
646: if (odd(r) && gome(y)) s= -s;
647: x=shifti(x,-r);
648: }
649: /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
650: if (modBIL(x) & modBIL(y) & 2) s = -s;
651: z=resii(y,x); y=x; x=z;
652: }
653: avma=av; return is_pm1(y)? s: 0;
654: }
655:
656: GEN
657: gkrogs(GEN x, long y)
658: {
659: return arith_proto2gs(krogs,x,y);
660: }
661:
662: long
663: krogs(GEN x, long y)
664: {
665: long av,r,s=1,x1,z;
666:
667: av=avma;
668: if (y<=0)
669: {
670: if (y) { y = -y; if (signe(x)<0) s = -1; }
671: else return is_pm1(x);
672: }
673: r=vals(y);
674: if (r)
675: {
676: if (mpodd(x))
677: {
678: if (odd(r) && gome(x)) s= -s;
679: y>>=r;
680: }
681: else return 0;
682: }
683: x1=smodis(x,y); avma = av;
684: while (x1)
685: {
686: r=vals(x1);
687: if (r)
688: {
689: if (odd(r) && ome(y)) s= -s;
690: x1>>=r;
691: }
692: if (x1 & y & 2) s= -s;
693: z=y%x1; y=x1; x1=z;
694: }
695: return (y==1)? s: 0;
696: }
697:
698: long
699: krosg(long s, GEN x)
700: {
701: long av = avma, y = kronecker(stoi(s),x);
702: avma = av; return y;
703: }
704:
705: long
706: kross(long x, long y)
707: {
708: long r,s=1,x1,z;
709:
710: if (y<=0)
711: {
712: if (y) { y= -y; if (x<0) s = -1; }
713: else return (labs(x)==1);
714: }
715: r=vals(y);
716: if (r)
717: {
718: if (odd(x))
719: {
720: if (odd(r) && ome(x)) s = -s;
721: y>>=r;
722: }
723: else return 0;
724: }
725: x1=x%y; if (x1<0) x1+=y;
726: while (x1)
727: {
728: r=vals(x1);
729: if (r)
730: {
731: if (odd(r) && ome(y)) s = -s;
732: x1>>=r;
733: }
734: if (x1 & y & 2) s = -s;
735: z=y%x1; y=x1; x1=z;
736: }
737: return (y==1)? s: 0;
738: }
739:
740: long
741: hil0(GEN x, GEN y, GEN p)
742: {
743: return hil(x,y, p? p: gzero);
744: }
745:
746: #define eps(t) (((signe(t)*(modBIL(t)))&3)==3)
747: long
748: hil(GEN x, GEN y, GEN p)
749: {
750: long a,b,av,tx,ty,z;
751: GEN p1,p2,u,v;
752:
753: if (gcmp0(x) || gcmp0(y)) return 0;
754: av=avma; tx=typ(x); ty=typ(y);
755: if (tx>ty) { p1=x; x=y; y=p1; a=tx,tx=ty; ty=a; }
756: switch(tx)
757: {
758: case t_INT:
759: switch(ty)
760: {
761: case t_INT:
762: if (signe(p)<=0)
763: return (signe(x)<0 && signe(y)<0)? -1: 1;
764: a = odd(pvaluation(x,p,&u));
765: b = odd(pvaluation(y,p,&v));
766: if (egalii(p,gdeux))
767: {
768: z = (eps(u) && eps(v))? -1: 1;
769: if (a && gome(v)) z= -z;
770: if (b && gome(u)) z= -z;
771: }
772: else
773: {
774: z = (a && b && eps(p))? -1: 1;
775: if (a && kronecker(v,p)<0) z= -z;
776: if (b && kronecker(u,p)<0) z= -z;
777: }
778: avma=av; return z;
779: case t_REAL:
780: return (signe(x)<0 && signe(y)<0)? -1: 1;
781: case t_INTMOD:
782: if (egalii(gdeux,(GEN)y[1])) err(hiler1);
783: return hil(x,(GEN)y[2],(GEN)y[1]);
784: case t_FRAC: case t_FRACN:
785: p1=mulii((GEN)y[1],(GEN)y[2]); z=hil(x,p1,p);
786: avma=av; return z;
787: case t_PADIC:
788: if (egalii(gdeux,(GEN)y[2]) && precp(y) <= 2) err(hiler1);
789: p1 = odd(valp(y))? mulii((GEN)y[2],(GEN)y[4]): (GEN)y[4];
790: z=hil(x,p1,(GEN)y[2]); avma=av; return z;
791: }
792: break;
793:
794: case t_REAL:
795: if (! is_frac_t(ty)) break;
796: if (signe(x) > 0) return 1;
797: return signe(y[1])*signe(y[2]);
798:
799: case t_INTMOD:
800: if (egalii(gdeux,(GEN)y[1])) err(hiler1);
801: switch(ty)
802: {
803: case t_INTMOD:
804: if (!egalii((GEN)x[1],(GEN)y[1])) break;
805: return hil((GEN)x[2],(GEN)y[2],(GEN)x[1]);
806: case t_FRAC: case t_FRACN:
807: return hil((GEN)x[2],y,(GEN)x[1]);
808: case t_PADIC:
809: if (!egalii((GEN)x[1],(GEN)y[2])) break;
810: return hil((GEN)x[2],y,(GEN)x[1]);
811: }
812: break;
813:
814: case t_FRAC: case t_FRACN:
815: p1=mulii((GEN)x[1],(GEN)x[2]);
816: switch(ty)
817: {
818: case t_FRAC: case t_FRACN:
819: p2=mulii((GEN)y[1],(GEN)y[2]);
820: z=hil(p1,p2,p); avma=av; return z;
821: case t_PADIC:
822: z=hil(p1,y,NULL); avma=av; return z;
823: }
824: break;
825:
826: case t_PADIC:
827: if (ty != t_PADIC || !egalii((GEN)x[2],(GEN)y[2])) break;
828: p1 = odd(valp(x))? mulii((GEN)x[2],(GEN)x[4]): (GEN)x[4];
829: p2 = odd(valp(y))? mulii((GEN)y[2],(GEN)y[4]): (GEN)y[4];
830: z=hil(p1,p2,(GEN)x[2]); avma=av; return z;
831: }
832: err(talker,"forbidden or incompatible types in hil");
833: return 0; /* not reached */
834: }
835: #undef eps
836: #undef ome
837: #undef gome
838:
839: /*******************************************************************/
840: /* */
841: /* SQUARE ROOT MODULO p */
842: /* */
843: /*******************************************************************/
844:
845: #define sqrmod(x,p) (resii(sqri(x),p))
846:
847: /* Assume p is prime. return NULL if (a,p) = -1 */
848: GEN
849: mpsqrtmod(GEN a, GEN p)
850: {
851: long av = avma, av1,lim,i,k,e;
852: GEN p1,q,v,y,w,m;
853:
854: if (typ(a) != t_INT || typ(p) != t_INT) err(arither1);
855: if (signe(p) <= 0 || is_pm1(p)) err(talker,"not a prime in mpsqrtmod");
856: p1 = addsi(-1,p); e = vali(p1);
857: if (e == 0) /* p = 2 */
858: {
859: avma = av;
860: if (!egalii(p, gdeux))
861: err(talker,"composite modulus in mpsqrtmod: %Z",p);
862: if (!signe(a) || !mod2(a)) return gzero;
863: return gun;
864: }
865: q = shifti(p1,-e); /* q = (p-1)/2^oo is odd */
866: if (e == 1) y = p1;
867: else /* look for an odd power of a primitive root */
868: for (k=2; ; k++)
869: { /* loop terminates for k < p (even if p composite) */
870:
871: i = krosg(k,p);
872: if (i >= 0)
873: {
874: if (i) continue;
875: err(talker,"composite modulus in mpsqrtmod: %Z",p);
876: }
877: av1 = avma;
878: y = m = powmodulo(stoi(k),q,p);
879: for (i=1; i<e; i++)
880: if (gcmp1(m = sqrmod(m,p))) break;
881: if (i == e) break; /* success */
882: avma = av1;
883: }
884:
885: p1 = powmodulo(a, shifti(q,-1), p); /* a ^ [(q-1)/2] */
886: if (!signe(p1)) { avma=av; return gzero; }
887: v = modii(mulii(a, p1), p);
888: w = modii(mulii(v, p1), p);
889: lim = stack_lim(av,1);
890: while (!gcmp1(w))
891: { /* a*w = v^2, y primitive 2^e-th root of 1
892: a square --> w even power of y, hence w^(2^(e-1)) = 1 */
893: p1 = sqrmod(w,p);
894: for (k=1; !gcmp1(p1) && k < e; k++) p1 = sqrmod(p1,p);
895: if (k == e) { avma=av; return NULL; } /* p composite or (a/p) != 1 */
896: /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
897: p1 = y;
898: for (i=1; i < e-k; i++) p1 = sqrmod(p1,p);
899: y = sqrmod(p1, p); e = k;
900: w = modii(mulii(y, w), p);
901: v = modii(mulii(v, p1), p);
902: if (low_stack(lim, stack_lim(av,1)))
903: {
904: GEN *gptr[3]; gptr[0]=&y; gptr[1]=&w; gptr[2]=&v;
905: if(DEBUGMEM>1) err(warnmem,"mpsqrtmod");
906: gerepilemany(av,gptr,3);
907: }
908: }
909: av1 = avma;
910: p1 = subii(p,v); if (cmpii(v,p1) > 0) v = p1; else avma = av1;
911: return gerepileuptoint(av, v);
912: }
913:
914: /*******************************************************************/
915: /* */
916: /* n-th ROOT MODULO p */
917: /* */
918: /*******************************************************************/
919:
920: /* UNCLEAN
921: * p-1=l^e*r
922: * l doit etre premier
923: * q=p-1
924: * q=(l^e)*r
925: * e>=1
926: * pgcd(r,l)=1
927: * return a non l-power residue and set zeta to a primitive l root of unity.
928: */
929:
930: static GEN
931: mplgenmod(GEN l, long e, GEN r,GEN p,GEN *zeta)
932: {
933: ulong av1;
934: GEN m,m1;
935: long k,i;
936: av1 = avma;
937: for (k=1; ; k++)
938: {
939: m1 = m = powmodulo(stoi(k+1),r,p);
940: for (i=1; i<e; i++)
941: if (gcmp1(m=powmodulo(m,l,p))) break;
942: if (i==e) break;
943: avma = av1;
944: }
945: *zeta=m;
946: return m1;
947: }
948: /* resoud x^l=a mod (p)
949: * l doit etre premier
950: * q=p-1
951: * q=(l^e)*r
952: * e>=1
953: * pgcd(r,l)=1
954: * m=y^(q/l)
955: * y n'est pas une puissance l-ieme
956: * m!=1
957: * ouf!
958: */
959: GEN Fp_shanks(GEN x,GEN g0,GEN p, GEN q);
960: static GEN
961: mpsqrtlmod(GEN a, GEN l, GEN p, GEN q,long e, GEN r, GEN y, GEN m)
962: {
963: ulong av = avma, tetpil,lim;
964: long k;
965: GEN p1,u1,u2,v,w,z,dl;
966: /* y contient un generateur de (Z/pZ)^* eleve a la puis (p-1)/(l^e) */
967: bezout(r,l,&u1,&u2);
968: v=powmodulo(a,u2,p);
969: w=powmodulo(a,modii(mulii(negi(u1),r),q),p);
970: lim = stack_lim(av,1);
971: while (!gcmp1(w))
972: {
973: /* if p is not prime, next loop will not end */
974: k=0;
975: p1=w;
976: do
977: {
978: z=p1;
979: p1=powmodulo(p1,l,p);
980: k++;
981: }while(!gcmp1(p1));
982: if (k==e) { avma=av; return NULL; }
983: dl=Fp_shanks(mpinvmod(z,p),m,p,l);
984: p1=powmodulo(y,modii(mulii(dl,gpowgs(l,e-k-1)),q),p);
985: m=powmodulo(m,dl,p);
986: e = k;
987: v = modii(mulii(p1,v),p);
988: y = powmodulo(p1,l,p);
989: w = modii(mulii(y,w),p);
990: if (low_stack(lim, stack_lim(av,1)))
991: {
992: GEN *gptr[4];
993: if(DEBUGMEM>1) err(warnmem,"mpsqrtlmod");
994: gptr[0]=&y; gptr[1]=&v; gptr[2]=&w; gptr[3]=&m;
995: gerepilemany(av,gptr,4);
996: }
997: }
998: tetpil=avma; return gerepile(av,tetpil,icopy(v));
999: }
1000: /* a and n are integers, p is prime
1001:
1002: return a solution of
1003:
1004: x^n=a mod p
1005:
1006: 1)If there is no solution return NULL and if zetan is not NULL set zetan to gzero.
1007:
1008: 2) If there is solution there are exactly m=gcd(p-1,n) of them.
1009:
1010: If zetan is not NULL, zetan is set to a primitive mth root of unity so that
1011: the set of solutions is {x*zetan^k;k=0 to m-1}
1012:
1013: If a=0 ,return 0 and if zetan is not NULL zetan is set to gun
1014: */
1015: GEN
1016: mpsqrtnmod(GEN a, GEN n, GEN p, GEN *zetan)
1017: {
1018: ulong ltop=avma,lbot=0,av1,lim;
1019: GEN m,u1,u2;
1020: GEN q,z;
1021: GEN *gptr[2];
1022: if (typ(a) != t_INT || typ(n) != t_INT || typ(p)!=t_INT)
1023: err(typeer,"mpsqrtnmod");
1024: if(!signe(n))
1025: err(talker,"1/0 exponent in mpsqrtnmod");
1026: if(gcmp1(n)) {if (zetan) *zetan=gun;return gcopy(a);}
1027: a=modii(a,p);
1028: if(gcmp0(a)) {if (zetan) *zetan=gun;avma=ltop;return gzero;}
1029: q=addsi(-1,p);
1030: m=bezout(n,q,&u1,&u2);
1031: if (zetan) z=gun;
1032: lim=stack_lim(ltop,1);
1033: if (!gcmp1(m))
1034: {
1035: GEN F=decomp(m);
1036: long i,j,e;
1037: GEN r,zeta,y,l;
1038: av1=avma;
1039: for (i = lg(F[1])-1; i; i--)
1040: {
1041: l=gcoeff(F,i,1); j=itos(gcoeff(F,i,2));
1042: e=pvaluation(q,l,&r);
1043: y=mplgenmod(l,e,r,p,&zeta);
1044: if (zetan) z=modii(mulii(z,powmodulo(y,gpowgs(l,e-j),p)),p);
1045: do
1046: {
1047: lbot=avma;
1048: if (!is_pm1(a) || signe(a)<0)
1049: a=mpsqrtlmod(a,l,p,q,e,r,y,zeta);
1050: else
1051: a=icopy(a);
1052: if (!a){avma=ltop;if (zetan) *zetan=gzero;return NULL;}
1053: j--;
1054: }while (j);
1055: if (low_stack(lim, stack_lim(ltop,1)))
1056: /* n can have lots of prime factors*/
1057: {
1058: if(DEBUGMEM>1) err(warnmem,"mpsqrtnmod");
1059: if (zetan)
1060: {
1061: z=gcopy(z);
1062: gptr[0]=&a;gptr[1]=&z;
1063: gerepilemanysp(av1,lbot,gptr,2);
1064: }
1065: else
1066: a=gerepile(av1,lbot,a);
1067: lbot=av1;
1068: }
1069: }
1070: }
1071: if (cmpii(m,n))
1072: {
1073: GEN b=modii(u1,q);
1074: lbot=avma;
1075: a=powmodulo(a,b,p);
1076: }
1077: if (zetan)
1078: {
1079: *zetan=gcopy(z);
1080: gptr[0]=&a;gptr[1]=zetan;
1081: gerepilemanysp(ltop,lbot,gptr,2);
1082: }
1083: else
1084: a=gerepile(ltop,lbot,a);
1085: return a;
1086: }
1087:
1088: /*********************************************************************/
1089: /** **/
1090: /** GCD & BEZOUT **/
1091: /** **/
1092: /*********************************************************************/
1093:
1094: /* Ultra-fast private ulong gcd for trial divisions. Called with y odd;
1095: x can be arbitrary (but will most of the time be smaller than y).
1096: Will also be used from inside ifactor2.c, so it's `semi-private' really.
1097: --GN */
1098:
1099: /* Gotos are Harmful, and Programming is a Science. E.W.Dijkstra. */
1100: ulong
1101: ugcd(ulong x, ulong y) /* assume y&1==1, y > 1 */
1102: {
1103: if (!x) return y;
1104: /* fix up x */
1105: while (!(x&1)) x>>=1;
1106: if (x==1) return 1;
1107: if (x==y) return y;
1108: else if (x>y) goto xislarger;/* will be rare, given how we'll use this */
1109: /* loop invariants: x,y odd and distinct. */
1110: yislarger:
1111: if ((x^y)&2) /* ...01, ...11 or vice versa */
1112: y=(x>>2)+(y>>2)+1; /* ==(x+y)>>2 except it can't overflow */
1113: else /* ...01,...01 or ...11,...11 */
1114: y=(y-x)>>2; /* now y!=0 in either case */
1115: while (!(y&1)) y>>=1; /* kill any windfall-gained powers of 2 */
1116: if (y==1) return 1; /* comparand == return value... */
1117: if (x==y) return y; /* this and the next is just one comparison */
1118: else if (x<y) goto yislarger;/* else fall through to xislarger */
1119:
1120: xislarger: /* same as above, seen through a mirror */
1121: if ((x^y)&2)
1122: x=(x>>2)+(y>>2)+1;
1123: else
1124: x=(x-y)>>2; /* x!=0 */
1125: while (!(x&1)) x>>=1;
1126: if (x==1) return 1;
1127: if (x==y) return y;
1128: else if (x>y) goto xislarger;/* again, a decent [g]cc should notice that
1129: it can re-use the comparison */
1130: goto yislarger;
1131: }
1132: /* Gotos are useful, and Programming is an Art. D.E.Knuth. */
1133: /* PS: Of course written with Dijkstra's lessons firmly in mind... --GN */
1134:
1135: /* modified right shift binary algorithm with at most one division */
1136: long
1137: cgcd(long a,long b)
1138: {
1139: long v;
1140: a=labs(a); if (!b) return a;
1141: b=labs(b); if (!a) return b;
1142: if (a>b) { a %= b; if (!a) return b; }
1143: else { b %= a; if (!b) return a; }
1144: v=vals(a|b); a>>=v; b>>=v;
1145: if (a==1 || b==1) return 1L<<v;
1146: if (b&1)
1147: return ((long)ugcd((ulong)a, (ulong)b)) << v;
1148: else
1149: return ((long)ugcd((ulong)b, (ulong)a)) << v;
1150: }
1151: long
1152: clcm(long a,long b)
1153: {
1154: long d;
1155: if(!a) return 0;
1156: d=cgcd(a,b);
1157: if(d!=1) return a*(b/d);
1158: return a*b;
1159: }
1160:
1161: /* assume a>b>0, return gcd(a,b) as a GEN (for mppgcd) */
1162: static GEN
1163: mppgcd_cgcd(ulong a, ulong b)
1164: {
1165: GEN r = cgeti(3);
1166: long v;
1167:
1168: r[1] = evalsigne(1)|evallgefint(3);
1169: a %= b; if (!a) { r[2]=(long)b; return r; }
1170: v=vals(a|b); a>>=v; b>>=v;
1171: if (a==1 || b==1) { r[2]=(long)(1UL<<v); return r; }
1172: if (b&1)
1173: r[2] = (long)(ugcd((ulong)a, (ulong)b) << v);
1174: else
1175: r[2] = (long)(ugcd((ulong)b, (ulong)a) << v);
1176: return r;
1177: }
1178:
1179: void mppgcd_plus_minus(GEN x, GEN y, GEN t);
1180: ulong mppgcd_resiu(GEN y, ulong x);
1181:
1182: /* uses modified right-shift binary algorithm now --GN 1998Jul23 */
1183: GEN
1184: mppgcd(GEN a, GEN b)
1185: {
1186: long av,v,w;
1187: GEN t, p1;
1188:
1189: if (typ(a) != t_INT || typ(b) != t_INT) err(arither1);
1190: switch (absi_cmp(a,b))
1191: {
1192: case 0: return absi(a);
1193: case -1: t=b; b=a; a=t;
1194: }
1195: if (!signe(b)) return absi(a);
1196: /* here |a|>|b|>0. Try single precision first */
1197: if (lgefint(a)==3)
1198: return mppgcd_cgcd((ulong)a[2], (ulong)b[2]);
1199: if (lgefint(b)==3)
1200: {
1201: ulong u = mppgcd_resiu(a,(ulong)b[2]);
1202: if (!u) return absi(b);
1203: return mppgcd_cgcd((ulong)b[2], u);
1204: }
1205:
1206: /* larger than gcd: "avma=av" gerepile (erasing t) is valid */
1207: av = avma; new_chunk(lgefint(b)); /* HACK */
1208: t = resii(a,b);
1209: if (!signe(t)) { avma=av; return absi(b); }
1210:
1211: a = b; b = t;
1212: v = vali(a); a = shifti(a,-v); setsigne(a,1);
1213: w = vali(b); b = shifti(b,-w); setsigne(b,1);
1214: if (w < v) v = w;
1215: switch(absi_cmp(a,b))
1216: {
1217: case 0: avma=av; a=shifti(a,v); return a;
1218: case -1: p1=b; b=a; a=p1;
1219: }
1220: if (is_pm1(b)) { avma=av; return shifti(gun,v); }
1221:
1222: /* we have three consecutive memory locations: a,b,t.
1223: * All computations are done in place */
1224:
1225: /* a and b are odd now, and a>b>1 */
1226: while (lgefint(a) > 3)
1227: {
1228: /* if a=b mod 4 set t=a-b, otherwise t=a+b, then strip powers of 2 */
1229: /* so that t <= (a+b)/4 < a/2 */
1230: mppgcd_plus_minus(a,b, t);
1231: if (is_pm1(t)) { avma=av; return shifti(gun,v); }
1232: switch(absi_cmp(t,b))
1233: {
1234: case -1: p1 = a; a = b; b = t; t = p1; break;
1235: case 1: p1 = a; a = t; t = p1; break;
1236: case 0: avma = av; b=shifti(b,v); return b;
1237: }
1238: }
1239: {
1240: long r[] = {evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3), 0};
1241: r[2] = (long) ugcd((ulong)b[2], (ulong)a[2]);
1242: avma = av; return shifti(r,v);
1243: }
1244: }
1245:
1246: GEN
1247: mpppcm(GEN x, GEN y)
1248: {
1249: ulong av;
1250: GEN p1,p2;
1251: if (typ(x) != t_INT || typ(y) != t_INT) err(arither1);
1252: if (!signe(x)) return gzero;
1253: av = avma;
1254: p1 = mppgcd(x,y); if (!is_pm1(p1)) y = divii(y,p1);
1255: p2 = mulii(x,y);
1256: if (signe(p2) < 0) setsigne(p2,1);
1257: return gerepileuptoint(av, p2);
1258: }
1259:
1260: /*********************************************************************/
1261: /** **/
1262: /** CHINESE REMAINDERS **/
1263: /** **/
1264: /*********************************************************************/
1265:
1266: /* P.M. & M.H.
1267: *
1268: * Chinese Remainder Theorem. x and y must have the same type (integermod,
1269: * polymod, or polynomial/vector/matrix recursively constructed with these
1270: * as coefficients). Creates (with the same type) a z in the same residue
1271: * class as x and the same residue class as y, if it is possible.
1272: *
1273: * We also allow (during recursion) two identical objects even if they are
1274: * not integermod or polymod. For example, if
1275: *
1276: * x = [1. mod(5, 11), mod(X + mod(2, 7), X^2 + 1)]
1277: * y = [1, mod(7, 17), mod(X + mod(0, 3), X^2 + 1)],
1278: *
1279: * then chinois(x, y) returns
1280: *
1281: * [1, mod(16, 187), mod(X + mod(9, 21), X^2 + 1)]
1282: *
1283: * Someone else may want to allow power series, complex numbers, and
1284: * quadratic numbers.
1285: */
1286:
1287: GEN chinese(GEN x, GEN y)
1288: {
1289: return gassoc_proto(chinois,x,y);
1290: }
1291:
1292: GEN
1293: chinois(GEN x, GEN y)
1294: {
1295: ulong av,tetpil, tx = typ(x);
1296: long i,lx,vx;
1297: GEN z,p1,p2,d,u,v;
1298:
1299: if (gegal(x,y)) return gcopy(x);
1300: if (tx == typ(y)) switch(tx)
1301: {
1302: case t_POLMOD:
1303: if (gegal((GEN)x[1],(GEN)y[1])) /* same modulus */
1304: {
1305: z=cgetg(3,tx);
1306: z[1]=lcopy((GEN)x[1]);
1307: z[2]=(long)chinois((GEN)x[2],(GEN)y[2]);
1308: return z;
1309: } /* fall through */
1310: case t_INTMOD:
1311: z=cgetg(3,tx); av=avma;
1312: d=gbezout((GEN)x[1],(GEN)y[1],&u,&v);
1313: if (!gegal(gmod((GEN)x[2],d), gmod((GEN)y[2],d))) break;
1314: p1 = gdiv((GEN)x[1],d);
1315: p2 = gadd((GEN)x[2], gmul(gmul(u,p1), gadd((GEN)y[2],gneg((GEN)x[2]))));
1316:
1317: tetpil=avma; z[1]=lmul(p1,(GEN)y[1]); z[2]=lmod(p2,(GEN)z[1]);
1318: gerepilemanyvec(av,tetpil,z+1,2); return z;
1319:
1320: case t_POL:
1321: lx=lgef(x); vx=varn(x); z=cgetg(lx,tx);
1322: if (lx!=lgef(y) || vx!=varn(y)) break;
1323: z[1]=evalsigne(1)|evallgef(lx)|evalvarn(vx);
1324: for (i=2; i<lx; i++) z[i]=(long)chinois((GEN)x[i],(GEN)y[i]);
1325: return z;
1326:
1327: case t_VEC: case t_COL: case t_MAT:
1328: lx=lg(x); z=cgetg(lx,tx); if (lx!=lg(y)) break;
1329: for (i=1; i<lx; i++) z[i]=(long)chinois((GEN)x[i],(GEN)y[i]);
1330: return z;
1331: }
1332: err(talker,"incompatible arguments in chinois");
1333: return NULL; /* not reached */
1334: }
1335:
1336: /* return lift(chinois(x2 mod x1, y2 mod y1))
1337: * assume(x1,y1)=1, xi,yi integers and z1 = x1*y1
1338: */
1339: GEN
1340: chinois_int_coprime(GEN x2, GEN y2, GEN x1, GEN y1, GEN z1)
1341: {
1342: ulong av = avma;
1343: GEN ax,p1;
1344: (void)new_chunk((lgefint(z1)<<1)+lgefint(x1)+lgefint(y1)); /* HACK */
1345: ax = mulii(mpinvmod(x1,y1), x1);
1346: p1 = addii(x2, mulii(ax, subii(y2,x2)));
1347: avma = av; return modii(p1,z1);
1348: }
1349:
1350: /*********************************************************************/
1351: /** **/
1352: /** INVERSE MODULO b **/
1353: /** **/
1354: /*********************************************************************/
1355:
1356: GEN
1357: mpinvmod(GEN a, GEN m)
1358: {
1359: GEN res;
1360: if (! invmod(a,m,&res))
1361: err(invmoder,"%Z",gmodulcp(res,m));
1362: return res;
1363: }
1364:
1365: /*********************************************************************/
1366: /** **/
1367: /** POWERING MODULO (A^N mod M) **/
1368: /** **/
1369: /*********************************************************************/
1370: GEN resiimul(GEN x, GEN y);
1371: GEN resmod2n(GEN x, long y);
1372: GEN _resii(GEN x, GEN y) { return resii(x,y); }
1373:
1374: GEN
1375: init_remainder(GEN M)
1376: {
1377: GEN sM = cgetg(3, t_VEC);
1378: GEN Mr = cgetr(lgefint(M) + 1);
1379: affir(M, Mr);
1380: sM[1] = (long)M;
1381: sM[2] = (long)linv(Mr);
1382: return sM;
1383: }
1384:
1385: /* optimal on UltraSPARC */
1386: static long RESIIMUL_LIMIT = 150;
1387:
1388: void
1389: setresiilimit(long n) { RESIIMUL_LIMIT = n; }
1390:
1391: GEN
1392: powmodulo(GEN A, GEN N, GEN M)
1393: {
1394: GEN y;
1395: long av = avma,av1,lim,man,k,nb,s, *p;
1396: GEN (*mul)(GEN,GEN) = mulii;
1397: GEN (*res)(GEN,GEN) = _resii;
1398:
1399: if (typ(A) != t_INT || typ(N) != t_INT || typ(M) != t_INT) err(arither1);
1400: s = signe(N);
1401: if (!s)
1402: {
1403: k = signe(resii(A,M)); avma=av;
1404: return k? gun: gzero;
1405: }
1406: if (s < 0) { A = mpinvmod(A,M); N = absi(N); }
1407: else
1408: {
1409: A = modii(A,M);
1410: if (!signe(A)) { avma=av; return gzero; }
1411: }
1412: y = A;
1413: if (lgefint(y)==3) switch(y[2])
1414: {
1415: case 1: /* y = 1 */
1416: avma=av; return gun;
1417: case 2: /* y = 2, use shifti not mulii */
1418: mul = (GEN (*)(GEN,GEN))shifti; A = (GEN)1L;
1419: }
1420:
1421: /* TODO: Move this out of here and use for general modular computations */
1422: if ((k = vali(M)) && k == expi(M))
1423: { /* M is a power of 2 */
1424: M = (GEN)k;
1425: res = (GEN(*)(GEN,GEN))resmod2n;
1426: }
1427: else if (lgefint(M) > RESIIMUL_LIMIT && (lgefint(N) > 3 || N[2] > 4))
1428: { /* compute x % M using multiplication by 1./M */
1429: M = init_remainder(M);
1430: res = (GEN(*)(GEN,GEN))resiimul;
1431: }
1432:
1433: p = N+2; man = *p; /* see puissii */
1434: k = 1+bfffo(man); man<<=k; k = BITS_IN_LONG-k;
1435: av1=avma; lim=stack_lim(av1,1);
1436: for (nb=lgefint(N)-2;;)
1437: {
1438: for (; k; man<<=1,k--)
1439: {
1440: y = res(sqri(y), M);
1441: if (man < 0) y = res(mul(y,A), M);
1442: if (low_stack(lim, stack_lim(av1,1)))
1443: {
1444: if (DEBUGMEM>1) err(warnmem,"powmodulo");
1445: y = gerepileuptoint(av1,y);
1446: }
1447: }
1448: if (--nb == 0) break;
1449: man = *++p, k = BITS_IN_LONG;
1450: }
1451: return gerepileupto(av,y);
1452: }
1453:
1454: /*********************************************************************/
1455: /** **/
1456: /** NEXT / PRECEDING (PSEUDO) PRIME **/
1457: /** **/
1458: /*********************************************************************/
1459: GEN
1460: gnextprime(GEN n)
1461: {
1462: return garith_proto(nextprime,n,0);
1463: }
1464:
1465: GEN
1466: gprecprime(GEN n)
1467: {
1468: return garith_proto(precprime,n,0);
1469: }
1470:
1471: GEN
1472: gisprime(GEN x, long flag)
1473: {
1474: switch (flag)
1475: {
1476: case 0:
1477: return arith_proto(isprime,x,1);
1478: case 1:
1479: return garith_proto2gs(plisprime,x,0);
1480: case 2:
1481: return garith_proto2gs(plisprime,x,1);
1482: }
1483: err(flagerr,"gisprime");
1484: return 0;
1485: }
1486:
1487: long
1488: isprime(GEN x)
1489: {
1490: return millerrabin(x,10);
1491: }
1492:
1493: GEN
1494: gispsp(GEN x)
1495: {
1496: return arith_proto(ispsp,x,1);
1497: }
1498:
1499: long
1500: ispsp(GEN x)
1501: {
1502: return millerrabin(x,1);
1503: }
1504:
1505: GEN
1506: gmillerrabin(GEN n, long k)
1507: {
1508: return arith_proto2gs(millerrabin,n,k);
1509: }
1510:
1511: /*********************************************************************/
1512: /** **/
1513: /** FUNDAMENTAL DISCRIMINANTS **/
1514: /** **/
1515: /*********************************************************************/
1516:
1517: /* gissquarefree, issquarefree moved into arith2.c with the other
1518: basic arithmetic functions (issquarefree is the square of the
1519: Moebius mu function, after all...) --GN */
1520:
1521: GEN
1522: gisfundamental(GEN x)
1523: {
1524: return arith_proto(isfundamental,x,1);
1525: }
1526:
1527: long
1528: isfundamental(GEN x)
1529: {
1530: long av,r;
1531: GEN p1;
1532:
1533: if (gcmp0(x)) return 0;
1534: r=mod4(x);
1535: if (!r)
1536: {
1537: av=avma; p1=shifti(x,-2);
1538: r=mod4(p1); if (!r) return 0;
1539: if (signe(x)<0) r=4-r;
1540: r = (r==1)? 0: issquarefree(p1);
1541: avma=av; return r;
1542: }
1543: if (signe(x)<0) r=4-r;
1544: return (r==1) ? issquarefree(x) : 0;
1545: }
1546:
1547: GEN
1548: quaddisc(GEN x)
1549: {
1550: long av=avma,tetpil=avma,i,r,tx=typ(x);
1551: GEN p1,p2,f,s;
1552:
1553: if (tx != t_INT && ! is_frac_t(tx)) err(arither1);
1554: f=factor(x); p1=(GEN)f[1]; p2=(GEN)f[2];
1555: s=gun;
1556: for (i=1; i<lg(p1); i++)
1557: if ( odd(mael(p2,i,2)) ) { tetpil=avma; s=gmul(s,(GEN)p1[i]); }
1558: r=mod4(s); if (gsigne(x)<0) r=4-r;
1559: if (r>1) { tetpil=avma; s=shifti(s,2); }
1560: return gerepile(av,tetpil,s);
1561: }
1562:
1563: /*********************************************************************/
1564: /** **/
1565: /** FACTORIAL **/
1566: /** **/
1567: /*********************************************************************/
1568: GEN muluu(ulong x, ulong y);
1569:
1570: GEN
1571: mpfact(long n)
1572: {
1573: long lx,k,l, av = avma;
1574: GEN x;
1575:
1576: if (n<2)
1577: {
1578: if (n<0) err(facter);
1579: return gun;
1580: }
1581: if (n < 60)
1582: {
1583: x = gdeux;
1584: for (k=3; k<=n; k++) x = mulsi(k,x);
1585: return gerepileuptoint(av,x);
1586: }
1587: lx = 1; x = cgetg(1 + n/2, t_VEC);
1588: for (k=2;; k++)
1589: {
1590: l = n+2-k; if (l <= k) break;
1591: x[lx++] = (long)muluu(k,l);
1592: }
1593: if (l == k) x[lx++] = lstoi(k);
1594: setlg(x, lx);
1595: return gerepileuptoint(av, divide_conquer_prod(x, mulii));
1596: }
1597:
1598: GEN
1599: mpfactr(long n, long prec)
1600: {
1601: long av = avma, lim,k;
1602: GEN f = cgetr(prec);
1603:
1604: affsr(1,f);
1605: if (n<2)
1606: {
1607: if (n<0) err(facter);
1608: return f;
1609: }
1610: lim = stack_lim(av,1);
1611: for (k=2; k<=n; k++)
1612: {
1613: f = mulsr(k,f);
1614: if (low_stack(lim, stack_lim(av,1)))
1615: {
1616: if(DEBUGMEM>1) err(warnmem,"mpfactr k=%ld",k);
1617: f = gerepileuptoleaf(av,f);
1618: }
1619: }
1620: return gerepileuptoleaf(av,f);
1621: }
1622:
1623: /*******************************************************************/
1624: /** **/
1625: /** LUCAS & FIBONACCI **/
1626: /** **/
1627: /*******************************************************************/
1628:
1629: void
1630: lucas(long n, GEN *ln, GEN *ln1)
1631: {
1632: long taille,av;
1633: GEN z,t;
1634:
1635: if (!n) { *ln=stoi(2); *ln1=stoi(1); return; }
1636:
1637: taille=(long)(pariC3*(1+labs(n))+3);
1638: *ln=cgeti(taille); *ln1=cgeti(taille);
1639: av=avma; lucas(n/2,&z,&t);
1640: switch(n % 4)
1641: {
1642: case -3:
1643: addsiz(2,sqri(z),*ln1);
1644: subiiz(addsi(1,mulii(z,t)),*ln1,*ln); break;
1645: case -2:
1646: addsiz(2,sqri(z),*ln); addsiz(1,mulii(z,t),*ln1); break;
1647: case -1:
1648: subisz(sqri(z),2,*ln1);
1649: subiiz(subis(mulii(z,t),1),*ln1,*ln); break;
1650: case 0: subisz(sqri(z),2,*ln); subisz(mulii(z,t),1,*ln1); break;
1651: case 1: subisz(mulii(z,t),1,*ln); addsiz(2,sqri(t),*ln1); break;
1652: case 2: addsiz(2,sqri(z),*ln); addsiz(1,mulii(z,t),*ln1); break;
1653: case 3: addsiz(1,mulii(z,t),*ln); subisz(sqri(t),2,*ln1);
1654: }
1655: avma=av;
1656: }
1657:
1658: GEN
1659: fibo(long n)
1660: {
1661: long av = avma;
1662: GEN ln,ln1;
1663:
1664: lucas(n-1,&ln,&ln1);
1665: return gerepileupto(av, divis(addii(shifti(ln,1),ln1), 5));
1666: }
1667:
1668: /*******************************************************************/
1669: /* */
1670: /* CONTINUED FRACTIONS */
1671: /* */
1672: /*******************************************************************/
1673: static GEN sfcont2(GEN b, GEN x, long k);
1674:
1675: GEN
1676: gcf(GEN x)
1677: {
1678: return sfcont(x,x,0);
1679: }
1680:
1681: GEN
1682: gcf2(GEN b, GEN x)
1683: {
1684: return contfrac0(x,b,0);
1685: }
1686:
1687: GEN
1688: gboundcf(GEN x, long k)
1689: {
1690: return sfcont(x,x,k);
1691: }
1692:
1693: GEN
1694: contfrac0(GEN x, GEN b, long flag)
1695: {
1696: long lb,tb,i;
1697: GEN y;
1698:
1699: if (!b || gcmp0(b)) return sfcont(x,x,flag);
1700: tb = typ(b);
1701: if (tb == t_INT) return sfcont(x,x,itos(b));
1702: if (! is_matvec_t(tb)) err(typeer,"contfrac0");
1703:
1704: lb=lg(b); if (lb==1) return cgetg(1,t_VEC);
1705: if (tb != t_MAT) return sfcont2(b,x,flag);
1706: if (lg(b[1])==1) return sfcont(x,x,flag);
1707: y = (GEN) gpmalloc(lb*sizeof(long));
1708: for (i=1; i<lb; i++) y[i]=coeff(b,1,i);
1709: x = sfcont2(y,x,flag); free(y); return x;
1710: }
1711:
1712: GEN
1713: sfcont(GEN x, GEN x1, long k)
1714: {
1715: long av,lx=lg(x),tx=typ(x),e,i,j,l,l1,lx1,tetpil,f;
1716: GEN y,p1,p2,p3,p4,yp;
1717:
1718: if (is_scalar_t(tx))
1719: {
1720: if (gcmp0(x)) { y=cgetg(2,t_VEC); y[1]=zero; return y; }
1721: switch(tx)
1722: {
1723: case t_INT:
1724: y=cgetg(2,t_VEC); y[1]=lcopy(x); return y;
1725: case t_REAL:
1726: l=avma; p1=cgetg(3,t_FRACN);
1727: p2 = rcopy(x); settyp(p2,t_INT); setlgefint(p2,lx);
1728: p1[1] = (long) p2; e = bit_accuracy(lx)-1-expo(x);
1729: if (e<0) err(talker,"integral part not significant in scfont");
1730:
1731: l1 = (e>>TWOPOTBITS_IN_LONG)+3; p2=cgeti(l1);
1732: p2[1]= evalsigne(1)|evallgefint(l1);
1733: p2[2]= (1L<<(e&(BITS_IN_LONG-1)));
1734: for (i=3; i<l1; i++) p2[i]=0;
1735: p1[2]=(long) p2; p3=cgetg(3,t_FRACN);
1736: p3[2]=lcopy(p2);
1737: p3[1]=laddsi(signe(x),(GEN)p1[1]);
1738: p1=sfcont(p1,p1,k); tetpil=avma;
1739: return gerepile(l,tetpil,sfcont(p3,p1,k));
1740:
1741: case t_FRAC: case t_FRACN:
1742: l=avma; lx1=lgefint(x[2]);
1743: l1 = (long) ((double) BYTES_IN_LONG/4.0*46.093443*(lx1-2)+3);
1744: if (k>0 && ++k > 0 && l1 > k) l1 = k; /* beware overflow */
1745: if ((ulong)l1>LGBITS) l1=LGBITS;
1746: if (lgefint(x[1]) >= lx1) p1=icopy((GEN)x[1]);
1747: else affii((GEN)x[1],p1=cgeti(lx1));
1748: p2=icopy((GEN)x[2]); lx1=lg(x1);
1749: y=cgetg(l1,t_VEC); f=(x!=x1); i=0;
1750: while (!gcmp0(p2) && i<=l1-2)
1751: {
1752: i++; y[i]=ldvmdii(p1,p2,&p3);
1753: if (signe(p3)>=0) affii(p3,p1);
1754: else
1755: {
1756: p4=addii(p3,p2); affii(p4,p1); cgiv(p4);
1757: y[1]=laddsi(-1,(GEN)y[1]);
1758: }
1759: cgiv(p3); p4=p1; p1=p2; p2=p4;
1760: if (f)
1761: {
1762: if (i>=lx1) { i--; break; }
1763: if (!egalii((GEN)y[i],(GEN)x1[i]))
1764: {
1765: av=avma;
1766: if (gcmp1(absi(subii((GEN)x1[i],(GEN)y[i]))))
1767: {
1768: if (i>=lx1-1 || !gcmp1((GEN)x1[i+1]))
1769: affii((GEN)x1[i],(GEN)y[i]);
1770: }
1771: else i--;
1772: avma=av; break;
1773: }
1774: }
1775: }
1776: if (i>=2 && gcmp1((GEN)y[i]))
1777: {
1778: cgiv((GEN)y[i]); --i; cgiv((GEN)y[i]);
1779: y[i]=laddsi(1,(GEN)y[i]);
1780: }
1781: setlg(y,i+1); return gerepileupto(l, y);
1782: }
1783: err(typeer,"sfcont");
1784: }
1785:
1786: switch(tx)
1787: {
1788: case t_POL:
1789: y=cgetg(2,t_VEC); y[1]=lcopy(x); break;
1790: case t_SER:
1791: av=avma; p1=gtrunc(x); tetpil=avma;
1792: y=gerepile(av,tetpil,sfcont(p1,p1,k)); break;
1793: case t_RFRAC:
1794: case t_RFRACN:
1795: av=avma; l1=lgef(x[1]); if (lgef(x[2])>l1) l1=lgef(x[2]);
1796: if (k>0 && l1>k+1) l1=k+1;
1797: yp=cgetg(l1,t_VEC); p1=(GEN)x[1]; i=0; p2=(GEN)x[2];
1798: while (!gcmp0(p2) && i<=l1-2)
1799: {
1800: i++; yp[i]=ldivres(p1,p2,&p3);
1801: p1=p2; p2=p3;
1802: }
1803: tetpil=avma; y=cgetg(i+1,t_VEC);
1804: for (j=1; j<=i; j++) y[j]=lcopy((GEN)yp[j]);
1805: y=gerepile(av,tetpil,y); break;
1806: default: err(typeer,"sfcont");
1807: return NULL; /* not reached */
1808: }
1809: return y;
1810: }
1811:
1812: static GEN
1813: sfcont2(GEN b, GEN x, long k)
1814: {
1815: ulong av = avma;
1816: long l1 = lg(b), tx = typ(x), i;
1817: GEN y,p1;
1818:
1819: if (k)
1820: {
1821: if (k>=l1) err(typeer,"sfcont2");
1822: l1 = k+1;
1823: }
1824: y=cgetg(l1,t_VEC);
1825: if (l1==1) return y;
1826: if (is_scalar_t(tx))
1827: {
1828: if (!is_intreal_t(tx) && !is_frac_t(tx)) err(typeer,"sfcont2");
1829: }
1830: else if (tx == t_SER) x = gtrunc(x);
1831:
1832: if (!gcmp1((GEN)b[1])) x = gmul((GEN)b[1],x);
1833: i = 2; y[1] = lfloor(x); p1 = gsub(x,(GEN)y[1]);
1834: for ( ; i<l1 && !gcmp0(p1); i++)
1835: {
1836: x = gdiv((GEN)b[i],p1);
1837: if (tx == t_REAL)
1838: {
1839: long e = expo(x);
1840: if (e>0 && (e>>TWOPOTBITS_IN_LONG)+3 > lg(x)) break;
1841: }
1842: y[i] = lfloor(x);
1843: p1 = gsub(x,(GEN)y[i]);
1844: }
1845: setlg(y,i);
1846: return gerepilecopy(av,y);
1847: }
1848:
1849: GEN
1850: pnqn(GEN x)
1851: {
1852: long av=avma,tetpil,lx,ly,tx=typ(x),i;
1853: GEN y,p0,p1,q0,q1,a,b,p2,q2;
1854:
1855: if (! is_matvec_t(tx)) err(typeer,"pnqn");
1856: lx=lg(x); if (lx==1) return idmat(2);
1857: p0=gun; q0=gzero;
1858: if (tx != t_MAT)
1859: {
1860: p1=(GEN)x[1]; q1=gun;
1861: for (i=2; i<lx; i++)
1862: {
1863: a=(GEN)x[i];
1864: p2=gadd(gmul(a,p1),p0); p0=p1; p1=p2;
1865: q2=gadd(gmul(a,q1),q0); q0=q1; q1=q2;
1866: }
1867: }
1868: else
1869: {
1870: ly=lg(x[1]);
1871: if (ly==2)
1872: {
1873: p1=cgetg(lx,t_VEC); for (i=1; i<lx; i++) p1[i]=mael(x,i,1);
1874: tetpil=avma; return gerepile(av,tetpil,pnqn(p1));
1875: }
1876: if (ly!=3) err(talker,"incorrect size in pnqn");
1877: p1=gcoeff(x,2,1); q1=gcoeff(x,1,1);
1878: for (i=2; i<lx; i++)
1879: {
1880: a=gcoeff(x,2,i); b=gcoeff(x,1,i);
1881: p2=gadd(gmul(a,p1),gmul(b,p0)); p0=p1; p1=p2;
1882: q2=gadd(gmul(a,q1),gmul(b,q0)); q0=q1; q1=q2;
1883: }
1884: }
1885: tetpil=avma; y=cgetg(3,t_MAT);
1886: p2=cgetg(3,t_COL); y[1]=(long)p2; p2[1]=lcopy(p1); p2[2]=lcopy(q1);
1887: p2=cgetg(3,t_COL); y[2]=(long)p2; p2[1]=lcopy(p0); p2[2]=lcopy(q0);
1888: return gerepile(av,tetpil,y);
1889: }
1890:
1891: /* x t_INTMOD. Look for coprime integers a<=A and b<=B, such that a/b = x */
1892: GEN
1893: bestappr_mod(GEN x, GEN A, GEN B)
1894: {
1895: long i,lx,tx;
1896: GEN y;
1897: tx = typ(x);
1898: switch(tx)
1899: {
1900: case t_INTMOD:
1901: {
1902: ulong av = avma;
1903: GEN a,b,d, t = cgetg(3, t_FRAC);
1904: if (! ratlift((GEN)x[2], (GEN)x[1], &a,&b,A,B)) return NULL;
1905: if (is_pm1(b)) return icopy_av(a, (GEN)av);
1906: d = mppgcd(a,b);
1907: if (!is_pm1(d)) { avma = av; return NULL; }
1908: cgiv(d);
1909: t[1] = (long)a;
1910: t[2] = (long)b; return t;
1911: }
1912: case t_COMPLEX: case t_POL: case t_SER: case t_RFRAC:
1913: case t_RFRACN: case t_VEC: case t_COL: case t_MAT:
1914: lx = (tx==t_POL)? lgef(x): lg(x); y=cgetg(lx,tx);
1915: for (i=1; i<lontyp[tx]; i++) y[i]=x[i];
1916: for ( ; i<lx; i++)
1917: {
1918: GEN t = bestappr_mod((GEN)x[i],A,B);
1919: if (!t) return NULL;
1920: y[i] = (long)t;
1921: }
1922: return y;
1923: }
1924: err(typeer,"bestappr0");
1925: return NULL; /* not reached */
1926: }
1927:
1928: GEN
1929: bestappr(GEN x, GEN k)
1930: {
1931: ulong av = avma, tetpil;
1932: long tx,tk=typ(k),lx,i,e;
1933: GEN p0,p1,p,q0,q1,q,a,y;
1934:
1935: if (tk != t_INT)
1936: {
1937: if (tk != t_REAL && !is_frac_t(tk))
1938: err(talker,"incorrect bound type in bestappr");
1939: k = gcvtoi(k,&e);
1940: }
1941: if (signe(k) <= 0) k=gun;
1942: tx=typ(x); if (tx == t_FRACN) x = gred(x);
1943: switch(tx)
1944: {
1945: case t_INT:
1946: if (av==avma) return icopy(x);
1947: tetpil=avma; return gerepile(av,tetpil,icopy(x));
1948:
1949: case t_FRAC:
1950: if (cmpii((GEN)x[2],k) <= 0)
1951: {
1952: if (av==avma) return gcopy(x);
1953: return gerepilecopy(av,x);
1954: }
1955:
1956: case t_REAL:
1957: p1=gun; p0=gfloor(x); q1=gzero; q0=gun; a=p0;
1958: while (gcmp(q0,k)<=0)
1959: {
1960: x = gsub(x,a);
1961: if (gcmp0(x)) { p1=p0; q1=q0; break; }
1962:
1963: x = ginv(x); a = (gcmp(x,k)<0)? gfloor(x): k;
1964: p = addii(mulii(a,p0),p1); p1=p0; p0=p;
1965: q = addii(mulii(a,q0),q1); q1=q0; q0=q;
1966: }
1967: tetpil=avma; return gerepile(av,tetpil,gdiv(p1,q1));
1968:
1969: case t_COMPLEX: case t_POL: case t_SER: case t_RFRAC:
1970: case t_RFRACN: case t_VEC: case t_COL: case t_MAT:
1971: lx = (tx==t_POL)? lgef(x): lg(x); y=cgetg(lx,tx);
1972: for (i=1; i<lontyp[tx]; i++) y[i]=x[i];
1973: for ( ; i<lx; i++) y[i]=(long)bestappr((GEN)x[i],k);
1974: return y;
1975: }
1976: err(typeer,"bestappr");
1977: return NULL; /* not reached */
1978: }
1979:
1980: GEN
1981: bestappr0(GEN x, GEN a, GEN b)
1982: {
1983: ulong av;
1984: GEN t;
1985: if (!b) return bestappr(x,a);
1986: av = avma;
1987: t = bestappr_mod(x,a,b);
1988: if (!t) { avma = av; return stoi(-1); }
1989: return t;
1990: }
1991:
1992: /***********************************************************************/
1993: /** **/
1994: /** FUNDAMENTAL UNIT AND REGULATOR (QUADRATIC FIELDS) **/
1995: /** **/
1996: /***********************************************************************/
1997:
1998: GEN
1999: gfundunit(GEN x)
2000: {
2001: return garith_proto(fundunit,x,1);
2002: }
2003:
2004: static GEN
2005: get_quad(GEN f, GEN pol, long r)
2006: {
2007: GEN y, c=(GEN)f[2], p1=(GEN)c[1], q1=(GEN)c[2];
2008: y=cgetg(4,t_QUAD); y[1]=(long)pol;
2009: y[2]=r? lsubii(p1,q1): (long)p1;
2010: y[3]=(long)q1; return y;
2011: }
2012:
2013: /* replace f by f * [a,1; 1,0] */
2014: static void
2015: update_f(GEN f, GEN a)
2016: {
2017: GEN p1;
2018: p1=gcoeff(f,1,1);
2019: coeff(f,1,1)=laddii(mulii(a,p1), gcoeff(f,1,2));
2020: coeff(f,1,2)=(long)p1;
2021:
2022: p1=gcoeff(f,2,1);
2023: coeff(f,2,1)=laddii(mulii(a,p1), gcoeff(f,2,2));
2024: coeff(f,2,2)=(long)p1;
2025: }
2026:
2027: GEN
2028: fundunit(GEN x)
2029: {
2030: long av = avma, av2,tetpil,lim,r,flp,flq;
2031: GEN pol,y,a,u,v,u1,v1,sqd,f;
2032:
2033: if (typ(x) != t_INT) err(arither1);
2034: if (signe(x)<=0) err(arither2);
2035: r=mod4(x); if (r==2 || r==3) err(funder2,"fundunit");
2036:
2037: sqd=racine(x); av2=avma; lim=stack_lim(av2,2);
2038: a = shifti(addsi(r,sqd),-1);
2039: f=cgetg(3,t_MAT); f[1]=lgetg(3,t_COL); f[2]=lgetg(3,t_COL);
2040: coeff(f,1,1)=(long)a; coeff(f,1,2)=un;
2041: coeff(f,2,1)=un; coeff(f,2,2)=zero;
2042: v = gdeux; u = stoi(r);
2043: for(;;)
2044: {
2045: u1=subii(mulii(a,v),u); flp=egalii(u,u1); u=u1;
2046: v1=divii(subii(x,sqri(u)),v); flq=egalii(v,v1); v=v1;
2047: if (flq) break; a = divii(addii(sqd,u),v);
2048: if (flp) break; update_f(f,a);
2049: if (low_stack(lim, stack_lim(av2,2)))
2050: {
2051: GEN *gptr[4]; gptr[0]=&a; gptr[1]=&f; gptr[2]=&u; gptr[3]=&v;
2052: if(DEBUGMEM>1) err(warnmem,"fundunit");
2053: gerepilemany(av2,gptr,4);
2054: }
2055: }
2056: pol = quadpoly(x); y = get_quad(f,pol,r);
2057: if (!flq) v1 = y; else { update_f(f,a); v1 = get_quad(f,pol,r); }
2058:
2059: u1=gconj(y); tetpil=avma; y=gdiv(v1,u1);
2060: if (signe(y[3]) < 0) { tetpil=avma; y=gneg(y); }
2061: return gerepile(av,tetpil,y);
2062: }
2063:
2064: GEN
2065: gregula(GEN x, long prec)
2066: {
2067: return garith_proto2gs(regula,x,prec);
2068: }
2069:
2070: GEN
2071: regula(GEN x, long prec)
2072: {
2073: long av = avma,av2,lim,r,fl,rexp;
2074: GEN reg,rsqd,y,u,v,u1,v1, sqd = racine(x);
2075:
2076: if (signe(x)<=0) err(arither2);
2077: r=mod4(x); if (r==2 || r==3) err(funder2,"regula");
2078:
2079: rsqd = gsqrt(x,prec);
2080: if (egalii(sqri(sqd),x)) err(talker,"square argument in regula");
2081:
2082: rexp=0; reg=cgetr(prec); affsr(2,reg);
2083: av2=avma; lim = stack_lim(av2,2);
2084: u = stoi(r); v = gdeux;
2085: for(;;)
2086: {
2087: u1 = subii(mulii(divii(addii(u,sqd),v), v), u);
2088: v1 = divii(subii(x,sqri(u1)),v); fl = egalii(v,v1);
2089: if (fl || egalii(u,u1)) break;
2090: reg = mulrr(reg, divri(addir(u1,rsqd),v));
2091: rexp += expo(reg); setexpo(reg,0);
2092: u = u1; v = v1;
2093: if (rexp & ~EXPOBITS) err(muler4);
2094: if (low_stack(lim, stack_lim(av2,2)))
2095: {
2096: GEN *gptr[3]; gptr[0]=® gptr[1]=&u; gptr[2]=&v;
2097: if(DEBUGMEM>1) err(warnmem,"regula");
2098: gerepilemany(av2,gptr,3);
2099: }
2100: }
2101: reg = gsqr(reg); setexpo(reg,expo(reg)-1);
2102: if (fl) reg = mulrr(reg, divri(addir(u1,rsqd),v));
2103: y = mplog(divri(reg,v));
2104: if (rexp)
2105: {
2106: u1 = mulsr(rexp, glog(gdeux, prec));
2107: setexpo(u1, expo(u1)+1);
2108: y = addrr(y,u1);
2109: }
2110: return gerepileupto(av, y);
2111: }
2112:
2113: /*************************************************************************/
2114: /** **/
2115: /** CLASS NUMBER **/
2116: /** **/
2117: /*************************************************************************/
2118:
2119: static GEN
2120: gclassno(GEN x)
2121: {
2122: return garith_proto(classno,x,1);
2123: }
2124:
2125: static GEN
2126: gclassno2(GEN x)
2127: {
2128: return garith_proto(classno2,x,1);
2129: }
2130:
2131: GEN
2132: qfbclassno0(GEN x,long flag)
2133: {
2134: switch(flag)
2135: {
2136: case 0: return gclassno(x);
2137: case 1: return gclassno2(x);
2138: default: err(flagerr,"qfbclassno");
2139: }
2140: return NULL; /* not reached */
2141: }
2142:
2143: /* f^h = 1, return order(f) */
2144: static GEN
2145: find_order(GEN f, GEN h)
2146: {
2147: GEN fh, p,e;
2148: long i,j,lim;
2149:
2150: p = decomp(h);
2151: e =(GEN)p[2];
2152: p =(GEN)p[1];
2153: for (i=1; i<lg(p); i++)
2154: {
2155: lim = itos((GEN)e[i]);
2156: for (j=1; j<=lim; j++)
2157: {
2158: GEN p1 = divii(h,(GEN)p[i]);
2159: fh = powgi(f,p1);
2160: if (!is_pm1(fh[1])) break;
2161: h = p1;
2162: }
2163: }
2164: return h;
2165: }
2166:
2167: static GEN
2168: end_classno(GEN h, GEN hin, GEN forms, long lform)
2169: {
2170: long i,com;
2171: GEN a,b,p1,q,fh,fg, f = (GEN)forms[0];
2172:
2173: h = find_order(f,h); /* H = <f> */
2174: q = ground(gdiv(hin,h)); /* approximate order of G/H */
2175: for (i=1; i < lform; i++)
2176: {
2177: long av = avma;
2178: fg = powgi((GEN)forms[i], h);
2179: fh = powgi(fg, q);
2180: a = (GEN)fh[1];
2181: if (is_pm1(a)) continue;
2182: b = (GEN)fh[2]; p1 = fg;
2183: for (com=1; ; com++)
2184: {
2185: if (egalii((GEN)p1[1], a) && absi_equal((GEN)p1[2], b)) break;
2186: p1 = gmul(p1,fg);
2187: }
2188: if (signe(p1[2]) == signe(b)) com = -com;
2189: /* f_i ^ h(q+com) = 1 */
2190: avma = av; q = addsi(com,q);
2191: }
2192: return mulii(q,h);
2193: }
2194:
2195: static GEN
2196: to_form(GEN x, long f)
2197: {
2198: return redimag(primeform(x, stoi(f), 0));
2199: }
2200:
2201: static long
2202: two_rank(GEN x)
2203: {
2204: GEN p = (GEN)decomp(absi(x))[1];
2205: long l = lg(p)-1;
2206: #if 0 /* positive disc not needed */
2207: if (signe(x) > 0)
2208: {
2209: long i;
2210: for (i=1; i<=l; i++)
2211: if (mod4((GEN)p[i]) == 3) { l--; break; }
2212: }
2213: #endif
2214: return l-1;
2215: }
2216:
2217: #define MAXFORM 11
2218: #define _low(x) (__x=(GEN)x, modBIL(__x))
2219:
2220: /* h(x) for x<0 using Baby Step/Giant Step.
2221: * Assumes G is not too far from being cyclic.
2222: *
2223: * Compute G^2 instead of G so as to kill most of the non-cyclicity */
2224: GEN
2225: classno(GEN x)
2226: {
2227: long av = avma, av2,c,lforms,k,l,i,j,j1,j2,lim,com,s, forms[MAXFORM];
2228: long r2;
2229: GEN a,b,count,index,tabla,tablb,hash,p1,p2,hin,h,f,fh,fg,ftest;
2230: GEN __x;
2231: byteptr p = diffptr;
2232:
2233: if (typ(x) != t_INT) err(arither1);
2234: s=signe(x); if (s>=0) return classno2(x);
2235:
2236: k=mod4(x); if (k==1 || k==2) err(funder2,"classno");
2237: if (gcmpgs(x,-12) >= 0) return gun;
2238:
2239: p2 = gsqrt(absi(x),DEFAULTPREC);
2240: p1 = divrr(p2,mppi(DEFAULTPREC));
2241: p2 = gtrunc(shiftr(gsqrt(p2,DEFAULTPREC),1));
2242: s = 1000;
2243: if (signe(p2))
2244: {
2245: if (is_bigint(p2)) err(talker,"discriminant too big in classno");
2246: s = itos(p2); if (s < 1000) s = 1000;
2247: }
2248: r2 = two_rank(x);
2249: c=lforms=0;
2250: while (c <= s && *p)
2251: {
2252: c += *p++; k = krogs(x,c);
2253: if (!k) continue;
2254:
2255: av2 = avma;
2256: if (k > 0)
2257: {
2258: divrsz(mulsr(c,p1),c-1, p1);
2259: if (lforms < MAXFORM && c > 2) forms[lforms++] = c;
2260: }
2261: else divrsz(mulsr(c,p1),c+1, p1);
2262: avma = av2;
2263: }
2264: h = hin = shifti(ground(p1), -r2);
2265: s = 2*itos(gceil(gpui(p1,ginv(stoi(4)),DEFAULTPREC)));
2266: if (s>=10000) s=10000;
2267:
2268: count = new_chunk(256); for (i=0; i<=255; i++) count[i]=0;
2269: index = new_chunk(257);
2270: tabla = new_chunk(10000);
2271: tablb = new_chunk(10000);
2272: hash = new_chunk(10000);
2273: f = gsqr(to_form(x, forms[0]));
2274: p1 = fh = powgi(f, h);
2275: for (i=0; i<s; i++)
2276: {
2277: tabla[i] = _low(p1[1]);
2278: tablb[i] = _low(p1[2]);
2279: count[tabla[i]&255]++; p1=compimag(p1,f);
2280: }
2281: index[0]=0; for (i=0; i<=254; i++) index[i+1] = index[i]+count[i];
2282: for (i=0; i<s; i++) hash[index[tabla[i]&255]++]=i;
2283: index[0]=0; for (i=0; i<=255; i++) index[i+1] = index[i]+count[i];
2284:
2285: fg=gpuigs(f,s); av2=avma; lim=stack_lim(av2,2);
2286: ftest=gpuigs(p1,0);
2287: for (com=0; ; com++)
2288: {
2289: a = (GEN)ftest[1]; k = _low(a);
2290: b = (GEN)ftest[2]; l = _low(b); j = k&255;
2291: for (j1=index[j]; j1 < index[j+1]; j1++)
2292: {
2293: j2 = hash[j1];
2294: if (tabla[j2] == k && tablb[j2] == l)
2295: {
2296: p1 = gmul(gpuigs(f,j2),fh);
2297: if (egalii((GEN)p1[1], a) && absi_equal((GEN)p1[2], b))
2298: { /* p1 = ftest or ftest^(-1), we are done */
2299: if (signe(p1[2]) == signe(b)) com = -com;
2300: h = addii(addis(h,j2), mulss(s,com));
2301: forms[0] = (long)f;
2302: for (i=1; i<lforms; i++)
2303: forms[i] = (long)gsqr(to_form(x, forms[i]));
2304: h = end_classno(h,hin,forms,lforms);
2305: return gerepileuptoint(av, shifti(h, r2));
2306: }
2307: }
2308: }
2309: ftest = gmul(ftest,fg);
2310: if (is_pm1(ftest[1])) err(impl,"classno with too small order");
2311: if (low_stack(lim, stack_lim(av2,2))) ftest = gerepileupto(av2,ftest);
2312: }
2313: }
2314:
2315: /* use Euler products */
2316: GEN
2317: classno2(GEN x)
2318: {
2319: long av=avma,tetpil,n,i,k,s=signe(x),fl2;
2320: GEN p1,p2,p3,p4,p5,p7,p8,pi4,reg,logd,d,fd;
2321:
2322: if (typ(x) != t_INT) err(arither1);
2323: if (!s) err(arither2);
2324: if (s<0 && gcmpgs(x,-12) >= 0) return gun;
2325:
2326: p1=auxdecomp(absi(x),1); p2=(GEN)p1[2]; p1=(GEN)p1[1];
2327: n=lg(p1); d=gun;
2328: for (i=1; i<n; i++)
2329: if (mod2((GEN)p2[i])) d=mulii(d,(GEN)p1[i]);
2330: fl2 = (mod4(x)==0); /* 4 | x */
2331: if (mod4(d) != 2-s)
2332: {
2333: if (!fl2) err(funder2,"classno2");
2334: d = shifti(d,2);
2335: }
2336: else fl2=0;
2337: p8=gun; fd = (s<0)? negi(d): d; /* d = abs(fd) */
2338: for (i=1; i<n; i++)
2339: {
2340: k=itos((GEN)p2[i]); p4=(GEN)p1[i];
2341: if (fl2 && i==1) k -= 2; /* p4 = 2 */
2342: if (k>=2)
2343: {
2344: p8=mulii(p8,subis(p4,kronecker(fd,p4)));
2345: if (k>=4) p8=mulii(p8,gpuigs(p4,(k>>1)-1));
2346: }
2347: }
2348: if (s<0 && lgefint(d)==3)
2349: {
2350: switch(d[2])
2351: {
2352: case 4: p8=gdivgs(p8,2); break;
2353: case 3: p8=gdivgs(p8,3); break;
2354: }
2355: if (d[2] < 12) /* |fd| < 12*/
2356: { tetpil=avma; return gerepile(av,tetpil,icopy(p8)); }
2357: }
2358:
2359: pi4 = mppi(DEFAULTPREC); logd = glog(d,DEFAULTPREC);
2360: p1 = gsqrt(gdiv(gmul(d,logd),gmul2n(pi4,1)),DEFAULTPREC);
2361: p4 = divri(pi4,d); p7 = ginv(mpsqrt(pi4));
2362: if (s > 0)
2363: {
2364: reg = regula(d,DEFAULTPREC);
2365: p2 = gsubsg(1, gmul2n(gdiv(glog(reg,DEFAULTPREC),logd),1));
2366: p3 = gsqrt(gdivsg(2,logd),DEFAULTPREC);
2367: if (gcmp(p2,p3)>=0) p1 = gmul(p2,p1);
2368: }
2369: else reg = NULL; /* for gcc -Wall */
2370: p1 = gtrunc(p1); n=p1[2];
2371: if (lgefint(p1)!=3 || n<0)
2372: err(talker,"discriminant too large in classno");
2373:
2374: p1 = gsqrt(d,DEFAULTPREC); p3 = gzero;
2375: if (s > 0)
2376: {
2377: for (i=1; i<=n; i++)
2378: {
2379: k=krogs(fd,i);
2380: if (k)
2381: {
2382: p2 = mulir(mulss(i,i),p4);
2383: p5 = subsr(1,mulrr(p7,incgam3(ghalf,p2,DEFAULTPREC)));
2384: p5 = addrr(divrs(mulrr(p1,p5),i),eint1(p2,DEFAULTPREC));
2385: p3 = (k>0)? addrr(p3,p5): subrr(p3,p5);
2386: }
2387: }
2388: p3 = shiftr(divrr(p3,reg),-1);
2389: if (!egalii(x,fd)) /* x != p3 */
2390: p8 = gdiv(p8,ground(gdiv(regula(x,DEFAULTPREC),reg)));
2391: }
2392: else
2393: {
2394: p1 = gdiv(p1,pi4);
2395: for (i=1; i<=n; i++)
2396: {
2397: k=krogs(fd,i);
2398: if (k)
2399: {
2400: p2=mulir(mulss(i,i),p4);
2401: p5=subsr(1,mulrr(p7,incgam3(ghalf,p2,DEFAULTPREC)));
2402: p5=addrr(p5,divrr(divrs(p1,i),mpexp(p2)));
2403: p3 = (k>0)? addrr(p3,p5): subrr(p3,p5);
2404: }
2405: }
2406: }
2407: p3=ground(p3); tetpil=avma;
2408: return gerepile(av,tetpil,mulii(p8,p3));
2409: }
2410:
2411: GEN
2412: hclassno(GEN x)
2413: {
2414: long d,a,b,h,b2,f;
2415:
2416: d= -itos(x); if (d>0 || (d & 3) > 1) return gzero;
2417: if (!d) return gdivgs(gun,-12);
2418: if (-d > (VERYBIGINT>>1))
2419: err(talker,"discriminant too big in hclassno. Use quadclassunit");
2420: h=0; b=d&1; b2=(b-d)>>2; f=0;
2421: if (!b)
2422: {
2423: for (a=1; a*a<b2; a++)
2424: if (b2%a == 0) h++;
2425: f = (a*a==b2); b=2; b2=(4-d)>>2;
2426: }
2427: while (b2*3+d<0)
2428: {
2429: if (b2%b == 0) h++;
2430: for (a=b+1; a*a<b2; a++)
2431: if (b2%a == 0) h+=2;
2432: if (a*a==b2) h++;
2433: b+=2; b2=(b*b-d)>>2;
2434: }
2435: if (b2*3+d==0)
2436: {
2437: GEN y = cgetg(3,t_FRAC);
2438: y[1]=lstoi(3*h+1); y[2]=lstoi(3); return y;
2439: }
2440: if (f) return gaddsg(h,ghalf);
2441: return stoi(h);
2442: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to arith1.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 / basemath