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