![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to gmpq.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine|
Return to gmpq.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 ! noro 1: #include "ca.h"
! 2: #include "gmp.h"
! 3: #include "base.h"
! 4: #include "inline.h"
! 5:
! 6: mpz_t ONEMPZ;
! 7: GZ ONEGZ;
! 8:
! 9: void isqrtgz(GZ a,GZ *r);
! 10: void bshiftgz(GZ a,int n,GZ *r);
! 11:
! 12: void *gc_realloc(void *p,size_t osize,size_t nsize)
! 13: {
! 14: return (void *)GC_realloc(p,nsize);
! 15: }
! 16:
! 17: void gc_free(void *p,size_t size)
! 18: {
! 19: GC_free(p);
! 20: }
! 21:
! 22: void init_gmpq()
! 23: {
! 24: mp_set_memory_functions(GC_malloc_atomic,gc_realloc,gc_free);
! 25:
! 26: mpz_init(ONEMPZ); mpz_set_ui(ONEMPZ,1); MPZTOGZ(ONEMPZ,ONEGZ);
! 27: }
! 28:
! 29: GZ utogz(unsigned int u)
! 30: {
! 31: mpz_t z;
! 32: GZ r;
! 33:
! 34: if ( !u ) return 0;
! 35: mpz_init(z); mpz_set_ui(z,u); MPZTOGZ(z,r); return r;
! 36: }
! 37:
! 38: GZ stogz(int s)
! 39: {
! 40: mpz_t z;
! 41: GZ r;
! 42:
! 43: if ( !s ) return 0;
! 44: mpz_init(z); mpz_set_si(z,s); MPZTOGZ(z,r); return r;
! 45: }
! 46:
! 47: GQ mpqtogzq(mpq_t a)
! 48: {
! 49: GZ z;
! 50: GQ q;
! 51:
! 52: if ( INTMPQ(a) ) {
! 53: MPZTOGZ(mpq_numref(a),z); return (GQ)z;
! 54: } else {
! 55: MPQTOGQ(a,q); return q;
! 56: }
! 57: }
! 58:
! 59: GZ ztogz(Q a)
! 60: {
! 61: mpz_t z;
! 62: mpq_t b;
! 63: N nm;
! 64: GZ s;
! 65:
! 66: if ( !a ) return 0;
! 67: nm = NM(a);
! 68: mpz_init(z);
! 69: mpz_import(z,PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));
! 70: if ( SGN(a)<0 ) mpz_neg(z,z);
! 71: MPZTOGZ(z,s);
! 72: return s;
! 73: }
! 74:
! 75: Q gztoz(GZ a)
! 76: {
! 77: N nm;
! 78: Q q;
! 79: int sgn;
! 80: size_t len;
! 81:
! 82: if ( !a ) return 0;
! 83: len = WORDSIZE_IN_N(BDY((GZ)a)); nm = NALLOC(len);
! 84: mpz_export(BD(nm),&len,-1,sizeof(int),0,0,BDY((GZ)a));
! 85: PL(nm) = len;
! 86: sgn = mpz_sgn(BDY((GZ)a)); NTOQ(nm,sgn,q);
! 87: return q;
! 88: }
! 89:
! 90: int n_bits_gz(GZ a)
! 91: {
! 92: return a ? mpz_sizeinbase(BDY(a),2) : 0;
! 93: }
! 94:
! 95: GQ qtogq(Q a)
! 96: {
! 97: mpz_t z;
! 98: mpq_t b;
! 99: N nm,dn;
! 100: GZ s;
! 101: GQ r;
! 102:
! 103: if ( !a ) return 0;
! 104: if ( INT(a) ) {
! 105: nm = NM(a);
! 106: mpz_init(z);
! 107: mpz_import(z,PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));
! 108: if ( SGN(a)<0 ) mpz_neg(z,z);
! 109: MPZTOGZ(z,s);
! 110: return (GQ)s;
! 111: } else {
! 112: nm = NM(a); dn = DN(a);
! 113: mpq_init(b);
! 114: mpz_import(mpq_numref(b),PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));
! 115: mpz_import(mpq_denref(b),PL(dn),-1,sizeof(BD(dn)[0]),0,0,BD(dn));
! 116: if ( SGN(a)<0 ) mpq_neg(b,b);
! 117: MPQTOGQ(b,r);
! 118: return r;
! 119: }
! 120: }
! 121:
! 122: Q gqtoq(GQ a)
! 123: {
! 124: N nm,dn;
! 125: Q q;
! 126: int sgn;
! 127: size_t len;
! 128:
! 129: if ( !a ) return 0;
! 130: if ( NID(a) == N_GZ ) {
! 131: len = WORDSIZE_IN_N(BDY((GZ)a)); nm = NALLOC(len);
! 132: mpz_export(BD(nm),&len,-1,sizeof(int),0,0,BDY((GZ)a));
! 133: PL(nm) = len;
! 134: sgn = mpz_sgn(BDY((GZ)a)); NTOQ(nm,sgn,q);
! 135: } else {
! 136: len = WORDSIZE_IN_N(mpq_numref(BDY(a))); nm = NALLOC(len);
! 137: mpz_export(BD(nm),&len,-1,sizeof(int),0,0,mpq_numref(BDY(a)));
! 138: PL(nm) = len;
! 139: len = WORDSIZE_IN_N(mpq_denref(BDY(a))); dn = NALLOC(len);
! 140: mpz_export(BD(dn),&len,-1,sizeof(int),0,0,mpq_denref(BDY(a)));
! 141: PL(dn) = len;
! 142: sgn = mpz_sgn(mpq_numref(BDY(a))); NDTOQ(nm,dn,sgn,q);
! 143: }
! 144: return q;
! 145: }
! 146:
! 147: P ptogp(P a)
! 148: {
! 149: DCP dc,dcr,dcr0;
! 150: P b;
! 151:
! 152: if ( !a ) return 0;
! 153: if ( NUM(a) ) return (P)qtogq((Q)a);
! 154: for ( dc = DC(a), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 155: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = (P)ptogp(COEF(dc));
! 156: }
! 157: NEXT(dcr) = 0; MKP(VR(a),dcr0,b);
! 158: return b;
! 159: }
! 160:
! 161: P gptop(P a)
! 162: {
! 163: DCP dc,dcr,dcr0;
! 164: P b;
! 165:
! 166: if ( !a ) return 0;
! 167: if ( NUM(a) ) b = (P)gqtoq((GQ)a);
! 168: else {
! 169: for ( dc = DC(a), dcr0 = 0; dc; dc = NEXT(dc) ) {
! 170: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc);
! 171: COEF(dcr) = (P)gptop(COEF(dc));
! 172: }
! 173: NEXT(dcr) = 0; MKP(VR(a),dcr0,b);
! 174: }
! 175: return b;
! 176: }
! 177:
! 178: void addgz(GZ n1,GZ n2,GZ *nr)
! 179: {
! 180: mpz_t t;
! 181: int s1,s2;
! 182:
! 183: if ( !n1 ) *nr = n2;
! 184: else if ( !n2 ) *nr = n1;
! 185: else {
! 186: mpz_init(t); mpz_add(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nr);
! 187: }
! 188: }
! 189:
! 190: void subgz(GZ n1,GZ n2,GZ *nr)
! 191: {
! 192: mpz_t t;
! 193:
! 194: if ( !n1 )
! 195: if ( !n2 )
! 196: *nr = 0;
! 197: else {
! 198: t[0] = BDY(n2)[0]; mpz_neg(t,t); MPZTOGZ(t,*nr);
! 199: }
! 200: else if ( !n2 )
! 201: *nr = n1;
! 202: else if ( n1 == n2 )
! 203: *nr = 0;
! 204: else {
! 205: mpz_init(t); mpz_sub(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nr);
! 206: }
! 207: }
! 208:
! 209: void mulgz(GZ n1,GZ n2,GZ *nr)
! 210: {
! 211: mpz_t t;
! 212:
! 213: if ( !n1 || !n2 ) *nr = 0;
! 214: else if ( UNIGZ(n1) ) *nr = n2;
! 215: else if ( UNIGZ(n2) ) *nr = n1;
! 216: else if ( MUNIGZ(n1) ) chsgngz(n2,nr);
! 217: else if ( MUNIGZ(n2) ) chsgngz(n1,nr);
! 218: else {
! 219: mpz_init(t); mpz_mul(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nr);
! 220: }
! 221: }
! 222:
! 223: void divgz(GZ n1,GZ n2,GZ *nq)
! 224: {
! 225: mpz_t t;
! 226: mpq_t a, b, q;
! 227:
! 228: if ( !n2 ) {
! 229: error("division by 0");
! 230: *nq = 0;
! 231: } else if ( !n1 )
! 232: *nq = 0;
! 233: else if ( n1 == n2 ) {
! 234: mpz_init(t); mpz_set_ui(t,1); MPZTOGZ(t,*nq);
! 235: } else {
! 236: MPZTOMPQ(BDY(n1),a); MPZTOMPQ(BDY(n2),b);
! 237: mpq_init(q); mpq_div(q,a,b); *nq = (GZ)mpqtogzq(q);
! 238: }
! 239: }
! 240:
! 241: void divqrgz(GZ n1,GZ n2,GZ *nq,GZ *nr)
! 242: {
! 243: mpz_t t, a, b, q, r;
! 244:
! 245: if ( !n2 ) {
! 246: error("division by 0");
! 247: *nq = 0; *nr = 0;
! 248: } else if ( !n1 ) {
! 249: *nq = 0; *nr = 0;
! 250: } else if ( n1 == n2 ) {
! 251: mpz_init(t); mpz_set_ui(t,1); MPZTOGZ(t,*nq); *nr = 0;
! 252: } else {
! 253: mpz_init(q); mpz_init(r);
! 254: mpz_fdiv_qr(q,r,BDY(n1),BDY(n2));
! 255: if ( !mpz_sgn(q) ) *nq = 0;
! 256: else MPZTOGZ(q,*nq);
! 257: if ( !mpz_sgn(r) ) *nr = 0;
! 258: else MPZTOGZ(r,*nr);
! 259: }
! 260: }
! 261:
! 262: void divsgz(GZ n1,GZ n2,GZ *nq)
! 263: {
! 264: mpz_t t;
! 265: mpq_t a, b, q;
! 266:
! 267: if ( !n2 ) {
! 268: error("division by 0");
! 269: *nq = 0;
! 270: } else if ( !n1 )
! 271: *nq = 0;
! 272: else if ( n1 == n2 ) {
! 273: mpz_init(t); mpz_set_ui(t,1); MPZTOGZ(t,*nq);
! 274: } else {
! 275: mpz_init(t); mpz_divexact(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nq);
! 276: }
! 277: }
! 278:
! 279: void chsgngz(GZ n,GZ *nr)
! 280: {
! 281: mpz_t t;
! 282:
! 283: if ( !n )
! 284: *nr = 0;
! 285: else {
! 286: t[0] = BDY(n)[0]; mpz_neg(t,t); MPZTOGZ(t,*nr);
! 287: }
! 288: }
! 289:
! 290: void pwrgz(GZ n1,Q n,GZ *nr)
! 291: {
! 292: mpq_t t,q;
! 293: mpz_t z;
! 294: GQ p,r;
! 295:
! 296: if ( !n || UNIGZ(n1) ) *nr = ONEGZ;
! 297: else if ( !n1 ) *nr = 0;
! 298: else if ( !INT(n) ) {
! 299: error("can't calculate fractional power."); *nr = 0;
! 300: } else if ( MUNIGZ(n1) ) *nr = BD(NM(n))[0]%2 ? n1 : ONEGZ;
! 301: else if ( PL(NM(n)) > 1 ) {
! 302: error("exponent too big."); *nr = 0;
! 303: } else if ( NID(n1)==N_GZ && SGN(n)>0 ) {
! 304: mpz_init(z); mpz_pow_ui(z,BDY(n1),QTOS(n)); MPZTOGZ(z,*nr);
! 305: } else {
! 306: MPZTOMPQ(BDY(n1),q); MPQTOGQ(q,r);
! 307: pwrgq(r,n,&p); *nr = (GZ)p;
! 308: }
! 309: }
! 310:
! 311: int cmpgz(GZ q1,GZ q2)
! 312: {
! 313: int sgn;
! 314:
! 315: if ( !q1 )
! 316: if ( !q2 )
! 317: return 0;
! 318: else
! 319: return -mpz_sgn(BDY(q2));
! 320: else if ( !q2 )
! 321: return mpz_sgn(BDY(q1));
! 322: else if ( (sgn = mpz_sgn(BDY(q1))) != mpz_sgn(BDY(q2)) )
! 323: return sgn;
! 324: else {
! 325: sgn = mpz_cmp(BDY(q1),BDY(q2));
! 326: if ( sgn > 0 ) return 1;
! 327: else if ( sgn < 0 ) return -1;
! 328: else return 0;
! 329: }
! 330: }
! 331:
! 332: void gcdgz(GZ n1,GZ n2,GZ *nq)
! 333: {
! 334: mpz_t t;
! 335:
! 336: if ( !n1 ) *nq = n2;
! 337: else if ( !n2 ) *nq = n1;
! 338: else {
! 339: mpz_init(t); mpz_gcd(t,BDY(n1),BDY(n2));
! 340: MPZTOGZ(t,*nq);
! 341: }
! 342: }
! 343:
! 344: void lcmgz(GZ n1,GZ n2,GZ *nq)
! 345: {
! 346: GZ g,t;
! 347:
! 348: if ( !n1 || !n2 ) *nq = 0;
! 349: else {
! 350: gcdgz(n1,n2,&g); divsgz(n1,g,&t);
! 351: mulgz(n2,t,nq);
! 352: }
! 353: }
! 354:
! 355: void gcdvgz(VECT v,GZ *q)
! 356: {
! 357: int n,i;
! 358: GZ *b;
! 359: GZ g,g1;
! 360:
! 361: n = v->len;
! 362: b = (GZ *)v->body;
! 363: g = b[0];
! 364: for ( i = 1; i < n; i++ ) {
! 365: gcdgz(g,b[i],&g1); g = g1;
! 366: }
! 367: *q = g;
! 368: }
! 369:
! 370: void gcdvgz_estimate(VECT v,GZ *q)
! 371: {
! 372: int n,m,i;
! 373: GZ s,t,u;
! 374: GZ *b;
! 375:
! 376: n = v->len;
! 377: b = (GZ *)v->body;
! 378: if ( n == 1 ) {
! 379: if ( mpz_sgn(BDY(b[0]))<0 ) chsgngz(b[0],q);
! 380: else *q = b[0];
! 381: }
! 382: m = n/2;
! 383: for ( i = 0, s = 0; i < m; i++ ) {
! 384: if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subgz(s,b[i],&u);
! 385: else addgz(s,b[i],&u);
! 386: s = u;
! 387: }
! 388: for ( i = 0, t = 0; i < n; i++ ) {
! 389: if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subgz(t,b[i],&u);
! 390: else addgz(t,b[i],&u);
! 391: t = u;
! 392: }
! 393: gcdgz(s,t,q);
! 394: }
! 395:
! 396: void addgq(GQ n1,GQ n2,GQ *nr)
! 397: {
! 398: mpq_t q1,q2,t;
! 399:
! 400: if ( !n1 ) *nr = n2;
! 401: else if ( !n2 ) *nr = n1;
! 402: else {
! 403: if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);
! 404: else q1[0] = BDY(n1)[0];
! 405: if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);
! 406: else q2[0] = BDY(n2)[0];
! 407: mpq_init(t); mpq_add(t,q1,q2); *nr = mpqtogzq(t);
! 408: }
! 409: }
! 410:
! 411: void subgq(GQ n1,GQ n2,GQ *nr)
! 412: {
! 413: mpq_t q1,q2,t;
! 414:
! 415: if ( !n1 )
! 416: if ( !n2 ) *nr = 0;
! 417: else {
! 418: if ( NID(n1) == N_GZ ) chsgngz((GZ)n1,(GZ *)nr);
! 419: else {
! 420: mpq_init(t); mpq_neg(t,BDY(n2)); MPQTOGQ(t,*nr);
! 421: }
! 422: }
! 423: else if ( !n2 ) *nr = n1;
! 424: else if ( n1 == n2 ) *nr = 0;
! 425: else {
! 426: if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);
! 427: else q1[0] = BDY(n1)[0];
! 428: if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);
! 429: else q2[0] = BDY(n2)[0];
! 430: mpq_init(t); mpq_sub(t,q1,q2); *nr = mpqtogzq(t);
! 431: }
! 432: }
! 433:
! 434: void mulgq(GQ n1,GQ n2,GQ *nr)
! 435: {
! 436: mpq_t t,q1,q2;
! 437:
! 438: if ( !n1 || !n2 ) *nr = 0;
! 439: else {
! 440: if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);
! 441: else q1[0] = BDY(n1)[0];
! 442: if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);
! 443: else q2[0] = BDY(n2)[0];
! 444: mpq_init(t); mpq_mul(t,q1,q2); *nr = mpqtogzq(t);
! 445: }
! 446: }
! 447:
! 448: void divgq(GQ n1,GQ n2,GQ *nq)
! 449: {
! 450: mpq_t t,q1,q2;
! 451:
! 452: if ( !n2 ) {
! 453: error("division by 0");
! 454: *nq = 0;
! 455: return;
! 456: } else if ( !n1 ) *nq = 0;
! 457: else if ( n1 == n2 ) *nq = (GQ)ONEGZ;
! 458: else {
! 459: if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);
! 460: else q1[0] = BDY(n1)[0];
! 461: if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);
! 462: else q2[0] = BDY(n2)[0];
! 463: mpq_init(t); mpq_div(t,q1,q2); *nq = mpqtogzq(t);
! 464: }
! 465: }
! 466:
! 467: void chsgngq(GQ n,GQ *nr)
! 468: {
! 469: mpq_t t;
! 470:
! 471: if ( !n ) *nr = 0;
! 472: else if ( NID(n) == N_GZ ) chsgngz((GZ)n,(GZ *)nr);
! 473: else {
! 474: mpq_init(t); mpq_neg(t,BDY(n)); MPQTOGQ(t,*nr);
! 475: }
! 476: }
! 477:
! 478: void pwrgq(GQ n1,Q n,GQ *nr)
! 479: {
! 480: int e;
! 481: mpz_t nm,dn;
! 482: mpq_t t;
! 483:
! 484: if ( !n || UNIGZ((GZ)n1) || UNIGQ(n1) ) *nr = (GQ)ONEGZ;
! 485: else if ( !n1 ) *nr = 0;
! 486: else if ( !INT(n) ) {
! 487: error("can't calculate fractional power."); *nr = 0;
! 488: } else if ( PL(NM(n)) > 1 ) {
! 489: error("exponent too big."); *nr = 0;
! 490: } else {
! 491: e = QTOS(n);
! 492: if ( e < 0 ) {
! 493: e = -e;
! 494: if ( NID(n1)==N_GZ ) {
! 495: nm[0] = ONEMPZ[0];
! 496: dn[0] = BDY((GZ)n1)[0];
! 497: } else {
! 498: nm[0] = mpq_denref(BDY(n1))[0]; dn[0] = mpq_numref(BDY(n1))[0];
! 499: }
! 500: } else {
! 501: if ( NID(n1)==N_GZ ) {
! 502: nm[0] = BDY((GZ)n1)[0]; dn[0] = ONEMPZ[0];
! 503: } else {
! 504: nm[0] = mpq_numref(BDY(n1))[0]; dn[0] = mpq_denref(BDY(n1))[0];
! 505: }
! 506: }
! 507: mpq_init(t);
! 508: mpz_pow_ui(mpq_numref(t),nm,e); mpz_pow_ui(mpq_denref(t),dn,e);
! 509: *nr = mpqtogzq(t);
! 510: }
! 511: }
! 512:
! 513: int cmpgq(GQ n1,GQ n2)
! 514: {
! 515: mpq_t q1,q2;
! 516: int sgn;
! 517:
! 518: if ( !n1 )
! 519: if ( !n2 ) return 0;
! 520: else return (NID(n2)==N_GZ) ? -mpz_sgn(BDY((GZ)n2)) : -mpq_sgn(BDY(n2));
! 521: if ( !n2 ) return (NID(n1)==N_GZ) ? mpz_sgn(BDY((GZ)n1)) : mpq_sgn(BDY(n1));
! 522: else if ( (sgn = mpq_sgn(BDY(n1))) != mpq_sgn(BDY(n2)) ) return sgn;
! 523: else {
! 524: if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);
! 525: else q1[0] = BDY(n1)[0];
! 526: if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);
! 527: else q2[0] = BDY(n2)[0];
! 528: sgn = mpq_cmp(q1,q2);
! 529: if ( sgn > 0 ) return 1;
! 530: else if ( sgn < 0 ) return -1;
! 531: else return 0;
! 532: }
! 533: }
! 534:
! 535: void mkgwc(int k,int l,GZ *t)
! 536: {
! 537: mpz_t a,b,q,nm,z,u;
! 538: int i,n;
! 539:
! 540: n = MIN(k,l);
! 541: mpz_init_set_ui(z,1);
! 542: mpz_init(u); mpz_set(u,z); MPZTOGZ(u,t[0]);
! 543: mpz_init(a); mpz_init(b); mpz_init(nm);
! 544: for ( i = 1; i <= n; i++ ) {
! 545: mpz_set_ui(a,k-i+1); mpz_set_ui(b,l-i+1); mpz_mul(nm,a,b);
! 546: mpz_mul(z,BDY(t[i-1]),nm); mpz_fdiv_q_ui(z,z,i);
! 547: mpz_init(u); mpz_set(u,z); MPZTOGZ(u,t[i]);
! 548: }
! 549: }
! 550:
! 551: void gz_lgp(P p,GZ *g,GZ *l);
! 552:
! 553: void gz_ptozp(P p,int sgn,GQ *c,P *pr)
! 554: {
! 555: GZ nm,dn;
! 556:
! 557: if ( !p ) {
! 558: *c = 0; *pr = 0;
! 559: } else {
! 560: gz_lgp(p,&nm,&dn);
! 561: divgz(nm,dn,(GZ *)c);
! 562: divsp(CO,p,(P)*c,pr);
! 563: }
! 564: }
! 565:
! 566: void gz_lgp(P p,GZ *g,GZ *l)
! 567: {
! 568: DCP dc;
! 569: GZ g1,g2,l1,l2,l3,l4;
! 570:
! 571: if ( NUM(p) ) {
! 572: if ( NID((GZ)p)==N_GZ ) {
! 573: MPZTOGZ(BDY((GZ)p),*g);
! 574: *l = ONEGZ;
! 575: } else {
! 576: MPZTOGZ(mpq_numref(BDY((GQ)p)),*g);
! 577: MPZTOGZ(mpq_denref(BDY((GQ)p)),*l);
! 578: }
! 579: } else {
! 580: dc = DC(p); gz_lgp(COEF(dc),g,l);
! 581: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
! 582: gz_lgp(COEF(dc),&g1,&l1); gcdgz(*g,g1,&g2); *g = g2;
! 583: gcdgz(*l,l1,&l2); mulgz(*l,l1,&l3); divgz(l3,l2,l);
! 584: }
! 585: }
! 586: }
! 587:
! 588: void gz_qltozl(GQ *w,int n,GZ *dvr)
! 589: {
! 590: GZ nm,dn;
! 591: GZ g,g1,l1,l2,l3;
! 592: GQ c;
! 593: int i;
! 594: struct oVECT v;
! 595:
! 596: for ( i = 0; i < n; i++ )
! 597: if ( w[i] && NID(w[i])==N_GQ )
! 598: break;
! 599: if ( i == n ) {
! 600: v.id = O_VECT; v.len = n; v.body = (pointer *)w;
! 601: gcdvgz(&v,dvr); return;
! 602: }
! 603: for ( i = 0; !w[i]; i++ );
! 604: c = w[i];
! 605: if ( NID(c)==N_GQ ) {
! 606: MPZTOGZ(mpq_numref(BDY(c)),nm); MPZTOGZ(mpq_denref(BDY(c)),dn);
! 607: } else {
! 608: MPZTOGZ(BDY((GZ)c),nm); dn = ONEGZ;
! 609: }
! 610: for ( i++; i < n; i++ ) {
! 611: c = w[i];
! 612: if ( !c ) continue;
! 613: if ( NID(c)==N_GQ ) {
! 614: MPZTOGZ(mpq_numref(BDY(c)),g1); MPZTOGZ(mpq_denref(BDY(c)),l1);
! 615: } else {
! 616: MPZTOGZ(BDY((GZ)c),g1); l1 = ONEGZ;
! 617: }
! 618: gcdgz(nm,g1,&g); nm = g;
! 619: gcdgz(dn,l1,&l2); mulgz(dn,l1,&l3); divgz(l3,l2,&dn);
! 620: }
! 621: divgz(nm,dn,dvr);
! 622: }
! 623:
! 624: int gz_bits(GQ q)
! 625: {
! 626: if ( !q ) return 0;
! 627: else if ( NID(q)==N_Q )
! 628: return n_bits(NM((Q)q))+(INT((Q)q)?0:n_bits(DN((Q)q)));
! 629: else if ( NID(q)==N_GZ ) return mpz_sizeinbase(BDY((GZ)q),2);
! 630: else
! 631: return mpz_sizeinbase(mpq_numref(BDY(q)),2)
! 632: + mpz_sizeinbase(mpq_denref(BDY(q)),2);
! 633: }
! 634:
! 635: int gzp_mag(P p)
! 636: {
! 637: int s;
! 638: DCP dc;
! 639:
! 640: if ( !p ) return 0;
! 641: else if ( OID(p) == O_N ) return gz_bits((GQ)p);
! 642: else {
! 643: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) s += gzp_mag(COEF(dc));
! 644: return s;
! 645: }
! 646: }
! 647:
! 648: void makesubstgz(VL v,NODE *s)
! 649: {
! 650: NODE r,r0;
! 651: GZ q;
! 652: unsigned int n;
! 653:
! 654: for ( r0 = 0; v; v = NEXT(v) ) {
! 655: NEXTNODE(r0,r); BDY(r) = (pointer)v->v;
! 656: #if defined(_PA_RISC1_1)
! 657: n = mrand48()&BMASK; q = utogz(n);
! 658: #else
! 659: n = random(); q = utogz(n);
! 660: #endif
! 661: NEXTNODE(r0,r); BDY(r) = (pointer)q;
! 662: }
! 663: if ( r0 ) NEXT(r) = 0;
! 664: *s = r0;
! 665: }
! 666:
! 667: unsigned int remgq(GQ a,unsigned int mod)
! 668: {
! 669: unsigned int c,nm,dn;
! 670: mpz_t r;
! 671:
! 672: if ( !a ) return 0;
! 673: else if ( NID(a)==N_GZ ) {
! 674: mpz_init(r);
! 675: c = mpz_fdiv_r_ui(r,BDY((GZ)a),mod);
! 676: } else {
! 677: mpz_init(r);
! 678: nm = mpz_fdiv_r_ui(r,mpq_numref(BDY(a)),mod);
! 679: dn = mpz_fdiv_r_ui(r,mpq_denref(BDY(a)),mod);
! 680: dn = invm(dn,mod);
! 681: DMAR(nm,dn,0,mod,c);
! 682: }
! 683: return c;
! 684: }
! 685:
! 686: extern int DP_Print;
! 687:
! 688: #define GZ_F4_INTRAT_PERIOD 8
! 689:
! 690: int gz_generic_gauss_elim(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp)
! 691: {
! 692: int **wmat;
! 693: GZ **bmat,**tmat,*bmi,*tmi;
! 694: GZ q,m1,m2,m3,s,u;
! 695: int *wmi,*colstat,*wcolstat,*rind,*cind;
! 696: int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv;
! 697: MAT r,crmat;
! 698: int ret;
! 699:
! 700: bmat = (GZ **)mat->body;
! 701: row = mat->row; col = mat->col;
! 702: wmat = (int **)almat(row,col);
! 703: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
! 704: wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
! 705: for ( ind = 0; ; ind++ ) {
! 706: if ( DP_Print ) {
! 707: fprintf(asir_out,"."); fflush(asir_out);
! 708: }
! 709: md = get_lprime(ind);
! 710: for ( i = 0; i < row; i++ )
! 711: for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
! 712: wmi[j] = remgq((GQ)bmi[j],md);
! 713: rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);
! 714: if ( !ind ) {
! 715: RESET:
! 716: m1 = utogz(md);
! 717: rank0 = rank;
! 718: bcopy(wcolstat,colstat,col*sizeof(int));
! 719: MKMAT(crmat,rank,col-rank);
! 720: MKMAT(r,rank,col-rank); *nm = r;
! 721: tmat = (GZ **)crmat->body;
! 722: for ( i = 0; i < rank; i++ )
! 723: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
! 724: if ( !colstat[j] ) tmi[k++] = utogz(wmi[j]);
! 725: } else {
! 726: if ( rank < rank0 ) {
! 727: if ( DP_Print ) {
! 728: fprintf(asir_out,"lower rank matrix; continuing...\n");
! 729: fflush(asir_out);
! 730: }
! 731: continue;
! 732: } else if ( rank > rank0 ) {
! 733: if ( DP_Print ) {
! 734: fprintf(asir_out,"higher rank matrix; resetting...\n");
! 735: fflush(asir_out);
! 736: }
! 737: goto RESET;
! 738: } else {
! 739: for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
! 740: if ( j < col ) {
! 741: if ( DP_Print ) {
! 742: fprintf(asir_out,"inconsitent colstat; resetting...\n");
! 743: fflush(asir_out);
! 744: }
! 745: goto RESET;
! 746: }
! 747: }
! 748:
! 749: inv = invm(remgq((GQ)m1,md),md);
! 750: m2 = utogz(md); mulgz(m1,m2,&m3);
! 751: for ( i = 0; i < rank; i++ )
! 752: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
! 753: if ( !colstat[j] ) {
! 754: if ( tmi[k] ) {
! 755: /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
! 756: t = remgq((GQ)tmi[k],md);
! 757: if ( wmi[j] >= t ) t = wmi[j]-t;
! 758: else t = md-(t-wmi[j]);
! 759: DMAR(t,inv,0,md,t1)
! 760: u = utogz(t1); mulgz(m1,u,&s);
! 761: addgz(tmi[k],s,&u); tmi[k] = u;
! 762: } else if ( wmi[j] ) {
! 763: /* f3 = m1*(m1 mod m2)^(-1)*f2 */
! 764: DMAR(wmi[j],inv,0,md,t)
! 765: u = utogz(t); mulgz(m1,u,&s); tmi[k] = s;
! 766: }
! 767: k++;
! 768: }
! 769: m1 = m3;
! 770: if ( ind % GZ_F4_INTRAT_PERIOD )
! 771: ret = 0;
! 772: else
! 773: ret = gz_intmtoratm(crmat,m1,*nm,dn);
! 774: if ( ret ) {
! 775: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
! 776: *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
! 777: for ( j = k = l = 0; j < col; j++ )
! 778: if ( colstat[j] ) rind[k++] = j;
! 779: else cind[l++] = j;
! 780: if ( gz_gensolve_check(mat,*nm,*dn,rind,cind) )
! 781: return rank;
! 782: }
! 783: }
! 784: }
! 785: }
! 786:
! 787: int gz_generic_gauss_elim2(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp)
! 788: {
! 789:
! 790: MAT full;
! 791: GZ **bmat,**b;
! 792: GZ *bmi;
! 793: GZ dn0;
! 794: int row,col,md,i,j,rank,ret;
! 795: int **wmat;
! 796: int *wmi;
! 797: int *colstat,*rowstat;
! 798:
! 799: bmat = (GZ **)mat->body;
! 800: row = mat->row; col = mat->col;
! 801: wmat = (int **)almat(row,col);
! 802: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
! 803: rowstat = (int *)MALLOC_ATOMIC(row*sizeof(int));
! 804: /* XXX */
! 805: md = get_lprime(0);
! 806: for ( i = 0; i < row; i++ )
! 807: for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
! 808: wmi[j] = remgq((GQ)bmi[j],md);
! 809: rank = generic_gauss_elim_mod2(wmat,row,col,md,colstat,rowstat);
! 810: b = (GZ **)MALLOC(rank*sizeof(GZ));
! 811: for ( i = 0; i < rank; i++ ) b[i] = bmat[rowstat[i]];
! 812: NEWMAT(full); full->row = rank; full->col = col; full->body = (pointer **)b;
! 813: ret = gz_generic_gauss_elim_full(full,nm,dn,rindp,cindp);
! 814: if ( !ret ) {
! 815: rank = gz_generic_gauss_elim(mat,nm,&dn0,rindp,cindp);
! 816: for ( i = 0; i < rank; i++ ) dn[i] = dn0;
! 817: }
! 818: return rank;
! 819: }
! 820:
! 821: int gz_generic_gauss_elim_full(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp)
! 822: {
! 823: int **wmat;
! 824: int *stat;
! 825: GZ **bmat,**tmat,*bmi,*tmi,*ri;
! 826: GZ q,m1,m2,m3,s,u;
! 827: int *wmi,*colstat,*wcolstat,*rind,*cind;
! 828: int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv,h;
! 829: MAT r,crmat;
! 830: int ret,initialized,done;
! 831:
! 832: initialized = 0;
! 833: bmat = (GZ **)mat->body;
! 834: row = mat->row; col = mat->col;
! 835: wmat = (int **)almat(row,col);
! 836: stat = (int *)MALLOC_ATOMIC(row*sizeof(int));
! 837: for ( i = 0; i < row; i++ ) stat[i] = 0;
! 838: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
! 839: wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
! 840: for ( ind = 0; ; ind++ ) {
! 841: if ( DP_Print ) {
! 842: fprintf(asir_out,"."); fflush(asir_out);
! 843: }
! 844: md = get_lprime(ind);
! 845: for ( i = 0; i < row; i++ )
! 846: for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
! 847: wmi[j] = remgq((GQ)bmi[j],md);
! 848: rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);
! 849: if ( rank < row ) continue;
! 850: if ( !initialized ) {
! 851: m1 = utogz(md);
! 852: bcopy(wcolstat,colstat,col*sizeof(int));
! 853: MKMAT(crmat,row,col-row);
! 854: MKMAT(r,row,col-row); *nm = r;
! 855: tmat = (GZ **)crmat->body;
! 856: for ( i = 0; i < row; i++ )
! 857: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
! 858: if ( !colstat[j] ) tmi[k++] = utogz(wmi[j]);
! 859: initialized = 1;
! 860: } else {
! 861: for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
! 862: if ( j < col ) continue;
! 863:
! 864: inv = invm(remgq((GQ)m1,md),md);
! 865: m2 = utogz(md); mulgz(m1,m2,&m3);
! 866: for ( i = 0; i < row; i++ )
! 867: switch ( stat[i] ) {
! 868: case 1:
! 869: /* consistency check */
! 870: ri = (GZ *)BDY(r)[i]; wmi = wmat[i];
! 871: for ( j = 0; j < col; j++ ) if ( colstat[j] ) break;
! 872: h = md-remgq((GQ)dn[i],md);
! 873: for ( j++, k = 0; j < col; j++ )
! 874: if ( !colstat[j] ) {
! 875: t = remgq((GQ)ri[k],md);
! 876: DMAR(wmi[i],h,t,md,t1);
! 877: if ( t1 ) break;
! 878: }
! 879: if ( j == col ) { stat[i]++; break; }
! 880: else {
! 881: /* fall to the case 0 */
! 882: stat[i] = 0;
! 883: }
! 884: case 0:
! 885: tmi = tmat[i]; wmi = wmat[i];
! 886: for ( j = k = 0; j < col; j++ )
! 887: if ( !colstat[j] ) {
! 888: if ( tmi[k] ) {
! 889: /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
! 890: t = remgq((GQ)tmi[k],md);
! 891: if ( wmi[j] >= t ) t = wmi[j]-t;
! 892: else t = md-(t-wmi[j]);
! 893: DMAR(t,inv,0,md,t1)
! 894: u = utogz(t1); mulgz(m1,u,&s);
! 895: addgz(tmi[k],s,&u); tmi[k] = u;
! 896: } else if ( wmi[j] ) {
! 897: /* f3 = m1*(m1 mod m2)^(-1)*f2 */
! 898: DMAR(wmi[j],inv,0,md,t)
! 899: u = utogz(t); mulgz(m1,u,&s); tmi[k] = s;
! 900: }
! 901: k++;
! 902: }
! 903: break;
! 904: case 2: default:
! 905: break;
! 906: }
! 907: m1 = m3;
! 908: if ( ind % 4 )
! 909: ret = 0;
! 910: else
! 911: ret = gz_intmtoratm2(crmat,m1,*nm,dn,stat);
! 912: if ( ret ) {
! 913: *rindp = rind = (int *)MALLOC_ATOMIC(row*sizeof(int));
! 914: *cindp = cind = (int *)MALLOC_ATOMIC((col-row)*sizeof(int));
! 915: for ( j = k = l = 0; j < col; j++ )
! 916: if ( colstat[j] ) rind[k++] = j;
! 917: else cind[l++] = j;
! 918: return gz_gensolve_check2(mat,*nm,dn,rind,cind);
! 919: }
! 920: }
! 921: }
! 922: }
! 923:
! 924: int gz_generic_gauss_elim_direct(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp){
! 925: GZ **bmat,*s;
! 926: GZ u,v,w,x,d,t,y;
! 927: int row,col,i,j,k,l,m,rank;
! 928: int *colstat,*colpos,*cind;
! 929: MAT r,in;
! 930:
! 931: row = mat->row; col = mat->col;
! 932: MKMAT(in,row,col);
! 933: for ( i = 0; i < row; i++ )
! 934: for ( j = 0; j < col; j++ ) in->body[i][j] = mat->body[i][j];
! 935: bmat = (GZ **)in->body;
! 936: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
! 937: *rindp = colpos = (int *)MALLOC_ATOMIC(row*sizeof(int));
! 938: for ( j = 0, rank = 0, d = ONEGZ; j < col; j++ ) {
! 939: for ( i = rank; i < row && !bmat[i][j]; i++ );
! 940: if ( i == row ) { colstat[j] = 0; continue; }
! 941: else { colstat[j] = 1; colpos[rank] = j; }
! 942: if ( i != rank )
! 943: for ( k = j; k < col; k++ ) {
! 944: t = bmat[i][k]; bmat[i][k] = bmat[rank][k]; bmat[rank][k] = t;
! 945: }
! 946: for ( i = rank+1, v = bmat[rank][j]; i < row; i++ )
! 947: for ( k = j, u = bmat[i][j]; k < col; k++ ) {
! 948: mulgz(bmat[i][k],v,&w); mulgz(bmat[rank][k],u,&x);
! 949: subgz(w,x,&y); divsgz(y,d,&bmat[i][k]);
! 950: }
! 951: d = v; rank++;
! 952: }
! 953: *dn = d;
! 954: s = (GZ *)MALLOC(col*sizeof(GZ));
! 955: for ( i = rank-1; i >= 0; i-- ) {
! 956: for ( k = colpos[i]; k < col; k++ ) mulgz(bmat[i][k],d,&s[k]);
! 957: for ( m = rank-1; m > i; m-- ) {
! 958: for ( k = colpos[m], u = bmat[i][k]; k < col; k++ ) {
! 959: mulgz(bmat[m][k],u,&w); subgz(s[k],w,&x); s[k] = x;
! 960: }
! 961: }
! 962: for ( k = colpos[i], u = bmat[i][k]; k < col; k++ )
! 963: divgz(s[k],u,&bmat[i][k]);
! 964: }
! 965: *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
! 966: MKMAT(r,rank,col-rank); *nm = r;
! 967: for ( j = 0, k = 0; j < col; j++ )
! 968: if ( !colstat[j] ) {
! 969: cind[k] = j;
! 970: for ( i = 0; i < rank; i++ ) r->body[i][k] = bmat[i][j];
! 971: k++;
! 972: }
! 973: return rank;
! 974: }
! 975:
! 976: int gz_intmtoratm(MAT mat,GZ md,MAT nm,GZ *dn)
! 977: {
! 978: GZ t,s,b,dn0,dn1,nm1,q,u,unm,udn,dmy;
! 979: int i,j,k,l,row,col,sgn,ret;
! 980: GZ **rmat,**tmat,*tmi,*nmk;
! 981:
! 982: if ( UNIGZ(md) )
! 983: return 0;
! 984: row = mat->row; col = mat->col;
! 985: bshiftgz(md,1,&t);
! 986: isqrtgz(t,&s);
! 987: bshiftgz(s,64,&b);
! 988: if ( !b ) b = ONEGZ;
! 989: dn0 = ONEGZ;
! 990: tmat = (GZ **)mat->body;
! 991: rmat = (GZ **)nm->body;
! 992: for ( i = 0; i < row; i++ )
! 993: for ( j = 0, tmi = tmat[i]; j < col; j++ )
! 994: if ( tmi[j] ) {
! 995: mulgz(tmi[j],dn0,&s);
! 996: divqrgz(s,md,&dmy,&u);
! 997: ret = gz_inttorat(u,md,b,&sgn,&unm,&udn);
! 998: if ( !ret ) return 0;
! 999: else {
! 1000: if ( sgn < 0 ) chsgngz(unm,&nm1);
! 1001: else nm1 = unm;
! 1002: dn1 = udn;
! 1003: if ( !UNIGZ(dn1) ) {
! 1004: for ( k = 0; k < i; k++ )
! 1005: for ( l = 0, nmk = rmat[k]; l < col; l++ ) {
! 1006: mulgz(nmk[l],dn1,&q); nmk[l] = q;
! 1007: }
! 1008: for ( l = 0, nmk = rmat[i]; l < j; l++ ) {
! 1009: mulgz(nmk[l],dn1,&q); nmk[l] = q;
! 1010: }
! 1011: }
! 1012: rmat[i][j] = nm1;
! 1013: mulgz(dn0,dn1,&q); dn0 = q;
! 1014: }
! 1015: }
! 1016: *dn = dn0;
! 1017: return 1;
! 1018: }
! 1019:
! 1020: int gz_intmtoratm2(MAT mat,GZ md,MAT nm,GZ *dn,int *stat)
! 1021: {
! 1022: int row,col,i,j,ret;
! 1023: GZ dn0,dn1,t,s,b;
! 1024: GZ *w,*tmi;
! 1025: GZ **tmat;
! 1026:
! 1027: bshiftgz(md,1,&t);
! 1028: isqrtgz(t,&s);
! 1029: bshiftgz(s,64,&b);
! 1030: tmat = (GZ **)mat->body;
! 1031: if ( UNIGZ(md) ) return 0;
! 1032: row = mat->row; col = mat->col;
! 1033: dn0 = ONEGZ;
! 1034: for ( i = 0; i < row; i++ )
! 1035: if ( cmpgz(dn[i],dn0) > 0 ) dn0 = dn[i];
! 1036: w = (GZ *)MALLOC(col*sizeof(GZ));
! 1037: for ( i = 0; i < row; i++ )
! 1038: if ( stat[i] == 0 ) {
! 1039: for ( j = 0, tmi = tmat[i]; j < col; j++ )
! 1040: mulgz(tmi[j],dn0,&w[j]);
! 1041: ret = gz_intvtoratv(w,col,md,b,BDY(nm)[i],&dn[i]);
! 1042: if ( ret ) {
! 1043: stat[i] = 1;
! 1044: mulgz(dn0,dn[i],&t); dn[i] = t; dn0 = t;
! 1045: }
! 1046: }
! 1047: for ( i = 0; i < row; i++ ) if ( !stat[i] ) break;
! 1048: if ( i == row ) return 1;
! 1049: else return 0;
! 1050: }
! 1051:
! 1052: int gz_intvtoratv(GZ *v,int n,GZ md,GZ b,GZ *nm,GZ *dn)
! 1053: {
! 1054: GZ dn0,dn1,q,s,u,nm1,unm,udn,dmy;
! 1055: GZ *nmk;
! 1056: int j,l,col,ret,sgn;
! 1057:
! 1058: for ( j = 0; j < n; j++ ) nm[j] = 0;
! 1059: dn0 = ONEGZ;
! 1060: for ( j = 0; j < n; j++ ) {
! 1061: if ( !v[j] ) continue;
! 1062: mulgz(v[j],dn0,&s);
! 1063: divqrgz(s,md,&dmy,&u);
! 1064: ret = gz_inttorat(u,md,b,&sgn,&unm,&udn);
! 1065: if ( !ret ) return 0;
! 1066: if ( sgn < 0 ) chsgngz(unm,&nm1);
! 1067: else nm1 = unm;
! 1068: dn1 = udn;
! 1069: if ( !UNIGZ(dn1) )
! 1070: for ( l = 0; l < j; l++ ) {
! 1071: mulgz(nm[l],dn1,&q); nm[l] = q;
! 1072: }
! 1073: nm[j] = nm1;
! 1074: mulgz(dn0,dn1,&q); dn0 = q;
! 1075: }
! 1076: *dn = dn0;
! 1077: return 1;
! 1078: }
! 1079:
! 1080: /* assuming 0 < c < m */
! 1081:
! 1082: int gz_inttorat(GZ c,GZ m,GZ b,int *sgnp,GZ *nmp,GZ *dnp)
! 1083: {
! 1084: GZ qq,t,u1,v1,r1;
! 1085: GZ q,u2,v2,r2;
! 1086:
! 1087: u1 = 0; v1 = ONEGZ; u2 = m; v2 = c;
! 1088: while ( cmpgz(v2,b) >= 0 ) {
! 1089: divqrgz(u2,v2,&q,&r2); u2 = v2; v2 = r2;
! 1090: mulgz(q,v1,&t); subgz(u1,t,&r1); u1 = v1; v1 = r1;
! 1091: }
! 1092: if ( cmpgz(v1,b) >= 0 ) return 0;
! 1093: else {
! 1094: *nmp = v2;
! 1095: if ( mpz_sgn(BDY(v1))<0 ) {
! 1096: *sgnp = -1; chsgngz(v1,dnp);
! 1097: } else {
! 1098: *sgnp = 1; *dnp = v1;
! 1099: }
! 1100: return 1;
! 1101: }
! 1102: }
! 1103:
! 1104: extern int f4_nocheck;
! 1105:
! 1106: int gz_gensolve_check(MAT mat,MAT nm,GZ dn,int *rind,int *cind)
! 1107: {
! 1108: int row,col,rank,clen,i,j,k,l;
! 1109: GZ s,t;
! 1110: GZ *w;
! 1111: GZ *mati,*nmk;
! 1112:
! 1113: if ( f4_nocheck ) return 1;
! 1114: row = mat->row; col = mat->col; rank = nm->row; clen = nm->col;
! 1115: w = (GZ *)MALLOC(clen*sizeof(GZ));
! 1116: for ( i = 0; i < row; i++ ) {
! 1117: mati = (GZ *)mat->body[i];
! 1118: bzero(w,clen*sizeof(GZ));
! 1119: for ( k = 0; k < rank; k++ )
! 1120: for ( l = 0, nmk = (GZ *)nm->body[k]; l < clen; l++ ) {
! 1121: mulgz(mati[rind[k]],nmk[l],&t); addgz(w[l],t,&s); w[l] = s;
! 1122: }
! 1123: for ( j = 0; j < clen; j++ ) {
! 1124: mulgz(dn,mati[cind[j]],&t);
! 1125: if ( cmpgz(w[j],t) ) break;
! 1126: }
! 1127: if ( j != clen ) break;
! 1128: }
! 1129: if ( i != row ) return 0;
! 1130: else return 1;
! 1131: }
! 1132:
! 1133: int gz_gensolve_check2(MAT mat,MAT nm,GZ *dn,int *rind,int *cind)
! 1134: {
! 1135: int row,col,rank,clen,i,j,k,l;
! 1136: GZ s,t,u,d;
! 1137: GZ *w,*m;
! 1138: GZ *mati,*nmk;
! 1139:
! 1140: if ( f4_nocheck ) return 1;
! 1141: row = mat->row; col = mat->col; rank = nm->row; clen = nm->col;
! 1142: w = (GZ *)MALLOC(clen*sizeof(GZ));
! 1143: m = (GZ *)MALLOC(clen*sizeof(GZ));
! 1144: for ( d = dn[0], i = 1; i < rank; i++ ) {
! 1145: lcmgz(d,dn[i],&t); d = t;
! 1146: }
! 1147: for ( i = 0; i < rank; i++ ) divsgz(d,dn[i],&m[i]);
! 1148: for ( i = 0; i < row; i++ ) {
! 1149: mati = (GZ *)mat->body[i];
! 1150: bzero(w,clen*sizeof(GZ));
! 1151: for ( k = 0; k < rank; k++ ) {
! 1152: mulgz(mati[rind[k]],m[k],&u);
! 1153: for ( l = 0, nmk = (GZ *)nm->body[k]; l < clen; l++ ) {
! 1154: mulgz(u,nmk[l],&t); addgz(w[l],t,&s); w[l] = s;
! 1155: }
! 1156: }
! 1157: for ( j = 0; j < clen; j++ ) {
! 1158: mulgz(d,mati[cind[j]],&t);
! 1159: if ( cmpgz(w[j],t) ) break;
! 1160: }
! 1161: if ( j != clen ) break;
! 1162: }
! 1163: if ( i != row ) return 0;
! 1164: else return 1;
! 1165: }
! 1166:
! 1167: void isqrtgz(GZ a,GZ *r)
! 1168: {
! 1169: int k;
! 1170: GZ x,t,x2,xh,quo,rem;
! 1171: Q two;
! 1172:
! 1173: if ( !a ) *r = 0;
! 1174: else if ( UNIGZ(a) ) *r = ONEGZ;
! 1175: else {
! 1176: k = gz_bits((GQ)a); /* a <= 2^k-1 */
! 1177: bshiftgz(ONEGZ,-((k>>1)+(k&1)),&x); /* a <= x^2 */
! 1178: STOQ(2,two);
! 1179: while ( 1 ) {
! 1180: pwrgz(x,two,&t);
! 1181: if ( cmpgz(t,a) <= 0 ) {
! 1182: *r = x; return;
! 1183: } else {
! 1184: if ( mpz_tstbit(BDY(x),0) ) addgz(x,a,&t);
! 1185: else t = a;
! 1186: bshiftgz(x,-1,&x2); divqrgz(t,x2,&quo,&rem);
! 1187: bshiftgz(x,1,&xh); addgz(quo,xh,&x);
! 1188: }
! 1189: }
! 1190: }
! 1191: }
! 1192:
! 1193: void bshiftgz(GZ a,int n,GZ *r)
! 1194: {
! 1195: mpz_t t;
! 1196:
! 1197: if ( !a ) *r = 0;
! 1198: else if ( n == 0 ) *r = a;
! 1199: else if ( n < 0 ) {
! 1200: mpz_init(t); mpz_mul_2exp(t,BDY(a),-n); MPZTOGZ(t,*r);
! 1201: } else {
! 1202: mpz_init(t); mpz_fdiv_q_2exp(t,BDY(a),n);
! 1203: if ( !mpz_sgn(t) ) *r = 0;
! 1204: else MPZTOGZ(t,*r);
! 1205: }
! 1206: }
! 1207: