Annotation of OpenXM_contrib/pari-2.2/src/basemath/alglin1.c, Revision 1.1
1.1 ! noro 1: /* $Id: alglin1.c,v 1.57 2001/10/01 12:11:27 karim Exp $
! 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:
! 16: /********************************************************************/
! 17: /** **/
! 18: /** LINEAR ALGEBRA **/
! 19: /** (first part) **/
! 20: /** **/
! 21: /********************************************************************/
! 22: #include "pari.h"
! 23:
! 24: /* for linear algebra mod p */
! 25: #ifdef LONG_IS_64BIT
! 26: # define MASK (0x7fffffff00000000UL)
! 27: #else
! 28: # define MASK (0x7fff0000UL)
! 29: #endif
! 30:
! 31: /* 2p^2 < 2^BITS_IN_LONG */
! 32: #ifdef LONG_IS_64BIT
! 33: # define u_OK_ULONG(p) ((ulong)p <= 3037000493UL)
! 34: #else
! 35: # define u_OK_ULONG(p) ((ulong)p <= 46337UL)
! 36: #endif
! 37: #define OK_ULONG(p) (lgefint(p) == 3 && u_OK_ULONG(p[2]))
! 38: extern GEN ZM_init_CRT(GEN Hp, ulong p);
! 39: extern ulong u_invmod(ulong x, ulong p);
! 40:
! 41: /*******************************************************************/
! 42: /* */
! 43: /* TRANSPOSE */
! 44: /* */
! 45: /*******************************************************************/
! 46:
! 47: /* No copy*/
! 48: GEN
! 49: gtrans_i(GEN x)
! 50: {
! 51: long i,j,lx,dx, tx=typ(x);
! 52: GEN y,p1;
! 53:
! 54: if (! is_matvec_t(tx)) err(typeer,"gtrans_i");
! 55: switch(tx)
! 56: {
! 57: case t_VEC:
! 58: y=dummycopy(x); settyp(y,t_COL); break;
! 59:
! 60: case t_COL:
! 61: y=dummycopy(x); settyp(y,t_VEC); break;
! 62:
! 63: case t_MAT:
! 64: lx=lg(x); if (lx==1) return cgetg(1,t_MAT);
! 65: dx=lg(x[1]); y=cgetg(dx,tx);
! 66: for (i=1; i<dx; i++)
! 67: {
! 68: p1=cgetg(lx,t_COL); y[i]=(long)p1;
! 69: for (j=1; j<lx; j++) p1[j]=coeff(x,i,j);
! 70: }
! 71: break;
! 72:
! 73: default: y=x; break;
! 74: }
! 75: return y;
! 76: }
! 77:
! 78:
! 79: GEN
! 80: gtrans(GEN x)
! 81: {
! 82: long i,j,lx,dx, tx=typ(x);
! 83: GEN y,p1;
! 84:
! 85: if (! is_matvec_t(tx)) err(typeer,"gtrans");
! 86: switch(tx)
! 87: {
! 88: case t_VEC:
! 89: y=gcopy(x); settyp(y,t_COL); break;
! 90:
! 91: case t_COL:
! 92: y=gcopy(x); settyp(y,t_VEC); break;
! 93:
! 94: case t_MAT:
! 95: lx=lg(x); if (lx==1) return cgetg(1,t_MAT);
! 96: dx=lg(x[1]); y=cgetg(dx,tx);
! 97: for (i=1; i<dx; i++)
! 98: {
! 99: p1=cgetg(lx,t_COL); y[i]=(long)p1;
! 100: for (j=1; j<lx; j++) p1[j]=lcopy(gcoeff(x,i,j));
! 101: }
! 102: break;
! 103:
! 104: default: y=gcopy(x); break;
! 105: }
! 106: return y;
! 107: }
! 108:
! 109: /*******************************************************************/
! 110: /* */
! 111: /* CONCATENATION & EXTRACTION */
! 112: /* */
! 113: /*******************************************************************/
! 114:
! 115: static GEN
! 116: strconcat(GEN x, GEN y)
! 117: {
! 118: long flx=0,fly=0,l;
! 119: char *sx,*sy,*str;
! 120:
! 121: if (typ(x)==t_STR) sx = GSTR(x); else { flx=1; sx = GENtostr(x); }
! 122: if (typ(y)==t_STR) sy = GSTR(y); else { fly=1; sy = GENtostr(y); }
! 123: l = strlen(sx) + strlen(sy) + 1;
! 124: l = (l+BYTES_IN_LONG) >> TWOPOTBYTES_IN_LONG;
! 125: x = cgetg(l+1, t_STR); str = GSTR(x);
! 126: strcpy(str,sx);
! 127: strcat(str,sy);
! 128: if (flx) free(sx);
! 129: if (fly) free(sy);
! 130: return x;
! 131: }
! 132:
! 133: /* concat 3 matrices. Internal */
! 134: GEN
! 135: concatsp3(GEN x, GEN y, GEN z)
! 136: {
! 137: long i, lx = lg(x), ly = lg(y), lz = lg(z);
! 138: GEN r = cgetg(lx+ly+lz-2, t_MAT), t = r;
! 139: for (i=1; i<lx; i++) *++t = *++x;
! 140: for (i=1; i<ly; i++) *++t = *++y;
! 141: for (i=1; i<lz; i++) *++t = *++z;
! 142: return r;
! 143: }
! 144:
! 145: /* concat A and B vertically. Internal */
! 146: GEN
! 147: vconcat(GEN A, GEN B)
! 148: {
! 149: long la,ha,hb,hc,i,j;
! 150: GEN M,a,b,c;
! 151:
! 152: la = lg(A); if (la==1) return A;
! 153: ha = lg(A[1]); M = cgetg(la,t_MAT);
! 154: hb = lg(B[1]); hc = ha+hb-1;
! 155: for (j=1; j<la; j++)
! 156: {
! 157: c = cgetg(hc,t_COL); M[j] = (long)c; a = (GEN)A[j]; b = (GEN)B[j];
! 158: for (i=1; i<ha; i++) *++c = *++a;
! 159: for (i=1; i<hb; i++) *++c = *++b;
! 160: }
! 161: return M;
! 162: }
! 163:
! 164: GEN
! 165: _vec(GEN x) { GEN v = cgetg(2, t_VEC); v[1] = (long)x; return v; }
! 166: GEN
! 167: _col(GEN x) { GEN v = cgetg(2, t_COL); v[1] = (long)x; return v; }
! 168:
! 169: GEN
! 170: concatsp(GEN x, GEN y)
! 171: {
! 172: long tx=typ(x),ty=typ(y),lx=lg(x),ly=lg(y),i;
! 173: GEN z,p1;
! 174:
! 175: if (tx==t_LIST || ty==t_LIST) return listconcat(x,y);
! 176: if (tx==t_STR || ty==t_STR) return strconcat(x,y);
! 177:
! 178: if (tx==t_MAT && lx==1)
! 179: {
! 180: if (ty!=t_VEC || ly==1) return gtomat(y);
! 181: err(concater);
! 182: }
! 183: if (ty==t_MAT && ly==1)
! 184: {
! 185: if (tx!=t_VEC || lx==1) return gtomat(x);
! 186: err(concater);
! 187: }
! 188:
! 189: if (! is_matvec_t(tx))
! 190: {
! 191: if (! is_matvec_t(ty))
! 192: {
! 193: z=cgetg(3,t_VEC); z[1]=(long)x; z[2]=(long)y;
! 194: return z;
! 195: }
! 196: z=cgetg(ly+1,ty);
! 197: if (ty != t_MAT) p1 = x;
! 198: else
! 199: {
! 200: if (lg(y[1])!=2) err(concater);
! 201: p1=cgetg(2,t_COL); p1[1]=(long)x;
! 202: }
! 203: for (i=2; i<=ly; i++) z[i]=y[i-1];
! 204: z[1]=(long)p1; return z;
! 205: }
! 206: if (! is_matvec_t(ty))
! 207: {
! 208: z=cgetg(lx+1,tx);
! 209: if (tx != t_MAT) p1 = y;
! 210: else
! 211: {
! 212: if (lg(x[1])!=2) err(concater);
! 213: p1=cgetg(2,t_COL); p1[1]=(long)y;
! 214: }
! 215: for (i=1; i<lx; i++) z[i]=x[i];
! 216: z[lx]=(long)p1; return z;
! 217: }
! 218:
! 219: if (tx == ty)
! 220: {
! 221: if (tx == t_MAT && lg(x[1]) != lg(y[1])) err(concater);
! 222: z=cgetg(lx+ly-1,tx);
! 223: for (i=1; i<lx; i++) z[i]=x[i];
! 224: for (i=1; i<ly; i++) z[lx+i-1]=y[i];
! 225: return z;
! 226: }
! 227:
! 228: switch(tx)
! 229: {
! 230: case t_VEC:
! 231: switch(ty)
! 232: {
! 233: case t_COL:
! 234: if (lx<=2) return (lx==1)? y: concatsp((GEN) x[1],y);
! 235: if (ly>=3) break;
! 236: return (ly==1)? x: concatsp(x,(GEN) y[1]);
! 237: case t_MAT:
! 238: z=cgetg(ly,ty); if (lx != ly) break;
! 239: for (i=1; i<ly; i++) z[i]=(long)concatsp((GEN) x[i],(GEN) y[i]);
! 240: return z;
! 241: }
! 242: break;
! 243:
! 244: case t_COL:
! 245: switch(ty)
! 246: {
! 247: case t_VEC:
! 248: if (lx<=2) return (lx==1)? y: concatsp((GEN) x[1],y);
! 249: if (ly>=3) break;
! 250: return (ly==1)? x: concatsp(x,(GEN) y[1]);
! 251: case t_MAT:
! 252: if (lx != lg(y[1])) break;
! 253: z=cgetg(ly+1,ty); z[1]=(long)x;
! 254: for (i=2; i<=ly; i++) z[i]=y[i-1];
! 255: return z;
! 256: }
! 257: break;
! 258:
! 259: case t_MAT:
! 260: switch(ty)
! 261: {
! 262: case t_VEC:
! 263: z=cgetg(lx,tx); if (ly != lx) break;
! 264: for (i=1; i<lx; i++) z[i]=(long)concatsp((GEN) x[i],(GEN) y[i]);
! 265: return z;
! 266: case t_COL:
! 267: if (ly != lg(x[1])) break;
! 268: z=cgetg(lx+1,tx); z[lx]=(long)y;
! 269: for (i=1; i<lx; i++) z[i]=x[i];
! 270: return z;
! 271: }
! 272: break;
! 273: }
! 274: err(concater);
! 275: return NULL; /* not reached */
! 276: }
! 277:
! 278: GEN
! 279: concat(GEN x, GEN y)
! 280: {
! 281: long tx = typ(x), lx,ty,ly,i;
! 282: GEN z,p1;
! 283:
! 284: if (!y)
! 285: {
! 286: ulong av = avma;
! 287: if (tx == t_LIST)
! 288: { lx = lgef(x); i = 2; }
! 289: else if (tx == t_VEC)
! 290: { lx = lg(x); i = 1; }
! 291: else
! 292: {
! 293: err(concater);
! 294: return NULL; /* not reached */
! 295: }
! 296: if (i>=lx) err(talker,"trying to concat elements of an empty vector");
! 297: y = (GEN)x[i++];
! 298: for (; i<lx; i++) y = concatsp(y, (GEN)x[i]);
! 299: return gerepilecopy(av,y);
! 300: }
! 301: ty = typ(y);
! 302: if (tx==t_LIST || ty==t_LIST) return listconcat(x,y);
! 303: if (tx==t_STR || ty==t_STR) return strconcat(x,y);
! 304: lx=lg(x); ly=lg(y);
! 305:
! 306: if (tx==t_MAT && lx==1)
! 307: {
! 308: if (ty!=t_VEC || ly==1) return gtomat(y);
! 309: err(concater);
! 310: }
! 311: if (ty==t_MAT && ly==1)
! 312: {
! 313: if (tx!=t_VEC || lx==1) return gtomat(x);
! 314: err(concater);
! 315: }
! 316:
! 317: if (! is_matvec_t(tx))
! 318: {
! 319: if (! is_matvec_t(ty))
! 320: {
! 321: z=cgetg(3,t_VEC); z[1]=lcopy(x); z[2]=lcopy(y);
! 322: return z;
! 323: }
! 324: z=cgetg(ly+1,ty);
! 325: if (ty != t_MAT) p1 = gcopy(x);
! 326: else
! 327: {
! 328: if (lg(y[1])!=2) err(concater);
! 329: p1=cgetg(2,t_COL); p1[1]=lcopy(x);
! 330: }
! 331: for (i=2; i<=ly; i++) z[i]=lcopy((GEN) y[i-1]);
! 332: z[1]=(long)p1; return z;
! 333: }
! 334: if (! is_matvec_t(ty))
! 335: {
! 336: z=cgetg(lx+1,tx);
! 337: if (tx != t_MAT) p1 = gcopy(y);
! 338: else
! 339: {
! 340: if (lg(x[1])!=2) err(concater);
! 341: p1=cgetg(2,t_COL); p1[1]=lcopy(y);
! 342: }
! 343: for (i=1; i<lx; i++) z[i]=lcopy((GEN) x[i]);
! 344: z[lx]=(long)p1; return z;
! 345: }
! 346:
! 347: if (tx == ty)
! 348: {
! 349: if (tx == t_MAT && lg(x[1]) != lg(y[1])) err(concater);
! 350: z=cgetg(lx+ly-1,tx);
! 351: for (i=1; i<lx; i++) z[i]=lcopy((GEN) x[i]);
! 352: for (i=1; i<ly; i++) z[lx+i-1]=lcopy((GEN) y[i]);
! 353: return z;
! 354: }
! 355:
! 356: switch(tx)
! 357: {
! 358: case t_VEC:
! 359: switch(ty)
! 360: {
! 361: case t_COL:
! 362: if (lx<=2) return (lx==1)? gcopy(y): concat((GEN) x[1],y);
! 363: if (ly>=3) break;
! 364: return (ly==1)? gcopy(x): concat(x,(GEN) y[1]);
! 365: case t_MAT:
! 366: z=cgetg(ly,ty); if (lx != ly) break;
! 367: for (i=1; i<ly; i++) z[i]=lconcat((GEN) x[i],(GEN) y[i]);
! 368: return z;
! 369: }
! 370: break;
! 371:
! 372: case t_COL:
! 373: switch(ty)
! 374: {
! 375: case t_VEC:
! 376: if (lx<=2) return (lx==1)? gcopy(y): concat((GEN) x[1],y);
! 377: if (ly>=3) break;
! 378: return (ly==1)? gcopy(x): concat(x,(GEN) y[1]);
! 379: case t_MAT:
! 380: if (lx != lg(y[1])) break;
! 381: z=cgetg(ly+1,ty); z[1]=lcopy(x);
! 382: for (i=2; i<=ly; i++) z[i]=lcopy((GEN) y[i-1]);
! 383: return z;
! 384: }
! 385: break;
! 386:
! 387: case t_MAT:
! 388: switch(ty)
! 389: {
! 390: case t_VEC:
! 391: z=cgetg(lx,tx); if (ly != lx) break;
! 392: for (i=1; i<lx; i++) z[i]=lconcat((GEN) x[i],(GEN) y[i]);
! 393: return z;
! 394: case t_COL:
! 395: if (ly != lg(x[1])) break;
! 396: z=cgetg(lx+1,tx); z[lx]=lcopy(y);
! 397: for (i=1; i<lx; i++) z[i]=lcopy((GEN) x[i]);
! 398: return z;
! 399: }
! 400: break;
! 401: }
! 402: err(concater);
! 403: return NULL; /* not reached */
! 404: }
! 405:
! 406: static long
! 407: str_to_long(char *s, char **pt)
! 408: {
! 409: long a = atol(s);
! 410: while (isspace((int)*s)) s++;
! 411: if (*s == '-' || *s == '+') s++;
! 412: while (isdigit((int)*s) || isspace((int)*s)) s++;
! 413: *pt = s; return a;
! 414: }
! 415:
! 416: static int
! 417: get_range(char *s, long *a, long *b, long *cmpl, long lx)
! 418: {
! 419: long max = lx - 1;
! 420:
! 421: *a = 1; *b = max;
! 422: if (*s == '^') { *cmpl = 1; s++; } else *cmpl = 0;
! 423: if (*s == 0) return 0;
! 424: if (*s != '.')
! 425: {
! 426: *a = str_to_long(s, &s);
! 427: if (*a < 0) *a += lx;
! 428: if (*a<1 || *a>max) return 0;
! 429: }
! 430: if (*s == '.')
! 431: {
! 432: s++; if (*s != '.') return 0;
! 433: do s++; while (isspace((int)*s));
! 434: if (*s)
! 435: {
! 436: *b = str_to_long(s, &s);
! 437: if (*b < 0) *b += lx;
! 438: if (*b<1 || *b>max || *s) return 0;
! 439: }
! 440: return 1;
! 441: }
! 442: if (*s) return 0;
! 443: *b = *a; return 1;
! 444: }
! 445:
! 446: GEN
! 447: vecextract_i(GEN A, long y1, long y2)
! 448: {
! 449: long i,lB = y2 - y1 + 2;
! 450: GEN B = cgetg(lB, typ(A));
! 451: for (i=1; i<lB; i++) B[i] = A[y1-1+i];
! 452: return B;
! 453: }
! 454:
! 455: GEN
! 456: rowextract_i(GEN A, long x1, long x2)
! 457: {
! 458: long i, lB = lg(A);
! 459: GEN B = cgetg(lB, typ(A));
! 460: for (i=1; i<lB; i++) B[i] = (long)vecextract_i((GEN)A[i],x1,x2);
! 461: return B;
! 462: }
! 463:
! 464: GEN
! 465: vecextract_p(GEN A, GEN p)
! 466: {
! 467: long i,lB = lg(p);
! 468: GEN B = cgetg(lB, typ(A));
! 469: for (i=1; i<lB; i++) B[i] = A[p[i]];
! 470: return B;
! 471: }
! 472:
! 473: GEN
! 474: rowextract_p(GEN A, GEN p)
! 475: {
! 476: long i, lB = lg(A);
! 477: GEN B = cgetg(lB, typ(A));
! 478: for (i=1; i<lB; i++) B[i] = (long)vecextract_p((GEN)A[i],p);
! 479: return B;
! 480: }
! 481:
! 482: void
! 483: vecselect_p(GEN A, GEN B, GEN p, long init, long lB)
! 484: {
! 485: long i; setlg(B, lB);
! 486: for (i=init; i<lB; i++) B[i] = A[p[i]];
! 487: }
! 488:
! 489: /* B := p . A = row selection according to permutation p. Treat only lower
! 490: * right corner init x init */
! 491: void
! 492: rowselect_p(GEN A, GEN B, GEN p, long init)
! 493: {
! 494: long i, lB = lg(A), lp = lg(p);
! 495: for (i=1; i<init; i++) setlg(B[i],lp);
! 496: for ( ; i<lB; i++) vecselect_p((GEN)A[i],(GEN)B[i],p,init,lp);
! 497: }
! 498:
! 499: GEN
! 500: extract(GEN x, GEN l)
! 501: {
! 502: long av,i,j, tl = typ(l), tx = typ(x), lx = lg(x);
! 503: GEN y;
! 504:
! 505: if (! is_matvec_t(tx)) err(typeer,"extract");
! 506: if (tl==t_INT)
! 507: {
! 508: /* extract components of x as per the bits of mask l */
! 509: if (!signe(l)) return cgetg(1,tx);
! 510: av=avma; y = (GEN) gpmalloc(lx*sizeof(long));
! 511: i = j = 1; while (!mpodd(l)) { l=shifti(l,-1); i++; }
! 512: while (signe(l) && i<lx)
! 513: {
! 514: if (mod2(l)) y[j++] = x[i];
! 515: i++; l=shifti(l,-1);
! 516: }
! 517: if (signe(l)) err(talker,"mask too large in vecextract");
! 518: y[0] = evaltyp(tx) | evallg(j);
! 519: avma=av; x = gcopy(y); free(y); return x;
! 520: }
! 521: if (tl==t_STR)
! 522: {
! 523: char *s = GSTR(l);
! 524: long first, last, cmpl;
! 525: if (! get_range(s, &first, &last, &cmpl, lx))
! 526: err(talker, "incorrect range in extract");
! 527: if (lx == 1) return gcopy(x);
! 528: if (cmpl)
! 529: {
! 530: if (first <= last)
! 531: {
! 532: y = cgetg(lx - (last - first + 1),tx);
! 533: for (j=1; j<first; j++) y[j] = lcopy((GEN)x[j]);
! 534: for (i=last+1; i<lx; i++,j++) y[j] = lcopy((GEN)x[i]);
! 535: }
! 536: else
! 537: {
! 538: y = cgetg(lx - (first - last + 1),tx);
! 539: for (j=1,i=lx-1; i>first; i--,j++) y[j] = lcopy((GEN)x[i]);
! 540: for (i=last-1; i>0; i--,j++) y[j] = lcopy((GEN)x[i]);
! 541: }
! 542: }
! 543: else
! 544: {
! 545: if (first <= last)
! 546: {
! 547: y = cgetg(last-first+2,tx);
! 548: for (i=first,j=1; i<=last; i++,j++) y[j] = lcopy((GEN)x[i]);
! 549: }
! 550: else
! 551: {
! 552: y = cgetg(first-last+2,tx);
! 553: for (i=first,j=1; i>=last; i--,j++) y[j] = lcopy((GEN)x[i]);
! 554: }
! 555: }
! 556: return y;
! 557: }
! 558:
! 559: if (is_vec_t(tl))
! 560: {
! 561: long ll=lg(l); y=cgetg(ll,tx);
! 562: for (i=1; i<ll; i++)
! 563: {
! 564: j = itos((GEN) l[i]);
! 565: if (j>=lx || j<=0) err(talker,"no such component in vecextract");
! 566: y[i] = lcopy((GEN) x[j]);
! 567: }
! 568: return y;
! 569: }
! 570: if (tl == t_VECSMALL)
! 571: {
! 572: long ll=lg(l); y=cgetg(ll,tx);
! 573: for (i=1; i<ll; i++)
! 574: {
! 575: j = l[i];
! 576: if (j>=lx || j<=0) err(talker,"no such component in vecextract");
! 577: y[i] = lcopy((GEN) x[j]);
! 578: }
! 579: return y;
! 580: }
! 581: err(talker,"incorrect mask in vecextract");
! 582: return NULL; /* not reached */
! 583: }
! 584:
! 585: GEN
! 586: matextract(GEN x, GEN l1, GEN l2)
! 587: {
! 588: long av = avma, tetpil;
! 589:
! 590: if (typ(x)!=t_MAT) err(typeer,"matextract");
! 591: x = extract(gtrans(extract(x,l2)),l1); tetpil=avma;
! 592: return gerepile(av,tetpil, gtrans(x));
! 593: }
! 594:
! 595: GEN
! 596: extract0(GEN x, GEN l1, GEN l2)
! 597: {
! 598: if (! l2) return extract(x,l1);
! 599: return matextract(x,l1,l2);
! 600: }
! 601:
! 602: /*******************************************************************/
! 603: /* */
! 604: /* SCALAR-MATRIX OPERATIONS */
! 605: /* */
! 606: /*******************************************************************/
! 607:
! 608: /* create the square nxn matrix equal to z*Id */
! 609: static GEN
! 610: gscalmat_proto(GEN z, GEN myzero, long n, int flag)
! 611: {
! 612: long i,j;
! 613: GEN y = cgetg(n+1,t_MAT);
! 614: if (n < 0) err(talker,"negative size in scalmat");
! 615: if (flag) z = (flag==1)? stoi((long)z): gcopy(z);
! 616: for (i=1; i<=n; i++)
! 617: {
! 618: y[i]=lgetg(n+1,t_COL);
! 619: for (j=1; j<=n; j++)
! 620: coeff(y,j,i) = (i==j)? (long)z: (long)myzero;
! 621: }
! 622: return y;
! 623: }
! 624:
! 625: GEN
! 626: gscalmat(GEN x, long n) { return gscalmat_proto(x,gzero,n,2); }
! 627:
! 628: GEN
! 629: gscalsmat(long x, long n) { return gscalmat_proto((GEN)x,gzero,n,1); }
! 630:
! 631: GEN
! 632: idmat(long n) { return gscalmat_proto(gun,gzero,n,0); }
! 633:
! 634: GEN
! 635: idmat_intern(long n,GEN myun,GEN z) { return gscalmat_proto(myun,z,n,0); }
! 636:
! 637: GEN
! 638: gscalcol_proto(GEN z, GEN myzero, long n)
! 639: {
! 640: GEN y = cgetg(n+1,t_COL);
! 641: long i;
! 642:
! 643: if (n)
! 644: {
! 645: y[1]=(long)z;
! 646: for (i=2; i<=n; i++) y[i]=(long)myzero;
! 647: }
! 648: return y;
! 649: }
! 650:
! 651: GEN
! 652: zerocol(long n)
! 653: {
! 654: GEN y = cgetg(n+1,t_COL);
! 655: long i; for (i=1; i<=n; i++) y[i]=zero;
! 656: return y;
! 657: }
! 658:
! 659: GEN
! 660: zerovec(long n)
! 661: {
! 662: GEN y = cgetg(n+1,t_VEC);
! 663: long i; for (i=1; i<=n; i++) y[i]=zero;
! 664: return y;
! 665: }
! 666:
! 667: GEN
! 668: zeromat(long m, long n)
! 669: {
! 670: GEN y = cgetg(n+1,t_MAT);
! 671: long i; for (i=1; i<=n; i++) y[i]=(long)zerocol(m);
! 672: return y;
! 673: }
! 674:
! 675: GEN
! 676: gscalcol(GEN x, long n)
! 677: {
! 678: GEN y=gscalcol_proto(gzero,gzero,n);
! 679: if (n) y[1]=lcopy(x);
! 680: return y;
! 681: }
! 682:
! 683: GEN
! 684: gscalcol_i(GEN x, long n) { return gscalcol_proto(x,gzero,n); }
! 685:
! 686: GEN
! 687: gtomat(GEN x)
! 688: {
! 689: long tx,lx,i;
! 690: GEN y,p1;
! 691:
! 692: if (!x) return cgetg(1, t_MAT);
! 693: tx = typ(x);
! 694: if (! is_matvec_t(tx))
! 695: {
! 696: y=cgetg(2,t_MAT); p1=cgetg(2,t_COL); y[1]=(long)p1;
! 697: p1[1]=lcopy(x); return y;
! 698: }
! 699: switch(tx)
! 700: {
! 701: case t_VEC:
! 702: lx=lg(x); y=cgetg(lx,t_MAT);
! 703: for (i=1; i<lx; i++)
! 704: {
! 705: p1=cgetg(2,t_COL); y[i]=(long)p1;
! 706: p1[1]=lcopy((GEN) x[i]);
! 707: }
! 708: break;
! 709: case t_COL:
! 710: y=cgetg(2,t_MAT); y[1]=lcopy(x); break;
! 711: default: /* case t_MAT: */
! 712: y=gcopy(x); break;
! 713: }
! 714: return y;
! 715: }
! 716:
! 717: long
! 718: isscalarmat(GEN x, GEN s)
! 719: {
! 720: long nco,i,j;
! 721:
! 722: if (typ(x)!=t_MAT) err(typeer,"isdiagonal");
! 723: nco=lg(x)-1; if (!nco) return 1;
! 724: if (nco != lg(x[1])-1) return 0;
! 725:
! 726: for (j=1; j<=nco; j++)
! 727: {
! 728: GEN *col = (GEN*) x[j];
! 729: for (i=1; i<=nco; i++)
! 730: if (i == j)
! 731: {
! 732: if (!gegal(col[i],s)) return 0;
! 733: }
! 734: else if (!gcmp0(col[i])) return 0;
! 735: }
! 736: return 1;
! 737: }
! 738:
! 739: long
! 740: isdiagonal(GEN x)
! 741: {
! 742: long nco,i,j;
! 743:
! 744: if (typ(x)!=t_MAT) err(typeer,"isdiagonal");
! 745: nco=lg(x)-1; if (!nco) return 1;
! 746: if (nco != lg(x[1])-1) return 0;
! 747:
! 748: for (j=1; j<=nco; j++)
! 749: {
! 750: GEN *col = (GEN*) x[j];
! 751: for (i=1; i<=nco; i++)
! 752: if (i!=j && !gcmp0(col[i])) return 0;
! 753: }
! 754: return 1;
! 755: }
! 756:
! 757: /* create the diagonal matrix, whose diagonal is given by x */
! 758: GEN
! 759: diagonal(GEN x)
! 760: {
! 761: long i,j,lx,tx=typ(x);
! 762: GEN y,p1;
! 763:
! 764: if (! is_matvec_t(tx)) return gscalmat(x,1);
! 765: if (tx==t_MAT)
! 766: {
! 767: if (isdiagonal(x)) return gcopy(x);
! 768: err(talker,"incorrect object in diagonal");
! 769: }
! 770: lx=lg(x); y=cgetg(lx,t_MAT);
! 771: for (j=1; j<lx; j++)
! 772: {
! 773: p1=cgetg(lx,t_COL); y[j]=(long)p1;
! 774: for (i=1; i<lx; i++)
! 775: p1[i] = (i==j)? lcopy((GEN) x[i]): zero;
! 776: }
! 777: return y;
! 778: }
! 779:
! 780: /* compute m*diagonal(d) */
! 781: GEN
! 782: matmuldiagonal(GEN m, GEN d)
! 783: {
! 784: long j=typ(d),lx=lg(m);
! 785: GEN y;
! 786:
! 787: if (typ(m)!=t_MAT) err(typeer,"matmuldiagonal");
! 788: if (! is_vec_t(j) || lg(d)!=lx)
! 789: err(talker,"incorrect vector in matmuldiagonal");
! 790: y=cgetg(lx,t_MAT);
! 791: for (j=1; j<lx; j++) y[j] = lmul((GEN) d[j],(GEN) m[j]);
! 792: return y;
! 793: }
! 794:
! 795: /* compute m*n assuming the result is a diagonal matrix */
! 796: GEN
! 797: matmultodiagonal(GEN m, GEN n)
! 798: {
! 799: long lx,i,j;
! 800: GEN s,y;
! 801:
! 802: if (typ(m)!=t_MAT || typ(n)!=t_MAT) err(typeer,"matmultodiagonal");
! 803: lx=lg(n); y=idmat(lx-1);
! 804: if (lx == 1)
! 805: { if (lg(m) != 1) err(consister,"matmultodiagonal"); }
! 806: else
! 807: { if (lg(m) != lg(n[1])) err(consister,"matmultodiagonal"); }
! 808: for (i=1; i<lx; i++)
! 809: {
! 810: s = gzero;
! 811: for (j=1; j<lx; j++)
! 812: s = gadd(s,gmul(gcoeff(m,i,j),gcoeff(n,j,i)));
! 813: coeff(y,i,i) = (long)s;
! 814: }
! 815: return y;
! 816: }
! 817:
! 818: /* [m[1,1], ..., m[l,l]], internal */
! 819: GEN
! 820: mattodiagonal_i(GEN m)
! 821: {
! 822: long i, lx = lg(m);
! 823: GEN y = cgetg(lx,t_VEC);
! 824: for (i=1; i<lx; i++) y[i] = coeff(m,i,i);
! 825: return y;
! 826: }
! 827:
! 828: /* same, public function */
! 829: GEN
! 830: mattodiagonal(GEN m)
! 831: {
! 832: long i, lx = lg(m);
! 833: GEN y = cgetg(lx,t_VEC);
! 834:
! 835: if (typ(m)!=t_MAT) err(typeer,"mattodiagonal");
! 836: for (i=1; i<lx; i++) y[i] = lcopy(gcoeff(m,i,i));
! 837: return y;
! 838: }
! 839:
! 840: /*******************************************************************/
! 841: /* */
! 842: /* ADDITION SCALAR + MATRIX */
! 843: /* */
! 844: /*******************************************************************/
! 845:
! 846: /* create the square matrix x*Id + y */
! 847: GEN
! 848: gaddmat(GEN x, GEN y)
! 849: {
! 850: long ly,dy,i,j;
! 851: GEN z;
! 852:
! 853: ly=lg(y); if (ly==1) return cgetg(1,t_MAT);
! 854: dy=lg(y[1]);
! 855: if (typ(y)!=t_MAT || ly!=dy) err(mattype1,"gaddmat");
! 856: z=cgetg(ly,t_MAT);
! 857: for (i=1; i<ly; i++)
! 858: {
! 859: z[i]=lgetg(dy,t_COL);
! 860: for (j=1; j<dy; j++)
! 861: coeff(z,j,i) = i==j? ladd(x,gcoeff(y,j,i)): lcopy(gcoeff(y,j,i));
! 862: }
! 863: return z;
! 864: }
! 865:
! 866: /*******************************************************************/
! 867: /* */
! 868: /* Solve A*X=B (Gauss pivot) */
! 869: /* */
! 870: /*******************************************************************/
! 871: #define swap(x,y) { long _t=x; x=y; y=_t; }
! 872:
! 873: /* Assume x is a non-empty matrix. Return 0 if maximal pivot should not be
! 874: * used, and the matrix precision (min real precision of coeffs) otherwise.
! 875: */
! 876: static long
! 877: matprec(GEN x)
! 878: {
! 879: long tx,i,j,l, k = VERYBIGINT, lx = lg(x), ly = lg(x[1]);
! 880: GEN p1;
! 881: for (i=1; i<lx; i++)
! 882: for (j=1; j<ly; j++)
! 883: {
! 884: p1 = gmael(x,i,j); tx = typ(p1);
! 885: if (!is_scalar_t(tx)) return 0;
! 886: l = precision(p1); if (l && l<k) k = l;
! 887: }
! 888: return (k==VERYBIGINT)? 0: k;
! 889: }
! 890:
! 891: /* As above, returning 1 if the precision would be non-zero, 0 otherwise */
! 892: static long
! 893: use_maximal_pivot(GEN x)
! 894: {
! 895: long tx,i,j, lx = lg(x), ly = lg(x[1]);
! 896: GEN p1;
! 897: for (i=1; i<lx; i++)
! 898: for (j=1; j<ly; j++)
! 899: {
! 900: p1 = gmael(x,i,j); tx = typ(p1);
! 901: if (!is_scalar_t(tx)) return 0;
! 902: if (precision(p1)) return 1;
! 903: }
! 904: return 0;
! 905: }
! 906:
! 907: static GEN
! 908: check_b(GEN b, long nbli)
! 909: {
! 910: GEN col;
! 911: if (!b) return idmat(nbli);
! 912: col = (typ(b) == t_MAT)? (GEN)b[1]: b;
! 913: if (nbli != lg(col)-1) err(talker,"incompatible matrix dimensions in gauss");
! 914: return dummycopy(b);
! 915: }
! 916:
! 917: /* C = A^(-1)(tB) where A, B, C are integral, A is upper triangular, t t_INT */
! 918: GEN
! 919: gauss_triangle_i(GEN A, GEN B, GEN t)
! 920: {
! 921: long n = lg(A)-1, i,j,k;
! 922: GEN m, c = cgetg(n+1,t_MAT);
! 923:
! 924: if (!n) return c;
! 925: for (k=1; k<=n; k++)
! 926: {
! 927: GEN u = cgetg(n+1, t_COL), b = (GEN)B[k];
! 928: ulong av = avma;
! 929: c[k] = (long)u; m = mulii((GEN)b[n],t);
! 930: u[n] = (long)gerepileuptoint(av, divii(m, gcoeff(A,n,n)));
! 931: for (i=n-1; i>0; i--)
! 932: {
! 933: av = avma; m = mulii(negi((GEN)b[i]),t);
! 934: for (j=i+1; j<=n; j++)
! 935: m = addii(m, mulii(gcoeff(A,i,j),(GEN) u[j]));
! 936: u[i] = (long)gerepileuptoint(av, divii(negi(m), gcoeff(A,i,i)));
! 937: }
! 938: }
! 939: return c;
! 940: }
! 941:
! 942: GEN
! 943: gauss_get_col(GEN a, GEN b, GEN p, long li)
! 944: {
! 945: GEN m, u=cgetg(li+1,t_COL);
! 946: long i,j;
! 947:
! 948: u[li] = ldiv((GEN)b[li], p);
! 949: for (i=li-1; i>0; i--)
! 950: {
! 951: m = gneg_i((GEN)b[i]);
! 952: for (j=i+1; j<=li; j++)
! 953: m = gadd(m, gmul(gcoeff(a,i,j), (GEN)u[j]));
! 954: u[i] = ldiv(gneg_i(m), gcoeff(a,i,i));
! 955: }
! 956: return u;
! 957: }
! 958:
! 959: GEN
! 960: Fp_gauss_get_col(GEN a, GEN b, GEN piv, long li, GEN p)
! 961: {
! 962: GEN m, u=cgetg(li+1,t_COL);
! 963: long i,j;
! 964:
! 965: u[li] = lresii(mulii((GEN)b[li], mpinvmod(piv,p)), p);
! 966: for (i=li-1; i>0; i--)
! 967: {
! 968: m = (GEN)b[i];
! 969: for (j=i+1; j<=li; j++)
! 970: m = subii(m, mulii(gcoeff(a,i,j), (GEN)u[j]));
! 971: m = resii(m, p);
! 972: u[i] = lresii(mulii(m, mpinvmod(gcoeff(a,i,i), p)), p);
! 973: }
! 974: return u;
! 975: }
! 976:
! 977: /* assume -p < a < p, return 1/a mod p */
! 978: static long
! 979: u_Fp_inv(long a, long p)
! 980: {
! 981: if (a < 0) a = p + a; /* pb with ulongs < 0 */
! 982: return u_invmod(a,p);
! 983: }
! 984:
! 985: GEN
! 986: u_Fp_gauss_get_col(GEN a, GEN b, long piv, long li, long p)
! 987: {
! 988: GEN u=cgetg(li+1,t_VECSMALL);
! 989: long i,j, m;
! 990:
! 991: u[li] = (b[li] * u_Fp_inv(piv,p)) % p;
! 992: if (u[li] < 0) u[li] += p;
! 993: for (i=li-1; i>0; i--)
! 994: {
! 995: m = b[i];
! 996: for (j=i+1; j<=li; j++) { m -= coeff(a,i,j) * u[j]; if (m & MASK) m %= p; }
! 997: u[i] = ((m%p) * u_Fp_inv(coeff(a,i,i), p)) % p;
! 998: if (u[i] < 0) u[i] += p;
! 999: }
! 1000: return u;
! 1001: }
! 1002:
! 1003: /* bk += m * bi */
! 1004: static void
! 1005: _addmul(GEN b, long k, long i, GEN m)
! 1006: {
! 1007: b[k] = ladd((GEN)b[k], gmul(m, (GEN)b[i]));
! 1008: }
! 1009:
! 1010: /* same, reduce mod p */
! 1011: static void
! 1012: _Fp_addmul(GEN b, long k, long i, GEN m, GEN p)
! 1013: {
! 1014: if (lgefint(b[i]) > lgefint(p)) b[i] = lresii((GEN)b[i], p);
! 1015: b[k] = laddii((GEN)b[k], mulii(m, (GEN)b[i]));
! 1016: }
! 1017:
! 1018: static void
! 1019: _u_Fp_addmul(GEN b, long k, long i, long m, long p)
! 1020: {
! 1021: long a;
! 1022: if (b[i] & MASK) b[i] %= p;
! 1023: a = b[k] + m * b[i];
! 1024: if (a & MASK) a %= p;
! 1025: b[k] = a;
! 1026: }
! 1027:
! 1028: /* Gaussan Elimination. Compute a^(-1)*b
! 1029: * b is a matrix or column vector, NULL meaning: take the identity matrix
! 1030: * If a and b are empty, the result is the empty matrix.
! 1031: *
! 1032: * li: nb lines of a and b
! 1033: * aco: nb columns of a
! 1034: * bco: nb columns of b (if matrix)
! 1035: *
! 1036: * li > aco is allowed if b = NULL, in which case return c such that c a = Id */
! 1037: GEN
! 1038: gauss_intern(GEN a, GEN b)
! 1039: {
! 1040: long inexact,iscol,i,j,k,av,lim,li,bco, aco = lg(a)-1;
! 1041: GEN p,m,u;
! 1042:
! 1043: if (typ(a)!=t_MAT) err(mattype1,"gauss");
! 1044: if (b && typ(b)!=t_COL && typ(b)!=t_MAT) err(typeer,"gauss");
! 1045: if (!aco)
! 1046: {
! 1047: if (b && lg(b)!=1) err(consister,"gauss");
! 1048: if (DEBUGLEVEL) err(warner,"in Gauss lg(a)=1 lg(b)=%ld", b?1:-1);
! 1049: return cgetg(1,t_MAT);
! 1050: }
! 1051: av=avma; lim=stack_lim(av,1);
! 1052: li = lg(a[1])-1;
! 1053: if (li != aco && (li < aco || b)) err(mattype1,"gauss");
! 1054: a = dummycopy(a);
! 1055: b = check_b(b,li); bco = lg(b)-1;
! 1056: inexact = use_maximal_pivot(a);
! 1057: iscol = (typ(b)==t_COL);
! 1058: if(DEBUGLEVEL>4)
! 1059: fprintferr("Entering gauss with inexact=%ld iscol=%ld\n",inexact,iscol);
! 1060:
! 1061: p = NULL; /* gcc -Wall */
! 1062: for (i=1; i<=aco; i++)
! 1063: {
! 1064: /* k is the line where we find the pivot */
! 1065: p = gcoeff(a,i,i); k = i;
! 1066: if (inexact) /* maximal pivot */
! 1067: {
! 1068: long e, ex = gexpo(p);
! 1069: for (j=i+1; j<=li; j++)
! 1070: {
! 1071: e = gexpo(gcoeff(a,j,i));
! 1072: if (e > ex) { ex=e; k=j; }
! 1073: }
! 1074: if (gcmp0(gcoeff(a,k,i))) return NULL;
! 1075: }
! 1076: else if (gcmp0(p)) /* first non-zero pivot */
! 1077: {
! 1078: do k++; while (k<=li && gcmp0(gcoeff(a,k,i)));
! 1079: if (k>li) return NULL;
! 1080: }
! 1081:
! 1082: /* if (k!=i), exchange the lines s.t. k = i */
! 1083: if (k != i)
! 1084: {
! 1085: for (j=i; j<=aco; j++) swap(coeff(a,i,j), coeff(a,k,j));
! 1086: if (iscol) { swap(b[i],b[k]); }
! 1087: else
! 1088: for (j=1; j<=bco; j++) swap(coeff(b,i,j), coeff(b,k,j));
! 1089: p = gcoeff(a,i,i);
! 1090: }
! 1091: if (i == aco) break;
! 1092:
! 1093: for (k=i+1; k<=li; k++)
! 1094: {
! 1095: m=gcoeff(a,k,i);
! 1096: if (!gcmp0(m))
! 1097: {
! 1098: m = gneg_i(gdiv(m,p));
! 1099: for (j=i+1; j<=aco; j++) _addmul((GEN)a[j],k,i,m);
! 1100: if (iscol) _addmul(b,k,i,m);
! 1101: else
! 1102: for (j=1; j<=bco; j++) _addmul((GEN)b[j],k,i,m);
! 1103: }
! 1104: }
! 1105: if (low_stack(lim, stack_lim(av,1)))
! 1106: {
! 1107: GEN *gptr[2]; gptr[0]=&a; gptr[1]=&b;
! 1108: if(DEBUGMEM>1) err(warnmem,"gauss. i=%ld",i);
! 1109: gerepilemany(av,gptr,2);
! 1110: }
! 1111: }
! 1112:
! 1113: if(DEBUGLEVEL>4) fprintferr("Solving the triangular system\n");
! 1114: if (iscol) u = gauss_get_col(a,b,p,aco);
! 1115: else
! 1116: {
! 1117: long av1 = avma;
! 1118: lim = stack_lim(av1,1); u=cgetg(bco+1,t_MAT);
! 1119: for (j=2; j<=bco; j++) u[j] = zero; /* dummy */
! 1120: for (j=1; j<=bco; j++)
! 1121: {
! 1122: u[j] = (long)gauss_get_col(a,(GEN)b[j],p,aco);
! 1123: if (low_stack(lim, stack_lim(av1,1)))
! 1124: {
! 1125: if(DEBUGMEM>1) err(warnmem,"gauss[2]. j=%ld", j);
! 1126: u = gerepilecopy(av1, u);
! 1127: }
! 1128: }
! 1129: }
! 1130: return gerepilecopy(av, u);
! 1131: }
! 1132:
! 1133: GEN
! 1134: gauss(GEN a, GEN b)
! 1135: {
! 1136: GEN z = gauss_intern(a,b);
! 1137: if (!z) err(matinv1);
! 1138: return z;
! 1139: }
! 1140:
! 1141: GEN
! 1142: u_FpM_gauss(GEN a, GEN b, long p)
! 1143: {
! 1144: long piv,m,iscol,i,j,k,li,bco, aco = lg(a)-1;
! 1145: GEN u;
! 1146:
! 1147: if (!aco) return cgetg(1,t_MAT);
! 1148: li = lg(a[1])-1;
! 1149: bco = lg(b)-1;
! 1150: iscol = (typ(b)==t_COL);
! 1151: piv = 0; /* gcc -Wall */
! 1152: for (i=1; i<=aco; i++)
! 1153: {
! 1154: /* k is the line where we find the pivot */
! 1155: coeff(a,i,i) = coeff(a,i,i) % p;
! 1156: piv = coeff(a,i,i); k = i;
! 1157: if (!piv)
! 1158: {
! 1159: for (k++; k <= li; k++)
! 1160: {
! 1161: coeff(a,k,i) %= p;
! 1162: if (coeff(a,k,i)) break;
! 1163: }
! 1164: if (k>li) return NULL;
! 1165: }
! 1166:
! 1167: /* if (k!=i), exchange the lines s.t. k = i */
! 1168: if (k != i)
! 1169: {
! 1170: for (j=i; j<=aco; j++) swap(coeff(a,i,j), coeff(a,k,j));
! 1171: if (iscol) { swap(b[i],b[k]); }
! 1172: else
! 1173: for (j=1; j<=bco; j++) swap(coeff(b,i,j), coeff(b,k,j));
! 1174: piv = coeff(a,i,i);
! 1175: }
! 1176: if (i == aco) break;
! 1177:
! 1178: for (k=i+1; k<=li; k++)
! 1179: {
! 1180: coeff(a,k,i) %= p;
! 1181: m = coeff(a,k,i); coeff(a,k,i) = 0;
! 1182: if (m)
! 1183: {
! 1184: m = - (m * u_Fp_inv(piv,p)) % p;
! 1185: for (j=i+1; j<=aco; j++) _u_Fp_addmul((GEN)a[j],k,i,m, p);
! 1186: if (iscol) _u_Fp_addmul(b,k,i,m, p);
! 1187: else
! 1188: for (j=1; j<=bco; j++) _u_Fp_addmul((GEN)b[j],k,i,m, p);
! 1189: }
! 1190: }
! 1191: }
! 1192: if (iscol) u = u_Fp_gauss_get_col(a,b,piv,aco,p);
! 1193: else
! 1194: {
! 1195: u=cgetg(bco+1,t_MAT);
! 1196: for (j=1; j<=bco; j++)
! 1197: u[j] = (long)u_Fp_gauss_get_col(a,(GEN)b[j],piv,aco,p);
! 1198: }
! 1199: return u;
! 1200: }
! 1201:
! 1202: GEN
! 1203: FpM_gauss(GEN a, GEN b, GEN p)
! 1204: {
! 1205: long iscol,i,j,k,av,lim,li,bco, aco = lg(a)-1;
! 1206: GEN piv,m,u;
! 1207:
! 1208: if (typ(a)!=t_MAT) err(mattype1,"gauss");
! 1209: if (b && typ(b)!=t_COL && typ(b)!=t_MAT) err(typeer,"gauss");
! 1210: if (!aco)
! 1211: {
! 1212: if (b && lg(b)!=1) err(consister,"gauss");
! 1213: if (DEBUGLEVEL) err(warner,"in Gauss lg(a)=1 lg(b)=%ld", b?1:-1);
! 1214: return cgetg(1,t_MAT);
! 1215: }
! 1216: li = lg(a[1])-1;
! 1217: if (li != aco && (li < aco || b)) err(mattype1,"gauss");
! 1218: b = check_b(b,li); av = avma;
! 1219: if (OK_ULONG(p))
! 1220: {
! 1221: ulong pp=p[2];
! 1222: a = u_Fp_FpM(a, pp);
! 1223: b = u_Fp_FpM(b, pp);
! 1224: u = u_FpM_gauss(a,b, pp);
! 1225: return gerepileupto(av, small_to_mat(u));
! 1226: }
! 1227: lim = stack_lim(av,1);
! 1228: a = dummycopy(a);
! 1229: bco = lg(b)-1;
! 1230: iscol = (typ(b)==t_COL);
! 1231: piv = NULL; /* gcc -Wall */
! 1232: for (i=1; i<=aco; i++)
! 1233: {
! 1234: /* k is the line where we find the pivot */
! 1235: coeff(a,i,i) = lresii(gcoeff(a,i,i), p);
! 1236: piv = gcoeff(a,i,i); k = i;
! 1237: if (!signe(piv))
! 1238: {
! 1239: for (k++; k <= li; k++)
! 1240: {
! 1241: coeff(a,k,i) = lresii(gcoeff(a,k,i), p);
! 1242: if (signe(coeff(a,k,i))) break;
! 1243: }
! 1244: if (k>li) return NULL;
! 1245: }
! 1246:
! 1247: /* if (k!=i), exchange the lines s.t. k = i */
! 1248: if (k != i)
! 1249: {
! 1250: for (j=i; j<=aco; j++) swap(coeff(a,i,j), coeff(a,k,j));
! 1251: if (iscol) { swap(b[i],b[k]); }
! 1252: else
! 1253: for (j=1; j<=bco; j++) swap(coeff(b,i,j), coeff(b,k,j));
! 1254: piv = gcoeff(a,i,i);
! 1255: }
! 1256: if (i == aco) break;
! 1257:
! 1258: for (k=i+1; k<=li; k++)
! 1259: {
! 1260: coeff(a,k,i) = lresii(gcoeff(a,k,i), p);
! 1261: m = gcoeff(a,k,i); coeff(a,k,i) = zero;
! 1262: if (signe(m))
! 1263: {
! 1264: m = mulii(m, mpinvmod(piv,p));
! 1265: m = resii(negi(m), p);
! 1266: for (j=i+1; j<=aco; j++) _Fp_addmul((GEN)a[j],k,i,m, p);
! 1267: if (iscol) _Fp_addmul(b,k,i,m, p);
! 1268: else
! 1269: for (j=1; j<=bco; j++) _Fp_addmul((GEN)b[j],k,i,m, p);
! 1270: }
! 1271: }
! 1272: if (low_stack(lim, stack_lim(av,1)))
! 1273: {
! 1274: GEN *gptr[2]; gptr[0]=&a; gptr[1]=&b;
! 1275: if(DEBUGMEM>1) err(warnmem,"FpM_gauss. i=%ld",i);
! 1276: gerepilemany(av,gptr,2);
! 1277: }
! 1278: }
! 1279:
! 1280: if(DEBUGLEVEL>4) fprintferr("Solving the triangular system\n");
! 1281: if (iscol) u = Fp_gauss_get_col(a,b,piv,aco,p);
! 1282: else
! 1283: {
! 1284: long av1 = avma;
! 1285: lim = stack_lim(av1,1); u=cgetg(bco+1,t_MAT);
! 1286: for (j=2; j<=bco; j++) u[j] = zero; /* dummy */
! 1287: for (j=1; j<=bco; j++)
! 1288: {
! 1289: u[j] = (long)Fp_gauss_get_col(a,(GEN)b[j],piv,aco,p);
! 1290: if (low_stack(lim, stack_lim(av1,1)))
! 1291: {
! 1292: if(DEBUGMEM>1) err(warnmem,"gauss[2]. j=%ld", j);
! 1293: u = gerepilecopy(av1, u);
! 1294: }
! 1295: }
! 1296: }
! 1297: return gerepilecopy(av, u);
! 1298: }
! 1299:
! 1300: GEN
! 1301: FpM_inv(GEN x, GEN p) { return FpM_gauss(x, NULL, p); }
! 1302:
! 1303: static GEN
! 1304: u_idmat(long n)
! 1305: {
! 1306: GEN y = cgetg(n+1,t_MAT);
! 1307: long i,j;
! 1308: if (n < 0) err(talker,"negative size in u_idmat");
! 1309: for (i=1; i<=n; i++)
! 1310: {
! 1311: y[i]=lgetg(n+1,t_VECSMALL);
! 1312: for (j=1; j<=n; j++) coeff(y,j,i) = (i==j)? 1: 0;
! 1313: }
! 1314: return y;
! 1315: }
! 1316:
! 1317: /* set y = x * y */
! 1318: #define HIGHWORD(a) ((a) >> BITS_IN_HALFULONG)
! 1319: static GEN
! 1320: u_FpM_Fp_mul_ip(GEN y, long x, long p)
! 1321: {
! 1322: int i,j, m = lg(y[1]), l = lg(y);
! 1323: if (HIGHWORD(x | p))
! 1324: for(j=1; j<l; j++)
! 1325: for(i=1; i<m; i++)
! 1326: coeff(y,i,j) = (long)mulssmod(coeff(y,i,j), x, p);
! 1327: else
! 1328: for(j=1; j<l; j++)
! 1329: for(i=1; i<m; i++)
! 1330: coeff(y,i,j) = (coeff(y,i,j) * x) % p;
! 1331: return y;
! 1332: }
! 1333:
! 1334: /* M integral, dM such that M' = dM M^-1 is integral [e.g det(M)]. Return M' */
! 1335: GEN
! 1336: ZM_inv(GEN M, GEN dM)
! 1337: {
! 1338: GEN ID,Hp,q,H;
! 1339: ulong p,dMp, av2, av = avma, lim = stack_lim(av,1);
! 1340: byteptr d = diffptr;
! 1341: long lM = lg(M), stable = 0;
! 1342:
! 1343: if (lM == 1) return cgetg(1,t_MAT);
! 1344: if (!dM) dM = det(M);
! 1345:
! 1346: av2 = avma;
! 1347: H = NULL;
! 1348: d += 3000; /* 27449 = prime(3000) */
! 1349: for(p = 27449; ; p+= *d++)
! 1350: {
! 1351: if (!*d) err(primer1);
! 1352: dMp = umodiu(dM,p);
! 1353: if (!dMp) continue;
! 1354: ID = u_idmat(lM-1);
! 1355: Hp = u_FpM_gauss(u_Fp_FpM(M,p), ID, p);
! 1356: if (dMp != 1) Hp = u_FpM_Fp_mul_ip(Hp, dMp, p);
! 1357:
! 1358: if (!H)
! 1359: {
! 1360: H = ZM_init_CRT(Hp, p);
! 1361: q = utoi(p);
! 1362: }
! 1363: else
! 1364: {
! 1365: GEN qp = muliu(q,p);
! 1366: stable = ZM_incremental_CRT(H, Hp, q,qp, p);
! 1367: q = qp;
! 1368: }
! 1369: if (DEBUGLEVEL>5) msgtimer("inverse mod %ld (stable=%ld)", p,stable);
! 1370: if (stable && isscalarmat(gmul(M, H), dM)) break; /* DONE */
! 1371:
! 1372: if (low_stack(lim, stack_lim(av,2)))
! 1373: {
! 1374: GEN *gptr[2]; gptr[0] = &H; gptr[1] = &q;
! 1375: if (DEBUGMEM>1) err(warnmem,"ZM_inv");
! 1376: gerepilemany(av2,gptr, 2);
! 1377: }
! 1378:
! 1379: }
! 1380: if (DEBUGLEVEL>5) msgtimer("ZM_inv done");
! 1381: return gerepilecopy(av, H);
! 1382: }
! 1383:
! 1384: /* same as above, M rational */
! 1385: GEN
! 1386: QM_inv(GEN M, GEN dM)
! 1387: {
! 1388: ulong av = avma;
! 1389: GEN cM = content(M);
! 1390: if (is_pm1(cM)) { avma = av; return ZM_inv(M,dM); }
! 1391: return gerepileupto(av, ZM_inv(gdiv(M,cM), gdiv(dM,cM)));
! 1392: }
! 1393:
! 1394: /* x a matrix with integer coefficients. Return a multiple of the determinant
! 1395: * of the lattice generated by the columns of x (to be used with hnfmod)
! 1396: */
! 1397: GEN
! 1398: detint(GEN x)
! 1399: {
! 1400: GEN pass,c,v,det1,piv,pivprec,vi,p1;
! 1401: long i,j,k,rg,n,m,m1,av=avma,av1,lim;
! 1402:
! 1403: if (typ(x)!=t_MAT) err(typeer,"detint");
! 1404: n=lg(x)-1; if (!n) return gun;
! 1405: m1=lg(x[1]); m=m1-1; lim=stack_lim(av,1);
! 1406: c=new_chunk(m1); for (k=1; k<=m; k++) c[k]=0;
! 1407: av1=avma; pass=cgetg(m1,t_MAT);
! 1408: for (j=1; j<=m; j++)
! 1409: {
! 1410: p1=cgetg(m1,t_COL); pass[j]=(long)p1;
! 1411: for (i=1; i<=m; i++) p1[i]=zero;
! 1412: }
! 1413: for (k=1; k<=n; k++)
! 1414: for (j=1; j<=m; j++)
! 1415: if (typ(gcoeff(x,j,k)) != t_INT)
! 1416: err(talker,"not an integer matrix in detint");
! 1417: v=cgetg(m1,t_COL);
! 1418: det1=gzero; piv=pivprec=gun;
! 1419: for (rg=0,k=1; k<=n; k++)
! 1420: {
! 1421: long t = 0;
! 1422: for (i=1; i<=m; i++)
! 1423: if (!c[i])
! 1424: {
! 1425: vi=mulii(piv,gcoeff(x,i,k));
! 1426: for (j=1; j<=m; j++)
! 1427: if (c[j]) vi=addii(vi,mulii(gcoeff(pass,i,j),gcoeff(x,j,k)));
! 1428: v[i]=(long)vi; if (!t && signe(vi)) t=i;
! 1429: }
! 1430: if (t)
! 1431: {
! 1432: if (rg == m-1)
! 1433: { det1=mppgcd((GEN)v[t],det1); c[t]=0; }
! 1434: else
! 1435: {
! 1436: rg++; pivprec = piv; piv=(GEN)v[t]; c[t]=k;
! 1437: for (i=1; i<=m; i++)
! 1438: if (!c[i])
! 1439: {
! 1440: GEN p2 = negi((GEN)v[i]);
! 1441: for (j=1; j<=m; j++)
! 1442: if (c[j] && j!=t)
! 1443: {
! 1444: p1 = addii(mulii(gcoeff(pass,i,j), piv),
! 1445: mulii(gcoeff(pass,t,j), p2));
! 1446: if (rg>1) p1 = divii(p1,pivprec);
! 1447: coeff(pass,i,j) = (long)p1;
! 1448: }
! 1449: coeff(pass,i,t) = (long)p2;
! 1450: }
! 1451: }
! 1452: }
! 1453: if (low_stack(lim, stack_lim(av,1)))
! 1454: {
! 1455: GEN *gptr[5];
! 1456: if(DEBUGMEM>1) err(warnmem,"detint. k=%ld",k);
! 1457: gptr[0]=&det1; gptr[1]=ϖ gptr[2]=&pivprec;
! 1458: gptr[3]=&pass; gptr[4]=&v; gerepilemany(av1,gptr,5);
! 1459: }
! 1460: }
! 1461: return gerepileupto(av, absi(det1));
! 1462: }
! 1463:
! 1464: static void
! 1465: gerepile_gauss_ker(GEN x, long m, long n, long k, long t, ulong av)
! 1466: {
! 1467: ulong tetpil = avma,A;
! 1468: long dec,u,i;
! 1469:
! 1470: if (DEBUGMEM > 1) err(warnmem,"gauss_pivot_ker. k=%ld, n=%ld",k,n);
! 1471: for (u=t+1; u<=m; u++) copyifstack(coeff(x,u,k), coeff(x,u,k));
! 1472: for (i=k+1; i<=n; i++)
! 1473: for (u=1; u<=m; u++) copyifstack(coeff(x,u,i), coeff(x,u,i));
! 1474:
! 1475: (void)gerepile(av,tetpil,NULL); dec = av-tetpil;
! 1476: for (u=t+1; u<=m; u++)
! 1477: {
! 1478: A=coeff(x,u,k);
! 1479: if (A<av && A>=bot) coeff(x,u,k)+=dec;
! 1480: }
! 1481: for (i=k+1; i<=n; i++)
! 1482: for (u=1; u<=m; u++)
! 1483: {
! 1484: A=coeff(x,u,i);
! 1485: if (A<av && A>=bot) coeff(x,u,i)+=dec;
! 1486: }
! 1487: }
! 1488:
! 1489: static void
! 1490: gerepile_gauss_FpM_ker(GEN x, GEN p, long m, long n, long k, long t, ulong av)
! 1491: {
! 1492: ulong tetpil = avma,A;
! 1493: long dec,u,i;
! 1494:
! 1495: if (DEBUGMEM > 1) err(warnmem,"gauss_pivot_ker. k=%ld, n=%ld",k,n);
! 1496: for (u=t+1; u<=m; u++)
! 1497: if (isonstack(coeff(x,u,k))) coeff(x,u,k) = lmodii(gcoeff(x,u,k),p);
! 1498: for (i=k+1; i<=n; i++)
! 1499: for (u=1; u<=m; u++)
! 1500: if (isonstack(coeff(x,u,i))) coeff(x,u,i) = lmodii(gcoeff(x,u,i),p);
! 1501:
! 1502: (void)gerepile(av,tetpil,NULL); dec = av-tetpil;
! 1503: for (u=t+1; u<=m; u++)
! 1504: {
! 1505: A=coeff(x,u,k);
! 1506: if (A<av && A>=bot) coeff(x,u,k)+=dec;
! 1507: }
! 1508: for (i=k+1; i<=n; i++)
! 1509: for (u=1; u<=m; u++)
! 1510: {
! 1511: A=coeff(x,u,i);
! 1512: if (A<av && A>=bot) coeff(x,u,i)+=dec;
! 1513: }
! 1514: }
! 1515:
! 1516: /* special gerepile for huge matrices */
! 1517:
! 1518: static void
! 1519: gerepile_gauss(GEN x,long m,long n,long k,long t,ulong av, long j, GEN c)
! 1520: {
! 1521: ulong tetpil = avma,A;
! 1522: long dec,u,i;
! 1523:
! 1524: if (DEBUGMEM > 1) err(warnmem,"gauss_pivot. k=%ld, n=%ld",k,n);
! 1525: for (u=t+1; u<=m; u++)
! 1526: if (u==j || !c[u]) copyifstack(coeff(x,u,k), coeff(x,u,k));
! 1527: for (u=1; u<=m; u++)
! 1528: if (u==j || !c[u])
! 1529: for (i=k+1; i<=n; i++) copyifstack(coeff(x,u,i), coeff(x,u,i));
! 1530:
! 1531: (void)gerepile(av,tetpil,NULL); dec = av-tetpil;
! 1532: for (u=t+1; u<=m; u++)
! 1533: if (u==j || !c[u])
! 1534: {
! 1535: A=coeff(x,u,k);
! 1536: if (A<av && A>=bot) coeff(x,u,k)+=dec;
! 1537: }
! 1538: for (u=1; u<=m; u++)
! 1539: if (u==j || !c[u])
! 1540: for (i=k+1; i<=n; i++)
! 1541: {
! 1542: A=coeff(x,u,i);
! 1543: if (A<av && A>=bot) coeff(x,u,i)+=dec;
! 1544: }
! 1545: }
! 1546:
! 1547: /*******************************************************************/
! 1548: /* */
! 1549: /* KERNEL of an m x n matrix */
! 1550: /* return n - rk(x) linearly independant vectors */
! 1551: /* */
! 1552: /*******************************************************************/
! 1553:
! 1554: /* x has INTEGER coefficients */
! 1555: GEN
! 1556: keri(GEN x)
! 1557: {
! 1558: GEN c,d,y,p,pp;
! 1559: long i,j,k,r,t,n,m,av,av0,tetpil,lim;
! 1560:
! 1561: if (typ(x)!=t_MAT) err(typeer,"keri");
! 1562: n=lg(x)-1; if (!n) return cgetg(1,t_MAT);
! 1563:
! 1564: av0=avma; m=lg(x[1])-1; r=0;
! 1565: pp=cgetg(n+1,t_COL);
! 1566: x=dummycopy(x); p=gun;
! 1567: c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0;
! 1568: d=new_chunk(n+1); av=avma; lim=stack_lim(av,1);
! 1569: for (k=1; k<=n; k++)
! 1570: {
! 1571: j=1;
! 1572: while (j<=m && (c[j] || !signe(gcoeff(x,j,k))) ) j++;
! 1573: if (j>m)
! 1574: {
! 1575: r++; d[k]=0;
! 1576: for(j=1; j<k; j++)
! 1577: if (d[j]) coeff(x,d[j],k) = lclone(gcoeff(x,d[j],k));
! 1578: pp[k]=lclone(p);
! 1579: }
! 1580: else
! 1581: {
! 1582: GEN p0 = p;
! 1583: long av1;
! 1584:
! 1585: c[j]=k; d[k]=j; p = gcoeff(x,j,k);
! 1586:
! 1587: for (t=1; t<=m; t++)
! 1588: if (t!=j)
! 1589: {
! 1590: GEN q=gcoeff(x,t,k), p1,p2;
! 1591: for (i=k+1; i<=n; i++)
! 1592: {
! 1593: av1=avma; (void)new_chunk(lgefint(p0));
! 1594: p1=mulii(q,gcoeff(x,j,i));
! 1595: p2=mulii(p,gcoeff(x,t,i));
! 1596: p1=subii(p2,p1); avma=av1;
! 1597: coeff(x,t,i) = ldivii(p1,p0);
! 1598: }
! 1599: if (low_stack(lim, stack_lim(av,1)))
! 1600: {
! 1601: p1 = gclone(p);
! 1602: gerepile_gauss_ker(x,m,n,k,t,av);
! 1603: p = gcopy(p1); gunclone(p1);
! 1604: }
! 1605: }
! 1606: }
! 1607: }
! 1608: if (!r) { avma=av0; y=cgetg(1,t_MAT); return y; }
! 1609:
! 1610: /* non trivial kernel */
! 1611: tetpil=avma; y=cgetg(r+1,t_MAT);
! 1612: for (j=k=1; j<=r; j++,k++)
! 1613: {
! 1614: p = cgetg(n+1, t_COL);
! 1615: y[j]=(long)p; while (d[k]) k++;
! 1616: for (i=1; i<k; i++)
! 1617: if (d[i])
! 1618: {
! 1619: c=gcoeff(x,d[i],k);
! 1620: p[i] = (long) forcecopy(c); gunclone(c);
! 1621: }
! 1622: else
! 1623: p[i] = zero;
! 1624: p[k]=lnegi((GEN)pp[k]); gunclone((GEN)pp[k]);
! 1625: for (i=k+1; i<=n; i++) p[i]=zero;
! 1626: }
! 1627: return gerepile(av0,tetpil,y);
! 1628: }
! 1629:
! 1630: GEN
! 1631: deplin(GEN x)
! 1632: {
! 1633: long i,j,k,t,nc,nl, av=avma;
! 1634: GEN y,q,c,l,d;
! 1635:
! 1636: if (typ(x) != t_MAT) err(typeer,"deplin");
! 1637: nc=lg(x)-1; if (!nc) err(talker,"empty matrix in deplin");
! 1638: nl=lg(x[1])-1;
! 1639: c=new_chunk(nl+1);
! 1640: l=new_chunk(nc+1);
! 1641: d=cgetg(nl+1,t_VEC);
! 1642: for (i=1; i<=nl; i++) { d[i]=un; c[i]=0; }
! 1643: k=1; t=1;
! 1644: while (t<=nl && k<=nc)
! 1645: {
! 1646: for (j=1; j<k; j++)
! 1647: for (i=1; i<=nl; i++)
! 1648: if (i!=l[j])
! 1649: coeff(x,i,k)=lsub(gmul((GEN) d[j],gcoeff(x,i,k)),
! 1650: gmul(gcoeff(x,i,j),gcoeff(x,l[j],k)));
! 1651: t=1;
! 1652: while ( t<=nl && (c[t] || gcmp0(gcoeff(x,t,k))) ) t++;
! 1653: if (t<=nl)
! 1654: {
! 1655: d[k]=coeff(x,t,k);
! 1656: c[t]=k; l[k++]=t;
! 1657: }
! 1658: }
! 1659: if (k>nc)
! 1660: {
! 1661: avma=av; y=cgetg(nc+1,t_COL);
! 1662: for (j=1; j<=nc; j++) y[j]=zero;
! 1663: return y;
! 1664: }
! 1665: y=cgetg(nc+1,t_COL);
! 1666: y[1]=(k>1)? coeff(x,l[1],k): un;
! 1667: for (q=gun,j=2; j<k; j++)
! 1668: {
! 1669: q=gmul(q,(GEN) d[j-1]);
! 1670: y[j]=lmul(gcoeff(x,l[j],k),q);
! 1671: }
! 1672: if (k>1) y[k]=lneg(gmul(q,(GEN) d[k-1]));
! 1673: for (j=k+1; j<=nc; j++) y[j]=zero;
! 1674: d=content(y); t=avma;
! 1675: return gerepile(av,t,gdiv(y,d));
! 1676: }
! 1677:
! 1678: /*******************************************************************/
! 1679: /* */
! 1680: /* GAUSS REDUCTION OF MATRICES (m lines x n cols) */
! 1681: /* (kernel, image, complementary image, rank) */
! 1682: /* */
! 1683: /*******************************************************************/
! 1684: static long gauss_ex;
! 1685: static int (*gauss_is_zero)(GEN);
! 1686:
! 1687: static int
! 1688: real0(GEN x)
! 1689: {
! 1690: return gcmp0(x) || (gexpo(x) < gauss_ex);
! 1691: }
! 1692:
! 1693: static void
! 1694: gauss_get_prec(GEN x, long prec)
! 1695: {
! 1696: long pr = matprec(x);
! 1697:
! 1698: if (!pr) { gauss_is_zero = &gcmp0; return; }
! 1699: if (pr > prec) prec = pr;
! 1700:
! 1701: gauss_ex = - (long)(0.85 * bit_accuracy(prec));
! 1702: gauss_is_zero = &real0;
! 1703: }
! 1704:
! 1705: static long
! 1706: gauss_get_pivot_NZ(GEN x, GEN x0/* unused */, GEN c, long i0)
! 1707: {
! 1708: long i,lx = lg(x);
! 1709: if (c)
! 1710: for (i=i0; i<lx; i++)
! 1711: {
! 1712: if (c[i]) continue;
! 1713: if (!gcmp0((GEN)x[i])) break;
! 1714: }
! 1715: else
! 1716: for (i=i0; i<lx; i++)
! 1717: if (!gcmp0((GEN)x[i])) break;
! 1718: return i;
! 1719:
! 1720: }
! 1721:
! 1722: /* ~ 0 compared to reference y */
! 1723: static int
! 1724: approx_0(GEN x, GEN y)
! 1725: {
! 1726: long tx = typ(x);
! 1727: if (tx == t_COMPLEX)
! 1728: return approx_0((GEN)x[1], y) && approx_0((GEN)x[2], y);
! 1729: return gcmp0(x) ||
! 1730: (tx == t_REAL && gexpo(y) - gexpo(x) > bit_accuracy(lg(x)));
! 1731: }
! 1732:
! 1733: static long
! 1734: gauss_get_pivot_max(GEN x, GEN x0, GEN c, long i0)
! 1735: {
! 1736: long i,e, k = i0, ex = -HIGHEXPOBIT, lx = lg(x);
! 1737: GEN p,r;
! 1738: if (c)
! 1739: for (i=i0; i<lx; i++)
! 1740: {
! 1741: if (c[i]) continue;
! 1742: e = gexpo((GEN)x[i]);
! 1743: if (e > ex) { ex=e; k=i; }
! 1744: }
! 1745: else
! 1746: for (i=i0; i<lx; i++)
! 1747: {
! 1748: e = gexpo((GEN)x[i]);
! 1749: if (e > ex) { ex=e; k=i; }
! 1750: }
! 1751: p = (GEN)x[k];
! 1752: r = (GEN)x0[k]; if (isexactzero(r)) r = x0;
! 1753: return approx_0(p, r)? i: k;
! 1754: }
! 1755:
! 1756: /* return the transform of x under a standard Gauss pivot. r = dim ker(x).
! 1757: * d[k] contains the index of the first non-zero pivot in column k
! 1758: * If a != NULL, use x - a Id instead (for eigen)
! 1759: */
! 1760: static GEN
! 1761: gauss_pivot_ker(GEN x0, GEN a, GEN *dd, long *rr)
! 1762: {
! 1763: GEN x,c,d,p,mun;
! 1764: long i,j,k,r,t,n,m,av,lim;
! 1765: long (*get_pivot)(GEN,GEN,GEN,long);
! 1766:
! 1767: if (typ(x0)!=t_MAT) err(typeer,"gauss_pivot");
! 1768: n=lg(x0)-1; if (!n) { *dd=NULL; *rr=0; return cgetg(1,t_MAT); }
! 1769: m=lg(x0[1])-1; r=0;
! 1770:
! 1771: x = dummycopy(x0); mun = negi(gun);
! 1772: if (a)
! 1773: {
! 1774: if (n != m) err(consister,"gauss_pivot_ker");
! 1775: for (k=1; k<=n; k++) coeff(x,k,k) = lsub(gcoeff(x,k,k), a);
! 1776: }
! 1777: get_pivot = use_maximal_pivot(x)? gauss_get_pivot_max: gauss_get_pivot_NZ;
! 1778: c=cgetg(m+1,t_VECSMALL); for (k=1; k<=m; k++) c[k]=0;
! 1779: d=cgetg(n+1,t_VECSMALL);
! 1780: av=avma; lim=stack_lim(av,1);
! 1781: for (k=1; k<=n; k++)
! 1782: {
! 1783: j = get_pivot((GEN)x[k], (GEN)x0[k], c, 1);
! 1784: if (j > m)
! 1785: {
! 1786: r++; d[k]=0;
! 1787: for(j=1; j<k; j++)
! 1788: if (d[j]) coeff(x,d[j],k) = lclone(gcoeff(x,d[j],k));
! 1789: }
! 1790: else
! 1791: { /* pivot for column k on row j */
! 1792: c[j]=k; d[k]=j; p = gdiv(mun,gcoeff(x,j,k));
! 1793: coeff(x,j,k) = (long)mun;
! 1794: /* x[j,] /= - x[j,k] */
! 1795: for (i=k+1; i<=n; i++)
! 1796: coeff(x,j,i) = lmul(p,gcoeff(x,j,i));
! 1797: for (t=1; t<=m; t++)
! 1798: if (t!=j)
! 1799: { /* x[t,] -= 1 / x[j,k] x[j,] */
! 1800: p=gcoeff(x,t,k); coeff(x,t,k)=zero;
! 1801: for (i=k+1; i<=n; i++)
! 1802: coeff(x,t,i) = ladd(gcoeff(x,t,i),gmul(p,gcoeff(x,j,i)));
! 1803: if (low_stack(lim, stack_lim(av,1)))
! 1804: gerepile_gauss_ker(x,m,n,k,t,av);
! 1805: }
! 1806: }
! 1807: }
! 1808: *dd=d; *rr=r; return x;
! 1809: }
! 1810:
! 1811: /* r = dim ker(x).
! 1812: * d[k] contains the index of the first non-zero pivot in column k
! 1813: */
! 1814: static void
! 1815: gauss_pivot(GEN x0, GEN *dd, long *rr)
! 1816: {
! 1817: GEN x,c,d,d0,mun,p;
! 1818: long i,j,k,r,t,n,m,av,lim;
! 1819: long (*get_pivot)(GEN,GEN,GEN,long);
! 1820:
! 1821: if (typ(x0)!=t_MAT) err(typeer,"gauss_pivot");
! 1822: n=lg(x0)-1; if (!n) { *dd=NULL; *rr=0; return; }
! 1823:
! 1824: d0 = cgetg(n+1, t_VECSMALL);
! 1825: if (use_maximal_pivot(x0))
! 1826: { /* put exact columns first, then largest inexact ones */
! 1827: get_pivot = gauss_get_pivot_max;
! 1828: for (k=1; k<=n; k++)
! 1829: d0[k] = isinexactreal((GEN)x0[k])? -gexpo((GEN)x0[k]): -VERYBIGINT;
! 1830: d0 = gen_sort(d0, cmp_C|cmp_IND, NULL);
! 1831: x0 = vecextract_p(x0, d0);
! 1832: }
! 1833: else
! 1834: {
! 1835: get_pivot = gauss_get_pivot_NZ;
! 1836: for (k=1; k<=n; k++) d0[k] = k;
! 1837: }
! 1838: x = dummycopy(x0); mun = negi(gun);
! 1839: m=lg(x[1])-1; r=0;
! 1840: c=cgetg(m+1, t_VECSMALL); for (k=1; k<=m; k++) c[k]=0;
! 1841: d=(GEN)gpmalloc((n+1)*sizeof(long)); av=avma; lim=stack_lim(av,1);
! 1842: for (k=1; k<=n; k++)
! 1843: {
! 1844: j = get_pivot((GEN)x[k], (GEN)x0[k], c, 1);
! 1845: if (j>m) { r++; d[d0[k]]=0; }
! 1846: else
! 1847: {
! 1848: c[j]=k; d[d0[k]]=j; p = gdiv(mun,gcoeff(x,j,k));
! 1849: for (i=k+1; i<=n; i++)
! 1850: coeff(x,j,i) = lmul(p,gcoeff(x,j,i));
! 1851:
! 1852: for (t=1; t<=m; t++)
! 1853: if (!c[t]) /* no pivot on that line yet */
! 1854: {
! 1855: p=gcoeff(x,t,k); coeff(x,t,k)=zero;
! 1856: for (i=k+1; i<=n; i++)
! 1857: coeff(x,t,i) = ladd(gcoeff(x,t,i), gmul(p,gcoeff(x,j,i)));
! 1858: if (low_stack(lim, stack_lim(av,1)))
! 1859: gerepile_gauss(x,m,n,k,t,av,j,c);
! 1860: }
! 1861: for (i=k; i<=n; i++) coeff(x,j,i) = zero; /* dummy */
! 1862: }
! 1863: }
! 1864: *dd=d; *rr=r;
! 1865: }
! 1866:
! 1867: /* compute ker(x - aI) */
! 1868: static GEN
! 1869: ker0(GEN x, GEN a)
! 1870: {
! 1871: GEN d,y;
! 1872: long i,j,k,r,n, av = avma, tetpil;
! 1873:
! 1874: x = gauss_pivot_ker(x,a,&d,&r);
! 1875: if (!r) { avma=av; return cgetg(1,t_MAT); }
! 1876: n = lg(x)-1; tetpil=avma; y=cgetg(r+1,t_MAT);
! 1877: for (j=k=1; j<=r; j++,k++)
! 1878: {
! 1879: GEN p = cgetg(n+1,t_COL);
! 1880:
! 1881: y[j]=(long)p; while (d[k]) k++;
! 1882: for (i=1; i<k; i++)
! 1883: if (d[i])
! 1884: {
! 1885: GEN p1=gcoeff(x,d[i],k);
! 1886: p[i] = (long)forcecopy(p1); gunclone(p1);
! 1887: }
! 1888: else
! 1889: p[i] = zero;
! 1890: p[k]=un; for (i=k+1; i<=n; i++) p[i]=zero;
! 1891: }
! 1892: return gerepile(av,tetpil,y);
! 1893: }
! 1894:
! 1895: GEN
! 1896: ker(GEN x) /* Programme pour types exacts */
! 1897: {
! 1898: return ker0(x,NULL);
! 1899: }
! 1900:
! 1901: GEN
! 1902: matker0(GEN x,long flag)
! 1903: {
! 1904: return flag? keri(x): ker(x);
! 1905: }
! 1906:
! 1907: GEN
! 1908: image(GEN x)
! 1909: {
! 1910: GEN d,y;
! 1911: long j,k,r, av = avma;
! 1912:
! 1913: gauss_pivot(x,&d,&r);
! 1914:
! 1915: /* r = dim ker(x) */
! 1916: if (!r) { avma=av; if (d) free(d); return gcopy(x); }
! 1917:
! 1918: /* r = dim Im(x) */
! 1919: r = lg(x)-1 - r; avma=av;
! 1920: y=cgetg(r+1,t_MAT);
! 1921: for (j=k=1; j<=r; k++)
! 1922: if (d[k]) y[j++] = lcopy((GEN)x[k]);
! 1923: free(d); return y;
! 1924: }
! 1925:
! 1926: GEN
! 1927: imagecompl(GEN x)
! 1928: {
! 1929: GEN d,y;
! 1930: long j,i,r,av = avma;
! 1931:
! 1932: gauss_pivot(x,&d,&r);
! 1933: avma=av; y=cgetg(r+1,t_VEC);
! 1934: for (i=j=1; j<=r; i++)
! 1935: if (!d[i]) y[j++]=lstoi(i);
! 1936: if (d) free(d); return y;
! 1937: }
! 1938:
! 1939: /* for hnfspec: imagecompl(trans(x)) + image(trans(x)) */
! 1940: GEN
! 1941: imagecomplspec(GEN x, long *nlze)
! 1942: {
! 1943: GEN d,y;
! 1944: long i,j,k,l,r,av = avma;
! 1945:
! 1946: x = gtrans(x); l = lg(x);
! 1947: gauss_pivot(x,&d,&r);
! 1948: avma=av; y = cgetg(l,t_VECSMALL);
! 1949: for (i=j=1, k=r+1; i<l; i++)
! 1950: if (d[i]) y[k++]=i; else y[j++]=i;
! 1951: *nlze = r;
! 1952: if (d) free(d); return y;
! 1953: }
! 1954:
! 1955: static GEN
! 1956: sinverseimage(GEN mat, GEN y)
! 1957: {
! 1958: long av=avma,tetpil,i, nbcol = lg(mat);
! 1959: GEN col,p1 = cgetg(nbcol+1,t_MAT);
! 1960:
! 1961: if (nbcol==1) return NULL;
! 1962: if (lg(y) != lg(mat[1])) err(consister,"inverseimage");
! 1963:
! 1964: p1[nbcol] = (long)y;
! 1965: for (i=1; i<nbcol; i++) p1[i]=mat[i];
! 1966: p1 = ker(p1); i=lg(p1)-1;
! 1967: if (!i) return NULL;
! 1968:
! 1969: col = (GEN)p1[i]; p1 = (GEN) col[nbcol];
! 1970: if (gcmp0(p1)) return NULL;
! 1971:
! 1972: p1 = gneg_i(p1); setlg(col,nbcol); tetpil=avma;
! 1973: return gerepile(av,tetpil, gdiv(col, p1));
! 1974: }
! 1975:
! 1976: /* Calcule l'image reciproque de v par m */
! 1977: GEN
! 1978: inverseimage(GEN m,GEN v)
! 1979: {
! 1980: long av=avma,j,lv,tv=typ(v);
! 1981: GEN y,p1;
! 1982:
! 1983: if (typ(m)!=t_MAT) err(typeer,"inverseimage");
! 1984: if (tv==t_COL)
! 1985: {
! 1986: p1 = sinverseimage(m,v);
! 1987: if (p1) return p1;
! 1988: avma = av; return cgetg(1,t_MAT);
! 1989: }
! 1990: if (tv!=t_MAT) err(typeer,"inverseimage");
! 1991:
! 1992: lv=lg(v)-1; y=cgetg(lv+1,t_MAT);
! 1993: for (j=1; j<=lv; j++)
! 1994: {
! 1995: p1 = sinverseimage(m,(GEN)v[j]);
! 1996: if (!p1) { avma = av; return cgetg(1,t_MAT); }
! 1997: y[j] = (long)p1;
! 1998: }
! 1999: return y;
! 2000: }
! 2001:
! 2002: /* x is an n x k matrix, rank(x) = k <= n. Return an invertible n x n matrix
! 2003: * whose first k columns are given by x. If rank(x)<k, the result may be wrong
! 2004: */
! 2005: GEN
! 2006: suppl_intern(GEN x, GEN myid)
! 2007: {
! 2008: long av = avma, lx = lg(x), n,i,j;
! 2009: GEN y,p1;
! 2010: stackzone *zone;
! 2011:
! 2012: if (typ(x) != t_MAT) err(typeer,"suppl");
! 2013: if (lx==1) err(talker,"empty matrix in suppl");
! 2014: n=lg(x[1]); if (lx>n) err(suppler2);
! 2015: if (lx == n) return gcopy(x);
! 2016:
! 2017: zone = switch_stack(NULL, n*n);
! 2018: switch_stack(zone,1);
! 2019: y = myid? dummycopy(myid): idmat(n-1);
! 2020: switch_stack(zone,0);
! 2021: gauss_get_prec(x,0);
! 2022: for (i=1; i<lx; i++)
! 2023: {
! 2024: p1 = gauss(y,(GEN)x[i]); j=i;
! 2025: while (j<n && gauss_is_zero((GEN)p1[j])) j++;
! 2026: if (j>=n) err(suppler2);
! 2027: p1=(GEN)y[i]; y[i]=x[i]; if (i!=j) y[j]=(long)p1;
! 2028: }
! 2029: avma = av; y = gcopy(y);
! 2030: free(zone); return y;
! 2031: }
! 2032:
! 2033: GEN
! 2034: suppl(GEN x)
! 2035: {
! 2036: return suppl_intern(x,NULL);
! 2037: }
! 2038:
! 2039: GEN
! 2040: image2(GEN x)
! 2041: {
! 2042: long av=avma,tetpil,k,n,i;
! 2043: GEN p1,p2;
! 2044:
! 2045: if (typ(x)!=t_MAT) err(typeer,"image2");
! 2046: k=lg(x)-1; if (!k) return gcopy(x);
! 2047: n=lg(x[1])-1; p1=ker(x); k=lg(p1)-1;
! 2048: if (k) { p1=suppl(p1); n=lg(p1)-1; }
! 2049: else p1=idmat(n);
! 2050:
! 2051: tetpil=avma; p2=cgetg(n-k+1,t_MAT);
! 2052: for (i=k+1; i<=n; i++) p2[i-k]=lmul(x,(GEN) p1[i]);
! 2053: return gerepile(av,tetpil,p2);
! 2054: }
! 2055:
! 2056: GEN
! 2057: matimage0(GEN x,long flag)
! 2058: {
! 2059: switch(flag)
! 2060: {
! 2061: case 0: return image(x);
! 2062: case 1: return image2(x);
! 2063: default: err(flagerr,"matimage");
! 2064: }
! 2065: return NULL; /* not reached */
! 2066: }
! 2067:
! 2068: long
! 2069: rank(GEN x)
! 2070: {
! 2071: long av = avma, r;
! 2072: GEN d;
! 2073:
! 2074: gauss_pivot(x,&d,&r);
! 2075: /* yield r = dim ker(x) */
! 2076:
! 2077: avma=av; if (d) free(d);
! 2078: return lg(x)-1 - r;
! 2079: }
! 2080:
! 2081: static void FpM_gauss_pivot(GEN x, GEN p, GEN *dd, long *rr);
! 2082:
! 2083: /* if p != NULL, assume x integral and compute rank over Fp */
! 2084: static GEN
! 2085: indexrank0(GEN x, GEN p, int small)
! 2086: {
! 2087: long av = avma, i,j,n,r;
! 2088: GEN res,d,p1,p2;
! 2089:
! 2090: /* yield r = dim ker(x) */
! 2091: FpM_gauss_pivot(x,p,&d,&r);
! 2092:
! 2093: /* now r = dim Im(x) */
! 2094: n = lg(x)-1; r = n - r;
! 2095:
! 2096: avma=av; res=cgetg(3,t_VEC);
! 2097: p1=cgetg(r+1,small? t_VECSMALL: t_VEC); res[1]=(long)p1;
! 2098: p2=cgetg(r+1,small? t_VECSMALL: t_VEC); res[2]=(long)p2;
! 2099: if (d)
! 2100: {
! 2101: for (i=0,j=1; j<=n; j++)
! 2102: if (d[j]) { i++; p1[i]=d[j]; p2[i]=j; }
! 2103: free(d);
! 2104: qsort(p1+1,r,sizeof(long),(QSCOMP)pari_compare_long);
! 2105: }
! 2106: if (!small)
! 2107: for (i=1;i<=r;i++) { p1[i]=lstoi(p1[i]); p2[i]=lstoi(p2[i]); }
! 2108: return res;
! 2109: }
! 2110:
! 2111: GEN
! 2112: indexrank(GEN x) { return indexrank0(x,NULL,0); }
! 2113:
! 2114: GEN
! 2115: sindexrank(GEN x) { return indexrank0(x,NULL,1); }
! 2116:
! 2117: GEN
! 2118: FpM_sindexrank(GEN x, GEN p) { return indexrank0(x,p,1); }
! 2119:
! 2120: /*******************************************************************/
! 2121: /* */
! 2122: /* LINEAR ALGEBRA MODULO P */
! 2123: /* */
! 2124: /*******************************************************************/
! 2125:
! 2126: /*If p is NULL no reduction is performed.*/
! 2127: GEN
! 2128: FpM_mul(GEN x, GEN y, GEN p)
! 2129: {
! 2130: long i,j,l,lx=lg(x), ly=lg(y);
! 2131: GEN z;
! 2132: if (ly==1) return cgetg(ly,t_MAT);
! 2133: if (lx != lg(y[1])) err(operi,"* [mod p]",t_MAT,t_MAT);
! 2134: z=cgetg(ly,t_MAT);
! 2135: if (lx==1)
! 2136: {
! 2137: for (i=1; i<ly; i++) z[i]=lgetg(1,t_COL);
! 2138: return z;
! 2139: }
! 2140: l=lg(x[1]);
! 2141: for (j=1; j<ly; j++)
! 2142: {
! 2143: z[j] = lgetg(l,t_COL);
! 2144: for (i=1; i<l; i++)
! 2145: {
! 2146: ulong av;
! 2147: GEN p1,p2;
! 2148: int k;
! 2149: p1=gzero; av=avma;
! 2150: for (k=1; k<lx; k++)
! 2151: {
! 2152: p2=mulii(gcoeff(x,i,k),gcoeff(y,k,j));
! 2153: p1=addii(p1,p2);
! 2154: }
! 2155: coeff(z,i,j)=lpileupto(av,p?modii(p1,p):p1);
! 2156: }
! 2157: }
! 2158: return z;
! 2159: }
! 2160:
! 2161: /*If p is NULL no reduction is performed.*/
! 2162: GEN
! 2163: FpM_FpV_mul(GEN x, GEN y, GEN p)
! 2164: {
! 2165: long i,l,lx=lg(x), ly=lg(y);
! 2166: GEN z;
! 2167: if (lx != ly) err(operi,"* [mod p]",t_MAT,t_COL);
! 2168: z=cgetg(ly,t_COL);
! 2169: if (lx==1) return cgetg(1,t_COL);
! 2170: l=lg(x[1]);
! 2171: for (i=1; i<l; i++)
! 2172: {
! 2173: ulong av;
! 2174: GEN p1,p2;
! 2175: int k;
! 2176: p1=gzero; av=avma;
! 2177: for (k=1; k<lx; k++)
! 2178: {
! 2179: p2=mulii(gcoeff(x,i,k),(GEN)y[k]);
! 2180: p1=addii(p1,p2);
! 2181: }
! 2182: z[i]=lpileupto(av,p?modii(p1,p):p1);
! 2183: }
! 2184: return z;
! 2185: }
! 2186:
! 2187: static GEN
! 2188: u_FpM_ker(GEN x, GEN pp, long nontriv)
! 2189: {
! 2190: GEN y,c,d;
! 2191: long a,i,j,k,r,t,n,m,av0,tetpil, p = pp[2], piv;
! 2192:
! 2193: n = lg(x)-1;
! 2194: m=lg(x[1])-1; r=0; av0 = avma;
! 2195: x = dummycopy(x);
! 2196: for (i=1; i<=n; i++)
! 2197: {
! 2198: GEN p1 = (GEN)x[i];
! 2199: for (j=1; j<=m; j++) p1[j] = itos((GEN)p1[j]);
! 2200: }
! 2201: c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0;
! 2202: d=new_chunk(n+1);
! 2203: a=0; /* for gcc -Wall */
! 2204: for (k=1; k<=n; k++)
! 2205: {
! 2206: for (j=1; j<=m; j++)
! 2207: if (!c[j])
! 2208: {
! 2209: a = coeff(x,j,k) % p;
! 2210: if (a) break;
! 2211: }
! 2212: if (j>m)
! 2213: {
! 2214: if (nontriv) { avma=av0; return NULL; }
! 2215: r++; d[k]=0;
! 2216: }
! 2217: else
! 2218: {
! 2219: c[j]=k; d[k]=j;
! 2220: piv = - u_Fp_inv(a, p);
! 2221: coeff(x,j,k) = -1;
! 2222: for (i=k+1; i<=n; i++)
! 2223: coeff(x,j,i) = (piv * coeff(x,j,i)) % p;
! 2224: for (t=1; t<=m; t++)
! 2225: if (t!=j)
! 2226: {
! 2227: piv = coeff(x,t,k) % p;
! 2228: if (piv)
! 2229: {
! 2230: coeff(x,t,k) = 0;
! 2231: for (i=k+1; i<=n; i++) _u_Fp_addmul((GEN)x[i],t,j,piv,p);
! 2232: }
! 2233: }
! 2234: }
! 2235: }
! 2236:
! 2237: tetpil=avma; y=cgetg(r+1,t_MAT);
! 2238: for (j=k=1; j<=r; j++,k++)
! 2239: {
! 2240: GEN c = cgetg(n+1,t_COL);
! 2241:
! 2242: y[j]=(long)c; while (d[k]) k++;
! 2243: for (i=1; i<k; i++)
! 2244: if (d[i])
! 2245: {
! 2246: long a = coeff(x,d[i],k) % p;
! 2247: if (a < 0) a += p;
! 2248: c[i] = lstoi(a);
! 2249: }
! 2250: else
! 2251: c[i] = zero;
! 2252: c[k]=un; for (i=k+1; i<=n; i++) c[i]=zero;
! 2253: }
! 2254: return gerepile(av0,tetpil,y);
! 2255: }
! 2256:
! 2257: static GEN
! 2258: FpM_ker_i(GEN x, GEN p, long nontriv)
! 2259: {
! 2260: GEN y,c,d,piv,mun;
! 2261: long i,j,k,r,t,n,m,av0,av,lim,tetpil;
! 2262:
! 2263: if (typ(x)!=t_MAT) err(typeer,"FpM_ker");
! 2264: n=lg(x)-1; if (!n) return cgetg(1,t_MAT);
! 2265: if (OK_ULONG(p)) return u_FpM_ker(x, p, nontriv);
! 2266:
! 2267: m=lg(x[1])-1; r=0; av0 = avma;
! 2268: x=dummycopy(x); mun=negi(gun);
! 2269: c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0;
! 2270: d=new_chunk(n+1);
! 2271: av=avma; lim=stack_lim(av,1);
! 2272: for (k=1; k<=n; k++)
! 2273: {
! 2274: for (j=1; j<=m; j++)
! 2275: if (!c[j])
! 2276: {
! 2277: coeff(x,j,k) = lmodii(gcoeff(x,j,k), p);
! 2278: if (signe(coeff(x,j,k))) break;
! 2279: }
! 2280: if (j>m)
! 2281: {
! 2282: if (nontriv) { avma = av0; return NULL; }
! 2283: r++; d[k]=0;
! 2284: for(j=1; j<k; j++)
! 2285: if (d[j]) coeff(x,d[j],k) = lclone(gcoeff(x,d[j],k));
! 2286: }
! 2287: else
! 2288: {
! 2289: c[j]=k; d[k]=j; piv = negi(mpinvmod(gcoeff(x,j,k), p));
! 2290: coeff(x,j,k) = (long)mun;
! 2291: for (i=k+1; i<=n; i++)
! 2292: coeff(x,j,i) = lmodii(mulii(piv,gcoeff(x,j,i)), p);
! 2293: for (t=1; t<=m; t++)
! 2294: if (t!=j)
! 2295: {
! 2296: piv = modii(gcoeff(x,t,k), p);
! 2297: if (signe(piv))
! 2298: {
! 2299: coeff(x,t,k)=zero;
! 2300: for (i=k+1; i<=n; i++)
! 2301: coeff(x,t,i) = laddii(gcoeff(x,t,i),mulii(piv,gcoeff(x,j,i)));
! 2302: if (low_stack(lim, stack_lim(av,1)))
! 2303: gerepile_gauss_FpM_ker(x,p,m,n,k,t,av);
! 2304: }
! 2305: }
! 2306: }
! 2307: }
! 2308:
! 2309: tetpil=avma; y=cgetg(r+1,t_MAT);
! 2310: for (j=k=1; j<=r; j++,k++)
! 2311: {
! 2312: GEN c = cgetg(n+1,t_COL);
! 2313:
! 2314: y[j]=(long)c; while (d[k]) k++;
! 2315: for (i=1; i<k; i++)
! 2316: if (d[i])
! 2317: {
! 2318: GEN p1=gcoeff(x,d[i],k);
! 2319: c[i] = lmodii(p1, p); gunclone(p1);
! 2320: }
! 2321: else
! 2322: c[i] = zero;
! 2323: c[k]=un; for (i=k+1; i<=n; i++) c[i]=zero;
! 2324: }
! 2325: return gerepile(av0,tetpil,y);
! 2326: }
! 2327:
! 2328: static void
! 2329: FpM_gauss_pivot(GEN x, GEN p, GEN *dd, long *rr)
! 2330: {
! 2331: GEN c,d,piv;
! 2332: long i,j,k,r,t,n,m,av,lim;
! 2333:
! 2334: if (!p) return gauss_pivot(x,dd,rr);
! 2335: if (typ(x)!=t_MAT) err(typeer,"FpM_gauss_pivot");
! 2336: n=lg(x)-1; if (!n) { *dd=NULL; *rr=0; return; }
! 2337:
! 2338: m=lg(x[1])-1; r=0;
! 2339: x=dummycopy(x);
! 2340: c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0;
! 2341: d=(GEN)gpmalloc((n+1)*sizeof(long)); av=avma; lim=stack_lim(av,1);
! 2342: for (k=1; k<=n; k++)
! 2343: {
! 2344: for (j=1; j<=m; j++)
! 2345: if (!c[j])
! 2346: {
! 2347: coeff(x,j,k) = lmodii(gcoeff(x,j,k), p);
! 2348: if (signe(coeff(x,j,k))) break;
! 2349: }
! 2350: if (j>m) { r++; d[k]=0; }
! 2351: else
! 2352: {
! 2353: c[j]=k; d[k]=j; piv = negi(mpinvmod(gcoeff(x,j,k), p));
! 2354: for (i=k+1; i<=n; i++)
! 2355: coeff(x,j,i) = lmodii(mulii(piv,gcoeff(x,j,i)), p);
! 2356: for (t=1; t<=m; t++)
! 2357: if (!c[t]) /* no pivot on that line yet */
! 2358: {
! 2359: piv=gcoeff(x,t,k);
! 2360: if (signe(piv))
! 2361: {
! 2362: coeff(x,t,k)=zero;
! 2363: for (i=k+1; i<=n; i++)
! 2364: coeff(x,t,i) = laddii(gcoeff(x,t,i), mulii(piv,gcoeff(x,j,i)));
! 2365: if (low_stack(lim, stack_lim(av,1)))
! 2366: gerepile_gauss(x,m,n,k,t,av,j,c);
! 2367: }
! 2368: }
! 2369: for (i=k; i<=n; i++) coeff(x,j,i) = zero; /* dummy */
! 2370: }
! 2371: }
! 2372: *dd=d; *rr=r;
! 2373: }
! 2374:
! 2375: GEN
! 2376: FpM_ker(GEN x, GEN p)
! 2377: {
! 2378: return FpM_ker_i(x, p, 0);
! 2379: }
! 2380:
! 2381: int
! 2382: ker_trivial_mod_p(GEN x, GEN p)
! 2383: {
! 2384: return FpM_ker_i(x, p, 1)==NULL;
! 2385: }
! 2386:
! 2387: GEN
! 2388: FpM_image(GEN x, GEN p)
! 2389: {
! 2390: GEN d,y;
! 2391: long j,k,r, av = avma;
! 2392:
! 2393: FpM_gauss_pivot(x,p,&d,&r);
! 2394:
! 2395: /* r = dim ker(x) */
! 2396: if (!r) { avma=av; if (d) free(d); return gcopy(x); }
! 2397:
! 2398: /* r = dim Im(x) */
! 2399: r = lg(x)-1 - r; avma=av;
! 2400: y=cgetg(r+1,t_MAT);
! 2401: for (j=k=1; j<=r; k++)
! 2402: if (d[k]) y[j++] = lcopy((GEN)x[k]);
! 2403: free(d); return y;
! 2404: }
! 2405:
! 2406: static GEN
! 2407: sFpM_invimage(GEN mat, GEN y, GEN p)
! 2408: {
! 2409: long av=avma,i, nbcol = lg(mat);
! 2410: GEN col,p1 = cgetg(nbcol+1,t_MAT),res;
! 2411:
! 2412: if (nbcol==1) return NULL;
! 2413: if (lg(y) != lg(mat[1])) err(consister,"FpM_invimage");
! 2414:
! 2415: p1[nbcol] = (long)y;
! 2416: for (i=1; i<nbcol; i++) p1[i]=mat[i];
! 2417: p1 = FpM_ker(p1,p); i=lg(p1)-1;
! 2418: if (!i) return NULL;
! 2419:
! 2420: col = (GEN)p1[i]; p1 = (GEN) col[nbcol];
! 2421: if (gcmp0(p1)) return NULL;
! 2422:
! 2423: p1 = mpinvmod(negi(p1),p);
! 2424: setlg(col,nbcol);
! 2425: for(i=1;i<nbcol;i++)
! 2426: col[i]=lmulii((GEN)col[i],p1);
! 2427: res=cgetg(nbcol,t_COL);
! 2428: for(i=1;i<nbcol;i++)
! 2429: res[i]=lmodii((GEN)col[i],p);
! 2430: return gerepileupto(av,res);
! 2431: }
! 2432:
! 2433: /* Calcule l'image reciproque de v par m */
! 2434: GEN
! 2435: FpM_invimage(GEN m, GEN v, GEN p)
! 2436: {
! 2437: long av=avma,j,lv,tv=typ(v);
! 2438: GEN y,p1;
! 2439:
! 2440: if (typ(m)!=t_MAT) err(typeer,"inverseimage");
! 2441: if (tv==t_COL)
! 2442: {
! 2443: p1 = sFpM_invimage(m,v,p);
! 2444: if (p1) return p1;
! 2445: avma = av; return cgetg(1,t_MAT);
! 2446: }
! 2447: if (tv!=t_MAT) err(typeer,"inverseimage");
! 2448:
! 2449: lv=lg(v)-1; y=cgetg(lv+1,t_MAT);
! 2450: for (j=1; j<=lv; j++)
! 2451: {
! 2452: p1 = sFpM_invimage(m,(GEN)v[j],p);
! 2453: if (!p1) { avma = av; return cgetg(1,t_MAT); }
! 2454: y[j] = (long)p1;
! 2455: }
! 2456: return y;
! 2457: }
! 2458: /**************************************************************
! 2459: Rather stupid implementation of linear algebra in finite fields
! 2460: **************************************************************/
! 2461: static GEN
! 2462: Fq_add(GEN x, GEN y, GEN T/*unused*/, GEN p)
! 2463: {
! 2464: switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
! 2465: {
! 2466: case 0: return modii(addii(x,y),p);
! 2467: case 1: return FpX_Fp_add(x,y,p);
! 2468: case 2: return FpX_Fp_add(y,x,p);
! 2469: case 3: return FpX_add(x,y,p);
! 2470: }
! 2471: return NULL;
! 2472: }
! 2473: #if 0
! 2474: /* this function is really for FpV_roots_to_pol in polarit1.c
! 2475: * For consistency we write the code there.
! 2476: * To avoid having to remove static status, we rewrite it in polarit1.c
! 2477: */
! 2478: static GEN
! 2479: Fq_neg(GEN x, GEN T, GEN p)
! 2480: {
! 2481: switch(typ(x)==t_POL)
! 2482: {
! 2483: case 0: return signe(x)?subii(p,x):gzero;
! 2484: case 1: return FpX_neg(x,p);
! 2485: }
! 2486: return NULL;
! 2487: }
! 2488: #endif
! 2489:
! 2490: static GEN
! 2491: Fq_mul(GEN x, GEN y, GEN T, GEN p)
! 2492: {
! 2493: switch((typ(x)==t_POL)|((typ(y)==t_POL)<<1))
! 2494: {
! 2495: case 0: return modii(mulii(x,y),p);
! 2496: case 1: return FpX_Fp_mul(x,y,p);
! 2497: case 2: return FpX_Fp_mul(y,x,p);
! 2498: case 3: return FpXQ_mul(x,y,T,p);
! 2499: }
! 2500: return NULL;
! 2501: }
! 2502:
! 2503: static GEN
! 2504: Fq_neg_inv(GEN x, GEN T, GEN p)
! 2505: {
! 2506: switch(typ(x)==t_POL)
! 2507: {
! 2508: case 0: return mpinvmod(negi(x),p);
! 2509: case 1: return FpXQ_inv(FpX_neg(x,p),T,p);
! 2510: }
! 2511: return NULL;
! 2512: }
! 2513:
! 2514: static GEN
! 2515: Fq_res(GEN x, GEN T, GEN p)
! 2516: {
! 2517: ulong ltop=avma;
! 2518: switch(typ(x)==t_POL)
! 2519: {
! 2520: case 0: return modii(x,p);
! 2521: case 1: return gerepileupto(ltop,FpX_res(FpX_red(x,p),T,p));
! 2522: }
! 2523: return NULL;
! 2524: }
! 2525:
! 2526: static void
! 2527: Fq_gerepile_gauss_ker(GEN x, GEN T, GEN p, long m, long n, long k, long t,
! 2528: ulong av)
! 2529: {
! 2530: ulong tetpil = avma,A;
! 2531: long dec,u,i;
! 2532:
! 2533: if (DEBUGMEM > 1) err(warnmem,"gauss_pivot_ker. k=%ld, n=%ld",k,n);
! 2534: for (u=t+1; u<=m; u++)
! 2535: if (isonstack(coeff(x,u,k))) coeff(x,u,k) = (long) Fq_res(gcoeff(x,u,k),T,p);
! 2536: for (i=k+1; i<=n; i++)
! 2537: for (u=1; u<=m; u++)
! 2538: if (isonstack(coeff(x,u,i))) coeff(x,u,i) = (long) Fq_res(gcoeff(x,u,i),T,p);
! 2539:
! 2540: (void)gerepile(av,tetpil,NULL); dec = av-tetpil;
! 2541: for (u=t+1; u<=m; u++)
! 2542: {
! 2543: A=coeff(x,u,k);
! 2544: if (A<av && A>=bot) coeff(x,u,k)+=dec;
! 2545: }
! 2546: for (i=k+1; i<=n; i++)
! 2547: for (u=1; u<=m; u++)
! 2548: {
! 2549: A=coeff(x,u,i);
! 2550: if (A<av && A>=bot) coeff(x,u,i)+=dec;
! 2551: }
! 2552: }
! 2553:
! 2554: static GEN
! 2555: FqM_ker_i(GEN x, GEN T, GEN p, long nontriv)
! 2556: {
! 2557: GEN y,c,d,piv,mun;
! 2558: long i,j,k,r,t,n,m,av0,av,lim,tetpil;
! 2559:
! 2560: if (typ(x)!=t_MAT) err(typeer,"FpM_ker");
! 2561: n=lg(x)-1; if (!n) return cgetg(1,t_MAT);
! 2562:
! 2563: m=lg(x[1])-1; r=0; av0 = avma;
! 2564: x=dummycopy(x); mun=negi(gun);
! 2565: c=new_chunk(m+1); for (k=1; k<=m; k++) c[k]=0;
! 2566: d=new_chunk(n+1);
! 2567: av=avma; lim=stack_lim(av,1);
! 2568: for (k=1; k<=n; k++)
! 2569: {
! 2570: for (j=1; j<=m; j++)
! 2571: if (!c[j])
! 2572: {
! 2573: coeff(x,j,k) = (long) Fq_res(gcoeff(x,j,k), T, p);
! 2574: if (signe(coeff(x,j,k))) break;
! 2575: }
! 2576: if (j>m)
! 2577: {
! 2578: if (nontriv) { avma = av0; return NULL; }
! 2579: r++; d[k]=0;
! 2580: for(j=1; j<k; j++)
! 2581: if (d[j]) coeff(x,d[j],k) = lclone(gcoeff(x,d[j],k));
! 2582: }
! 2583: else
! 2584: {
! 2585: c[j]=k; d[k]=j; piv = Fq_neg_inv(gcoeff(x,j,k), T, p);
! 2586: coeff(x,j,k) = (long)mun;
! 2587: for (i=k+1; i<=n; i++)
! 2588: coeff(x,j,i) = (long) Fq_mul(piv,gcoeff(x,j,i), T, p);
! 2589: for (t=1; t<=m; t++)
! 2590: if (t!=j)
! 2591: {
! 2592: piv = Fq_res(gcoeff(x,t,k), T, p);
! 2593: if (signe(piv))
! 2594: {
! 2595: coeff(x,t,k)=zero;
! 2596: for (i=k+1; i<=n; i++)
! 2597: coeff(x,t,i) = (long)Fq_add(gcoeff(x,t,i),Fq_mul(piv,gcoeff(x,j,i), T, p), T, p);
! 2598: if (low_stack(lim, stack_lim(av,1)))
! 2599: Fq_gerepile_gauss_ker(x,T,p,m,n,k,t,av);
! 2600: }
! 2601: }
! 2602: }
! 2603: }
! 2604:
! 2605: tetpil=avma; y=cgetg(r+1,t_MAT);
! 2606: for (j=k=1; j<=r; j++,k++)
! 2607: {
! 2608: GEN c = cgetg(n+1,t_COL);
! 2609:
! 2610: y[j]=(long)c; while (d[k]) k++;
! 2611: for (i=1; i<k; i++)
! 2612: if (d[i])
! 2613: {
! 2614: GEN p1=gcoeff(x,d[i],k);
! 2615: c[i] = (long) Fq_res(p1, T, p); gunclone(p1);
! 2616: }
! 2617: else
! 2618: c[i] = zero;
! 2619: c[k]=un; for (i=k+1; i<=n; i++) c[i]=zero;
! 2620: }
! 2621: return gerepile(av0,tetpil,y);
! 2622: }
! 2623: GEN
! 2624: FqM_ker(GEN x, GEN T, GEN p)
! 2625: {
! 2626: return FqM_ker_i(x, T, p, 0);
! 2627: }
! 2628:
! 2629: /*******************************************************************/
! 2630: /* */
! 2631: /* EIGENVECTORS */
! 2632: /* (independent eigenvectors, sorted by increasing eigenvalue) */
! 2633: /* */
! 2634: /*******************************************************************/
! 2635:
! 2636: GEN
! 2637: eigen(GEN x, long prec)
! 2638: {
! 2639: GEN y,rr,p,ssesp,r1,r2,r3;
! 2640: long e,i,k,l,ly,ex, n = lg(x);
! 2641: ulong av = avma;
! 2642:
! 2643: if (typ(x)!=t_MAT) err(typeer,"eigen");
! 2644: if (n != 1 && n != lg(x[1])) err(mattype1,"eigen");
! 2645: if (n<=2) return gcopy(x);
! 2646:
! 2647: ex = 16 - bit_accuracy(prec);
! 2648: y=cgetg(n,t_MAT);
! 2649: p=caradj(x,0,NULL); rr=roots(p,prec);
! 2650: for (i=1; i<n; i++)
! 2651: {
! 2652: GEN p1 = (GEN)rr[i];
! 2653: if (!signe(p1[2]) || gexpo((GEN)p1[2]) - gexpo(p1) < ex) rr[i] = p1[1];
! 2654: }
! 2655: ly=1; k=2; r2=(GEN)rr[1];
! 2656: for(;;)
! 2657: {
! 2658: r3 = grndtoi(r2, &e); if (e < ex) r2 = r3;
! 2659: ssesp = ker0(x,r2); l = lg(ssesp);
! 2660: if (l == 1 || ly + (l-1) > n)
! 2661: err(talker2, "missing eigenspace. Compute the matrix to higher accuracy, then restart eigen at the current precision",NULL,NULL);
! 2662: for (i=1; i<l; ) y[ly++]=ssesp[i++]; /* done with this eigenspace */
! 2663:
! 2664: r1=r2; /* try to find a different eigenvalue */
! 2665: do
! 2666: {
! 2667: if (k == n || ly == n)
! 2668: {
! 2669: setlg(y,ly); /* x may not be diagonalizable */
! 2670: return gerepilecopy(av,y);
! 2671: }
! 2672: r2 = (GEN)rr[k++];
! 2673: r3 = gsub(r1,r2);
! 2674: }
! 2675: while (gcmp0(r3) || gexpo(r3) < ex);
! 2676: }
! 2677: }
! 2678:
! 2679: /*******************************************************************/
! 2680: /* */
! 2681: /* DETERMINANT */
! 2682: /* */
! 2683: /*******************************************************************/
! 2684:
! 2685: GEN
! 2686: det0(GEN a,long flag)
! 2687: {
! 2688: switch(flag)
! 2689: {
! 2690: case 0: return det(a);
! 2691: case 1: return det2(a);
! 2692: default: err(flagerr,"matdet");
! 2693: }
! 2694: return NULL; /* not reached */
! 2695: }
! 2696:
! 2697: /* Exact types: choose the first non-zero pivot. Otherwise: maximal pivot */
! 2698: static GEN
! 2699: det_simple_gauss(GEN a, long inexact)
! 2700: {
! 2701: long i,j,k,av,av1,s, nbco = lg(a)-1;
! 2702: GEN x,p;
! 2703:
! 2704: av=avma; s=1; x=gun; a=dummycopy(a);
! 2705: for (i=1; i<nbco; i++)
! 2706: {
! 2707: p=gcoeff(a,i,i); k=i;
! 2708: if (inexact)
! 2709: {
! 2710: long e, ex = gexpo(p);
! 2711: GEN p1;
! 2712:
! 2713: for (j=i+1; j<=nbco; j++)
! 2714: {
! 2715: e = gexpo(gcoeff(a,i,j));
! 2716: if (e > ex) { ex=e; k=j; }
! 2717: }
! 2718: p1 = gcoeff(a,i,k);
! 2719: if (gcmp0(p1)) return gerepilecopy(av, p1);
! 2720: }
! 2721: else if (gcmp0(p))
! 2722: {
! 2723: do k++; while(k<=nbco && gcmp0(gcoeff(a,i,k)));
! 2724: if (k>nbco) return gerepilecopy(av, p);
! 2725: }
! 2726: if (k != i)
! 2727: {
! 2728: swap(a[i],a[k]); s = -s;
! 2729: p = gcoeff(a,i,i);
! 2730: }
! 2731:
! 2732: x = gmul(x,p);
! 2733: for (k=i+1; k<=nbco; k++)
! 2734: {
! 2735: GEN m = gcoeff(a,i,k);
! 2736: if (!gcmp0(m))
! 2737: {
! 2738: m = gneg_i(gdiv(m,p));
! 2739: for (j=i+1; j<=nbco; j++)
! 2740: coeff(a,j,k) = ladd(gcoeff(a,j,k), gmul(m,gcoeff(a,j,i)));
! 2741: }
! 2742: }
! 2743: }
! 2744: if (s<0) x = gneg_i(x);
! 2745: av1=avma; return gerepile(av,av1,gmul(x,gcoeff(a,nbco,nbco)));
! 2746: }
! 2747:
! 2748: /* a has integer entries, N = P^n */
! 2749: GEN
! 2750: det_mod_P_n(GEN a, GEN N, GEN P)
! 2751: {
! 2752: long va,i,j,k,s, av = avma, nbco = lg(a)-1;
! 2753: GEN x,p;
! 2754:
! 2755: s=1; va=0; x=gun; a=dummycopy(a);
! 2756: p = NULL; /* for gcc -Wall */
! 2757: for (i=1; i<nbco; i++)
! 2758: {
! 2759: long fl = 0;
! 2760: for(;;)
! 2761: {
! 2762: for (k=i; k<=nbco; k++)
! 2763: {
! 2764: p=gcoeff(a,i,k);
! 2765: if (signe(p))
! 2766: {
! 2767: fl = 1;
! 2768: if (resii(p,P) != gzero) break;
! 2769: }
! 2770: }
! 2771: if (k <= nbco) break;
! 2772: va++; N = divii(N, P);
! 2773: if (!fl || is_pm1(N)) { avma=av; return gzero; }
! 2774: for (k=i; k<=nbco; k++) coeff(a,i,k) = ldivii(gcoeff(a,i,k), P);
! 2775: }
! 2776: if (k != i) { swap(a[i],a[k]); s = -s; }
! 2777:
! 2778: x = gmul(x,p); p = mpinvmod(p,N);
! 2779: for (k=i+1; k<=nbco; k++)
! 2780: {
! 2781: GEN m = resii(gcoeff(a,i,k), N);
! 2782: coeff(a,i,k) = zero;
! 2783: if (signe(m))
! 2784: {
! 2785: m = negi(resii(mulii(m,p), N));
! 2786: for (j=i+1; j<=nbco; j++)
! 2787: coeff(a,j,k) = lresii(addii(gcoeff(a,j,k),
! 2788: mulii(m,gcoeff(a,j,i))), N);
! 2789: }
! 2790: }
! 2791: }
! 2792: if (s<0) x = negi(x);
! 2793: x = resii(mulii(x,gcoeff(a,nbco,nbco)), N);
! 2794: return gerepileuptoint(av, mulii(x, gpowgs(P,va)));
! 2795: }
! 2796:
! 2797: GEN
! 2798: det2(GEN a)
! 2799: {
! 2800: long nbco = lg(a)-1;
! 2801: if (typ(a)!=t_MAT) err(mattype1,"det2");
! 2802: if (!nbco) return gun;
! 2803: if (nbco != lg(a[1])-1) err(mattype1,"det2");
! 2804: return det_simple_gauss(a,use_maximal_pivot(a));
! 2805: }
! 2806:
! 2807: static GEN
! 2808: mydiv(GEN x, GEN y)
! 2809: {
! 2810: long tx = typ(x), ty = typ(y);
! 2811: if (tx == ty && tx == t_POL && varn(x) == varn(y))
! 2812: return gdeuc(x,y);
! 2813: return gdiv(x,y);
! 2814: }
! 2815:
! 2816: /* determinant in a ring A: all computations are done within A
! 2817: * (Gauss-Bareiss algorithm) */
! 2818: GEN
! 2819: det(GEN a)
! 2820: {
! 2821: ulong av, lim;
! 2822: long nbco = lg(a)-1,i,j,k,s;
! 2823: GEN p,pprec;
! 2824:
! 2825: if (typ(a)!=t_MAT) err(mattype1,"det");
! 2826: if (!nbco) return gun;
! 2827: if (nbco != lg(a[1])-1) err(mattype1,"det");
! 2828: if (use_maximal_pivot(a)) return det_simple_gauss(a,1);
! 2829: if (DEBUGLEVEL > 7) timer2();
! 2830:
! 2831: av = avma; lim = stack_lim(av,2);
! 2832: a = dummycopy(a); s = 1;
! 2833: for (pprec=gun,i=1; i<nbco; i++,pprec=p)
! 2834: {
! 2835: GEN *ci, *ck, m, p1;
! 2836: int diveuc = (gcmp1(pprec)==0);
! 2837:
! 2838: p = gcoeff(a,i,i);
! 2839: if (gcmp0(p))
! 2840: {
! 2841: k=i+1; while (k<=nbco && gcmp0(gcoeff(a,i,k))) k++;
! 2842: if (k>nbco) return gerepilecopy(av, p);
! 2843: swap(a[k], a[i]); s = -s;
! 2844: p = gcoeff(a,i,i);
! 2845: }
! 2846: ci = (GEN*)a[i];
! 2847: for (k=i+1; k<=nbco; k++)
! 2848: {
! 2849: ck = (GEN*)a[k]; m = (GEN)ck[i];
! 2850: if (gcmp0(m))
! 2851: {
! 2852: if (gcmp1(p))
! 2853: {
! 2854: if (diveuc)
! 2855: a[k] = (long)mydiv((GEN)a[k], pprec);
! 2856: }
! 2857: else
! 2858: for (j=i+1; j<=nbco; j++)
! 2859: {
! 2860: p1 = gmul(p,ck[j]);
! 2861: if (diveuc) p1 = mydiv(p1,pprec);
! 2862: ck[j] = p1;
! 2863: }
! 2864: }
! 2865: else
! 2866: {
! 2867: m = gneg_i(m);
! 2868: for (j=i+1; j<=nbco; j++)
! 2869: {
! 2870: p1 = gadd(gmul(p,ck[j]), gmul(m,ci[j]));
! 2871: if (diveuc) p1 = mydiv(p1,pprec);
! 2872: ck[j] = p1;
! 2873: }
! 2874: }
! 2875: if (low_stack(lim,stack_lim(av,2)))
! 2876: {
! 2877: GEN *gptr[2]; gptr[0]=&a; gptr[1]=&pprec;
! 2878: if(DEBUGMEM>1) err(warnmem,"det. col = %ld",i);
! 2879: gerepilemany(av,gptr,2); p = gcoeff(a,i,i); ci = (GEN*)a[i];
! 2880: }
! 2881: }
! 2882: if (DEBUGLEVEL > 7) msgtimer("det, col %ld / %ld",i,nbco-1);
! 2883: }
! 2884: p = gcoeff(a,nbco,nbco);
! 2885: if (s < 0) p = gneg(p); else p = gcopy(p);
! 2886: return gerepileupto(av, p);
! 2887: }
! 2888:
! 2889: /* return a solution of congruence system sum M_{ i,j } X_j = Y_i mod D_i
! 2890: * If ptu1 != NULL, put in *ptu1 a Z-basis of the homogeneous system
! 2891: */
! 2892: static GEN
! 2893: gaussmoduloall(GEN M, GEN D, GEN Y, GEN *ptu1)
! 2894: {
! 2895: long n,m,i,j,lM,av=avma,tetpil;
! 2896: GEN p1,delta,H,U,u1,u2,x;
! 2897:
! 2898: if (typ(M)!=t_MAT) err(typeer,"gaussmodulo");
! 2899: lM = lg(M); m = lM-1;
! 2900: if (!m)
! 2901: {
! 2902: if ((typ(Y)!=t_INT && lg(Y)!=1)
! 2903: || (typ(D)!=t_INT && lg(D)!=1)) err(consister,"gaussmodulo");
! 2904: return gzero;
! 2905: }
! 2906: n = lg(M[1])-1;
! 2907: switch(typ(D))
! 2908: {
! 2909: case t_VEC:
! 2910: case t_COL: delta=diagonal(D); break;
! 2911: case t_INT: delta=gscalmat(D,n); break;
! 2912: default: err(typeer,"gaussmodulo");
! 2913: return NULL; /* not reached */
! 2914: }
! 2915: if (typ(Y) == t_INT)
! 2916: {
! 2917: p1 = cgetg(n+1,t_COL);
! 2918: for (i=1; i<=n; i++) p1[i]=(long)Y;
! 2919: Y = p1;
! 2920: }
! 2921: p1 = hnfall(concatsp(M,delta));
! 2922: H = (GEN)p1[1]; U = (GEN)p1[2];
! 2923: Y = gauss(H,Y);
! 2924: if (!gcmp1(denom(Y))) return gzero;
! 2925: u1 = cgetg(m+1,t_MAT);
! 2926: u2 = cgetg(n+1,t_MAT);
! 2927: for (j=1; j<=m; j++)
! 2928: {
! 2929: p1 = (GEN)U[j]; setlg(p1,lM);
! 2930: u1[j] = (long)p1;
! 2931: }
! 2932: U += m;
! 2933: for (j=1; j<=n; j++)
! 2934: {
! 2935: p1 = (GEN)U[j]; setlg(p1,lM);
! 2936: u2[j] = (long)p1;
! 2937: }
! 2938: x = gmul(u2,Y);
! 2939: tetpil=avma; x=lllreducemodmatrix(x,u1);
! 2940: if (!ptu1) x = gerepile(av,tetpil,x);
! 2941: else
! 2942: {
! 2943: GEN *gptr[2];
! 2944: *ptu1=gcopy(u1); gptr[0]=ptu1; gptr[1]=&x;
! 2945: gerepilemanysp(av,tetpil,gptr,2);
! 2946: }
! 2947: return x;
! 2948: }
! 2949:
! 2950: GEN
! 2951: matsolvemod0(GEN M, GEN D, GEN Y, long flag)
! 2952: {
! 2953: long av;
! 2954: GEN p1,y;
! 2955:
! 2956: if (!flag) return gaussmoduloall(M,D,Y,NULL);
! 2957:
! 2958: av=avma; y = cgetg(3,t_VEC);
! 2959: p1 = gaussmoduloall(M,D,Y, (GEN*)y+2);
! 2960: if (p1==gzero) { avma=av; return gzero; }
! 2961: y[1] = (long)p1; return y;
! 2962: }
! 2963:
! 2964: GEN
! 2965: gaussmodulo2(GEN M, GEN D, GEN Y)
! 2966: {
! 2967: return matsolvemod0(M,D,Y,1);
! 2968: }
! 2969:
! 2970: GEN
! 2971: gaussmodulo(GEN M, GEN D, GEN Y)
! 2972: {
! 2973: return matsolvemod0(M,D,Y,0);
! 2974: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to alglin1.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 / basemath