Loic Grenie on Sun, 19 Nov 2006 17:27:29 +0100

 rnfkummer modification

```    Hi again,

I've asked in the pari-users mailing list help for determining a
"basis" of linearly indenpendant Kummer extensions of a base field
(or bnr). I've devised a modification of _rnfkummer and rndkummersimple
that does the following:

- rnfkummer(bnr, subgroup, -ell) will output a list of linearly
independant extensions for the given subgroup.

--> - rnfkummer(bnr, subgroup, -1) will output the same thing as
-->   before, rnfkummer(bnr, subgroup)) i.e. the equation of a Kummer
-->   extension for the given subgroup but using an algorithm which
-->   should be, most of the time,

-->   (sorry for shouting), when going through rnfkummersimple() (it
-->   might apply also to _rnfkummer but I don't know and I've
-->   preferred not to do something that is beyond my limited knowledge).

While the second feature is not very useful to me, the first one is
nearly vital (well, I can always apply the patch myself !). The second
is a incentive for the developpers to look benevolently at my patch !
(I confess that I'm very proud of it, though).

Full set of disclaimers:

- I don't know the theorem that are behind the algorithm. I've just
used the F_ell-vector space structure of the
Gal(compositum(Kummer extensions of deg ell of K))/K (when K
contains all ell-th roots of unity) + a trick to improve the
probability that ok_congruence() will return 1.
- I've tried to check it only on a small set of examples.
- I'm a very bad programmer.

Note: I've not indented two do {...} while (firstpass--); loops
to minimize the set of differences and make clearer what I've changed.

Loïc```
```Index: src/modules/kummer.c
===================================================================
RCS file: /home/cvs/pari/src/modules/kummer.c,v
retrieving revision 1.130
diff -c -r1.130 kummer.c
*** src/modules/kummer.c	30 Sep 2006 23:43:03 -0000	1.130
--- src/modules/kummer.c	19 Nov 2006 15:47:15 -0000
***************
*** 340,348 ****
get_gell(GEN bnr, GEN subgp, long all)
{
GEN gell;
!   if (all)        gell = stoi(all);
!   else if (subgp) gell = det(subgp);
!   else            gell = det(diagonal_i(gmael(bnr,5,2)));
if (typ(gell) != t_INT) pari_err(arither1);
return gell;
}
--- 340,348 ----
get_gell(GEN bnr, GEN subgp, long all)
{
GEN gell;
!   if (all && all != -1) gell = stoi(all);
!   else if (subgp)       gell = det(subgp);
!   else                  gell = det(diagonal_i(gmael(bnr,5,2)));
if (typ(gell) != t_INT) pari_err(arither1);
return gell;
}
***************
*** 509,514 ****
--- 509,537 ----
return x;
}

+ /* A column vector representing a subgroup of prime index */
+ /* !! Not Stack Clean: return pointers to subgroup's elements !! */
+ static GEN
+ grptocol(GEN subgroup)
+ {
+   long l, i, j;
+   GEN col;
+
+   l = lg(subgroup);
+   col = cgetg(l, t_COL);
+   for (i = 1; i < l; i++)
+     if (gcmp1(gcoeff(subgroup, i, i)))
+       gel(col, i) = gen_0;
+     else
+     {
+       gel(col, i) = gen_m1;
+       break;
+     }
+   for (j=i; ++j < l; )
+     gel(col, j) = gcoeff(subgroup, i, j);
+   return col;
+ }
+

/* if all!=0, give all equations of degree 'all'. Assume bnr modulus is the
* conductor */
***************
*** 518,528 ****
long ell, i, j, degK, dK;
long lSml2, lSl2, lSp, rc, lW;
long prec;

GEN bnf,nf,bid,ideal,arch,cycgen;
GEN clgp,cyc;
GEN Sp,listprSp,matP;
!   GEN res,u,M,K,y,vecMsup,vecW,vecWB,vecBp,msign;
primlist L;

bnf = gel(bnr,1);
--- 541,554 ----
long ell, i, j, degK, dK;
long lSml2, lSl2, lSp, rc, lW;
long prec;
+   long rk=0, ncyc=0;
+   GEN mat=NULL, matgrp=NULL, xell;
+   long firstpass = all<0;

GEN bnf,nf,bid,ideal,arch,cycgen;
GEN clgp,cyc;
GEN Sp,listprSp,matP;
!   GEN res=NULL,u,M,K,y,vecMsup,vecW,vecWB,vecBp,msign;
primlist L;

bnf = gel(bnr,1);
***************
*** 551,560 ****
matP  = cgetg(lSp+1, t_MAT);
for (j=1; j<=lSp; j++)
{
!     GEN e, a, L;
L = isprincipalell(bnf,gel(Sp,j), cycgen,u,gell,rc);
!     e = gel(L,1); gel(matP,j) = e;
!     a = gel(L,2); gel(vecBp,j) = a;
}
vecWB = shallowconcat(vecW, vecBp);

--- 577,586 ----
matP  = cgetg(lSp+1, t_MAT);
for (j=1; j<=lSp; j++)
{
!     GEN L;
L = isprincipalell(bnf,gel(Sp,j), cycgen,u,gell,rc);
!     gel( matP,j) = gel(L,1);
!     gel(vecBp,j) = gel(L,2);
}
vecWB = shallowconcat(vecW, vecBp);

***************
*** 580,615 ****
lW = lg(vecW);
M = vconcat(M, shallowconcat(zeromat(rc,lW-1), matP));

!   K = FpM_ker(M, gell);
!   dK = lg(K)-1;
y = cgetg(dK+1,t_VECSMALL);
!   res = cgetg(1,t_VEC); /* in case all = 1 */
while (dK)
{
for (i=1; i<dK; i++) y[i] = 0;
y[i] = 1; /* y = [0,...,0,1,0,...,0], 1 at dK'th position */
do
{
pari_sp av = avma;
!       GEN be, P, X = FpC_red(ZM_zc_mul(K, y), gell);
if (ok_congruence(X, gell, lW, vecMsup) && ok_sign(X, msign, arch))
{/* be satisfies all congruences, x^ell - be is irreducible, signature
!         * and relative discriminant are correct */
!         be = compute_beta(X, vecWB, gell, bnf);
!         be = lift_if_rational(coltoalg(nf, be));
!         P = gsub(monomial(gen_1, ell, 0), be);
!         if (all) res = shallowconcat(res, gerepileupto(av, P));
!         else
!         {
!           if (gequal(rnfnormgroup(bnr,P),subgroup)) return P; /*DONE*/
!           avma = av;
!         }
}
else avma = av;
} while (increment(y, dK, ell));
y[dK--] = 0;
}
!   if (all) return res;
return gen_0;
}

--- 606,854 ----
lW = lg(vecW);
M = vconcat(M, shallowconcat(zeromat(rc,lW-1), matP));

!   if (all < 0)
!   {
!     pari_sp av2 = avma;
!     GEN Kidx;
!     long idx, ffree;
!     /* Reorganize the kernel so that the tests of ok_congruence can be ok
!      * for y[ncyc]=1 and y[1..ncyc]=1 */
!     K = ker(gmodulo(M, gell));
!     dK = lg(K)-1;
!     Kidx = cgetg(dK+1, t_VECSMALL);
!     /* First step: Gauss elimination on vectors lW...lg(M) */
!     for (idx = lg(K), i=lg(M); --i >= lW; )
!     {
!       GEN tmp, Ki;
!       for (j=dK; j > 0; j--) if (!gcmp0(gcoeff(K, i, j))) break;
!       if (!j)
!       {
! 	/* Do our best to ensure that K[dK,i]!=0 */
! 	if (gcmp0(gcoeff(K, i, dK)))
! 	{
! 	  for (j = idx; j < dK; j++)
! 	    if (!gcmp0(gcoeff(K, i, j)) &&
! 		  itos(gcoeff(K, Kidx[j], dK)) < ell - 1)
! 	  gel(K, dK) = gadd(gel(K, dK), gel(K, j));
! 	}
! 	continue;
!       }
!       if (j != --idx)
!       {
! 	tmp = gel(K, j);
! 	gel(K, j) = gel(K, idx);
! 	gel(K, idx) = tmp;
!       }
!       Kidx[idx] = i;
!       if (!gcmp1(gcoeff(K, idx, idx)))
! 	gel(K, idx) = gdiv(gel(K, idx), gcoeff(K, i, idx));
!       Ki = gel(K, idx);
!       if (!gcmp1(gcoeff(K, i, dK)))
! 	gel(K, dK) = gsub(gel(K, dK), gmul(gsubgs(gcoeff(K, i, dK), 1), Ki));
!       for (j = dK; --j > 0; )
!       {
! 	if (j == idx) continue;
! 	if (!gcmp0(gcoeff(K, i, j)))
! 	  gel(K, j) = gsub(gel(K, j), gmul(gcoeff(K, i, j), Ki));
!       }
!     }
!     /* ffree = first vector that is not "free" for the scalar products */
!     ffree = idx;
!     /* Second step: for each hyperplace equation in vecMsup, do the same
!      * thing as before. */
!     for (i=1; i < lg(vecMsup); i++)
!     {
!       GEN tmp, Ki;
!       long scalarprod;
!       for (j=ffree; --j > 0; )
! 	if (lg(gel(vecMsup, i)) == 2 &&
! 	    (scalarprod = itos(lift(gmul(gel(vecMsup, i), gel(K, j))))))
! 	  break;
!       if (!j)
!       {
! 	/* Do our best to ensure that vecMsup.K[dK] != 0 */
! 	if (!signe(lift(gmul(gel(vecMsup, i), gel(K, dK)))))
! 	{
! 	  for (j = ffree; j < dK; j++)
! 	    if (itos(lift(gmul(gel(vecMsup, i), gel(K, j)))) &&
! 		  itos(gcoeff(K, Kidx[j], dK)) < ell - 1)
! 	  gel(K, dK) = gadd(gel(K, dK), gel(K, j));
! 	}
! 	continue;
!       }
!       if (j != --ffree)
!       {
! 	tmp = gel(K, j);
! 	gel(K, j) = gel(K, ffree);
! 	gel(K, ffree) = tmp;
!       }
!       if (!signe(lift(gmul(gel(vecMsup, i), gel(K, dK)))))
! 	gel(K, ffree) = gdiv(gel(K, ffree), gcoeff(K, i, ffree));
!       Ki = gel(K, ffree);
!       if (!gcmp1(gcoeff(K, i, dK)))
! 	gel(K, dK) = gsub(gel(K, dK), gmul(gsubgs(gcoeff(K, i, dK), 1), Ki));
!       for (j = dK; --j > 0; )
!       {
! 	if (j == ffree) continue;
! 	if (!gcmp0(gcoeff(K, i, j)))
! 	  gel(K, j) = gsub(gel(K, j), gmul(gcoeff(K, i, j), Ki));
!       }
!     }
!     if (ell == 2)
!     {
!       for (i = ffree, j = ffree - 1; i < dK; i++, j--)
!       {
! 	GEN tmp = gel(K, i);
! 	gel(K, i) = gel(K, j);
! 	gel(K, j) = tmp;
!       }
!     }
!     K = gerepileupto(av2, lift(K));
!   }
!   else
!   {
!     K = FpM_ker(M, gell);
!     dK = lg(K)-1;
!   }
y = cgetg(dK+1,t_VECSMALL);
!   if (all) res = cgetg(1,t_VEC); /* in case all = 1 */
!   if (all<0)
!   {
!     GEN col;
!     ncyc = dK;
!     col = cgetg(lg(M)+1, t_COL);
!     for (i=1; i <= lg(M); i++)
!       gel(col, i) = gen_0;
!     mat = cgetg(ncyc+1, t_MAT);
!     for (i = 1; i <= ncyc; i++)
!       gel(mat, i) = col;
!     if (all == -1)
!     {
!       long l = lg(gmael(bnr, 5, 2));
!       col = cgetg(l, t_COL);
!       for (i = 1; i < l; i++)
! 	gel(col, i) = gen_0;
!       matgrp = cgetg(ncyc+2, t_MAT);
!       gel(matgrp, 1) = grptocol(subgroup);
!       for (i = 2; i <= ncyc+1; i++)
! 	gel(matgrp, i) = col;
!     }
!     rk = 0;
!   }
!   xell = monomial(gen_1, ell, 0);
!   do {
!   dK = lg(K)-1;
while (dK)
{
for (i=1; i<dK; i++) y[i] = 0;
y[i] = 1; /* y = [0,...,0,1,0,...,0], 1 at dK'th position */
+     if (firstpass) y[ncyc] = 1;
do
{
pari_sp av = avma;
!       GEN be, P=NULL, X = FpC_red(ZM_zc_mul(K, y), gell);
if (ok_congruence(X, gell, lW, vecMsup) && ok_sign(X, msign, arch))
{/* be satisfies all congruences, x^ell - be is irreducible, signature
! 	* and relative discriminant are correct */
! 	if (all<0)
! 	{
! 	  gel(mat, rk+1) = X;
! 	  if (rank(mat) > rk)
! 	  {
! 	    rk++;
! 	    gel(mat, rk) = gclone(X);
! 	  }
! 	  else
! 	    continue;
! 	}
! 	be = compute_beta(X, vecWB, gell, bnf);
! 	be = lift_if_rational(coltoalg(nf, be));
! 	if (all != -1) P = gsub(xell, be);
! 	if (all)
! 	{
! 	  res = shallowconcat(res, gerepileupto(av, (all+1) ? P:gcopy(be)));
! 	}
! 	else
! 	{
! 	  if (gequal(rnfnormgroup(bnr,P),subgroup)) return P; /*DONE*/
! 	  avma = av;
! 	  continue;
! 	}
! 	if (all == -1)
! 	{
! 	  pari_sp av2 = avma;
! 	  GEN Kgrp, colgrp;
! 	  long dKgrp;
! 	  colgrp = grptocol(rnfnormgroup(bnr, gsub(xell, gel(res, rk))));
! 	  if (ell != 2 && rk != 1)
! 	  {
! 	    /* Compute the pesky scalar */
! 	    GEN mat2, K2, multiplicator;
! 	    mat2 = cgetg(4, t_MAT);
! 	    gel(mat2, 1) = gel(matgrp, 2);
! 	    gel(mat2, 2) = colgrp;
! 	    gel(mat2, 3) = grptocol(rnfnormgroup(bnr, gsub(xell, gmul(gel(res, 1), gel(res, rk)))));
! 	    K2 = FpM_ker(mat2, gell);
! 	    if (lg(K2) != 2)
! 	      pari_err(talker, "L'algèbre linéaire m'a tuer");
! 	    K2 = gel(K2, 1);
! 	    if (!gequal(gel(K2, 1), gel(K2, 2)))
! 	    {
! 	      multiplicator = gdiv(gel(K2, 2), gmodulo(gel(K2, 1), gell));
! 	      colgrp = lift(gmul(multiplicator, colgrp));
! 	    }
! 	  }
! 	  gel(matgrp, rk+1) = gclone(colgrp);
! 	  setlg(matgrp, rk+2);
! 	  Kgrp = FpM_ker(matgrp, gell);
! 	  setlg(matgrp, ncyc+2);
! 	  dKgrp = lg(Kgrp)-1;
! 	  if (dKgrp /* == 1 */)
! 	  {
! 	    long l;
! 	    Kgrp = gel(Kgrp, 1);
! 	    l = lg(Kgrp);
! 	    for (i = 2; i < l; i++) {
! 	      GEN pow = gel(Kgrp, i);
!
! 	      if (signe(pow))
! 	      {
! 		if (pow[2] != 1)
! 		  gel(res, i-1) = gpowgs(gel(res, i-1), pow[2]);
! 	      }
! 	      else
! 		  gel(res, i-1) = gen_1;
! 	    }
! 	    be = divide_conquer_prod(res, gmul);
! 	    be = reducebeta(bnf, be, gell);
! 	    be = coltoalg(nf, be);
! 	    res = gsub(xell, be);
! 	    goto DONE2;
! 	  }
! 	  avma = av2;
! 	}
! 	if (all < 0 && rk == ncyc)
! 	  goto DONE2;
! 	if (all < 0) break;
}
else avma = av;
} while (increment(y, dK, ell));
y[dK--] = 0;
}
!   } while (firstpass--);
!   if (all)
!   {
!     if (all<0)
!     {
! DONE2:
!       for (i=1;i<=rk;i++)
! 	gunclone(gel(mat, i));
!       if (all == -1)
! 	for (i=1;i<=rk;i++)
! 	  gunclone(gel(matgrp, i+1));
!     }
!     return res;
!   }
return gen_0;
}

***************
*** 882,894 ****
GEN polnf,bnf,nf,bnfz,nfz,bid,ideal,cycgen,gell,p1,p2,wk,U,vselmer;
GEN clgp,cyc,gen;
GEN Q,idealz,gothf;
!   GEN res,u,M,K,y,vecMsup,vecW,vecWA,vecWB,vecB,vecC,vecAp,vecBp;
GEN matP,Sp,listprSp,Tc,Tv,P;
primlist L;
toK_s T;
tau_s _tau, *tau;
compo_s COMPO;
pari_timer t;

if (DEBUGLEVEL) TIMERstart(&t);
checkbnrgen(bnr);
--- 1121,1136 ----
GEN polnf,bnf,nf,bnfz,nfz,bid,ideal,cycgen,gell,p1,p2,wk,U,vselmer;
GEN clgp,cyc,gen;
GEN Q,idealz,gothf;
!   GEN res=NULL,u,M,K,y,vecMsup,vecW,vecWA,vecWB,vecB,vecC,vecAp,vecBp;
GEN matP,Sp,listprSp,Tc,Tv,P;
primlist L;
toK_s T;
tau_s _tau, *tau;
compo_s COMPO;
pari_timer t;
+   long rk=0, ncyc=0;
+   GEN mat=NULL;
+   long firstpass = all<0;

if (DEBUGLEVEL) TIMERstart(&t);
checkbnrgen(bnr);
***************
*** 902,913 ****
if (DEBUGLEVEL) msgTIMER(&t, "[rnfkummer] conductor");
bnr      = gel(p1,2);
subgroup = gel(p1,3);
!   gell = get_gell(bnr,subgroup,all);
ell = itos(gell);
if (ell == 1) return pol_x(vnf);
if (!uisprime(ell)) pari_err(impl,"kummer for composite relative degree");
if (!umodiu(wk,ell)) return rnfkummersimple(bnr, subgroup, gell, all);

bid = gel(bnr,2);
ideal = gmael(bid,1,1);
/* step 1 of alg 5.3.5. */
--- 1144,1156 ----
if (DEBUGLEVEL) msgTIMER(&t, "[rnfkummer] conductor");
bnr      = gel(p1,2);
subgroup = gel(p1,3);
!   gell = get_gell(bnr,subgroup,all<-1?-all:all);
ell = itos(gell);
if (ell == 1) return pol_x(vnf);
if (!uisprime(ell)) pari_err(impl,"kummer for composite relative degree");
if (!umodiu(wk,ell)) return rnfkummersimple(bnr, subgroup, gell, all);

+   if (all == -1) all = 0;
bid = gel(bnr,2);
ideal = gmael(bid,1,1);
/* step 1 of alg 5.3.5. */
***************
*** 1073,1080 ****
if (DEBUGLEVEL>2) fprintferr("Step 18\n");
dK = lg(K)-1;
y = cgetg(dK+1,t_VECSMALL);
!   res = cgetg(1, t_VEC); /* in case all = 1 */
if (DEBUGLEVEL) msgTIMER(&t, "[rnfkummer] candidate list");
while (dK)
{
for (i=1; i<dK; i++) y[i] = 0;
--- 1316,1337 ----
if (DEBUGLEVEL>2) fprintferr("Step 18\n");
dK = lg(K)-1;
y = cgetg(dK+1,t_VECSMALL);
!   if (all) res = cgetg(1, t_VEC);
if (DEBUGLEVEL) msgTIMER(&t, "[rnfkummer] candidate list");
+   if (all<0)
+   {
+     GEN col;
+     ncyc = dK;
+     col = cgetg(lg(M)+1, t_COL);
+     for (i = 1; i <= lg(M); i++)
+       gel(col, i) = gen_0;
+     mat = cgetg(ncyc+1, t_MAT);
+     for (i = 1; i <= ncyc; i++)
+       gel(mat, i) = col;
+     rk = 0;
+   }
+   do {
+   dK = lg(K)-1;
while (dK)
{
for (i=1; i<dK; i++) y[i] = 0;
***************
*** 1084,1105 ****
GEN H, be, P, X = FpC_red(ZM_zc_mul(K, y), gell);
if (ok_congruence(X, gell, lW, vecMsup))
{
!         be = compute_beta(X, vecWB, gell, bnfz);
!         P = compute_polrel(nfz, &T, be, g, ell);
!         P = lift_if_rational(P);
!         if (DEBUGLEVEL>1) fprintferr("polrel(beta) = %Z\n", P);
!         H = rnfnormgroup(bnr, P);
!         if (!all) {
!           if (gequal(subgroup, H)) return P; /* DONE */
!         } else {
!           if (!gequal(subgroup,H) && conductor(bnr, H, -1) == gen_0) continue;
!         }
!         res = shallowconcat(res, P);
}
} while (increment(y, dK, ell));
y[dK--] = 0;
}
!   if (all) return res;
return gen_0; /* FAIL */
}

--- 1341,1387 ----
GEN H, be, P, X = FpC_red(ZM_zc_mul(K, y), gell);
if (ok_congruence(X, gell, lW, vecMsup))
{
! 	if (all<0)
! 	{
! 	  gel(mat, rk+1) = gclone(X);
! 	  if (rank(mat) > rk)
! 	  {
! 	    rk++;
! 	    gel(mat, rk) = gclone(X);
! 	  }
! 	  else
! 	    continue;
! 	}
! 	be = compute_beta(X, vecWB, gell, bnfz);
! 	P = compute_polrel(nfz, &T, be, g, ell);
! 	P = lift_if_rational(P);
! 	if (DEBUGLEVEL>1) fprintferr("polrel(beta) = %Z\n", P);
! 	H = rnfnormgroup(bnr, P);
! 	if (!all) {
! 	  if (gequal(subgroup, H)) return P; /* DONE */
! 	  continue;
! 	} else {
! 	  if (!gequal(subgroup,H) && conductor(bnr, H, -1) == gen_0) {
! 	    continue;
! 	  }
! 	}
! 	res = shallowconcat(res, P);
! 	if (all < 0 && rk == ncyc)
! 	  goto DONE2;
! 	if (all < 0) break;
}
} while (increment(y, dK, ell));
y[dK--] = 0;
}
!   } while (firstpass--);
!   if (all)
!   {
!     if (all<0)
! DONE2:
!       for (i=1;i<=rk;i++)
! 	gunclone(gel(mat, i));
!     return res;
!   }
return gen_0; /* FAIL */
}

```