Loic Grenie on Wed, 22 Nov 2006 19:49:41 +0100

 Re: rnfkummer modification

```    I've made a few errors in the last submission. Let's hope this one is
better.

I discuss the algorithm and the trick:

- I've not changed the core part of the computations.
- I modify _rnfkummer and rnfkummersimple so that they give only
linearly independant extensions in all = -ell.
- I modify rnfkummersimple so that if all = -1, it follows a
different strategy: it looks only for linearly independant
extensions. Then if it finds that the subgroup can be written
as a linea combination of these extensions, it does the
corresponding product of the beta's.
The problem is that rnfnormgroup gives a subgroup, which can
be identified with a point in the projective space. To have
the vector in the vector space, I use a second rnfnormgroup
for <the first be>*be. Since <the first be> and be are linearly
independant, I can adjust the scalar. This means that I need
(at most) n^2-n rnfnormgroup(), with n=#bnr.cyc.

- The trick: I've tried to put the basis of K so that ok_congruence
returns 1 for as many vectors of the basis as possible. This is
done before the main loop (it starts with a
K = gmodulo(FpM_ker(M), gell);). (It's mostly a glorified Gauss
elimination). Then I do a first pass in the main loop and try to
pick a quick basis. It that fails, I do a second pass enumerating
all the vectors.

On a Pentium M 1.73GHz, under WinXP+cygwin, with a bnr such that

bnr.pol == y^16 - y^14 + y^12 + 4*y^10 - 3*y^8 + 2*y^6 + 4*y^4 - 2*y^2 + 1
bnr.mod == 2^2*3^3*5^2*11^2

the command

ext = rnfkummer(bnr, matdiagonal(vector(#bnr.cyc, i, 1+2*(i==19))), -1)

finishes in 1mn, 7,265 ms while

rnfnormgroup(bnr, ext)

requires 3,484 ms. (I think the bnfinit has been done with setrand(1)
but I can't remember; this precise diagonal matrix is one of the 3
that have the largest ray class with 21 generators, which means that
rnfkummer(bnr, this_matrix) might need up to (3^21-1)/2
rnfnormgroup()).

I hope I've not made too many errors. I am rather confident that the
basic ideas are correct. I'd like to see that included in Pari if it
is possible.

Loïc```
```Index: src/modules/kummer.c
===================================================================
RCS file: /home/cvs/pari/src/modules/kummer.c,v
retrieving revision 1.130
diff -u -r1.130 kummer.c
--- src/modules/kummer.c	30 Sep 2006 23:43:03 -0000	1.130
+++ src/modules/kummer.c	22 Nov 2006 18:46:09 -0000
@@ -340,9 +340,9 @@
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 (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,6 +509,41 @@
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, ell;
+  GEN col;
+
+  l = lg(subgroup);
+  col = cgetg(l, t_COL);
+  for (i = 1; i < l; i++)
+    /* subgroup is in HNF, diagonal elements are > 0 */
+    if ((ell = gcoeff(subgroup, i, i)[2]) == 1)
+      gel(col, i) = gen_0;
+    else
+    {
+      switch (ell)
+      {
+	case 2:
+	  gel(col, i) = gen_1;
+	  break;
+	case 3:
+	  gel(col, i) = gen_2;
+	  break;
+	default:
+	  gel(col, i) = stoi(ell-1);
+	  break;
+      }
+      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,11 +553,14 @@
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,u,M,K,y,vecMsup,vecW,vecWB,vecBp,msign;
+  GEN res=NULL,u,M,K,y,vecMsup,vecW,vecWB,vecBp,msign;
primlist L;

bnf = gel(bnr,1);
@@ -551,10 +589,10 @@
matP  = cgetg(lSp+1, t_MAT);
for (j=1; j<=lSp; j++)
{
-    GEN e, a, L;
+    GEN 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;
+    gel( matP,j) = gel(L,1);
+    gel(vecBp,j) = gel(L,2);
}
vecWB = shallowconcat(vecW, vecBp);

@@ -580,10 +618,143 @@
lW = lg(vecW);
M = vconcat(M, shallowconcat(zeromat(rc,lW-1), matP));

-  K = FpM_ker(M, gell);
-  dK = lg(K)-1;
+  if (all < 0)
+  {
+    pari_sp av2 = avma;
+    GEN Kidx, Kl;
+    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 = gmodulo(FpM_ker(M, gell), gell); /* Computation mod ell is slow but easier. */
+    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, i, idx)))
+	gel(K, idx) = gdiv(gel(K, idx), gcoeff(K, i, idx));
+      Ki = gel(K, idx);
+      if (gcmp0(gcoeff(K, i, dK)))
+	gel(K, dK) = gadd(gel(K, dK), 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;
+      GEN scalarprod = NULL;
+      if (lg(gmael(vecMsup, i, 1)) != 2)
+	continue;
+      for (j=ffree; --j > 0; )
+	if (!gcmp0(scalarprod = gel(gmul(gel(vecMsup, i), gel(K, j)), 1)))
+	  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-1; 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 (!gcmp1(scalarprod))
+	gel(K, ffree) = gdiv(gel(K, ffree), scalarprod);
+      Ki = gel(K, ffree);
+      if (gcmp0(gmul(gel(vecMsup, i), gel(K, dK))))
+	gel(K, dK) = gadd(gel(K, dK), Ki);
+      for (j = dK; --j > 0; )
+      {
+	if (j == ffree) continue;
+	if (!gcmp0(scalarprod = gel(gmul(gel(vecMsup, i), gel(K, j)), 1)))
+	  gel(K, j) = gsub(gel(K, j), gmul(scalarprod, Ki));
+      }
+    }
+    if (ell == 2)
+    {
+      for (i = ffree, j = ffree - 1; i <= dK && j; i++, j--)
+      {
+	GEN tmp = gel(K, i);
+	gel(K, i) = gel(K, j);
+	gel(K, j) = tmp;
+      }
+    }
+    /* Try to ensure that bzero(y); y[i]=1; gives a good candidate */
+    Kl = gel(K, dK);
+    for (i = 1; i < dK; i++)
+      gel(K, i) = gadd(gel(K, i), Kl);
+    K = gerepileupto(av2, lift(K));
+  }
+  else
+  {
+    K = FpM_ker(M, gell);
+    dK = lg(K)-1;
+  }
y = cgetg(dK+1,t_VECSMALL);
-  res = cgetg(1,t_VEC); /* in case all = 1 */
+  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;
@@ -591,25 +762,111 @@
do
{
pari_sp av = avma;
-      GEN be, P, X = FpC_red(ZM_zc_mul(K, y), gell);
+      GEN be, P=NULL, X;
+      if (all<0)
+      {
+	gel(mat, rk+1) = y;
+	if (Flm_rank(mat, ell) <= rk)
+	  continue;
+      }
+      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;
-        }
+	* and relative discriminant are correct */
+	if (all < 0)
+	{
+	  rk++;
+	  gel(mat, rk) = gclone(y);
+	}
+	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 (firstpass) break;
}
else avma = av;
} while (increment(y, dK, ell));
y[dK--] = 0;
}
-  if (all) return res;
+  } 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,13 +1139,16 @@
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 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,12 +1162,13 @@
if (DEBUGLEVEL) msgTIMER(&t, "[rnfkummer] conductor");
bnr      = gel(p1,2);
subgroup = gel(p1,3);
-  gell = get_gell(bnr,subgroup,all);
+  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. */
@@ -1068,13 +1329,134 @@
if (!M) M = zeromat(1, lSp + lW - 1);
/* step 16 */
if (DEBUGLEVEL>2) fprintferr("Step 16\n");
-  K = FpM_ker(M, gell);
+  if (all < 0)
+  {
+    pari_sp av2 = avma;
+    GEN Kidx, Kl;
+    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 = gmodulo(FpM_ker(M, gell), gell); /* Computation mod ell is slow but easier. */
+    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, i, 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;
+      GEN scalarprod = NULL;
+      if (lg(gmael(vecMsup, i, 1)) != 2)
+	continue;
+      for (j=ffree; --j > 0; )
+	if (!gcmp0(scalarprod = gel(gmul(gel(vecMsup, i), gel(K, j)), 1)))
+	  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-1; 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 (!gcmp1(scalarprod))
+	gel(K, ffree) = gdiv(gel(K, ffree), scalarprod);
+      Ki = gel(K, ffree);
+      if (!gcmp1(scalarprod = gel(gmul(gel(vecMsup, i), gel(K, dK)), 1)))
+	gel(K, dK) = gsub(gel(K, dK), gmul(gsubgs(scalarprod, 1), Ki));
+      for (j = dK; --j > 0; )
+      {
+	if (j == ffree) continue;
+	if (!gcmp0(scalarprod = gel(gmul(gel(vecMsup, i), gel(K, j)), 1)))
+	  gel(K, j) = gsub(gel(K, j), gmul(scalarprod, Ki));
+      }
+    }
+    if (ell == 2)
+    {
+      for (i = ffree, j = ffree - 1; i <= dK && j; i++, j--)
+      {
+	GEN tmp = gel(K, i);
+	gel(K, i) = gel(K, j);
+	gel(K, j) = tmp;
+      }
+    }
+    /* Try to ensure that bzero(y); y[i]=1; gives a good candidate */
+    Kl = gel(K, dK);
+    for (i = 1; i < dK; i++)
+      gel(K, i) = gadd(gel(K, i), Kl);
+    K = gerepileupto(av2, lift(K));
+  }
+  else
+  {
+    K = FpM_ker(M, gell);
+    dK = lg(K)-1;
+  }
/* step 18 & ff */
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 (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,22 +1466,47 @@
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);
+	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 (firstpass) break;
}
} while (increment(y, dK, ell));
y[dK--] = 0;
}
-  if (all) return res;
+  } while (firstpass--);
+  if (all)
+  {
+    if (all<0)
+DONE2:
+      for (i=1;i<=rk;i++)
+	gunclone(gel(mat, i));
+    return res;
+  }
return gen_0; /* FAIL */
}

```