Annotation of OpenXM_contrib2/asir2000/engine/pari-mp.c, Revision 1.3
1.1 saito 1: /*
1.3 ! ohara 2: * $OpenXM: OpenXM_contrib2/asir2000/engine/pari-mp.c,v 1.2 2000/12/26 06:12:36 noro Exp $
1.1 saito 3: */
4: /* for f-itv.c */
5:
1.3 ! ohara 6: #if defined(PARI)
1.2 noro 7:
1.1 saito 8: #include "genpari.h"
9: #include "itv-pari.h"
10:
11:
12: GEN
13: PariAddDown(GEN x, GEN y)
14: {
15: if(typ(x)==1) return (typ(y)==1) ? addii(x,y) : PariAddirDown(y,x);
16: return (typ(y)==1) ? PariAddirDown(y,x) : PariAddrrDown(x,y);
17: }
18:
19: GEN
20: PariAddirDown(GEN x, GEN y)
21: {
22: long l,e,ly,av,i,l1;
23: GEN z;
24:
25: if(!signe(x)) return rcopy(y);
26: if(!signe(y))
27: {
28: l=lgef(x)-(expo(y)>>TWOPOTBITS_IN_LONG);if((l<3)||(l>32767)) err(adder3);
29: z=cgetr(l);affir(x,z);return z;
30: }
31: else
32: {
33: e=expo(y)-expi(x);ly=lg(y);
34: if(e>0)
35: {
36: l=ly-(e>>TWOPOTBITS_IN_LONG);if(l<=2) return rcopy(y);
37: }
38: else
39: {
40: l=ly+((-e)>>TWOPOTBITS_IN_LONG)+1;if(l>32767) err(adder3);
41: }
42: av=avma;z=cgetr(l);affir(x,z);l1=av-avma;l=l1>>TWOPOTBYTES_IN_LONG;
43: z=PariAddrrDown(z,y);
44: for(i=lg(z)-1;i>=0;i--) z[i+l]=z[i];z+=l;avma+=l1;
45: }
46: return z;
47: }
48:
49: GEN
50: PariAddrrDown(GEN x, GEN y)
51: {
52: long sx=signe(x),sy=signe(y),lx=lg(x),ly=lg(y),lz,ex=expo(x),ey=expo(y),sz;
53: long av0=avma,e,l,i,d,m,flag,lp1,lp2,av,k,j,cex,f2;
54: GEN z,p1,p2;
55:
56: if(!sy)
57: {
58: if(!sx) {e=(ex>=ey)?ex:ey;z=cgetr(3);z[2]=0;z[1]=e+HIGHEXPOBIT;return z;}
59: e=ex-ey;
60: if(e<0) {z=cgetr(3);z[2]=0;z[1]=ey+HIGHEXPOBIT;return z;}
61: l=(e>>TWOPOTBITS_IN_LONG)+3;if(l>lx) l=lx;z=cgetr(l);
62: for(i=1;i<l;i++) z[i]=x[i];return z;
63: }
64: e=ey-ex;
65: if(!sx)
66: {
67: if(e<0) {z=cgetr(3);z[2]=0;z[1]=ex+HIGHEXPOBIT;return z;}
68: l=(e>>TWOPOTBITS_IN_LONG)+3;if(l>ly) l=ly;z=cgetr(l);
69: for(i=1;i<l;i++) z[i]=y[i];return z;
70: }
71: if(e)
72: {
73: if(e<0) {z=x;x=y;y=z;lz=lx;lx=ly;ly=lz;ey=ex;e= -e;sz=sx;sx=sy;sy=sz;}
74: d=(e>>TWOPOTBITS_IN_LONG);m=e&(BITS_IN_LONG-1);
75: if(d>=ly-2) return rcopy(y);
76: l=d+lx;
77: if(l>=ly)
78: {
79: flag=1;p1=cgetr(ly);lp1=ly;lp2=ly-d;
80: }
81: else
82: {
83: flag=0;p1=cgetr(l+1);lp2=lx+1;lp1=l+1;
84: }
85: av=avma;
86: if(m)
87: {
88: p2=cgetr(lp2);m=BITS_IN_LONG-m;
89: if(flag) {shiftl(x[lp2-1],m);k=hiremainder;}
90: else k=0;
91: for(i=lp2-1;i>=3;i--)
92: {
93: p2[i]=shiftl(x[i-1],m)+k;k=hiremainder;
94: }
95: p2[2]=k;
96: }
97: else p2=x;
98: }
99: else
100: {
101: l=(lx>ly)?ly:lx;p1=cgetr(l);av=avma;lp2=lp1=l;flag=2;p2=x;m=0;
102: }
103: if(sx==sy)
104: {
105: overflow=0;
106: if(m+flag) for(i=lp1-1,j=lp2-1;j>=2;i--,j--) p1[i]=addllx(p2[j],y[i]);
107: else
108: {
109: p1[lp1-1]=y[lp1-1];
110: for(i=lp1-2,j=lp2-2;j>=2;i--,j--) p1[i]=addllx(p2[j],y[i]);
111: }
112: if(overflow)
113: {
114: for(;(i>=2)&&(y[i]==(long)MAXULONG);i--) p1[i]=0;
115: if(i>=2) {cex=0;p1[i]=y[i]+1;while(i>=3) {i--;p1[i]=y[i];}}
116: else
117: {
118: cex=1;k=HIGHBIT;if(ey==(HIGHEXPOBIT-1)) err(adder4);
119: for(i=2;i<lp1;i++) {p1[i]=shiftlr(p1[i],1)+k;k=hiremainder;}
120:
121: if ( sx < 0 && hiremainder ) { /* |p1| + 1 */
122: overflow=1;
123: for(j=lg(p1);(overflow) && j>=2;j--) p1[j]=addllx(p1[j],0);
124: }
125:
126: }
127: }
128: else {cex=0;for(;i>=2;i--) p1[i]=y[i];}
129: p1[1]=evalsigne(sx)+ey+cex+HIGHEXPOBIT;
130: avma=av;return p1;
131: }
132: else
133: {
134: if(!e)
135: {
136: for(i=2;(i<l)&&(p2[i]==y[i]);i++);
137: if(i==l)
138: {
139: e=ex-((l-2)<<TWOPOTBITS_IN_LONG)+HIGHEXPOBIT;if(e<0) err(adder5);
140: if(e>EXPOBITS) err(adder4);
141: avma=av0;z=cgetr(3);z[2]=0;z[1]=e;return z;
142: }
143: else
144: {
145: f2=(((ulong)y[i])>((ulong)p2[i]))?1:0;
146: }
147: }
148: else f2=1;
149: if(f2)
150: {
151: overflow=0;
152: if(m+flag) for(i=lp1-1,j=lp2-1;j>=2;i--,j--) p1[i]=subllx(y[i],p2[j]);
153: else
154: {
155: p1[lp1-1]=y[lp1-1];
156: for(i=lp1-2,j=lp2-2;j>=2;i--,j--) p1[i]=subllx(y[i],p2[j]);
157: }
158: if(overflow)
159: {
160: for(;(i>=2)&&(!y[i]);i--) p1[i]=(long)MAXULONG;
161: p1[i]=y[i]-1;while(i>=3) {i--;p1[i]=y[i];}
162: }
163: else for(;i>=2;i--) p1[i]=y[i];
164: }
165: else
166: {
167: overflow=0;
168: if(m+flag) for(i=lp1-1;i>=2;i--) p1[i]=subllx(p2[i],y[i]);
169: else
170: {
171: p1[lp1-1]=subllx(0,y[lp1-1]);
172: for(i=lp1-2;i>=2;i--) p1[i]=subllx(p2[i],y[i]);
173: }
174: }
175: for(i=2;!p1[i];i++);j=i-2;avma=av+(j<<TWOPOTBYTES_IN_LONG);p1[j]=p1[0]-j;p1+=j;
176: m=bfffo(p1[2]);e=ey-(j<<TWOPOTBITS_IN_LONG)-m+HIGHEXPOBIT;
177: if(e<0) err(adder5);
178: p1[1]=f2 ? evalsigne(sy)+e : evalsigne(sx)+e;
179: if(m)
180: {
181: k=0;for(i=lp1-1-j;i>=2;i--) {p1[i]=shiftl(p1[i],m)+k;k=hiremainder;}
182: }
183: return p1;
184: }
185: }
186:
187: GEN
188: PariSubDown(GEN x, GEN y)
189: {
190: if(typ(x)==1) return (typ(y)==1) ? subii(x,y) : PariSubirDown(x,y);
191: return (typ(y)==1) ? PariSubriDown(x,y) : PariSubrrDown(x,y);
192: }
193:
194: GEN
195: PariSubrrDown(GEN x, GEN y)
196: {
197: long s=signe(y);
198: GEN z;
199:
200: if(x==y)
201: {
202: z=cgetr(3);z[2]=0;z[1]=HIGHEXPOBIT-(lg(x)<<TWOPOTBITS_IN_LONG);return z;
203: }
204: setsigne(y,-s);z=PariAddrrDown(x,y);setsigne(y,s);return z;
205: }
206:
207: GEN
208: PariSubirDown(GEN x, GEN y)
209: {
210: long s=signe(y);
211: GEN z;
212:
213: setsigne(y,-s);z=PariAddirDown(x,y);setsigne(y,s);return z;
214: }
215:
216: GEN
217: PariSubriDown(GEN x, GEN y)
218: {
219: long s=signe(y);
220: GEN z;
221:
222: setsigne(y,-s);z=PariAddirDown(y,x);setsigne(y,s);return z;
223: }
224:
225: GEN
226: PariMulDown(GEN x, GEN y)
227: {
228: if(typ(x)==1) return (typ(y)==1) ? mulii(x,y) : PariMulirDown(x,y);
229: return (typ(y)==1) ? PariMulirDown(y,x) : PariMulrrDown(x,y);
230: }
231:
232: GEN
233: PariMulrrDown(GEN x, GEN y)
234: {
235: long i,j,lx=lg(x),ly=lg(y),sx=signe(x),sy=signe(y),ex=expo(x),ey=expo(y);
236: long e,flag,garde,p1,p2,lz;
237: GEN z;
238:
239: e=ex+ey+HIGHEXPOBIT;if(e>=EXPOBITS) err(muler4);
240: if(e<0) err(muler5);
241: if((!sx)||(!sy)) {z=cgetr(3);z[2]=0;z[1]=e;return z;}
242: if(sy<0) sx= -sx;
243: if(lx>ly) {lz=ly;z=x;x=y;y=z;flag=1;} else {lz=lx;flag=(lx!=ly);}
244: z=cgetr(lz);z[1]=evalsigne(sx)+e;
245: if(flag) mulll(x[2],y[lz]);else hiremainder=0;
246: if(lz==3)
247: {
248: garde=flag ? addmul(x[2],y[2]) : mulll(x[2],y[2]);
249: if((long)hiremainder<0) {z[2]=hiremainder;z[1]++;}
250: else {z[2]=(garde<0)?(hiremainder<<1)+1:(hiremainder<<1);}
251: return z;
252: }
253: p1=x[lz-1];garde=hiremainder;
254: if(p1)
255: {
256: mulll(p1,y[3]);p2=addmul(p1,y[2]);
257: garde=addll(p2,garde);z[lz-1]=overflow+hiremainder;
258: }
259: else z[lz-1]=0;
260: for(j=lz-2;j>=3;j--)
261: {
262: p1=x[j];
263: if(p1)
264: {
265: mulll(p1,y[lz+2-j]);
266: p2=addmul(p1,y[lz+1-j]);
267: garde=addll(p2,garde);hiremainder+=overflow;
268: for(i=lz-j;i>=2;i--)
269: {
270: p2=addmul(p1,y[i]);
271: z[i+j-1]=addll(p2,z[i+j-1]);hiremainder+=overflow;
272: }
273: z[j]=hiremainder;
274: }
275: else z[j]=0;
276: }
277: p1=x[2];p2=mulll(p1,y[lz-1]);
278: garde=addll(p2,garde);hiremainder+=overflow;
279: for(i=lz-2;i>=2;i--)
280: {
281: p2=addmul(p1,y[i]);
282: z[i+1]=addll(p2,z[i+1]);hiremainder+=overflow;
283: }
284: z[2]=hiremainder;
285: if((long)hiremainder>0)
286: {
287: overflow=(garde<0)?1:0;
288: for(i=lz-1;i>=2;i--) {p1=z[i];z[i]=addllx(p1,p1);}
289: }
290: else z[1]++;
291: return z;
292: }
293:
294: GEN
295: PariMulirDown(GEN x, GEN y)
296: {
297: long sx=signe(x),sy,av,lz,ey,e,garde,p1,p2,i,j;
298: GEN z,temp;
299:
300: if(!sx) return gzero;
301: sy=signe(y);ey=expo(y);
302: if(!sy)
303: {
304: z=cgetr(3);z[2]=0;e=expi(x)+ey+HIGHEXPOBIT;if(e>EXPOBITS) err(muler6);
305: z[1]=e;return z;
306: }
307: lz=lg(y);if(sy<0) sx= -sx;
308: z=cgetr(lz);av=avma;
309: temp=cgetr(lz+1);affir(x,temp);x=y;y=temp;
310: e=expo(y)+ey+HIGHEXPOBIT;if(e>=EXPOBITS) err(muler4);
311: if(e<0) err(muler5);
312: z[1]=evalsigne(sx)+e;
313: mulll(x[2],y[lz]);
314: if(lz==3)
315: {
316: garde=addmul(x[2],y[2]);
317: if((long)hiremainder<0) {z[2]=hiremainder;z[1]++;}
318: else {z[2]=(garde<0)?(hiremainder<<1)+1:(hiremainder<<1);}
319: avma=av;return z;
320: }
321: garde=hiremainder;
322: p1=x[lz-1];mulll(p1,y[3]);p2=addmul(p1,y[2]);
323: garde=addll(p2,garde);z[lz-1]=overflow+hiremainder;
324: for(j=lz-2;j>=3;j--)
325: {
326: p1=x[j];mulll(p1,y[lz+2-j]);
327: p2=addmul(p1,y[lz+1-j]);
328: garde=addll(p2,garde);hiremainder+=overflow;
329: for(i=lz-j;i>=2;i--)
330: {
331: p2=addmul(p1,y[i]);
332: z[i+j-1]=addll(p2,z[i+j-1]);hiremainder+=overflow;
333: }
334: z[j]=hiremainder;
335: }
336: p1=x[2];p2=mulll(p1,y[lz-1]);
337: garde=addll(p2,garde);hiremainder+=overflow;
338: for(i=lz-2;i>=2;i--)
339: {
340: p2=addmul(p1,y[i]);
341: z[i+1]=addll(p2,z[i+1]);hiremainder+=overflow;
342: }
343: z[2]=hiremainder;
344: if((long)hiremainder>0)
345: {
346: overflow=(garde<0)?1:0;
347: for(i=lz-1;i>=2;i--) {p1=z[i];z[i]=addllx(p1,p1);}
348: }
349: else z[1]++;
350: avma=av;return z;
351: }
352:
1.3 ! ohara 353: #if defined(LONG_IS_32BIT)
1.1 saito 354: #define DIVCONVI 14
355: #endif
356:
1.3 ! ohara 357: #if defined(LONG_IS_64BIT)
1.1 saito 358: #define DIVCONVI 7
359: #endif
360:
361: GEN
362: PariDivDown(GEN x, GEN y)
363: {
364: if(typ(x)==1) return (typ(y)==1) ? divii(x,y) : PariDivirDown(x,y);
365: return (typ(y)==1) ? PariDivriDown(x,y) : PariDivrrDown(x,y);
366: }
367:
368: GEN
369: PariDivirDown(GEN x, GEN y)
370: {
371: GEN xr,z;
372: long av,ly;
373:
374: if(!signe(y)) err(diver5);
375: if(!signe(x)) return gzero;
376: ly=lg(y);z=cgetr(ly);av=avma;affir(x,xr=cgetr(ly+1));
377: xr=PariDivrrDown(xr,y);affrr(xr,z);avma=av;return z;
378: }
379:
380: GEN
381: PariDivriDown(GEN x, GEN y)
382: {
383: GEN yr,z;
384: long av,lx,ex,s=signe(y);
385:
386: if(!s) err(diver8);
387: if(!signe(x))
388: {
389: ex=expo(x)-expi(y)+HIGHEXPOBIT;
390: if(ex<0) err(diver12);
391: z=cgetr(3);z[1]=ex;z[2]=0;return z;
392: }
393: if((lg(y)==3)&&(y[2]>0)) return (s>0) ? divrs(x,y[2]) : divrs(x,-y[2]);
394: lx=lg(x);z=cgetr(lx);av=avma;affir(y,yr=cgetr(lx+1));
395: yr=PariDivrrDown(x,yr);affrr(yr,z);avma=av;return z;
396: }
397:
398: GEN
399: PariDivrrDown(GEN x, GEN y)
400: {
401: long sx=signe(x),sy=signe(y),lx,ly,lz,ex,ex1,i,z0;
402: ulong ldif,y0,y1,si,saux,qp,k,k3,k4,j;
403: GEN z;
404: if(!sy) err(diver9);
405: ex=expo(x)-expo(y)+HIGHEXPOBIT;
406: if(ex<=0) err(diver10);
407: if(ex>EXPOBITS) err(diver11);
408: if(!sx)
409: {
410: z=cgetr(3);z[1]=ex;z[2]=0;return z;
411: }
412: lx=lg(x);ly=lg(y);lz=(lx<=ly)?lx:ly;
413: z=cgetr(lz);if(sy<0) sx= -sx;
414: ex1=evalsigne(sx)+ex;
415: if(ly==3)
416: {
417: i=x[2];si=(lx>3)?x[3]:0;
418: if((ulong)i<(ulong)y[2])
419: {
420: hiremainder=i;z[2]=divll(si,y[2]);
421: z[1]=ex1-1;return z;
422: }
423: else
424: {
425: hiremainder=((ulong)i)>>1;
426: z[2]=(i&1)?divll((((ulong)si)>>1)|(HIGHBIT),y[2]):divll(((ulong)si)>>1,y[2]);
427: z[1]=ex1;return z;
428: }
429: }
430: z0= *z;*z=0;
431: for(i=2;i<=lz-1;i++) z[i-1]=x[i];
432: z[lz-1]=(lx>lz) ? x[lz] : 0;
433: ldif=ly-lz;if(!ldif) {y0=y[lz];y[lz]=0;}
434: if(ldif<=1) {y1=y[lz+1];y[lz+1]=0;}
435: si=y[2];saux=y[3];
436: for(i=0;i<lz-1;i++)
437: {
438: if(z[i]!=si)
439: {
440: if((ulong)z[i]>si)
441: {
442: overflow=0;
443: for(j=lz-i+1;j>=2;j--) z[i+j-2]=subllx(z[i+j-2],y[j]);
444: {z[i-1]++;for(j=i-1;j&&(!z[j]);j--) z[j-1]++;}
445: }
446: hiremainder=z[i];qp=divll(z[i+1],si);
447: overflow=0;k=hiremainder;
448: }
449: else
450: {
451: qp=(long)MAXULONG;k=addll(si,z[i+1]);
452: }
453: if(!overflow)
454: {
455: k3=subll(mulll(qp,saux),z[i+2]);k4=subllx(hiremainder,k);
456: while((!overflow)&&k4) {qp--;k3=subll(k3,saux);k4=subllx(k4,si);}
457: }
458: mulll(qp,y[lz+1-i]);
459: for(j=lz-i;j>=2;j--)
460: {
461: z[i+j-1]=subll(z[i+j-1],addmul(qp,y[j]));hiremainder+=overflow;
462: }
463: if((ulong)z[i]!=(ulong)hiremainder)
464: {
465: if((ulong)z[i]<(ulong)hiremainder)
466: {
467: overflow=0;qp--;
468: for(j=lz-i;j>=2;j--) z[i+j-1]=addllx(z[i+j-1],y[j]);
469: }
470: else
471: {
472: z[i]-=hiremainder;
473: while(z[i])
474: {
475: overflow=0;qp++;
476: if(!qp) {z[i-1]++;for(j=i-1;j&&(!z[j]);j--) z[j-1]++;}
477: for(j=lz-i;j>=2;j--) z[i+j-1]=subllx(z[i+j-1],y[j]);
478: z[i]-=overflow;
479: }
480: }
481: }
482: z[i]=qp;
483: }
484: if(!ldif) y[lz]=y0;if(ldif<=1) y[lz+1]=y1;
485: for(j=lz-1;j>=2;j--) z[j]=z[j-1];
486: if(*z)
487: {
488: k=HIGHBIT;
489: for(j=2;j<lz;j++) {z[j]=shiftlr(z[j],1)+k;k=hiremainder;}
490: }
491: else ex1--;
492: z[1]=ex1;*z=z0;return z;
493: }
494:
1.2 noro 495: #endif /* PARI */
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to pari-mp.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / engine