![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to gmpq.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine|
Return to gmpq.c CVS log | Up to [local] / OpenXM_contrib2 / asir2000 / engine |
1.1 noro 1: #include "ca.h"
2: #include "gmp.h"
3: #include "base.h"
4: #include "inline.h"
5:
6: mpz_t ONEMPZ;
7: GZ ONEGZ;
1.5 noro 8: int lf_lazy;
9: GZ current_mod_lf;
1.6 noro 10: int current_mod_lf_size;
1.1 noro 11:
12: void isqrtgz(GZ a,GZ *r);
13: void bshiftgz(GZ a,int n,GZ *r);
14:
15: void *gc_realloc(void *p,size_t osize,size_t nsize)
16: {
1.2 noro 17: return (void *)Risa_GC_realloc(p,nsize);
1.1 noro 18: }
19:
20: void gc_free(void *p,size_t size)
21: {
1.2 noro 22: Risa_GC_free(p);
1.1 noro 23: }
24:
25: void init_gmpq()
26: {
1.4 noro 27: mp_set_memory_functions(Risa_GC_malloc,gc_realloc,gc_free);
1.1 noro 28:
29: mpz_init(ONEMPZ); mpz_set_ui(ONEMPZ,1); MPZTOGZ(ONEMPZ,ONEGZ);
30: }
31:
32: GZ utogz(unsigned int u)
33: {
34: mpz_t z;
35: GZ r;
36:
37: if ( !u ) return 0;
38: mpz_init(z); mpz_set_ui(z,u); MPZTOGZ(z,r); return r;
39: }
40:
41: GZ stogz(int s)
42: {
43: mpz_t z;
44: GZ r;
45:
46: if ( !s ) return 0;
47: mpz_init(z); mpz_set_si(z,s); MPZTOGZ(z,r); return r;
48: }
49:
50: GQ mpqtogzq(mpq_t a)
51: {
52: GZ z;
53: GQ q;
54:
55: if ( INTMPQ(a) ) {
56: MPZTOGZ(mpq_numref(a),z); return (GQ)z;
57: } else {
58: MPQTOGQ(a,q); return q;
59: }
60: }
61:
62: GZ ztogz(Q a)
63: {
64: mpz_t z;
65: mpq_t b;
66: N nm;
67: GZ s;
68:
69: if ( !a ) return 0;
70: nm = NM(a);
71: mpz_init(z);
72: mpz_import(z,PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));
73: if ( SGN(a)<0 ) mpz_neg(z,z);
74: MPZTOGZ(z,s);
75: return s;
76: }
77:
78: Q gztoz(GZ a)
79: {
80: N nm;
81: Q q;
82: int sgn;
83: size_t len;
84:
85: if ( !a ) return 0;
86: len = WORDSIZE_IN_N(BDY((GZ)a)); nm = NALLOC(len);
1.8 ! noro 87: fprintf(stderr,"%d ",len);
1.1 noro 88: mpz_export(BD(nm),&len,-1,sizeof(int),0,0,BDY((GZ)a));
89: PL(nm) = len;
90: sgn = mpz_sgn(BDY((GZ)a)); NTOQ(nm,sgn,q);
91: return q;
92: }
93:
1.5 noro 94: void dupgz(GZ a,GZ *b)
95: {
96: mpz_t t;
97:
98: if ( !a ) *b = a;
99: else {
100: mpz_init(t); mpz_set(t,BDY(a)); MPZTOGZ(t,*b);
101: }
102: }
103:
1.1 noro 104: int n_bits_gz(GZ a)
105: {
106: return a ? mpz_sizeinbase(BDY(a),2) : 0;
107: }
108:
109: GQ qtogq(Q a)
110: {
111: mpz_t z;
112: mpq_t b;
113: N nm,dn;
114: GZ s;
115: GQ r;
116:
117: if ( !a ) return 0;
118: if ( INT(a) ) {
119: nm = NM(a);
120: mpz_init(z);
121: mpz_import(z,PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));
122: if ( SGN(a)<0 ) mpz_neg(z,z);
123: MPZTOGZ(z,s);
124: return (GQ)s;
125: } else {
126: nm = NM(a); dn = DN(a);
127: mpq_init(b);
128: mpz_import(mpq_numref(b),PL(nm),-1,sizeof(BD(nm)[0]),0,0,BD(nm));
129: mpz_import(mpq_denref(b),PL(dn),-1,sizeof(BD(dn)[0]),0,0,BD(dn));
130: if ( SGN(a)<0 ) mpq_neg(b,b);
131: MPQTOGQ(b,r);
132: return r;
133: }
134: }
135:
136: Q gqtoq(GQ a)
137: {
138: N nm,dn;
139: Q q;
140: int sgn;
141: size_t len;
142:
143: if ( !a ) return 0;
144: if ( NID(a) == N_GZ ) {
145: len = WORDSIZE_IN_N(BDY((GZ)a)); nm = NALLOC(len);
146: mpz_export(BD(nm),&len,-1,sizeof(int),0,0,BDY((GZ)a));
147: PL(nm) = len;
148: sgn = mpz_sgn(BDY((GZ)a)); NTOQ(nm,sgn,q);
149: } else {
150: len = WORDSIZE_IN_N(mpq_numref(BDY(a))); nm = NALLOC(len);
151: mpz_export(BD(nm),&len,-1,sizeof(int),0,0,mpq_numref(BDY(a)));
152: PL(nm) = len;
153: len = WORDSIZE_IN_N(mpq_denref(BDY(a))); dn = NALLOC(len);
154: mpz_export(BD(dn),&len,-1,sizeof(int),0,0,mpq_denref(BDY(a)));
155: PL(dn) = len;
156: sgn = mpz_sgn(mpq_numref(BDY(a))); NDTOQ(nm,dn,sgn,q);
157: }
158: return q;
159: }
160:
161: P ptogp(P a)
162: {
163: DCP dc,dcr,dcr0;
164: P b;
165:
166: if ( !a ) return 0;
167: if ( NUM(a) ) return (P)qtogq((Q)a);
168: for ( dc = DC(a), dcr0 = 0; dc; dc = NEXT(dc) ) {
169: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc); COEF(dcr) = (P)ptogp(COEF(dc));
170: }
171: NEXT(dcr) = 0; MKP(VR(a),dcr0,b);
172: return b;
173: }
174:
175: P gptop(P a)
176: {
177: DCP dc,dcr,dcr0;
178: P b;
179:
180: if ( !a ) return 0;
181: if ( NUM(a) ) b = (P)gqtoq((GQ)a);
182: else {
183: for ( dc = DC(a), dcr0 = 0; dc; dc = NEXT(dc) ) {
184: NEXTDC(dcr0,dcr); DEG(dcr) = DEG(dc);
185: COEF(dcr) = (P)gptop(COEF(dc));
186: }
187: NEXT(dcr) = 0; MKP(VR(a),dcr0,b);
188: }
189: return b;
190: }
191:
192: void addgz(GZ n1,GZ n2,GZ *nr)
193: {
194: mpz_t t;
195: int s1,s2;
196:
197: if ( !n1 ) *nr = n2;
198: else if ( !n2 ) *nr = n1;
199: else {
200: mpz_init(t); mpz_add(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nr);
201: }
202: }
203:
204: void subgz(GZ n1,GZ n2,GZ *nr)
205: {
206: mpz_t t;
207:
208: if ( !n1 )
209: if ( !n2 )
210: *nr = 0;
211: else {
212: t[0] = BDY(n2)[0]; mpz_neg(t,t); MPZTOGZ(t,*nr);
213: }
214: else if ( !n2 )
215: *nr = n1;
216: else if ( n1 == n2 )
217: *nr = 0;
218: else {
219: mpz_init(t); mpz_sub(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nr);
220: }
221: }
222:
223: void mulgz(GZ n1,GZ n2,GZ *nr)
224: {
225: mpz_t t;
226:
227: if ( !n1 || !n2 ) *nr = 0;
1.5 noro 228: #if 1
1.1 noro 229: else if ( UNIGZ(n1) ) *nr = n2;
230: else if ( UNIGZ(n2) ) *nr = n1;
231: else if ( MUNIGZ(n1) ) chsgngz(n2,nr);
232: else if ( MUNIGZ(n2) ) chsgngz(n1,nr);
1.5 noro 233: #endif
1.1 noro 234: else {
235: mpz_init(t); mpz_mul(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nr);
236: }
237: }
238:
1.5 noro 239: /* nr += n1*n2 */
240:
241: void muladdtogz(GZ n1,GZ n2,GZ *nr)
242: {
243: GZ t;
244:
245: if ( n1 && n2 ) {
246: if ( !(*nr) ) {
247: NEWGZ(t); mpz_init(BDY(t)); *nr = t;
248: }
249: mpz_addmul(BDY(*nr),BDY(n1),BDY(n2));
250: }
251: }
252:
1.3 noro 253: void mul1gz(GZ n1,int n2,GZ *nr)
254: {
255: mpz_t t;
256:
257: if ( !n1 || !n2 ) *nr = 0;
258: else {
259: mpz_init(t); mpz_mul_ui(t,BDY(n1),(long)n2); MPZTOGZ(t,*nr);
260: }
261: }
262:
1.1 noro 263: void divgz(GZ n1,GZ n2,GZ *nq)
264: {
265: mpz_t t;
266: mpq_t a, b, q;
267:
268: if ( !n2 ) {
269: error("division by 0");
270: *nq = 0;
271: } else if ( !n1 )
272: *nq = 0;
273: else if ( n1 == n2 ) {
274: mpz_init(t); mpz_set_ui(t,1); MPZTOGZ(t,*nq);
275: } else {
276: MPZTOMPQ(BDY(n1),a); MPZTOMPQ(BDY(n2),b);
277: mpq_init(q); mpq_div(q,a,b); *nq = (GZ)mpqtogzq(q);
278: }
279: }
280:
1.5 noro 281: void remgz(GZ n1,GZ n2,GZ *nr)
282: {
283: mpz_t r;
284:
285: if ( !n2 ) {
286: error("division by 0");
287: *nr = 0;
288: } else if ( !n1 || n1 == n2 )
289: *nr = 0;
290: else {
291: mpz_init(r);
292: mpz_mod(r,BDY(n1),BDY(n2));
293: if ( !mpz_sgn(r) ) *nr = 0;
294: else MPZTOGZ(r,*nr);
295: }
296: }
297:
1.1 noro 298: void divqrgz(GZ n1,GZ n2,GZ *nq,GZ *nr)
299: {
300: mpz_t t, a, b, q, r;
301:
302: if ( !n2 ) {
303: error("division by 0");
304: *nq = 0; *nr = 0;
305: } else if ( !n1 ) {
306: *nq = 0; *nr = 0;
307: } else if ( n1 == n2 ) {
308: mpz_init(t); mpz_set_ui(t,1); MPZTOGZ(t,*nq); *nr = 0;
309: } else {
310: mpz_init(q); mpz_init(r);
311: mpz_fdiv_qr(q,r,BDY(n1),BDY(n2));
312: if ( !mpz_sgn(q) ) *nq = 0;
313: else MPZTOGZ(q,*nq);
314: if ( !mpz_sgn(r) ) *nr = 0;
315: else MPZTOGZ(r,*nr);
316: }
317: }
318:
319: void divsgz(GZ n1,GZ n2,GZ *nq)
320: {
321: mpz_t t;
322: mpq_t a, b, q;
323:
324: if ( !n2 ) {
325: error("division by 0");
326: *nq = 0;
327: } else if ( !n1 )
328: *nq = 0;
329: else if ( n1 == n2 ) {
330: mpz_init(t); mpz_set_ui(t,1); MPZTOGZ(t,*nq);
331: } else {
332: mpz_init(t); mpz_divexact(t,BDY(n1),BDY(n2)); MPZTOGZ(t,*nq);
333: }
334: }
335:
336: void chsgngz(GZ n,GZ *nr)
337: {
338: mpz_t t;
339:
340: if ( !n )
341: *nr = 0;
342: else {
343: t[0] = BDY(n)[0]; mpz_neg(t,t); MPZTOGZ(t,*nr);
344: }
345: }
346:
347: void pwrgz(GZ n1,Q n,GZ *nr)
348: {
349: mpq_t t,q;
350: mpz_t z;
351: GQ p,r;
352:
353: if ( !n || UNIGZ(n1) ) *nr = ONEGZ;
354: else if ( !n1 ) *nr = 0;
355: else if ( !INT(n) ) {
356: error("can't calculate fractional power."); *nr = 0;
357: } else if ( MUNIGZ(n1) ) *nr = BD(NM(n))[0]%2 ? n1 : ONEGZ;
358: else if ( PL(NM(n)) > 1 ) {
359: error("exponent too big."); *nr = 0;
360: } else if ( NID(n1)==N_GZ && SGN(n)>0 ) {
361: mpz_init(z); mpz_pow_ui(z,BDY(n1),QTOS(n)); MPZTOGZ(z,*nr);
362: } else {
363: MPZTOMPQ(BDY(n1),q); MPQTOGQ(q,r);
364: pwrgq(r,n,&p); *nr = (GZ)p;
365: }
366: }
367:
368: int cmpgz(GZ q1,GZ 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 ( (sgn = mpz_sgn(BDY(q1))) != mpz_sgn(BDY(q2)) )
380: return sgn;
381: else {
382: sgn = mpz_cmp(BDY(q1),BDY(q2));
383: if ( sgn > 0 ) return 1;
384: else if ( sgn < 0 ) return -1;
385: else return 0;
386: }
387: }
388:
389: void gcdgz(GZ n1,GZ n2,GZ *nq)
390: {
391: mpz_t t;
392:
393: if ( !n1 ) *nq = n2;
394: else if ( !n2 ) *nq = n1;
395: else {
396: mpz_init(t); mpz_gcd(t,BDY(n1),BDY(n2));
397: MPZTOGZ(t,*nq);
398: }
399: }
400:
1.5 noro 401: void invgz(GZ n1,GZ *nq)
402: {
403: mpz_t t;
404:
405: mpz_init(t); mpz_invert(t,BDY(n1),BDY(current_mod_lf));
406: MPZTOGZ(t,*nq);
407: }
408:
1.1 noro 409: void lcmgz(GZ n1,GZ n2,GZ *nq)
410: {
411: GZ g,t;
412:
413: if ( !n1 || !n2 ) *nq = 0;
414: else {
415: gcdgz(n1,n2,&g); divsgz(n1,g,&t);
416: mulgz(n2,t,nq);
417: }
418: }
419:
420: void gcdvgz(VECT v,GZ *q)
421: {
422: int n,i;
423: GZ *b;
424: GZ g,g1;
425:
426: n = v->len;
427: b = (GZ *)v->body;
428: g = b[0];
429: for ( i = 1; i < n; i++ ) {
430: gcdgz(g,b[i],&g1); g = g1;
431: }
432: *q = g;
433: }
434:
435: void gcdvgz_estimate(VECT v,GZ *q)
436: {
437: int n,m,i;
438: GZ s,t,u;
439: GZ *b;
440:
441: n = v->len;
442: b = (GZ *)v->body;
443: if ( n == 1 ) {
444: if ( mpz_sgn(BDY(b[0]))<0 ) chsgngz(b[0],q);
445: else *q = b[0];
446: }
447: m = n/2;
448: for ( i = 0, s = 0; i < m; i++ ) {
449: if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subgz(s,b[i],&u);
450: else addgz(s,b[i],&u);
451: s = u;
452: }
453: for ( i = 0, t = 0; i < n; i++ ) {
454: if ( b[i] && mpz_sgn(BDY(b[i]))<0 ) subgz(t,b[i],&u);
455: else addgz(t,b[i],&u);
456: t = u;
457: }
458: gcdgz(s,t,q);
459: }
460:
461: void addgq(GQ n1,GQ n2,GQ *nr)
462: {
463: mpq_t q1,q2,t;
464:
465: if ( !n1 ) *nr = n2;
466: else if ( !n2 ) *nr = n1;
467: else {
468: if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);
469: else q1[0] = BDY(n1)[0];
470: if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);
471: else q2[0] = BDY(n2)[0];
472: mpq_init(t); mpq_add(t,q1,q2); *nr = mpqtogzq(t);
473: }
474: }
475:
476: void subgq(GQ n1,GQ n2,GQ *nr)
477: {
478: mpq_t q1,q2,t;
479:
480: if ( !n1 )
481: if ( !n2 ) *nr = 0;
482: else {
483: if ( NID(n1) == N_GZ ) chsgngz((GZ)n1,(GZ *)nr);
484: else {
485: mpq_init(t); mpq_neg(t,BDY(n2)); MPQTOGQ(t,*nr);
486: }
487: }
488: else if ( !n2 ) *nr = n1;
489: else if ( n1 == n2 ) *nr = 0;
490: else {
491: if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);
492: else q1[0] = BDY(n1)[0];
493: if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);
494: else q2[0] = BDY(n2)[0];
495: mpq_init(t); mpq_sub(t,q1,q2); *nr = mpqtogzq(t);
496: }
497: }
498:
499: void mulgq(GQ n1,GQ n2,GQ *nr)
500: {
501: mpq_t t,q1,q2;
502:
503: if ( !n1 || !n2 ) *nr = 0;
504: else {
505: if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);
506: else q1[0] = BDY(n1)[0];
507: if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);
508: else q2[0] = BDY(n2)[0];
509: mpq_init(t); mpq_mul(t,q1,q2); *nr = mpqtogzq(t);
510: }
511: }
512:
513: void divgq(GQ n1,GQ n2,GQ *nq)
514: {
515: mpq_t t,q1,q2;
516:
517: if ( !n2 ) {
518: error("division by 0");
519: *nq = 0;
520: return;
521: } else if ( !n1 ) *nq = 0;
522: else if ( n1 == n2 ) *nq = (GQ)ONEGZ;
523: else {
524: if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);
525: else q1[0] = BDY(n1)[0];
526: if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);
527: else q2[0] = BDY(n2)[0];
528: mpq_init(t); mpq_div(t,q1,q2); *nq = mpqtogzq(t);
529: }
530: }
531:
532: void chsgngq(GQ n,GQ *nr)
533: {
534: mpq_t t;
535:
536: if ( !n ) *nr = 0;
537: else if ( NID(n) == N_GZ ) chsgngz((GZ)n,(GZ *)nr);
538: else {
539: mpq_init(t); mpq_neg(t,BDY(n)); MPQTOGQ(t,*nr);
540: }
541: }
542:
543: void pwrgq(GQ n1,Q n,GQ *nr)
544: {
545: int e;
546: mpz_t nm,dn;
547: mpq_t t;
548:
549: if ( !n || UNIGZ((GZ)n1) || UNIGQ(n1) ) *nr = (GQ)ONEGZ;
550: else if ( !n1 ) *nr = 0;
551: else if ( !INT(n) ) {
552: error("can't calculate fractional power."); *nr = 0;
553: } else if ( PL(NM(n)) > 1 ) {
554: error("exponent too big."); *nr = 0;
555: } else {
556: e = QTOS(n);
557: if ( e < 0 ) {
558: e = -e;
559: if ( NID(n1)==N_GZ ) {
560: nm[0] = ONEMPZ[0];
561: dn[0] = BDY((GZ)n1)[0];
562: } else {
563: nm[0] = mpq_denref(BDY(n1))[0]; dn[0] = mpq_numref(BDY(n1))[0];
564: }
565: } else {
566: if ( NID(n1)==N_GZ ) {
567: nm[0] = BDY((GZ)n1)[0]; dn[0] = ONEMPZ[0];
568: } else {
569: nm[0] = mpq_numref(BDY(n1))[0]; dn[0] = mpq_denref(BDY(n1))[0];
570: }
571: }
572: mpq_init(t);
573: mpz_pow_ui(mpq_numref(t),nm,e); mpz_pow_ui(mpq_denref(t),dn,e);
574: *nr = mpqtogzq(t);
575: }
576: }
577:
578: int cmpgq(GQ n1,GQ n2)
579: {
580: mpq_t q1,q2;
581: int sgn;
582:
583: if ( !n1 )
584: if ( !n2 ) return 0;
585: else return (NID(n2)==N_GZ) ? -mpz_sgn(BDY((GZ)n2)) : -mpq_sgn(BDY(n2));
586: if ( !n2 ) return (NID(n1)==N_GZ) ? mpz_sgn(BDY((GZ)n1)) : mpq_sgn(BDY(n1));
587: else if ( (sgn = mpq_sgn(BDY(n1))) != mpq_sgn(BDY(n2)) ) return sgn;
588: else {
589: if ( NID(n1) == N_GZ ) MPZTOMPQ(BDY((GZ)n1),q1);
590: else q1[0] = BDY(n1)[0];
591: if ( NID(n2) == N_GZ ) MPZTOMPQ(BDY((GZ)n2),q2);
592: else q2[0] = BDY(n2)[0];
593: sgn = mpq_cmp(q1,q2);
594: if ( sgn > 0 ) return 1;
595: else if ( sgn < 0 ) return -1;
596: else return 0;
597: }
598: }
599:
600: void mkgwc(int k,int l,GZ *t)
601: {
602: mpz_t a,b,q,nm,z,u;
603: int i,n;
604:
605: n = MIN(k,l);
606: mpz_init_set_ui(z,1);
607: mpz_init(u); mpz_set(u,z); MPZTOGZ(u,t[0]);
608: mpz_init(a); mpz_init(b); mpz_init(nm);
609: for ( i = 1; i <= n; i++ ) {
610: mpz_set_ui(a,k-i+1); mpz_set_ui(b,l-i+1); mpz_mul(nm,a,b);
611: mpz_mul(z,BDY(t[i-1]),nm); mpz_fdiv_q_ui(z,z,i);
612: mpz_init(u); mpz_set(u,z); MPZTOGZ(u,t[i]);
613: }
614: }
615:
616: void gz_lgp(P p,GZ *g,GZ *l);
617:
618: void gz_ptozp(P p,int sgn,GQ *c,P *pr)
619: {
620: GZ nm,dn;
621:
622: if ( !p ) {
623: *c = 0; *pr = 0;
624: } else {
625: gz_lgp(p,&nm,&dn);
626: divgz(nm,dn,(GZ *)c);
627: divsp(CO,p,(P)*c,pr);
628: }
629: }
630:
631: void gz_lgp(P p,GZ *g,GZ *l)
632: {
633: DCP dc;
634: GZ g1,g2,l1,l2,l3,l4;
635:
636: if ( NUM(p) ) {
637: if ( NID((GZ)p)==N_GZ ) {
638: MPZTOGZ(BDY((GZ)p),*g);
639: *l = ONEGZ;
640: } else {
641: MPZTOGZ(mpq_numref(BDY((GQ)p)),*g);
642: MPZTOGZ(mpq_denref(BDY((GQ)p)),*l);
643: }
644: } else {
645: dc = DC(p); gz_lgp(COEF(dc),g,l);
646: for ( dc = NEXT(dc); dc; dc = NEXT(dc) ) {
647: gz_lgp(COEF(dc),&g1,&l1); gcdgz(*g,g1,&g2); *g = g2;
648: gcdgz(*l,l1,&l2); mulgz(*l,l1,&l3); divgz(l3,l2,l);
649: }
650: }
651: }
652:
653: void gz_qltozl(GQ *w,int n,GZ *dvr)
654: {
655: GZ nm,dn;
656: GZ g,g1,l1,l2,l3;
657: GQ c;
658: int i;
659: struct oVECT v;
660:
661: for ( i = 0; i < n; i++ )
662: if ( w[i] && NID(w[i])==N_GQ )
663: break;
664: if ( i == n ) {
665: v.id = O_VECT; v.len = n; v.body = (pointer *)w;
666: gcdvgz(&v,dvr); return;
667: }
668: for ( i = 0; !w[i]; i++ );
669: c = w[i];
670: if ( NID(c)==N_GQ ) {
671: MPZTOGZ(mpq_numref(BDY(c)),nm); MPZTOGZ(mpq_denref(BDY(c)),dn);
672: } else {
673: MPZTOGZ(BDY((GZ)c),nm); dn = ONEGZ;
674: }
675: for ( i++; i < n; i++ ) {
676: c = w[i];
677: if ( !c ) continue;
678: if ( NID(c)==N_GQ ) {
679: MPZTOGZ(mpq_numref(BDY(c)),g1); MPZTOGZ(mpq_denref(BDY(c)),l1);
680: } else {
681: MPZTOGZ(BDY((GZ)c),g1); l1 = ONEGZ;
682: }
683: gcdgz(nm,g1,&g); nm = g;
684: gcdgz(dn,l1,&l2); mulgz(dn,l1,&l3); divgz(l3,l2,&dn);
685: }
686: divgz(nm,dn,dvr);
687: }
688:
689: int gz_bits(GQ q)
690: {
691: if ( !q ) return 0;
692: else if ( NID(q)==N_Q )
693: return n_bits(NM((Q)q))+(INT((Q)q)?0:n_bits(DN((Q)q)));
694: else if ( NID(q)==N_GZ ) return mpz_sizeinbase(BDY((GZ)q),2);
695: else
696: return mpz_sizeinbase(mpq_numref(BDY(q)),2)
697: + mpz_sizeinbase(mpq_denref(BDY(q)),2);
698: }
699:
700: int gzp_mag(P p)
701: {
702: int s;
703: DCP dc;
704:
705: if ( !p ) return 0;
706: else if ( OID(p) == O_N ) return gz_bits((GQ)p);
707: else {
708: for ( dc = DC(p), s = 0; dc; dc = NEXT(dc) ) s += gzp_mag(COEF(dc));
709: return s;
710: }
711: }
712:
713: void makesubstgz(VL v,NODE *s)
714: {
715: NODE r,r0;
716: GZ q;
717: unsigned int n;
718:
719: for ( r0 = 0; v; v = NEXT(v) ) {
720: NEXTNODE(r0,r); BDY(r) = (pointer)v->v;
721: #if defined(_PA_RISC1_1)
722: n = mrand48()&BMASK; q = utogz(n);
723: #else
724: n = random(); q = utogz(n);
725: #endif
726: NEXTNODE(r0,r); BDY(r) = (pointer)q;
727: }
728: if ( r0 ) NEXT(r) = 0;
729: *s = r0;
730: }
731:
732: unsigned int remgq(GQ a,unsigned int mod)
733: {
734: unsigned int c,nm,dn;
735: mpz_t r;
736:
737: if ( !a ) return 0;
738: else if ( NID(a)==N_GZ ) {
739: mpz_init(r);
740: c = mpz_fdiv_r_ui(r,BDY((GZ)a),mod);
741: } else {
742: mpz_init(r);
743: nm = mpz_fdiv_r_ui(r,mpq_numref(BDY(a)),mod);
744: dn = mpz_fdiv_r_ui(r,mpq_denref(BDY(a)),mod);
745: dn = invm(dn,mod);
746: DMAR(nm,dn,0,mod,c);
747: }
748: return c;
749: }
750:
751: extern int DP_Print;
752:
753: #define GZ_F4_INTRAT_PERIOD 8
754:
755: int gz_generic_gauss_elim(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp)
756: {
757: int **wmat;
758: GZ **bmat,**tmat,*bmi,*tmi;
759: GZ q,m1,m2,m3,s,u;
760: int *wmi,*colstat,*wcolstat,*rind,*cind;
761: int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv;
762: MAT r,crmat;
763: int ret;
764:
765: bmat = (GZ **)mat->body;
766: row = mat->row; col = mat->col;
767: wmat = (int **)almat(row,col);
768: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
769: wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
770: for ( ind = 0; ; ind++ ) {
771: if ( DP_Print ) {
772: fprintf(asir_out,"."); fflush(asir_out);
773: }
774: md = get_lprime(ind);
775: for ( i = 0; i < row; i++ )
776: for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
777: wmi[j] = remgq((GQ)bmi[j],md);
778: rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);
779: if ( !ind ) {
780: RESET:
781: m1 = utogz(md);
782: rank0 = rank;
783: bcopy(wcolstat,colstat,col*sizeof(int));
784: MKMAT(crmat,rank,col-rank);
785: MKMAT(r,rank,col-rank); *nm = r;
786: tmat = (GZ **)crmat->body;
787: for ( i = 0; i < rank; i++ )
788: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
789: if ( !colstat[j] ) tmi[k++] = utogz(wmi[j]);
790: } else {
791: if ( rank < rank0 ) {
792: if ( DP_Print ) {
793: fprintf(asir_out,"lower rank matrix; continuing...\n");
794: fflush(asir_out);
795: }
796: continue;
797: } else if ( rank > rank0 ) {
798: if ( DP_Print ) {
799: fprintf(asir_out,"higher rank matrix; resetting...\n");
800: fflush(asir_out);
801: }
802: goto RESET;
803: } else {
804: for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
805: if ( j < col ) {
806: if ( DP_Print ) {
807: fprintf(asir_out,"inconsitent colstat; resetting...\n");
808: fflush(asir_out);
809: }
810: goto RESET;
811: }
812: }
813:
814: inv = invm(remgq((GQ)m1,md),md);
815: m2 = utogz(md); mulgz(m1,m2,&m3);
816: for ( i = 0; i < rank; i++ )
817: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
818: if ( !colstat[j] ) {
819: if ( tmi[k] ) {
820: /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
821: t = remgq((GQ)tmi[k],md);
822: if ( wmi[j] >= t ) t = wmi[j]-t;
823: else t = md-(t-wmi[j]);
824: DMAR(t,inv,0,md,t1)
825: u = utogz(t1); mulgz(m1,u,&s);
826: addgz(tmi[k],s,&u); tmi[k] = u;
827: } else if ( wmi[j] ) {
828: /* f3 = m1*(m1 mod m2)^(-1)*f2 */
829: DMAR(wmi[j],inv,0,md,t)
830: u = utogz(t); mulgz(m1,u,&s); tmi[k] = s;
831: }
832: k++;
833: }
834: m1 = m3;
835: if ( ind % GZ_F4_INTRAT_PERIOD )
836: ret = 0;
837: else
838: ret = gz_intmtoratm(crmat,m1,*nm,dn);
839: if ( ret ) {
840: *rindp = rind = (int *)MALLOC_ATOMIC(rank*sizeof(int));
841: *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
842: for ( j = k = l = 0; j < col; j++ )
843: if ( colstat[j] ) rind[k++] = j;
844: else cind[l++] = j;
845: if ( gz_gensolve_check(mat,*nm,*dn,rind,cind) )
846: return rank;
847: }
848: }
849: }
850: }
851:
852: int gz_generic_gauss_elim2(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp)
853: {
854:
855: MAT full;
856: GZ **bmat,**b;
857: GZ *bmi;
858: GZ dn0;
859: int row,col,md,i,j,rank,ret;
860: int **wmat;
861: int *wmi;
862: int *colstat,*rowstat;
863:
864: bmat = (GZ **)mat->body;
865: row = mat->row; col = mat->col;
866: wmat = (int **)almat(row,col);
867: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
868: rowstat = (int *)MALLOC_ATOMIC(row*sizeof(int));
869: /* XXX */
870: md = get_lprime(0);
871: for ( i = 0; i < row; i++ )
872: for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
873: wmi[j] = remgq((GQ)bmi[j],md);
874: rank = generic_gauss_elim_mod2(wmat,row,col,md,colstat,rowstat);
875: b = (GZ **)MALLOC(rank*sizeof(GZ));
876: for ( i = 0; i < rank; i++ ) b[i] = bmat[rowstat[i]];
877: NEWMAT(full); full->row = rank; full->col = col; full->body = (pointer **)b;
878: ret = gz_generic_gauss_elim_full(full,nm,dn,rindp,cindp);
879: if ( !ret ) {
880: rank = gz_generic_gauss_elim(mat,nm,&dn0,rindp,cindp);
881: for ( i = 0; i < rank; i++ ) dn[i] = dn0;
882: }
883: return rank;
884: }
885:
886: int gz_generic_gauss_elim_full(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp)
887: {
888: int **wmat;
889: int *stat;
890: GZ **bmat,**tmat,*bmi,*tmi,*ri;
891: GZ q,m1,m2,m3,s,u;
892: int *wmi,*colstat,*wcolstat,*rind,*cind;
893: int row,col,ind,md,i,j,k,l,t,t1,rank,rank0,inv,h;
894: MAT r,crmat;
895: int ret,initialized,done;
896:
897: initialized = 0;
898: bmat = (GZ **)mat->body;
899: row = mat->row; col = mat->col;
900: wmat = (int **)almat(row,col);
901: stat = (int *)MALLOC_ATOMIC(row*sizeof(int));
902: for ( i = 0; i < row; i++ ) stat[i] = 0;
903: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
904: wcolstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
905: for ( ind = 0; ; ind++ ) {
906: if ( DP_Print ) {
907: fprintf(asir_out,"."); fflush(asir_out);
908: }
909: md = get_lprime(ind);
910: for ( i = 0; i < row; i++ )
911: for ( j = 0, bmi = bmat[i], wmi = wmat[i]; j < col; j++ )
912: wmi[j] = remgq((GQ)bmi[j],md);
913: rank = generic_gauss_elim_mod(wmat,row,col,md,wcolstat);
914: if ( rank < row ) continue;
915: if ( !initialized ) {
916: m1 = utogz(md);
917: bcopy(wcolstat,colstat,col*sizeof(int));
918: MKMAT(crmat,row,col-row);
919: MKMAT(r,row,col-row); *nm = r;
920: tmat = (GZ **)crmat->body;
921: for ( i = 0; i < row; i++ )
922: for ( j = k = 0, tmi = tmat[i], wmi = wmat[i]; j < col; j++ )
923: if ( !colstat[j] ) tmi[k++] = utogz(wmi[j]);
924: initialized = 1;
925: } else {
926: for ( j = 0; (j<col) && (colstat[j]==wcolstat[j]); j++ );
927: if ( j < col ) continue;
928:
929: inv = invm(remgq((GQ)m1,md),md);
930: m2 = utogz(md); mulgz(m1,m2,&m3);
931: for ( i = 0; i < row; i++ )
932: switch ( stat[i] ) {
933: case 1:
934: /* consistency check */
935: ri = (GZ *)BDY(r)[i]; wmi = wmat[i];
936: for ( j = 0; j < col; j++ ) if ( colstat[j] ) break;
937: h = md-remgq((GQ)dn[i],md);
938: for ( j++, k = 0; j < col; j++ )
939: if ( !colstat[j] ) {
940: t = remgq((GQ)ri[k],md);
941: DMAR(wmi[i],h,t,md,t1);
942: if ( t1 ) break;
943: }
944: if ( j == col ) { stat[i]++; break; }
945: else {
946: /* fall to the case 0 */
947: stat[i] = 0;
948: }
949: case 0:
950: tmi = tmat[i]; wmi = wmat[i];
951: for ( j = k = 0; j < col; j++ )
952: if ( !colstat[j] ) {
953: if ( tmi[k] ) {
954: /* f3 = f1+m1*(m1 mod m2)^(-1)*(f2 - f1 mod m2) */
955: t = remgq((GQ)tmi[k],md);
956: if ( wmi[j] >= t ) t = wmi[j]-t;
957: else t = md-(t-wmi[j]);
958: DMAR(t,inv,0,md,t1)
959: u = utogz(t1); mulgz(m1,u,&s);
960: addgz(tmi[k],s,&u); tmi[k] = u;
961: } else if ( wmi[j] ) {
962: /* f3 = m1*(m1 mod m2)^(-1)*f2 */
963: DMAR(wmi[j],inv,0,md,t)
964: u = utogz(t); mulgz(m1,u,&s); tmi[k] = s;
965: }
966: k++;
967: }
968: break;
969: case 2: default:
970: break;
971: }
972: m1 = m3;
973: if ( ind % 4 )
974: ret = 0;
975: else
976: ret = gz_intmtoratm2(crmat,m1,*nm,dn,stat);
977: if ( ret ) {
978: *rindp = rind = (int *)MALLOC_ATOMIC(row*sizeof(int));
979: *cindp = cind = (int *)MALLOC_ATOMIC((col-row)*sizeof(int));
980: for ( j = k = l = 0; j < col; j++ )
981: if ( colstat[j] ) rind[k++] = j;
982: else cind[l++] = j;
983: return gz_gensolve_check2(mat,*nm,dn,rind,cind);
984: }
985: }
986: }
987: }
988:
989: int gz_generic_gauss_elim_direct(MAT mat,MAT *nm,GZ *dn,int **rindp,int **cindp){
990: GZ **bmat,*s;
991: GZ u,v,w,x,d,t,y;
992: int row,col,i,j,k,l,m,rank;
993: int *colstat,*colpos,*cind;
994: MAT r,in;
995:
996: row = mat->row; col = mat->col;
997: MKMAT(in,row,col);
998: for ( i = 0; i < row; i++ )
999: for ( j = 0; j < col; j++ ) in->body[i][j] = mat->body[i][j];
1000: bmat = (GZ **)in->body;
1001: colstat = (int *)MALLOC_ATOMIC(col*sizeof(int));
1002: *rindp = colpos = (int *)MALLOC_ATOMIC(row*sizeof(int));
1003: for ( j = 0, rank = 0, d = ONEGZ; j < col; j++ ) {
1004: for ( i = rank; i < row && !bmat[i][j]; i++ );
1005: if ( i == row ) { colstat[j] = 0; continue; }
1006: else { colstat[j] = 1; colpos[rank] = j; }
1007: if ( i != rank )
1008: for ( k = j; k < col; k++ ) {
1009: t = bmat[i][k]; bmat[i][k] = bmat[rank][k]; bmat[rank][k] = t;
1010: }
1011: for ( i = rank+1, v = bmat[rank][j]; i < row; i++ )
1012: for ( k = j, u = bmat[i][j]; k < col; k++ ) {
1013: mulgz(bmat[i][k],v,&w); mulgz(bmat[rank][k],u,&x);
1014: subgz(w,x,&y); divsgz(y,d,&bmat[i][k]);
1015: }
1016: d = v; rank++;
1017: }
1018: *dn = d;
1019: s = (GZ *)MALLOC(col*sizeof(GZ));
1020: for ( i = rank-1; i >= 0; i-- ) {
1021: for ( k = colpos[i]; k < col; k++ ) mulgz(bmat[i][k],d,&s[k]);
1022: for ( m = rank-1; m > i; m-- ) {
1023: for ( k = colpos[m], u = bmat[i][k]; k < col; k++ ) {
1024: mulgz(bmat[m][k],u,&w); subgz(s[k],w,&x); s[k] = x;
1025: }
1026: }
1027: for ( k = colpos[i], u = bmat[i][k]; k < col; k++ )
1028: divgz(s[k],u,&bmat[i][k]);
1029: }
1030: *cindp = cind = (int *)MALLOC_ATOMIC((col-rank)*sizeof(int));
1031: MKMAT(r,rank,col-rank); *nm = r;
1032: for ( j = 0, k = 0; j < col; j++ )
1033: if ( !colstat[j] ) {
1034: cind[k] = j;
1035: for ( i = 0; i < rank; i++ ) r->body[i][k] = bmat[i][j];
1036: k++;
1037: }
1038: return rank;
1039: }
1040:
1041: int gz_intmtoratm(MAT mat,GZ md,MAT nm,GZ *dn)
1042: {
1043: GZ t,s,b,dn0,dn1,nm1,q,u,unm,udn,dmy;
1044: int i,j,k,l,row,col,sgn,ret;
1045: GZ **rmat,**tmat,*tmi,*nmk;
1046:
1047: if ( UNIGZ(md) )
1048: return 0;
1049: row = mat->row; col = mat->col;
1050: bshiftgz(md,1,&t);
1051: isqrtgz(t,&s);
1052: bshiftgz(s,64,&b);
1053: if ( !b ) b = ONEGZ;
1054: dn0 = ONEGZ;
1055: tmat = (GZ **)mat->body;
1056: rmat = (GZ **)nm->body;
1057: for ( i = 0; i < row; i++ )
1058: for ( j = 0, tmi = tmat[i]; j < col; j++ )
1059: if ( tmi[j] ) {
1060: mulgz(tmi[j],dn0,&s);
1061: divqrgz(s,md,&dmy,&u);
1062: ret = gz_inttorat(u,md,b,&sgn,&unm,&udn);
1063: if ( !ret ) return 0;
1064: else {
1065: if ( sgn < 0 ) chsgngz(unm,&nm1);
1066: else nm1 = unm;
1067: dn1 = udn;
1068: if ( !UNIGZ(dn1) ) {
1069: for ( k = 0; k < i; k++ )
1070: for ( l = 0, nmk = rmat[k]; l < col; l++ ) {
1071: mulgz(nmk[l],dn1,&q); nmk[l] = q;
1072: }
1073: for ( l = 0, nmk = rmat[i]; l < j; l++ ) {
1074: mulgz(nmk[l],dn1,&q); nmk[l] = q;
1075: }
1076: }
1077: rmat[i][j] = nm1;
1078: mulgz(dn0,dn1,&q); dn0 = q;
1079: }
1080: }
1081: *dn = dn0;
1082: return 1;
1083: }
1084:
1085: int gz_intmtoratm2(MAT mat,GZ md,MAT nm,GZ *dn,int *stat)
1086: {
1087: int row,col,i,j,ret;
1088: GZ dn0,dn1,t,s,b;
1089: GZ *w,*tmi;
1090: GZ **tmat;
1091:
1092: bshiftgz(md,1,&t);
1093: isqrtgz(t,&s);
1094: bshiftgz(s,64,&b);
1095: tmat = (GZ **)mat->body;
1096: if ( UNIGZ(md) ) return 0;
1097: row = mat->row; col = mat->col;
1098: dn0 = ONEGZ;
1099: for ( i = 0; i < row; i++ )
1100: if ( cmpgz(dn[i],dn0) > 0 ) dn0 = dn[i];
1101: w = (GZ *)MALLOC(col*sizeof(GZ));
1102: for ( i = 0; i < row; i++ )
1103: if ( stat[i] == 0 ) {
1104: for ( j = 0, tmi = tmat[i]; j < col; j++ )
1105: mulgz(tmi[j],dn0,&w[j]);
1106: ret = gz_intvtoratv(w,col,md,b,BDY(nm)[i],&dn[i]);
1107: if ( ret ) {
1108: stat[i] = 1;
1109: mulgz(dn0,dn[i],&t); dn[i] = t; dn0 = t;
1110: }
1111: }
1112: for ( i = 0; i < row; i++ ) if ( !stat[i] ) break;
1113: if ( i == row ) return 1;
1114: else return 0;
1115: }
1116:
1117: int gz_intvtoratv(GZ *v,int n,GZ md,GZ b,GZ *nm,GZ *dn)
1118: {
1119: GZ dn0,dn1,q,s,u,nm1,unm,udn,dmy;
1120: GZ *nmk;
1121: int j,l,col,ret,sgn;
1122:
1123: for ( j = 0; j < n; j++ ) nm[j] = 0;
1124: dn0 = ONEGZ;
1125: for ( j = 0; j < n; j++ ) {
1126: if ( !v[j] ) continue;
1127: mulgz(v[j],dn0,&s);
1128: divqrgz(s,md,&dmy,&u);
1129: ret = gz_inttorat(u,md,b,&sgn,&unm,&udn);
1130: if ( !ret ) return 0;
1131: if ( sgn < 0 ) chsgngz(unm,&nm1);
1132: else nm1 = unm;
1133: dn1 = udn;
1134: if ( !UNIGZ(dn1) )
1135: for ( l = 0; l < j; l++ ) {
1136: mulgz(nm[l],dn1,&q); nm[l] = q;
1137: }
1138: nm[j] = nm1;
1139: mulgz(dn0,dn1,&q); dn0 = q;
1140: }
1141: *dn = dn0;
1142: return 1;
1143: }
1144:
1145: /* assuming 0 < c < m */
1146:
1147: int gz_inttorat(GZ c,GZ m,GZ b,int *sgnp,GZ *nmp,GZ *dnp)
1148: {
1149: GZ qq,t,u1,v1,r1;
1150: GZ q,u2,v2,r2;
1151:
1152: u1 = 0; v1 = ONEGZ; u2 = m; v2 = c;
1153: while ( cmpgz(v2,b) >= 0 ) {
1154: divqrgz(u2,v2,&q,&r2); u2 = v2; v2 = r2;
1155: mulgz(q,v1,&t); subgz(u1,t,&r1); u1 = v1; v1 = r1;
1156: }
1157: if ( cmpgz(v1,b) >= 0 ) return 0;
1158: else {
1159: *nmp = v2;
1160: if ( mpz_sgn(BDY(v1))<0 ) {
1161: *sgnp = -1; chsgngz(v1,dnp);
1162: } else {
1163: *sgnp = 1; *dnp = v1;
1164: }
1165: return 1;
1166: }
1167: }
1168:
1169: extern int f4_nocheck;
1170:
1171: int gz_gensolve_check(MAT mat,MAT nm,GZ dn,int *rind,int *cind)
1172: {
1173: int row,col,rank,clen,i,j,k,l;
1174: GZ s,t;
1175: GZ *w;
1176: GZ *mati,*nmk;
1177:
1178: if ( f4_nocheck ) return 1;
1179: row = mat->row; col = mat->col; rank = nm->row; clen = nm->col;
1180: w = (GZ *)MALLOC(clen*sizeof(GZ));
1181: for ( i = 0; i < row; i++ ) {
1182: mati = (GZ *)mat->body[i];
1183: bzero(w,clen*sizeof(GZ));
1184: for ( k = 0; k < rank; k++ )
1185: for ( l = 0, nmk = (GZ *)nm->body[k]; l < clen; l++ ) {
1186: mulgz(mati[rind[k]],nmk[l],&t); addgz(w[l],t,&s); w[l] = s;
1187: }
1188: for ( j = 0; j < clen; j++ ) {
1189: mulgz(dn,mati[cind[j]],&t);
1190: if ( cmpgz(w[j],t) ) break;
1191: }
1192: if ( j != clen ) break;
1193: }
1194: if ( i != row ) return 0;
1195: else return 1;
1196: }
1197:
1198: int gz_gensolve_check2(MAT mat,MAT nm,GZ *dn,int *rind,int *cind)
1199: {
1200: int row,col,rank,clen,i,j,k,l;
1201: GZ s,t,u,d;
1202: GZ *w,*m;
1203: GZ *mati,*nmk;
1204:
1205: if ( f4_nocheck ) return 1;
1206: row = mat->row; col = mat->col; rank = nm->row; clen = nm->col;
1207: w = (GZ *)MALLOC(clen*sizeof(GZ));
1208: m = (GZ *)MALLOC(clen*sizeof(GZ));
1209: for ( d = dn[0], i = 1; i < rank; i++ ) {
1210: lcmgz(d,dn[i],&t); d = t;
1211: }
1212: for ( i = 0; i < rank; i++ ) divsgz(d,dn[i],&m[i]);
1213: for ( i = 0; i < row; i++ ) {
1214: mati = (GZ *)mat->body[i];
1215: bzero(w,clen*sizeof(GZ));
1216: for ( k = 0; k < rank; k++ ) {
1217: mulgz(mati[rind[k]],m[k],&u);
1218: for ( l = 0, nmk = (GZ *)nm->body[k]; l < clen; l++ ) {
1219: mulgz(u,nmk[l],&t); addgz(w[l],t,&s); w[l] = s;
1220: }
1221: }
1222: for ( j = 0; j < clen; j++ ) {
1223: mulgz(d,mati[cind[j]],&t);
1224: if ( cmpgz(w[j],t) ) break;
1225: }
1226: if ( j != clen ) break;
1227: }
1228: if ( i != row ) return 0;
1229: else return 1;
1230: }
1231:
1232: void isqrtgz(GZ a,GZ *r)
1233: {
1234: int k;
1235: GZ x,t,x2,xh,quo,rem;
1236: Q two;
1237:
1238: if ( !a ) *r = 0;
1239: else if ( UNIGZ(a) ) *r = ONEGZ;
1240: else {
1241: k = gz_bits((GQ)a); /* a <= 2^k-1 */
1242: bshiftgz(ONEGZ,-((k>>1)+(k&1)),&x); /* a <= x^2 */
1243: STOQ(2,two);
1244: while ( 1 ) {
1245: pwrgz(x,two,&t);
1246: if ( cmpgz(t,a) <= 0 ) {
1247: *r = x; return;
1248: } else {
1249: if ( mpz_tstbit(BDY(x),0) ) addgz(x,a,&t);
1250: else t = a;
1251: bshiftgz(x,-1,&x2); divqrgz(t,x2,&quo,&rem);
1252: bshiftgz(x,1,&xh); addgz(quo,xh,&x);
1253: }
1254: }
1255: }
1256: }
1257:
1258: void bshiftgz(GZ a,int n,GZ *r)
1259: {
1260: mpz_t t;
1261:
1262: if ( !a ) *r = 0;
1263: else if ( n == 0 ) *r = a;
1264: else if ( n < 0 ) {
1265: mpz_init(t); mpz_mul_2exp(t,BDY(a),-n); MPZTOGZ(t,*r);
1266: } else {
1267: mpz_init(t); mpz_fdiv_q_2exp(t,BDY(a),n);
1268: if ( !mpz_sgn(t) ) *r = 0;
1269: else MPZTOGZ(t,*r);
1270: }
1271: }
1272:
1.5 noro 1273: void addlf(GZ a,GZ b,GZ *c)
1274: {
1275: GZ t;
1276:
1277: addgz(a,b,c);
1278: if ( !lf_lazy ) {
1279: remgz(*c,current_mod_lf,&t); *c = t;
1280: }
1281: }
1282:
1283: void sublf(GZ a,GZ b,GZ *c)
1284: {
1285: GZ t;
1286:
1287: subgz(a,b,c);
1288: if ( !lf_lazy ) {
1289: remgz(*c,current_mod_lf,&t); *c = t;
1290: }
1291: }
1292:
1293: void mullf(GZ a,GZ b,GZ *c)
1294: {
1295: GZ t;
1296:
1297: mulgz(a,b,c);
1298: if ( !lf_lazy ) {
1299: remgz(*c,current_mod_lf,&t); *c = t;
1300: }
1301: }
1302:
1303: void divlf(GZ a,GZ b,GZ *c)
1304: {
1305: GZ t,inv;
1306:
1307: invgz(b,&inv);
1308: mulgz(a,inv,c);
1309: if ( !lf_lazy ) {
1310: remgz(*c,current_mod_lf,&t); *c = t;
1311: }
1312: }
1313:
1314: void chsgnlf(GZ a,GZ *c)
1315: {
1316: GZ t;
1317:
1318: chsgngz(a,c);
1319: if ( !lf_lazy ) {
1320: remgz(*c,current_mod_lf,&t); *c = t;
1321: }
1322: }
1323:
1324: void lmtolf(LM a,GZ *b)
1325: {
1326: Q q;
1327:
1.7 noro 1328: if ( !a ) *b = 0;
1329: else {
1330: NTOQ(BDY(a),1,q); *b = ztogz(q);
1331: }
1.5 noro 1332: }
1333:
1334: void setmod_lf(N p)
1335: {
1336: Q q;
1337:
1338: NTOQ(p,1,q); current_mod_lf = ztogz(q);
1.6 noro 1339: current_mod_lf_size = mpz_size(BDY(current_mod_lf))+1;
1.5 noro 1340: }
1341:
1342: void simplf_force(GZ a,GZ *b)
1343: {
1344: GZ t;
1345:
1346: remgz(a,current_mod_lf,&t); *b = t;
1347: }