Annotation of OpenXM_contrib/pari-2.2/src/modules/thue.c, Revision 1.2
1.2 ! noro 1: /* $Id: thue.c,v 1.19 2002/07/17 20:23:16 karim Exp $
1.1 noro 2:
3: Copyright (C) 2000 The PARI group.
4:
5: This file is part of the PARI/GP package.
6:
7: PARI/GP is free software; you can redistribute it and/or modify it under the
8: terms of the GNU General Public License as published by the Free Software
9: Foundation. It is distributed in the hope that it will be useful, but WITHOUT
10: ANY WARRANTY WHATSOEVER.
11:
12: Check the License for details. You should have received a copy of it, along
13: with the package; see the file 'COPYING'. If not, write to the Free Software
14: Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15:
1.2 ! noro 16: #include "pari.h"
! 17:
1.1 noro 18: /* Thue equation solver. In all the forthcoming remarks, "paper"
19: * designs the paper "Thue Equations of High Degree", by Yu. Bilu and
20: * G. Hanrot, J. Number Theory (1996). Note that numbering of the constants
21: * is that of Hanrot's thesis rather than that of the paper.
1.2 ! noro 22: * The last part of the program (bnfisintnorm) was written by K. Belabas. */
1.1 noro 23:
24: /* Check whether tnf is a valid structure */
25: static int
26: checktnf(GEN tnf)
27: {
1.2 ! noro 28: int n, R, S, T;
1.1 noro 29: if (typ(tnf)!=t_VEC || (lg(tnf)!=8 && lg(tnf)!=3)) return 0;
1.2 ! noro 30: if (typ(tnf[1]) != t_POL) return 0;
! 31: if (lg(tnf) != 8) return 1; /* S=0 */
1.1 noro 32:
1.2 ! noro 33: n = degpol(tnf[1]);
! 34: if (n <= 2) err(talker,"invalid polynomial in thue (need n>2)");
! 35: S = sturm((GEN)tnf[1]); T = (n-S)>>1; R = S+T-1;
1.1 noro 36: (void)checkbnf((GEN)tnf[2]);
1.2 ! noro 37: if (typ(tnf[3]) != t_COL || lg(tnf[3]) != n+1) return 0;
! 38: if (typ(tnf[4]) != t_COL || lg(tnf[4]) != R+1) return 0;
! 39: if (typ(tnf[5]) != t_MAT || lg(tnf[5])!=R+1 || lg(gmael(tnf,5,1)) != n+1) return 0;
! 40: if (typ(tnf[6]) != t_MAT || lg(tnf[6])!=R+1 || lg(gmael(tnf,6,1)) != R+1) return 0;
1.1 noro 41: if (typ(tnf[7]) != t_VEC || lg(tnf[7]) != 7) return 0;
42: return 1;
43: }
44:
1.2 ! noro 45: static GEN
! 46: distoZ(GEN z)
1.1 noro 47: {
1.2 ! noro 48: GEN p1 = gfrac(z);
1.1 noro 49: return gmin(p1, gsub(gun,p1));
50: }
51:
1.2 ! noro 52: /* Compensates rounding errors for computation/display of the constants.
! 53: * Round up if dir > 0, down otherwise */
1.1 noro 54: static GEN
1.2 ! noro 55: myround(GEN x, long dir)
1.1 noro 56: {
1.2 ! noro 57: GEN eps = gpowgs(stoi(dir > 0? 10: -10), -10);
! 58: return gmul(x, gadd(gun, eps));
1.1 noro 59: }
60:
61: /* Returns the index of the largest element in a vector */
1.2 ! noro 62: static GEN
! 63: _Vecmax(GEN Vec, int *ind)
1.1 noro 64: {
1.2 ! noro 65: int k, tno = 1, l = lg(Vec);
1.1 noro 66: GEN tmax = (GEN)Vec[1];
1.2 ! noro 67: for (k = 2; k < l; k++)
! 68: if (gcmp((GEN)Vec[k],tmax) > 0) { tmax = (GEN)Vec[k]; tno = k; }
! 69: if (ind) *ind = tno;
! 70: return tmax;
1.1 noro 71: }
72:
1.2 ! noro 73: static GEN
! 74: Vecmax(GEN v) { return _Vecmax(v, NULL); }
1.1 noro 75:
1.2 ! noro 76: static int
! 77: Vecmaxind(GEN v) { int i; (void)_Vecmax(v, &i); return i; }
1.1 noro 78:
1.2 ! noro 79: static GEN
! 80: tnf_get_roots(GEN poly, long prec, int S, int T)
! 81: {
! 82: GEN R0 = roots(poly, prec), R = cgetg(lg(R0), t_COL);
! 83: int k;
1.1 noro 84:
1.2 ! noro 85: for (k=1; k<=S; k++) R[k] = lreal((GEN)R0[k]);
! 86: /* swap roots to get the usual order */
! 87: for (k=1; k<=T; k++)
1.1 noro 88: {
1.2 ! noro 89: R[k+S] = R0[2*k+S-1];
! 90: R[k+S+T]= R0[2*k+S];
1.1 noro 91: }
1.2 ! noro 92: return R;
! 93: }
1.1 noro 94:
1.2 ! noro 95: /* Computation of the logarithmic height of x (given by embeddings) */
! 96: static GEN
! 97: LogHeight(GEN x, long prec)
! 98: {
! 99: int i, n = lg(x)-1;
! 100: GEN LH = gun;
! 101: for (i=1; i<=n; i++) LH = gmul(LH, gmax(gun, gabs((GEN)x[i], prec)));
! 102: return gdivgs(glog(LH,prec), n);
! 103: }
1.1 noro 104:
1.2 ! noro 105: /* x^(1/n) */
! 106: static GEN
! 107: mpsqrtn(GEN x, long n)
! 108: {
! 109: if (typ(x) != t_REAL) err(typeer,"mpsqrtn");
! 110: return mpexp(divrs(mplog(x), n));
! 111: }
1.1 noro 112:
1.2 ! noro 113: /* |x|^(1/n), x t_INT */
! 114: static GEN
! 115: absisqrtn(GEN x, long n, long prec) {
! 116: GEN r = itor(x,prec);
! 117: if (signe(r) < 0) r = negr(r);
! 118: return mpsqrtn(r, n);
1.1 noro 119: }
120:
1.2 ! noro 121: static GEN
! 122: get_emb(GEN x, GEN r)
1.1 noro 123: {
1.2 ! noro 124: int l = lg(r), i, tx;
! 125: GEN e, y = cgetg(l, t_COL);
1.1 noro 126:
1.2 ! noro 127: tx = typ(x);
! 128: if (tx != t_INT && tx != t_POL) err(typeer,"get_emb");
! 129: for (i=1; i<l; i++)
! 130: {
! 131: e = poleval(x, (GEN)r[i]);
! 132: if (gcmp0(e)) return NULL;
! 133: y[i] = (long)e;
! 134: }
! 135: return y;
1.1 noro 136: }
137:
1.2 ! noro 138: /* Computation of the conjugates (sigma_i(v_j)), and log. heights, of elts of v */
! 139: static GEN
! 140: Conj_LH(GEN v, GEN *H, GEN r, long prec)
1.1 noro 141: {
1.2 ! noro 142: int j, l = lg(v);
! 143: GEN e, M = (GEN)cgetg(l,t_MAT);
1.1 noro 144:
1.2 ! noro 145: (*H) = cgetg(l,t_COL);
! 146: for (j = 1; j < l; j++)
1.1 noro 147: {
1.2 ! noro 148: if (! (e = get_emb((GEN)v[j], r)) ) return NULL; /* FAIL */
! 149: M[j] = (long)e;
! 150: (*H)[j]= (long)LogHeight(e, prec);
1.1 noro 151: }
1.2 ! noro 152: return M;
1.1 noro 153: }
154:
1.2 ! noro 155: static GEN abslog(GEN x, long prec) { return gabs(glog(x,prec), prec); }
! 156: static GEN logabs(GEN x, long prec) { return glog(gabs(x,prec), prec); }
! 157:
! 158: /* Computation of M, its inverse A and precision check (see paper) */
! 159: static GEN
! 160: T_A_Matrices(GEN MatFU, int r, GEN *eps5, long prec)
1.1 noro 161: {
1.2 ! noro 162: GEN A, p1, m1, IntM, nia, eps3, eps2, eps1 = shifti(gun, bit_accuracy(prec));
1.1 noro 163:
1.2 ! noro 164: m1 = rowextract_i(vecextract_i(MatFU, 1,r), 1,r); /* minor order r */
! 165: m1 = logabs(m1,prec);
1.1 noro 166:
1.2 ! noro 167: A = invmat(m1);
! 168: IntM = gsub(gmul(A,m1), idmat(r));
1.1 noro 169:
1.2 ! noro 170: eps2 = gadd(vecmax(gabs(IntM,prec)), ginv(eps1));
! 171: nia = vecmax(gabs(A,prec));
1.1 noro 172:
1.2 ! noro 173: /* Check for the precision in matrix inversion. See paper, Lemma 2.4.2. */
! 174: p1 = gadd(gmulsg(r, gmul(nia,eps1)), eps2);
! 175: if (gcmp(gmulgs(p1, 2*r), gun) < 0) err(precer,"thue");
1.1 noro 176:
1.2 ! noro 177: p1 = gadd(gmulsg(r, gdiv(nia,eps1)), eps2);
! 178: eps3 = gmul(gmulsg(2*r*r,nia), p1);
! 179: eps3 = myround(eps3, 1);
1.1 noro 180:
1.2 ! noro 181: if (DEBUGLEVEL>1) fprintferr("epsilon_3 -> %Z\n",eps3);
! 182: *eps5 = mulsr(r, eps3); return A;
1.1 noro 183: }
184:
1.2 ! noro 185: /* Performs basic computations concerning the equation.
! 186: * Returns a "tnf" structure containing
! 187: * 1) the polynomial
! 188: * 2) the bnf (used to solve the norm equation)
! 189: * 3) roots, with presumably enough precision
! 190: * 4) The logarithmic heights of units
! 191: * 5) The matrix of conjugates of units
! 192: * 6) its inverse
! 193: * 7) a few technical constants */
1.1 noro 194: static GEN
1.2 ! noro 195: inithue(GEN P, GEN bnf, long flag, long prec)
1.1 noro 196: {
1.2 ! noro 197: GEN MatFU, x0, tnf, tmp, gpmin, dP, csts, ALH, eps5, ro, c1, c2;
! 198: int k,j, s,t, n = degpol(P);
1.1 noro 199:
1.2 ! noro 200: if (!bnf)
! 201: {
! 202: bnf = bnfinit0(P,1,NULL,prec);
! 203: if (bnf != checkbnf_discard(bnf)) err(talker,"non-monic polynomial in thue");
! 204: if (flag) (void)certifybuchall(bnf);
! 205: }
! 206: nf_get_sign(checknf(bnf), (long*)&s, (long*)&t);
! 207: ro = tnf_get_roots(P, prec, s, t);
! 208: MatFU = Conj_LH(gmael(bnf,8,5), &ALH, ro, prec);
! 209: if (!MatFU) return NULL; /* FAIL */
1.1 noro 210:
1.2 ! noro 211: dP = derivpol(P);
! 212: c1 = NULL; /* min |P'(r_i)|, i <= s */
! 213: for (k=1; k<=s; k++)
1.1 noro 214: {
1.2 ! noro 215: tmp = gabs(poleval(dP,(GEN)ro[k]),prec);
! 216: if (!c1 || gcmp(tmp,c1) < 0) c1 = tmp;
1.1 noro 217: }
1.2 ! noro 218: c1 = gdiv(shifti(gun,n-1), c1);
! 219: c1 = gprec_w(myround(c1, 1), DEFAULTPREC);
1.1 noro 220:
1.2 ! noro 221: c2 = NULL; /* max |r_i - r_j|, i!=j */
! 222: for (k=1; k<=n; k++)
! 223: for (j=k+1; j<=n; j++)
! 224: {
! 225: tmp = gabs(gsub((GEN)ro[j],(GEN)ro[k]), prec);
! 226: if (!c2 || gcmp(c2,tmp) > 0) c2 = tmp;
! 227: }
! 228: c2 = gprec_w(myround(c2, -1), DEFAULTPREC);
1.1 noro 229:
1.2 ! noro 230: if (t==0) x0 = gun;
! 231: else
1.1 noro 232: {
1.2 ! noro 233: gpmin = NULL; /* min |P'(r_i)|, i > s */
! 234: for (k=1; k<=t; k++)
1.1 noro 235: {
1.2 ! noro 236: tmp = gabs(poleval(dP,(GEN)ro[s+k]), prec);
! 237: if (!gpmin || gcmp(tmp,gpmin) < 0) gpmin = tmp;
1.1 noro 238: }
1.2 ! noro 239: gpmin = gprec_w(gpmin, DEFAULTPREC);
! 240:
! 241: /* Compute x0. See paper, Prop. 2.2.1 */
! 242: x0 = gmul(gpmin, Vecmax(gabs(gimag(ro), prec)));
! 243: x0 = mpsqrtn(gdiv(shifti(gun,n-1), x0), n);
1.1 noro 244: }
1.2 ! noro 245: if (DEBUGLEVEL>1) fprintferr("c1 = %Z\nc2 = %Z\n", c1, c2);
! 246:
! 247: ALH = gmul2n(ALH, 1);
! 248: tnf = cgetg(8,t_VEC); csts = cgetg(7,t_VEC);
! 249: tnf[1] = (long)P;
! 250: tnf[2] = (long)bnf;
! 251: tnf[3] = (long)ro;
! 252: tnf[4] = (long)ALH;
! 253: tnf[5] = (long)MatFU;
! 254: tnf[6] = (long)T_A_Matrices(MatFU, s+t-1, &eps5, prec);
! 255: tnf[7] = (long)csts;
! 256: csts[1] = (long)c1; csts[2] = (long)c2; csts[3] = (long)LogHeight(ro, prec);
! 257: csts[4] = (long)x0; csts[5] = (long)eps5; csts[6] = (long)stoi(prec);
! 258: return tnf;
! 259: }
! 260:
! 261: typedef struct {
! 262: GEN c10, c11, c13, c15, bak, NE, ALH, hal, MatFU, ro, Hmu;
! 263: GEN delta, lambda, errdelta;
! 264: int r, iroot;
! 265: } baker_s;
1.1 noro 266:
1.2 ! noro 267: /* Compute Baker's bound c9 and B_0, the bound for the b_i's. See Thm 2.3.1 */
! 268: static GEN
! 269: Baker(baker_s *BS)
1.1 noro 270: {
1.2 ! noro 271: const long prec = DEFAULTPREC;
! 272: GEN tmp, B0, hb0, c9 = gun, ro = BS->ro, ro0 = (GEN)ro[BS->iroot];
! 273: int k, i1, i2, r = BS->r;
! 274:
! 275: switch (BS->iroot) {
! 276: case 1: i1=2; i2=3; break;
! 277: case 2: i1=1; i2=3; break;
! 278: default: i1=1; i2=2; break;
1.1 noro 279: }
280:
281: /* Compute h_1....h_r */
282: for (k=1; k<=r; k++)
283: {
1.2 ! noro 284: tmp = gdiv(gcoeff(BS->MatFU,i1,k), gcoeff(BS->MatFU,i2,k));
! 285: tmp = gmax(gun, abslog(tmp,prec));
! 286: c9 = gmul(c9, gmax((GEN)BS->ALH[k], gdiv(tmp, BS->bak)));
1.1 noro 287: }
288:
289: /* Compute a bound for the h_0 */
1.2 ! noro 290: hb0 = gadd(gmul2n(BS->hal,2), gmul(gdeux, gadd(BS->Hmu,mplog2(prec))));
! 291: tmp = gdiv(gmul(gsub(ro0, (GEN)ro[i2]), (GEN)BS->NE[i1]),
! 292: gmul(gsub(ro0, (GEN)ro[i1]), (GEN)BS->NE[i2]));
! 293: tmp = gmax(gun, abslog(tmp, prec));
! 294: hb0 = gmax(hb0, gdiv(tmp, BS->bak));
! 295: c9 = gmul(c9,hb0);
1.1 noro 296: /* Multiply c9 by the "constant" factor */
1.2 ! noro 297: c9 = gmul(c9, gmul(mulri(mulsr(18,mppi(prec)), gpowgs(stoi(32),4+r)),
! 298: gmul(gmul(mpfact(r+3), gpowgs(mulis(BS->bak,r+2), r+3)),
! 299: glog(mulis(BS->bak,2*(r+2)),prec))));
! 300: c9 = gprec_w(myround(c9, 1), DEFAULTPREC);
1.1 noro 301: /* Compute B0 according to Lemma 2.3.3 */
1.2 ! noro 302: B0 = mulsr(2, divrr(addrr(mulrr(c9,mplog(divrr(c9,BS->c10))), mplog(BS->c11)),
! 303: BS->c10));
! 304: B0 = gmax(B0, dbltor(2.71828183));
1.1 noro 305:
1.2 ! noro 306: if (DEBUGLEVEL>1) {
! 307: fprintferr(" B0 = %Z\n",B0);
! 308: fprintferr(" Baker = %Z\n",c9);
! 309: }
! 310: return B0;
1.1 noro 311: }
312:
1.2 ! noro 313: /* || x d ||, x t_REAL, d t_INT */
! 314: static GEN
! 315: errnum(GEN x, GEN d)
1.1 noro 316: {
1.2 ! noro 317: GEN dx = mulir(d, x);
! 318: return mpabs(subri(dx, ground(dx)));
1.1 noro 319: }
320:
321: /* Try to reduce the bound through continued fractions; see paper. */
322: static int
1.2 ! noro 323: CF_1stPass(GEN *B0, GEN kappa, baker_s *BS)
1.1 noro 324: {
1.2 ! noro 325: GEN q, ql, qd, l0, denbound = mulri(*B0, kappa);
! 326:
! 327: if (gcmp(gmul(dbltor(0.1),gsqr(denbound)), ginv(BS->errdelta)) > 0) return -1;
1.1 noro 328:
1.2 ! noro 329: q = denom( bestappr(BS->delta, denbound) );
! 330: qd = errnum(BS->delta, q);
! 331: ql = errnum(BS->lambda,q);
! 332:
! 333: l0 = subrr(ql, addrr(mulrr(qd, *B0), divri(dbltor(0.1),kappa)));
! 334: if (signe(l0) <= 0) return 0;
1.1 noro 335:
1.2 ! noro 336: if (BS->r > 1)
! 337: *B0 = divrr(mplog(divrr(mulir(q,BS->c15), l0)), BS->c13);
! 338: else
1.1 noro 339: {
1.2 ! noro 340: l0 = mulrr(l0,Pi2n(1, DEFAULTPREC));
! 341: *B0 = divrr(mplog(divrr(mulir(q,BS->c11), l0)), BS->c10);
1.1 noro 342: }
1.2 ! noro 343: if (DEBUGLEVEL>1) fprintferr(" B0 -> %Z\n",*B0);
! 344: return 1;
! 345: }
1.1 noro 346:
1.2 ! noro 347: /* Check whether a solution has already been found */
! 348: static int
! 349: new_sol(GEN z, GEN S)
! 350: {
! 351: int i, l = lg(S);
! 352: for (i=1; i<l; i++)
! 353: if (gegal(z,(GEN)S[i])) return 0;
! 354: return 1;
! 355: }
1.1 noro 356:
1.2 ! noro 357: static void
! 358: add_sol(GEN *pS, GEN x, GEN y)
! 359: {
! 360: GEN u = cgetg(3,t_VEC); u[1] = (long)x; u[2] = (long)y;
! 361: if (new_sol(u, *pS)) *pS = concatsp(*pS, _vec(u));
! 362: }
1.1 noro 363:
364: static void
1.2 ! noro 365: check_sol(GEN x, GEN y, GEN P, GEN rhs, GEN *pS)
1.1 noro 366: {
1.2 ! noro 367: if (gcmp0(y))
! 368: { if (gegal(gpowgs(x,degpol(P)), rhs)) add_sol(pS, x, gzero); }
! 369: else
! 370: { if (gegal(poleval(rescale_pol(P,y),x), rhs)) add_sol(pS, x, y); }
! 371: }
1.1 noro 372:
1.2 ! noro 373: /* Check whether a potential solution is a true solution. Return 0 if
! 374: * truncation error (increase precision) */
! 375: static int
! 376: CheckSol(GEN *pS, GEN z1, GEN z2, GEN P, GEN rhs, GEN ro)
! 377: {
! 378: GEN x, y, ro1 = (GEN)ro[1], ro2 = (GEN)ro[2];
! 379: long e;
1.1 noro 380:
1.2 ! noro 381: y = grndtoi(greal(gdiv(gsub(z2,z1), gsub(ro1,ro2))), &e);
! 382: if (e > 0) return 0;
! 383: x = gadd(z1, gmul(ro1, y));
! 384: x = grndtoi(greal(x), &e);
! 385: if (e > 0) return 0;
! 386: if (e <= -13)
1.1 noro 387: {
1.2 ! noro 388: check_sol( x , y, P,rhs,pS);
! 389: check_sol(gneg(x), gneg(y), P,rhs,pS);
1.1 noro 390: }
1.2 ! noro 391: return 1;
1.1 noro 392: }
393:
1.2 ! noro 394: /* find q1,q2,q3 st q1 + b q2 + c q3 ~ 0 */
1.1 noro 395: static GEN
1.2 ! noro 396: GuessQi(GEN b, GEN c, GEN *eps)
1.1 noro 397: {
1.2 ! noro 398: GEN Q, Lat, C = gpowgs(stoi(10), 10);
1.1 noro 399:
1.2 ! noro 400: Lat = idmat(3);
! 401: mael(Lat,1,3) = (long)C;
! 402: mael(Lat,2,3) = lround(gmul(C,b));
! 403: mael(Lat,3,3) = lround(gmul(C,c));
1.1 noro 404:
1.2 ! noro 405: Q = (GEN)lllint(Lat)[1];
! 406: if (gcmp0((GEN)Q[3])) return NULL; /* FAIL */
1.1 noro 407:
1.2 ! noro 408: *eps = gadd(gadd((GEN)Q[1], gmul((GEN)Q[2],b)), gmul((GEN)Q[3],c));
! 409: *eps = mpabs(*eps); return Q;
1.1 noro 410: }
411:
412: /* Check for solutions under a small bound (see paper) */
1.2 ! noro 413: static GEN
! 414: SmallSols(GEN S, int Bx, GEN poly, GEN rhs, GEN ro)
1.1 noro 415: {
1.2 ! noro 416: gpmem_t av = avma, lim = stack_lim(av, 1);
! 417: const long prec = DEFAULTPREC;
! 418: GEN Hpoly, interm, X, Xn, Xnm1, Y, sqrtnRHS;
! 419: int x, y, j, By, n = degpol(poly);
1.1 noro 420: double bndyx;
421:
1.2 ! noro 422: if (DEBUGLEVEL>1) fprintferr("* Checking for small solutions\n");
! 423:
! 424: sqrtnRHS = absisqrtn(rhs, n, prec);
! 425: bndyx = gtodouble(gadd(sqrtnRHS, Vecmax(gabs(ro,prec))));
! 426:
! 427: /* x = 0 first */
! 428: Y = ground(sqrtnRHS);
! 429: if (gegal(gpowgs(Y,n), rhs)) add_sol(&S, Y, gzero);
! 430: Y = negi(Y);
! 431: if (gegal(gpowgs(Y,n), rhs)) add_sol(&S, Y, gzero);
! 432: /* x != 0 */
! 433: for (x = -Bx; x <= Bx; x++)
! 434: {
! 435: if (!x) continue;
! 436:
! 437: X = stoi(x); Xn = gpowgs(X,n);
! 438: interm = subii(rhs, mulii(Xn, (GEN)poly[2]));
! 439: Xnm1 = Xn; j = 2;
! 440: while (gcmp0(interm))
! 441: {
! 442: Xnm1 = divis(Xnm1, x);
! 443: interm = mulii((GEN)poly[++j], Xnm1);
! 444: }
! 445: /* y | interm */
1.1 noro 446:
1.2 ! noro 447: Hpoly = rescale_pol(poly,X); /* homogenize */
! 448: By = (int)(x>0? bndyx*x: -bndyx*x);
! 449: if (gegal(gmul((GEN)poly[2],Xn),rhs)) add_sol(&S, gzero, X); /* y = 0 */
! 450: for(y = -By; y <= By; y++)
! 451: {
! 452: if (!y || smodis(interm, y)) continue;
! 453: Y = stoi(y);
! 454: if (gegal(poleval(Hpoly, Y), rhs)) add_sol(&S, Y, X);
! 455: }
1.1 noro 456:
457: if (low_stack(lim,stack_lim(av,1)))
458: {
459: if(DEBUGMEM>1) err(warnmem,"Check_small");
1.2 ! noro 460: S = gerepilecopy(av, S);
1.1 noro 461: }
462: }
1.2 ! noro 463: return S;
1.1 noro 464: }
465:
466: /* Computes [x]! */
467: static double
468: fact(double x)
469: {
1.2 ! noro 470: double ft = 1.0;
! 471: x = (int)x; while (x>1) { ft *= x; x--; }
1.1 noro 472: return ft ;
473: }
474:
475: /* From a polynomial and optionally a system of fundamental units for the
1.2 ! noro 476: * field defined by poly, computes all relevant constants needed to solve
! 477: * the equation P(x,y)=a given the solutions of N_{K/Q}(x)=a (see inithue). */
1.1 noro 478: GEN
479: thueinit(GEN poly, long flag, long prec)
480: {
1.2 ! noro 481: GEN tnf, bnf = NULL;
! 482: gpmem_t av = avma;
! 483: long k, s;
! 484:
! 485: if (checktnf(poly)) { bnf = checkbnf((GEN)poly[2]); poly = (GEN)poly[1]; }
! 486: if (typ(poly)!=t_POL) err(notpoler,"thueinit");
! 487: if (degpol(poly) <= 2) err(talker,"invalid polynomial in thue (need deg>2)");
! 488:
! 489: if (!gisirreducible(poly)) err(redpoler,"thueinit");
! 490: s = sturm(poly);
! 491: if (s)
! 492: {
! 493: long PREC, n = degpol(poly);
! 494: double d, dr, dn = (double)n;
1.1 noro 495:
1.2 ! noro 496: dr = (double)((s+n-2)>>1); /* s+t-1 */
! 497: d = dn*(dn-1)*(dn-2);
! 498: /* Guess precision by approximating Baker's bound. The guess is most of
! 499: * the time not sharp, ie 10 to 30 decimal digits above what is _really_
! 500: * necessary. Note that the limiting step is the reduction. See paper. */
! 501: PREC = 3 + (long)((5.83 + (dr+4)*5 + log(fact(dr+3)) + (dr+3)*log(dr+2) +
! 502: (dr+3)*log(d) + log(log(2*d*(dr+2))) + (dr+1)) / 10.);
! 503: if (PREC < prec) PREC = prec;
! 504: for (;;)
! 505: {
! 506: if (( tnf = inithue(poly, bnf, flag, PREC) )) break;
! 507: PREC = (PREC<<1)-2;
! 508: if (DEBUGLEVEL>1) err(warnprec,"thueinit",PREC);
! 509: bnf = NULL; avma = av;
! 510: }
! 511: }
1.1 noro 512: else
1.2 ! noro 513: {
! 514: GEN c0 = gun, ro = roots(poly, DEFAULTPREC);
! 515: for (k=1; k<lg(ro); k++) c0 = gmul(c0, gimag((GEN)ro[k]));
! 516: c0 = ginv( mpabs(c0) );
! 517: tnf = cgetg(3,t_VEC);
! 518: tnf[1] = (long)poly;
! 519: tnf[2] = (long)c0;
! 520: }
! 521: return gerepilecopy(av,tnf);
! 522: }
1.1 noro 523:
1.2 ! noro 524: static GEN
! 525: get_B0(int i1, GEN Delta, GEN Lambda, GEN eps5, long prec, baker_s *BS)
! 526: {
! 527: GEN delta, lambda, errdelta, B0 = Baker(BS);
! 528: int i2, r = BS->r;
! 529:
! 530: i2 = (i1 == 1)? 2: 1;
! 531: for(;;) /* i2 from 1 to r unless r = 1 [then i2 = 2] */
1.1 noro 532: {
1.2 ! noro 533: if (r > 1)
! 534: {
! 535: delta = divrr((GEN)Delta[i2],(GEN)Delta[i1]);
! 536: lambda = gdiv(gsub(gmul((GEN)Delta[i2],(GEN)Lambda[i1]),
! 537: gmul((GEN)Delta[i1],(GEN)Lambda[i2])),
! 538: (GEN)Delta[i1]);
! 539: errdelta = mulrr(addsr(1,delta),
! 540: divrr(eps5, subrr(mpabs((GEN)Delta[i1]),eps5)));
! 541: }
! 542: else
! 543: { /* r == 1, single fundamental unit (i1 = s = t = 1) */
! 544: GEN p1, Pi2 = Pi2n(1, prec);
! 545: GEN fu = (GEN)BS->MatFU[1], ro = BS->ro;
! 546:
! 547: p1 = gdiv((GEN)fu[2], (GEN)fu[3]);
! 548: delta = divrr(garg(p1,prec), Pi2);
! 549:
! 550: p1 = gmul(gdiv(gsub((GEN)ro[1], (GEN)ro[2]),
! 551: gsub((GEN)ro[1], (GEN)ro[3])),
! 552: gdiv((GEN)BS->NE[3], (GEN)BS->NE[2]));
! 553: lambda = divrr(garg(p1,prec), Pi2);
! 554:
! 555: errdelta = gdiv(gmul2n(gun, 1 - bit_accuracy(prec)),
! 556: gabs((GEN)fu[2],prec));
! 557: }
! 558: BS->delta = delta;
! 559: BS->lambda= lambda; BS->errdelta = errdelta;
! 560: if (DEBUGLEVEL>1) fprintferr(" errdelta = %Z\n",errdelta);
! 561:
! 562: if (DEBUGLEVEL>1) fprintferr(" Entering CF...\n");
! 563: /* Reduce B0 as long as we make progress: newB0 < oldB0 - 0.1 */
! 564: for (;;)
! 565: {
! 566: GEN oldB0 = B0, kappa = stoi(10);
! 567: int cf;
! 568:
! 569: for (cf = 0; cf < 10; cf++, kappa = mulis(kappa,10))
! 570: {
! 571: int res = CF_1stPass(&B0, kappa, BS);
! 572: if (res < 0) return NULL; /* prec problem */
! 573: if (res) break;
! 574: if (DEBUGLEVEL>1) fprintferr("CF failed. Increasing kappa\n");
! 575: }
! 576: if (cf == 10)
! 577: { /* Semirational or totally rational case */
! 578: GEN Q, ep, q, l0, denbound;
! 579:
! 580: if (! (Q = GuessQi(delta, lambda, &ep)) ) break;
! 581:
! 582: denbound = gadd(B0, absi((GEN)Q[2]));
! 583: q = denom( bestappr(delta, denbound) );
! 584: l0 = subrr(errnum(delta, q), ep);
! 585: if (signe(l0) <= 0) break;
! 586:
! 587: B0 = divrr(mplog(divrr(mulir((GEN)Q[3], BS->c15), l0)), BS->c13);
! 588: if (DEBUGLEVEL>1) fprintferr("Semirat. reduction: B0 -> %Z\n",B0);
! 589: }
! 590: /* if no progress, stop */
! 591: if (gcmp(oldB0, gadd(B0,dbltor(0.1))) <= 0) return gmin(oldB0, B0);
! 592: }
! 593: i2++; if (i2 == i1) i2++;
! 594: if (i2 > r) break;
! 595: }
! 596: err(bugparier,"thue (totally rational case)");
! 597: return NULL; /* not reached */
1.1 noro 598: }
599:
1.2 ! noro 600: static GEN
! 601: LargeSols(GEN tnf, GEN rhs, GEN ne, GEN *pro, GEN *pS)
1.1 noro 602: {
1.2 ! noro 603: GEN Vect, P, ro, bnf, MatFU, A, csts, dP, vecdP;
! 604: GEN c1,c2,c3,c4,c10,c11,c13,c14,c15, x0, x1, x2, x3, b, zp1, tmp, eps5;
! 605: int iroot, ine, n, i, r,s,t;
! 606: long upb, bi1, Prec, prec;
! 607: baker_s BS;
! 608: gpmem_t av = avma;
! 609:
! 610: bnf = (GEN)tnf[2];
! 611: if (!ne) ne = bnfisintnorm(bnf, rhs);
! 612: if (lg(ne)==1) return NULL;
! 613: nf_get_sign(checknf(bnf), (long*)&s, (long*)&t);
! 614: BS.r = r = s+t-1;
! 615: P = (GEN)tnf[1]; n = degpol(P);
! 616: ro = (GEN)tnf[3];
! 617: BS.ALH = (GEN)tnf[4];
! 618: MatFU = (GEN)tnf[5];
! 619: A = (GEN)tnf[6];
! 620: csts = (GEN)tnf[7];
! 621: c1 = (GEN)csts[1]; c1 = gmul(absi(rhs), c1);
! 622: c2 = (GEN)csts[2];
! 623: BS.hal = (GEN)csts[3];
! 624: x0 = (GEN)csts[4];
! 625: eps5 = (GEN)csts[5];
! 626: Prec = gtolong((GEN)csts[6]);
! 627: BS.MatFU = MatFU;
! 628: BS.bak = mulss(n, (n-1)*(n-2)); /* safe */
! 629: *pS = cgetg(1, t_VEC);
! 630:
! 631: if (t) x0 = gmul(x0, absisqrtn(rhs, n, Prec));
! 632: tmp = divrr(c1,c2);
! 633: c3 = mulrr(dbltor(1.39), tmp);
! 634: c4 = mulsr(n-1, c3);
! 635: x1 = gmax(x0, mpsqrtn(mulsr(2,tmp),n));
! 636:
! 637: Vect = cgetg(r+1,t_COL); for (i=1; i<=r; i++) Vect[i]=un;
! 638: Vect = gmul(gabs(A,DEFAULTPREC), Vect);
! 639: c14 = mulrr(c4, Vecmax(Vect));
! 640: x2 = gmax(x1, mpsqrtn(mulsr(10,c14), n));
! 641: if (DEBUGLEVEL>1) {
! 642: fprintferr("x1 -> %Z\n",x1);
! 643: fprintferr("x2 -> %Z\n",x2);
! 644: fprintferr("c14 = %Z\n",c14);
! 645: }
1.1 noro 646:
1.2 ! noro 647: dP = derivpol(P);
! 648: vecdP = cgetg(s+1, t_VEC);
! 649: for (i=1; i<=s; i++) vecdP[i] = (long)poleval(dP, (GEN)ro[i]);
! 650:
! 651: zp1 = dbltor(0.01);
! 652: x3 = gmax(x2, mpsqrtn(mulsr(2,divrr(c14,zp1)),n));
! 653:
! 654: b = cgetg(r+1,t_COL);
! 655: for (iroot=1; iroot<=s; iroot++)
! 656: {
! 657: GEN Delta, MatNE, Hmu, c5, c7;
! 658:
! 659: Vect = cgetg(r+1,t_COL); for (i=1; i<=r; i++) Vect[i] = un;
! 660: if (iroot <= r) Vect[iroot] = lstoi(1-n);
! 661: Delta = gmul(A,Vect);
! 662:
! 663: c5 = Vecmax(gabs(Delta,Prec));
! 664: c5 = myround(gprec_w(c5,DEFAULTPREC), 1);
! 665: c7 = mulsr(r,c5);
! 666: c10 = divsr(n,c7); BS.c10 = c10;
! 667: c13 = divsr(n,c5); BS.c13 = c13;
! 668: if (DEBUGLEVEL>1) {
! 669: fprintferr("* real root no %ld/%ld\n", iroot,s);
! 670: fprintferr(" c10 = %Z\n",c10);
! 671: fprintferr(" c13 = %Z\n",c13);
! 672: }
! 673:
! 674: prec = Prec;
! 675: for (;;)
! 676: {
! 677: if (( MatNE = Conj_LH(ne, &Hmu, ro, prec) )) break;
! 678: prec = (prec<<1)-2;
! 679: if (DEBUGLEVEL>1) err(warnprec,"thue",prec);
! 680: ro = tnf_get_roots(P, prec, s, t);
! 681: }
! 682: BS.ro = ro;
! 683: BS.iroot = iroot;
1.1 noro 684:
1.2 ! noro 685: for (ine=1; ine<lg(ne); ine++)
1.1 noro 686: {
1.2 ! noro 687: GEN Lambda, B0, c6, c8;
! 688: GEN NE = (GEN)MatNE[ine], Vect2 = cgetg(r+1,t_COL);
! 689: int k, i1;
! 690:
! 691: if (DEBUGLEVEL>1) fprintferr(" - norm sol. no %ld/%ld\n",ine,lg(ne)-1);
! 692: for (k=1; k<=r; k++)
1.1 noro 693: {
1.2 ! noro 694: if (k == iroot)
! 695: tmp = gdiv(rhs, gmul((GEN)vecdP[k], (GEN)NE[k]));
! 696: else
! 697: tmp = gdiv(gsub((GEN)ro[iroot],(GEN)ro[k]), (GEN)NE[k]);
! 698: Vect2[k] = llog(gabs(tmp,prec), prec);
! 699: }
! 700: Lambda = gmul(A,Vect2);
! 701:
! 702: c6 = addrr(dbltor(0.1), Vecmax(gabs(Lambda,DEFAULTPREC)));
! 703: c6 = myround(c6, 1);
! 704: c8 = addrr(dbltor(1.23), mulsr(r,c6));
! 705: c11= mulrr(mulsr(2,c3) , mpexp(divrr(mulsr(n,c8),c7)));
! 706: c15= mulrr(mulsr(2,c14), mpexp(divrr(mulsr(n,c6),c5)));
! 707:
! 708: if (DEBUGLEVEL>1) {
! 709: fprintferr(" c6 = %Z\n",c6);
! 710: fprintferr(" c8 = %Z\n",c8);
! 711: fprintferr(" c11 = %Z\n",c11);
! 712: fprintferr(" c15 = %Z\n",c15);
! 713: }
! 714: BS.c11 = c11;
! 715: BS.c15 = c15;
! 716: BS.NE = NE;
! 717: BS.Hmu = (GEN)Hmu[ine];
! 718:
! 719: i1 = Vecmaxind(gabs(Delta,prec));
! 720: if (! (B0 = get_B0(i1, Delta, Lambda, eps5, prec, &BS)) ) goto PRECPB;
! 721:
! 722: /* For each possible value of b_i1, compute the b_i's
! 723: * and 2 conjugates of z = x - alpha y. Then check. */
! 724: upb = gtolong(gceil(B0));
! 725: for (bi1=-upb; bi1<=upb; bi1++)
! 726: {
! 727: GEN z1, z2;
! 728: for (i=1; i<=r; i++)
! 729: {
! 730: b[i] = ldiv(gsub(gmul((GEN)Delta[i], stoi(bi1)),
! 731: gsub(gmul((GEN)Delta[i],(GEN)Lambda[i1]),
! 732: gmul((GEN)Delta[i1],(GEN)Lambda[i]))),
! 733: (GEN)Delta[i1]);
! 734: if (gcmp(distoZ((GEN)b[i]), zp1) > 0) break;
! 735: }
! 736: if (i <= r) continue;
! 737:
! 738: z1 = z2 = gun;
! 739: for(i=1; i<=r; i++)
! 740: {
! 741: GEN c = ground((GEN)b[i]);
! 742: z1 = gmul(z1, powgi(gcoeff(MatFU,1,i), c));
! 743: z2 = gmul(z2, powgi(gcoeff(MatFU,2,i), c));
! 744: }
! 745: z1 = gmul(z1, (GEN)NE[1]);
! 746: z2 = gmul(z2, (GEN)NE[2]);
! 747: if (!CheckSol(pS, z1,z2,P,rhs,ro)) goto PRECPB;
1.1 noro 748: }
749: }
750: }
1.2 ! noro 751: *pro = ro; return x3;
1.1 noro 752:
1.2 ! noro 753: PRECPB:
! 754: ne = gerepilecopy(av, ne);
! 755: prec += 5 * (DEFAULTPREC-2);
! 756: if (DEBUGLEVEL>1) err(warnprec,"thue",prec);
! 757: tnf = inithue(P, bnf, 0, prec);
! 758: return LargeSols(tnf, rhs, ne, pro, pS);
! 759: }
! 760:
! 761: /* Given a tnf structure as returned by thueinit, a RHS and
! 762: * optionally the solutions to the norm equation, returns the solutions to
! 763: * the Thue equation F(x,y)=a
! 764: */
! 765: GEN
! 766: thue(GEN tnf, GEN rhs, GEN ne)
! 767: {
! 768: gpmem_t av = avma;
! 769: GEN P, ro, x3, S;
! 770:
! 771: if (!checktnf(tnf)) err(talker,"not a tnf in thue");
! 772: if (typ(rhs) != t_INT) err(typeer,"thue");
! 773:
! 774: P = (GEN)tnf[1];
! 775: if (lg(tnf) == 8)
! 776: {
! 777: x3 = LargeSols(tnf, rhs, ne, &ro, &S);
! 778: if (!x3) { avma = av; return cgetg(1,t_VEC); }
! 779: }
! 780: else
! 781: { /* Case s=0. All solutions are "small". */
! 782: GEN c0 = (GEN)tnf[2]; /* t_REAL */
! 783: S = cgetg(1,t_VEC);
! 784: ro = roots(P, DEFAULTPREC);
! 785: x3 = mpsqrtn(mulir(constant_term(P), divir(absi(rhs),c0)), degpol(P));
! 786: }
! 787: return gerepilecopy(av, SmallSols(S, gtolong(x3), P, rhs, ro));
1.1 noro 788: }
789:
790: static GEN *Relations; /* primes above a, expressed on generators of Cl(K) */
791: static GEN *Partial; /* list of vvectors, Partial[i] = Relations[1..i] * u[1..i] */
792: static GEN *gen_ord; /* orders of generators of Cl(K) given in bnf */
793:
794: static long *f; /* f[i] = f(Primes[i]/p), inertial degree */
795: static long *n; /* a = prod p^{ n_p }. n[i]=n_p if Primes[i] divides p */
796: static long *inext; /* index of first P above next p, 0 if p is last */
797: static long *S; /* S[i] = n[i] - sum_{ 1<=k<=i } f[k].u[k] */
798: static long *u; /* We want principal ideals I = prod Primes[i]^u[i] */
799: static long **normsol; /* lists of copies of the u[] which are solutions */
800:
801: static long Nprimes; /* length(Relations) = #{max ideal above divisors of a} */
802: static long sindex, smax; /* current index in normsol; max. index */
803:
804: /* u[1..i] has been filled. Norm(u) is correct.
805: * Check relations in class group then save it.
806: */
807: static void
808: test_sol(long i)
809: {
810: long k,*sol;
811:
812: if (Partial)
813: {
1.2 ! noro 814: gpmem_t av=avma;
1.1 noro 815: for (k=1; k<lg(Partial[1]); k++)
816: if ( signe(modii( (GEN)Partial[i][k], gen_ord[k] )) )
817: { avma=av; return; }
818: avma=av;
819: }
820: if (sindex == smax)
821: {
822: long new_smax = smax << 1;
823: long **new_normsol = (long **) new_chunk(new_smax+1);
824:
825: for (k=1; k<=smax; k++) new_normsol[k] = normsol[k];
826: normsol = new_normsol; smax = new_smax;
827: }
828: sol = cgetg(Nprimes+1,t_VECSMALL);
829: normsol[++sindex] = sol;
830:
831: for (k=1; k<=i; k++) sol[k]=u[k];
832: for ( ; k<=Nprimes; k++) sol[k]=0;
1.2 ! noro 833: if (DEBUGLEVEL>2)
1.1 noro 834: {
835: fprintferr("sol = %Z\n",sol);
836: if (Partial) fprintferr("Partial = %Z\n",Partial);
837: flusherr();
838: }
839: }
840: static void
841: fix_Partial(long i)
842: {
1.2 ! noro 843: long k;
! 844: gpmem_t av = avma;
1.1 noro 845: for (k=1; k<lg(Partial[1]); k++)
846: addiiz(
847: (GEN) Partial[i-1][k],
848: mulis((GEN) Relations[i][k], u[i]),
849: (GEN) Partial[i][k]
850: );
851: avma=av;
852: }
853:
854: /* Recursive loop. Suppose u[1..i] has been filled
855: * Find possible solutions u such that, Norm(prod Prime[i]^u[i]) = a, taking
856: * into account:
857: * 1) the relations in the class group if need be.
858: * 2) the factorization of a.
859: */
860: static void
861: isintnorm_loop(long i)
862: {
863: if (S[i] == 0) /* sum u[i].f[i] = n[i], do another prime */
864: {
865: int k;
866: if (inext[i] == 0) { test_sol(i); return; }
867:
868: /* there are some primes left */
869: if (Partial) gaffect(Partial[i], Partial[inext[i]-1]);
870: for (k=i+1; k < inext[i]; k++) u[k]=0;
871: i=inext[i]-1;
872: }
873: else if (i == inext[i]-2 || i == Nprimes-1)
874: {
875: /* only one Prime left above prime; change prime, fix u[i+1] */
876: if (S[i] % f[i+1]) return;
877: i++; u[i] = S[i-1] / f[i];
878: if (Partial) fix_Partial(i);
879: if (inext[i]==0) { test_sol(i); return; }
880: }
881:
882: i++; u[i]=0;
883: if (Partial) gaffect(Partial[i-1], Partial[i]);
884: if (i == inext[i-1])
885: { /* change prime */
886: if (inext[i] == i+1 || i == Nprimes) /* only one Prime above p */
887: {
888: S[i]=0; u[i] = n[i] / f[i]; /* we already know this is exact */
889: if (Partial) fix_Partial(i);
890: }
891: else S[i] = n[i];
892: }
893: else S[i] = S[i-1]; /* same prime, different Prime */
894: for(;;)
895: {
896: isintnorm_loop(i); S[i]-=f[i]; if (S[i]<0) break;
897: if (Partial)
898: gaddz(Partial[i], Relations[i], Partial[i]);
899: u[i]++;
900: }
901: }
902:
903: static void
904: get_sol_abs(GEN bnf, GEN a, GEN **ptPrimes)
905: {
906: GEN dec, fact, primes, *Primes, *Fact;
907: long *gcdlist, gcd,nprimes,Ngen,i,j;
908:
909: if (gcmp1(a))
910: {
911: long *sol = new_chunk(Nprimes+1);
912: sindex = 1; normsol = (long**) new_chunk(2);
913: normsol[1] = sol; for (i=1; i<=Nprimes; i++) sol[i]=0;
914: return;
915: }
916:
917: fact=factor(a); primes=(GEN)fact[1];
918: nprimes=lg(primes)-1; sindex = 0;
919: gen_ord = (GEN *) mael3(bnf,8,1,2); Ngen = lg(gen_ord)-1;
920:
921: Fact = (GEN*) new_chunk(1+nprimes);
922: gcdlist = new_chunk(1+nprimes);
923:
924: Nprimes=0;
925: for (i=1; i<=nprimes; i++)
926: {
927: long ldec;
928:
929: dec = primedec(bnf,(GEN)primes[i]); ldec = lg(dec)-1;
930: gcd = itos(gmael(dec,1,4));
931:
932: /* check that gcd_{P|p} f_P | n_p */
933: for (j=2; gcd>1 && j<=ldec; j++)
934: gcd = cgcd(gcd,itos(gmael(dec,j,4)));
935:
936: gcdlist[i]=gcd;
937:
938: if (gcd != 1 && smodis(gmael(fact,2,i),gcd))
939: {
940: if (DEBUGLEVEL>2)
941: { fprintferr("gcd f_P does not divide n_p\n"); flusherr(); }
942: return;
943: }
944: Fact[i] = dec; Nprimes += ldec;
945: }
946:
947: f = new_chunk(1+Nprimes); u = new_chunk(1+Nprimes);
948: n = new_chunk(1+Nprimes); S = new_chunk(1+Nprimes);
949: inext = new_chunk(1+Nprimes);
950: Primes = (GEN*) new_chunk(1+Nprimes);
951: *ptPrimes = Primes;
952:
953: if (Ngen)
954: {
955: Partial = (GEN*) new_chunk(1+Nprimes);
956: Relations = (GEN*) new_chunk(1+Nprimes);
957: }
958: else /* trivial class group, no relations to check */
959: Relations = Partial = NULL;
960:
961: j=0;
962: for (i=1; i<=nprimes; i++)
963: {
964: long k,lim,v, vn=itos(gmael(fact,2,i));
965:
966: gcd = gcdlist[i];
967: dec = Fact[i]; lim = lg(dec);
968: v = (i==nprimes)? 0: j + lim;
969: for (k=1; k < lim; k++)
970: {
971: j++; Primes[j] = (GEN)dec[k];
972: f[j] = itos(gmael(dec,k,4)) / gcd;
973: n[j] = vn / gcd; inext[j] = v;
974: if (Partial)
975: Relations[j] = isprincipal(bnf,Primes[j]);
976: }
977: }
978: if (Partial)
979: {
980: for (i=1; i <= Nprimes; i++)
981: if (!gcmp0(Relations[i])) break;
982: if (i > Nprimes) Partial = NULL; /* all ideals dividing a are principal */
983: }
984: if (Partial)
985: for (i=0; i<=Nprimes; i++) /* Partial[0] needs to be initialized */
986: {
987: Partial[i]=cgetg(Ngen+1,t_COL);
988: for (j=1; j<=Ngen; j++)
989: {
990: Partial[i][j]=lgeti(4);
991: gaffect(gzero,(GEN) Partial[i][j]);
992: }
993: }
994: smax=511; normsol = (long**) new_chunk(smax+1);
995: S[0]=n[1]; inext[0]=1; isintnorm_loop(0);
996: }
997:
998: static long
999: get_unit_1(GEN bnf, GEN pol, GEN *unit)
1000: {
1001: GEN fu;
1002: long j;
1003:
1004: if (DEBUGLEVEL > 2)
1005: fprintferr("looking for a fundamental unit of norm -1\n");
1006:
1007: *unit = gmael3(bnf,8,4,2); /* torsion unit */
1008: if (signe( gnorm(gmodulcp(*unit,pol)) ) < 0) return 1;
1009:
1010: fu = gmael(bnf,8,5); /* fundamental units */
1011: for (j=1; j<lg(fu); j++)
1012: {
1013: *unit = (GEN)fu[j];
1014: if (signe( gnorm(gmodulcp(*unit,pol)) ) < 0) return 1;
1015: }
1016: return 0;
1017: }
1018:
1019: GEN
1020: bnfisintnorm(GEN bnf, GEN a)
1021: {
1022: GEN nf,pol,res,unit,x,id, *Primes;
1023: long sa,i,j,norm_1;
1.2 ! noro 1024: gpmem_t av = avma;
1.1 noro 1025:
1026: bnf = checkbnf(bnf); nf = (GEN)bnf[7]; pol = (GEN)nf[1];
1027: if (typ(a)!=t_INT)
1028: err(talker,"expected an integer in bnfisintnorm");
1029: sa = signe(a);
1.2 ! noro 1030: if (sa == 0) return _vec(gzero);
! 1031: if (gcmp1(a)) return _vec(gun);
1.1 noro 1032:
1033: get_sol_abs(bnf, absi(a), &Primes);
1034:
1035: res = cgetg(1,t_VEC); unit = NULL;
1036: norm_1 = 0; /* gcc -Wall */
1037: for (i=1; i<=sindex; i++)
1038: {
1.2 ! noro 1039: x = normsol[i];
! 1040: id = idealpow(nf,Primes[1], stoi(x[1]));
! 1041: for (j=2; j<=Nprimes; j++) /* compute prod Primes[i]^u[i] */
! 1042: id = idealmulpowprime(nf,id, Primes[j], stoi(x[j]));
1.1 noro 1043: x = (GEN) isprincipalgenforce(bnf,id)[2];
1044: x = gmul((GEN)nf[7],x); /* x possible solution */
1045: if (signe(gnorm(gmodulcp(x,(GEN)nf[1]))) != sa)
1046: {
1047: if (! unit) norm_1 = get_unit_1(bnf,pol,&unit);
1048: if (norm_1) x = gmul(unit,x);
1049: else
1050: {
1.2 ! noro 1051: if (DEBUGLEVEL > 2) fprintferr("%Z eliminated because of sign\n",x);
1.1 noro 1052: x = NULL;
1053: }
1054: }
1055: if (x) res = concatsp(res, gmod(x,pol));
1056: }
1057: return gerepilecopy(av,res);
1058: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to thue.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 / modules