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