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