Annotation of OpenXM_contrib/pari/src/basemath/base3.c, Revision 1.1
1.1 ! maekawa 1: /*******************************************************************/
! 2: /* */
! 3: /* BASIC NF OPERATIONS */
! 4: /* */
! 5: /*******************************************************************/
! 6: /* $Id: base3.c,v 1.1.1.1 1999/09/16 13:47:20 karim Exp $ */
! 7: #include "pari.h"
! 8: int ker_trivial_mod_p(GEN x, GEN p);
! 9: long ideal_is_zk(GEN ideal,long N);
! 10: GEN idealaddtoone_i(GEN nf, GEN x, GEN y);
! 11:
! 12: /*******************************************************************/
! 13: /* */
! 14: /* OPERATIONS OVER NUMBER FIELD ELEMENTS. */
! 15: /* These are always represented as column vectors over the */
! 16: /* integral basis nf[7] */
! 17: /* */
! 18: /*******************************************************************/
! 19:
! 20: int
! 21: isnfscalar(GEN x)
! 22: {
! 23: long lx=lg(x),i;
! 24:
! 25: for (i=2; i<lx; i++)
! 26: if (!gcmp0((GEN)x[i])) return 0;
! 27: return 1;
! 28: }
! 29:
! 30: static GEN
! 31: checknfelt_mod(GEN nf, GEN x)
! 32: {
! 33: if (gegal((GEN)x[1],(GEN)nf[1])) return (GEN) x[2];
! 34: err(talker,"not the same polynomial in number field operation");
! 35: return NULL; /* not reached */
! 36: }
! 37:
! 38: static GEN
! 39: scal_mul(GEN nf, GEN x, GEN y, long ty)
! 40: {
! 41: long av=avma, tetpil;
! 42: GEN p1;
! 43:
! 44: if (!is_extscalar_t(ty))
! 45: {
! 46: if (ty!=t_COL) err(typeer,"nfmul");
! 47: y = gmul((GEN)nf[7],y);
! 48: }
! 49: p1 = gmul(x,y); tetpil=avma;
! 50: return gerepile(av,tetpil,algtobasis(nf,p1));
! 51: }
! 52:
! 53: /* product of x and y in nf */
! 54: GEN
! 55: element_mul(GEN nf, GEN x, GEN y)
! 56: {
! 57: long av,i,j,k,N,tx,ty;
! 58: GEN s,v,c,p1,tab;
! 59:
! 60: if (x == y) return element_sqr(nf,x);
! 61:
! 62: tx=typ(x); ty=typ(y);
! 63: nf=checknf(nf); tab = (GEN)nf[9]; N=lgef(nf[1])-3;
! 64: if (tx==t_POLMOD) x=checknfelt_mod(nf,x);
! 65: if (ty==t_POLMOD) y=checknfelt_mod(nf,y);
! 66: if (is_extscalar_t(tx)) return scal_mul(nf,x,y,ty);
! 67: if (is_extscalar_t(ty)) return scal_mul(nf,y,x,tx);
! 68: if (isnfscalar(x)) return gmul((GEN)x[1],y);
! 69: if (isnfscalar(y)) return gmul((GEN)y[1],x);
! 70:
! 71: v=cgetg(N+1,t_COL); av=avma;
! 72: for (k=1; k<=N; k++)
! 73: {
! 74: if (k == 1)
! 75: s = gmul((GEN)x[1],(GEN)y[1]);
! 76: else
! 77: s = gadd(gmul((GEN)x[1],(GEN)y[k]),
! 78: gmul((GEN)x[k],(GEN)y[1]));
! 79: for (i=2; i<=N; i++)
! 80: {
! 81: c=gcoeff(tab,k,(i-1)*N+i);
! 82: if (signe(c))
! 83: {
! 84: p1 = gmul((GEN)x[i],(GEN)y[i]);
! 85: if (!gcmp1(c)) p1 = gmul(p1,c);
! 86: s = gadd(s, p1);
! 87: }
! 88: for (j=i+1; j<=N; j++)
! 89: {
! 90: c=gcoeff(tab,k,(i-1)*N+j);
! 91: if (signe(c))
! 92: {
! 93: p1 = gadd(gmul((GEN)x[i],(GEN)y[j]),
! 94: gmul((GEN)x[j],(GEN)y[i]));
! 95: if (!gcmp1(c)) p1 = gmul(p1,c);
! 96: s = gadd(s, p1);
! 97: }
! 98: }
! 99: }
! 100: v[k]=(long)gerepileupto(av,s); av=avma;
! 101: }
! 102: return v;
! 103: }
! 104:
! 105: /* inverse of x in nf */
! 106: GEN
! 107: element_inv(GEN nf, GEN x)
! 108: {
! 109: long av=avma,tetpil,flx,i,N,tx=typ(x);
! 110: GEN p1,p;
! 111:
! 112: nf=checknf(nf); N=lgef(nf[1])-3;
! 113: if (is_extscalar_t(tx))
! 114: {
! 115: if (tx==t_POLMOD) checknfelt_mod(nf,x);
! 116: else if (tx==t_POL) x=gmodulcp(x,(GEN)nf[1]);
! 117: p1=ginv(x); tetpil=avma;
! 118: return gerepile(av,tetpil,algtobasis(nf,p1));
! 119: }
! 120: if (isnfscalar(x))
! 121: {
! 122: p1=cgetg(N+1,t_COL); p1[1]=linv((GEN)x[1]);
! 123: for (i=2; i<=N; i++) p1[i]=lcopy((GEN)x[i]);
! 124: return p1;
! 125: }
! 126: flx=1;
! 127: for (i=1; i<=N; i++)
! 128: if (typ(x[i])==t_INTMOD)
! 129: {
! 130: p=gmael(x,i,1); x=lift(x);
! 131: flx=0; break;
! 132: }
! 133: p1 = ginvmod(gmul((GEN)nf[7],x), (GEN)nf[1]);
! 134: p1 = algtobasis_intern(nf,p1);
! 135:
! 136: if (flx) return gerepileupto(av,p1);
! 137: tetpil=avma; return gerepile(av,tetpil,Fp_vec(p1, p));
! 138: }
! 139:
! 140: /* quotient of x and y in nf */
! 141: GEN
! 142: element_div(GEN nf, GEN x, GEN y)
! 143: {
! 144: long av=avma,tetpil,flx,i,N,tx=typ(x),ty=typ(y);
! 145: GEN p1,p;
! 146:
! 147: nf=checknf(nf); N=lgef(nf[1])-3;
! 148: if (tx==t_POLMOD) checknfelt_mod(nf,x);
! 149: else if (tx==t_POL) x=gmodulcp(x,(GEN)nf[1]);
! 150:
! 151: if (ty==t_POLMOD) checknfelt_mod(nf,y);
! 152: else if (ty==t_POL) y=gmodulcp(y,(GEN)nf[1]);
! 153:
! 154: if (is_extscalar_t(tx))
! 155: {
! 156: if (is_extscalar_t(ty)) p1=gdiv(x,y);
! 157: else
! 158: {
! 159: if (ty!=t_COL) err(typeer,"nfdiv");
! 160: p1=gdiv(x,gmodulcp(gmul((GEN)nf[7],y),(GEN)nf[1]));
! 161: }
! 162: tetpil=avma; return gerepile(av,tetpil,algtobasis(nf,p1));
! 163: }
! 164: if (is_extscalar_t(ty))
! 165: {
! 166: if (tx!=t_COL) err(typeer,"nfdiv");
! 167: p1=gdiv(gmodulcp(gmul((GEN)nf[7],x),(GEN)nf[1]),y);
! 168: tetpil=avma; return gerepile(av,tetpil,algtobasis(nf,p1));
! 169: }
! 170:
! 171: if (isnfscalar(y)) return gdiv(x,(GEN)y[1]);
! 172: if (isnfscalar(x))
! 173: {
! 174: p1=element_inv(nf,y); tetpil=avma;
! 175: return gerepile(av,tetpil,gmul((GEN)x[1],p1));
! 176: }
! 177:
! 178: flx=1;
! 179: for (i=1; i<=N; i++)
! 180: if (typ(x[i])==t_INTMOD)
! 181: {
! 182: p=gmael(x,i,1); x=lift(x);
! 183: flx=0; break;
! 184: }
! 185: for (i=1; i<=N; i++)
! 186: if (typ(y[i])==t_INTMOD)
! 187: {
! 188: p=gmael(y,i,1); y=lift(y);
! 189: flx=0; break;
! 190: }
! 191:
! 192: p1 = gmul(gmul((GEN)nf[7],x), ginvmod(gmul((GEN)nf[7],y), (GEN)nf[1]));
! 193: p1 = algtobasis_intern(nf, gres(p1, (GEN)nf[1]));
! 194: if (flx) return gerepileupto(av,p1);
! 195: tetpil=avma; return gerepile(av,tetpil, Fp_vec(p1,p));
! 196: }
! 197:
! 198: /* product of INTEGERS (i.e vectors with integral coeffs) x and y in nf */
! 199: GEN
! 200: element_muli(GEN nf, GEN x, GEN y)
! 201: {
! 202: long av,i,j,k,N=lgef(nf[1])-3;
! 203: GEN p1,s,v,c,tab = (GEN)nf[9];
! 204:
! 205: v=cgetg(N+1,t_COL); av=avma;
! 206: for (k=1; k<=N; k++)
! 207: {
! 208: if (k == 1)
! 209: s = mulii((GEN)x[1],(GEN)y[1]);
! 210: else
! 211: s = addii(mulii((GEN)x[1],(GEN)y[k]),
! 212: mulii((GEN)x[k],(GEN)y[1]));
! 213: for (i=2; i<=N; i++)
! 214: {
! 215: c=gcoeff(tab,k,(i-1)*N+i);
! 216: if (signe(c))
! 217: {
! 218: p1 = mulii((GEN)x[i],(GEN)y[i]);
! 219: if (!gcmp1(c)) p1 = mulii(p1,c);
! 220: s = addii(s,p1);
! 221: }
! 222: for (j=i+1; j<=N; j++)
! 223: {
! 224: c=gcoeff(tab,k,(i-1)*N+j);
! 225: if (signe(c))
! 226: {
! 227: p1 = addii(mulii((GEN)x[i],(GEN)y[j]),
! 228: mulii((GEN)x[j],(GEN)y[i]));
! 229: if (!gcmp1(c)) p1 = mulii(p1,c);
! 230: s = addii(s,p1);
! 231: }
! 232: }
! 233: }
! 234: v[k]=(long) gerepileuptoint(av,s); av=avma;
! 235: }
! 236: return v;
! 237: }
! 238:
! 239: /* product of INTEGERS (i.e vectors with integral coeffs) x and y in nf */
! 240: GEN
! 241: element_sqri(GEN nf, GEN x)
! 242: {
! 243: long av,i,j,k,N=lgef(nf[1])-3;
! 244: GEN p1,s,v,c,tab = (GEN)nf[9];
! 245:
! 246: v=cgetg(N+1,t_COL); av=avma;
! 247: for (k=1; k<=N; k++)
! 248: {
! 249: if (k == 1)
! 250: s = sqri((GEN)x[1]);
! 251: else
! 252: s = shifti(mulii((GEN)x[1],(GEN)x[k]), 1);
! 253: for (i=2; i<=N; i++)
! 254: {
! 255: c=gcoeff(tab,k,(i-1)*N+i);
! 256: if (signe(c))
! 257: {
! 258: p1 = sqri((GEN)x[i]);
! 259: if (!gcmp1(c)) p1 = mulii(p1,c);
! 260: s = addii(s,p1);
! 261: }
! 262: for (j=i+1; j<=N; j++)
! 263: {
! 264: c=gcoeff(tab,k,(i-1)*N+j);
! 265: if (signe(c))
! 266: {
! 267: p1 = shifti(mulii((GEN)x[i],(GEN)x[j]),1);
! 268: if (!gcmp1(c)) p1 = mulii(p1,c);
! 269: s = addii(s,p1);
! 270: }
! 271: }
! 272: }
! 273: v[k]=(long) gerepileuptoint(av,s); av=avma;
! 274: }
! 275: return v;
! 276: }
! 277:
! 278: /* square of x in nf */
! 279: GEN
! 280: element_sqr(GEN nf, GEN x)
! 281: {
! 282: long av,i,j,k,N=lgef(nf[1])-3;
! 283: GEN p1,s,v,c, tab = (GEN)nf[9];
! 284:
! 285: if (isnfscalar(x))
! 286: {
! 287: s=cgetg(N+1,t_COL); s[1]=lsqr((GEN)x[1]);
! 288: for (i=2; i<=N; i++) s[i]=lcopy((GEN)x[i]);
! 289: return s;
! 290: }
! 291: v=cgetg(N+1,t_COL); av=avma;
! 292: for (k=1; k<=N; k++)
! 293: {
! 294: if (k == 1)
! 295: s = gsqr((GEN)x[1]);
! 296: else
! 297: s = gmul2n(gmul((GEN)x[1],(GEN)x[k]), 1);
! 298: for (i=2; i<=N; i++)
! 299: {
! 300: c=gcoeff(tab,k,(i-1)*N+i);
! 301: if (signe(c))
! 302: {
! 303: p1 = gsqr((GEN)x[i]);
! 304: if (!gcmp1(c)) p1 = gmul(p1,c);
! 305: s = gadd(s,p1);
! 306: }
! 307: for (j=i+1; j<=N; j++)
! 308: {
! 309: c=gcoeff(tab,k,(i-1)*N+j);
! 310: if (signe(c))
! 311: {
! 312: p1 = gmul((GEN)x[i],(GEN)x[j]);
! 313: p1 = gcmp1(c)? gmul2n(p1,1): gmul(p1,shifti(c,1));
! 314: s = gadd(s,p1);
! 315: }
! 316: }
! 317: }
! 318: v[k]=(long)gerepileupto(av,s); av=avma;
! 319: }
! 320: return v;
! 321: }
! 322:
! 323: /* Compute x^n in nf, left-shift binary powering */
! 324: GEN
! 325: element_pow(GEN nf, GEN x, GEN n)
! 326: {
! 327: long s,av=avma,N,i,j,m;
! 328: GEN y,p1;
! 329:
! 330: if (typ(n)!=t_INT) err(talker,"not an integer exponent in nfpow");
! 331: nf=checknf(nf); N=lgef(nf[1])-3;
! 332: s=signe(n); if (!s) return gscalcol_i(gun,N);
! 333: if (typ(x)!=t_COL) x=algtobasis(nf,x);
! 334:
! 335: if (isnfscalar(x))
! 336: {
! 337: y = gscalcol_i(gun,N);
! 338: y[1] = (long)powgi((GEN)x[1],n); return y;
! 339: }
! 340: p1 = n+2; m = *p1;
! 341: y=x; j=1+bfffo(m); m<<=j; j = BITS_IN_LONG-j;
! 342: for (i=lgefint(n)-2;;)
! 343: {
! 344: for (; j; m<<=1,j--)
! 345: {
! 346: y = element_sqr(nf, y);
! 347: if (m<0) y=element_mul(nf, y, x);
! 348: }
! 349: if (--i == 0) break;
! 350: m = *++p1, j = BITS_IN_LONG;
! 351: }
! 352: if (s<0) y = element_inv(nf, y);
! 353: return av==avma? gcopy(y): gerepileupto(av,y);
! 354: }
! 355:
! 356: /* x in Z^n, compute lift(x^n mod p) */
! 357: GEN
! 358: element_pow_mod_p(GEN nf, GEN x, GEN n, GEN p)
! 359: {
! 360: long s,av=avma,N,i,j,m;
! 361: GEN y,p1;
! 362:
! 363: if (typ(n)!=t_INT) err(talker,"not an integer exponent in nfpow");
! 364: nf=checknf(nf); N=lgef(nf[1])-3;
! 365: s=signe(n); if (!s) return gscalcol_i(gun,N);
! 366: if (typ(x)!=t_COL) x=algtobasis(nf,x);
! 367:
! 368: if (isnfscalar(x))
! 369: {
! 370: y = gscalcol_i(gun,N);
! 371: y[1] = (long)powmodulo((GEN)x[1],n,p); return y;
! 372: }
! 373: p1 = n+2; m = *p1;
! 374: y=x; j=1+bfffo(m); m<<=j; j = BITS_IN_LONG-j;
! 375: for (i=lgefint(n)-2;;)
! 376: {
! 377: for (; j; m<<=1,j--)
! 378: {
! 379: y = element_sqri(nf, y);
! 380: if (m<0) y=element_muli(nf, y, x);
! 381: y = Fp_vec_red(y, p);
! 382: }
! 383: if (--i == 0) break;
! 384: m = *++p1, j = BITS_IN_LONG;
! 385: }
! 386: if (s<0) y = Fp_vec_red(element_inv(nf,y), p);
! 387: return av==avma? gcopy(y): gerepileupto(av,y);
! 388: }
! 389:
! 390: /* x = I-th vector of the Z-basis of Z_K, in Z^n, compute lift(x^n mod p) */
! 391: GEN
! 392: element_powid_mod_p(GEN nf, long I, GEN n, GEN p)
! 393: {
! 394: long s,av=avma,N,i,j,m;
! 395: GEN y,p1;
! 396:
! 397: if (typ(n)!=t_INT) err(talker,"not an integer exponent in nfpow");
! 398: nf=checknf(nf); N=lgef(nf[1])-3;
! 399: s=signe(n);
! 400: if (!s || I == 1) return gscalcol_i(gun,N);
! 401: p1 = n+2; m = *p1;
! 402: y = zerocol(N); y[I] = un;
! 403: j=1+bfffo(m); m<<=j; j = BITS_IN_LONG-j;
! 404: for (i=lgefint(n)-2;;)
! 405: {
! 406: for (; j; m<<=1,j--)
! 407: {
! 408: y = element_sqri(nf, y);
! 409: if (m<0) y=element_mulid(nf, y, I);
! 410: y = Fp_vec_red(y, p);
! 411: }
! 412: if (--i == 0) break;
! 413: m = *++p1, j = BITS_IN_LONG;
! 414: }
! 415: if (s<0) y = Fp_vec_red(element_inv(nf,y), p);
! 416: return av==avma? gcopy(y): gerepileupto(av,y);
! 417: }
! 418:
! 419: /* Outputs x.w_i, where w_i is the i-th elt of the integral basis */
! 420: GEN
! 421: element_mulid(GEN nf, GEN x, long i)
! 422: {
! 423: long av,j,k,N;
! 424: GEN s,v,c,p1, tab;
! 425:
! 426: if (i==1) return gcopy(x);
! 427: N = lgef(nf[1])-3;
! 428: tab = (GEN)nf[9]; tab += (i-1)*N;
! 429: v=cgetg(N+1,t_COL); av=avma;
! 430: for (k=1; k<=N; k++)
! 431: {
! 432: s = gzero;
! 433: for (j=1; j<=N; j++)
! 434: if (signe(c = gcoeff(tab,k,j)) && !gcmp0(p1 = (GEN)x[j]))
! 435: {
! 436: if (!gcmp1(c)) p1 = gmul(p1,c);
! 437: s = gadd(s,p1);
! 438: }
! 439: v[k]=lpileupto(av,s); av=avma;
! 440: }
! 441: return v;
! 442: }
! 443:
! 444: /* valuation of integer x, with resp. to prime ideal P above p.
! 445: * p.P^(-1) = b Z_K, v <= val_p(norm(x)), and N = deg(nf)
! 446: */
! 447: long
! 448: int_elt_val(GEN nf, GEN x, GEN p, GEN b, long v)
! 449: {
! 450: long i,k,w, N = lgef(nf[1])-3;
! 451: GEN r,a,y, mul = cgetg(N+1,t_MAT);
! 452:
! 453: for (i=1; i<=N; i++) mul[i] = (long)element_mulid(nf,b,i);
! 454: y = cgetg(N+1, t_COL); /* will hold the new x */
! 455: x = dummycopy(x);
! 456: for(w=0; w<=v; w++)
! 457: {
! 458: for (i=1; i<=N; i++)
! 459: { /* compute (x.b)_i */
! 460: a = mulii((GEN)x[1], gcoeff(mul,i,1));
! 461: for (k=2; k<=N; k++) a = addii(a, mulii((GEN)x[k], gcoeff(mul,i,k)));
! 462: /* is it divisible by p ? */
! 463: y[i] = ldvmdii(a,p,&r);
! 464: if (signe(r)) return w;
! 465: }
! 466: r=x; x=y; y=r;
! 467: }
! 468: return w;
! 469: }
! 470:
! 471: long
! 472: element_val(GEN nf, GEN x, GEN vp)
! 473: {
! 474: long av = avma,N,w,vcx,e;
! 475: GEN cx,p;
! 476:
! 477: if (gcmp0(x)) return VERYBIGINT;
! 478: nf=checknf(nf); N=lgef(nf[1])-3;
! 479: checkprimeid(vp); p=(GEN)vp[1]; e=itos((GEN)vp[3]);
! 480: switch(typ(x))
! 481: {
! 482: case t_INT: case t_FRAC: case t_FRACN:
! 483: return ggval(x,p)*e;
! 484: case t_POLMOD: x = (GEN)x[2]; /* fall through */
! 485: case t_POL:
! 486: x = algtobasis_intern(nf,x); break;
! 487: case t_COL:
! 488: if (lg(x)==N+1) break;
! 489: default: err(typeer,"element_val");
! 490: }
! 491: if (isnfscalar(x)) return ggval((GEN)x[1],p)*e;
! 492:
! 493: cx = content(x);
! 494: if (gcmp1(cx)) vcx=0; else { x = gdiv(x,cx); vcx = ggval(cx,p); }
! 495: w = int_elt_val(nf,x,p,(GEN)vp[5],VERYBIGINT);
! 496: avma=av; return w + vcx*e;
! 497: }
! 498:
! 499: /* d = a multiple of norm(x), x integral */
! 500: long
! 501: element_val2(GEN nf, GEN x, GEN d, GEN vp)
! 502: {
! 503: GEN p = (GEN)vp[1];
! 504: long av, v = ggval(d,p);
! 505:
! 506: if (!v) return 0;
! 507: av=avma;
! 508: v = int_elt_val(nf,x,p,(GEN)vp[5],v);
! 509: avma=av; return v;
! 510: }
! 511:
! 512: /* polegal without comparing variables */
! 513: long
! 514: polegal_spec(GEN x, GEN y)
! 515: {
! 516: long i = lgef(x);
! 517:
! 518: if (i != lgef(y)) return 0;
! 519: for (i--; i > 1; i--)
! 520: if (!gegal((GEN)x[i],(GEN)y[i])) return 0;
! 521: return 1;
! 522: }
! 523:
! 524: GEN
! 525: basistoalg(GEN nf, GEN x)
! 526: {
! 527: long tx=typ(x),lx=lg(x),i;
! 528: GEN z;
! 529:
! 530: nf=checknf(nf);
! 531: switch(tx)
! 532: {
! 533: case t_COL:
! 534: for (i=1; i<lx; i++)
! 535: {
! 536: long t = typ(x[i]);
! 537: if (is_matvec_t(t)) break;
! 538: }
! 539: if (i==lx)
! 540: {
! 541: z = cgetg(3,t_POLMOD);
! 542: z[1] = lcopy((GEN)nf[1]);
! 543: z[2] = lmul((GEN)nf[7],x); return z;
! 544: }
! 545: /* fall through */
! 546:
! 547: case t_VEC: case t_MAT: z=cgetg(lx,tx);
! 548: for (i=1; i<lx; i++) z[i]=(long)basistoalg(nf,(GEN)x[i]);
! 549: return z;
! 550:
! 551: case t_POLMOD:
! 552: if (!polegal_spec((GEN)nf[1],(GEN)x[1]))
! 553: err(talker,"not the same number field in basistoalg");
! 554: return gcopy(x);
! 555: default: z=cgetg(3,t_POLMOD);
! 556: z[1]=lcopy((GEN)nf[1]);
! 557: z[2]=lmul(x,polun[varn(nf[1])]); return z;
! 558: }
! 559: }
! 560:
! 561: /* return the (N-dimensional) vector of coeffs of p */
! 562: GEN
! 563: pol_to_vec(GEN x, long N)
! 564: {
! 565: long i,l=lgef(x)-1;
! 566: GEN z = cgetg(N+1,t_COL);
! 567: x++;
! 568: for (i=1; i<l ; i++) z[i]=x[i];
! 569: for ( ; i<=N; i++) z[i]=zero;
! 570: return z;
! 571: }
! 572:
! 573: /* valid for scalars and polynomial, degree less than N.
! 574: * No garbage collecting. No check (SEGV for vectors).
! 575: */
! 576: GEN
! 577: algtobasis_intern(GEN nf,GEN x)
! 578: {
! 579: long i,tx=typ(x),N=lgef(nf[1])-3;
! 580: GEN z;
! 581:
! 582: if (tx==t_POLMOD) { x=(GEN)x[2]; tx=typ(x); }
! 583: if (tx==t_POL)
! 584: {
! 585: if (lgef(x)-3 >= N) x=gres(x,(GEN)nf[1]);
! 586: return gmul((GEN)nf[8], pol_to_vec(x, N));
! 587: }
! 588: z = cgetg(N+1,t_COL);
! 589: z[1]=lcopy(x); for (i=2; i<=N; i++) z[i]=zero;
! 590: return z;
! 591: }
! 592:
! 593: GEN
! 594: algtobasis(GEN nf, GEN x)
! 595: {
! 596: long tx=typ(x),lx=lg(x),av=avma,i,N;
! 597: GEN z;
! 598:
! 599: nf=checknf(nf);
! 600: switch(tx)
! 601: {
! 602: case t_VEC: case t_COL: case t_MAT:
! 603: z=cgetg(lx,tx);
! 604: for (i=1; i<lx; i++) z[i]=(long)algtobasis(nf,(GEN)x[i]);
! 605: return z;
! 606: case t_POLMOD:
! 607: if (!polegal_spec((GEN)nf[1],(GEN)x[1]))
! 608: err(talker,"not the same number field in algtobasis");
! 609: x = (GEN)x[2]; /* fall through */
! 610: case t_POL:
! 611: return gerepileupto(av,algtobasis_intern(nf,x));
! 612:
! 613: default: N=lgef(nf[1])-3; return gscalcol(x,N);
! 614: }
! 615: }
! 616:
! 617: /* Given a and b in nf, gives an algebraic integer y in nf such that a-b.y
! 618: * is "small"
! 619: */
! 620: GEN
! 621: nfdiveuc(GEN nf, GEN a, GEN b)
! 622: {
! 623: long av=avma, tetpil;
! 624: a = element_div(nf,a,b); tetpil=avma;
! 625: return gerepile(av,tetpil,ground(a));
! 626: }
! 627:
! 628: /* Given a and b in nf, gives a "small" algebraic integer r in nf
! 629: * of the form a-b.y
! 630: */
! 631: GEN
! 632: nfmod(GEN nf, GEN a, GEN b)
! 633: {
! 634: long av=avma,tetpil;
! 635: GEN p1=gneg_i(element_mul(nf,b,ground(element_div(nf,a,b))));
! 636: tetpil=avma; return gerepile(av,tetpil,gadd(a,p1));
! 637: }
! 638:
! 639: /* Given a and b in nf, gives a two-component vector [y,r] in nf such
! 640: * that r=a-b.y is "small".
! 641: */
! 642: GEN
! 643: nfdivres(GEN nf, GEN a, GEN b)
! 644: {
! 645: long av=avma,tetpil;
! 646: GEN p1,z, y = ground(element_div(nf,a,b));
! 647:
! 648: p1=gneg_i(element_mul(nf,b,y)); tetpil=avma;
! 649: z=cgetg(3,t_VEC); z[1]=lcopy(y); z[2]=ladd(a,p1);
! 650: return gerepile(av,tetpil,z);
! 651: }
! 652:
! 653: /*************************************************************************/
! 654: /** **/
! 655: /** (Z_K/I)^* **/
! 656: /** **/
! 657: /*************************************************************************/
! 658:
! 659: /* return (column) vector of R1 signatures of x (coeff modulo 2)
! 660: * if arch = NULL, assume arch = [0,..0]
! 661: */
! 662: GEN
! 663: zsigne(GEN nf,GEN x,GEN arch)
! 664: {
! 665: GEN _0,_1, vecsign, rac = (GEN)nf[6];
! 666: long av, i,j,l;
! 667:
! 668: if (!arch) return cgetg(1,t_COL);
! 669: switch(typ(x))
! 670: {
! 671: case t_COL: x = gmul((GEN)nf[7],x); break;
! 672: case t_POLMOD: x = (GEN)x[2];
! 673: }
! 674: if (gcmp0(x)) err(talker,"zero element in zsigne");
! 675:
! 676: l = lg(arch); vecsign = cgetg(l,t_COL);
! 677: _0 = gmodulss(0,2);
! 678: _1 = gmodulss(1,2); av = avma;
! 679: for (j=1,i=1; i<l; i++)
! 680: if (signe(arch[i]))
! 681: vecsign[j++] = (gsigne(poleval(x,(GEN)rac[i])) > 0)? (long)_0
! 682: : (long)_1;
! 683: avma = av; setlg(vecsign,j); return vecsign;
! 684: }
! 685:
! 686: /* For internal use. Reduce x modulo ideal. We want a non-zero result */
! 687: GEN
! 688: nfreducemodideal(GEN nf,GEN x,GEN ideal)
! 689: {
! 690: long N = lg(x)-1, do_copy = 1, i;
! 691: GEN p1,q;
! 692:
! 693: ideal=idealhermite(nf,ideal);
! 694: for (i=N; i>=1; i--)
! 695: {
! 696: p1=gcoeff(ideal,i,i); q=gdivround((GEN)x[i],p1);
! 697: if (signe(q)) { x=gsub(x,gmul(q,(GEN)ideal[i])); do_copy=0; }
! 698: }
! 699: if (gcmp0(x)) return (GEN) ideal[1];
! 700: return do_copy? gcopy(x) : x;
! 701: }
! 702:
! 703: /* a usage interne...Reduction de la colonne x modulo la matrice y inversible
! 704: utilisant LLL */
! 705: GEN
! 706: lllreducemodmatrix(GEN x,GEN y)
! 707: {
! 708: long av = avma,tetpil;
! 709: GEN z = gmul(y,lllint(y));
! 710:
! 711: z = gneg(gmul(z, ground(gauss(z,x))));
! 712: tetpil=avma; return gerepile(av,tetpil,gadd(x,z));
! 713: }
! 714:
! 715: /* Reduce column x modulo y in HNF */
! 716: static GEN
! 717: colreducemodmat(GEN x,GEN y)
! 718: {
! 719: long av = avma, i;
! 720: GEN q;
! 721:
! 722: for (i=lg(x)-1; i>=1; i--)
! 723: {
! 724: q = gdivround((GEN)x[i], gcoeff(y,i,i));
! 725: if (signe(q)) { q = gneg(q); x = gadd(x, gmul(q,(GEN)y[i])); }
! 726: }
! 727: return avma==av? gcopy(x): gerepileupto(av,x);
! 728: }
! 729:
! 730: /* for internal use...Reduce matrix x modulo matrix y */
! 731: GEN
! 732: reducemodmatrix(GEN x, GEN y)
! 733: {
! 734: long i, lx = lg(x);
! 735: GEN z = cgetg(lx,t_MAT);
! 736:
! 737: if (DEBUGLEVEL>=8)
! 738: {
! 739: fprintferr("entering reducemodmatrix; avma-bot = %ld\n",avma-bot);
! 740: flusherr();
! 741: }
! 742: y = hnfmod(y,detint(y));
! 743: for (i=1; i<lx; i++)
! 744: {
! 745: if(DEBUGLEVEL>=8) { fprintferr("%ld ",i); flusherr(); }
! 746: z[i]=(long)colreducemodmat((GEN)x[i],y);
! 747: }
! 748: if(DEBUGLEVEL>=8) { fprintferr("\n"); flusherr(); }
! 749: return z;
! 750: }
! 751:
! 752: /* compute elt = x mod idele, with sign(elt) = sign(x) at arch */
! 753: GEN
! 754: nfreducemodidele(GEN nf,GEN x,GEN idele,GEN sarch)
! 755: {
! 756: GEN p1,p2,arch,gen;
! 757: long nba,i;
! 758:
! 759: if (gcmp0(x)) return gcopy(x);
! 760: if (!sarch || typ(idele)!=t_VEC || lg(idele)!=3)
! 761: return nfreducemodideal(nf,x,idele);
! 762:
! 763: arch =(GEN)idele[2];
! 764: nba = lg(sarch[1]);
! 765: gen =(GEN)sarch[2];
! 766: p1 = nfreducemodideal(nf,x,(GEN)idele[1]);
! 767: p2 = gadd(zsigne(nf,p1,arch), zsigne(nf,x,arch));
! 768: p2 = lift_intern(gmul((GEN)sarch[3],p2));
! 769: for (i=1; i<nba; i++)
! 770: if (signe(p2[i])) p1 = element_mul(nf,p1,(GEN)gen[i]);
! 771: return (gcmp(gnorml2(p1),gnorml2(x)) > 0)? x: p1;
! 772: }
! 773:
! 774: GEN
! 775: element_powmodideal(GEN nf,GEN x,GEN k,GEN ideal)
! 776: {
! 777: GEN y = gscalcol_i(gun,lgef(nf[1])-3);
! 778: for(;;)
! 779: {
! 780: if (mpodd(k)) y=element_mulmodideal(nf,x,y,ideal);
! 781: k=shifti(k,-1); if (!signe(k)) return y;
! 782: x = element_sqrmodideal(nf,x,ideal);
! 783: }
! 784: }
! 785:
! 786: GEN
! 787: element_powmodidele(GEN nf,GEN x,GEN k,GEN idele,GEN sarch)
! 788: {
! 789: GEN y = gscalcol_i(gun,lgef(nf[1])-3);
! 790: for(;;)
! 791: {
! 792: if (mpodd(k)) y=element_mulmodidele(nf,x,y,idele,sarch);
! 793: k=shifti(k,-1); if (!signe(k)) return y;
! 794: x = element_sqrmodidele(nf,x,idele,sarch);
! 795: }
! 796: }
! 797:
! 798: /* given 2 integral ideals x, y in HNF s.t x|y|x^2, compute the quotient
! 799: (1+x)/(1+y) in the form [[cyc],[gen],ux^-1]. */
! 800: static GEN
! 801: zidealij(GEN x, GEN y)
! 802: {
! 803: GEN p1,p2,p3,p4,d,z,x1;
! 804: long j,N,c;
! 805:
! 806: if(DEBUGLEVEL>=6)
! 807: {fprintferr("entree dans zidealij; avma-bot = %ld\n",avma-bot); flusherr();}
! 808: x1 = ginv(x);
! 809: p1 = gmul(x1,y);
! 810: if(DEBUGLEVEL>=6)
! 811: {fprintferr("p1 = y/x; avma-bot = %ld\n",avma-bot); flusherr();}
! 812: p2 = smith2(p1);
! 813: if(DEBUGLEVEL>=6)
! 814: {fprintferr("p2 = smith2(p1); avma-bot = %ld\n",avma-bot); flusherr();}
! 815: p3 = ginv((GEN)p2[1]);
! 816: if(DEBUGLEVEL>=6)
! 817: {fprintferr("p3 = 1/p2[1]; avma-bot = %ld\n",avma-bot); flusherr();}
! 818: p3 = reducemodmatrix(p3,p1);
! 819: p3 = gmul(x,p3); N=lg(p3)-1;
! 820: if(DEBUGLEVEL>=6)
! 821: {fprintferr("p3 = x.p3 mod p1; avma-bot = %ld\n",avma-bot); flusherr();}
! 822: for (j=1; j<=N; j++) coeff(p3,1,j)=laddsi(1,gcoeff(p3,1,j));
! 823: p4 = smithclean(p2); d=(GEN)p4[3];
! 824:
! 825: z=cgetg(4,t_VEC); c=lg(d); p1=cgetg(c,t_VEC);
! 826: /* transform p3 in a vector (gen) */
! 827: p3[0] = evaltyp(t_VEC) | evallg(c);
! 828: for (j=1; j<c; j++) p1[j] = coeff(d,j,j);
! 829: z[1]=(long)p1;
! 830: z[2]=(long)p3;
! 831: z[3] = lmul((GEN)p4[1],x1); return z;
! 832: }
! 833:
! 834: #if 0
! 835: /* un element g generateur d'un p^k divisant x etant donne, calcule
! 836: un element congru a g modulo p^k et a 1 modulo x/p^k et de plus
! 837: positif aux places de arch */
! 838: GEN
! 839: zconvert(GEN nf,GEN uv,GEN x,GEN arch,GEN structarch,GEN g)
! 840: {
! 841: long i,nba;
! 842: GEN p1,p2,generator;
! 843:
! 844: p1=nfreducemodideal(nf,gadd((GEN)uv[1],element_mul(nf,g,(GEN)uv[2])),x);
! 845: nba=lg(structarch[1])-1; generator=(GEN)structarch[2];
! 846: p2 = zsigne(nf,p1,arch);
! 847: for (i=1; i<=nba; i++)
! 848: if (signe(p2[i])) p1 = element_mul(nf,p1,(GEN)generator[i]);
! 849: return p1;
! 850: }
! 851: #endif
! 852:
! 853: static GEN
! 854: get_multab(GEN nf, GEN x)
! 855: {
! 856: long lx = lg(x), i;
! 857: GEN multab = cgetg(lx, t_MAT);
! 858: for (i=1; i<lx; i++)
! 859: multab[i]=(long)element_mulid(nf,x,i);
! 860: return multab;
! 861: }
! 862:
! 863: /* return mat * y mod prh */
! 864: static GEN
! 865: mul_matvec_mod_pr(GEN mat, GEN y, GEN prh)
! 866: {
! 867: long av,i,j, lx = lg(mat);
! 868: GEN v, res = cgetg(lx,t_COL), p = gcoeff(prh,1,1);
! 869:
! 870: av = avma; (void)new_chunk((lx-1) * lgefint(p));
! 871: v = zerocol(lx-1);
! 872: for (i=lx-1; i; i--)
! 873: {
! 874: GEN p1 = (GEN)v[i], t = (GEN)prh[i];
! 875: for (j=1; j<lx; j++)
! 876: p1 = addii(p1, mulii(gcoeff(mat,i,j), (GEN)y[j]));
! 877: p1 = modii(p1, p);
! 878: if (p1 != gzero && is_pm1(t[i]))
! 879: {
! 880: for (j=1; j<i; j++)
! 881: v[j] = lsubii((GEN)v[j], mulii(p1, (GEN)t[j]));
! 882: p1 = gzero;
! 883: }
! 884: if (p1 == gzero) /* intended */
! 885: res[i] = zero;
! 886: else
! 887: av = res[i] = (long)icopy_av(p1, (GEN)av);
! 888: }
! 889: avma = av; return res;
! 890: }
! 891:
! 892: /* smallest integer n such that g0^n=x modulo p prime */
! 893: static GEN
! 894: Fp_shanks(GEN x,GEN g0,GEN p)
! 895: {
! 896: long av=avma,av1,lim,lbaby,i,k,c;
! 897: GEN p1,smalltable,giant,perm,v,g0inv;
! 898:
! 899: x = modii(x,p);
! 900: if (is_pm1(x) || egalii(p,gdeux)) { avma = av; return gzero; }
! 901: p1 = addsi(-1, p);
! 902: if (egalii(p1,x)) { avma = av; return shifti(p,-1); }
! 903: p1 = racine(p1);
! 904: if (cmpis(p1,LGBITS) >= 0) err(talker,"module too large in Fp_shanks");
! 905: lbaby=itos(p1)+1; smalltable=cgetg(lbaby+1,t_VEC);
! 906: g0inv = mpinvmod(g0, p); p1 = x;
! 907:
! 908: c = 3 * lgefint(p);
! 909: for (i=1;;i++)
! 910: {
! 911: av1 = avma;
! 912: if (is_pm1(p1)) { avma=av; return stoi(i-1); }
! 913: smalltable[i]=(long)p1; if (i==lbaby) break;
! 914: new_chunk(c); p1 = mulii(p1,g0inv);
! 915: avma = av1; p1 = resii(p1, p);
! 916: }
! 917: giant = resii(mulii(x, mpinvmod(p1,p)), p);
! 918: p1=cgetg(lbaby+1,t_VEC);
! 919: perm = gen_sort(smalltable, cmp_IND | cmp_C, cmpii);
! 920: for (i=1; i<=lbaby; i++) p1[i]=smalltable[perm[i]];
! 921: smalltable=p1; p1=giant;
! 922:
! 923: av1 = avma; lim=stack_lim(av1,2);
! 924: for (k=1;;k++)
! 925: {
! 926: i=tablesearch(smalltable,p1,cmpii);
! 927: if (i)
! 928: {
! 929: v=addis(mulss(lbaby-1,k),perm[i]);
! 930: return gerepileuptoint(av,addsi(-1,v));
! 931: }
! 932: p1 = resii(mulii(p1,giant), p);
! 933:
! 934: if (low_stack(lim, stack_lim(av1,2)))
! 935: {
! 936: if(DEBUGMEM>1) err(warnmem,"Fp_shanks, k = %ld", k);
! 937: p1 = gerepileuptoint(av1, p1);
! 938: }
! 939: }
! 940: }
! 941:
! 942: /* smallest integer n such that g0^n=x modulo pr, assume g0 reduced */
! 943: /* TODO: should be done in F_(p^f), not in Z_k mod pr (done for f=1) */
! 944: GEN
! 945: nfshanks(GEN nf,GEN x,GEN g0,GEN pr,GEN prhall)
! 946: {
! 947: long av=avma,av1,lim,lbaby,i,k, f = itos((GEN)pr[4]);
! 948: GEN p1,smalltable,giant,perm,v,g0inv,prh = (GEN)prhall[1];
! 949: GEN multab, p = (GEN)pr[1];
! 950:
! 951: x = lift_intern(nfreducemodpr(nf,x,prhall));
! 952: if (f == 1) return gerepileuptoint(av, Fp_shanks((GEN)x[1],(GEN)g0[1],p));
! 953: p1 = addsi(-1, gpowgs(p,f));
! 954: if (isnfscalar(x))
! 955: {
! 956: x = (GEN)x[1];
! 957: if (gcmp1(x) || egalii((GEN)pr[1], gdeux)) { avma = av; return gzero; }
! 958: if (egalii(x, p1)) return gerepileuptoint(av,shifti(p1,-1));
! 959: v = divii(p1, addsi(-1,p));
! 960: g0 = lift_intern((GEN)element_powmodpr(nf,g0,v,prhall)[1]);
! 961: return gerepileuptoint(av, mulii(v, Fp_shanks(x,g0,p)));
! 962: }
! 963: p1 = racine(p1);
! 964: if (cmpis(p1,LGBITS) >= 0) err(talker,"module too large in nfshanks");
! 965: lbaby=itos(p1)+1; smalltable=cgetg(lbaby+1,t_VEC);
! 966: g0inv = lift_intern(element_invmodpr(nf,g0,prhall));
! 967: p1 = x;
! 968:
! 969: multab = get_multab(nf, g0inv);
! 970: for (i=lg(multab)-1; i; i--)
! 971: multab[i]=(long)Fp_vec_red((GEN)multab[i], p);
! 972:
! 973: for (i=1;;i++)
! 974: {
! 975: if (isnfscalar(p1) && gcmp1((GEN)p1[1])) { avma=av; return stoi(i-1); }
! 976:
! 977: smalltable[i]=(long)p1; if (i==lbaby) break;
! 978: p1 = mul_matvec_mod_pr(multab,p1,prh);
! 979: }
! 980: giant=lift_intern(element_divmodpr(nf,x,p1,prhall));
! 981: p1=cgetg(lbaby+1,t_VEC);
! 982: perm = gen_sort(smalltable, cmp_IND | cmp_C, cmp_vecint);
! 983: for (i=1; i<=lbaby; i++) p1[i]=smalltable[perm[i]];
! 984: smalltable=p1; p1=giant;
! 985:
! 986: multab = get_multab(nf, giant);
! 987: for (i=lg(multab)-1; i; i--)
! 988: multab[i]=(long)Fp_vec_red((GEN)multab[i], p);
! 989:
! 990: av1 = avma; lim=stack_lim(av1,2);
! 991: for (k=1;;k++)
! 992: {
! 993: i=tablesearch(smalltable,p1,cmp_vecint);
! 994: if (i)
! 995: {
! 996: v=addis(mulss(lbaby-1,k),perm[i]);
! 997: return gerepileuptoint(av,addsi(-1,v));
! 998: }
! 999: p1 = mul_matvec_mod_pr(multab,p1,prh);
! 1000:
! 1001: if (low_stack(lim, stack_lim(av1,2)))
! 1002: {
! 1003: if(DEBUGMEM>1) err(warnmem,"nfshanks");
! 1004: p1 = gerepileupto(av1, p1);
! 1005: }
! 1006: }
! 1007: }
! 1008:
! 1009: GEN
! 1010: znlog(GEN x, GEN g0)
! 1011: {
! 1012: long av=avma;
! 1013: GEN p = (GEN)g0[1];
! 1014: if (typ(g0) != t_INTMOD)
! 1015: err(talker,"not an element of (Z/pZ)* in znlog");
! 1016: switch(typ(x))
! 1017: {
! 1018: case t_INT: break;
! 1019: default: x = gmul(x, gmodulsg(1,p));
! 1020: if (typ(x) != t_INTMOD)
! 1021: err(talker,"not an element of (Z/pZ)* in znlog");
! 1022: case t_INTMOD: x = (GEN)x[2]; break;
! 1023: }
! 1024: return gerepileuptoint(av, Fp_shanks(x,(GEN)g0[2],p));
! 1025: }
! 1026:
! 1027: GEN
! 1028: dethnf(GEN mat)
! 1029: {
! 1030: long av,i,l = lg(mat);
! 1031: GEN s;
! 1032:
! 1033: if (l<3) return l<2? gun: icopy(gcoeff(mat,1,1));
! 1034: av = avma; s = gcoeff(mat,1,1);
! 1035: for (i=2; i<l; i++) s = gmul(s,gcoeff(mat,i,i));
! 1036: return av==avma? gcopy(s): gerepileupto(av,s);
! 1037: }
! 1038:
! 1039: GEN
! 1040: dethnf_i(GEN mat)
! 1041: {
! 1042: long av,i,l = lg(mat);
! 1043: GEN s;
! 1044:
! 1045: if (l<3) return l<2? gun: icopy(gcoeff(mat,1,1));
! 1046: av = avma; s = gcoeff(mat,1,1);
! 1047: for (i=2; i<l; i++) s = mulii(s,gcoeff(mat,i,i));
! 1048: return gerepileuptoint(av,s);
! 1049: }
! 1050:
! 1051: static GEN
! 1052: makeprimetoideal(GEN nf,GEN id,GEN uv,GEN x)
! 1053: {
! 1054: GEN p1 = gadd((GEN)uv[1], element_mul(nf,x,(GEN)uv[2]));
! 1055: return nfreducemodideal(nf,p1,id);
! 1056: }
! 1057:
! 1058: static GEN
! 1059: makeprimetoidealvec(GEN nf,GEN ideal,GEN uv,GEN listgen)
! 1060: {
! 1061: long i, lx = lg(listgen);
! 1062: GEN y = cgetg(lx,t_VEC);
! 1063:
! 1064: for (i=1; i<lx; i++)
! 1065: y[i] = (long)makeprimetoideal(nf,ideal,uv,(GEN)listgen[i]);
! 1066: return y;
! 1067: }
! 1068:
! 1069: /* Given an ideal pr^ep, and an integral ideal x (in HNF form) compute a list
! 1070: * of vectors, each with 5 components as follows :
! 1071: * [[clh],[gen1],[gen2],[signat2],U.X^-1]. Each component corresponds to
! 1072: * d_i,g_i,g'_i,s_i. Generators g_i are not necessarily prime to x, the
! 1073: * generators g'_i are. signat2 is the (horizontal) vector made of the
! 1074: * signatures (column vectors) of the g'_i. If x = NULL, the original ideal
! 1075: * was a prime power
! 1076: */
! 1077: static GEN
! 1078: zprimestar(GEN nf,GEN pr,GEN ep,GEN x,GEN arch)
! 1079: {
! 1080: long av=avma,av1,N,f,nbp,j,n,m,tetpil,i,e,a,b;
! 1081: GEN prh,p,pf_1,list,v,p1,p2,p3,p4,prk,uv,g0,newgen,pra,prb;
! 1082: GEN *gptr[2];
! 1083:
! 1084: if(DEBUGLEVEL>=4)
! 1085: { fprintferr("treating pr = %Z ^ %Z\n",pr,ep); flusherr(); }
! 1086: prh=prime_to_ideal(nf,pr); N=lg(prh)-1;
! 1087: f=itos((GEN)pr[4]); p=(GEN)pr[1];
! 1088:
! 1089: pf_1 = addis(gpowgs(p,f), -1);
! 1090: v = zerocol(N);
! 1091: if (f==1) v[1]=gener(p)[2];
! 1092: else
! 1093: {
! 1094: GEN prhall = cgetg(3,t_VEC);
! 1095: long psim;
! 1096: if (is_bigint(p)) err(talker,"prime too big in zprimestar");
! 1097: psim = itos(p);
! 1098: list = (GEN)factor(pf_1)[1]; nbp=lg(list)-1;
! 1099: prhall[1]=(long)prh; prhall[2]=zero;
! 1100: for (n=psim; n>=0; n++)
! 1101: {
! 1102: m=n;
! 1103: for (i=1; i<=N; i++)
! 1104: if (!gcmp1(gcoeff(prh,i,i))) { v[i]=lstoi(m%psim); m/=psim; }
! 1105: for (j=1; j<=nbp; j++)
! 1106: {
! 1107: p1 = divii(pf_1,(GEN)list[j]);
! 1108: p1 = lift_intern(element_powmodpr(nf,v,p1,prhall));
! 1109: if (isnfscalar(p1) && gcmp1((GEN)p1[1])) break;
! 1110: }
! 1111: if (j>nbp) break;
! 1112: }
! 1113: if (n < 0) err(talker,"prime too big in zprimestar");
! 1114: }
! 1115: /* v generates (Z_K / pr)^* */
! 1116: if(DEBUGLEVEL>=4) {fprintferr("v calcule\n");flusherr();}
! 1117: e = itos(ep); prk=(e==1)? pr: idealpow(nf,pr,ep);
! 1118: if(DEBUGLEVEL>=4) {fprintferr("prk calcule\n");flusherr();}
! 1119: g0 = v;
! 1120: if (x)
! 1121: {
! 1122: uv = idealaddtoone(nf,prk,idealdivexact(nf,x,prk));
! 1123: g0 = makeprimetoideal(nf,x,uv,v);
! 1124: }
! 1125: if(DEBUGLEVEL>=4) {fprintferr("g0 calcule\n");flusherr();}
! 1126:
! 1127: list=cgetg(2,t_VEC);
! 1128: p1=cgetg(6,t_VEC); list[1]=(long)p1; p1[5]=un;
! 1129: p2=cgetg(2,t_VEC); p1[1]=(long)p2; p2[1]=(long)pf_1;
! 1130: p2=cgetg(2,t_VEC); p1[2]=(long)p2; p2[1]=(long)v;
! 1131: p2=cgetg(2,t_VEC); p1[3]=(long)p2; p2[1]=(long)g0;
! 1132: p2=cgetg(2,t_VEC); p1[4]=(long)p2; p2[1]=(long)zsigne(nf,g0,arch);
! 1133: if (e==1)
! 1134: {
! 1135: tetpil=avma; return gerepile(av,tetpil,gcopy(list));
! 1136: }
! 1137:
! 1138: a=1; b=2; av1=avma;
! 1139: pra = prh; prb = (e==2)? prk: idealpow(nf,pr,gdeux);
! 1140: for(;;)
! 1141: {
! 1142: if(DEBUGLEVEL>=4)
! 1143: {fprintferr("on traite a = %ld, b = %ld\n",a,b); flusherr();}
! 1144: p1 = zidealij(pra,prb);
! 1145: newgen = dummycopy((GEN)p1[2]);
! 1146: p3 = cgetg(lg(newgen),t_VEC);
! 1147: if(DEBUGLEVEL>=4) {fprintferr("zidealij fait\n"); flusherr();}
! 1148: for (i=1; i<lg(newgen); i++)
! 1149: {
! 1150: if (x) newgen[i]=(long)makeprimetoideal(nf,x,uv,(GEN)newgen[i]);
! 1151: p3[i]=(long)zsigne(nf,(GEN)newgen[i],arch);
! 1152: }
! 1153: p2=cgetg(2,t_VEC); p4=cgetg(6,t_VEC); p2[1]=(long)p4;
! 1154: p4[1] = p1[1];
! 1155: p4[2] = p1[2];
! 1156: p4[3] = (long)newgen;
! 1157: p4[4] = (long)p3;
! 1158: p4[5] = p1[3];
! 1159:
! 1160: a=b; b=min(e,b<<1); tetpil = avma;
! 1161: list = concat(list,p2);
! 1162: if (a==b) return gerepile(av,tetpil,list);
! 1163:
! 1164: pra = gcopy(prb);
! 1165: gptr[0]=&pra; gptr[1]=&list;
! 1166: gerepilemanysp(av1,tetpil,gptr,2);
! 1167: prb = (b==(a<<1))? idealpow(nf,pra,gdeux): prk;
! 1168: }
! 1169: }
! 1170:
! 1171: /* x ideal, compute elements in 1+x whose sign matrix is invertible */
! 1172: GEN
! 1173: zarchstar(GEN nf,GEN x,GEN arch,long nba)
! 1174: {
! 1175: long av,N,i,j,k,r,rr,limr,zk,lgmat;
! 1176: GEN p1,y,bas,genarch,mat,lambda;
! 1177:
! 1178: if (!nba)
! 1179: {
! 1180: y=cgetg(4,t_VEC);
! 1181: y[1]=lgetg(1,t_VEC);
! 1182: y[2]=lgetg(1,t_VEC);
! 1183: y[3]=lgetg(1,t_MAT); return y;
! 1184: }
! 1185: N=lgef(nf[1])-3; y=cgetg(4,t_VEC);
! 1186: p1=cgetg(nba+1,t_VEC); for (i=1; i<=nba; i++) p1[i]=deux;
! 1187: y[1]=(long)p1; av = avma;
! 1188: if (N==1)
! 1189: {
! 1190: p1 = scalarpol(subsi(1, shifti(gcoeff(x,1,1),1)), varn(nf[1]));
! 1191: p1 = gerepileupto(av, p1);
! 1192: y[2]=lgetg(2,t_VEC); mael(y,2,1) = (long)p1;
! 1193: y[3]=(long)gscalmat(gun,1); return y;
! 1194: }
! 1195: zk = ideal_is_zk(x,N);
! 1196: x = gmul(x,lllintpartial(x));
! 1197: genarch = cgetg(nba+1,t_VEC);
! 1198: mat = cgetg(nba+1,t_MAT); setlg(mat,2);
! 1199: for (i=1; i<=nba; i++) mat[i] = lgetg(nba+1, t_MAT);
! 1200: lambda = new_chunk(N+1);
! 1201:
! 1202: bas = dummycopy(gmael(nf,5,1)); r = lg(arch);
! 1203: for (k=1,i=1; i<r; i++)
! 1204: if (signe(arch[i]))
! 1205: {
! 1206: if (k < i)
! 1207: for (j=1; j<=N; j++) coeff(bas,k,j) = coeff(bas,i,j);
! 1208: k++;
! 1209: }
! 1210: r = nba+1;
! 1211: for (i=1; i<=N; i++) setlg(bas[i], r);
! 1212: bas = gmul(bas, x);
! 1213:
! 1214: for (lgmat=1,r=1, rr=3; ; r<<=1, rr=(r<<1)+1)
! 1215: {
! 1216: if (DEBUGLEVEL>5) fprintferr("zarchstar: r = %ld\n",r);
! 1217: p1 = gpuigs(stoi(rr),N);
! 1218: limr = (cmpis(p1,BIGINT) > 0)? BIGINT: p1[2];
! 1219: limr = (limr-1)>>1;
! 1220: k = zk? rr: 1; /* to avoid 0 */
! 1221: for ( ; k<=limr; k++)
! 1222: {
! 1223: long av1=avma, kk = k;
! 1224: GEN alpha = gzero;
! 1225: for (i=1; i<=N; i++)
! 1226: {
! 1227: lambda[i] = (kk+r)%rr - r; kk/=rr;
! 1228: if (lambda[i])
! 1229: alpha = gadd(alpha, gmulsg(lambda[i],(GEN)bas[i]));
! 1230: }
! 1231: for (i=1; i<=nba; i++)
! 1232: alpha[i] = ladd((GEN)alpha[i], gun);
! 1233: p1 = (GEN)mat[lgmat];
! 1234: for (i=1; i<=nba; i++)
! 1235: p1[i] = (gsigne((GEN)alpha[i]) > 0)? zero: un;
! 1236:
! 1237: if (ker_trivial_mod_p(mat, gdeux)) avma = av1;
! 1238: else
! 1239: { /* new vector indep. of previous ones */
! 1240: avma = av1; alpha = gzero;
! 1241: for (i=1; i<=N; i++)
! 1242: if (lambda[i])
! 1243: alpha = gadd(alpha, gmulsg(lambda[i],(GEN)x[i]));
! 1244: alpha[1] = ladd((GEN)alpha[1], gun);
! 1245: genarch[lgmat++] = (long)alpha;
! 1246: if (lgmat > nba)
! 1247: {
! 1248: GEN *gptr[2];
! 1249: mat = ginv(mat); gptr[0]=&genarch; gptr[1]=&mat;
! 1250: gerepilemany(av,gptr,2);
! 1251: y[2]=(long)genarch;
! 1252: y[3]=(long)mat; return y;
! 1253: }
! 1254: setlg(mat,lgmat+1);
! 1255: }
! 1256: }
! 1257: }
! 1258: }
! 1259:
! 1260: /* Retourne la decomposition de a sur les nbgen generateurs successifs
! 1261: * contenus dans list_set et si index !=0 on ne fait ce calcul que pour
! 1262: * l'ideal premier correspondant a cet index en mettant a 0 les autres
! 1263: * composantes
! 1264: */
! 1265: static GEN
! 1266: zinternallog(GEN nf,GEN list_set,long nbgen,GEN arch,GEN fa,GEN a,long index)
! 1267: {
! 1268: GEN prlist,ep,y,ainit,list,pr,prk,cyc,gen,psigne,p1,p2,p3;
! 1269: long av,nbp,cp,i,j,k;
! 1270:
! 1271: y = cgetg(nbgen+1,t_COL); cp=0; av=avma;
! 1272: prlist=(GEN)fa[1]; ep=(GEN)fa[2]; nbp=lg(ep)-1;
! 1273: i=typ(a); if (is_extscalar_t(i)) a = algtobasis(nf,a);
! 1274: if (DEBUGLEVEL>3)
! 1275: {
! 1276: fprintferr("entering zinternallog, "); flusherr();
! 1277: if (DEBUGLEVEL>5) fprintferr("with a = %Z\n",a);
! 1278: }
! 1279: ainit = a; psigne = zsigne(nf,ainit,arch);
! 1280: for (k=1; k<=nbp; k++)
! 1281: {
! 1282: list=(GEN)list_set[k];
! 1283: if (index && index!=k)
! 1284: {
! 1285: for (j=1; j<lg(list); j++)
! 1286: {
! 1287: cyc = gmael(list,j,1);
! 1288: for (i=1; i<lg(cyc); i++) y[++cp]=zero;
! 1289: }
! 1290: continue;
! 1291: }
! 1292: pr=(GEN)prlist[k]; prk=idealpow(nf,pr,(GEN)ep[k]);
! 1293: for (j=1; j<lg(list); j++)
! 1294: {
! 1295: p1 = (GEN)list[j]; cyc=(GEN)p1[1]; gen=(GEN)p1[2];
! 1296: if (j==1)
! 1297: {
! 1298: if (DEBUGLEVEL>5) { fprintferr("do nfshanks\n"); flusherr(); }
! 1299: a=ainit; p3=nfmodprinit(nf,pr);
! 1300: p3 = nfshanks(nf,a,(GEN)gen[1],pr,p3);
! 1301: }
! 1302: else
! 1303: {
! 1304: p3 = (GEN)a[1]; a[1] = laddsi(-1,(GEN)a[1]);
! 1305: p2 = gmul((GEN)p1[5],a); a[1] = (long)p3;
! 1306: if (lg(p2)!=lg(cyc)) err(bugparier,"zinternallog");
! 1307: p3 = (GEN)p2[1];
! 1308: }
! 1309: for(i=1;;)
! 1310: {
! 1311: p3 = modii(negi(p3), (GEN)cyc[i]);
! 1312: y[++cp] = lnegi(p3);
! 1313: if (signe(p3))
! 1314: {
! 1315: if (mpodd((GEN)y[cp])) psigne = gadd(psigne,gmael(p1,4,i));
! 1316: if (DEBUGLEVEL>5) fprintferr("do element_powmodideal\n");
! 1317: p3 = element_powmodideal(nf,(GEN)gen[i],p3,prk);
! 1318: a = element_mulmodideal(nf,a,p3,prk);
! 1319: }
! 1320: i++; if (i==lg(cyc)) break;
! 1321: p3 = (GEN)p2[i];
! 1322: }
! 1323: }
! 1324: }
! 1325: p1=lift_intern(gmul(gmael(list_set,nbp+1,3), psigne));
! 1326: avma=av; for (i=1; i<lg(p1); i++) y[++cp] = p1[i];
! 1327: if (DEBUGLEVEL>3) { fprintferr("leaving\n"); flusherr(); }
! 1328: for (i=1; i<=nbgen; i++) y[i] = licopy((GEN)y[i]);
! 1329: return y;
! 1330: }
! 1331:
! 1332: GEN
! 1333: compute_gen(GEN nf, GEN u1, GEN met, GEN gen, GEN module, long nbp, GEN sarch)
! 1334: {
! 1335: long i,j,s,nba, c = lg(met), lgen = lg(gen);
! 1336: GEN basecl=cgetg(c,t_VEC), unnf = gscalcol_i(gun, lgef(nf[1])-3);
! 1337: GEN arch,x,p1,p2,p3,p4,p5,generator;
! 1338:
! 1339: if (sarch)
! 1340: {
! 1341: arch = (GEN)module[2];
! 1342: x = (GEN)module[1];
! 1343: generator = (GEN)sarch[2];
! 1344: nba = lg(generator)-1;
! 1345: }
! 1346: else
! 1347: {
! 1348: x = (typ(module) == t_MAT)? module: (GEN)module[1];
! 1349: nba = 0;
! 1350: }
! 1351: for (j=1; j<c; j++)
! 1352: {
! 1353: GEN *op, plus = unnf, minus = unnf;
! 1354:
! 1355: for (i=1; i<lgen; i++)
! 1356: {
! 1357: p1 = gcoeff(u1,i,j);
! 1358: if (!(s = signe(p1))) continue;
! 1359:
! 1360: if (s > 0) op = + else { op = − p1 = negi(p1); }
! 1361: p1 = element_powmodidele(nf,(GEN)gen[i],p1,module,sarch);
! 1362: *op = *op==unnf? p1: element_mulmodidele(nf,*op,p1,module,sarch);
! 1363: }
! 1364: if (nbp)
! 1365: {
! 1366: p5=idealaddtoone_i(nf,minus,x);
! 1367: p4=element_div(nf,p5,minus); /* = minus^(-1) mod x */
! 1368: p3=element_mulmodideal(nf,plus,p4,x);
! 1369: }
! 1370: else p3=unnf;
! 1371:
! 1372: if (nba)
! 1373: { /* correct so that sign(p3) == sign(plus/minus) */
! 1374: p4=gadd(gadd(zsigne(nf,p3,arch),
! 1375: zsigne(nf,plus,arch)),
! 1376: zsigne(nf,minus,arch));
! 1377: p2=lift_intern(gmul((GEN)sarch[3],p4));
! 1378: for (i=1; i<=nba; i++)
! 1379: if (signe(p2[i])) p3=element_mul(nf,p3,(GEN)generator[i]);
! 1380: }
! 1381: basecl[j]=(long)p3;
! 1382: }
! 1383: return basecl;
! 1384: }
! 1385:
! 1386: /* Compute [[ideal,arch], [h,[cyc],[gen]], idealfact, [liste], U]
! 1387: gen not computed unless add_gen != 0 */
! 1388: GEN
! 1389: zidealstarinitall(GEN nf, GEN ideal,long add_gen)
! 1390: {
! 1391: long av=avma,tetpil,i,j,k,nba,nbp,R1,nbgen,cp;
! 1392: GEN p1,p2,p3,p4,y,h,met,u1,u1u2,u1u2clean,fa,fa2,ep,x,arch,gen,list,sarch;
! 1393:
! 1394: nf = checknf(nf); R1=itos(gmael(nf,2,1));
! 1395: if (typ(ideal)==t_VEC && lg(ideal)==3)
! 1396: {
! 1397: arch=(GEN)ideal[2]; ideal = (GEN)ideal[1];
! 1398: i = typ(arch);
! 1399: if (!is_vec_t(i) || lg(arch) != R1+1)
! 1400: err(talker,"incorrect archimedean component in zidealstarinit");
! 1401: for (nba=0,i=1; i<=R1; i++)
! 1402: if (signe(arch[i])) nba++;
! 1403: }
! 1404: else
! 1405: {
! 1406: arch=cgetg(R1+1,t_VEC);
! 1407: for (nba=0,i=1; i<=R1; i++) arch[i]=zero;
! 1408: }
! 1409: x = idealhermite(nf,ideal);
! 1410: if (!gcmp1(denom(x)))
! 1411: err(talker,"zidealstarinit needs an integral ideal. x =\n%Z\n",x);
! 1412: p1=cgetg(3,t_VEC); ideal=p1;
! 1413: p1[1]=(long)x;
! 1414: p1[2]=(long)arch;
! 1415:
! 1416: fa=idealfactor(nf,x); list=(GEN)fa[1]; ep=(GEN)fa[2];
! 1417: nbp=lg(list)-1; fa2=cgetg(nbp+2,t_VEC);
! 1418:
! 1419: /* compute them */
! 1420: gen = cgetg(1,t_VEC);
! 1421: p2 = (nbp==1)? (GEN)NULL: x;
! 1422: for (i=1; i<=nbp; i++)
! 1423: {
! 1424: p1 = zprimestar(nf,(GEN)list[i],(GEN)ep[i],p2,arch);
! 1425: fa2[i]=(long)p1;
! 1426: for (j=1; j<lg(p1); j++)
! 1427: gen = concatsp(gen,gmael(p1,j,3));
! 1428: }
! 1429: sarch = zarchstar(nf,x,arch,nba);
! 1430: fa2[i]=(long)sarch;
! 1431: gen = concatsp(gen,(GEN)sarch[2]);
! 1432: nbgen = lg(gen)-1;
! 1433: h=cgetg(nbgen+1,t_MAT); cp=0;
! 1434: for (i=1; i<=nbp; i++)
! 1435: {
! 1436: list=(GEN)fa2[i];
! 1437: for (j=1; j<lg(list); j++)
! 1438: {
! 1439: p1=(GEN)list[j]; p2=(GEN)p1[1]; p3=(GEN)p1[3];
! 1440: for (k=1; k<lg(p3); k++)
! 1441: {
! 1442: if (DEBUGLEVEL>=6)
! 1443: { fprintferr("entering element_powmodidele\n"); flusherr(); }
! 1444: p4 = element_powmodidele(nf,(GEN)p3[k],(GEN)p2[k],ideal,sarch);
! 1445: h[++cp] = lneg(zinternallog(nf,fa2,nbgen,arch,fa,p4,i));
! 1446: coeff(h,cp,cp) = p2[k];
! 1447: }
! 1448: }
! 1449: }
! 1450: for (j=1; j<=nba; j++)
! 1451: { h[++cp]=(long)zerocol(nbgen); coeff(h,cp,cp)=deux; }
! 1452: if (cp!=nbgen) err(talker,"bug in zidealstarinit");
! 1453: u1u2 = smith2(h); u1u2clean = smithclean(u1u2);
! 1454: met=(GEN)u1u2clean[3];
! 1455: if (add_gen)
! 1456: {
! 1457: u1 = reducemodmatrix(ginv((GEN)u1u2[1]),h);
! 1458: p1 = cgetg(4,t_VEC);
! 1459: p1[3] = (long)compute_gen(nf,u1,met,gen,ideal, nbp,sarch);
! 1460: }
! 1461: else p1 = cgetg(3,t_VEC);
! 1462: y = cgetg(6,t_VEC);
! 1463: y[1]=(long)ideal;
! 1464: y[2]=(long)p1;
! 1465: p1[1]=(long)dethnf(met);
! 1466: p1[2]=(long)mattodiagonal(met);
! 1467: y[3]=(long)fa;
! 1468: y[4]=(long)fa2;
! 1469: y[5] = u1u2clean[1];
! 1470: tetpil=avma; return gerepile(av,tetpil,gcopy(y));
! 1471: }
! 1472:
! 1473: GEN
! 1474: zidealstarinitgen(GEN nf, GEN ideal)
! 1475: {
! 1476: return zidealstarinitall(nf,ideal,1);
! 1477: }
! 1478:
! 1479: GEN
! 1480: zidealstarinit(GEN nf, GEN ideal)
! 1481: {
! 1482: return zidealstarinitall(nf,ideal,0);
! 1483: }
! 1484:
! 1485: GEN
! 1486: zidealstar(GEN nf, GEN ideal)
! 1487: {
! 1488: long av = avma,tetpil;
! 1489: GEN y = zidealstarinitall(nf,ideal,1);
! 1490: tetpil=avma; return gerepile(av,tetpil,gcopy((GEN)y[2]));
! 1491: }
! 1492:
! 1493: GEN
! 1494: idealstar0(GEN nf, GEN ideal,long flag)
! 1495: {
! 1496: switch(flag)
! 1497: {
! 1498: case 0: return zidealstar(nf,ideal);
! 1499: case 1: return zidealstarinit(nf,ideal);
! 1500: case 2: return zidealstarinitgen(nf,ideal);
! 1501: default: err(flagerr,"idealstar");
! 1502: }
! 1503: return NULL; /* not reached */
! 1504: }
! 1505:
! 1506: /* x is not integral, but check that v_p(x)=0 for all prime divisors of the
! 1507: * ideal. We need x = a/b with integral a and b, prime to ideal
! 1508: * denmat = den * id.
! 1509: */
! 1510: static GEN
! 1511: rat_zinternallog(GEN nf, GEN x, GEN bigideal, GEN denmat)
! 1512: {
! 1513: long nbp,i,v,k, N = lgef(nf[1])-3;
! 1514: GEN den,fa,list,ep,pr,p1,p2,p3,x1,dinv,ideal;
! 1515:
! 1516: ideal = (GEN)bigideal[1];
! 1517: if (lg(ideal) == 3) ideal = (GEN)ideal[1];
! 1518: fa=(GEN)bigideal[3]; list=(GEN)fa[1]; ep=(GEN)fa[2];
! 1519: den=gmael(denmat,1,1); k=1; nbp=lg(list)-1;
! 1520: for (i=1; i<=nbp; i++)
! 1521: {
! 1522: pr=(GEN)list[i];
! 1523: v = (ggval(den,(GEN)pr[1])*itos((GEN)pr[3])) / itos((GEN)ep[i]) + 1;
! 1524: if (v>k) k=v;
! 1525: }
! 1526: p3=idealpow(nf,ideal,stoi(k));
! 1527: p1=idealadd(nf,denmat,p3); dinv=idealinv(nf,p1);
! 1528: p2=idealmullll(nf,denmat,dinv);
! 1529: p3=idealmullll(nf,p3,dinv);
! 1530: x1=idealaddtoone_i(nf,p2,p3);
! 1531: if (gcmp0(x1)) x1 = (GEN)denmat[1];
! 1532: /* x = a/b is suitable, with a=x1*x, b=x1 */
! 1533: p1=element_mul(nf,x1,x);
! 1534: /* x1 is prime to ideal. Check that x1 * x also */
! 1535: p2=idealadd(nf,p1,ideal);
! 1536: if (! ideal_is_zk(p2,N))
! 1537: err(talker,"element is not coprime to ideal in zideallog");
! 1538: p1=zideallog(nf,p1,bigideal);
! 1539: p2=zideallog(nf,x1,bigideal);
! 1540: return gsub(p1,p2);
! 1541: }
! 1542:
! 1543: /* Given x (not necessarily integral), and bid as output by zidealstarinit,
! 1544: * compute the vector of components on the generators bid[2]
! 1545: */
! 1546: GEN
! 1547: zideallog(GEN nf, GEN x, GEN bid)
! 1548: {
! 1549: long av,l,i,N,c;
! 1550: GEN fa,fa2,ideal,arch,den,p1,cyc,y;
! 1551:
! 1552: nf=checknf(nf); checkbid(bid);
! 1553: cyc=gmael(bid,2,2); c=lg(cyc);
! 1554: y=cgetg(c,t_COL); av=avma;
! 1555: N = lgef(nf[1])-3; ideal = (GEN) bid[1];
! 1556: if (typ(ideal)==t_VEC && lg(ideal)==3)
! 1557: arch = (GEN)ideal[2];
! 1558: else
! 1559: arch = NULL;
! 1560: switch(typ(x))
! 1561: {
! 1562: case t_INT: case t_FRAC: case t_FRACN:
! 1563: x = gscalcol_i(x,N); break;
! 1564: case t_POLMOD: case t_POL:
! 1565: x = algtobasis(nf,x); break;
! 1566: case t_COL: break;
! 1567: default: err(talker,"not an element in zideallog");
! 1568: }
! 1569: if (lg(x) != N+1) err(talker,"not an element in zideallog");
! 1570:
! 1571: den=denom(x);
! 1572: if (!gcmp1(den))
! 1573: p1 = rat_zinternallog(nf,x,bid, gscalmat(den,N));
! 1574: else
! 1575: {
! 1576: l=lg(bid[5])-1; fa=(GEN)bid[3]; fa2=(GEN)bid[4];
! 1577: p1 = zinternallog(nf,fa2,l,arch,fa,x,0);
! 1578: p1 = gmul((GEN)bid[5],p1); /* apply smith */
! 1579: }
! 1580: if (lg(p1)!=c) err(bugparier,"zideallog");
! 1581: for (i=1; i<c; i++)
! 1582: y[i] = lmodii((GEN)p1[i],(GEN)cyc[i]);
! 1583: avma=av; /* following line does a gerepile ! */
! 1584: for (i=1; i<c; i++)
! 1585: y[i] = (long)icopy((GEN)y[i]);
! 1586: return y;
! 1587: }
! 1588:
! 1589: /* Etant donnes bid1, bid2 resultats de zidealstarinit pour deux modules m1
! 1590: * et m2 premiers entre eux sans partie archimedienne, calcule le
! 1591: * zidealstarinit [[ideal,arch],[h,[cyc],[gen]],idealfact,[liste],U] du
! 1592: * produit
! 1593: */
! 1594: static GEN
! 1595: zidealstarinitjoinall(GEN nf, GEN bid1, GEN bid2, long add_gen)
! 1596: {
! 1597: long av=avma,i,j,nbp,nbgen,lx1,lx2,llx1,llx2,lx,llx;
! 1598: GEN module1,module2,struct1,struct2,fact1,fact2,liste1,liste2,U1,U2;
! 1599: GEN module,liste,fact,U,cyc,ex1,ex2,uv;
! 1600: GEN p1,p2,y,met,u1;
! 1601: GEN fa1,fa2,gen,u1u2,u1u2clean;
! 1602:
! 1603: nf=checknf(nf); checkbid(bid1); checkbid(bid2);
! 1604: module1=(GEN)bid1[1]; struct1=(GEN)bid1[2]; fact1=(GEN)bid1[3];
! 1605: module2=(GEN)bid2[1]; struct2=(GEN)bid2[2]; fact2=(GEN)bid2[3];
! 1606: module=cgetg(3,t_VEC);
! 1607: module[1]=(long)idealmul(nf,(GEN)module1[1],(GEN)module2[1]);
! 1608: module[2]=ladd((GEN)module1[2],(GEN)module2[2]);
! 1609: if (gcmpgs(vecmax((GEN)module[2]),1)>=0)
! 1610: err(talker,"nontrivial Archimedian components in zidealstarinitjoin");
! 1611:
! 1612: fa1=(GEN)fact1[1]; ex1=(GEN)fact1[2];
! 1613: fa2=(GEN)fact2[1]; ex2=(GEN)fact2[2];
! 1614: fact=cgetg(3,t_MAT);
! 1615: fact[1]=lconcat(fa1,fa2); lx1=lg(fa1);
! 1616: fact[2]=lconcat(ex1,ex2); lx2=lg(fa2);
! 1617: nbp=lx1+lx2-2;
! 1618: for (i=1; i<lx1; i++)
! 1619: if (isinvector(fa2,(GEN)fa1[i],lx2-1))
! 1620: err(talker,"noncoprime ideals in zidealstarinitjoin");
! 1621:
! 1622: liste1=(GEN)bid1[4]; lx1=lg(liste1);
! 1623: liste2=(GEN)bid2[4]; lx2=lg(liste2);
! 1624: lx=lx1+lx2-2; liste = cgetg(lx,t_VEC);
! 1625: for (i=1; i<lx1-1; i++) liste[i]=liste1[i];
! 1626: for ( ; i<lx; i++) liste[i]=liste2[i-lx1+2];
! 1627: U1=(GEN)bid1[5]; lx1=lg(U1);
! 1628: U2=(GEN)bid2[5]; lx2=lg(U2);
! 1629: lx=lx1+lx2-1; U = cgetg(lx,t_MAT);
! 1630:
! 1631: llx1=lg(struct1[2]);
! 1632: llx2=lg(struct2[2]);
! 1633: llx=llx1+llx2-1; nbgen=llx-1;
! 1634: cyc = diagonal(concatsp((GEN)struct1[2],(GEN)struct2[2]));
! 1635: u1u2=smith2(cyc); u1u2clean=smithclean(u1u2);
! 1636: met=(GEN)u1u2clean[3];
! 1637: if (nbgen)
! 1638: {
! 1639: for (j=1; j<lx1; j++)
! 1640: {
! 1641: p1=cgetg(llx,t_COL); U[j]=(long)p1;
! 1642: p2=(GEN)U1[j];
! 1643: for (i=1; i<llx1; i++) p1[i]=p2[i];
! 1644: for ( ; i<llx; i++) p1[i]=zero;
! 1645: }
! 1646: for ( ; j<lx; j++)
! 1647: {
! 1648: p1=cgetg(llx,t_COL); U[j]=(long)p1;
! 1649: p2=((GEN)U2[j-(lx1-1)]) + 1-llx1;
! 1650: for (i=1; i<llx1; i++) p1[i]=zero;
! 1651: for ( ; i<llx; i++) p1[i]=p2[i];
! 1652: }
! 1653: U = gmul((GEN)u1u2clean[1],U);
! 1654: }
! 1655: else
! 1656: for (j=1; j<lx; j++) U[j]=lgetg(1,t_COL);
! 1657:
! 1658: if (add_gen)
! 1659: {
! 1660: if (lg(struct1)<=3 || lg(struct2)<=3)
! 1661: err(talker,"please apply idealstar(,,2) and not idealstar(,,1)");
! 1662:
! 1663: uv = idealaddtoone(nf,(GEN)module1[1],(GEN)module2[1]);
! 1664: p1 = makeprimetoidealvec(nf,(GEN)module[1],uv,(GEN)struct1[3]);
! 1665: p2=(GEN)uv[1]; uv[1]=uv[2]; uv[2]=(long)p2;
! 1666: p2 = makeprimetoidealvec(nf,(GEN)module[1],uv,(GEN)struct2[3]);
! 1667: gen=concatsp(p1,p2);
! 1668:
! 1669: u1 = reducemodmatrix(ginv((GEN)u1u2[1]),cyc);
! 1670: p1 = cgetg(4,t_VEC);
! 1671: p1[3] = (long)compute_gen(nf,u1,met,gen,module, nbp,NULL);
! 1672: }
! 1673: else p1 = cgetg(3,t_VEC);
! 1674:
! 1675: y=cgetg(6,t_VEC);
! 1676: y[1]=(long)module;
! 1677: y[2]=(long)p1;
! 1678: p1[1]=(long)dethnf(met);
! 1679: p1[2]=(long)mattodiagonal(met);
! 1680: y[3]=(long)fact;
! 1681: y[4]=(long)liste;
! 1682: y[5]=(long)U;
! 1683: return gerepileupto(av,gcopy(y));
! 1684: }
! 1685:
! 1686: GEN
! 1687: zidealstarinitjoin(GEN nf, GEN bid1, GEN bid2)
! 1688: {
! 1689: return zidealstarinitjoinall(nf,bid1,bid2,0);
! 1690: }
! 1691:
! 1692: GEN
! 1693: zidealstarinitjoingen(GEN nf, GEN bid1, GEN bid2)
! 1694: {
! 1695: return zidealstarinitjoinall(nf,bid1,bid2,1);
! 1696: }
! 1697:
! 1698: /* Etant donnes bid1 resultat de zidealstarinit pour un module m1 sans partie
! 1699: * archimedienne et une partie archimedienne arch, calcule le zidealstarinit
! 1700: * [[ideal,arch],[h,[cyc],[gen]],idealfact,[liste],U] du produit
! 1701: */
! 1702: static GEN
! 1703: zidealstarinitjoinarchall(GEN nf, GEN bid1, GEN arch, long nba, long add_gen)
! 1704: {
! 1705: long av=avma,i,nbp,lx1,lx;
! 1706: GEN module1,struct1,fact1,liste1,U1;
! 1707: GEN module,liste,h,p1,y,met,u1,x,gen,sarch,u1u2,u1u2clean;
! 1708:
! 1709: nf=checknf(nf); checkbid(bid1);
! 1710: module1=(GEN)bid1[1]; struct1=(GEN)bid1[2]; fact1=(GEN)bid1[3];
! 1711: nbp=lg((GEN)fact1[1])-1;
! 1712: x=(GEN)module1[1];
! 1713: sarch = zarchstar(nf,x,arch,nba);
! 1714: module=cgetg(3,t_VEC);
! 1715: module[1]=(long)x;
! 1716: module[2]=(long)arch;
! 1717: if (gcmpgs(vecmax((GEN)module1[2]),1)>=0)
! 1718: err(talker,"nontrivial Archimedian components in zidealstarinitjoinarchall");
! 1719: liste1=(GEN)bid1[4]; lx=lg(liste1);
! 1720: liste=cgetg(lx,t_VEC); for (i=1; i<lx-1; i++) liste[i]=liste1[i];
! 1721: liste[lx-1]=(long)sarch;
! 1722: U1=(GEN)bid1[5]; lx1=lg(U1);
! 1723: lx=lx1+nba;
! 1724: h = diagonal(concatsp((GEN)struct1[2],(GEN)sarch[1]));
! 1725: u1u2 = smith2(h); u1u2clean = smithclean(u1u2);
! 1726: met=(GEN)u1u2clean[3];
! 1727:
! 1728: if (add_gen)
! 1729: {
! 1730: if (lg(struct1)<=3)
! 1731: err(talker,"please apply idealstar(,,2) and not idealstar(,,1)");
! 1732: gen=concatsp((GEN)struct1[3],(GEN)sarch[2]);
! 1733:
! 1734: u1 = reducemodmatrix(ginv((GEN)u1u2[1]),h);
! 1735: p1 = cgetg(4,t_VEC);
! 1736: p1[3] = (long)compute_gen(nf,u1,met,gen,module, nbp,sarch);
! 1737: }
! 1738: else p1 = cgetg(3,t_VEC);
! 1739:
! 1740: y=cgetg(6,t_VEC);
! 1741: y[1]=(long)module;
! 1742: y[2]=(long)p1;
! 1743: p1[1]=(long)dethnf(met);
! 1744: p1[2]=(long)mattodiagonal(met);
! 1745: y[3]=(long)fact1;
! 1746: y[4]=(long)liste;
! 1747: y[5] = u1u2clean[1];
! 1748: return gerepileupto(av,gcopy(y));
! 1749: }
! 1750:
! 1751: GEN
! 1752: zidealstarinitjoinarch(GEN nf, GEN bid1, GEN arch, long nba)
! 1753: {
! 1754: return zidealstarinitjoinarchall(nf,bid1,arch,nba,0);
! 1755: }
! 1756:
! 1757: GEN
! 1758: zidealstarinitjoinarchgen(GEN nf, GEN bid1, GEN arch, long nba)
! 1759: {
! 1760: return zidealstarinitjoinarchall(nf,bid1,arch,nba,1);
! 1761: }
! 1762:
! 1763: /* calcule la matrice des zinternallog des unites */
! 1764: GEN
! 1765: logunitmatrix(GEN nf,GEN funits,GEN racunit,GEN bid)
! 1766: {
! 1767: long R,j,sizeh;
! 1768: GEN m,fa2,fa,arch;
! 1769:
! 1770: R=lg(funits)-1; m=cgetg(R+2,t_MAT);
! 1771: fa2=(GEN)bid[4]; sizeh=lg(bid[5])-1; arch=gmael(bid,1,2);
! 1772: fa=(GEN)bid[3];
! 1773: m[1]=(long)zinternallog(nf,fa2,sizeh,arch,fa,racunit,0);
! 1774: for (j=2; j<=R+1; j++)
! 1775: m[j]=(long)zinternallog(nf,fa2,sizeh,arch,fa,(GEN)funits[j-1],0);
! 1776: return m;
! 1777: }
! 1778:
! 1779: /* calcule la matrice des zinternallog des unites */
! 1780: static GEN
! 1781: logunitmatrixarch(GEN nf,GEN funits,GEN racunit,GEN bid)
! 1782: {
! 1783: long R,j;
! 1784: GEN m,liste,structarch,arch;
! 1785:
! 1786: R=lg(funits)-1; m=cgetg(R+2,t_MAT); arch=gmael(bid,1,2);
! 1787: liste=(GEN)bid[4]; structarch=(GEN)liste[lg(liste)-1];
! 1788: m[1]=(long)zsigne(nf,racunit,arch);
! 1789: for (j=2; j<=R+1; j++)
! 1790: m[j]=(long)zsigne(nf,(GEN)funits[j-1],arch);
! 1791: return lift_intern(gmul((GEN)structarch[3],m));
! 1792: }
! 1793:
! 1794: /* concatenation verticale de Q1 et Q2. Ne cree pas le resultat. */
! 1795: GEN
! 1796: vconcat(GEN Q1, GEN Q2)
! 1797: {
! 1798: long lc,lr,lx1,lx2,i,j;
! 1799: GEN m,p1,p2,p3;
! 1800:
! 1801: lc=lg(Q1); if (lc==1) return Q1;
! 1802: lx1=lg(Q1[1]); lx2=lg(Q2[1]); lr=lx1+lx2-1;
! 1803: m=cgetg(lc,t_MAT);
! 1804: for (j=1; j<lc; j++)
! 1805: {
! 1806: p1=cgetg(lr,t_COL); m[j]=(long)p1; p2=(GEN)Q1[j]; p3=(GEN)Q2[j];
! 1807: for (i=1; i<lx1; i++) p1[i]=p2[i];
! 1808: for ( ; i<lr; i++) p1[i]=p3[i-lx1+1];
! 1809: }
! 1810: return m;
! 1811: }
! 1812:
! 1813: static void
! 1814: init_units(GEN bnf, GEN *funits, GEN *racunit)
! 1815: {
! 1816: GEN p1;
! 1817: bnf = checkbnf(bnf); p1=(GEN)bnf[8];
! 1818: if (lg(p1)==5) *funits=(GEN)buchfu(bnf)[1];
! 1819: else
! 1820: {
! 1821: if (lg(p1)!=7) err(talker,"incorrect big number field");
! 1822: *funits=(GEN)p1[5];
! 1823: }
! 1824: *racunit=gmael(p1,4,2);
! 1825: }
! 1826:
! 1827: /* flag &1 : generateurs, sinon non
! 1828: * flag &2 : unites, sinon pas.
! 1829: * flag &4 : ideaux, sinon zidealstar.
! 1830: */
! 1831: static GEN
! 1832: ideallistzstarall(GEN bnf,long bound,long flag)
! 1833: {
! 1834: byteptr ptdif=diffptr;
! 1835: long lim,av0=avma,av,i,j,k,l,q2,lp1,q;
! 1836: long do_gen = flag & 1, do_units = flag & 2, big_id = !(flag & 4);
! 1837: GEN y,nf,p,z,z2,p1,p2,p3,fa,pr,ideal,lu,lu2,funits,racunit,embunit;
! 1838:
! 1839: nf=checknf(bnf);
! 1840: z = cgetg(bound+1,t_VEC);
! 1841: for (i=2; i<=bound; i++) z[i] = lgetg(1,t_VEC);
! 1842: if (do_units)
! 1843: {
! 1844: init_units(bnf,&funits,&racunit);
! 1845: lu = cgetg(bound+1,t_VEC);
! 1846: for (i=2; i<=bound; i++) lu[i]=lgetg(1,t_VEC);
! 1847: }
! 1848:
! 1849: ideal = idmat(lgef(nf[1])-3);
! 1850: if (big_id) ideal = zidealstarinitall(nf,ideal,do_gen);
! 1851: p2 = cgetg(2,t_VEC); z[1] = (long)p2;
! 1852: p2[1] = (long)ideal;
! 1853: if (do_units)
! 1854: {
! 1855: p2 = cgetg(2,t_VEC); lu[1] = (long)p2;
! 1856: p2[1] = (long)logunitmatrix(nf,funits,racunit,ideal);
! 1857: }
! 1858:
! 1859: p=cgeti(3); p[1]=evalsigne(1) | evallgefint(3);
! 1860: av=avma; lim=stack_lim(av,1);
! 1861: for (p[2]=0; p[2]<=bound; )
! 1862: {
! 1863: if (!*ptdif) err(primer1);
! 1864: p[2] += *ptdif++;
! 1865: if (DEBUGLEVEL>1) { fprintferr("%ld ",p[2]); flusherr(); }
! 1866: fa = primedec(nf,p);
! 1867: for (j=1; j<lg(fa); j++)
! 1868: {
! 1869: pr = (GEN)fa[j]; p1 = powgi(p,(GEN)pr[4]);
! 1870: if (is_bigint(p1) || (q = itos(p1)) > bound) continue;
! 1871:
! 1872: q2=q; ideal=pr; z2=dummycopy(z);
! 1873: if (do_units) lu2=dummycopy(lu);
! 1874: for (l=2; ;l++)
! 1875: {
! 1876: if (big_id) ideal = zidealstarinitall(nf,ideal,do_gen);
! 1877: if (do_units) embunit = logunitmatrix(nf,funits,racunit,ideal);
! 1878: for (i=q; i<=bound; i+=q)
! 1879: {
! 1880: p1 = (GEN)z[i/q];
! 1881: if ((lp1 = lg(p1)) == 1) continue;
! 1882:
! 1883: p2 = cgetg(lp1,t_VEC);
! 1884: for (k=1; k<lp1; k++)
! 1885: if (big_id)
! 1886: p2[k] = (long)zidealstarinitjoinall(nf,(GEN)p1[k],ideal,do_gen);
! 1887: else
! 1888: p2[k] = (long)idealmul(nf,(GEN)p1[k],ideal);
! 1889: z2[i] = (long)concatsp((GEN)z2[i],p2);
! 1890: if (do_units)
! 1891: {
! 1892: p1 = (GEN)lu[i/q];
! 1893: p2 = cgetg(lp1,t_VEC);
! 1894: for (k=1; k<lp1; k++)
! 1895: p2[k] = (long)vconcat((GEN)p1[k],embunit);
! 1896: lu2[i] = (long)concatsp((GEN)lu2[i],p2);
! 1897: }
! 1898: }
! 1899: q *= q2; if ((ulong)q > (ulong)bound) break;
! 1900: ideal = idealpows(nf,pr,l);
! 1901: }
! 1902: z = z2; if (do_units) lu = lu2;
! 1903: }
! 1904: if (low_stack(lim, stack_lim(av,1)))
! 1905: {
! 1906: GEN *gptr[2]; gptr[0]=&z; gptr[1]=&lu;
! 1907: if(DEBUGMEM>1) err(warnmem,"ideallistzstarall");
! 1908: gerepilemany(av,gptr,do_units?2:1);
! 1909: }
! 1910: }
! 1911: if (!do_units) return gerepileupto(av0, gcopy(z));
! 1912: y = cgetg(3,t_VEC);
! 1913: y[1] = lcopy(z);
! 1914: lu2 = cgetg(lg(z),t_VEC);
! 1915: for (i=1; i<lg(z); i++)
! 1916: {
! 1917: p1=(GEN)z[i]; p2=(GEN)lu[i]; lp1=lg(p1);
! 1918: p3=cgetg(lp1,t_VEC); lu2[i]=(long)p3;
! 1919: for (j=1; j<lp1; j++) p3[j] = lmul(gmael(p1,j,5),(GEN)p2[j]);
! 1920: }
! 1921: y[2]=(long)lu2; return gerepileupto(av0, y);
! 1922: }
! 1923:
! 1924: GEN
! 1925: ideallist0(GEN bnf,long bound, long flag)
! 1926: {
! 1927: if (flag<0 || flag>4) err(flagerr,"ideallist");
! 1928: return ideallistzstarall(bnf,bound,flag);
! 1929: }
! 1930:
! 1931: GEN
! 1932: ideallistzstar(GEN nf,long bound)
! 1933: {
! 1934: return ideallistzstarall(nf,bound,0);
! 1935: }
! 1936:
! 1937: GEN
! 1938: ideallistzstargen(GEN nf,long bound)
! 1939: {
! 1940: return ideallistzstarall(nf,bound,1);
! 1941: }
! 1942:
! 1943: GEN
! 1944: ideallistunit(GEN nf,long bound)
! 1945: {
! 1946: return ideallistzstarall(nf,bound,2);
! 1947: }
! 1948:
! 1949: GEN
! 1950: ideallistunitgen(GEN nf,long bound)
! 1951: {
! 1952: return ideallistzstarall(nf,bound,3);
! 1953: }
! 1954:
! 1955: GEN
! 1956: ideallist(GEN bnf,long bound)
! 1957: {
! 1958: return ideallistzstarall(bnf,bound,4);
! 1959: }
! 1960:
! 1961: static GEN
! 1962: ideallist_arch(GEN nf,GEN list,GEN arch,long flun)
! 1963: {
! 1964: long nba,i,j,lx,ly;
! 1965: GEN p1,z,p2;
! 1966:
! 1967: nba=0; for (i=1; i<lg(arch); i++) if (signe(arch[i])) nba++;
! 1968: lx=lg(list); z=cgetg(lx,t_VEC);
! 1969: for (i=1; i<lx; i++)
! 1970: {
! 1971: p2=(GEN)list[i]; checkbid(p2); ly=lg(p2);
! 1972: p1=cgetg(ly,t_VEC); z[i]=(long)p1;
! 1973: for (j=1; j<ly; j++)
! 1974: p1[j]=(long)zidealstarinitjoinarchall(nf,(GEN)p2[j],arch,nba,flun);
! 1975: }
! 1976: return z;
! 1977: }
! 1978:
! 1979: static GEN
! 1980: ideallistarchall(GEN bnf,GEN list,GEN arch,long flag)
! 1981: {
! 1982: long av,tetpil,i,j,lp1;
! 1983: long do_units = flag & 2;
! 1984: GEN nf = checknf(bnf), p1,p2,p3,racunit,funits,lu2,lu,embunit,z,y;
! 1985:
! 1986: if (typ(list) != t_VEC || (do_units && lg(list) != 3))
! 1987: err(typeer, "ideallistarch");
! 1988: if (lg(list) == 1) return cgetg(1,t_VEC);
! 1989: if (do_units)
! 1990: {
! 1991: y = cgetg(3,t_VEC);
! 1992: z = (GEN)list[1];
! 1993: lu= (GEN)list[2]; if (typ(lu) != t_VEC) err(typeer, "ideallistarch");
! 1994: }
! 1995: else z = list;
! 1996: if (typ(z) != t_VEC) err(typeer, "ideallistarch");
! 1997: z = ideallist_arch(nf,z,arch, flag & 1);
! 1998: if (!do_units) return z;
! 1999: y[1]=(long)z; av=avma;
! 2000: init_units(bnf,&funits,&racunit);
! 2001: lu2=cgetg(lg(z),t_VEC);
! 2002: for (i=1; i<lg(z); i++)
! 2003: {
! 2004: p1=(GEN)z[i]; p2=(GEN)lu[i]; lp1=lg(p1);
! 2005: p3=cgetg(lp1,t_VEC); lu2[i]=(long)p3;
! 2006: for (j=1; j<lp1; j++)
! 2007: {
! 2008: embunit = logunitmatrixarch(nf,funits,racunit,(GEN)p1[j]);
! 2009: p3[j] = lmul(gmael(p1,j,5), vconcat((GEN)p2[j],embunit));
! 2010: }
! 2011: }
! 2012: tetpil=avma; y[2]=lpile(av,tetpil,gcopy(lu2)); return y;
! 2013: }
! 2014:
! 2015: GEN
! 2016: ideallistarch(GEN nf, GEN list, GEN arch)
! 2017: {
! 2018: return ideallistarchall(nf,list,arch,0);
! 2019: }
! 2020:
! 2021: GEN
! 2022: ideallistarchgen(GEN nf, GEN list, GEN arch)
! 2023: {
! 2024: return ideallistarchall(nf,list,arch,1);
! 2025: }
! 2026:
! 2027: GEN
! 2028: ideallistunitarch(GEN bnf,GEN list,GEN arch)
! 2029: {
! 2030: return ideallistarchall(bnf,list,arch,2);
! 2031: }
! 2032:
! 2033: GEN
! 2034: ideallistunitarchgen(GEN bnf,GEN list,GEN arch)
! 2035: {
! 2036: return ideallistarchall(bnf,list,arch,3);
! 2037: }
! 2038:
! 2039: GEN
! 2040: ideallistarch0(GEN nf, GEN list, GEN arch,long flag)
! 2041: {
! 2042: if (!arch) arch=cgetg(1,t_VEC);
! 2043: if (flag<0 || flag>3) err(flagerr,"ideallistarch");
! 2044: return ideallistarchall(nf,list,arch,flag);
! 2045: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to base3.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