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