Annotation of OpenXM_contrib/pari/src/basemath/arith1.c, Revision 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