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