Annotation of OpenXM_contrib2/asir2018/engine/Q.c, Revision 1.1
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: Z ONE;
! 8: int lf_lazy;
! 9: Z current_mod_lf;
! 10: int current_mod_lf_size;
! 11: gmp_randstate_t GMP_RAND;
! 12:
! 13: void isqrtz(Z a,Z *r);
! 14: void bshiftz(Z a,int n,Z *r);
! 15:
! 16: void *gc_realloc(void *p,size_t osize,size_t nsize)
! 17: {
! 18: return (void *)Risa_GC_realloc(p,nsize);
! 19: }
! 20:
! 21: void gc_free(void *p,size_t size)
! 22: {
! 23: Risa_GC_free(p);
! 24: }
! 25:
! 26: void init_gmpq()
! 27: {
! 28: mp_set_memory_functions(Risa_GC_malloc_atomic,gc_realloc,gc_free);
! 29:
! 30: mpz_init(ONEMPZ); mpz_set_ui(ONEMPZ,1); MPZTOZ(ONEMPZ,ONE);
! 31: gmp_randinit_default(GMP_RAND);
! 32: }
! 33:
! 34: Z utoz(unsigned int u)
! 35: {
! 36: mpz_t z;
! 37: Z r;
! 38:
! 39: if ( !u ) return 0;
! 40: mpz_init(z); mpz_set_ui(z,u); MPZTOZ(z,r); return r;
! 41: }
! 42:
! 43: Z stoz(int s)
! 44: {
! 45: mpz_t z;
! 46: Z r;
! 47:
! 48: if ( !s ) return 0;
! 49: mpz_init(z); mpz_set_si(z,s); MPZTOZ(z,r); return r;
! 50: }
! 51:
! 52: int sgnz(Z z)
! 53: {
! 54: if ( !z ) return 0;
! 55: else return mpz_sgn(BDY(z));
! 56: }
! 57:
! 58: void nmq(Q q,Z *r)
! 59: {
! 60: if ( !q ) *r = 0;
! 61: else if ( INT(q) ) *r = (Z)q;
! 62: else {
! 63: MPZTOZ(mpq_numref(BDY(q)),*r);
! 64: }
! 65: }
! 66:
! 67: void dnq(Q q,Z *r)
! 68: {
! 69: if ( !q ) *r = 0;
! 70: else if ( INT(q) ) *r = ONE;
! 71: else {
! 72: MPZTOZ(mpq_denref(BDY(q)),*r);
! 73: }
! 74: }
! 75:
! 76: int sgnq(Q q)
! 77: {
! 78: if ( !q ) return 0;
! 79: else if ( q->z ) return mpz_sgn(BDY((Z)q));
! 80: else return mpz_sgn(mpq_numref(BDY(q)));
! 81: }
! 82:
! 83: Q mpqtozq(mpq_t a)
! 84: {
! 85: Z z;
! 86: Q q;
! 87:
! 88: if ( INTMPQ(a) ) {
! 89: MPZTOZ(mpq_numref(a),z); return (Q)z;
! 90: } else {
! 91: MPQTOQ(a,q); return q;
! 92: }
! 93: }
! 94:
! 95: void dupz(Z a,Z *b)
! 96: {
! 97: mpz_t t;
! 98:
! 99: if ( !a ) *b = a;
! 100: else {
! 101: mpz_init(t); mpz_set(t,BDY(a)); MPZTOZ(t,*b);
! 102: }
! 103: }
! 104:
! 105: int n_bits_z(Z a)
! 106: {
! 107: return a ? mpz_sizeinbase(BDY(a),2) : 0;
! 108: }
! 109:
! 110: void addz(Z n1,Z n2,Z *nr)
! 111: {
! 112: mpz_t t;
! 113: int s1,s2;
! 114:
! 115: if ( !n1 ) *nr = n2;
! 116: else if ( !n2 ) *nr = n1;
! 117: else if ( !n1->z || !n2->z )
! 118: error("addz : invalid argument");
! 119: else {
! 120: mpz_init(t); mpz_add(t,BDY(n1),BDY(n2)); MPZTOZ(t,*nr);
! 121: }
! 122: }
! 123:
! 124: void subz(Z n1,Z n2,Z *nr)
! 125: {
! 126: mpz_t t;
! 127:
! 128: if ( !n1 ) {
! 129: if ( !n2 )
! 130: *nr = 0;
! 131: else
! 132: chsgnz(n2,nr);
! 133: } else if ( !n2 )
! 134: *nr = n1;
! 135: else if ( n1 == n2 )
! 136: *nr = 0;
! 137: else if ( !n1->z || !n2->z )
! 138: error("subz : invalid argument");
! 139: else {
! 140: mpz_init(t); mpz_sub(t,BDY(n1),BDY(n2)); MPZTOZ(t,*nr);
! 141: }
! 142: }
! 143:
! 144: void mulz(Z n1,Z n2,Z *nr)
! 145: {
! 146: mpz_t t;
! 147:
! 148: if ( !n1 || !n2 ) *nr = 0;
! 149: else if ( !n1->z || !n2->z )
! 150: error("mulz : invalid argument");
! 151: else if ( UNIQ(n1) ) *nr = n2;
! 152: else if ( UNIQ(n2) ) *nr = n1;
! 153: else if ( MUNIQ(n1) ) chsgnz(n2,nr);
! 154: else if ( MUNIQ(n2) ) chsgnz(n1,nr);
! 155: else {
! 156: mpz_init(t); mpz_mul(t,BDY(n1),BDY(n2)); MPZTOZ(t,*nr);
! 157: }
! 158: }
! 159:
! 160: /* nr += n1*n2 */
! 161:
! 162: void muladdtoz(Z n1,Z n2,Z *nr)
! 163: {
! 164: Z t;
! 165:
! 166: if ( n1 && n2 ) {
! 167: if ( !(*nr) ) {
! 168: NEWZ(t); mpz_init(BDY(t)); *nr = t;
! 169: }
! 170: mpz_addmul(BDY(*nr),BDY(n1),BDY(n2));
! 171: }
! 172: }
! 173:
! 174: /* nr += n1*u */
! 175:
! 176: void mul1addtoz(Z n1,long u,Z *nr)
! 177: {
! 178: Z t;
! 179:
! 180: if ( n1 && u ) {
! 181: if ( !(*nr) ) {
! 182: NEWZ(t); mpz_init(BDY(t)); *nr = t;
! 183: }
! 184: if ( u >= 0 )
! 185: mpz_addmul_ui(BDY(*nr),BDY(n1),(unsigned long)u);
! 186: else
! 187: mpz_submul_ui(BDY(*nr),BDY(n1),(unsigned long)(-u));
! 188: }
! 189: }
! 190:
! 191: void mul1z(Z n1,long n2,Z *nr)
! 192: {
! 193: mpz_t t;
! 194:
! 195: if ( !n1 || !n2 ) *nr = 0;
! 196: else {
! 197: mpz_init(t); mpz_mul_si(t,BDY(n1),n2); MPZTOZ(t,*nr);
! 198: }
! 199: }
! 200:
! 201: void divz(Z n1,Z n2,Z *nq)
! 202: {
! 203: mpz_t t;
! 204: mpq_t a, b, q;
! 205:
! 206: if ( !n2 ) {
! 207: error("division by 0");
! 208: *nq = 0;
! 209: } else if ( !n1 )
! 210: *nq = 0;
! 211: else if ( n1 == n2 ) {
! 212: mpz_init(t); mpz_set_ui(t,1); MPZTOZ(t,*nq);
! 213: } else {
! 214: MPZTOMPQ(BDY(n1),a); MPZTOMPQ(BDY(n2),b);
! 215: mpq_init(q); mpq_div(q,a,b); *nq = (Z)mpqtozq(q);
! 216: }
! 217: }
! 218:
! 219: void remz(Z n1,Z n2,Z *nr)
! 220: {
! 221: mpz_t r;
! 222:
! 223: if ( !n2 ) {
! 224: error("division by 0");
! 225: *nr = 0;
! 226: } else if ( !n1 || n1 == n2 )
! 227: *nr = 0;
! 228: else if ( !n1->z || !n2->z )
! 229: error("remz : invalid argument");
! 230: else {
! 231: mpz_init(r);
! 232: mpz_mod(r,BDY(n1),BDY(n2));
! 233: if ( !mpz_sgn(r) ) *nr = 0;
! 234: else MPZTOZ(r,*nr);
! 235: }
! 236: }
! 237:
! 238: void divqrz(Z n1,Z n2,Z *nq,Z *nr)
! 239: {
! 240: mpz_t t, a, b, q, r;
! 241:
! 242: if ( !n2 ) {
! 243: error("division by 0");
! 244: *nq = 0; *nr = 0;
! 245: } else if ( !n1 ) {
! 246: *nq = 0; *nr = 0;
! 247: } else if ( !n1->z || !n2->z )
! 248: error("divqrz : invalid argument");
! 249: else if ( n1 == n2 ) {
! 250: mpz_init(t); mpz_set_ui(t,1); MPZTOZ(t,*nq); *nr = 0;
! 251: } else {
! 252: mpz_init(q); mpz_init(r);
! 253: mpz_fdiv_qr(q,r,BDY(n1),BDY(n2));
! 254: if ( !mpz_sgn(q) ) *nq = 0;
! 255: else MPZTOZ(q,*nq);
! 256: if ( !mpz_sgn(r) ) *nr = 0;
! 257: else MPZTOZ(r,*nr);
! 258: }
! 259: }
! 260:
! 261: void divsz(Z n1,Z n2,Z *nq)
! 262: {
! 263: mpz_t t;
! 264: mpq_t a, b, q;
! 265:
! 266: if ( !n2 ) {
! 267: error("division by 0");
! 268: *nq = 0;
! 269: } else if ( !n1 )
! 270: *nq = 0;
! 271: else if ( !n1->z || !n2->z )
! 272: error("divsz : invalid argument");
! 273: else if ( n1 == n2 ) {
! 274: mpz_init(t); mpz_set_ui(t,1); MPZTOZ(t,*nq);
! 275: } else {
! 276: mpz_init(t); mpz_divexact(t,BDY(n1),BDY(n2)); MPZTOZ(t,*nq);
! 277: }
! 278: }
! 279:
! 280: void chsgnz(Z n,Z *nr)
! 281: {
! 282: mpz_t t;
! 283:
! 284: if ( !n )
! 285: *nr = 0;
! 286: else if ( !n->z )
! 287: error("chsgnz : invalid argument");
! 288: else {
! 289: t[0] = BDY(n)[0]; mpz_neg(t,t); MPZTOZ(t,*nr);
! 290: }
! 291: }
! 292:
! 293: void absz(Z n,Z *nr)
! 294: {
! 295: if ( !n ) *nr = 0;
! 296: else if ( !n->z )
! 297: error("absz : invalid argument");
! 298: else if ( sgnz(n) < 0 ) chsgnz(n,nr);
! 299: else *nr = n;
! 300: }
! 301:
! 302: int evenz(Z n)
! 303: {
! 304: return !n ? 1 : mpz_even_p(BDY(n));
! 305: }
! 306:
! 307: int smallz(Z n)
! 308: {
! 309: if ( !n ) return 1;
! 310: else if ( INT(n) && mpz_fits_sint_p(BDY(n)) ) return 1;
! 311: else return 0;
! 312: }
! 313:
! 314: void pwrz(Z n1,Z n,Z *nr)
! 315: {
! 316: mpq_t t,q;
! 317: mpz_t z;
! 318: Q p,r;
! 319:
! 320: if ( !n || UNIQ(n1) ) *nr = ONE;
! 321: else if ( !n1 ) *nr = 0;
! 322: else if ( !n->z || !n1->z )
! 323: error("pwrz : invalid argument");
! 324: else if ( MUNIQ(n1) ) {
! 325: if ( mpz_even_p(BDY((Z)n)) ) *nr = ONE;
! 326: else *nr = n1;
! 327: } else if ( !smallz(n) ) {
! 328: error("exponent too big."); *nr = 0;
! 329: } else if ( n1->z && mpz_sgn(BDY((Z)n))>0 ) {
! 330: mpz_init(z); mpz_pow_ui(z,BDY(n1),QTOS(n)); MPZTOZ(z,*nr);
! 331: } else {
! 332: MPZTOMPQ(BDY(n1),q); MPQTOQ(q,r);
! 333: pwrq(r,(Q)n,&p); *nr = (Z)p;
! 334: }
! 335: }
! 336:
! 337: int cmpz(Z q1,Z q2)
! 338: {
! 339: int sgn;
! 340:
! 341: if ( !q1 ) {
! 342: if ( !q2 )
! 343: return 0;
! 344: else
! 345: return -mpz_sgn(BDY(q2));
! 346: } else if ( !q2 )
! 347: return mpz_sgn(BDY(q1));
! 348: else if ( !q1->z || !q2->z )
! 349: error("mpqz : invalid argument");
! 350: else if ( (sgn = mpz_sgn(BDY(q1))) != mpz_sgn(BDY(q2)) )
! 351: return sgn;
! 352: else {
! 353: sgn = mpz_cmp(BDY(q1),BDY(q2));
! 354: if ( sgn > 0 ) return 1;
! 355: else if ( sgn < 0 ) return -1;
! 356: else return 0;
! 357: }
! 358: }
! 359:
! 360: void gcdz(Z n1,Z n2,Z *nq)
! 361: {
! 362: mpz_t t;
! 363:
! 364: if ( !n1 ) *nq = n2;
! 365: else if ( !n2 ) *nq = n1;
! 366: else if ( !n1->z || !n2->z )
! 367: error("gcdz : invalid argument");
! 368: else {
! 369: mpz_init(t); mpz_gcd(t,BDY(n1),BDY(n2));
! 370: MPZTOZ(t,*nq);
! 371: }
! 372: }
! 373:
! 374: void invz(Z n1,Z n2,Z *nq)
! 375: {
! 376: mpz_t t;
! 377:
! 378: if ( !n1 || !n2 || !n1->z || !n2->z )
! 379: error("invz : invalid argument");
! 380: mpz_init(t); mpz_invert(t,BDY(n1),BDY(n2));
! 381: MPZTOZ(t,*nq);
! 382: }
! 383:
! 384: void lcmz(Z n1,Z n2,Z *nq)
! 385: {
! 386: Z g,t;
! 387:
! 388: if ( !n1 || !n2 ) *nq = 0;
! 389: else if ( !n1->z || !n2->z )
! 390: error("lcmz : invalid argument");
! 391: else {
! 392: gcdz(n1,n2,&g); divsz(n1,g,&t);
! 393: mulz(n2,t,nq);
! 394: }
! 395: }
! 396:
! 397: void gcdvz(VECT v,Z *q)
! 398: {
! 399: int n,i;
! 400: Z *b;
! 401: Z g,g1;
! 402:
! 403: n = v->len;
! 404: b = (Z *)v->body;
! 405: g = b[0];
! 406: for ( i = 1; i < n; i++ ) {
! 407: gcdz(g,b[i],&g1); g = g1;
! 408: }
! 409: *q = g;
! 410: }
! 411:
! 412: void gcdvz_estimate(VECT v,Z *q)
! 413: {
! 414: int n,m,i;
! 415: Z s,t,u;
! 416: Z *b;
! 417:
! 418: n = v->len;
! 419: b = (Z *)v->body;
! 420: if ( n == 1 ) {
! 421: if ( mpz_sgn(BDY(b[0]))<0 ) chsgnz(b[0],q);
! 422: else *q = b[0];
! 423: }
! 424: m = n/2;
! 425: for ( i = 0, s = 0; i < m; i++ ) {
! 426: if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subz(s,b[i],&u);
! 427: else addz(s,b[i],&u);
! 428: s = u;
! 429: }
! 430: for ( i = 0, t = 0; i < n; i++ ) {
! 431: if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subz(t,b[i],&u);
! 432: else addz(t,b[i],&u);
! 433: t = u;
! 434: }
! 435: gcdz(s,t,q);
! 436: }
! 437:
! 438: void factorialz(unsigned int n,Z *nr)
! 439: {
! 440: mpz_t a;
! 441: mpz_init(a);
! 442: mpz_fac_ui(a,n);
! 443: MPZTOZ(a,*nr);
! 444: }
! 445:
! 446: void randomz(int blen,Z *nr)
! 447: {
! 448: mpz_t z;
! 449:
! 450: mpz_init(z);
! 451: mpz_urandomb(z,GMP_RAND,blen);
! 452: MPZTOZ(z,*nr);
! 453: }
! 454:
! 455: int tstbitz(Z n,int k)
! 456: {
! 457: if ( !n || !n->z )
! 458: error("tstbitz : invalid argument");
! 459: return !n ? 0 : mpz_tstbit(BDY(n),k);
! 460: }
! 461:
! 462: void addq(Q n1,Q n2,Q *nr)
! 463: {
! 464: mpq_t q1,q2,t;
! 465:
! 466: if ( !n1 ) *nr = n2;
! 467: else if ( !n2 ) *nr = n1;
! 468: else if ( n1->z && n2->z )
! 469: addz((Z)n1,(Z)n2,(Z *)nr);
! 470: else {
! 471: if ( n1->z ) MPZTOMPQ(BDY((Z)n1),q1);
! 472: else q1[0] = BDY(n1)[0];
! 473: if ( n2->z ) MPZTOMPQ(BDY((Z)n2),q2);
! 474: else q2[0] = BDY(n2)[0];
! 475: mpq_init(t); mpq_add(t,q1,q2); *nr = mpqtozq(t);
! 476: }
! 477: }
! 478:
! 479: void subq(Q n1,Q n2,Q *nr)
! 480: {
! 481: mpq_t q1,q2,t;
! 482:
! 483: if ( !n1 ) {
! 484: if ( !n2 ) *nr = 0;
! 485: else if ( n1->z ) chsgnz((Z)n1,(Z *)nr);
! 486: else {
! 487: mpq_init(t); mpq_neg(t,BDY(n2)); MPQTOQ(t,*nr);
! 488: }
! 489: } else if ( !n2 ) *nr = n1;
! 490: else if ( n1 == n2 ) *nr = 0;
! 491: else if ( n1->z && n2->z )
! 492: subz((Z)n1,(Z)n2,(Z *)nr);
! 493: else {
! 494: if ( n1->z ) MPZTOMPQ(BDY((Z)n1),q1);
! 495: else q1[0] = BDY(n1)[0];
! 496: if ( n2->z ) MPZTOMPQ(BDY((Z)n2),q2);
! 497: else q2[0] = BDY(n2)[0];
! 498: mpq_init(t); mpq_sub(t,q1,q2); *nr = mpqtozq(t);
! 499: }
! 500: }
! 501:
! 502: void mulq(Q n1,Q n2,Q *nr)
! 503: {
! 504: mpq_t t,q1,q2;
! 505:
! 506: if ( !n1 || !n2 ) *nr = 0;
! 507: else if ( n1->z && n2->z )
! 508: mulz((Z)n1,(Z)n2,(Z *)nr);
! 509: else {
! 510: if ( n1->z ) MPZTOMPQ(BDY((Z)n1),q1);
! 511: else q1[0] = BDY(n1)[0];
! 512: if ( n2->z ) MPZTOMPQ(BDY((Z)n2),q2);
! 513: else q2[0] = BDY(n2)[0];
! 514: mpq_init(t); mpq_mul(t,q1,q2); *nr = mpqtozq(t);
! 515: }
! 516: }
! 517:
! 518: void divq(Q n1,Q n2,Q *nq)
! 519: {
! 520: mpq_t t,q1,q2;
! 521:
! 522: if ( !n2 ) {
! 523: error("division by 0");
! 524: *nq = 0;
! 525: return;
! 526: } else if ( !n1 ) *nq = 0;
! 527: else if ( n1 == n2 ) *nq = (Q)ONE;
! 528: else {
! 529: if ( n1->z ) MPZTOMPQ(BDY((Z)n1),q1);
! 530: else q1[0] = BDY(n1)[0];
! 531: if ( n2->z ) MPZTOMPQ(BDY((Z)n2),q2);
! 532: else q2[0] = BDY(n2)[0];
! 533: mpq_init(t); mpq_div(t,q1,q2); *nq = mpqtozq(t);
! 534: }
! 535: }
! 536:
! 537: void invq(Q n,Q *nr)
! 538: {
! 539: Z nm,dn;
! 540:
! 541: if ( INT(n) )
! 542: divq((Q)ONE,n,nr);
! 543: else {
! 544: nmq(n,&nm);
! 545: dnq(n,&dn);
! 546: divq((Q)dn,(Q)nm,nr);
! 547: }
! 548: }
! 549:
! 550: void chsgnq(Q n,Q *nr)
! 551: {
! 552: mpq_t t;
! 553:
! 554: if ( !n ) *nr = 0;
! 555: else if (n->z ) chsgnz((Z)n,(Z *)nr);
! 556: else {
! 557: mpq_init(t); mpq_neg(t,BDY(n)); MPQTOQ(t,*nr);
! 558: }
! 559: }
! 560:
! 561: void absq(Q n,Q *nr)
! 562: {
! 563: if ( !n ) *nr = 0;
! 564: else if ( n->z ) absz((Z)n,(Z *)nr);
! 565: else if ( sgnq(n) < 0 ) chsgnq(n,nr);
! 566: else *nr = n;
! 567: }
! 568:
! 569: void pwrq(Q n1,Q n,Q *nr)
! 570: {
! 571: int e;
! 572: mpz_t nm,dn;
! 573: mpq_t t;
! 574:
! 575: if ( !n || UNIQ((Z)n1) || UNIQ(n1) ) *nr = (Q)ONE;
! 576: else if ( !n1 ) *nr = 0;
! 577: else if ( !INT(n) ) {
! 578: error("can't calculate fractional power."); *nr = 0;
! 579: } else if ( !smallz((Z)n) ) {
! 580: error("exponent too big."); *nr = 0;
! 581: } else {
! 582: e = QTOS(n);
! 583: if ( e < 0 ) {
! 584: e = -e;
! 585: if ( n1->z ) {
! 586: nm[0] = ONEMPZ[0];
! 587: dn[0] = BDY((Z)n1)[0];
! 588: } else {
! 589: nm[0] = mpq_denref(BDY(n1))[0];
! 590: dn[0] = mpq_numref(BDY(n1))[0];
! 591: }
! 592: } else {
! 593: if ( n1->z ) {
! 594: nm[0] = BDY((Z)n1)[0];
! 595: dn[0] = ONEMPZ[0];
! 596: } else {
! 597: nm[0] = mpq_numref(BDY(n1))[0];
! 598: dn[0] = mpq_denref(BDY(n1))[0];
! 599: }
! 600: }
! 601: mpq_init(t);
! 602: mpz_pow_ui(mpq_numref(t),nm,e); mpz_pow_ui(mpq_denref(t),dn,e);
! 603: *nr = mpqtozq(t);
! 604: }
! 605: }
! 606:
! 607: int cmpq(Q n1,Q n2)
! 608: {
! 609: mpq_t q1,q2;
! 610: int sgn;
! 611:
! 612: if ( !n1 ) {
! 613: if ( !n2 ) return 0;
! 614: else return (n2->z) ? -mpz_sgn(BDY((Z)n2)) : -mpq_sgn(BDY(n2));
! 615: } if ( !n2 ) return (n1->z) ? mpz_sgn(BDY((Z)n1)) : mpq_sgn(BDY(n1));
! 616: else if ( n1->z && n2->z )
! 617: return cmpz((Z)n1,(Z)n2);
! 618: else if ( (sgn = mpq_sgn(BDY(n1))) != mpq_sgn(BDY(n2)) ) return sgn;
! 619: else {
! 620: if ( n1->z ) MPZTOMPQ(BDY((Z)n1),q1);
! 621: else q1[0] = BDY(n1)[0];
! 622: if ( n2->z ) MPZTOMPQ(BDY((Z)n2),q2);
! 623: else q2[0] = BDY(n2)[0];
! 624: sgn = mpq_cmp(q1,q2);
! 625: if ( sgn > 0 ) return 1;
! 626: else if ( sgn < 0 ) return -1;
! 627: else return 0;
! 628: }
! 629: }
! 630:
! 631: /* t = [nC0 nC1 ... nCn] */
! 632:
! 633: void mkbc(int n,Z *t)
! 634: {
! 635: int i;
! 636: Z c,d,iq;
! 637:
! 638: for ( t[0] = ONE, i = 1; i <= n/2; i++ ) {
! 639: STOQ(n-i+1,c); mulz(t[i-1],c,&d);
! 640: STOQ(i,iq); divsz(d,iq,&t[i]);
! 641: }
! 642: for ( ; i <= n; i++ )
! 643: t[i] = t[n-i];
! 644: }
! 645:
! 646: /*
! 647: * Dx^k*x^l = W(k,l,0)*x^l*Dx^k+W(k,l,1)*x^(l-1)*x^(k-1)*+...
! 648: *
! 649: * t = [W(k,l,0) W(k,l,1) ... W(k,l,min(k,l)]
! 650: * where W(k,l,i) = i! * kCi * lCi
! 651: */
! 652:
! 653: /* mod m table */
! 654: /* XXX : should be optimized */
! 655:
! 656: void mkwcm(int k,int l,int m,int *t)
! 657: {
! 658: int i,n;
! 659: Z *s;
! 660:
! 661: n = MIN(k,l);
! 662: s = (Z *)ALLOCA((n+1)*sizeof(Q));
! 663: mkwc(k,l,s);
! 664: for ( i = 0; i <= n; i++ ) {
! 665: t[i] = remqi((Q)s[i],m);
! 666: }
! 667: }
! 668:
! 669: void mkwc(int k,int l,Z *t)
! 670: {
! 671: mpz_t a,b,q,nm,z,u;
! 672: int i,n;
! 673:
! 674: n = MIN(k,l);
! 675: mpz_init_set_ui(z,1);
! 676: mpz_init(u); mpz_set(u,z); MPZTOZ(u,t[0]);
! 677: mpz_init(a); mpz_init(b); mpz_init(nm);
! 678: for ( i = 1; i <= n; i++ ) {
! 679: mpz_set_ui(a,k-i+1); mpz_set_ui(b,l-i+1); mpz_mul(nm,a,b);
! 680: mpz_mul(z,BDY(t[i-1]),nm); mpz_fdiv_q_ui(z,z,i);
! 681: mpz_init(u); mpz_set(u,z); MPZTOZ(u,t[i]);
! 682: }
! 683: }
! 684:
! 685: void lgp(P p,Z *g,Z *l);
! 686:
! 687: void ptozp(P p,int sgn,Q *c,P *pr)
! 688: {
! 689: Z nm,dn;
! 690:
! 691: if ( !p ) {
! 692: *c = 0; *pr = 0;
! 693: } else {
! 694: lgp(p,&nm,&dn);
! 695: divz(nm,dn,(Z *)c);
! 696: divsp(CO,p,(P)*c,pr);
! 697: }
! 698: }
! 699:
! 700: void lgp(P p,Z *g,Z *l)
! 701: {
! 702: DCP dc;
! 703: Z g1,g2,l1,l2,l3,l4;
! 704:
! 705: if ( NUM(p) ) {
! 706: if ( ((Q)p)->z ) {
! 707: MPZTOZ(BDY((Z)p),*g);
! 708: *l = ONE;
! 709: } else {
! 710: MPZTOZ(mpq_numref(BDY((Q)p)),*g);
! 711: MPZTOZ(mpq_denref(BDY((Q)p)),*l);
! 712: }
! 713: } else {
! 714: dc = DC(p); lgp(COEF(dc),g,l);
! 715: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
! 716: lgp(COEF(dc),&g1,&l1); gcdz(*g,g1,&g2); *g = g2;
! 717: gcdz(*l,l1,&l2); mulz(*l,l1,&l3); divz(l3,l2,l);
! 718: }
! 719: }
! 720: }
! 721:
! 722: void qltozl(Q *w,int n,Z *dvr)
! 723: {
! 724: Z nm,dn;
! 725: Z g,g1,l1,l2,l3;
! 726: Q c;
! 727: int i;
! 728: struct oVECT v;
! 729:
! 730: for ( i = 0; i < n; i++ )
! 731: if ( w[i] && !w[i]->z )
! 732: break;
! 733: if ( i == n ) {
! 734: v.id = O_VECT; v.len = n; v.body = (pointer *)w;
! 735: gcdvz(&v,dvr); return;
! 736: }
! 737: for ( i = 0; !w[i]; i++ );
! 738: c = w[i];
! 739: if ( !c->z ) {
! 740: MPZTOZ(mpq_numref(BDY(c)),nm); MPZTOZ(mpq_denref(BDY(c)),dn);
! 741: } else {
! 742: MPZTOZ(BDY((Z)c),nm); dn = ONE;
! 743: }
! 744: for ( i++; i < n; i++ ) {
! 745: c = w[i];
! 746: if ( !c ) continue;
! 747: if ( !c->z ) {
! 748: MPZTOZ(mpq_numref(BDY(c)),g1); MPZTOZ(mpq_denref(BDY(c)),l1);
! 749: } else {
! 750: MPZTOZ(BDY((Z)c),g1); l1 = ONE;
! 751: }
! 752: gcdz(nm,g1,&g); nm = g;
! 753: gcdz(dn,l1,&l2); mulz(dn,l1,&l3); divz(l3,l2,&dn);
! 754: }
! 755: divz(nm,dn,dvr);
! 756: }
! 757:
! 758: int z_bits(Q q)
! 759: {
! 760: if ( !q ) return 0;
! 761: else if ( q->z ) return mpz_sizeinbase(BDY((Z)q),2);
! 762: else
! 763: return mpz_sizeinbase(mpq_numref(BDY(q)),2)
! 764: + mpz_sizeinbase(mpq_denref(BDY(q)),2);
! 765: }
! 766:
! 767: int zp_mag(P p)
! 768: {
! 769: int s;
! 770: DCP dc;
! 771:
! 772: if ( !p ) return 0;
! 773: else if ( OID(p) == O_N ) return z_bits((Q)p);
! 774: else {
! 775: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) s += zp_mag(COEF(dc));
! 776: return s;
! 777: }
! 778: }
! 779:
! 780: void makesubstz(VL v,NODE *s)
! 781: {
! 782: NODE r,r0;
! 783: Z q;
! 784: unsigned int n;
! 785:
! 786: for ( r0 = 0; v; v = NEXT(v) ) {
! 787: NEXTNODE(r0,r); BDY(r) = (pointer)v->v;
! 788: #if defined(_PA_RISC1_1)
! 789: n = mrand48()&BMASK; q = utoz(n);
! 790: #else
! 791: n = random(); q = utoz(n);
! 792: #endif
! 793: NEXTNODE(r0,r); BDY(r) = (pointer)q;
! 794: }
! 795: if ( r0 ) NEXT(r) = 0;
! 796: *s = r0;
! 797: }
! 798:
! 799: unsigned int remqi(Q a,unsigned int mod)
! 800: {
! 801: unsigned int c,nm,dn;
! 802: mpz_t r;
! 803:
! 804: if ( !a ) return 0;
! 805: else if ( a->z ) {
! 806: mpz_init(r);
! 807: c = mpz_fdiv_r_ui(r,BDY((Z)a),mod);
! 808: } else {
! 809: mpz_init(r);
! 810: nm = mpz_fdiv_r_ui(r,mpq_numref(BDY(a)),mod);
! 811: dn = mpz_fdiv_r_ui(r,mpq_denref(BDY(a)),mod);
! 812: dn = invm(dn,mod);
! 813: DMAR(nm,dn,0,mod,c);
! 814: }
! 815: return c;
! 816: }
! 817:
! 818: extern int DP_Print;
! 819:
! 820: #define F4_INTRAT_PERIOD 8
! 821:
! 822: int generic_gauss_elim(MAT mat,MAT *nm,Z *dn,int **rindp,int **cindp)
! 823: {
! 824: int **wmat;
! 825: Z **bmat,**tmat,*bmi,*tmi;
! 826: Z 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;
! 829: MAT r,crmat;
! 830: int ret;
! 831:
! 832: bmat = (Z **)mat->body;
! 833: row = mat->row; col = mat->col;
! 834: wmat = (int **)almat(row,col);
! 835: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
! 836: wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
! 837: for ( ind = 0; ; ind++ ) {
! 838: if ( DP_Print ) {
! 839: fprintf(asir_out,"."); fflush(asir_out);
! 840: }
! 841: md = get_lprime(ind);
! 842: for ( i = 0; i < row; i++ )
! 843: for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
! 844: wmi[j] = remqi((Q)bmi[j],md);
! 845: rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);
! 846: if ( !ind ) {
! 847: RESET:
! 848: m1 = utoz(md);
! 849: rank0 = rank;
! 850: bcopy(wcolstat,colstat,col*sizeof(int));
! 851: MKMAT(crmat,rank,col-rank);
! 852: MKMAT(r,rank,col-rank); *nm = r;
! 853: tmat = (Z **)crmat->body;
! 854: for ( i = 0; i < rank; i++ )
! 855: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
! 856: if ( !colstat[j] ) tmi[k++] = utoz(wmi[j]);
! 857: } else {
! 858: if ( rank < rank0 ) {
! 859: if ( DP_Print ) {
! 860: fprintf(asir_out,"lower rank matrix; continuing...\n");
! 861: fflush(asir_out);
! 862: }
! 863: continue;
! 864: } else if ( rank > rank0 ) {
! 865: if ( DP_Print ) {
! 866: fprintf(asir_out,"higher rank matrix; resetting...\n");
! 867: fflush(asir_out);
! 868: }
! 869: goto RESET;
! 870: } else {
! 871: for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
! 872: if ( j < col ) {
! 873: if ( DP_Print ) {
! 874: fprintf(asir_out,"inconsitent colstat; resetting...\n");
! 875: fflush(asir_out);
! 876: }
! 877: goto RESET;
! 878: }
! 879: }
! 880:
! 881: inv = invm(remqi((Q)m1,md),md);
! 882: m2 = utoz(md); mulz(m1,m2,&m3);
! 883: for ( i = 0; i < rank; i++ )
! 884: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
! 885: if ( !colstat[j] ) {
! 886: if ( tmi[k] ) {
! 887: /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
! 888: t = remqi((Q)tmi[k],md);
! 889: if ( wmi[j] >= t ) t = wmi[j]-t;
! 890: else t = md-(t-wmi[j]);
! 891: DMAR(t,inv,0,md,t1)
! 892: u = utoz(t1); mulz(m1,u,&s);
! 893: addz(tmi[k],s,&u); tmi[k] = u;
! 894: } else if ( wmi[j] ) {
! 895: /* f3 = m1*(m1 mod m2)^(-1)*f2 */
! 896: DMAR(wmi[j],inv,0,md,t)
! 897: u = utoz(t); mulz(m1,u,&s); tmi[k] = s;
! 898: }
! 899: k++;
! 900: }
! 901: m1 = m3;
! 902: if ( ind % F4_INTRAT_PERIOD )
! 903: ret = 0;
! 904: else
! 905: ret = intmtoratm(crmat,m1,*nm,dn);
! 906: if ( ret ) {
! 907: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
! 908: *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
! 909: for ( j = k = l = 0; j < col; j++ )
! 910: if ( colstat[j] ) rind[k++] = j;
! 911: else cind[l++] = j;
! 912: if ( gensolve_check(mat,*nm,*dn,rind,cind) )
! 913: return rank;
! 914: }
! 915: }
! 916: }
! 917: }
! 918:
! 919: int generic_gauss_elim2(MAT mat,MAT *nm,Z *dn,int **rindp,int **cindp)
! 920: {
! 921:
! 922: MAT full;
! 923: Z **bmat,**b;
! 924: Z *bmi;
! 925: Z dn0;
! 926: int row,col,md,i,j,rank,ret;
! 927: int **wmat;
! 928: int *wmi;
! 929: int *colstat,*rowstat;
! 930:
! 931: bmat = (Z **)mat->body;
! 932: row = mat->row; col = mat->col;
! 933: wmat = (int **)almat(row,col);
! 934: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
! 935: rowstat = (int *)MALLOC_ATOMIC(row*sizeof(int));
! 936: /* XXX */
! 937: md = get_lprime(0);
! 938: for ( i = 0; i < row; i++ )
! 939: for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
! 940: wmi[j] = remqi((Q)bmi[j],md);
! 941: rank = generic_gauss_elim_mod2(wmat,row,col,md,colstat,rowstat);
! 942: b = (Z **)MALLOC(rank*sizeof(Z));
! 943: for ( i = 0; i < rank; i++ ) b[i] = bmat[rowstat[i]];
! 944: NEWMAT(full); full->row = rank; full->col = col; full->body = (pointer **)b;
! 945: ret = generic_gauss_elim_full(full,nm,dn,rindp,cindp);
! 946: if ( !ret ) {
! 947: rank = generic_gauss_elim(mat,nm,&dn0,rindp,cindp);
! 948: for ( i = 0; i < rank; i++ ) dn[i] = dn0;
! 949: }
! 950: return rank;
! 951: }
! 952:
! 953: int generic_gauss_elim_full(MAT mat,MAT *nm,Z *dn,int **rindp,int **cindp)
! 954: {
! 955: int **wmat;
! 956: int *stat;
! 957: Z **bmat,**tmat,*bmi,*tmi,*ri;
! 958: Z q,m1,m2,m3,s,u;
! 959: int *wmi,*colstat,*wcolstat,*rind,*cind;
! 960: int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv,h;
! 961: MAT r,crmat;
! 962: int ret,initialized,done;
! 963:
! 964: initialized = 0;
! 965: bmat = (Z **)mat->body;
! 966: row = mat->row; col = mat->col;
! 967: wmat = (int **)almat(row,col);
! 968: stat = (int *)MALLOC_ATOMIC(row*sizeof(int));
! 969: for ( i = 0; i < row; i++ ) stat[i] = 0;
! 970: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
! 971: wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
! 972: for ( ind = 0; ; ind++ ) {
! 973: if ( DP_Print ) {
! 974: fprintf(asir_out,"."); fflush(asir_out);
! 975: }
! 976: md = get_lprime(ind);
! 977: for ( i = 0; i < row; i++ )
! 978: for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
! 979: wmi[j] = remqi((Q)bmi[j],md);
! 980: rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);
! 981: if ( rank < row ) continue;
! 982: if ( !initialized ) {
! 983: m1 = utoz(md);
! 984: bcopy(wcolstat,colstat,col*sizeof(int));
! 985: MKMAT(crmat,row,col-row);
! 986: MKMAT(r,row,col-row); *nm = r;
! 987: tmat = (Z **)crmat->body;
! 988: for ( i = 0; i < row; i++ )
! 989: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
! 990: if ( !colstat[j] ) tmi[k++] = utoz(wmi[j]);
! 991: initialized = 1;
! 992: } else {
! 993: for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
! 994: if ( j < col ) continue;
! 995:
! 996: inv = invm(remqi((Q)m1,md),md);
! 997: m2 = utoz(md); mulz(m1,m2,&m3);
! 998: for ( i = 0; i < row; i++ )
! 999: switch ( stat[i] ) {
! 1000: case 1:
! 1001: /* consistency check */
! 1002: ri = (Z *)BDY(r)[i]; wmi = wmat[i];
! 1003: for ( j = 0; j < col; j++ ) if ( colstat[j] ) break;
! 1004: h = md-remqi((Q)dn[i],md);
! 1005: for ( j++, k = 0; j < col; j++ )
! 1006: if ( !colstat[j] ) {
! 1007: t = remqi((Q)ri[k],md);
! 1008: DMAR(wmi[i],h,t,md,t1);
! 1009: if ( t1 ) break;
! 1010: }
! 1011: if ( j == col ) { stat[i]++; break; }
! 1012: else {
! 1013: /* fall to the case 0 */
! 1014: stat[i] = 0;
! 1015: }
! 1016: case 0:
! 1017: tmi = tmat[i]; wmi = wmat[i];
! 1018: for ( j = k = 0; j < col; j++ )
! 1019: if ( !colstat[j] ) {
! 1020: if ( tmi[k] ) {
! 1021: /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
! 1022: t = remqi((Q)tmi[k],md);
! 1023: if ( wmi[j] >= t ) t = wmi[j]-t;
! 1024: else t = md-(t-wmi[j]);
! 1025: DMAR(t,inv,0,md,t1)
! 1026: u = utoz(t1); mulz(m1,u,&s);
! 1027: addz(tmi[k],s,&u); tmi[k] = u;
! 1028: } else if ( wmi[j] ) {
! 1029: /* f3 = m1*(m1 mod m2)^(-1)*f2 */
! 1030: DMAR(wmi[j],inv,0,md,t)
! 1031: u = utoz(t); mulz(m1,u,&s); tmi[k] = s;
! 1032: }
! 1033: k++;
! 1034: }
! 1035: break;
! 1036: case 2: default:
! 1037: break;
! 1038: }
! 1039: m1 = m3;
! 1040: if ( ind % 4 )
! 1041: ret = 0;
! 1042: else
! 1043: ret = intmtoratm2(crmat,m1,*nm,dn,stat);
! 1044: if ( ret ) {
! 1045: *rindp = rind = (int *)MALLOC_ATOMIC(row*sizeof(int));
! 1046: *cindp = cind = (int *)MALLOC_ATOMIC((col-row)*sizeof(int));
! 1047: for ( j = k = l = 0; j < col; j++ )
! 1048: if ( colstat[j] ) rind[k++] = j;
! 1049: else cind[l++] = j;
! 1050: return gensolve_check2(mat,*nm,dn,rind,cind);
! 1051: }
! 1052: }
! 1053: }
! 1054: }
! 1055:
! 1056: int generic_gauss_elim_direct(MAT mat,MAT *nm,Z *dn,int **rindp,int **cindp){
! 1057: Z **bmat,*s;
! 1058: Z u,v,w,x,d,t,y;
! 1059: int row,col,i,j,k,l,m,rank;
! 1060: int *colstat,*colpos,*cind;
! 1061: MAT r,in;
! 1062:
! 1063: row = mat->row; col = mat->col;
! 1064: MKMAT(in,row,col);
! 1065: for ( i = 0; i < row; i++ )
! 1066: for ( j = 0; j < col; j++ ) in->body[i][j] = mat->body[i][j];
! 1067: bmat = (Z **)in->body;
! 1068: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
! 1069: *rindp = colpos = (int *)MALLOC_ATOMIC(row*sizeof(int));
! 1070: for ( j = 0, rank = 0, d = ONE; j < col; j++ ) {
! 1071: for ( i = rank; i < row && !bmat[i][j]; i++ );
! 1072: if ( i == row ) { colstat[j] = 0; continue; }
! 1073: else { colstat[j] = 1; colpos[rank] = j; }
! 1074: if ( i != rank )
! 1075: for ( k = j; k < col; k++ ) {
! 1076: t = bmat[i][k]; bmat[i][k] = bmat[rank][k]; bmat[rank][k] = t;
! 1077: }
! 1078: for ( i = rank+1, v = bmat[rank][j]; i < row; i++ )
! 1079: for ( k = j, u = bmat[i][j]; k < col; k++ ) {
! 1080: mulz(bmat[i][k],v,&w); mulz(bmat[rank][k],u,&x);
! 1081: subz(w,x,&y); divsz(y,d,&bmat[i][k]);
! 1082: }
! 1083: d = v; rank++;
! 1084: }
! 1085: *dn = d;
! 1086: s = (Z *)MALLOC(col*sizeof(Z));
! 1087: for ( i = rank-1; i >= 0; i-- ) {
! 1088: for ( k = colpos[i]; k < col; k++ ) mulz(bmat[i][k],d,&s[k]);
! 1089: for ( m = rank-1; m > i; m-- ) {
! 1090: for ( k = colpos[m], u = bmat[i][k]; k < col; k++ ) {
! 1091: mulz(bmat[m][k],u,&w); subz(s[k],w,&x); s[k] = x;
! 1092: }
! 1093: }
! 1094: for ( k = colpos[i], u = bmat[i][k]; k < col; k++ )
! 1095: divz(s[k],u,&bmat[i][k]);
! 1096: }
! 1097: *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
! 1098: MKMAT(r,rank,col-rank); *nm = r;
! 1099: for ( j = 0, k = 0; j < col; j++ )
! 1100: if ( !colstat[j] ) {
! 1101: cind[k] = j;
! 1102: for ( i = 0; i < rank; i++ ) r->body[i][k] = bmat[i][j];
! 1103: k++;
! 1104: }
! 1105: return rank;
! 1106: }
! 1107:
! 1108: int intmtoratm(MAT mat,Z md,MAT nm,Z *dn)
! 1109: {
! 1110: Z t,s,b,dn0,dn1,nm1,q,u,unm,udn,dmy;
! 1111: int i,j,k,l,row,col,sgn,ret;
! 1112: Z **rmat,**tmat,*tmi,*nmk;
! 1113:
! 1114: if ( UNIQ(md) )
! 1115: return 0;
! 1116: row = mat->row; col = mat->col;
! 1117: bshiftz(md,1,&t);
! 1118: isqrtz(t,&s);
! 1119: bshiftz(s,64,&b);
! 1120: if ( !b ) b = ONE;
! 1121: dn0 = ONE;
! 1122: tmat = (Z **)mat->body;
! 1123: rmat = (Z **)nm->body;
! 1124: for ( i = 0; i < row; i++ )
! 1125: for ( j = 0, tmi = tmat[i]; j < col; j++ )
! 1126: if ( tmi[j] ) {
! 1127: mulz(tmi[j],dn0,&s);
! 1128: divqrz(s,md,&dmy,&u);
! 1129: ret = inttorat(u,md,b,&nm1,&dn1);
! 1130: if ( !ret ) return 0;
! 1131: else {
! 1132: if ( !UNIQ(dn1) ) {
! 1133: for ( k = 0; k < i; k++ )
! 1134: for ( l = 0, nmk = rmat[k]; l < col; l++ ) {
! 1135: mulz(nmk[l],dn1,&q); nmk[l] = q;
! 1136: }
! 1137: for ( l = 0, nmk = rmat[i]; l < j; l++ ) {
! 1138: mulz(nmk[l],dn1,&q); nmk[l] = q;
! 1139: }
! 1140: }
! 1141: rmat[i][j] = nm1;
! 1142: mulz(dn0,dn1,&q); dn0 = q;
! 1143: }
! 1144: }
! 1145: *dn = dn0;
! 1146: return 1;
! 1147: }
! 1148:
! 1149: int intmtoratm2(MAT mat,Z md,MAT nm,Z *dn,int *stat)
! 1150: {
! 1151: int row,col,i,j,ret;
! 1152: Z dn0,dn1,t,s,b;
! 1153: Z *w,*tmi;
! 1154: Z **tmat;
! 1155:
! 1156: bshiftz(md,1,&t);
! 1157: isqrtz(t,&s);
! 1158: bshiftz(s,64,&b);
! 1159: tmat = (Z **)mat->body;
! 1160: if ( UNIQ(md) ) return 0;
! 1161: row = mat->row; col = mat->col;
! 1162: dn0 = ONE;
! 1163: for ( i = 0; i < row; i++ )
! 1164: if ( cmpz(dn[i],dn0) > 0 ) dn0 = dn[i];
! 1165: w = (Z *)MALLOC(col*sizeof(Z));
! 1166: for ( i = 0; i < row; i++ )
! 1167: if ( stat[i] == 0 ) {
! 1168: for ( j = 0, tmi = tmat[i]; j < col; j++ )
! 1169: mulz(tmi[j],dn0,&w[j]);
! 1170: ret = intvtoratv(w,col,md,b,(Z *)BDY(nm)[i],&dn[i]);
! 1171: if ( ret ) {
! 1172: stat[i] = 1;
! 1173: mulz(dn0,dn[i],&t); dn[i] = t; dn0 = t;
! 1174: }
! 1175: }
! 1176: for ( i = 0; i < row; i++ ) if ( !stat[i] ) break;
! 1177: if ( i == row ) return 1;
! 1178: else return 0;
! 1179: }
! 1180:
! 1181: int intvtoratv(Z *v,int n,Z md,Z b,Z *nm,Z *dn)
! 1182: {
! 1183: Z dn0,dn1,q,s,u,nm1,unm,udn,dmy;
! 1184: Z *nmk;
! 1185: int j,l,col,ret,sgn;
! 1186:
! 1187: for ( j = 0; j < n; j++ ) nm[j] = 0;
! 1188: dn0 = ONE;
! 1189: for ( j = 0; j < n; j++ ) {
! 1190: if ( !v[j] ) continue;
! 1191: mulz(v[j],dn0,&s);
! 1192: divqrz(s,md,&dmy,&u);
! 1193: ret = inttorat(u,md,b,&nm1,&dn1);
! 1194: if ( !ret ) return 0;
! 1195: if ( !UNIQ(dn1) )
! 1196: for ( l = 0; l < j; l++ ) {
! 1197: mulz(nm[l],dn1,&q); nm[l] = q;
! 1198: }
! 1199: nm[j] = nm1;
! 1200: mulz(dn0,dn1,&q); dn0 = q;
! 1201: }
! 1202: *dn = dn0;
! 1203: return 1;
! 1204: }
! 1205:
! 1206: /* assuming 0 < c < m */
! 1207:
! 1208: int inttorat(Z c,Z m,Z b,Z *nmp,Z *dnp)
! 1209: {
! 1210: Z qq,t,u1,v1,r1;
! 1211: Z q,u2,v2,r2;
! 1212:
! 1213: u1 = 0; v1 = ONE; u2 = m; v2 = c;
! 1214: while ( cmpz(v2,b) >= 0 ) {
! 1215: divqrz(u2,v2,&q,&r2); u2 = v2; v2 = r2;
! 1216: mulz(q,v1,&t); subz(u1,t,&r1); u1 = v1; v1 = r1;
! 1217: }
! 1218: if ( cmpz(v1,b) >= 0 ) return 0;
! 1219: else {
! 1220: if ( mpz_sgn(BDY(v1))<0 ) {
! 1221: chsgnz(v1,dnp); chsgnz(v2,nmp);
! 1222: } else {
! 1223: *dnp = v1; *nmp = v2;
! 1224: }
! 1225: return 1;
! 1226: }
! 1227: }
! 1228:
! 1229: extern int f4_nocheck;
! 1230:
! 1231: int gensolve_check(MAT mat,MAT nm,Z dn,int *rind,int *cind)
! 1232: {
! 1233: int row,col,rank,clen,i,j,k,l;
! 1234: Z s,t;
! 1235: Z *w;
! 1236: Z *mati,*nmk;
! 1237:
! 1238: if ( f4_nocheck ) return 1;
! 1239: row = mat->row; col = mat->col; rank = nm->row; clen = nm->col;
! 1240: w = (Z *)MALLOC(clen*sizeof(Z));
! 1241: for ( i = 0; i < row; i++ ) {
! 1242: mati = (Z *)mat->body[i];
! 1243: bzero(w,clen*sizeof(Z));
! 1244: for ( k = 0; k < rank; k++ )
! 1245: for ( l = 0, nmk = (Z *)nm->body[k]; l < clen; l++ ) {
! 1246: mulz(mati[rind[k]],nmk[l],&t); addz(w[l],t,&s); w[l] = s;
! 1247: }
! 1248: for ( j = 0; j < clen; j++ ) {
! 1249: mulz(dn,mati[cind[j]],&t);
! 1250: if ( cmpz(w[j],t) ) break;
! 1251: }
! 1252: if ( j != clen ) break;
! 1253: }
! 1254: if ( i != row ) return 0;
! 1255: else return 1;
! 1256: }
! 1257:
! 1258: int gensolve_check2(MAT mat,MAT nm,Z *dn,int *rind,int *cind)
! 1259: {
! 1260: int row,col,rank,clen,i,j,k,l;
! 1261: Z s,t,u,d;
! 1262: Z *w,*m;
! 1263: Z *mati,*nmk;
! 1264:
! 1265: if ( f4_nocheck ) return 1;
! 1266: row = mat->row; col = mat->col; rank = nm->row; clen = nm->col;
! 1267: w = (Z *)MALLOC(clen*sizeof(Z));
! 1268: m = (Z *)MALLOC(clen*sizeof(Z));
! 1269: for ( d = dn[0], i = 1; i < rank; i++ ) {
! 1270: lcmz(d,dn[i],&t); d = t;
! 1271: }
! 1272: for ( i = 0; i < rank; i++ ) divsz(d,dn[i],&m[i]);
! 1273: for ( i = 0; i < row; i++ ) {
! 1274: mati = (Z *)mat->body[i];
! 1275: bzero(w,clen*sizeof(Z));
! 1276: for ( k = 0; k < rank; k++ ) {
! 1277: mulz(mati[rind[k]],m[k],&u);
! 1278: for ( l = 0, nmk = (Z *)nm->body[k]; l < clen; l++ ) {
! 1279: mulz(u,nmk[l],&t); addz(w[l],t,&s); w[l] = s;
! 1280: }
! 1281: }
! 1282: for ( j = 0; j < clen; j++ ) {
! 1283: mulz(d,mati[cind[j]],&t);
! 1284: if ( cmpz(w[j],t) ) break;
! 1285: }
! 1286: if ( j != clen ) break;
! 1287: }
! 1288: if ( i != row ) return 0;
! 1289: else return 1;
! 1290: }
! 1291:
! 1292: void isqrtz(Z a,Z *r)
! 1293: {
! 1294: int k;
! 1295: Z x,t,x2,xh,quo,rem;
! 1296: Z two;
! 1297:
! 1298: if ( !a ) *r = 0;
! 1299: else if ( UNIQ(a) ) *r = ONE;
! 1300: else {
! 1301: k = z_bits((Q)a); /* a <= 2^k-1 */
! 1302: bshiftz(ONE,-((k>>1)+(k&1)),&x); /* a <= x^2 */
! 1303: STOQ(2,two);
! 1304: while ( 1 ) {
! 1305: pwrz(x,two,&t);
! 1306: if ( cmpz(t,a) <= 0 ) {
! 1307: *r = x; return;
! 1308: } else {
! 1309: if ( mpz_tstbit(BDY(x),0) ) addz(x,a,&t);
! 1310: else t = a;
! 1311: bshiftz(x,-1,&x2); divqrz(t,x2,&quo,&rem);
! 1312: bshiftz(x,1,&xh); addz(quo,xh,&x);
! 1313: }
! 1314: }
! 1315: }
! 1316: }
! 1317:
! 1318: void bshiftz(Z a,int n,Z *r)
! 1319: {
! 1320: mpz_t t;
! 1321:
! 1322: if ( !a ) *r = 0;
! 1323: else if ( n == 0 ) *r = a;
! 1324: else if ( n < 0 ) {
! 1325: mpz_init(t); mpz_mul_2exp(t,BDY(a),-n); MPZTOZ(t,*r);
! 1326: } else {
! 1327: mpz_init(t); mpz_fdiv_q_2exp(t,BDY(a),n);
! 1328: if ( !mpz_sgn(t) ) *r = 0;
! 1329: else MPZTOZ(t,*r);
! 1330: }
! 1331: }
! 1332:
! 1333: void addlf(Z a,Z b,Z *c)
! 1334: {
! 1335: addz(a,b,c);
! 1336: if ( !lf_lazy ) {
! 1337: if ( cmpz(*c,current_mod_lf) >= 0 ) {
! 1338: subz(*c,current_mod_lf,c);
! 1339: }
! 1340: }
! 1341: }
! 1342:
! 1343: void sublf(Z a,Z b,Z *c)
! 1344: {
! 1345: subz(a,b,c);
! 1346: if ( !lf_lazy ) {
! 1347: remz(*c,current_mod_lf,c);
! 1348: }
! 1349: }
! 1350:
! 1351: void mullf(Z a,Z b,Z *c)
! 1352: {
! 1353: mulz(a,b,c);
! 1354: if ( !lf_lazy ) {
! 1355: remz(*c,current_mod_lf,c);
! 1356: }
! 1357: }
! 1358:
! 1359: void divlf(Z a,Z b,Z *c)
! 1360: {
! 1361: Z inv;
! 1362:
! 1363: invz(b,current_mod_lf,&inv);
! 1364: mulz(a,inv,c);
! 1365: if ( !lf_lazy ) {
! 1366: remz(*c,current_mod_lf,c);
! 1367: }
! 1368: }
! 1369:
! 1370: void chsgnlf(Z a,Z *c)
! 1371: {
! 1372: chsgnz(a,c);
! 1373: if ( !lf_lazy ) {
! 1374: remz(*c,current_mod_lf,c);
! 1375: }
! 1376: }
! 1377:
! 1378: void lmtolf(LM a,Z *b)
! 1379: {
! 1380: if ( !a ) *b = 0;
! 1381: else {
! 1382: MPZTOZ(BDY(a),*b);
! 1383: }
! 1384: }
! 1385:
! 1386: void setmod_lf(Z p)
! 1387: {
! 1388: current_mod_lf = p;
! 1389: current_mod_lf_size = mpz_size(BDY(current_mod_lf))+1;
! 1390: }
! 1391:
! 1392: void simplf_force(Z a,Z *b)
! 1393: {
! 1394: remz(a,current_mod_lf,b);
! 1395: }
! 1396:
! 1397: int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn,int **rindp,int **cindp)
! 1398: {
! 1399: MAT bmat,xmat;
! 1400: Z **a0,**a,**b,**x,**nm;
! 1401: Z *ai,*bi,*xi;
! 1402: int row,col;
! 1403: int **w;
! 1404: int *wi;
! 1405: int **wc;
! 1406: Z mdq,q,s,u;
! 1407: Z tn;
! 1408: int ind,md,i,j,k,l,li,ri,rank;
! 1409: unsigned int t;
! 1410: int *cinfo,*rinfo;
! 1411: int *rind,*cind;
! 1412: int count;
! 1413: int ret;
! 1414: struct oEGT eg_mul,eg_inv,eg_intrat,eg_check,tmp0,tmp1;
! 1415: int period;
! 1416: int *wx,*ptr;
! 1417: int wxsize,nsize;
! 1418: Z wn;
! 1419: Z wq;
! 1420:
! 1421: a0 = (Z **)mat->body;
! 1422: row = mat->row; col = mat->col;
! 1423: w = (int **)almat(row,col);
! 1424: for ( ind = 0; ; ind++ ) {
! 1425: md = get_lprime(ind);
! 1426: STOQ(md,mdq);
! 1427: for ( i = 0; i < row; i++ )
! 1428: for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
! 1429: wi[j] = remqi((Q)ai[j],md);
! 1430:
! 1431: if ( DP_Print > 3 ) {
! 1432: fprintf(asir_out,"LU decomposition.."); fflush(asir_out);
! 1433: }
! 1434: rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo);
! 1435: if ( DP_Print > 3 ) {
! 1436: fprintf(asir_out,"done.\n"); fflush(asir_out);
! 1437: }
! 1438: a = (Z **)almat_pointer(rank,rank); /* lhs mat */
! 1439: MKMAT(bmat,rank,col-rank); b = (Z **)bmat->body; /* lhs mat */
! 1440: for ( j = li = ri = 0; j < col; j++ )
! 1441: if ( cinfo[j] ) {
! 1442: /* the column is in lhs */
! 1443: for ( i = 0; i < rank; i++ ) {
! 1444: w[i][li] = w[i][j];
! 1445: a[i][li] = a0[rinfo[i]][j];
! 1446: }
! 1447: li++;
! 1448: } else {
! 1449: /* the column is in rhs */
! 1450: for ( i = 0; i < rank; i++ )
! 1451: b[i][ri] = a0[rinfo[i]][j];
! 1452: ri++;
! 1453: }
! 1454:
! 1455: /* solve Ax=B; A: rank x rank, B: rank x ri */
! 1456: /* algorithm
! 1457: c <- B
! 1458: x <- 0
! 1459: q <- 1
! 1460: do
! 1461: t <- A^(-1)c mod p
! 1462: x <- x+qt
! 1463: c <- (c-At)/p
! 1464: q <- qp
! 1465: end do
! 1466: then Ax-B=0 mod q and b=(B-Ax)/q hold after "do".
! 1467: */
! 1468: MKMAT(xmat,rank,ri); x = (Z **)(xmat)->body;
! 1469: MKMAT(*nmmat,rank,ri); nm = (Z **)(*nmmat)->body;
! 1470: wc = (int **)almat(rank,ri);
! 1471: for ( i = 0; i < rank; i++ )
! 1472: wc[i] = w[i]+rank;
! 1473: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
! 1474: *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));
! 1475:
! 1476: period = F4_INTRAT_PERIOD;
! 1477: for ( q = ONE, count = 0; ; ) {
! 1478: if ( DP_Print > 3 )
! 1479: fprintf(stderr,"o");
! 1480: /* wc = b mod md */
! 1481: for ( i = 0; i < rank; i++ )
! 1482: for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ ) {
! 1483: wi[j] = remqi((Q)bi[j],md);
! 1484: if ( wi[j] && sgnz(bi[j]) < 0 )
! 1485: wi[j] = md-wi[j];
! 1486: }
! 1487: /* wc = A^(-1)wc; wc is normalized */
! 1488: solve_by_lu_mod(w,rank,md,wc,ri,1);
! 1489: /* x += q*wc */
! 1490: for ( i = 0; i < rank; i++ )
! 1491: for ( j = 0, wi = wc[i]; j < ri; j++ ) mul1addtoz(q,wi[j],&x[i][j]);
! 1492: /* b =(A*wc+b)/md */
! 1493: for ( i = 0; i < rank; i++ )
! 1494: for ( j = 0; j < ri; j++ ) {
! 1495: u = b[i][j];
! 1496: for ( k = 0; k < rank; k++ ) mul1addtoz(a[i][k],wc[k][j],&u);
! 1497: divsz(u,mdq,&b[i][j]);
! 1498: }
! 1499: count++;
! 1500: /* q = q*md */
! 1501: mulz(q,mdq,&u); q = u;
! 1502: if ( count == period ) {
! 1503: ret = intmtoratm(xmat,q,*nmmat,dn);
! 1504: if ( ret ) {
! 1505: for ( j = k = l = 0; j < col; j++ )
! 1506: if ( cinfo[j] )
! 1507: rind[k++] = j;
! 1508: else
! 1509: cind[l++] = j;
! 1510: ret = gensolve_check(mat,*nmmat,*dn,rind,cind);
! 1511: if ( ret ) {
! 1512: *rindp = rind;
! 1513: *cindp = cind;
! 1514: for ( j = k = 0; j < col; j++ )
! 1515: if ( !cinfo[j] )
! 1516: cind[k++] = j;
! 1517: return rank;
! 1518: }
! 1519: } else {
! 1520: period = period*3/2;
! 1521: count = 0;
! 1522: }
! 1523: }
! 1524: }
! 1525: }
! 1526: }
! 1527:
! 1528: /* for inv_or_split_dalg */
! 1529:
! 1530: int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT *nmmat,Z *dn,int **rindp,int **cindp)
! 1531: {
! 1532: MAT bmat,xmat;
! 1533: Z **a0,**a,**b,**x,**nm;
! 1534: Z *ai,*bi,*xi;
! 1535: int row,col;
! 1536: int **w;
! 1537: int *wi;
! 1538: int **wc;
! 1539: Z mdq,q,s,u;
! 1540: Z tn;
! 1541: int ind,md,i,j,k,l,li,ri,rank;
! 1542: unsigned int t;
! 1543: int *cinfo,*rinfo;
! 1544: int *rind,*cind;
! 1545: int count;
! 1546: int ret;
! 1547: struct oEGT eg_mul,eg_inv,eg_intrat,eg_check,tmp0,tmp1;
! 1548: int period;
! 1549: int *wx,*ptr;
! 1550: int wxsize,nsize;
! 1551: Z wn;
! 1552: Z wq;
! 1553: DP m;
! 1554:
! 1555: a0 = (Z **)mat->body;
! 1556: row = mat->row; col = mat->col;
! 1557: w = (int **)almat(row,col);
! 1558: for ( ind = 0; ; ind++ ) {
! 1559: md = get_lprime(ind);
! 1560: STOQ(md,mdq);
! 1561: for ( i = 0; i < row; i++ )
! 1562: for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
! 1563: wi[j] = remqi((Q)ai[j],md);
! 1564:
! 1565: if ( DP_Print > 3 ) {
! 1566: fprintf(asir_out,"LU decomposition.."); fflush(asir_out);
! 1567: }
! 1568: rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo);
! 1569: if ( DP_Print > 3 ) {
! 1570: fprintf(asir_out,"done.\n"); fflush(asir_out);
! 1571: }
! 1572:
! 1573: /* this part is added for inv_or_split_dalg */
! 1574: for ( i = 0; i < col-1; i++ ) {
! 1575: if ( !cinfo[i] ) {
! 1576: m = mb[i];
! 1577: for ( j = i+1; j < col-1; j++ )
! 1578: if ( dp_redble(mb[j],m) )
! 1579: cinfo[j] = -1;
! 1580: }
! 1581: }
! 1582:
! 1583: a = (Z **)almat_pointer(rank,rank); /* lhs mat */
! 1584: MKMAT(bmat,rank,col-rank); b = (Z **)bmat->body; /* lhs mat */
! 1585: for ( j = li = ri = 0; j < col; j++ )
! 1586: if ( cinfo[j] ) {
! 1587: /* the column is in lhs */
! 1588: for ( i = 0; i < rank; i++ ) {
! 1589: w[i][li] = w[i][j];
! 1590: a[i][li] = a0[rinfo[i]][j];
! 1591: }
! 1592: li++;
! 1593: } else {
! 1594: /* the column is in rhs */
! 1595: for ( i = 0; i < rank; i++ )
! 1596: b[i][ri] = a0[rinfo[i]][j];
! 1597: ri++;
! 1598: }
! 1599:
! 1600: /* solve Ax=B; A: rank x rank, B: rank x ri */
! 1601: /* algorithm
! 1602: c <- B
! 1603: x <- 0
! 1604: q <- 1
! 1605: do
! 1606: t <- A^(-1)c mod p
! 1607: x <- x+qt
! 1608: c <- (c-At)/p
! 1609: q <- qp
! 1610: end do
! 1611: then Ax-B=0 mod q and b=(B-Ax)/q hold after "do".
! 1612: */
! 1613: MKMAT(xmat,rank,ri); x = (Z **)(xmat)->body;
! 1614: MKMAT(*nmmat,rank,ri); nm = (Z **)(*nmmat)->body;
! 1615: wc = (int **)almat(rank,ri);
! 1616: for ( i = 0; i < rank; i++ )
! 1617: wc[i] = w[i]+rank;
! 1618: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
! 1619: *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));
! 1620:
! 1621: period = F4_INTRAT_PERIOD;
! 1622: for ( q = ONE, count = 0; ; ) {
! 1623: if ( DP_Print > 3 )
! 1624: fprintf(stderr,"o");
! 1625: /* wc = b mod md */
! 1626: for ( i = 0; i < rank; i++ )
! 1627: for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ ) {
! 1628: wi[j] = remqi((Q)bi[j],md);
! 1629: if ( wi[j] && sgnz(bi[j]) < 0 )
! 1630: wi[j] = md-wi[j];
! 1631: }
! 1632: /* wc = A^(-1)wc; wc is normalized */
! 1633: solve_by_lu_mod(w,rank,md,wc,ri,1);
! 1634: /* x += q*wc */
! 1635: for ( i = 0; i < rank; i++ )
! 1636: for ( j = 0, wi = wc[i]; j < ri; j++ ) mul1addtoz(q,wi[j],&x[i][j]);
! 1637: /* b =(A*wc+b)/md */
! 1638: for ( i = 0; i < rank; i++ )
! 1639: for ( j = 0; j < ri; j++ ) {
! 1640: u = b[i][j];
! 1641: for ( k = 0; k < rank; k++ ) mul1addtoz(a[i][k],wc[k][j],&u);
! 1642: divsz(u,mdq,&b[i][j]);
! 1643: }
! 1644: count++;
! 1645: /* q = q*md */
! 1646: mulz(q,mdq,&u); q = u;
! 1647: if ( count == period ) {
! 1648: ret = intmtoratm(xmat,q,*nmmat,dn);
! 1649: if ( ret ) {
! 1650: for ( j = k = l = 0; j < col; j++ )
! 1651: if ( cinfo[j] > 0 )
! 1652: rind[k++] = j;
! 1653: else if ( !cinfo[j] )
! 1654: cind[l++] = j;
! 1655: ret = gensolve_check(mat,*nmmat,*dn,rind,cind);
! 1656: if ( ret ) {
! 1657: *rindp = rind;
! 1658: *cindp = cind;
! 1659: for ( j = k = 0; j < col; j++ )
! 1660: if ( !cinfo[j] )
! 1661: cind[k++] = j;
! 1662: return rank;
! 1663: }
! 1664: } else {
! 1665: period = period*3/2;
! 1666: count = 0;
! 1667: }
! 1668: }
! 1669: }
! 1670: }
! 1671: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to Q.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2018 / engine