diff -ru src/src/basemath/base1.c b/src/basemath/base1.c --- src/src/basemath/base1.c 2013-05-06 16:01:56.000000000 +0200 +++ b/src/basemath/base1.c 2013-05-06 16:49:09.042406927 +0200 @@ -1675,40 +1675,42 @@ ZX_is_better(GEN y, GEN x, GEN *dx) { GEN d = ZX_disc(y); - long cmp = absi_cmp(d, *dx); + int cmp; + if (!*dx) *dx = ZX_disc(x); + cmp = absi_cmp(d, *dx); if (cmp < 0) { *dx = d; return 1; } if (cmp == 0) return cmp_abs_ZX(y, x) < 0; return 0; } -static GEN polred_aux(nfbasic_t *T, GEN *pro, long flag); +static void polredbest_aux(nfbasic_t *T, GEN *pro, GEN *px, GEN *pdx, GEN *pa); /* Seek a simpler, polynomial pol defining the same number field as * x (assumed to be monic at this point) */ static GEN nfpolred(nfbasic_t *T, GEN *pro) { - GEN x = T->x, dx = T->dx, a, z, rev, pow, dpow; + GEN x = T->x, dx, b, rev, pow, dpow; long i, n = degpol(x), v = varn(x); if (n == 1) { T->x = deg1pol_shallow(gen_1, gen_m1, v); *pro = NULL; return pol_1(v); } - z = polred_aux(T, pro, nf_ORIG | nf_RED); - if (typ(z) != t_VEC || !ZX_is_better(gel(z,1),x,&dx)) - return NULL; /* no improvement */ - - rev = QXQ_reverse(gel(z,2), x); - x = gel(z,1); if (DEBUGLEVEL>1) err_printf("xbest = %Ps\n",x); + polredbest_aux(T, pro, &x, &dx, &b); + if (x == T->x) return NULL; /* no improvement */ + if (DEBUGLEVEL>1) err_printf("xbest = %Ps\n",x); /* update T */ + rev = QXQ_reverse(b, T->x); pow = QXQ_powers(rev, n-1, x); pow = Q_remove_denom(pow, &dpow); - a = T->bas; - for (i=2; i<=n; i++) gel(a,i) = QX_ZXQV_eval(gel(a,i), pow, dpow); + for (i=2; i<=n; i++) gel(T->bas,i) = QX_ZXQV_eval(gel(T->bas,i), pow, dpow); (void)Z_issquareall(diviiexact(dx,T->dK), &(T->index)); - T->basden = get_bas_den(a); - T->dx = dx; T->x = x; *pro = NULL; return rev; + T->basden = get_bas_den(T->bas); + T->dx = dx; + T->x = x; + *pro = NULL; /* reset */ + return rev; } /* let bas a t_VEC of QX giving a Z-basis of O_K. Return the index of the @@ -1776,7 +1778,7 @@ x = Q_primpart(x); RgX_check_ZX(x, "nfinit"); if (!ZX_is_irred(x)) pari_err(redpoler, "nfinit"); - if (flag & nf_RED || !gequal1(gel(x,lg(x)-1))) + if (flag & nf_RED || !equali1(gel(x,lg(x)-1))) x = ZX_Q_normalize(x, &(T->lead)); nfmaxord(&S, x, flag, fa); index = S.index; @@ -1822,26 +1824,32 @@ nfinitall(GEN x, long flag, long prec) { const pari_sp av = avma; - GEN nf; + GEN nf, lead; nfbasic_t T; nfbasic_init(x, flag, NULL, &T); nfbasic_add_disc(&T); /* more expensive after set_LLL_basis */ - if (T.lead != gen_1 && !(flag & nf_RED)) + lead = T.lead; + if (lead != gen_1 && !(flag & nf_RED)) { pari_warn(warner,"non-monic polynomial. Result of the form [nf,c]"); flag |= nf_RED | nf_ORIG; } if (flag & nf_RED) { - GEN ro, rev = nfpolred(&T, &ro); + GEN ro, rev; + /* lie to polred: more efficient to update *after* modreverse, than to + * unscale in the polred subsystem */ + T.lead = gen_1; + rev = nfpolred(&T, &ro); nf = nfbasic_to_nf(&T, ro, prec); if (flag & nf_ORIG) { if (!rev) rev = pol_x(varn(T.x)); /* no improvement */ - if (T.lead != gen_1) rev = RgX_Rg_div(rev, T.lead); + if (lead != gen_1) rev = RgX_Rg_div(rev, lead); nf = mkvec2(nf, mkpolmod(rev, T.x)); } + T.lead = lead; /* restore */ } else { GEN ro; set_LLL_basis(&T, &ro, 0.99); nf = nfbasic_to_nf(&T, ro, prec); @@ -1948,7 +1956,7 @@ get_polchar(CG_data *d, GEN x) { return get_pol(d, RgM_RgC_mul(d->ZKembed,x)); } -/* return a defining polynomial for Q(w_i) */ +/* return a defining polynomial for Q(w_k) */ static GEN get_polmin_w(CG_data *d, long k) { @@ -1956,6 +1964,22 @@ if (g) (void)ZX_gcd_all(g, ZX_deriv(g), &g); return g; } +/* return a defining polynomial for Q(w_k+w_l) */ +static GEN +get_polmin_add2(CG_data *d, long k, long l) +{ + GEN g = get_pol(d, RgV_add(gel(d->ZKembed,k), gel(d->ZKembed,l))); + if (g) (void)ZX_gcd_all(g, ZX_deriv(g), &g); + return g; +} +/* return a defining polynomial for Q(w_k-w_l) */ +static GEN +get_polmin_sub2(CG_data *d, long k, long l) +{ + GEN g = get_pol(d, RgV_sub(gel(d->ZKembed,k), gel(d->ZKembed,l))); + if (g) (void)ZX_gcd_all(g, ZX_deriv(g), &g); + return g; +} /* does x generate the correct field ? */ static GEN @@ -2044,13 +2068,67 @@ d->v = varn(T->x); d->r1= T->r1; return prec; } +static void +update(GEN *pai, GEN *pch, nfbasic_t *T, long orig) +{ + GEN ch = *pch, ai = *pai; + if (!ch) + { /* accuracy too low, compute algebraically */ + ch = ZXQ_charpoly(ai, T->x, varn(T->x)); + (void)ZX_gcd_all(ch, ZX_deriv(ch), &ch); + } + if (ZX_canon_neg(ch) && orig) ai = RgX_neg(ai); + if (DEBUGLEVEL>3) err_printf("polred: generator %Ps\n", ch); + if (T->lead != gen_1 && orig) ai = RgX_unscale(ai, ginv(T->lead)); + *pch = ch; *pai = ai; +} +static GEN +findmindisc(GEN y, GEN *pa) +{ + GEN a = *pa, x = gel(y,1), b = gel(a,1), dx = NULL; + long i, l = lg(y); + for (i = 2; i < l; i++) + { + GEN yi = gel(y,i); + if (ZX_is_better(yi,x,&dx)) { x = yi; b = gel(a,i); } + } + *pa = b; return x; +} +/* filter [y,b] from polred_aux: keep a single polynomial of degree n in y + * [ the best wrt discriminant ordering ], but keep all non-primitive + * polynomials */ +static void +filter(GEN y, GEN b, long n) +{ + GEN x, a, dx; + long i, k = 1, l = lg(y); + a = x = dx = NULL; + for (i = 1; i < l; i++) + { + GEN yi = gel(y,i), ai = gel(b,i); + if (degpol(yi) == n) + { + if (dx && !ZX_is_better(yi,x,&dx)) continue; + if (!dx) dx = ZX_disc(yi); + x = yi; a = ai; continue; + } + gel(y,k) = yi; + gel(b,k) = ai; k++; + } + if (dx) + { + gel(y,k) = x; + gel(b,k) = a; k++; + } + setlg(y, k); + setlg(b, k); +} + static GEN -polred_aux(nfbasic_t *T, GEN *pro, long flag) +polred_aux(nfbasic_t *T, GEN *pro, long orig) { GEN b, y, x = T->x; - long i, v = varn(x), l = lg(T->bas); - const long orig = flag & nf_ORIG; - const long nfred = flag & nf_RED; + long maxi, i, j, k, v = varn(x), n = lg(T->bas)-1; nffp_t F; CG_data d; @@ -2058,27 +2136,41 @@ *pro = F.ro; d.ZKembed = F.M; - y = cgetg(l, t_VEC); - b = cgetg(l, t_COL); + /* n + 2 sum_{1 <= i <= n} n-i = n + n(n-1) = n*n */ + y = cgetg(n*n + 1, t_VEC); + b = cgetg(n*n + 1, t_COL); + /* i = 1 */ gel(y,1) = deg1pol_shallow(gen_1, gen_m1, v); gel(b,1) = gen_1; - for (i = 2; i < l; i++) + for (i = k = 2; i <= n; i++) { - GEN ch, ai = gel(T->bas,i); + GEN ch, ai; + ai = gel(T->bas,i); ch = get_polmin_w(&d, i); - /* if accuracy too low, compute algebraically */ - if (!ch) + update(&ai, &ch, T, orig); + gel(y,k) = ch; + gel(b,k) = ai; k++; + } + k = i; + maxi = minss(n, 3); + for (i = 1; i <= maxi; i++) + for (j = i+1; j <= n; j++) { - ch = ZXQ_charpoly(ai, x, v); - (void)ZX_gcd_all(ch, ZX_deriv(ch), &ch); + GEN ch, ai; + ai = gadd(gel(T->bas,i), gel(T->bas,j)); + ch = get_polmin_add2(&d, i, j); + update(&ai, &ch, T, orig); + gel(y,k) = ch; + gel(b,k) = ai; k++; + + ai = gsub(gel(T->bas,i), gel(T->bas,j)); + ch = get_polmin_sub2(&d, i, j); + update(&ai, &ch, T, orig); + gel(y,k) = ch; + gel(b,k) = ai; k++; } - if (ZX_canon_neg(ch) && orig) ai = RgX_neg(ai); - if (nfred && degpol(ch) == l-1) return mkvec2(ch, ai); - if (DEBUGLEVEL>3) err_printf("polred: generator %Ps\n", ch); - if (T->lead != gen_1 && orig) ai = RgX_unscale(ai, ginv(T->lead)); - gel(y,i) = ch; - gel(b,i) = ai; - } + setlg(y, k); + setlg(b, k); filter(y, b, n); if (!orig) return gen_sort_uniq(y, (void*)cmpii, &gen_cmp_RgX); (void)sort_factor_pol(mkmat2(y, b), cmpii); settyp(y, t_COL); return mkmat2(b, y); @@ -2089,10 +2181,58 @@ { pari_sp av = avma; GEN ro; - nfbasic_t T; nfbasic_init(x, flag & (nf_PARTIALFACT|nf_RED), fa, &T); + nfbasic_t T; nfbasic_init(x, flag & nf_PARTIALFACT, fa, &T); return gerepilecopy(av, polred_aux(&T, &ro, flag & nf_ORIG)); } +/* finds "best" polynomial in polred_aux list, defaulting to T->x if none of + * them is primitive. *px is the ZX, characteristic polynomial of *pb, *pdx + * its discriminant. + * Set *pro = polroots(T->x) [ NOT *px ], in case caller needs it. */ +static void +polredbest_aux(nfbasic_t *T, GEN *pro, GEN *px, GEN *pdx, GEN *pb) +{ + GEN a, v, y, x = T->x, b = pol_x(varn(x)); /* default values */ + long i, l, n = degpol(x); + v = polred_aux(T, pro, nf_ORIG); + *pdx = T->dx; + y = gel(v,2); + a = gel(v,1); l = lg(a); + for (i=1; i