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