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