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