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