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