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