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