=================================================================== RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/alglin1.c,v retrieving revision 1.1.1.1 retrieving revision 1.2 diff -u -p -r1.1.1.1 -r1.2 --- OpenXM_contrib/pari-2.2/src/basemath/Attic/alglin1.c 2001/10/02 11:17:00 1.1.1.1 +++ OpenXM_contrib/pari-2.2/src/basemath/Attic/alglin1.c 2002/09/11 07:26:48 1.2 @@ -1,4 +1,4 @@ -/* $Id: alglin1.c,v 1.1.1.1 2001/10/02 11:17:00 noro Exp $ +/* $Id: alglin1.c,v 1.2 2002/09/11 07:26:48 noro Exp $ Copyright (C) 2000 The PARI group. @@ -23,9 +23,9 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston, /* for linear algebra mod p */ #ifdef LONG_IS_64BIT -# define MASK (0x7fffffff00000000UL) +# define MASK (0xffffffff00000000UL) #else -# define MASK (0x7fff0000UL) +# define MASK (0xffff0000UL) #endif /* 2p^2 < 2^BITS_IN_LONG */ @@ -36,7 +36,7 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston, #endif #define OK_ULONG(p) (lgefint(p) == 3 && u_OK_ULONG(p[2])) extern GEN ZM_init_CRT(GEN Hp, ulong p); -extern ulong u_invmod(ulong x, ulong p); +extern int ZM_incremental_CRT(GEN H, GEN Hp, GEN q, GEN qp, ulong p); /*******************************************************************/ /* */ @@ -75,7 +75,6 @@ gtrans_i(GEN x) return y; } - GEN gtrans(GEN x) { @@ -135,7 +134,8 @@ GEN concatsp3(GEN x, GEN y, GEN z) { long i, lx = lg(x), ly = lg(y), lz = lg(z); - GEN r = cgetg(lx+ly+lz-2, t_MAT), t = r; + GEN t, r = cgetg(lx+ly+lz-2, t_MAT); + t = r; for (i=1; i b) return gzero; + if (b > lg(v)-1) err(talker,"incorrect length in sum"); + p = (GEN)v[a]; + for (i=a+1; i<=b; i++) p = gadd(p, (GEN)v[i]); + return p; +} + /*******************************************************************/ /* */ /* SCALAR-MATRIX OPERATIONS */ @@ -668,15 +719,16 @@ GEN zeromat(long m, long n) { GEN y = cgetg(n+1,t_MAT); - long i; for (i=1; i<=n; i++) y[i]=(long)zerocol(m); + GEN v = zerocol(m); + long i; for (i=1; i<=n; i++) y[i]=(long)v; return y; } GEN -gscalcol(GEN x, long n) -{ - GEN y=gscalcol_proto(gzero,gzero,n); - if (n) y[1]=lcopy(x); +gscalcol(GEN x, long n) +{ + GEN y=gscalcol_proto(gzero,gzero,n); + if (n) y[1]=lcopy(x); return y; } @@ -847,22 +899,42 @@ mattodiagonal(GEN m) GEN gaddmat(GEN x, GEN y) { - long ly,dy,i,j; - GEN z; + long l,d,i,j; + GEN z, cz,cy; - ly=lg(y); if (ly==1) return cgetg(1,t_MAT); - dy=lg(y[1]); - if (typ(y)!=t_MAT || ly!=dy) err(mattype1,"gaddmat"); - z=cgetg(ly,t_MAT); - for (i=1; i0; i--) { av = avma; m = mulii(negi((GEN)b[i]),t); for (j=i+1; j<=n; j++) m = addii(m, mulii(gcoeff(A,i,j),(GEN) u[j])); - u[i] = (long)gerepileuptoint(av, divii(negi(m), gcoeff(A,i,i))); + u[i] = lpileuptoint(av, divii(negi(m), gcoeff(A,i,i))); } } return c; } +/* A, B integral upper HNF. A^(-1) B integral ? */ +int +hnfdivide(GEN A, GEN B) +{ + gpmem_t av = avma; + long n = lg(A)-1, i,j,k; + GEN u, b, m, r; + + if (!n) return 1; + u = cgetg(n+1, t_COL); + for (k=1; k<=n; k++) + { + b = (GEN)B[k]; + m = (GEN)b[k]; + u[k] = ldvmdii(m, gcoeff(A,k,k), &r); + if (r != gzero) { avma = av; return 0; } + for (i=k-1; i>0; i--) + { + m = negi((GEN)b[i]); + for (j=i+1; j<=k; j++) + m = addii(m, mulii(gcoeff(A,i,j),(GEN) u[j])); + m = dvmdii(m, gcoeff(A,i,i), &r); + if (r != gzero) { avma = av; return 0; } + u[i] = lnegi(m); + } + } + avma = av; return 1; +} + +/* A upper HNF, b integral vector. Return A^(-1) b if integral, + * NULL otherwise. */ GEN +hnf_invimage(GEN A, GEN b) +{ + gpmem_t av = avma, av2; + long n = lg(A)-1, i,j; + GEN u, m, r; + + if (!n) return NULL; + u = cgetg(n+1, t_COL); + m = (GEN)b[n]; + u[n] = ldvmdii(m, gcoeff(A,n,n), &r); + if (r != gzero) { avma = av; return NULL; } + for (i=n-1; i>0; i--) + { + av2 = avma; + m = negi((GEN)b[i]); + for (j=i+1; j<=n; j++) + m = addii(m, mulii(gcoeff(A,i,j),(GEN) u[j])); + m = dvmdii(m, gcoeff(A,i,i), &r); + if (r != gzero) { avma = av; return NULL; } + u[i] = lpileuptoint(av2, negi(m)); + } + return u; +} + +GEN gauss_get_col(GEN a, GEN b, GEN p, long li) { GEN m, u=cgetg(li+1,t_COL); @@ -948,10 +1074,11 @@ gauss_get_col(GEN a, GEN b, GEN p, long li) u[li] = ldiv((GEN)b[li], p); for (i=li-1; i>0; i--) { + gpmem_t av = avma; m = gneg_i((GEN)b[i]); for (j=i+1; j<=li; j++) m = gadd(m, gmul(gcoeff(a,i,j), (GEN)u[j])); - u[i] = ldiv(gneg_i(m), gcoeff(a,i,i)); + u[i] = lpileupto(av, gdiv(gneg_i(m), gcoeff(a,i,i))); } return u; } @@ -965,38 +1092,79 @@ Fp_gauss_get_col(GEN a, GEN b, GEN piv, long li, GEN p u[li] = lresii(mulii((GEN)b[li], mpinvmod(piv,p)), p); for (i=li-1; i>0; i--) { + gpmem_t av = avma; m = (GEN)b[i]; for (j=i+1; j<=li; j++) m = subii(m, mulii(gcoeff(a,i,j), (GEN)u[j])); m = resii(m, p); - u[i] = lresii(mulii(m, mpinvmod(gcoeff(a,i,i), p)), p); + u[i] = lpileuptoint(av, resii(mulii(m, mpinvmod(gcoeff(a,i,i), p)), p)); } return u; } +GEN +Fq_gauss_get_col(GEN a, GEN b, GEN piv, long li, GEN T, GEN p) +{ + GEN m, u=cgetg(li+1,t_COL); + long i,j; + u[li] = (long)FpXQ_mul((GEN)b[li], FpXQ_inv(piv,T,p), T,p); + for (i=li-1; i>0; i--) + { + gpmem_t av = avma; + m = (GEN)b[i]; + for (j=i+1; j<=li; j++) + m = gsub(m, gmul(gcoeff(a,i,j), (GEN)u[j])); + m = FpX_res(m, T,p); + u[i] = lpileupto(av, FpXQ_mul(m, FpXQ_inv(gcoeff(a,i,i), T,p), T,p)); + } + return u; +} + /* assume -p < a < p, return 1/a mod p */ static long u_Fp_inv(long a, long p) { if (a < 0) a = p + a; /* pb with ulongs < 0 */ - return u_invmod(a,p); + return (long)u_invmod((ulong)a,(ulong)p); } +/* assume 0 <= a[i,i] < p */ GEN -u_Fp_gauss_get_col(GEN a, GEN b, long piv, long li, long p) +u_Fp_gauss_get_col(GEN a, GEN b, ulong piv, long li, ulong p) { - GEN u=cgetg(li+1,t_VECSMALL); - long i,j, m; + GEN u = cgetg(li+1,t_VECSMALL); + ulong m; + long i,j; - u[li] = (b[li] * u_Fp_inv(piv,p)) % p; - if (u[li] < 0) u[li] += p; - for (i=li-1; i>0; i--) + m = b[li] % p; + if (u_OK_ULONG(p)) { - m = b[i]; - for (j=i+1; j<=li; j++) { m -= coeff(a,i,j) * u[j]; if (m & MASK) m %= p; } - u[i] = ((m%p) * u_Fp_inv(coeff(a,i,i), p)) % p; - if (u[i] < 0) u[i] += p; + u[li] = (m * u_Fp_inv(piv,p)) % p; + for (i=li-1; i>0; i--) + { + m = p - b[i]%p; + for (j=i+1; j<=li; j++) { + if (m & HIGHBIT) m %= p; + m += coeff(a,i,j) * u[j]; + } + m %= p; + if (m) m = ((p-m) * u_invmod((ulong)coeff(a,i,i), p)) % p; + u[i] = m; + } } + else + { + u[li] = mulssmod(m, u_Fp_inv(piv,p), p); + for (i=li-1; i>0; i--) + { + m = p - b[i]%p; + for (j=i+1; j<=li; j++) + m += mulssmod(coeff(a,i,j), u[j], p); + m %= p; + if (m) m = mulssmod(p-m, u_invmod((ulong)coeff(a,i,i), p), p); + u[i] = m; + } + } return u; } @@ -1015,16 +1183,37 @@ _Fp_addmul(GEN b, long k, long i, GEN m, GEN p) b[k] = laddii((GEN)b[k], mulii(m, (GEN)b[i])); } +/* same, reduce mod (T,p) */ static void -_u_Fp_addmul(GEN b, long k, long i, long m, long p) +_Fq_addmul(GEN b, long k, long i, GEN m, GEN T, GEN p) { - long a; - if (b[i] & MASK) b[i] %= p; - a = b[k] + m * b[i]; - if (a & MASK) a %= p; - b[k] = a; + b[i] = (long)FpX_res((GEN)b[i], T,p); + b[k] = (long)ladd((GEN)b[k], gmul(m, (GEN)b[i])); } +/* assume m < p && u_OK_ULONG(p) && (! (b[i] & MASK)) */ +static void +_u_Fp_addmul_OK(GEN b, long k, long i, ulong m, ulong p) +{ + b[k] += m * b[i]; + if (b[k] & MASK) b[k] %= p; +} +/* assume m < p */ +static void +_u_Fp_addmul(GEN b, long k, long i, ulong m, ulong p) +{ + b[i] %= p; + b[k] += mulssmod(m, b[i], p); + if (b[k] & MASK) b[k] %= p; +} +/* same m = 1 */ +static void +_u_Fp_add(GEN b, long k, long i, ulong p) +{ + b[k] += b[i]; + if (b[k] & MASK) b[k] %= p; +} + /* Gaussan Elimination. Compute a^(-1)*b * b is a matrix or column vector, NULL meaning: take the identity matrix * If a and b are empty, the result is the empty matrix. @@ -1037,7 +1226,9 @@ _u_Fp_addmul(GEN b, long k, long i, long m, long p) GEN gauss_intern(GEN a, GEN b) { - long inexact,iscol,i,j,k,av,lim,li,bco, aco = lg(a)-1; + gpmem_t av, lim; + long i,j,k,li,bco, aco = lg(a)-1; + int inexact, iscol; GEN p,m,u; if (typ(a)!=t_MAT) err(mattype1,"gauss"); @@ -1048,15 +1239,13 @@ gauss_intern(GEN a, GEN b) if (DEBUGLEVEL) err(warner,"in Gauss lg(a)=1 lg(b)=%ld", b?1:-1); return cgetg(1,t_MAT); } - av=avma; lim=stack_lim(av,1); + av = avma; lim = stack_lim(av,1); li = lg(a[1])-1; if (li != aco && (li < aco || b)) err(mattype1,"gauss"); a = dummycopy(a); - b = check_b(b,li); bco = lg(b)-1; + b = check_b(b,li, &iscol); bco = lg(b)-1; inexact = use_maximal_pivot(a); - iscol = (typ(b)==t_COL); - if(DEBUGLEVEL>4) - fprintferr("Entering gauss with inexact=%ld iscol=%ld\n",inexact,iscol); + if(DEBUGLEVEL>4) fprintferr("Entering gauss with inexact=%ld\n",inexact); p = NULL; /* gcc -Wall */ for (i=1; i<=aco; i++) @@ -1083,51 +1272,32 @@ gauss_intern(GEN a, GEN b) if (k != i) { for (j=i; j<=aco; j++) swap(coeff(a,i,j), coeff(a,k,j)); - if (iscol) { swap(b[i],b[k]); } - else - for (j=1; j<=bco; j++) swap(coeff(b,i,j), coeff(b,k,j)); + for (j=1; j<=bco; j++) swap(coeff(b,i,j), coeff(b,k,j)); p = gcoeff(a,i,i); } if (i == aco) break; for (k=i+1; k<=li; k++) { - m=gcoeff(a,k,i); + m = gcoeff(a,k,i); if (!gcmp0(m)) { m = gneg_i(gdiv(m,p)); for (j=i+1; j<=aco; j++) _addmul((GEN)a[j],k,i,m); - if (iscol) _addmul(b,k,i,m); - else - for (j=1; j<=bco; j++) _addmul((GEN)b[j],k,i,m); + for (j=1; j<=bco; j++) _addmul((GEN)b[j],k,i,m); } } if (low_stack(lim, stack_lim(av,1))) { - GEN *gptr[2]; gptr[0]=&a; gptr[1]=&b; if(DEBUGMEM>1) err(warnmem,"gauss. i=%ld",i); - gerepilemany(av,gptr,2); + gerepileall(av,2, &a,&b); } } if(DEBUGLEVEL>4) fprintferr("Solving the triangular system\n"); - if (iscol) u = gauss_get_col(a,b,p,aco); - else - { - long av1 = avma; - lim = stack_lim(av1,1); u=cgetg(bco+1,t_MAT); - for (j=2; j<=bco; j++) u[j] = zero; /* dummy */ - for (j=1; j<=bco; j++) - { - u[j] = (long)gauss_get_col(a,(GEN)b[j],p,aco); - if (low_stack(lim, stack_lim(av1,1))) - { - if(DEBUGMEM>1) err(warnmem,"gauss[2]. j=%ld", j); - u = gerepilecopy(av1, u); - } - } - } - return gerepilecopy(av, u); + u = cgetg(bco+1,t_MAT); + for (j=1; j<=bco; j++) u[j] = (long)gauss_get_col(a,(GEN)b[j],p,aco); + return gerepilecopy(av, iscol? (GEN)u[1]: u); } GEN @@ -1138,29 +1308,32 @@ gauss(GEN a, GEN b) return z; } -GEN -u_FpM_gauss(GEN a, GEN b, long p) +/* destroy a, b */ +static GEN +u_FpM_gauss_sp(GEN a, GEN b, ulong p) { - long piv,m,iscol,i,j,k,li,bco, aco = lg(a)-1; + long iscol,i,j,k,li,bco, aco = lg(a)-1; + ulong piv, m; GEN u; if (!aco) return cgetg(1,t_MAT); li = lg(a[1])-1; bco = lg(b)-1; - iscol = (typ(b)==t_COL); + iscol = (typ(b)!=t_MAT); + if (iscol) b = col_to_mat(b); piv = 0; /* gcc -Wall */ for (i=1; i<=aco; i++) { /* k is the line where we find the pivot */ - coeff(a,i,i) = coeff(a,i,i) % p; - piv = coeff(a,i,i); k = i; + piv = coeff(a,i,i) % p; + coeff(a,i,i) = piv; k = i; if (!piv) { for (k++; k <= li; k++) { coeff(a,k,i) %= p; if (coeff(a,k,i)) break; - } + } if (k>li) return NULL; } @@ -1168,41 +1341,65 @@ u_FpM_gauss(GEN a, GEN b, long p) if (k != i) { for (j=i; j<=aco; j++) swap(coeff(a,i,j), coeff(a,k,j)); - if (iscol) { swap(b[i],b[k]); } - else - for (j=1; j<=bco; j++) swap(coeff(b,i,j), coeff(b,k,j)); + for (j=1; j<=bco; j++) swap(coeff(b,i,j), coeff(b,k,j)); piv = coeff(a,i,i); } if (i == aco) break; for (k=i+1; k<=li; k++) { - coeff(a,k,i) %= p; - m = coeff(a,k,i); coeff(a,k,i) = 0; + m = coeff(a,k,i) % p; coeff(a,k,i) = 0; if (m) { - m = - (m * u_Fp_inv(piv,p)) % p; - for (j=i+1; j<=aco; j++) _u_Fp_addmul((GEN)a[j],k,i,m, p); - if (iscol) _u_Fp_addmul(b,k,i,m, p); - else - for (j=1; j<=bco; j++) _u_Fp_addmul((GEN)b[j],k,i,m, p); + m = mulssmod(m, u_invmod(piv,p), p); + m = p - m; + if (u_OK_ULONG(p)) + { + for (j=i+1; j<=aco; j++) _u_Fp_addmul_OK((GEN)a[j],k,i,m, p); + for (j=1; j<=bco; j++) _u_Fp_addmul_OK((GEN)b[j],k,i,m, p); + } + else + { + for (j=i+1; j<=aco; j++) _u_Fp_addmul((GEN)a[j],k,i,m, p); + for (j=1; j<=bco; j++) _u_Fp_addmul((GEN)b[j],k,i,m, p); + } } } } - if (iscol) u = u_Fp_gauss_get_col(a,b,piv,aco,p); - else - { - u=cgetg(bco+1,t_MAT); - for (j=1; j<=bco; j++) - u[j] = (long)u_Fp_gauss_get_col(a,(GEN)b[j],piv,aco,p); - } - return u; + u = cgetg(bco+1,t_MAT); + for (j=1; j<=bco; j++) u[j] = (long)u_Fp_gauss_get_col(a,(GEN)b[j],piv,aco,p); + return iscol? (GEN)u[1]: u; } +static GEN +u_idmat(long n) +{ + GEN y = cgetg(n+1,t_MAT); + long i; + if (n < 0) err(talker,"negative size in u_idmat"); + for (i=1; i<=n; i++) { y[i] = (long)vecsmall_const(n, 0); coeff(y, i,i) = 1; } + return y; +} + GEN +u_FpM_gauss(GEN a, GEN b, ulong p) { + return u_FpM_gauss_sp(dummycopy(a), dummycopy(b), p); +} +static GEN +u_FpM_inv_sp(GEN a, ulong p) { + return u_FpM_gauss_sp(a, u_idmat(lg(a)-1), p); +} +GEN +u_FpM_inv(GEN a, ulong p) { + return u_FpM_inv_sp(dummycopy(a), p); +} + +GEN FpM_gauss(GEN a, GEN b, GEN p) { - long iscol,i,j,k,av,lim,li,bco, aco = lg(a)-1; + gpmem_t av,lim; + long i,j,k,li,bco, aco = lg(a)-1; + int iscol; GEN piv,m,u; if (typ(a)!=t_MAT) err(mattype1,"gauss"); @@ -1215,19 +1412,19 @@ FpM_gauss(GEN a, GEN b, GEN p) } li = lg(a[1])-1; if (li != aco && (li < aco || b)) err(mattype1,"gauss"); - b = check_b(b,li); av = avma; + b = check_b(b,li, &iscol); av = avma; if (OK_ULONG(p)) { ulong pp=p[2]; a = u_Fp_FpM(a, pp); b = u_Fp_FpM(b, pp); - u = u_FpM_gauss(a,b, pp); - return gerepileupto(av, small_to_mat(u)); + u = u_FpM_gauss_sp(a,b, pp); + u = iscol? small_to_col((GEN)u[1]): small_to_mat(u); + return gerepileupto(av, u); } lim = stack_lim(av,1); a = dummycopy(a); bco = lg(b)-1; - iscol = (typ(b)==t_COL); piv = NULL; /* gcc -Wall */ for (i=1; i<=aco; i++) { @@ -1240,7 +1437,7 @@ FpM_gauss(GEN a, GEN b, GEN p) { coeff(a,k,i) = lresii(gcoeff(a,k,i), p); if (signe(coeff(a,k,i))) break; - } + } if (k>li) return NULL; } @@ -1248,9 +1445,7 @@ FpM_gauss(GEN a, GEN b, GEN p) if (k != i) { for (j=i; j<=aco; j++) swap(coeff(a,i,j), coeff(a,k,j)); - if (iscol) { swap(b[i],b[k]); } - else - for (j=1; j<=bco; j++) swap(coeff(b,i,j), coeff(b,k,j)); + for (j=1; j<=bco; j++) swap(coeff(b,i,j), coeff(b,k,j)); piv = gcoeff(a,i,i); } if (i == aco) break; @@ -1264,9 +1459,7 @@ FpM_gauss(GEN a, GEN b, GEN p) m = mulii(m, mpinvmod(piv,p)); m = resii(negi(m), p); for (j=i+1; j<=aco; j++) _Fp_addmul((GEN)a[j],k,i,m, p); - if (iscol) _Fp_addmul(b,k,i,m, p); - else - for (j=1; j<=bco; j++) _Fp_addmul((GEN)b[j],k,i,m, p); + for (j=1 ; j<=bco; j++) _Fp_addmul((GEN)b[j],k,i,m, p); } } if (low_stack(lim, stack_lim(av,1))) @@ -1278,42 +1471,89 @@ FpM_gauss(GEN a, GEN b, GEN p) } if(DEBUGLEVEL>4) fprintferr("Solving the triangular system\n"); - if (iscol) u = Fp_gauss_get_col(a,b,piv,aco,p); - else + u = cgetg(bco+1,t_MAT); + for (j=1; j<=bco; j++) u[j] = (long)Fp_gauss_get_col(a,(GEN)b[j],piv,aco,p); + return gerepilecopy(av, iscol? (GEN)u[1]: u); +} +GEN +FqM_gauss(GEN a, GEN b, GEN T, GEN p) +{ + gpmem_t av,lim; + long i,j,k,li,bco, aco = lg(a)-1; + int iscol; + GEN piv,m,u; + + if (!T) return FpM_gauss(a,b,p); + + if (typ(a)!=t_MAT) err(mattype1,"gauss"); + if (b && typ(b)!=t_COL && typ(b)!=t_MAT) err(typeer,"gauss"); + if (!aco) { - long av1 = avma; - lim = stack_lim(av1,1); u=cgetg(bco+1,t_MAT); - for (j=2; j<=bco; j++) u[j] = zero; /* dummy */ - for (j=1; j<=bco; j++) + if (b && lg(b)!=1) err(consister,"gauss"); + if (DEBUGLEVEL) err(warner,"in Gauss lg(a)=1 lg(b)=%ld", b?1:-1); + return cgetg(1,t_MAT); + } + li = lg(a[1])-1; + if (li != aco && (li < aco || b)) err(mattype1,"gauss"); + b = check_b(b,li,&iscol); av = avma; + lim = stack_lim(av,1); + a = dummycopy(a); + bco = lg(b)-1; + piv = NULL; /* gcc -Wall */ + for (i=1; i<=aco; i++) + { + /* k is the line where we find the pivot */ + coeff(a,i,i) = (long)FpX_res(gcoeff(a,i,i), T,p); + piv = gcoeff(a,i,i); k = i; + if (!signe(piv)) { - u[j] = (long)Fp_gauss_get_col(a,(GEN)b[j],piv,aco,p); - if (low_stack(lim, stack_lim(av1,1))) + for (k++; k <= li; k++) { - if(DEBUGMEM>1) err(warnmem,"gauss[2]. j=%ld", j); - u = gerepilecopy(av1, u); + coeff(a,k,i) = (long)FpX_res(gcoeff(a,k,i), T,p); + if (signe(coeff(a,k,i))) break; } + if (k>li) return NULL; } + + /* if (k!=i), exchange the lines s.t. k = i */ + if (k != i) + { + for (j=i; j<=aco; j++) swap(coeff(a,i,j), coeff(a,k,j)); + for (j=1; j<=bco; j++) swap(coeff(b,i,j), coeff(b,k,j)); + piv = gcoeff(a,i,i); + } + if (i == aco) break; + + for (k=i+1; k<=li; k++) + { + coeff(a,k,i) = (long)FpX_res(gcoeff(a,k,i), T,p); + m = gcoeff(a,k,i); coeff(a,k,i) = zero; + if (signe(m)) + { + m = FpXQ_mul(m, FpXQ_inv(piv,T,p),T,p); + m = resii(negi(m), p); + for (j=i+1; j<=aco; j++) _Fq_addmul((GEN)a[j],k,i,m, T,p); + for (j=1; j<=bco; j++) _Fq_addmul((GEN)b[j],k,i,m, T,p); + } + } + if (low_stack(lim, stack_lim(av,1))) + { + GEN *gptr[2]; gptr[0]=&a; gptr[1]=&b; + if(DEBUGMEM>1) err(warnmem,"FpM_gauss. i=%ld",i); + gerepilemany(av,gptr,2); + } } - return gerepilecopy(av, u); + + if(DEBUGLEVEL>4) fprintferr("Solving the triangular system\n"); + u = cgetg(bco+1,t_MAT); + for (j=1; j<=bco; j++) u[j] = (long)Fq_gauss_get_col(a,(GEN)b[j],piv,aco,T,p); + return gerepilecopy(av, iscol? (GEN)u[1]: u); } + GEN FpM_inv(GEN x, GEN p) { return FpM_gauss(x, NULL, p); } -static GEN -u_idmat(long n) -{ - GEN y = cgetg(n+1,t_MAT); - long i,j; - if (n < 0) err(talker,"negative size in u_idmat"); - for (i=1; i<=n; i++) - { - y[i]=lgetg(n+1,t_VECSMALL); - for (j=1; j<=n; j++) coeff(y,j,i) = (i==j)? 1: 0; - } - return y; -} - /* set y = x * y */ #define HIGHWORD(a) ((a) >> BITS_IN_HALFULONG) static GEN @@ -1332,11 +1572,12 @@ u_FpM_Fp_mul_ip(GEN y, long x, long p) } /* M integral, dM such that M' = dM M^-1 is integral [e.g det(M)]. Return M' */ -GEN +GEN ZM_inv(GEN M, GEN dM) { - GEN ID,Hp,q,H; - ulong p,dMp, av2, av = avma, lim = stack_lim(av,1); + gpmem_t av2, av = avma, lim = stack_lim(av,1); + GEN Hp,q,H; + ulong p,dMp; byteptr d = diffptr; long lM = lg(M), stable = 0; @@ -1346,16 +1587,14 @@ ZM_inv(GEN M, GEN dM) av2 = avma; H = NULL; d += 3000; /* 27449 = prime(3000) */ - for(p = 27449; ; p+= *d++) + for(p = 27449; ; ) { - if (!*d) err(primer1); dMp = umodiu(dM,p); - if (!dMp) continue; - ID = u_idmat(lM-1); - Hp = u_FpM_gauss(u_Fp_FpM(M,p), ID, p); + if (!dMp) goto repeat; + Hp = u_FpM_inv_sp(u_Fp_FpM(M,p), p); if (dMp != 1) Hp = u_FpM_Fp_mul_ip(Hp, dMp, p); - if (!H) + if (!H) { H = ZM_init_CRT(Hp, p); q = utoi(p); @@ -1368,27 +1607,28 @@ ZM_inv(GEN M, GEN dM) } if (DEBUGLEVEL>5) msgtimer("inverse mod %ld (stable=%ld)", p,stable); if (stable && isscalarmat(gmul(M, H), dM)) break; /* DONE */ - + if (low_stack(lim, stack_lim(av,2))) { GEN *gptr[2]; gptr[0] = &H; gptr[1] = &q; if (DEBUGMEM>1) err(warnmem,"ZM_inv"); gerepilemany(av2,gptr, 2); } - + repeat: + NEXT_PRIME_VIADIFF_CHECK(p,d); } if (DEBUGLEVEL>5) msgtimer("ZM_inv done"); return gerepilecopy(av, H); } /* same as above, M rational */ -GEN +GEN QM_inv(GEN M, GEN dM) { - ulong av = avma; - GEN cM = content(M); - if (is_pm1(cM)) { avma = av; return ZM_inv(M,dM); } - return gerepileupto(av, ZM_inv(gdiv(M,cM), gdiv(dM,cM))); + gpmem_t av = avma; + GEN cM, pM = Q_primitive_part(M, &cM); + if (!cM) return ZM_inv(pM,dM); + return gerepileupto(av, ZM_inv(pM, gdiv(dM,cM))); } /* x a matrix with integer coefficients. Return a multiple of the determinant @@ -1397,8 +1637,9 @@ QM_inv(GEN M, GEN dM) GEN detint(GEN x) { + gpmem_t av = avma, av1, lim; GEN pass,c,v,det1,piv,pivprec,vi,p1; - long i,j,k,rg,n,m,m1,av=avma,av1,lim; + long i,j,k,rg,n,m,m1; if (typ(x)!=t_MAT) err(typeer,"detint"); n=lg(x)-1; if (!n) return gun; @@ -1462,10 +1703,11 @@ detint(GEN x) } static void -gerepile_gauss_ker(GEN x, long m, long n, long k, long t, ulong av) +gerepile_gauss_ker(GEN x, long k, long t, gpmem_t av) { - ulong tetpil = avma,A; - long dec,u,i; + gpmem_t tetpil = avma, A; + long u,i, n = lg(x)-1, m = n? lg(x[1])-1: 0; + size_t dec; if (DEBUGMEM > 1) err(warnmem,"gauss_pivot_ker. k=%ld, n=%ld",k,n); for (u=t+1; u<=m; u++) copyifstack(coeff(x,u,k), coeff(x,u,k)); @@ -1475,22 +1717,23 @@ gerepile_gauss_ker(GEN x, long m, long n, long k, long (void)gerepile(av,tetpil,NULL); dec = av-tetpil; for (u=t+1; u<=m; u++) { - A=coeff(x,u,k); + A=(gpmem_t)coeff(x,u,k); if (A=bot) coeff(x,u,k)+=dec; } for (i=k+1; i<=n; i++) for (u=1; u<=m; u++) { - A=coeff(x,u,i); + A=(gpmem_t)coeff(x,u,i); if (A=bot) coeff(x,u,i)+=dec; } } static void -gerepile_gauss_FpM_ker(GEN x, GEN p, long m, long n, long k, long t, ulong av) +gerepile_gauss_FpM_ker(GEN x, GEN p, long k, long t, gpmem_t av) { - ulong tetpil = avma,A; - long dec,u,i; + gpmem_t tetpil = avma, A; + long u,i, n = lg(x)-1, m = n? lg(x[1])-1: 0; + size_t dec; if (DEBUGMEM > 1) err(warnmem,"gauss_pivot_ker. k=%ld, n=%ld",k,n); for (u=t+1; u<=m; u++) @@ -1502,13 +1745,13 @@ gerepile_gauss_FpM_ker(GEN x, GEN p, long m, long n, l (void)gerepile(av,tetpil,NULL); dec = av-tetpil; for (u=t+1; u<=m; u++) { - A=coeff(x,u,k); + A=(gpmem_t)coeff(x,u,k); if (A=bot) coeff(x,u,k)+=dec; } for (i=k+1; i<=n; i++) for (u=1; u<=m; u++) { - A=coeff(x,u,i); + A=(gpmem_t)coeff(x,u,i); if (A=bot) coeff(x,u,i)+=dec; } } @@ -1516,10 +1759,11 @@ gerepile_gauss_FpM_ker(GEN x, GEN p, long m, long n, l /* special gerepile for huge matrices */ static void -gerepile_gauss(GEN x,long m,long n,long k,long t,ulong av, long j, GEN c) +gerepile_gauss(GEN x,long k,long t,gpmem_t av, long j, GEN c) { - ulong tetpil = avma,A; - long dec,u,i; + gpmem_t tetpil = avma, A; + long u,i, n = lg(x)-1, m = n? lg(x[1])-1: 0; + size_t dec; if (DEBUGMEM > 1) err(warnmem,"gauss_pivot. k=%ld, n=%ld",k,n); for (u=t+1; u<=m; u++) @@ -1532,14 +1776,14 @@ gerepile_gauss(GEN x,long m,long n,long k,long t,ulong for (u=t+1; u<=m; u++) if (u==j || !c[u]) { - A=coeff(x,u,k); + A=(gpmem_t)coeff(x,u,k); if (A=bot) coeff(x,u,k)+=dec; } for (u=1; u<=m; u++) if (u==j || !c[u]) for (i=k+1; i<=n; i++) { - A=coeff(x,u,i); + A=(gpmem_t)coeff(x,u,i); if (A=bot) coeff(x,u,i)+=dec; } } @@ -1550,13 +1794,65 @@ gerepile_gauss(GEN x,long m,long n,long k,long t,ulong /* return n - rk(x) linearly independant vectors */ /* */ /*******************************************************************/ +static long +gauss_get_pivot_NZ(GEN x, GEN x0/*unused*/, GEN c, long i0) +{ + long i,lx = lg(x); + (void)x0; + if (c) + for (i=i0; i bit_accuracy(lg(x))); +} + +static long +gauss_get_pivot_max(GEN x, GEN x0, GEN c, long i0) +{ + long i,e, k = i0, ex = - (long)HIGHEXPOBIT, lx = lg(x); + GEN p,r; + if (c) + for (i=i0; i ex) { ex=e; k=i; } + } + else + for (i=i0; i ex) { ex=e; k=i; } + } + p = (GEN)x[k]; + r = (GEN)x0[k]; if (isexactzero(r)) r = x0; + return approx_0(p, r)? i: k; +} + /* x has INTEGER coefficients */ GEN keri(GEN x) { - GEN c,d,y,p,pp; - long i,j,k,r,t,n,m,av,av0,tetpil,lim; + gpmem_t av, av0, tetpil, lim; + GEN c,l,y,p,pp; + long i,j,k,r,t,n,m; if (typ(x)!=t_MAT) err(typeer,"keri"); n=lg(x)-1; if (!n) return cgetg(1,t_MAT); @@ -1564,42 +1860,40 @@ keri(GEN x) av0=avma; m=lg(x[1])-1; r=0; pp=cgetg(n+1,t_COL); x=dummycopy(x); p=gun; - c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0; - d=new_chunk(n+1); av=avma; lim=stack_lim(av,1); + c=cgetg(m+1, t_VECSMALL); for (k=1; k<=m; k++) c[k]=0; + l=cgetg(n+1, t_VECSMALL); + av = avma; lim = stack_lim(av,1); for (k=1; k<=n; k++) { - j=1; - while (j<=m && (c[j] || !signe(gcoeff(x,j,k))) ) j++; - if (j>m) + j = 1; + while ( j <= m && (c[j] || !signe(gcoeff(x,j,k))) ) j++; + if (j > m) { - r++; d[k]=0; + r++; l[k] = 0; for(j=1; j nl) break; + + d[k] = ck[i]; + c[i] = k; l[k] = i; /* pivot d[k] in x[i,k] */ } - if (k>nc) + if (k > nc) { avma = av; return zerocol(nc); } + if (k == 1) { avma = av; return gscalcol_i(gun,nc); } + y = cgetg(nc+1,t_COL); + y[1] = (long)ck[ l[1] ]; + for (q=d[1],j=2; j1)? coeff(x,l[1],k): un; - for (q=gun,j=2; j1) y[k]=lneg(gmul(q,(GEN) d[k-1])); - for (j=k+1; j<=nc; j++) y[j]=zero; - d=content(y); t=avma; - return gerepile(av,t,gdiv(y,d)); + y[j] = lneg(q); + for (j++; j<=nc; j++) y[j] = zero; + return gerepileupto(av, gdiv(y,content(y))); } /*******************************************************************/ @@ -1681,78 +1978,6 @@ deplin(GEN x) /* (kernel, image, complementary image, rank) */ /* */ /*******************************************************************/ -static long gauss_ex; -static int (*gauss_is_zero)(GEN); - -static int -real0(GEN x) -{ - return gcmp0(x) || (gexpo(x) < gauss_ex); -} - -static void -gauss_get_prec(GEN x, long prec) -{ - long pr = matprec(x); - - if (!pr) { gauss_is_zero = &gcmp0; return; } - if (pr > prec) prec = pr; - - gauss_ex = - (long)(0.85 * bit_accuracy(prec)); - gauss_is_zero = &real0; -} - -static long -gauss_get_pivot_NZ(GEN x, GEN x0/* unused */, GEN c, long i0) -{ - long i,lx = lg(x); - if (c) - for (i=i0; i bit_accuracy(lg(x))); -} - -static long -gauss_get_pivot_max(GEN x, GEN x0, GEN c, long i0) -{ - long i,e, k = i0, ex = -HIGHEXPOBIT, lx = lg(x); - GEN p,r; - if (c) - for (i=i0; i ex) { ex=e; k=i; } - } - else - for (i=i0; i ex) { ex=e; k=i; } - } - p = (GEN)x[k]; - r = (GEN)x0[k]; if (isexactzero(r)) r = x0; - return approx_0(p, r)? i: k; -} - /* return the transform of x under a standard Gauss pivot. r = dim ker(x). * d[k] contains the index of the first non-zero pivot in column k * If a != NULL, use x - a Id instead (for eigen) @@ -1761,7 +1986,8 @@ static GEN gauss_pivot_ker(GEN x0, GEN a, GEN *dd, long *rr) { GEN x,c,d,p,mun; - long i,j,k,r,t,n,m,av,lim; + gpmem_t av, lim; + long i,j,k,r,t,n,m; long (*get_pivot)(GEN,GEN,GEN,long); if (typ(x0)!=t_MAT) err(typeer,"gauss_pivot"); @@ -1801,7 +2027,7 @@ gauss_pivot_ker(GEN x0, GEN a, GEN *dd, long *rr) for (i=k+1; i<=n; i++) coeff(x,t,i) = ladd(gcoeff(x,t,i),gmul(p,gcoeff(x,j,i))); if (low_stack(lim, stack_lim(av,1))) - gerepile_gauss_ker(x,m,n,k,t,av); + gerepile_gauss_ker(x,k,t,av); } } } @@ -1815,7 +2041,8 @@ static void gauss_pivot(GEN x0, GEN *dd, long *rr) { GEN x,c,d,d0,mun,p; - long i,j,k,r,t,n,m,av,lim; + long i, j, k, r, t, n, m; + gpmem_t av, lim; long (*get_pivot)(GEN,GEN,GEN,long); if (typ(x0)!=t_MAT) err(typeer,"gauss_pivot"); @@ -1856,7 +2083,7 @@ gauss_pivot(GEN x0, GEN *dd, long *rr) for (i=k+1; i<=n; i++) coeff(x,t,i) = ladd(gcoeff(x,t,i), gmul(p,gcoeff(x,j,i))); if (low_stack(lim, stack_lim(av,1))) - gerepile_gauss(x,m,n,k,t,av,j,c); + gerepile_gauss(x,k,t,av,j,c); } for (i=k; i<=n; i++) coeff(x,j,i) = zero; /* dummy */ } @@ -1868,8 +2095,9 @@ gauss_pivot(GEN x0, GEN *dd, long *rr) static GEN ker0(GEN x, GEN a) { + gpmem_t av = avma, tetpil; GEN d,y; - long i,j,k,r,n, av = avma, tetpil; + long i,j,k,r,n; x = gauss_pivot_ker(x,a,&d,&r); if (!r) { avma=av; return cgetg(1,t_MAT); } @@ -1907,8 +2135,9 @@ matker0(GEN x,long flag) GEN image(GEN x) { + gpmem_t av = avma; GEN d,y; - long j,k,r, av = avma; + long j,k,r; gauss_pivot(x,&d,&r); @@ -1926,8 +2155,9 @@ image(GEN x) GEN imagecompl(GEN x) { + gpmem_t av = avma; GEN d,y; - long j,i,r,av = avma; + long j,i,r; gauss_pivot(x,&d,&r); avma=av; y=cgetg(r+1,t_VEC); @@ -1940,10 +2170,11 @@ imagecompl(GEN x) GEN imagecomplspec(GEN x, long *nlze) { + gpmem_t av = avma; GEN d,y; - long i,j,k,l,r,av = avma; + long i,j,k,l,r; - x = gtrans(x); l = lg(x); + x = gtrans_i(x); l = lg(x); gauss_pivot(x,&d,&r); avma=av; y = cgetg(l,t_VECSMALL); for (i=j=1, k=r+1; in) err(suppler2); - if (lx == n) return gcopy(x); +/* NB: d freed */ +static GEN +get_suppl(GEN x, GEN d, long r) +{ + gpmem_t av; + GEN y,c; + long j,k,rx,n; - zone = switch_stack(NULL, n*n); - switch_stack(zone,1); - y = myid? dummycopy(myid): idmat(n-1); - switch_stack(zone,0); - gauss_get_prec(x,0); - for (i=1; i=n) err(suppler2); - p1=(GEN)y[i]; y[i]=x[i]; if (i!=j) y[j]=(long)p1; - } - avma = av; y = gcopy(y); - free(zone); return y; + rx = lg(x)-1; + if (!rx) err(talker,"empty matrix in suppl"); + n = lg(x[1])-1; + if (rx == n && r == 0) { free(d); return gcopy(x); } + y = cgetg(n+1, t_MAT); + av = avma; + c = vecsmall_const(n,0); + k = 1; + /* c = lines containing pivots (could get it from gauss_pivot, but cheap) + * In theory r = 0 and d[j] > 0 for all j, but why take chances? */ + for (j=1; j<=rx; j++) + if (d[j]) + { + c[ d[j] ] = 1; + y[k++] = x[j]; + } + for (j=1; j<=n; j++) + if (!c[j]) y[k++] = j; + avma = av; + + rx -= r; + for (j=1; j<=rx; j++) + y[j] = lcopy((GEN)y[j]); + for ( ; j<=n; j++) + y[j] = (long)_ei(n, y[j]); + free(d); return y; } +/* x is an n x k matrix, rank(x) = k <= n. Return an invertible n x n matrix + * whose first k columns are given by x. If rank(x) < k, undefined result. */ GEN suppl(GEN x) { - return suppl_intern(x,NULL); + gpmem_t av = avma; + GEN d; + long r; + + gauss_pivot(x,&d,&r); + avma = av; return get_suppl(x,d,r); } +static void FpM_gauss_pivot(GEN x, GEN p, GEN *dd, long *rr); +static void FqM_gauss_pivot(GEN x, GEN T, GEN p, GEN *dd, long *rr); + GEN +FpM_suppl(GEN x, GEN p) +{ + gpmem_t av = avma; + GEN d; + long r; + + FpM_gauss_pivot(x,p, &d,&r); + avma = av; return get_suppl(x,d,r); +} +GEN +FqM_suppl(GEN x, GEN T, GEN p) +{ + gpmem_t av = avma; + GEN d; + long r; + + if (!T) return FpM_suppl(x,p); + FqM_gauss_pivot(x,T,p, &d,&r); + avma = av; return get_suppl(x,d,r); +} + +GEN image2(GEN x) { - long av=avma,tetpil,k,n,i; + gpmem_t av=avma,tetpil; + long k,n,i; GEN p1,p2; if (typ(x)!=t_MAT) err(typeer,"image2"); @@ -2068,7 +2346,8 @@ matimage0(GEN x,long flag) long rank(GEN x) { - long av = avma, r; + gpmem_t av = avma; + long r; GEN d; gauss_pivot(x,&d,&r); @@ -2078,13 +2357,12 @@ rank(GEN x) return lg(x)-1 - r; } -static void FpM_gauss_pivot(GEN x, GEN p, GEN *dd, long *rr); - /* if p != NULL, assume x integral and compute rank over Fp */ static GEN indexrank0(GEN x, GEN p, int small) { - long av = avma, i,j,n,r; + gpmem_t av = avma; + long i,j,n,r; GEN res,d,p1,p2; /* yield r = dim ker(x) */ @@ -2108,13 +2386,13 @@ indexrank0(GEN x, GEN p, int small) return res; } -GEN +GEN indexrank(GEN x) { return indexrank0(x,NULL,0); } -GEN +GEN sindexrank(GEN x) { return indexrank0(x,NULL,1); } -GEN +GEN FpM_sindexrank(GEN x, GEN p) { return indexrank0(x,p,1); } /*******************************************************************/ @@ -2122,7 +2400,7 @@ FpM_sindexrank(GEN x, GEN p) { return indexrank0(x,p,1 /* LINEAR ALGEBRA MODULO P */ /* */ /*******************************************************************/ - + /*If p is NULL no reduction is performed.*/ GEN FpM_mul(GEN x, GEN y, GEN p) @@ -2130,7 +2408,7 @@ FpM_mul(GEN x, GEN y, GEN p) long i,j,l,lx=lg(x), ly=lg(y); GEN z; if (ly==1) return cgetg(ly,t_MAT); - if (lx != lg(y[1])) err(operi,"* [mod p]",t_MAT,t_MAT); + if (lx != lg(y[1])) err(operi,"* [mod p]",x,y); z=cgetg(ly,t_MAT); if (lx==1) { @@ -2142,8 +2420,8 @@ FpM_mul(GEN x, GEN y, GEN p) { z[j] = lgetg(l,t_COL); for (i=1; im) + if (j > m) { - if (nontriv) { avma=av0; return NULL; } - r++; d[k]=0; + if (deplin) { + c = cgetg(n+1, t_VECSMALL); + for (i=1; im) { - if (nontriv) { avma = av0; return NULL; } + if (deplin) { + c = cgetg(n+1, t_COL); + for (i=1; im) { r++; d[k]=0; } + else + { + c[j]=k; d[k]=j; piv = gneg(FpXQ_inv(gcoeff(x,j,k), T,p)); + for (i=k+1; i<=n; i++) + coeff(x,j,i) = (long)Fq_mul(piv,gcoeff(x,j,i), T, p); + for (t=1; t<=m; t++) + if (!c[t]) /* no pivot on that line yet */ + { + piv=gcoeff(x,t,k); + if (signe(piv)) + { + coeff(x,t,k)=zero; + for (i=k+1; i<=n; i++) + coeff(x,t,i) = ladd(gcoeff(x,t,i), gmul(piv,gcoeff(x,j,i))); + if (low_stack(lim, stack_lim(av,1))) + gerepile_gauss(x,k,t,av,j,c); + } + } + for (i=k; i<=n; i++) coeff(x,j,i) = zero; /* dummy */ + } + } + *dd=d; *rr=r; } GEN FpM_image(GEN x, GEN p) { + gpmem_t av = avma; GEN d,y; - long j,k,r, av = avma; + long j,k,r; FpM_gauss_pivot(x,p,&d,&r); @@ -2406,7 +2761,8 @@ FpM_image(GEN x, GEN p) static GEN sFpM_invimage(GEN mat, GEN y, GEN p) { - long av=avma,i, nbcol = lg(mat); + gpmem_t av=avma; + long i, nbcol = lg(mat); GEN col,p1 = cgetg(nbcol+1,t_MAT),res; if (nbcol==1) return NULL; @@ -2434,7 +2790,8 @@ sFpM_invimage(GEN mat, GEN y, GEN p) GEN FpM_invimage(GEN m, GEN v, GEN p) { - long av=avma,j,lv,tv=typ(v); + gpmem_t av=avma; + long j,lv,tv=typ(v); GEN y,p1; if (typ(m)!=t_MAT) err(typeer,"inverseimage"); @@ -2455,7 +2812,7 @@ FpM_invimage(GEN m, GEN v, GEN p) } return y; } -/************************************************************** +/************************************************************** Rather stupid implementation of linear algebra in finite fields **************************************************************/ static GEN @@ -2467,7 +2824,7 @@ Fq_add(GEN x, GEN y, GEN T/*unused*/, GEN p) case 1: return FpX_Fp_add(x,y,p); case 2: return FpX_Fp_add(y,x,p); case 3: return FpX_add(x,y,p); - } + } return NULL; } #if 0 @@ -2475,7 +2832,7 @@ Fq_add(GEN x, GEN y, GEN T/*unused*/, GEN p) * For consistency we write the code there. * To avoid having to remove static status, we rewrite it in polarit1.c */ -static GEN +static GEN Fq_neg(GEN x, GEN T, GEN p) { switch(typ(x)==t_POL) @@ -2487,7 +2844,7 @@ Fq_neg(GEN x, GEN T, GEN p) } #endif -static GEN +GEN Fq_mul(GEN x, GEN y, GEN T, GEN p) { switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1)) @@ -2514,7 +2871,7 @@ Fq_neg_inv(GEN x, GEN T, GEN p) static GEN Fq_res(GEN x, GEN T, GEN p) { - ulong ltop=avma; + gpmem_t ltop=avma; switch(typ(x)==t_POL) { case 0: return modii(x,p); @@ -2524,11 +2881,11 @@ Fq_res(GEN x, GEN T, GEN p) } static void -Fq_gerepile_gauss_ker(GEN x, GEN T, GEN p, long m, long n, long k, long t, - ulong av) +Fq_gerepile_gauss_ker(GEN x, GEN T, GEN p, long k, long t, gpmem_t av) { - ulong tetpil = avma,A; - long dec,u,i; + gpmem_t tetpil = avma,A; + long u,i, n = lg(x)-1, m = n? lg(x[1])-1: 0; + size_t dec; if (DEBUGMEM > 1) err(warnmem,"gauss_pivot_ker. k=%ld, n=%ld",k,n); for (u=t+1; u<=m; u++) @@ -2540,23 +2897,26 @@ Fq_gerepile_gauss_ker(GEN x, GEN T, GEN p, long m, lon (void)gerepile(av,tetpil,NULL); dec = av-tetpil; for (u=t+1; u<=m; u++) { - A=coeff(x,u,k); + A=(gpmem_t)coeff(x,u,k); if (A=bot) coeff(x,u,k)+=dec; } for (i=k+1; i<=n; i++) for (u=1; u<=m; u++) { - A=coeff(x,u,i); + A=(gpmem_t)coeff(x,u,i); if (A=bot) coeff(x,u,i)+=dec; } } static GEN -FqM_ker_i(GEN x, GEN T, GEN p, long nontriv) +FqM_ker_i(GEN x, GEN T, GEN p, long deplin) { + gpmem_t av0,av,lim,tetpil; GEN y,c,d,piv,mun; - long i,j,k,r,t,n,m,av0,av,lim,tetpil; + long i,j,k,r,t,n,m; + if (!T) return FpM_ker_i(x,p,deplin); + if (typ(x)!=t_MAT) err(typeer,"FpM_ker"); n=lg(x)-1; if (!n) return cgetg(1,t_MAT); @@ -2575,7 +2935,12 @@ FqM_ker_i(GEN x, GEN T, GEN p, long nontriv) } if (j>m) { - if (nontriv) { avma = av0; return NULL; } + if (deplin) { + c = cgetg(n+1, t_COL); + for (i=1; i 7) timer2(); + if (DEBUGLEVEL > 7) (void)timer2(); - av = avma; lim = stack_lim(av,2); + av = avma; lim = stack_lim(av,2); a = dummycopy(a); s = 1; for (pprec=gun,i=1; i