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