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