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