Annotation of OpenXM_contrib2/asir2018/engine/Q.c, Revision 1.16
1.16 ! noro 1: /* $OpenXM: OpenXM_contrib2/asir2018/engine/Q.c,v 1.15 2019/04/23 08:07:24 noro Exp $ */
1.1 noro 2: #include "ca.h"
3: #include "gmp.h"
4: #include "base.h"
5: #include "inline.h"
6:
7: mpz_t ONEMPZ;
8: Z ONE;
9: int lf_lazy;
10: Z current_mod_lf;
11: int current_mod_lf_size;
12: gmp_randstate_t GMP_RAND;
13:
1.9 noro 14: #define F4_INTRAT_PERIOD 4
1.6 noro 15:
16: extern int DP_Print;
17:
1.1 noro 18: void isqrtz(Z a,Z *r);
19: void bshiftz(Z a,int n,Z *r);
20:
21: void *gc_realloc(void *p,size_t osize,size_t nsize)
22: {
23: return (void *)Risa_GC_realloc(p,nsize);
24: }
25:
26: void gc_free(void *p,size_t size)
27: {
28: Risa_GC_free(p);
29: }
30:
31: void init_gmpq()
32: {
1.10 noro 33: mp_set_memory_functions(Risa_GC_malloc,gc_realloc,gc_free);
1.1 noro 34:
35: mpz_init(ONEMPZ); mpz_set_ui(ONEMPZ,1); MPZTOZ(ONEMPZ,ONE);
36: gmp_randinit_default(GMP_RAND);
37: }
38:
1.7 noro 39: void printexpr(VL,Obj);
40:
1.3 noro 41: void pmat(Z **a,int row,int col)
42: {
43: int i,j;
44:
45: for ( i = 0; i < row; i++, printf("\n") )
46: for ( j = 0; j < col; j++, printf(" ") )
1.7 noro 47: printexpr(CO,(Obj)a[i][j]);
1.3 noro 48: printf("\n");
49: }
50:
1.1 noro 51: Z utoz(unsigned int u)
52: {
53: mpz_t z;
54: Z r;
55:
56: if ( !u ) return 0;
57: mpz_init(z); mpz_set_ui(z,u); MPZTOZ(z,r); return r;
58: }
59:
60: Z stoz(int s)
61: {
62: mpz_t z;
63: Z r;
64:
65: if ( !s ) return 0;
66: mpz_init(z); mpz_set_si(z,s); MPZTOZ(z,r); return r;
67: }
68:
69: int sgnz(Z z)
70: {
71: if ( !z ) return 0;
72: else return mpz_sgn(BDY(z));
73: }
74:
75: void nmq(Q q,Z *r)
76: {
77: if ( !q ) *r = 0;
78: else if ( INT(q) ) *r = (Z)q;
79: else {
80: MPZTOZ(mpq_numref(BDY(q)),*r);
81: }
82: }
83:
84: void dnq(Q q,Z *r)
85: {
86: if ( !q ) *r = 0;
87: else if ( INT(q) ) *r = ONE;
88: else {
89: MPZTOZ(mpq_denref(BDY(q)),*r);
90: }
91: }
92:
93: int sgnq(Q q)
94: {
95: if ( !q ) return 0;
96: else if ( q->z ) return mpz_sgn(BDY((Z)q));
97: else return mpz_sgn(mpq_numref(BDY(q)));
98: }
99:
100: Q mpqtozq(mpq_t a)
101: {
102: Z z;
103: Q q;
104:
105: if ( INTMPQ(a) ) {
106: MPZTOZ(mpq_numref(a),z); return (Q)z;
107: } else {
108: MPQTOQ(a,q); return q;
109: }
110: }
111:
112: void dupz(Z a,Z *b)
113: {
114: mpz_t t;
115:
116: if ( !a ) *b = a;
117: else {
118: mpz_init(t); mpz_set(t,BDY(a)); MPZTOZ(t,*b);
119: }
120: }
121:
122: int n_bits_z(Z a)
123: {
124: return a ? mpz_sizeinbase(BDY(a),2) : 0;
125: }
126:
127: void addz(Z n1,Z n2,Z *nr)
128: {
129: mpz_t t;
130: int s1,s2;
131:
132: if ( !n1 ) *nr = n2;
133: else if ( !n2 ) *nr = n1;
134: else if ( !n1->z || !n2->z )
135: error("addz : invalid argument");
136: else {
137: mpz_init(t); mpz_add(t,BDY(n1),BDY(n2)); MPZTOZ(t,*nr);
138: }
139: }
140:
141: void subz(Z n1,Z n2,Z *nr)
142: {
143: mpz_t t;
144:
145: if ( !n1 ) {
146: if ( !n2 )
147: *nr = 0;
148: else
149: chsgnz(n2,nr);
150: } else if ( !n2 )
151: *nr = n1;
152: else if ( n1 == n2 )
153: *nr = 0;
154: else if ( !n1->z || !n2->z )
155: error("subz : invalid argument");
156: else {
157: mpz_init(t); mpz_sub(t,BDY(n1),BDY(n2)); MPZTOZ(t,*nr);
158: }
159: }
160:
161: void mulz(Z n1,Z n2,Z *nr)
162: {
163: mpz_t t;
164:
165: if ( !n1 || !n2 ) *nr = 0;
166: else if ( !n1->z || !n2->z )
167: error("mulz : invalid argument");
168: else if ( UNIQ(n1) ) *nr = n2;
169: else if ( UNIQ(n2) ) *nr = n1;
170: else if ( MUNIQ(n1) ) chsgnz(n2,nr);
171: else if ( MUNIQ(n2) ) chsgnz(n1,nr);
172: else {
173: mpz_init(t); mpz_mul(t,BDY(n1),BDY(n2)); MPZTOZ(t,*nr);
174: }
175: }
176:
177: /* nr += n1*n2 */
178:
179: void muladdtoz(Z n1,Z n2,Z *nr)
180: {
1.3 noro 181: #if 0
1.1 noro 182: Z t;
183:
184: if ( n1 && n2 ) {
185: if ( !(*nr) ) {
186: NEWZ(t); mpz_init(BDY(t)); *nr = t;
187: }
188: mpz_addmul(BDY(*nr),BDY(n1),BDY(n2));
1.2 noro 189: if ( !mpz_sgn(BDY(*nr)) )
190: *nr = 0;
1.3 noro 191: }
1.2 noro 192: #else
193: Z t,s;
194:
195: mulz(n1,n2,&t); addz(*nr,t,&s); *nr = s;
196: #endif
1.1 noro 197: }
198:
199: /* nr += n1*u */
200:
201: void mul1addtoz(Z n1,long u,Z *nr)
202: {
1.3 noro 203: #if 0
1.1 noro 204: Z t;
205:
206: if ( n1 && u ) {
207: if ( !(*nr) ) {
208: NEWZ(t); mpz_init(BDY(t)); *nr = t;
209: }
210: if ( u >= 0 )
211: mpz_addmul_ui(BDY(*nr),BDY(n1),(unsigned long)u);
212: else
213: mpz_submul_ui(BDY(*nr),BDY(n1),(unsigned long)(-u));
1.2 noro 214: if ( !mpz_sgn(BDY(*nr)) )
215: *nr = 0;
1.1 noro 216: }
1.3 noro 217: #else
218: Z t,s;
219:
220: mul1z(n1,u,&t); addz(*nr,t,&s); *nr = s;
221: #endif
1.1 noro 222: }
223:
224: void mul1z(Z n1,long n2,Z *nr)
225: {
226: mpz_t t;
227:
228: if ( !n1 || !n2 ) *nr = 0;
229: else {
230: mpz_init(t); mpz_mul_si(t,BDY(n1),n2); MPZTOZ(t,*nr);
231: }
232: }
233:
234: void divz(Z n1,Z n2,Z *nq)
235: {
236: mpz_t t;
237: mpq_t a, b, q;
238:
239: if ( !n2 ) {
240: error("division by 0");
241: *nq = 0;
242: } else if ( !n1 )
243: *nq = 0;
244: else if ( n1 == n2 ) {
245: mpz_init(t); mpz_set_ui(t,1); MPZTOZ(t,*nq);
246: } else {
247: MPZTOMPQ(BDY(n1),a); MPZTOMPQ(BDY(n2),b);
248: mpq_init(q); mpq_div(q,a,b); *nq = (Z)mpqtozq(q);
249: }
250: }
251:
252: void remz(Z n1,Z n2,Z *nr)
253: {
254: mpz_t r;
255:
256: if ( !n2 ) {
257: error("division by 0");
258: *nr = 0;
259: } else if ( !n1 || n1 == n2 )
260: *nr = 0;
261: else if ( !n1->z || !n2->z )
262: error("remz : invalid argument");
263: else {
264: mpz_init(r);
265: mpz_mod(r,BDY(n1),BDY(n2));
266: if ( !mpz_sgn(r) ) *nr = 0;
267: else MPZTOZ(r,*nr);
268: }
269: }
270:
271: void divqrz(Z n1,Z n2,Z *nq,Z *nr)
272: {
273: mpz_t t, a, b, q, r;
274:
275: if ( !n2 ) {
276: error("division by 0");
277: *nq = 0; *nr = 0;
278: } else if ( !n1 ) {
279: *nq = 0; *nr = 0;
280: } else if ( !n1->z || !n2->z )
281: error("divqrz : invalid argument");
282: else if ( n1 == n2 ) {
283: mpz_init(t); mpz_set_ui(t,1); MPZTOZ(t,*nq); *nr = 0;
284: } else {
285: mpz_init(q); mpz_init(r);
286: mpz_fdiv_qr(q,r,BDY(n1),BDY(n2));
287: if ( !mpz_sgn(q) ) *nq = 0;
288: else MPZTOZ(q,*nq);
289: if ( !mpz_sgn(r) ) *nr = 0;
290: else MPZTOZ(r,*nr);
291: }
292: }
293:
294: void divsz(Z n1,Z n2,Z *nq)
295: {
296: mpz_t t;
297: mpq_t a, b, q;
298:
299: if ( !n2 ) {
300: error("division by 0");
301: *nq = 0;
302: } else if ( !n1 )
303: *nq = 0;
304: else if ( !n1->z || !n2->z )
305: error("divsz : invalid argument");
306: else if ( n1 == n2 ) {
307: mpz_init(t); mpz_set_ui(t,1); MPZTOZ(t,*nq);
308: } else {
309: mpz_init(t); mpz_divexact(t,BDY(n1),BDY(n2)); MPZTOZ(t,*nq);
310: }
311: }
312:
313: void chsgnz(Z n,Z *nr)
314: {
315: mpz_t t;
316:
317: if ( !n )
318: *nr = 0;
319: else if ( !n->z )
320: error("chsgnz : invalid argument");
321: else {
322: t[0] = BDY(n)[0]; mpz_neg(t,t); MPZTOZ(t,*nr);
323: }
324: }
325:
326: void absz(Z n,Z *nr)
327: {
328: if ( !n ) *nr = 0;
329: else if ( !n->z )
330: error("absz : invalid argument");
331: else if ( sgnz(n) < 0 ) chsgnz(n,nr);
332: else *nr = n;
333: }
334:
335: int evenz(Z n)
336: {
337: return !n ? 1 : mpz_even_p(BDY(n));
338: }
339:
340: int smallz(Z n)
341: {
342: if ( !n ) return 1;
343: else if ( INT(n) && mpz_fits_sint_p(BDY(n)) ) return 1;
344: else return 0;
345: }
346:
347: void pwrz(Z n1,Z n,Z *nr)
348: {
349: mpq_t t,q;
350: mpz_t z;
351: Q p,r;
352:
353: if ( !n || UNIQ(n1) ) *nr = ONE;
354: else if ( !n1 ) *nr = 0;
355: else if ( !n->z || !n1->z )
356: error("pwrz : invalid argument");
357: else if ( MUNIQ(n1) ) {
358: if ( mpz_even_p(BDY((Z)n)) ) *nr = ONE;
359: else *nr = n1;
360: } else if ( !smallz(n) ) {
361: error("exponent too big."); *nr = 0;
362: } else if ( n1->z && mpz_sgn(BDY((Z)n))>0 ) {
1.5 noro 363: mpz_init(z); mpz_pow_ui(z,BDY(n1),ZTOS(n)); MPZTOZ(z,*nr);
1.1 noro 364: } else {
365: MPZTOMPQ(BDY(n1),q); MPQTOQ(q,r);
366: pwrq(r,(Q)n,&p); *nr = (Z)p;
367: }
368: }
369:
370: int cmpz(Z q1,Z q2)
371: {
372: int sgn;
373:
374: if ( !q1 ) {
375: if ( !q2 )
376: return 0;
377: else
378: return -mpz_sgn(BDY(q2));
379: } else if ( !q2 )
380: return mpz_sgn(BDY(q1));
381: else if ( !q1->z || !q2->z )
382: error("mpqz : invalid argument");
383: else if ( (sgn = mpz_sgn(BDY(q1))) != mpz_sgn(BDY(q2)) )
384: return sgn;
385: else {
386: sgn = mpz_cmp(BDY(q1),BDY(q2));
387: if ( sgn > 0 ) return 1;
388: else if ( sgn < 0 ) return -1;
389: else return 0;
390: }
391: }
392:
393: void gcdz(Z n1,Z n2,Z *nq)
394: {
395: mpz_t t;
396:
397: if ( !n1 ) *nq = n2;
398: else if ( !n2 ) *nq = n1;
399: else if ( !n1->z || !n2->z )
400: error("gcdz : invalid argument");
401: else {
402: mpz_init(t); mpz_gcd(t,BDY(n1),BDY(n2));
403: MPZTOZ(t,*nq);
404: }
405: }
406:
407: void invz(Z n1,Z n2,Z *nq)
408: {
409: mpz_t t;
410:
411: if ( !n1 || !n2 || !n1->z || !n2->z )
412: error("invz : invalid argument");
413: mpz_init(t); mpz_invert(t,BDY(n1),BDY(n2));
414: MPZTOZ(t,*nq);
415: }
416:
417: void lcmz(Z n1,Z n2,Z *nq)
418: {
419: Z g,t;
420:
421: if ( !n1 || !n2 ) *nq = 0;
422: else if ( !n1->z || !n2->z )
423: error("lcmz : invalid argument");
424: else {
425: gcdz(n1,n2,&g); divsz(n1,g,&t);
426: mulz(n2,t,nq);
427: }
428: }
429:
430: void gcdvz(VECT v,Z *q)
431: {
432: int n,i;
433: Z *b;
434: Z g,g1;
435:
436: n = v->len;
437: b = (Z *)v->body;
438: g = b[0];
439: for ( i = 1; i < n; i++ ) {
440: gcdz(g,b[i],&g1); g = g1;
441: }
442: *q = g;
443: }
444:
445: void gcdvz_estimate(VECT v,Z *q)
446: {
447: int n,m,i;
448: Z s,t,u;
449: Z *b;
450:
451: n = v->len;
452: b = (Z *)v->body;
453: if ( n == 1 ) {
454: if ( mpz_sgn(BDY(b[0]))<0 ) chsgnz(b[0],q);
455: else *q = b[0];
456: }
457: m = n/2;
458: for ( i = 0, s = 0; i < m; i++ ) {
459: if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subz(s,b[i],&u);
460: else addz(s,b[i],&u);
461: s = u;
462: }
1.4 noro 463: for ( t = 0; i < n; i++ ) {
1.1 noro 464: if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subz(t,b[i],&u);
465: else addz(t,b[i],&u);
466: t = u;
467: }
468: gcdz(s,t,q);
469: }
470:
1.4 noro 471: void gcdv_mpz_estimate(mpz_t g,mpz_t *b,int n)
472: {
473: int m,m2,i,j;
474: mpz_t s,t;
475:
476: mpz_init(g);
477: for ( i = 0, m = 0; i < n; i++ )
478: if ( mpz_sgn(b[i]) ) m++;
479: if ( !m ) {
480: mpz_set_ui(g,0);
481: return;
482: }
483: if ( m == 1 ) {
484: for ( i = 0, m = 0; i < n; i++ )
485: if ( mpz_sgn(b[i]) ) break;
486: if ( mpz_sgn(b[i])<0 ) mpz_neg(g,b[i]);
487: else mpz_set(g,b[i]);
488: return ;
489: }
490: m2 = m/2;
491: mpz_init_set_ui(s,0);
492: for ( i = j = 0; j < m2; i++ ) {
493: if ( mpz_sgn(b[i]) ) {
494: if ( mpz_sgn(b[i])<0 )
495: mpz_sub(s,s,b[i]);
496: else
497: mpz_add(s,s,b[i]);
498: j++;
499: }
500: }
501: mpz_init_set_ui(t,0);
502: for ( ; i < n; i++ ) {
503: if ( mpz_sgn(b[i]) ) {
504: if ( mpz_sgn(b[i])<0 )
505: mpz_sub(t,t,b[i]);
506: else
507: mpz_add(t,t,b[i]);
508: }
509: }
510: mpz_gcd(g,s,t);
511: }
512:
513:
1.1 noro 514: void factorialz(unsigned int n,Z *nr)
515: {
516: mpz_t a;
517: mpz_init(a);
1.13 noro 518: mpz_fac_ui(a,(unsigned long)n);
1.1 noro 519: MPZTOZ(a,*nr);
520: }
521:
522: void randomz(int blen,Z *nr)
523: {
524: mpz_t z;
525:
526: mpz_init(z);
527: mpz_urandomb(z,GMP_RAND,blen);
528: MPZTOZ(z,*nr);
529: }
530:
531: int tstbitz(Z n,int k)
532: {
533: if ( !n || !n->z )
534: error("tstbitz : invalid argument");
535: return !n ? 0 : mpz_tstbit(BDY(n),k);
536: }
537:
538: void addq(Q n1,Q n2,Q *nr)
539: {
540: mpq_t q1,q2,t;
541:
542: if ( !n1 ) *nr = n2;
543: else if ( !n2 ) *nr = n1;
544: else if ( n1->z && n2->z )
545: addz((Z)n1,(Z)n2,(Z *)nr);
546: else {
547: if ( n1->z ) MPZTOMPQ(BDY((Z)n1),q1);
548: else q1[0] = BDY(n1)[0];
549: if ( n2->z ) MPZTOMPQ(BDY((Z)n2),q2);
550: else q2[0] = BDY(n2)[0];
551: mpq_init(t); mpq_add(t,q1,q2); *nr = mpqtozq(t);
552: }
553: }
554:
555: void subq(Q n1,Q n2,Q *nr)
556: {
557: mpq_t q1,q2,t;
558:
559: if ( !n1 ) {
560: if ( !n2 ) *nr = 0;
1.15 noro 561: else if ( n2->z ) chsgnz((Z)n2,(Z *)nr);
1.1 noro 562: else {
563: mpq_init(t); mpq_neg(t,BDY(n2)); MPQTOQ(t,*nr);
564: }
565: } else if ( !n2 ) *nr = n1;
566: else if ( n1 == n2 ) *nr = 0;
567: else if ( n1->z && n2->z )
568: subz((Z)n1,(Z)n2,(Z *)nr);
569: else {
570: if ( n1->z ) MPZTOMPQ(BDY((Z)n1),q1);
571: else q1[0] = BDY(n1)[0];
572: if ( n2->z ) MPZTOMPQ(BDY((Z)n2),q2);
573: else q2[0] = BDY(n2)[0];
574: mpq_init(t); mpq_sub(t,q1,q2); *nr = mpqtozq(t);
575: }
576: }
577:
578: void mulq(Q n1,Q n2,Q *nr)
579: {
580: mpq_t t,q1,q2;
581:
582: if ( !n1 || !n2 ) *nr = 0;
583: else if ( n1->z && n2->z )
584: mulz((Z)n1,(Z)n2,(Z *)nr);
585: else {
586: if ( n1->z ) MPZTOMPQ(BDY((Z)n1),q1);
587: else q1[0] = BDY(n1)[0];
588: if ( n2->z ) MPZTOMPQ(BDY((Z)n2),q2);
589: else q2[0] = BDY(n2)[0];
590: mpq_init(t); mpq_mul(t,q1,q2); *nr = mpqtozq(t);
591: }
592: }
593:
594: void divq(Q n1,Q n2,Q *nq)
595: {
596: mpq_t t,q1,q2;
597:
598: if ( !n2 ) {
599: error("division by 0");
600: *nq = 0;
601: return;
602: } else if ( !n1 ) *nq = 0;
603: else if ( n1 == n2 ) *nq = (Q)ONE;
604: else {
605: if ( n1->z ) MPZTOMPQ(BDY((Z)n1),q1);
606: else q1[0] = BDY(n1)[0];
607: if ( n2->z ) MPZTOMPQ(BDY((Z)n2),q2);
608: else q2[0] = BDY(n2)[0];
609: mpq_init(t); mpq_div(t,q1,q2); *nq = mpqtozq(t);
610: }
611: }
612:
613: void invq(Q n,Q *nr)
614: {
615: Z nm,dn;
616:
617: if ( INT(n) )
618: divq((Q)ONE,n,nr);
619: else {
620: nmq(n,&nm);
621: dnq(n,&dn);
622: divq((Q)dn,(Q)nm,nr);
623: }
624: }
625:
626: void chsgnq(Q n,Q *nr)
627: {
628: mpq_t t;
629:
630: if ( !n ) *nr = 0;
631: else if (n->z ) chsgnz((Z)n,(Z *)nr);
632: else {
633: mpq_init(t); mpq_neg(t,BDY(n)); MPQTOQ(t,*nr);
634: }
635: }
636:
637: void absq(Q n,Q *nr)
638: {
639: if ( !n ) *nr = 0;
640: else if ( n->z ) absz((Z)n,(Z *)nr);
641: else if ( sgnq(n) < 0 ) chsgnq(n,nr);
642: else *nr = n;
643: }
644:
645: void pwrq(Q n1,Q n,Q *nr)
646: {
647: int e;
648: mpz_t nm,dn;
649: mpq_t t;
650:
651: if ( !n || UNIQ((Z)n1) || UNIQ(n1) ) *nr = (Q)ONE;
652: else if ( !n1 ) *nr = 0;
653: else if ( !INT(n) ) {
654: error("can't calculate fractional power."); *nr = 0;
655: } else if ( !smallz((Z)n) ) {
656: error("exponent too big."); *nr = 0;
657: } else {
1.5 noro 658: e = ZTOS(n);
1.1 noro 659: if ( e < 0 ) {
660: e = -e;
661: if ( n1->z ) {
662: nm[0] = ONEMPZ[0];
663: dn[0] = BDY((Z)n1)[0];
664: } else {
665: nm[0] = mpq_denref(BDY(n1))[0];
666: dn[0] = mpq_numref(BDY(n1))[0];
667: }
668: } else {
669: if ( n1->z ) {
670: nm[0] = BDY((Z)n1)[0];
671: dn[0] = ONEMPZ[0];
672: } else {
673: nm[0] = mpq_numref(BDY(n1))[0];
674: dn[0] = mpq_denref(BDY(n1))[0];
675: }
676: }
677: mpq_init(t);
678: mpz_pow_ui(mpq_numref(t),nm,e); mpz_pow_ui(mpq_denref(t),dn,e);
679: *nr = mpqtozq(t);
680: }
681: }
682:
683: int cmpq(Q n1,Q n2)
684: {
685: mpq_t q1,q2;
686: int sgn;
687:
688: if ( !n1 ) {
689: if ( !n2 ) return 0;
690: else return (n2->z) ? -mpz_sgn(BDY((Z)n2)) : -mpq_sgn(BDY(n2));
691: } if ( !n2 ) return (n1->z) ? mpz_sgn(BDY((Z)n1)) : mpq_sgn(BDY(n1));
692: else if ( n1->z && n2->z )
693: return cmpz((Z)n1,(Z)n2);
694: else if ( (sgn = mpq_sgn(BDY(n1))) != mpq_sgn(BDY(n2)) ) return sgn;
695: else {
696: if ( n1->z ) MPZTOMPQ(BDY((Z)n1),q1);
697: else q1[0] = BDY(n1)[0];
698: if ( n2->z ) MPZTOMPQ(BDY((Z)n2),q2);
699: else q2[0] = BDY(n2)[0];
700: sgn = mpq_cmp(q1,q2);
701: if ( sgn > 0 ) return 1;
702: else if ( sgn < 0 ) return -1;
703: else return 0;
704: }
705: }
706:
707: /* t = [nC0 nC1 ... nCn] */
708:
709: void mkbc(int n,Z *t)
710: {
711: int i;
712: Z c,d,iq;
713:
714: for ( t[0] = ONE, i = 1; i <= n/2; i++ ) {
1.5 noro 715: STOZ(n-i+1,c); mulz(t[i-1],c,&d);
716: STOZ(i,iq); divsz(d,iq,&t[i]);
1.1 noro 717: }
718: for ( ; i <= n; i++ )
719: t[i] = t[n-i];
720: }
721:
722: /*
723: * Dx^k*x^l = W(k,l,0)*x^l*Dx^k+W(k,l,1)*x^(l-1)*x^(k-1)*+...
724: *
725: * t = [W(k,l,0) W(k,l,1) ... W(k,l,min(k,l)]
726: * where W(k,l,i) = i! * kCi * lCi
727: */
728:
729: /* mod m table */
730: /* XXX : should be optimized */
731:
732: void mkwcm(int k,int l,int m,int *t)
733: {
734: int i,n;
735: Z *s;
736:
737: n = MIN(k,l);
738: s = (Z *)ALLOCA((n+1)*sizeof(Q));
739: mkwc(k,l,s);
740: for ( i = 0; i <= n; i++ ) {
741: t[i] = remqi((Q)s[i],m);
742: }
743: }
744:
745: void mkwc(int k,int l,Z *t)
746: {
747: mpz_t a,b,q,nm,z,u;
748: int i,n;
749:
750: n = MIN(k,l);
751: mpz_init_set_ui(z,1);
752: mpz_init(u); mpz_set(u,z); MPZTOZ(u,t[0]);
753: mpz_init(a); mpz_init(b); mpz_init(nm);
754: for ( i = 1; i <= n; i++ ) {
755: mpz_set_ui(a,k-i+1); mpz_set_ui(b,l-i+1); mpz_mul(nm,a,b);
756: mpz_mul(z,BDY(t[i-1]),nm); mpz_fdiv_q_ui(z,z,i);
757: mpz_init(u); mpz_set(u,z); MPZTOZ(u,t[i]);
758: }
759: }
760:
761: void lgp(P p,Z *g,Z *l);
762:
763: void ptozp(P p,int sgn,Q *c,P *pr)
764: {
1.16 ! noro 765: Z nm,dn,nm1;
1.1 noro 766:
767: if ( !p ) {
768: *c = 0; *pr = 0;
769: } else {
770: lgp(p,&nm,&dn);
1.16 ! noro 771: if ( sgn < 0 ) {
! 772: chsgnz(nm,&nm1); nm = nm1;
! 773: }
1.1 noro 774: divz(nm,dn,(Z *)c);
775: divsp(CO,p,(P)*c,pr);
776: }
777: }
778:
779: void lgp(P p,Z *g,Z *l)
780: {
781: DCP dc;
782: Z g1,g2,l1,l2,l3,l4;
783:
784: if ( NUM(p) ) {
785: if ( ((Q)p)->z ) {
786: MPZTOZ(BDY((Z)p),*g);
787: *l = ONE;
788: } else {
789: MPZTOZ(mpq_numref(BDY((Q)p)),*g);
790: MPZTOZ(mpq_denref(BDY((Q)p)),*l);
791: }
792: } else {
793: dc = DC(p); lgp(COEF(dc),g,l);
794: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
795: lgp(COEF(dc),&g1,&l1); gcdz(*g,g1,&g2); *g = g2;
796: gcdz(*l,l1,&l2); mulz(*l,l1,&l3); divz(l3,l2,l);
797: }
798: }
799: }
800:
801: void qltozl(Q *w,int n,Z *dvr)
802: {
803: Z nm,dn;
804: Z g,g1,l1,l2,l3;
805: Q c;
806: int i;
807: struct oVECT v;
808:
809: for ( i = 0; i < n; i++ )
810: if ( w[i] && !w[i]->z )
811: break;
812: if ( i == n ) {
813: v.id = O_VECT; v.len = n; v.body = (pointer *)w;
814: gcdvz(&v,dvr); return;
815: }
816: for ( i = 0; !w[i]; i++ );
817: c = w[i];
818: if ( !c->z ) {
819: MPZTOZ(mpq_numref(BDY(c)),nm); MPZTOZ(mpq_denref(BDY(c)),dn);
820: } else {
821: MPZTOZ(BDY((Z)c),nm); dn = ONE;
822: }
823: for ( i++; i < n; i++ ) {
824: c = w[i];
825: if ( !c ) continue;
826: if ( !c->z ) {
827: MPZTOZ(mpq_numref(BDY(c)),g1); MPZTOZ(mpq_denref(BDY(c)),l1);
828: } else {
829: MPZTOZ(BDY((Z)c),g1); l1 = ONE;
830: }
831: gcdz(nm,g1,&g); nm = g;
832: gcdz(dn,l1,&l2); mulz(dn,l1,&l3); divz(l3,l2,&dn);
833: }
834: divz(nm,dn,dvr);
835: }
836:
837: int z_bits(Q q)
838: {
839: if ( !q ) return 0;
840: else if ( q->z ) return mpz_sizeinbase(BDY((Z)q),2);
841: else
842: return mpz_sizeinbase(mpq_numref(BDY(q)),2)
843: + mpz_sizeinbase(mpq_denref(BDY(q)),2);
844: }
845:
846: int zp_mag(P p)
847: {
848: int s;
849: DCP dc;
850:
851: if ( !p ) return 0;
852: else if ( OID(p) == O_N ) return z_bits((Q)p);
853: else {
854: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) s += zp_mag(COEF(dc));
855: return s;
856: }
857: }
858:
859: void makesubstz(VL v,NODE *s)
860: {
861: NODE r,r0;
862: Z q;
863: unsigned int n;
864:
865: for ( r0 = 0; v; v = NEXT(v) ) {
866: NEXTNODE(r0,r); BDY(r) = (pointer)v->v;
867: #if defined(_PA_RISC1_1)
868: n = mrand48()&BMASK; q = utoz(n);
869: #else
870: n = random(); q = utoz(n);
871: #endif
872: NEXTNODE(r0,r); BDY(r) = (pointer)q;
873: }
874: if ( r0 ) NEXT(r) = 0;
875: *s = r0;
876: }
877:
878: unsigned int remqi(Q a,unsigned int mod)
879: {
880: unsigned int c,nm,dn;
881: mpz_t r;
882:
883: if ( !a ) return 0;
884: else if ( a->z ) {
885: mpz_init(r);
886: c = mpz_fdiv_r_ui(r,BDY((Z)a),mod);
887: } else {
888: mpz_init(r);
889: nm = mpz_fdiv_r_ui(r,mpq_numref(BDY(a)),mod);
890: dn = mpz_fdiv_r_ui(r,mpq_denref(BDY(a)),mod);
891: dn = invm(dn,mod);
892: DMAR(nm,dn,0,mod,c);
893: }
894: return c;
895: }
896:
897: int generic_gauss_elim(MAT mat,MAT *nm,Z *dn,int **rindp,int **cindp)
898: {
899: int **wmat;
900: Z **bmat,**tmat,*bmi,*tmi;
901: Z q,m1,m2,m3,s,u;
902: int *wmi,*colstat,*wcolstat,*rind,*cind;
903: int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv;
904: MAT r,crmat;
905: int ret;
1.8 noro 906: MAT mat2,nm2;
907: Z dn2;
908: int *rind2,*cind2;
909: int ret2;
1.1 noro 910:
1.6 noro 911: #if SIZEOF_LONG == 8
1.8 noro 912: ret = generic_gauss_elim64(mat,nm,dn,rindp,cindp);
913: return ret;
1.6 noro 914: #endif
1.1 noro 915: bmat = (Z **)mat->body;
916: row = mat->row; col = mat->col;
917: wmat = (int **)almat(row,col);
918: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
919: wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
920: for ( ind = 0; ; ind++ ) {
921: if ( DP_Print ) {
922: fprintf(asir_out,"."); fflush(asir_out);
923: }
924: md = get_lprime(ind);
925: for ( i = 0; i < row; i++ )
926: for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
927: wmi[j] = remqi((Q)bmi[j],md);
928: rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);
929: if ( !ind ) {
930: RESET:
931: m1 = utoz(md);
932: rank0 = rank;
933: bcopy(wcolstat,colstat,col*sizeof(int));
934: MKMAT(crmat,rank,col-rank);
935: MKMAT(r,rank,col-rank); *nm = r;
936: tmat = (Z **)crmat->body;
937: for ( i = 0; i < rank; i++ )
938: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
939: if ( !colstat[j] ) tmi[k++] = utoz(wmi[j]);
940: } else {
941: if ( rank < rank0 ) {
942: if ( DP_Print ) {
943: fprintf(asir_out,"lower rank matrix; continuing...\n");
944: fflush(asir_out);
945: }
946: continue;
947: } else if ( rank > rank0 ) {
948: if ( DP_Print ) {
949: fprintf(asir_out,"higher rank matrix; resetting...\n");
950: fflush(asir_out);
951: }
952: goto RESET;
953: } else {
954: for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
955: if ( j < col ) {
956: if ( DP_Print ) {
957: fprintf(asir_out,"inconsitent colstat; resetting...\n");
958: fflush(asir_out);
959: }
960: goto RESET;
961: }
962: }
963:
964: inv = invm(remqi((Q)m1,md),md);
965: m2 = utoz(md); mulz(m1,m2,&m3);
966: for ( i = 0; i < rank; i++ )
967: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
968: if ( !colstat[j] ) {
969: if ( tmi[k] ) {
970: /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
971: t = remqi((Q)tmi[k],md);
972: if ( wmi[j] >= t ) t = wmi[j]-t;
973: else t = md-(t-wmi[j]);
974: DMAR(t,inv,0,md,t1)
975: u = utoz(t1); mulz(m1,u,&s);
976: addz(tmi[k],s,&u); tmi[k] = u;
977: } else if ( wmi[j] ) {
978: /* f3 = m1*(m1 mod m2)^(-1)*f2 */
979: DMAR(wmi[j],inv,0,md,t)
980: u = utoz(t); mulz(m1,u,&s); tmi[k] = s;
981: }
982: k++;
983: }
984: m1 = m3;
985: if ( ind % F4_INTRAT_PERIOD )
986: ret = 0;
987: else
988: ret = intmtoratm(crmat,m1,*nm,dn);
989: if ( ret ) {
990: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
991: *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
992: for ( j = k = l = 0; j < col; j++ )
993: if ( colstat[j] ) rind[k++] = j;
994: else cind[l++] = j;
995: if ( gensolve_check(mat,*nm,*dn,rind,cind) )
996: return rank;
997: }
998: }
999: }
1000: }
1001:
1002: int generic_gauss_elim2(MAT mat,MAT *nm,Z *dn,int **rindp,int **cindp)
1003: {
1004:
1005: MAT full;
1006: Z **bmat,**b;
1007: Z *bmi;
1008: Z dn0;
1009: int row,col,md,i,j,rank,ret;
1010: int **wmat;
1011: int *wmi;
1012: int *colstat,*rowstat;
1013:
1014: bmat = (Z **)mat->body;
1015: row = mat->row; col = mat->col;
1016: wmat = (int **)almat(row,col);
1017: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
1018: rowstat = (int *)MALLOC_ATOMIC(row*sizeof(int));
1019: /* XXX */
1020: md = get_lprime(0);
1021: for ( i = 0; i < row; i++ )
1022: for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
1023: wmi[j] = remqi((Q)bmi[j],md);
1024: rank = generic_gauss_elim_mod2(wmat,row,col,md,colstat,rowstat);
1025: b = (Z **)MALLOC(rank*sizeof(Z));
1026: for ( i = 0; i < rank; i++ ) b[i] = bmat[rowstat[i]];
1027: NEWMAT(full); full->row = rank; full->col = col; full->body = (pointer **)b;
1028: ret = generic_gauss_elim_full(full,nm,dn,rindp,cindp);
1029: if ( !ret ) {
1030: rank = generic_gauss_elim(mat,nm,&dn0,rindp,cindp);
1031: for ( i = 0; i < rank; i++ ) dn[i] = dn0;
1032: }
1033: return rank;
1034: }
1035:
1036: int generic_gauss_elim_full(MAT mat,MAT *nm,Z *dn,int **rindp,int **cindp)
1037: {
1038: int **wmat;
1039: int *stat;
1040: Z **bmat,**tmat,*bmi,*tmi,*ri;
1041: Z q,m1,m2,m3,s,u;
1042: int *wmi,*colstat,*wcolstat,*rind,*cind;
1043: int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv,h;
1044: MAT r,crmat;
1045: int ret,initialized,done;
1046:
1047: initialized = 0;
1048: bmat = (Z **)mat->body;
1049: row = mat->row; col = mat->col;
1050: wmat = (int **)almat(row,col);
1051: stat = (int *)MALLOC_ATOMIC(row*sizeof(int));
1052: for ( i = 0; i < row; i++ ) stat[i] = 0;
1053: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
1054: wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
1055: for ( ind = 0; ; ind++ ) {
1056: if ( DP_Print ) {
1057: fprintf(asir_out,"."); fflush(asir_out);
1058: }
1059: md = get_lprime(ind);
1060: for ( i = 0; i < row; i++ )
1061: for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
1062: wmi[j] = remqi((Q)bmi[j],md);
1063: rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);
1064: if ( rank < row ) continue;
1065: if ( !initialized ) {
1066: m1 = utoz(md);
1067: bcopy(wcolstat,colstat,col*sizeof(int));
1068: MKMAT(crmat,row,col-row);
1069: MKMAT(r,row,col-row); *nm = r;
1070: tmat = (Z **)crmat->body;
1071: for ( i = 0; i < row; i++ )
1072: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
1073: if ( !colstat[j] ) tmi[k++] = utoz(wmi[j]);
1074: initialized = 1;
1075: } else {
1076: for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
1077: if ( j < col ) continue;
1078:
1079: inv = invm(remqi((Q)m1,md),md);
1080: m2 = utoz(md); mulz(m1,m2,&m3);
1081: for ( i = 0; i < row; i++ )
1082: switch ( stat[i] ) {
1083: case 1:
1084: /* consistency check */
1085: ri = (Z *)BDY(r)[i]; wmi = wmat[i];
1086: for ( j = 0; j < col; j++ ) if ( colstat[j] ) break;
1087: h = md-remqi((Q)dn[i],md);
1088: for ( j++, k = 0; j < col; j++ )
1089: if ( !colstat[j] ) {
1090: t = remqi((Q)ri[k],md);
1091: DMAR(wmi[i],h,t,md,t1);
1092: if ( t1 ) break;
1093: }
1094: if ( j == col ) { stat[i]++; break; }
1095: else {
1096: /* fall to the case 0 */
1097: stat[i] = 0;
1098: }
1099: case 0:
1100: tmi = tmat[i]; wmi = wmat[i];
1101: for ( j = k = 0; j < col; j++ )
1102: if ( !colstat[j] ) {
1103: if ( tmi[k] ) {
1104: /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
1105: t = remqi((Q)tmi[k],md);
1106: if ( wmi[j] >= t ) t = wmi[j]-t;
1107: else t = md-(t-wmi[j]);
1108: DMAR(t,inv,0,md,t1)
1109: u = utoz(t1); mulz(m1,u,&s);
1110: addz(tmi[k],s,&u); tmi[k] = u;
1111: } else if ( wmi[j] ) {
1112: /* f3 = m1*(m1 mod m2)^(-1)*f2 */
1113: DMAR(wmi[j],inv,0,md,t)
1114: u = utoz(t); mulz(m1,u,&s); tmi[k] = s;
1115: }
1116: k++;
1117: }
1118: break;
1119: case 2: default:
1120: break;
1121: }
1122: m1 = m3;
1123: if ( ind % 4 )
1124: ret = 0;
1125: else
1126: ret = intmtoratm2(crmat,m1,*nm,dn,stat);
1127: if ( ret ) {
1128: *rindp = rind = (int *)MALLOC_ATOMIC(row*sizeof(int));
1129: *cindp = cind = (int *)MALLOC_ATOMIC((col-row)*sizeof(int));
1130: for ( j = k = l = 0; j < col; j++ )
1131: if ( colstat[j] ) rind[k++] = j;
1132: else cind[l++] = j;
1133: return gensolve_check2(mat,*nm,dn,rind,cind);
1134: }
1135: }
1136: }
1137: }
1138:
1139: int generic_gauss_elim_direct(MAT mat,MAT *nm,Z *dn,int **rindp,int **cindp){
1140: Z **bmat,*s;
1141: Z u,v,w,x,d,t,y;
1142: int row,col,i,j,k,l,m,rank;
1143: int *colstat,*colpos,*cind;
1144: MAT r,in;
1145:
1146: row = mat->row; col = mat->col;
1147: MKMAT(in,row,col);
1148: for ( i = 0; i < row; i++ )
1149: for ( j = 0; j < col; j++ ) in->body[i][j] = mat->body[i][j];
1150: bmat = (Z **)in->body;
1151: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
1152: *rindp = colpos = (int *)MALLOC_ATOMIC(row*sizeof(int));
1153: for ( j = 0, rank = 0, d = ONE; j < col; j++ ) {
1154: for ( i = rank; i < row && !bmat[i][j]; i++ );
1155: if ( i == row ) { colstat[j] = 0; continue; }
1156: else { colstat[j] = 1; colpos[rank] = j; }
1157: if ( i != rank )
1158: for ( k = j; k < col; k++ ) {
1159: t = bmat[i][k]; bmat[i][k] = bmat[rank][k]; bmat[rank][k] = t;
1160: }
1161: for ( i = rank+1, v = bmat[rank][j]; i < row; i++ )
1162: for ( k = j, u = bmat[i][j]; k < col; k++ ) {
1163: mulz(bmat[i][k],v,&w); mulz(bmat[rank][k],u,&x);
1164: subz(w,x,&y); divsz(y,d,&bmat[i][k]);
1165: }
1166: d = v; rank++;
1167: }
1168: *dn = d;
1169: s = (Z *)MALLOC(col*sizeof(Z));
1170: for ( i = rank-1; i >= 0; i-- ) {
1171: for ( k = colpos[i]; k < col; k++ ) mulz(bmat[i][k],d,&s[k]);
1172: for ( m = rank-1; m > i; m-- ) {
1173: for ( k = colpos[m], u = bmat[i][k]; k < col; k++ ) {
1174: mulz(bmat[m][k],u,&w); subz(s[k],w,&x); s[k] = x;
1175: }
1176: }
1177: for ( k = colpos[i], u = bmat[i][k]; k < col; k++ )
1178: divz(s[k],u,&bmat[i][k]);
1179: }
1180: *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
1181: MKMAT(r,rank,col-rank); *nm = r;
1182: for ( j = 0, k = 0; j < col; j++ )
1183: if ( !colstat[j] ) {
1184: cind[k] = j;
1185: for ( i = 0; i < rank; i++ ) r->body[i][k] = bmat[i][j];
1186: k++;
1187: }
1188: return rank;
1189: }
1190:
1.8 noro 1191: int mpz_intmtoratm(mpz_t **mat,int row,int col,mpz_t md,mpz_t **nm,mpz_t dn)
1192: {
1193: mpz_t t,s,b,u,nm1,dn1;
1194: int i,j,k,l,ret;
1195: mpz_t *mi,*nmk;
1196:
1197: if ( UNIMPZ(md) )
1198: return 0;
1199: mpz_init(t); mpz_init(s); mpz_init(b); mpz_init(u);
1200: mpz_init(nm1); mpz_init(dn1);
1201: mpz_fdiv_q_2exp(t,md,1); mpz_sqrt(s,t); mpz_fdiv_q_2exp(b,s,64);
1202: if ( !mpz_sgn(b) ) mpz_set_ui(b,1);
1203: mpz_set_ui(dn,1);
1204: for ( i = 0; i < row; i++ )
1205: for ( j = 0, mi = mat[i]; j < col; j++ )
1206: if ( mpz_sgn(mi[j]) ) {
1207: mpz_mul(s,mi[j],dn);
1208: mpz_mod(u,s,md);
1209: ret = mpz_inttorat(u,md,b,nm1,dn1);
1210: if ( !ret )
1211: return 0;
1212: else {
1213: if ( !UNIMPZ(dn1) ) {
1214: for ( k = 0; k < i; k++ )
1215: for ( l = 0, nmk = nm[k]; l < col; l++ ) mpz_mul(nmk[l],nmk[l],dn1);
1216: for ( l = 0, nmk = nm[i]; l < j; l++ ) mpz_mul(nmk[l],nmk[l],dn1);
1217: }
1218: mpz_set(nm[i][j],nm1);
1219: mpz_mul(dn,dn,dn1);
1220: }
1221: }
1222: return 1;
1223: }
1224:
1.1 noro 1225: int intmtoratm(MAT mat,Z md,MAT nm,Z *dn)
1226: {
1227: Z t,s,b,dn0,dn1,nm1,q,u,unm,udn,dmy;
1228: int i,j,k,l,row,col,sgn,ret;
1229: Z **rmat,**tmat,*tmi,*nmk;
1230:
1231: if ( UNIQ(md) )
1232: return 0;
1233: row = mat->row; col = mat->col;
1234: bshiftz(md,1,&t);
1.11 noro 1235: isqrt(t,&s);
1.1 noro 1236: bshiftz(s,64,&b);
1237: if ( !b ) b = ONE;
1238: dn0 = ONE;
1239: tmat = (Z **)mat->body;
1240: rmat = (Z **)nm->body;
1241: for ( i = 0; i < row; i++ )
1242: for ( j = 0, tmi = tmat[i]; j < col; j++ )
1243: if ( tmi[j] ) {
1244: mulz(tmi[j],dn0,&s);
1245: divqrz(s,md,&dmy,&u);
1246: ret = inttorat(u,md,b,&nm1,&dn1);
1247: if ( !ret ) return 0;
1248: else {
1249: if ( !UNIQ(dn1) ) {
1250: for ( k = 0; k < i; k++ )
1251: for ( l = 0, nmk = rmat[k]; l < col; l++ ) {
1252: mulz(nmk[l],dn1,&q); nmk[l] = q;
1253: }
1254: for ( l = 0, nmk = rmat[i]; l < j; l++ ) {
1255: mulz(nmk[l],dn1,&q); nmk[l] = q;
1256: }
1257: }
1258: rmat[i][j] = nm1;
1259: mulz(dn0,dn1,&q); dn0 = q;
1260: }
1261: }
1262: *dn = dn0;
1263: return 1;
1264: }
1265:
1266: int intmtoratm2(MAT mat,Z md,MAT nm,Z *dn,int *stat)
1267: {
1268: int row,col,i,j,ret;
1269: Z dn0,dn1,t,s,b;
1270: Z *w,*tmi;
1271: Z **tmat;
1272:
1273: bshiftz(md,1,&t);
1274: isqrtz(t,&s);
1275: bshiftz(s,64,&b);
1276: tmat = (Z **)mat->body;
1277: if ( UNIQ(md) ) return 0;
1278: row = mat->row; col = mat->col;
1279: dn0 = ONE;
1280: for ( i = 0; i < row; i++ )
1281: if ( cmpz(dn[i],dn0) > 0 ) dn0 = dn[i];
1282: w = (Z *)MALLOC(col*sizeof(Z));
1283: for ( i = 0; i < row; i++ )
1284: if ( stat[i] == 0 ) {
1285: for ( j = 0, tmi = tmat[i]; j < col; j++ )
1286: mulz(tmi[j],dn0,&w[j]);
1287: ret = intvtoratv(w,col,md,b,(Z *)BDY(nm)[i],&dn[i]);
1288: if ( ret ) {
1289: stat[i] = 1;
1290: mulz(dn0,dn[i],&t); dn[i] = t; dn0 = t;
1291: }
1292: }
1293: for ( i = 0; i < row; i++ ) if ( !stat[i] ) break;
1294: if ( i == row ) return 1;
1295: else return 0;
1296: }
1297:
1298: int intvtoratv(Z *v,int n,Z md,Z b,Z *nm,Z *dn)
1299: {
1300: Z dn0,dn1,q,s,u,nm1,unm,udn,dmy;
1301: Z *nmk;
1302: int j,l,col,ret,sgn;
1303:
1304: for ( j = 0; j < n; j++ ) nm[j] = 0;
1305: dn0 = ONE;
1306: for ( j = 0; j < n; j++ ) {
1307: if ( !v[j] ) continue;
1308: mulz(v[j],dn0,&s);
1309: divqrz(s,md,&dmy,&u);
1310: ret = inttorat(u,md,b,&nm1,&dn1);
1311: if ( !ret ) return 0;
1312: if ( !UNIQ(dn1) )
1313: for ( l = 0; l < j; l++ ) {
1314: mulz(nm[l],dn1,&q); nm[l] = q;
1315: }
1316: nm[j] = nm1;
1317: mulz(dn0,dn1,&q); dn0 = q;
1318: }
1319: *dn = dn0;
1320: return 1;
1321: }
1322:
1323: /* assuming 0 < c < m */
1324:
1.8 noro 1325: int mpz_inttorat(mpz_t c,mpz_t m,mpz_t b,mpz_t nm,mpz_t dn)
1326: {
1327: mpz_t u1,v1,u2,v2,r1,r2;
1328: mpz_t q,t;
1329:
1330: mpz_init_set_ui(u1,0); mpz_init_set_ui(v1,1);
1331: mpz_init_set(u2,m); mpz_init_set(v2,c);
1332: mpz_init(q); mpz_init(t); mpz_init(r1); mpz_init(r2);
1333: while ( mpz_cmp(v2,b) >= 0 ) {
1334: /* r2 = u2-q*v2 */
1335: mpz_fdiv_qr(q,r2,u2,v2);
1336: mpz_set(u2,v2); mpz_set(v2,r2);
1337: /* r1 = u1-q*v1 */
1338: mpz_mul(t,q,v1); mpz_sub(r1,u1,t);
1339: mpz_set(u1,v1); mpz_set(v1,r1);
1340: }
1341: if ( mpz_cmp(v1,b) >= 0 ) return 0;
1342: else {
1.14 noro 1343: mpz_gcd(t,v1,v2);
1344: if ( UNIMPZ(t) )
1345: mpz_set_ui(r1,0);
1346: else {
1347: /* v1 /= t, v2 /= t, t=c*v1-v2, r1=t%m */
1348: mpz_divexact(v1,v1,t); mpz_divexact(v2,v2,t);
1349: mpz_mul(t,c,v1); mpz_sub(t,t,v2); mpz_mod(r1,t,m);
1350: }
1351: if ( mpz_sgn(r1) ) return 0;
1.8 noro 1352: if ( mpz_sgn(v1)<0 ) {
1353: mpz_neg(dn,v1); mpz_neg(nm,v2);
1354: } else {
1355: mpz_set(dn,v1); mpz_set(nm,v2);
1356: }
1357: return 1;
1358: }
1359: }
1360:
1.1 noro 1361: int inttorat(Z c,Z m,Z b,Z *nmp,Z *dnp)
1362: {
1.14 noro 1363: Z qq,t,s,r,u1,v1,r1;
1.1 noro 1364: Z q,u2,v2,r2;
1365:
1366: u1 = 0; v1 = ONE; u2 = m; v2 = c;
1367: while ( cmpz(v2,b) >= 0 ) {
1368: divqrz(u2,v2,&q,&r2); u2 = v2; v2 = r2;
1369: mulz(q,v1,&t); subz(u1,t,&r1); u1 = v1; v1 = r1;
1370: }
1371: if ( cmpz(v1,b) >= 0 ) return 0;
1372: else {
1.14 noro 1373: /* reduction and check */
1374: /* v2/v1 = u2/u1, c*u1-u2 = 0 mod m? */
1375: gcdz(v1,v2,&t);
1376: if ( UNIZ(t) ) {
1377: u1 = v1; u2 = v2; r = 0;
1378: } else {
1379: divsz(v1,t,&u1); divsz(v2,t,&u2);
1380: mulz(c,u1,&t); subz(t,u2,&s); remz(s,m,&r);
1381: }
1382: if ( r ) return 0;
1383: if ( mpz_sgn(BDY(u1))<0 ) {
1384: chsgnz(u1,dnp); chsgnz(u2,nmp);
1.1 noro 1385: } else {
1.14 noro 1386: *dnp = u1; *nmp = u2;
1.1 noro 1387: }
1388: return 1;
1389: }
1390: }
1391:
1392: extern int f4_nocheck;
1393:
1.12 noro 1394: int mpz_gensolve_check(MAT mat,mpz_t **nm,mpz_t dn,int rank,int clen,int *rind,int *cind)
1.8 noro 1395: {
1.12 noro 1396: int row,col,i,j,k,l;
1.8 noro 1397: mpz_t t;
1398: mpz_t *w;
1399: Z *mati;
1400: mpz_t *nmk;
1401:
1402: if ( f4_nocheck ) return 1;
1.12 noro 1403: row = mat->row; col = mat->col;
1.8 noro 1404: w = (mpz_t *)MALLOC(clen*sizeof(mpz_t));
1405: mpz_init(t);
1406: for ( i = 0; i < clen; i++ ) mpz_init(w[i]);
1407: for ( i = 0; i < row; i++ ) {
1408: mati = (Z *)mat->body[i];
1409: for ( l = 0; l < clen; l++ ) mpz_set_ui(w[l],0);
1410: for ( k = 0; k < rank; k++ )
1411: for ( l = 0, nmk = (mpz_t *)nm[k]; l < clen; l++ ) {
1412: /* w[l] += mati[rind[k]]*nmk[k] */
1413: if ( mati[rind[k]] ) mpz_addmul(w[l],BDY(mati[rind[k]]),nmk[l]);
1414: }
1415: for ( j = 0; j < clen; j++ ) {
1416: if ( mati[cind[j]] ) mpz_mul(t,dn,BDY(mati[cind[j]]));
1417: else mpz_set_ui(t,0);
1418: if ( mpz_cmp(w[j],t) ) break;
1419: }
1420: if ( j != clen ) break;
1421: }
1422: if ( i != row ) return 0;
1423: else return 1;
1424: }
1425:
1.1 noro 1426: int gensolve_check(MAT mat,MAT nm,Z dn,int *rind,int *cind)
1427: {
1428: int row,col,rank,clen,i,j,k,l;
1429: Z s,t;
1430: Z *w;
1431: Z *mati,*nmk;
1432:
1433: if ( f4_nocheck ) return 1;
1434: row = mat->row; col = mat->col; rank = nm->row; clen = nm->col;
1435: w = (Z *)MALLOC(clen*sizeof(Z));
1436: for ( i = 0; i < row; i++ ) {
1437: mati = (Z *)mat->body[i];
1438: bzero(w,clen*sizeof(Z));
1439: for ( k = 0; k < rank; k++ )
1440: for ( l = 0, nmk = (Z *)nm->body[k]; l < clen; l++ ) {
1441: mulz(mati[rind[k]],nmk[l],&t); addz(w[l],t,&s); w[l] = s;
1442: }
1443: for ( j = 0; j < clen; j++ ) {
1444: mulz(dn,mati[cind[j]],&t);
1445: if ( cmpz(w[j],t) ) break;
1446: }
1447: if ( j != clen ) break;
1448: }
1449: if ( i != row ) return 0;
1450: else return 1;
1451: }
1452:
1453: int gensolve_check2(MAT mat,MAT nm,Z *dn,int *rind,int *cind)
1454: {
1455: int row,col,rank,clen,i,j,k,l;
1456: Z s,t,u,d;
1457: Z *w,*m;
1458: Z *mati,*nmk;
1459:
1460: if ( f4_nocheck ) return 1;
1461: row = mat->row; col = mat->col; rank = nm->row; clen = nm->col;
1462: w = (Z *)MALLOC(clen*sizeof(Z));
1463: m = (Z *)MALLOC(clen*sizeof(Z));
1464: for ( d = dn[0], i = 1; i < rank; i++ ) {
1465: lcmz(d,dn[i],&t); d = t;
1466: }
1467: for ( i = 0; i < rank; i++ ) divsz(d,dn[i],&m[i]);
1468: for ( i = 0; i < row; i++ ) {
1469: mati = (Z *)mat->body[i];
1470: bzero(w,clen*sizeof(Z));
1471: for ( k = 0; k < rank; k++ ) {
1472: mulz(mati[rind[k]],m[k],&u);
1473: for ( l = 0, nmk = (Z *)nm->body[k]; l < clen; l++ ) {
1474: mulz(u,nmk[l],&t); addz(w[l],t,&s); w[l] = s;
1475: }
1476: }
1477: for ( j = 0; j < clen; j++ ) {
1478: mulz(d,mati[cind[j]],&t);
1479: if ( cmpz(w[j],t) ) break;
1480: }
1481: if ( j != clen ) break;
1482: }
1483: if ( i != row ) return 0;
1484: else return 1;
1485: }
1486:
1487: void isqrtz(Z a,Z *r)
1488: {
1489: int k;
1490: Z x,t,x2,xh,quo,rem;
1491: Z two;
1492:
1493: if ( !a ) *r = 0;
1.11 noro 1494: else if ( UNIZ(a) ) *r = ONE;
1.1 noro 1495: else {
1496: k = z_bits((Q)a); /* a <= 2^k-1 */
1497: bshiftz(ONE,-((k>>1)+(k&1)),&x); /* a <= x^2 */
1.5 noro 1498: STOZ(2,two);
1.1 noro 1499: while ( 1 ) {
1500: pwrz(x,two,&t);
1501: if ( cmpz(t,a) <= 0 ) {
1502: *r = x; return;
1503: } else {
1504: if ( mpz_tstbit(BDY(x),0) ) addz(x,a,&t);
1505: else t = a;
1506: bshiftz(x,-1,&x2); divqrz(t,x2,&quo,&rem);
1507: bshiftz(x,1,&xh); addz(quo,xh,&x);
1508: }
1509: }
1510: }
1511: }
1512:
1513: void bshiftz(Z a,int n,Z *r)
1514: {
1515: mpz_t t;
1516:
1517: if ( !a ) *r = 0;
1518: else if ( n == 0 ) *r = a;
1519: else if ( n < 0 ) {
1520: mpz_init(t); mpz_mul_2exp(t,BDY(a),-n); MPZTOZ(t,*r);
1521: } else {
1522: mpz_init(t); mpz_fdiv_q_2exp(t,BDY(a),n);
1523: if ( !mpz_sgn(t) ) *r = 0;
1524: else MPZTOZ(t,*r);
1525: }
1526: }
1527:
1528: void addlf(Z a,Z b,Z *c)
1529: {
1530: addz(a,b,c);
1531: if ( !lf_lazy ) {
1532: if ( cmpz(*c,current_mod_lf) >= 0 ) {
1533: subz(*c,current_mod_lf,c);
1534: }
1535: }
1536: }
1537:
1538: void sublf(Z a,Z b,Z *c)
1539: {
1540: subz(a,b,c);
1541: if ( !lf_lazy ) {
1542: remz(*c,current_mod_lf,c);
1543: }
1544: }
1545:
1546: void mullf(Z a,Z b,Z *c)
1547: {
1548: mulz(a,b,c);
1549: if ( !lf_lazy ) {
1550: remz(*c,current_mod_lf,c);
1551: }
1552: }
1553:
1554: void divlf(Z a,Z b,Z *c)
1555: {
1556: Z inv;
1557:
1558: invz(b,current_mod_lf,&inv);
1559: mulz(a,inv,c);
1560: if ( !lf_lazy ) {
1561: remz(*c,current_mod_lf,c);
1562: }
1563: }
1564:
1565: void chsgnlf(Z a,Z *c)
1566: {
1567: chsgnz(a,c);
1568: if ( !lf_lazy ) {
1569: remz(*c,current_mod_lf,c);
1570: }
1571: }
1572:
1573: void lmtolf(LM a,Z *b)
1574: {
1575: if ( !a ) *b = 0;
1576: else {
1577: MPZTOZ(BDY(a),*b);
1578: }
1579: }
1580:
1581: void setmod_lf(Z p)
1582: {
1583: current_mod_lf = p;
1584: current_mod_lf_size = mpz_size(BDY(current_mod_lf))+1;
1585: }
1586:
1587: void simplf_force(Z a,Z *b)
1588: {
1589: remz(a,current_mod_lf,b);
1590: }
1591:
1592: int generic_gauss_elim_hensel(MAT mat,MAT *nmmat,Z *dn,int **rindp,int **cindp)
1593: {
1594: MAT bmat,xmat;
1595: Z **a0,**a,**b,**x,**nm;
1596: Z *ai,*bi,*xi;
1597: int row,col;
1598: int **w;
1599: int *wi;
1600: int **wc;
1601: Z mdq,q,s,u;
1602: Z tn;
1603: int ind,md,i,j,k,l,li,ri,rank;
1604: unsigned int t;
1605: int *cinfo,*rinfo;
1606: int *rind,*cind;
1607: int count;
1608: int ret;
1.3 noro 1609: struct oEGT eg_mul1,eg_mul2,tmp0,tmp1,tmp2;
1.1 noro 1610: int period;
1611: int *wx,*ptr;
1612: int wxsize,nsize;
1613: Z wn;
1614: Z wq;
1615:
1.9 noro 1616: #if SIZEOF_LONG == 8
1.11 noro 1617: return generic_gauss_elim_hensel64(mat,nmmat,dn,rindp,cindp,0);
1.9 noro 1618: #endif
1.3 noro 1619: init_eg(&eg_mul1); init_eg(&eg_mul2);
1.1 noro 1620: a0 = (Z **)mat->body;
1621: row = mat->row; col = mat->col;
1622: w = (int **)almat(row,col);
1623: for ( ind = 0; ; ind++ ) {
1624: md = get_lprime(ind);
1.5 noro 1625: STOZ(md,mdq);
1.1 noro 1626: for ( i = 0; i < row; i++ )
1627: for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
1628: wi[j] = remqi((Q)ai[j],md);
1629:
1630: if ( DP_Print > 3 ) {
1631: fprintf(asir_out,"LU decomposition.."); fflush(asir_out);
1632: }
1633: rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo);
1634: if ( DP_Print > 3 ) {
1635: fprintf(asir_out,"done.\n"); fflush(asir_out);
1636: }
1637: a = (Z **)almat_pointer(rank,rank); /* lhs mat */
1638: MKMAT(bmat,rank,col-rank); b = (Z **)bmat->body; /* lhs mat */
1639: for ( j = li = ri = 0; j < col; j++ )
1640: if ( cinfo[j] ) {
1641: /* the column is in lhs */
1642: for ( i = 0; i < rank; i++ ) {
1643: w[i][li] = w[i][j];
1644: a[i][li] = a0[rinfo[i]][j];
1645: }
1646: li++;
1647: } else {
1648: /* the column is in rhs */
1649: for ( i = 0; i < rank; i++ )
1650: b[i][ri] = a0[rinfo[i]][j];
1651: ri++;
1652: }
1653:
1654: /* solve Ax=B; A: rank x rank, B: rank x ri */
1655: /* algorithm
1656: c <- B
1657: x <- 0
1658: q <- 1
1659: do
1660: t <- A^(-1)c mod p
1661: x <- x+qt
1662: c <- (c-At)/p
1663: q <- qp
1664: end do
1665: then Ax-B=0 mod q and b=(B-Ax)/q hold after "do".
1666: */
1667: MKMAT(xmat,rank,ri); x = (Z **)(xmat)->body;
1668: MKMAT(*nmmat,rank,ri); nm = (Z **)(*nmmat)->body;
1669: wc = (int **)almat(rank,ri);
1670: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
1671: *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));
1672:
1673: period = F4_INTRAT_PERIOD;
1674: for ( q = ONE, count = 0; ; ) {
1.3 noro 1675: /* check Ax=B mod q */
1.1 noro 1676: if ( DP_Print > 3 )
1677: fprintf(stderr,"o");
1678: /* wc = b mod md */
1679: for ( i = 0; i < rank; i++ )
1.3 noro 1680: for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )
1.1 noro 1681: wi[j] = remqi((Q)bi[j],md);
1.3 noro 1682: /* wc = A^(-1)wc; wc is not normalized */
1683: solve_by_lu_mod(w,rank,md,wc,ri,0);
1.1 noro 1684: /* x += q*wc */
1.3 noro 1685: get_eg(&tmp0);
1.1 noro 1686: for ( i = 0; i < rank; i++ )
1687: for ( j = 0, wi = wc[i]; j < ri; j++ ) mul1addtoz(q,wi[j],&x[i][j]);
1.3 noro 1688: /* b =(b-A*wc)/md */
1689: get_eg(&tmp1); add_eg(&eg_mul1,&tmp0,&tmp1);
1.1 noro 1690: for ( i = 0; i < rank; i++ )
1691: for ( j = 0; j < ri; j++ ) {
1.3 noro 1692: mpz_t uz;
1693:
1694: if ( b[i][j] )
1695: mpz_init_set(uz,BDY(b[i][j]));
1696: else
1697: mpz_init_set_ui(uz,0);
1698: for ( k = 0; k < rank; k++ ) {
1699: if ( a[i][k] && wc[k][j] ) {
1700: if ( wc[k][j] < 0 )
1701: mpz_addmul_ui(uz,BDY(a[i][k]),-wc[k][j]);
1702: else
1703: mpz_submul_ui(uz,BDY(a[i][k]),wc[k][j]);
1704: }
1705: }
1706: MPZTOZ(uz,u);
1.1 noro 1707: divsz(u,mdq,&b[i][j]);
1708: }
1.3 noro 1709: get_eg(&tmp2); add_eg(&eg_mul2,&tmp1,&tmp2);
1.1 noro 1710: count++;
1711: /* q = q*md */
1712: mulz(q,mdq,&u); q = u;
1713: if ( count == period ) {
1714: ret = intmtoratm(xmat,q,*nmmat,dn);
1715: if ( ret ) {
1.3 noro 1716: print_eg("MUL1",&eg_mul1);
1717: print_eg("MUL2",&eg_mul2);
1.1 noro 1718: for ( j = k = l = 0; j < col; j++ )
1719: if ( cinfo[j] )
1720: rind[k++] = j;
1721: else
1722: cind[l++] = j;
1723: ret = gensolve_check(mat,*nmmat,*dn,rind,cind);
1724: if ( ret ) {
1725: *rindp = rind;
1726: *cindp = cind;
1727: for ( j = k = 0; j < col; j++ )
1728: if ( !cinfo[j] )
1729: cind[k++] = j;
1730: return rank;
1.11 noro 1731: } else
1732: goto reset;
1.1 noro 1733: } else {
1.11 noro 1734: reset:
1.1 noro 1735: period = period*3/2;
1736: count = 0;
1737: }
1738: }
1739: }
1740: }
1741: }
1742:
1743: /* for inv_or_split_dalg */
1744:
1745: int generic_gauss_elim_hensel_dalg(MAT mat,DP *mb,MAT *nmmat,Z *dn,int **rindp,int **cindp)
1746: {
1747: MAT bmat,xmat;
1748: Z **a0,**a,**b,**x,**nm;
1749: Z *ai,*bi,*xi;
1750: int row,col;
1751: int **w;
1752: int *wi;
1753: int **wc;
1754: Z mdq,q,s,u;
1755: Z tn;
1756: int ind,md,i,j,k,l,li,ri,rank;
1757: unsigned int t;
1758: int *cinfo,*rinfo;
1759: int *rind,*cind;
1760: int count;
1761: int ret;
1762: struct oEGT eg_mul,eg_inv,eg_intrat,eg_check,tmp0,tmp1;
1763: int period;
1764: int *wx,*ptr;
1765: int wxsize,nsize;
1766: Z wn;
1767: Z wq;
1768: DP m;
1769:
1.11 noro 1770: #if SIZEOF_LONG == 8
1771: return generic_gauss_elim_hensel64(mat,nmmat,dn,rindp,cindp,mb);
1772: #endif
1.1 noro 1773: a0 = (Z **)mat->body;
1774: row = mat->row; col = mat->col;
1775: w = (int **)almat(row,col);
1776: for ( ind = 0; ; ind++ ) {
1777: md = get_lprime(ind);
1.5 noro 1778: STOZ(md,mdq);
1.1 noro 1779: for ( i = 0; i < row; i++ )
1780: for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
1781: wi[j] = remqi((Q)ai[j],md);
1782:
1783: if ( DP_Print > 3 ) {
1784: fprintf(asir_out,"LU decomposition.."); fflush(asir_out);
1785: }
1786: rank = find_lhs_and_lu_mod((unsigned int **)w,row,col,md,&rinfo,&cinfo);
1787: if ( DP_Print > 3 ) {
1788: fprintf(asir_out,"done.\n"); fflush(asir_out);
1789: }
1790:
1791: /* this part is added for inv_or_split_dalg */
1792: for ( i = 0; i < col-1; i++ ) {
1793: if ( !cinfo[i] ) {
1794: m = mb[i];
1795: for ( j = i+1; j < col-1; j++ )
1796: if ( dp_redble(mb[j],m) )
1797: cinfo[j] = -1;
1798: }
1799: }
1800:
1801: a = (Z **)almat_pointer(rank,rank); /* lhs mat */
1802: MKMAT(bmat,rank,col-rank); b = (Z **)bmat->body; /* lhs mat */
1803: for ( j = li = ri = 0; j < col; j++ )
1.4 noro 1804: if ( cinfo[j] > 0 ) {
1.1 noro 1805: /* the column is in lhs */
1806: for ( i = 0; i < rank; i++ ) {
1807: w[i][li] = w[i][j];
1808: a[i][li] = a0[rinfo[i]][j];
1809: }
1810: li++;
1.4 noro 1811: } else if ( !cinfo[j] ) {
1.1 noro 1812: /* the column is in rhs */
1813: for ( i = 0; i < rank; i++ )
1814: b[i][ri] = a0[rinfo[i]][j];
1815: ri++;
1816: }
1817:
1818: /* solve Ax=B; A: rank x rank, B: rank x ri */
1819: /* algorithm
1820: c <- B
1821: x <- 0
1822: q <- 1
1823: do
1824: t <- A^(-1)c mod p
1825: x <- x+qt
1826: c <- (c-At)/p
1827: q <- qp
1828: end do
1829: then Ax-B=0 mod q and b=(B-Ax)/q hold after "do".
1830: */
1831: MKMAT(xmat,rank,ri); x = (Z **)(xmat)->body;
1832: MKMAT(*nmmat,rank,ri); nm = (Z **)(*nmmat)->body;
1833: wc = (int **)almat(rank,ri);
1834: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
1835: *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));
1836:
1837: period = F4_INTRAT_PERIOD;
1838: for ( q = ONE, count = 0; ; ) {
1839: if ( DP_Print > 3 )
1840: fprintf(stderr,"o");
1841: /* wc = b mod md */
1842: for ( i = 0; i < rank; i++ )
1.3 noro 1843: for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )
1.1 noro 1844: wi[j] = remqi((Q)bi[j],md);
1.11 noro 1845: /* wc = A^(-1)wc; wc is not normalized */
1846: solve_by_lu_mod(w,rank,md,wc,ri,0);
1.1 noro 1847: /* x += q*wc */
1848: for ( i = 0; i < rank; i++ )
1849: for ( j = 0, wi = wc[i]; j < ri; j++ ) mul1addtoz(q,wi[j],&x[i][j]);
1.3 noro 1850: /* b =(b-A*wc)/md */
1.1 noro 1851: for ( i = 0; i < rank; i++ )
1852: for ( j = 0; j < ri; j++ ) {
1.3 noro 1853: mpz_t uz;
1854:
1855: if ( b[i][j] )
1856: mpz_init_set(uz,BDY(b[i][j]));
1857: else
1858: mpz_init_set_ui(uz,0);
1.11 noro 1859: for ( k = 0; k < rank; k++ )
1860: if ( a[i][k] && wc[k][j] )
1.3 noro 1861: mpz_submul_ui(uz,BDY(a[i][k]),wc[k][j]);
1862: MPZTOZ(uz,u);
1.1 noro 1863: divsz(u,mdq,&b[i][j]);
1864: }
1.11 noro 1865:
1.1 noro 1866: count++;
1867: /* q = q*md */
1868: mulz(q,mdq,&u); q = u;
1869: if ( count == period ) {
1870: ret = intmtoratm(xmat,q,*nmmat,dn);
1871: if ( ret ) {
1872: for ( j = k = l = 0; j < col; j++ )
1873: if ( cinfo[j] > 0 )
1874: rind[k++] = j;
1875: else if ( !cinfo[j] )
1876: cind[l++] = j;
1877: ret = gensolve_check(mat,*nmmat,*dn,rind,cind);
1878: if ( ret ) {
1879: *rindp = rind;
1880: *cindp = cind;
1881: for ( j = k = 0; j < col; j++ )
1882: if ( !cinfo[j] )
1883: cind[k++] = j;
1884: return rank;
1.11 noro 1885: } else
1886: goto reset;
1.1 noro 1887: } else {
1.11 noro 1888: reset:
1.1 noro 1889: period = period*3/2;
1890: count = 0;
1891: }
1892: }
1893: }
1894: }
1895: }
1.6 noro 1896:
1897: #if SIZEOF_LONG == 8
1898: mp_limb_t remqi64(Q a,mp_limb_t mod)
1899: {
1900: mp_limb_t c,nm,dn;
1901: mpz_t r;
1902:
1903: if ( !a ) return 0;
1904: else if ( a->z ) {
1905: mpz_init(r);
1906: c = mpz_fdiv_r_ui(r,BDY((Z)a),mod);
1907: } else {
1908: mpz_init(r);
1909: nm = mpz_fdiv_r_ui(r,mpq_numref(BDY(a)),mod);
1910: dn = mpz_fdiv_r_ui(r,mpq_denref(BDY(a)),mod);
1911: dn = invmod64(dn,mod);
1912: c = mulmod64(nm,dn,mod);
1913: }
1914: return c;
1915: }
1916:
1917: int generic_gauss_elim_mod64(mp_limb_t **mat,int row,int col,mp_limb_t md,int *colstat);
1918: mp_limb_t get_lprime64(int ind);
1919:
1.8 noro 1920: void mpz_print(mpz_t a)
1921: {
1922: mpz_out_str(stdout,10,a); printf("\n");
1923: }
1924:
1925: void mpz_printmat(mpz_t **a,int row,int col)
1926: {
1927: int i,j;
1928: for ( i = 0; i < row; i++ ) {
1929: for ( j = 0; j < col; j++ ) {
1930: mpz_out_str(stdout,10,a[i][j]); printf(" ");
1931: }
1932: printf("\n");
1933: }
1934: }
1935:
1936: mpz_t **mpz_allocmat(int row,int col)
1937: {
1938: mpz_t **p;
1939: int i,j;
1940:
1941: p = (mpz_t **)MALLOC(row*sizeof(mpz_t *));
1942: for ( i = 0; i < row; i++ ) {
1943: p[i] = (mpz_t *)MALLOC(col*sizeof(mpz_t));
1944: for ( j = 0; j < col; j++ ) mpz_init(p[i][j]);
1945: }
1946: return p;
1947: }
1948:
1949: #if 1
1950: int generic_gauss_elim64(MAT mat,MAT *nm,Z *dn,int **rindp,int **cindp)
1951: {
1952: mp_limb_t **wmat;
1953: mp_limb_t *wmi;
1954: mp_limb_t md,inv,t,t1;
1955: Z z;
1956: Z **bmat,*bmi;
1957: mpz_t **tmat,**num;
1958: mpz_t *tmi;
1959: mpz_t den;
1960: mpz_t q,m1,m3,s,u;
1961: int *colstat,*wcolstat,*rind,*cind;
1962: int row,col,ind,i,j,k,l,rank,rank0;
1963: MAT r;
1964: int ret;
1965:
1966: bmat = (Z **)mat->body;
1967: row = mat->row; col = mat->col;
1968: wmat = (mp_limb_t **)almat64(row,col);
1969: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
1970: wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
1971: mpz_init(m1); mpz_init(m3); mpz_init(den);
1972: for ( ind = 0; ; ind++ ) {
1973: if ( DP_Print ) {
1974: fprintf(asir_out,"."); fflush(asir_out);
1975: }
1976: md = get_lprime64(ind);
1977: for ( i = 0; i < row; i++ )
1978: for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
1979: wmi[j] = bmi[j]==0?0:mpz_fdiv_ui(BDY(bmi[j]),md);
1980: rank = generic_gauss_elim_mod64(wmat,row,col,md,wcolstat);
1981: if ( !ind ) {
1982: RESET:
1983: mpz_set_ui(m1,md);
1984: rank0 = rank;
1985: bcopy(wcolstat,colstat,col*sizeof(int));
1986: // crmat
1987: tmat = mpz_allocmat(rank,col-rank);
1988: //
1989: num = mpz_allocmat(rank,col-rank);
1990: for ( i = 0; i < rank; i++ )
1991: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
1992: if ( !colstat[j] ) { mpz_set_ui(tmi[k],wmi[j]); k++; }
1993: } else {
1994: if ( rank < rank0 ) {
1995: if ( DP_Print ) {
1996: fprintf(asir_out,"lower rank matrix; continuing...\n");
1997: fflush(asir_out);
1998: }
1999: continue;
2000: } else if ( rank > rank0 ) {
2001: if ( DP_Print ) {
2002: fprintf(asir_out,"higher rank matrix; resetting...\n");
2003: fflush(asir_out);
2004: }
2005: goto RESET;
2006: } else {
2007: for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
2008: if ( j < col ) {
2009: if ( DP_Print ) {
2010: fprintf(asir_out,"inconsitent colstat; resetting...\n");
2011: fflush(asir_out);
2012: }
2013: goto RESET;
2014: }
2015: }
2016:
2017: inv = invmod64(mpz_fdiv_ui(m1,md),md);
2018: mpz_mul_ui(m3,m1,md);
2019: for ( i = 0; i < rank; i++ )
2020: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
2021: if ( !colstat[j] ) {
2022: if ( mpz_sgn(tmi[k]) ) {
2023: /* f3 = f1+m1*(m1 mod md)^(-1)*(f2 - f1 mod md) */
2024: t = mpz_fdiv_ui(tmi[k],md);
2025: if ( wmi[j] >= t ) t = wmi[j]-t;
2026: else t = md-(t-wmi[j]);
2027: mpz_addmul_ui(tmi[k],m1,mulmod64(t,inv,md));
2028: } else if ( wmi[j] ) {
2029: /* f3 = m1*(m1 mod m2)^(-1)*f2 */
2030: mpz_mul_ui(tmi[k],m1,mulmod64(wmi[j],inv,md));
2031: }
2032: k++;
2033: }
2034: mpz_set(m1,m3);
2035: if ( ind % F4_INTRAT_PERIOD )
2036: ret = 0;
2037: else
2038: ret = mpz_intmtoratm(tmat,rank,col-rank,m1,num,den);
2039: if ( ret ) {
2040: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
2041: *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
2042: for ( j = k = l = 0; j < col; j++ )
2043: if ( colstat[j] ) rind[k++] = j;
2044: else cind[l++] = j;
1.12 noro 2045: if ( mpz_gensolve_check(mat,num,den,rank,col-rank,rind,cind) ) {
1.8 noro 2046: MKMAT(r,rank,col-rank); *nm = r;
2047: for ( i = 0; i < rank; i++ )
2048: for ( j = 0; j < col-rank; j++ ) {
2049: MPZTOZ(num[i][j],z); BDY(r)[i][j] = z;
2050: }
2051: MPZTOZ(den,*dn);
2052: return rank;
2053: }
2054: }
2055: }
2056: }
2057: }
2058: #else
1.6 noro 2059: int generic_gauss_elim64(MAT mat,MAT *nm,Z *dn,int **rindp,int **cindp)
2060: {
2061: mp_limb_t **wmat;
2062: mp_limb_t *wmi;
2063: mp_limb_t md,inv,t,t1;
2064: Z **bmat,**tmat,*bmi,*tmi;
2065: Z q,m1,m2,m3,s,u;
2066: int *colstat,*wcolstat,*rind,*cind;
2067: int row,col,ind,i,j,k,l,rank,rank0;
2068: MAT r,crmat;
2069: int ret;
2070:
2071: bmat = (Z **)mat->body;
2072: row = mat->row; col = mat->col;
2073: wmat = (mp_limb_t **)almat64(row,col);
2074: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
2075: wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
2076: for ( ind = 0; ; ind++ ) {
2077: if ( DP_Print ) {
2078: fprintf(asir_out,"."); fflush(asir_out);
2079: }
2080: md = get_lprime64(ind);
2081: for ( i = 0; i < row; i++ )
2082: for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
2083: wmi[j] = remqi64((Q)bmi[j],md);
2084: rank = generic_gauss_elim_mod64(wmat,row,col,md,wcolstat);
2085: if ( !ind ) {
2086: RESET:
2087: UTOZ(md,m1);
2088: rank0 = rank;
2089: bcopy(wcolstat,colstat,col*sizeof(int));
2090: MKMAT(crmat,rank,col-rank);
2091: MKMAT(r,rank,col-rank); *nm = r;
2092: tmat = (Z **)crmat->body;
2093: for ( i = 0; i < rank; i++ )
2094: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
2095: if ( !colstat[j] ) { UTOZ(wmi[j],tmi[k]); k++; }
2096: } else {
2097: if ( rank < rank0 ) {
2098: if ( DP_Print ) {
2099: fprintf(asir_out,"lower rank matrix; continuing...\n");
2100: fflush(asir_out);
2101: }
2102: continue;
2103: } else if ( rank > rank0 ) {
2104: if ( DP_Print ) {
2105: fprintf(asir_out,"higher rank matrix; resetting...\n");
2106: fflush(asir_out);
2107: }
2108: goto RESET;
2109: } else {
2110: for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
2111: if ( j < col ) {
2112: if ( DP_Print ) {
2113: fprintf(asir_out,"inconsitent colstat; resetting...\n");
2114: fflush(asir_out);
2115: }
2116: goto RESET;
2117: }
2118: }
2119:
2120: inv = invmod64(remqi64((Q)m1,md),md);
2121: UTOZ(md,m2); mulz(m1,m2,&m3);
2122: for ( i = 0; i < rank; i++ )
2123: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
2124: if ( !colstat[j] ) {
2125: if ( tmi[k] ) {
2126: /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
2127: t = remqi64((Q)tmi[k],md);
2128: if ( wmi[j] >= t ) t = wmi[j]-t;
2129: else t = md-(t-wmi[j]);
2130: t1 = mulmod64(t,inv,md);
2131: UTOZ(t1,u); mulz(m1,u,&s);
2132: addz(tmi[k],s,&u); tmi[k] = u;
2133: } else if ( wmi[j] ) {
2134: /* f3 = m1*(m1 mod m2)^(-1)*f2 */
2135: t = mulmod64(wmi[j],inv,md);
2136: UTOZ(t,u); mulz(m1,u,&s); tmi[k] = s;
2137: }
2138: k++;
2139: }
2140: m1 = m3;
2141: if ( ind % F4_INTRAT_PERIOD )
2142: ret = 0;
2143: else
2144: ret = intmtoratm(crmat,m1,*nm,dn);
2145: if ( ret ) {
2146: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
2147: *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
2148: for ( j = k = l = 0; j < col; j++ )
2149: if ( colstat[j] ) rind[k++] = j;
2150: else cind[l++] = j;
2151: if ( gensolve_check(mat,*nm,*dn,rind,cind) )
2152: return rank;
2153: }
2154: }
2155: }
2156: }
2157: #endif
1.8 noro 2158:
1.11 noro 2159: int generic_gauss_elim_hensel64(MAT mat,MAT *nmmat,Z *dn,int **rindp,int **cindp,DP *mb)
1.9 noro 2160: {
2161: MAT r;
2162: Z z;
2163: Z **a0;
2164: Z *ai;
2165: mpz_t **a,**b,**x,**nm;
2166: mpz_t *bi,*xi;
2167: mpz_t q,u,den;
2168: mp_limb_t **w;
2169: mp_limb_t *wi;
2170: mp_limb_t **wc;
2171: mp_limb_t md;
2172: int row,col;
2173: int ind,i,j,k,l,li,ri,rank;
2174: int *cinfo,*rinfo;
2175: int *rind,*cind;
2176: int count;
2177: int ret;
2178: int period;
1.11 noro 2179: DP m;
1.9 noro 2180:
2181: a0 = (Z **)mat->body;
2182: row = mat->row; col = mat->col;
2183: w = (mp_limb_t **)almat64(row,col);
2184: for ( ind = 0; ; ind++ ) {
2185: md = get_lprime64(ind);
2186: for ( i = 0; i < row; i++ )
2187: for ( j = 0, ai = a0[i], wi = w[i]; j < col; j++ )
2188: wi[j] = remqi64((Q)ai[j],md);
2189:
2190: if ( DP_Print > 3 ) {
2191: fprintf(asir_out,"LU decomposition.."); fflush(asir_out);
2192: }
2193: rank = find_lhs_and_lu_mod64(w,row,col,md,&rinfo,&cinfo);
2194: if ( DP_Print > 3 ) {
2195: fprintf(asir_out,"done.\n"); fflush(asir_out);
2196: }
1.11 noro 2197:
2198: if ( mb ) {
2199: /* this part is added for inv_or_split_dalg */
2200: for ( i = 0; i < col-1; i++ ) {
2201: if ( !cinfo[i] ) {
2202: m = mb[i];
2203: for ( j = i+1; j < col-1; j++ )
2204: if ( dp_redble(mb[j],m) )
2205: cinfo[j] = -1;
2206: }
2207: }
2208: }
2209:
1.9 noro 2210: a = (mpz_t **)mpz_allocmat(rank,rank); /* lhs mat */
2211: b = (mpz_t **)mpz_allocmat(rank,col-rank);
2212: for ( j = li = ri = 0; j < col; j++ )
1.11 noro 2213: if ( cinfo[j] > 0 ) {
1.9 noro 2214: /* the column is in lhs */
2215: for ( i = 0; i < rank; i++ ) {
2216: w[i][li] = w[i][j];
2217: if ( a0[rinfo[i]][j] )
2218: mpz_set(a[i][li],BDY(a0[rinfo[i]][j]));
2219: else
2220: mpz_set_ui(a[i][li],0);
2221: }
2222: li++;
1.11 noro 2223: } else if ( !cinfo[j] ) {
1.9 noro 2224: /* the column is in rhs */
2225: for ( i = 0; i < rank; i++ ) {
2226: if ( a0[rinfo[i]][j] )
2227: mpz_set(b[i][ri],BDY(a0[rinfo[i]][j]));
2228: else
2229: mpz_set_ui(b[i][ri],0);
2230: }
2231: ri++;
2232: }
2233:
2234: /* solve Ax=B; A: rank x rank, B: rank x ri */
2235: /* algorithm
2236: c <- B
2237: x <- 0
2238: q <- 1
2239: do
2240: t <- A^(-1)c mod p
2241: x <- x+qt
2242: c <- (c-At)/p
2243: q <- qp
2244: end do
2245: then Ax-B=0 mod q and b=(B-Ax)/q hold after "do".
2246: */
2247: x = (mpz_t **)mpz_allocmat(rank,ri);
2248: nm = (mpz_t **)mpz_allocmat(rank,ri);
2249: wc = (mp_limb_t **)almat64(rank,ri);
2250: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
2251: *cindp = cind = (int *)MALLOC_ATOMIC((ri)*sizeof(int));
2252:
2253: period = F4_INTRAT_PERIOD;
2254: mpz_init_set_ui(q,1);
2255: mpz_init(u);
2256: mpz_init(den);
2257: for ( count = 0; ; ) {
2258: /* check Ax=B mod q */
2259: if ( DP_Print > 3 )
2260: fprintf(stderr,"o");
2261: /* wc = b mod md */
2262: for ( i = 0; i < rank; i++ )
2263: for ( j = 0, bi = b[i], wi = wc[i]; j < ri; j++ )
2264: wi[j] = mpz_fdiv_ui(bi[j],md);
2265: /* wc = A^(-1)wc; wc is not normalized */
2266: solve_by_lu_mod64(w,rank,md,wc,ri,0);
2267: /* x += q*wc */
2268: for ( i = 0; i < rank; i++ )
2269: for ( j = 0, wi = wc[i]; j < ri; j++ )
2270: if ( wi[j] > 0 )
2271: mpz_addmul_ui(x[i][j],q,wi[j]);
2272: else if ( wi[j] < 0 )
2273: mpz_submul_ui(x[i][j],q,-wi[j]);
2274: /* b =(b-A*wc)/md */
2275: for ( i = 0; i < rank; i++ )
2276: for ( j = 0; j < ri; j++ ) {
2277: mpz_set(u,b[i][j]);
2278: for ( k = 0; k < rank; k++ ) {
2279: if ( a[i][k] && wc[k][j] ) {
2280: if ( wc[k][j] < 0 )
2281: mpz_addmul_ui(u,a[i][k],-wc[k][j]);
2282: else
2283: mpz_submul_ui(u,a[i][k],wc[k][j]);
2284: }
2285: }
2286: mpz_divexact_ui(b[i][j],u,md);
2287: }
2288: count++;
2289: /* q = q*md */
2290: mpz_mul_ui(q,q,md);
2291: if ( count == period ) {
2292: ret = mpz_intmtoratm(x,rank,ri,q,nm,den);
2293: if ( ret ) {
2294: for ( j = k = l = 0; j < col; j++ )
1.11 noro 2295: if ( cinfo[j] > 0 )
1.9 noro 2296: rind[k++] = j;
1.11 noro 2297: else if ( !cinfo[j] )
1.9 noro 2298: cind[l++] = j;
1.12 noro 2299: ret = mpz_gensolve_check(mat,nm,den,rank,ri,rind,cind);
1.9 noro 2300: if ( ret ) {
2301: *rindp = rind;
2302: *cindp = cind;
2303: for ( j = k = 0; j < col; j++ )
2304: if ( !cinfo[j] )
2305: cind[k++] = j;
2306: MKMAT(r,rank,ri); *nmmat = r;
2307: for ( i = 0; i < rank; i++ )
2308: for ( j = 0; j < ri; j++ ) {
2309: MPZTOZ(nm[i][j],z); BDY(r)[i][j] = z;
2310: }
2311: MPZTOZ(den,*dn);
2312: return rank;
1.11 noro 2313: } else
2314: goto reset;
1.9 noro 2315: } else {
1.11 noro 2316: reset:
1.9 noro 2317: fprintf(stderr,"F");
2318: period = period*3/2;
2319: count = 0;
2320: }
2321: }
2322: }
2323: }
2324: }
2325:
1.8 noro 2326: #endif
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