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