Annotation of OpenXM_contrib/pari/src/language/sumiter.c, Revision 1.1.1.1
1.1 maekawa 1: /********************************************************************/
2: /** **/
3: /** BIBLIOTHEQUE MATHEMATIQUE **/
4: /** (sommes, produits, iterations) **/
5: /** **/
6: /********************************************************************/
7: /* $Id: sumiter.c,v 1.3 1999/09/23 17:50:57 karim Exp $ */
8: #include "pari.h"
9: #include "anal.h"
10:
11: /********************************************************************/
12: /** **/
13: /** ITERATIONS **/
14: /** **/
15: /********************************************************************/
16:
17: void
18: forpari(entree *ep, GEN a, GEN b, char *ch)
19: {
20: long av,av0 = avma, lim;
21:
22: b = gcopy(b); av=avma; lim = stack_lim(av,1);
23: /* gcopy nedeed in case b gets overwritten in ch, as in
24: * b=10; for(a=1,b, print(a);b=1)
25: */
26: push_val(ep, a);
27: while (gcmp(a,b) <= 0)
28: {
29: long av1=avma; (void)lisseq(ch); avma=av1;
30: if (loop_break()) break;
31: a = (GEN) ep->value; a = gadd(a,gun);
32: if (low_stack(lim, stack_lim(av,1)))
33: {
34: if (DEBUGMEM>1) err(warnmem,"forpari");
35: a = gerepileupto(av,a);
36: }
37: ep->value = (void*)a;
38: }
39: pop_val(ep); avma = av0;
40: }
41:
42: static int negcmp(GEN x, GEN y) { return gcmp(y,x); }
43:
44: void
45: forstep(entree *ep, GEN a, GEN b, GEN s, char *ch)
46: {
47: long ss, av,av0 = avma, lim, i, lgv;
48: GEN v = NULL;
49: int (*cmp)(GEN,GEN);
50:
51: b = gcopy(b); av=avma; lim = stack_lim(av,1);
52: push_val(ep, a);
53: if (is_vec_t(typ(s)))
54: {
55: v=s; lgv=lg(v)-1; s=gzero;
56: for (i=1; i<=lgv; i++) s = gadd(s,(GEN)v[i]);
57: }
58: ss = gsigne(s);
59: if (!ss) err(talker, "step equal to zero in forstep");
60: cmp = (ss > 0)? gcmp: negcmp;
61: i = 0;
62: while (cmp(a,b) <= 0)
63: {
64: long av1=avma; (void)lisseq(ch); avma=av1;
65: if (loop_break()) break;
66: if (v) { if (++i > lgv) i=1; s = (GEN)v[i]; }
67: a = (GEN) ep->value; a = gadd(a,s);
68:
69: if (low_stack(lim, stack_lim(av,1)))
70: {
71: if (DEBUGMEM>1) err(warnmem,"forstep");
72: a = gerepileupto(av,a);
73: }
74: ep->value = (void*)a;
75: }
76: pop_val(ep); avma = av0;
77: }
78:
79: void
80: forprime(entree *ep, GEN ga, GEN gb, char *ch)
81: {
82: long av = avma, a,b;
83: long prime[] = {evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3), 0};
84: byteptr p=diffptr;
85:
86: ga = gceil(ga); gb = gfloor(gb);
87: if (is_bigint(ga) || is_bigint(gb)) err(primer1);
88: a = itos(ga); if (a<=0) a = 1;
89: b = 0; while (*p && a > b) b += *p++;
90: prime[2] = b; b = itos(gb);
91: avma = av; push_val(ep, prime);
92: while (prime[2] <= b)
93: {
94: if (!*p) err(primer1);
95: (void)lisseq(ch); if (loop_break()) break;
96: avma = av; prime[2] += *p++;
97: }
98: pop_val(ep);
99: }
100:
101: void
102: fordiv(GEN a, entree *ep, char *ch)
103: {
104: long i,av2,l, av = avma;
105: GEN t = divisors(a);
106:
107: push_val(ep, NULL); l=lg(t); av2 = avma;
108: for (i=1; i<l; i++)
109: {
110: ep->value = (void*) t[i];
111: (void)lisseq(ch); if (loop_break()) break;
112: avma = av2;
113: }
114: pop_val(ep); avma=av;
115: }
116:
117: /* Embedded for loops:
118: * fl = 0: execute ch (a), where a = (ai) runs through all n-uplets in
119: * [m1,M1] x ... x [mn,Mn]
120: * fl = 1: impose a1 <= ... <= an
121: * fl = 2: a1 < ... < an
122: */
123: static GEN *fv_a, *fv_m, *fv_M;
124: static long fv_n, fv_fl;
125: static char *fv_ch;
126:
127: /* general case */
128: static void
129: fvloop(long i)
130: {
131: fv_a[i] = fv_m[i];
132: if (fv_fl && i > 1)
133: {
134: GEN p1 = gsub(fv_a[i], fv_a[i-1]);
135: if (gsigne(p1) < 0)
136: fv_a[i] = gadd(fv_a[i], gceil(gneg_i(p1)));
137: if (fv_fl == 2 && gegal(fv_a[i], fv_a[i-1]))
138: fv_a[i] = gadd(fv_a[i], gun);
139: }
140: if (i+1 == fv_n)
141: while (gcmp(fv_a[i], fv_M[i]) <= 0)
142: {
143: long av = avma; (void)lisseq(fv_ch); avma = av;
144: if (loop_break()) { fv_n = 0; return; }
145: fv_a[i] = gadd(fv_a[i], gun);
146: }
147: else
148: while (gcmp(fv_a[i], fv_M[i]) <= 0)
149: {
150: long av = avma; fvloop(i+1); avma = av;
151: if (!fv_n) return;
152: fv_a[i] = gadd(fv_a[i], gun);
153: }
154: }
155:
156: /* we run through integers */
157: static void
158: fvloop_i(long i)
159: {
160: fv_a[i] = setloop(fv_m[i]);
161: if (fv_fl && i > 1)
162: {
163: int c = cmpii(fv_a[i], fv_a[i-1]);
164: if (c < 0)
165: {
166: fv_a[i] = setloop(fv_a[i-1]);
167: c = 0;
168: }
169: if (c == 0 && fv_fl == 2)
170: fv_a[i] = incloop(fv_a[i]);
171: }
172: if (i+1 == fv_n)
173: while (gcmp(fv_a[i], fv_M[i]) <= 0)
174: {
175: long av = avma; (void)lisseq(fv_ch); avma = av;
176: if (loop_break()) { fv_n = 0; return; }
177: fv_a[i] = incloop(fv_a[i]);
178: }
179: else
180: while (gcmp(fv_a[i], fv_M[i]) <= 0)
181: {
182: long av = avma; fvloop_i(i+1); avma = av;
183: if (!fv_n) return;
184: fv_a[i] = incloop(fv_a[i]);
185: }
186: }
187:
188: void
189: forvec(entree *ep, GEN x, char *c, long flag)
190: {
191: long i, av = avma, tx = typ(x), n = fv_n, fl = fv_fl;
192: GEN *a = fv_a, *M = fv_M, *m = fv_m;
193: char *ch = fv_ch;
194: void (*fv_fun)(long) = fvloop_i;
195:
196: if (!is_vec_t(tx)) err(talker,"not a vector in forvec");
197: if (fl<0 || fl>2) err(flagerr);
198: fv_n = lg(x);
199: fv_ch = c;
200: fv_fl = flag;
201: fv_a = (GEN*)cgetg(fv_n,t_VEC); push_val(ep, (GEN)fv_a);
202: fv_m = (GEN*)cgetg(fv_n,t_VEC);
203: fv_M = (GEN*)cgetg(fv_n,t_VEC);
204: if (fv_n == 1) lisseq(fv_ch);
205: else
206: {
207: for (i=1; i<fv_n; i++)
208: {
209: GEN *c = (GEN*) x[i];
210: tx = typ(c);
211: if (! is_vec_t(tx) || lg(c)!=3)
212: err(talker,"not a vector of two-component vectors in forvec");
213: if (gcmp(c[1],c[2]) > 0) fv_n = 0;
214: /* in case x is an ep->value and c destroys it, we have to copy */
215: if (typ(c[1]) != t_INT) fv_fun = fvloop;
216: fv_m[i] = gcopy(c[1]);
217: fv_M[i] = gcopy(c[2]);
218: }
219: fv_fun(1);
220: }
221: pop_val(ep);
222: fv_n = n;
223: fv_ch = ch;
224: fv_fl = fl;
225: fv_a = a;
226: fv_m = m;
227: fv_M = M; avma = av;
228: }
229:
230: /********************************************************************/
231: /** **/
232: /** SUMS **/
233: /** **/
234: /********************************************************************/
235:
236: GEN
237: somme(entree *ep, GEN a, GEN b, char *ch, GEN x)
238: {
239: long av,av0 = avma, lim;
240: GEN p1;
241:
242: if (typ(a) != t_INT) err(talker,"non integral index in sum");
243: if (!x) x = gzero;
244: if (gcmp(b,a) < 0) return gcopy(x);
245:
246: b = gfloor(b);
247: a = setloop(a);
248: av=avma; lim = stack_lim(av,1);
249: push_val(ep, a);
250: for(;;)
251: {
252: p1 = lisexpr(ch); if (did_break()) err(breaker,"sum");
253: x=gadd(x,p1); if (cmpii(a,b) >= 0) break;
254: a = incloop(a);
255: if (low_stack(lim, stack_lim(av,1)))
256: {
257: if (DEBUGMEM>1) err(warnmem,"sum");
258: x = gerepileupto(av,x);
259: }
260: ep->value = (void*) a;
261: }
262: pop_val(ep); return gerepileupto(av0,x);
263: }
264:
265: GEN
266: suminf(entree *ep, GEN a, char *ch, long prec)
267: {
268: long fl,G,tetpil, av0 = avma, av,lim;
269: GEN p1,x = realun(prec);
270:
271: if (typ(a) != t_INT) err(talker,"non integral index in suminf");
272: a = setloop(a);
273: av = avma; lim = stack_lim(av,1);
274: push_val(ep, a);
275: fl=0; G = bit_accuracy(prec) + 5;
276: for(;;)
277: {
278: p1 = lisexpr(ch); if (did_break()) err(breaker,"suminf");
279: x = gadd(x,p1); a = incloop(a);
280: if (gcmp0(p1) || gexpo(p1) <= gexpo(x)-G)
281: { if (++fl==3) break; }
282: else
283: fl=0;
284: if (low_stack(lim, stack_lim(av,1)))
285: {
286: if (DEBUGMEM>1) err(warnmem,"suminf");
287: x = gerepileupto(av,x);
288: }
289: ep->value = (void*)a;
290: }
291: pop_val(ep); tetpil=avma;
292: return gerepile(av0,tetpil,gsub(x,gun));
293: }
294:
295: GEN
296: divsum(GEN num, entree *ep, char *ch)
297: {
298: long av=avma,d,n,d2;
299: GEN y,z, p1 = icopy(gun);
300:
301: push_val(ep, p1); n=itos(num); /* provisoire */
302: y=gzero;
303: for (d=d2=1; d2 < n; d++, d2 += d+d-1)
304: if (n%d == 0)
305: {
306: p1[2]=d; y=gadd(y, lisexpr(ch));
307: if (did_break()) err(breaker,"divsum");
308: p1[2]=n/d; z = lisexpr(ch);
309: if (did_break()) err(breaker,"divsum");
310: y=gadd(y,z);
311: }
312: if (d2 == n)
313: {
314: p1[2]=d; z = lisexpr(ch);
315: if (did_break()) err(breaker,"divsum");
316: y=gadd(y,z);
317: }
318: pop_val(ep); return gerepileupto(av,y);
319: }
320:
321: /********************************************************************/
322: /** **/
323: /** PRODUCTS **/
324: /** **/
325: /********************************************************************/
326:
327: GEN
328: produit(entree *ep, GEN a, GEN b, char *ch, GEN x)
329: {
330: long av,av0 = avma, lim;
331: GEN p1;
332:
333: if (typ(a) != t_INT) err(talker,"non integral index in sum");
334: if (!x) x = gun;
335: if (gcmp(b,a) < 0) return gcopy(x);
336:
337: b = gfloor(b);
338: a = setloop(a);
339: av=avma; lim = stack_lim(av,1);
340: push_val(ep, a);
341: for(;;)
342: {
343: p1 = lisexpr(ch); if (did_break()) err(breaker,"prod");
344: x = gmul(x,p1); if (cmpii(a,b) >= 0) break;
345: a = incloop(a);
346: if (low_stack(lim, stack_lim(av,1)))
347: {
348: if (DEBUGMEM>1) err(warnmem,"prod");
349: x = gerepileupto(av,x);
350: }
351: ep->value = (void*) a;
352: }
353: pop_val(ep); return gerepileupto(av0,x);
354: }
355:
356: GEN
357: prodinf0(entree *ep, GEN a, char *ch, long flag, long prec)
358: {
359: switch(flag)
360: {
361: case 0: return prodinf(ep,a,ch,prec);
362: case 1: return prodinf1(ep,a,ch,prec);
363: }
364: err(flagerr);
365: return NULL; /* not reached */
366: }
367:
368: GEN
369: prodinf(entree *ep, GEN a, char *ch, long prec)
370: {
371: long fl,G,tetpil, av0 = avma, av,lim;
372: GEN p1,x = realun(prec);
373:
374: if (typ(a) != t_INT) err(talker,"non integral index in prodinf");
375: a = setloop(a);
376: av = avma; lim = stack_lim(av,1);
377: push_val(ep, a);
378: fl=0; G = -bit_accuracy(prec)-5;
379: for(;;)
380: {
381: p1 = lisexpr(ch); if (did_break()) err(breaker,"prodinf");
382: x=gmul(x,p1); a = incloop(a);
383: if (gexpo(gsubgs(p1,1)) <= G) { if (++fl==3) break; } else fl=0;
384: if (low_stack(lim, stack_lim(av,1)))
385: {
386: if (DEBUGMEM>1) err(warnmem,"prodinf");
387: x = gerepileupto(av,x);
388: }
389: ep->value = (void*)a;
390: }
391: pop_val(ep); tetpil=avma;
392: return gerepile(av0,tetpil,gcopy(x));
393: }
394:
395: GEN
396: prodinf1(entree *ep, GEN a, char *ch, long prec)
397: {
398: long fl,G,tetpil, av0 = avma, av,lim;
399: GEN p1,p2,x = realun(prec);
400:
401: if (typ(a) != t_INT) err(talker,"non integral index in prodinf1");
402: a = setloop(a);
403: av = avma; lim = stack_lim(av,1);
404: push_val(ep, a);
405: fl=0; G = -bit_accuracy(prec)-5;
406: for(;;)
407: {
408: p2 = lisexpr(ch); if (did_break()) err(breaker,"prodinf1");
409: p1 = gadd(p2,gun); x=gmul(x,p1); a = incloop(a);
410: if (gcmp0(p1) || gexpo(p2) <= G) { if (++fl==3) break; } else fl=0;
411: if (low_stack(lim, stack_lim(av,1)))
412: {
413: if (DEBUGMEM>1) err(warnmem,"prodinf1");
414: x = gerepileupto(av,x);
415: }
416: ep->value = (void*)a;
417: }
418: pop_val(ep); tetpil=avma;
419: return gerepile(av0,tetpil,gcopy(x));
420: }
421:
422: GEN
423: prodeuler(entree *ep, GEN a, GEN b, char *ch, long prec)
424: {
425: long prime,tetpil, av,av0 = avma, lim;
426: GEN p1,x = realun(prec);
427: byteptr p=diffptr;
428:
429: prime = 0;
430: b = gfloor(b); a = gceil(a);
431: while (*p && cmpis(a,prime)>0) prime += *p++;
432: if (cmpsi(prime,b) > 0)
433: {
434: av=avma; if (gcmp1(subii(a,b))) { avma=av; return x; }
435: err(talker,"incorrect indices in prodeuler");
436: }
437: av=avma; lim = stack_lim(avma,1);
438: a = stoi(prime); push_val(ep, a);
439: do
440: {
441: if (!*p) err(primer1);
442: p1 = lisexpr(ch); if (did_break()) err(breaker,"prodeuler");
443: x=gmul(x,p1); a = addsi(*p++,a);
444: if (low_stack(lim, stack_lim(av,1)))
445: {
446: GEN *gptr[2]; gptr[0]=&x; gptr[1]=&a;
447: if (DEBUGMEM>1) err(warnmem,"prodeuler");
448: gerepilemany(av,gptr,2);
449: }
450: ep->value = (void*)a;
451: }
452: while (gcmp(a,b)<=0);
453: pop_val(ep); tetpil=avma;
454: return gerepile(av0,tetpil,gcopy(x));
455: }
456:
457: GEN
458: direuler(entree *ep, GEN a, GEN b, char *ch)
459: {
460: GEN p1,x,x1,s,polnum,polden,c0;
461: long av0 = avma,av,tetpil,lim = (av0+bot)>>1, prime = 0,q,n,i,j,k,k1,tx,lx;
462: byteptr p = diffptr;
463:
464: if (typ(a) != t_INT) err(talker,"non integral index in direuler");
465: if (gcmpgs(b,2)<0) { x=cgetg(2,t_VEC); x[1]=un; return x; }
466: if (gcmpgs(a,2) < 0) a = gdeux;
467: if (gcmpgs(b, MAXHALFULONG-1) > 0) b = stoi(MAXHALFULONG-1);
468: n = itos(b);
469:
470: x1=cgetg(n+1,t_VEC); b = gcopy(b); av=avma;
471: x=cgetg(n+1,t_VEC); x[1]=un; for (i=2; i<=n; i++) x[i]=zero;
472:
473: while (*p && gcmpgs(a,prime) > 0) prime += *p++;
474: a = stoi(prime); push_val(ep, a);
475: while (gcmp(a,b)<=0)
476: {
477: if (!*p) err(primer1);
478: p1 = lisexpr(ch); if (did_break()) err(breaker,"direuler");
479: polnum=numer(p1); polden=denom(p1);
480: tx = typ(polnum);
481: if (is_scalar_t(tx))
482: {
483: if (!gcmp1(polnum))
484: err(talker,"constant term not equal to 1 in direuler");
485: }
486: else
487: {
488: if (tx != t_POL) err(typeer,"direuler");
489: c0 = truecoeff(polnum,0);
490: if (!gcmp1(c0)) err(talker,"constant term not equal to 1 in direuler");
491: for (i=1; i<=n; i++) x1[i]=x[i];
492: prime=itos(a); q=prime; j=1; lx=lgef(polnum)-3;
493: while (q<=n && j<=lx)
494: {
495: c0=(GEN)polnum[j+2];
496: if (!gcmp0(c0))
497: for (k=1,k1=q; k1<=n; k1+=q,k++)
498: x[k1] = ladd((GEN)x[k1], gmul(c0,(GEN)x1[k]));
499: q*=prime; j++;
500: }
501: }
502: tx=typ(polden);
503: if (is_scalar_t(tx))
504: {
505: if (!gcmp1(polden))
506: err(talker,"constant term not equal to 1 in direuler");
507: }
508: else
509: {
510: if (tx != t_POL) err(typeer,"direuler");
511: c0 = truecoeff(polden,0);
512: if (!gcmp1(c0)) err(talker,"constant term not equal to 1 in direuler");
513: prime=itos(a); lx=lgef(polden)-3;
514: for (i=prime; i<=n; i+=prime)
515: {
516: s=gzero; k=i; j=1;
517: while (!(k%prime) && j<=lx)
518: {
519: c0=(GEN)polden[j+2]; k/=prime; j++;
520: if (!gcmp0(c0)) s=gadd(s,gmul(c0,(GEN)x[k]));
521: }
522: x[i]=lsub((GEN)x[i],s);
523: }
524: }
525: if (low_stack(lim, stack_lim(av,1)))
526: {
527: GEN *gptr[2]; gptr[0]=&x; gptr[1]=&a;
528: if (DEBUGMEM>1) err(warnmem,"direuler");
529: gerepilemany(av,gptr,2);
530: }
531: a = addsi(*p++,a); ep->value = (void*) a;
532: }
533: pop_val(ep); tetpil=avma;
534: return gerepile(av0,tetpil,gcopy(x));
535: }
536:
537: /********************************************************************/
538: /** **/
539: /** VECTORS & MATRICES **/
540: /** **/
541: /********************************************************************/
542:
543: GEN
544: vecteur(GEN nmax, entree *ep, char *ch)
545: {
546: GEN y,p1;
547: long i,m;
548: long c[]={evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3), 0};
549:
550: if (typ(nmax)!=t_INT || signe(nmax) < 0)
551: err(talker,"bad number of components in vector");
552: m=itos(nmax); y=cgetg(m+1,t_VEC);
553: if (!ep || !ch)
554: {
555: for (i=1; i<=m; i++) y[i] = zero;
556: return y;
557: }
558: push_val(ep, c);
559: for (i=1; i<=m; i++)
560: {
561: c[2]=i; p1 = lisseq(ch);
562: if (did_break()) err(breaker,"vector");
563: y[i] = isonstack(p1)? (long)p1 : (long)forcecopy(p1);
564: }
565: pop_val(ep); return y;
566: }
567:
568: GEN
569: vvecteur(GEN nmax, entree *ep, char *ch)
570: {
571: GEN y = vecteur(nmax,ep,ch);
572: settyp(y,t_COL); return y;
573: }
574:
575: GEN
576: matrice(GEN nlig, GEN ncol,entree *ep1, entree *ep2, char *ch)
577: {
578: GEN y,z,p1;
579: long i,j,m,n,s;
580: long c1[]={evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3), 1};
581: long c2[]={evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3), 1};
582:
583: s=signe(ncol);
584: if (typ(ncol)!=t_INT || s<0) err(talker,"bad number of columns in matrix");
585: if (!s) return cgetg(1,t_MAT);
586:
587: s=signe(nlig);
588: if (typ(nlig)!=t_INT || s<0) err(talker,"bad number of rows in matrix");
589: m=itos(ncol)+1;
590: n=itos(nlig)+1; y=cgetg(m,t_MAT);
591: if (!s)
592: {
593: for (i=1; i<m; i++) y[i]=lgetg(1,t_COL);
594: return y;
595: }
596: if (!ep1 || !ep2 || !ch)
597: {
598: for (i=1; i<m; i++)
599: {
600: z=cgetg(n,t_COL); y[i]=(long)z;
601: for (j=1; j<n; j++) z[j] = zero;
602: }
603: return y;
604: }
605: push_val(ep1, c1);
606: push_val(ep2, c2);
607: for (i=1; i<m; i++)
608: {
609: c2[2]=i; z=cgetg(n,t_COL); y[i]=(long)z;
610: for (j=1; j<n; j++)
611: {
612: c1[2]=j; p1=lisseq(ch);
613: if (did_break()) err(breaker,"matrix");
614: z[j] = isonstack(p1)? (long)p1 : (long)forcecopy(p1);
615: }
616: }
617: pop_val(ep2);
618: pop_val(ep1); return y;
619: }
620:
621: /********************************************************************/
622: /** **/
623: /** SOMMATION DE SERIES **/
624: /** **/
625: /********************************************************************/
626:
627: GEN
628: sumalt0(entree *ep, GEN a, char *ch, long flag, long prec)
629: {
630: switch(flag)
631: {
632: case 0: return sumalt(ep,a,ch,prec);
633: case 1: return sumalt2(ep,a,ch,prec);
634: default: err(flagerr);
635: }
636: return NULL; /* not reached */
637: }
638:
639: GEN
640: sumpos0(entree *ep, GEN a, char *ch, long flag, long prec)
641: {
642: switch(flag)
643: {
644: case 0: return sumpos(ep,a,ch,prec);
645: case 1: return sumpos2(ep,a,ch,prec);
646: default: err(flagerr);
647: }
648: return NULL; /* not reached */
649: }
650:
651: GEN
652: sumalt(entree *ep, GEN a, char *ch, long prec)
653: {
654: long av = avma, tetpil,k,N;
655: GEN s,az,c,x,e1,d;
656:
657: if (typ(a) != t_INT) err(talker,"non integral index in sumalt");
658:
659: push_val(ep, a);
660: e1=addsr(3,gsqrt(stoi(8),prec));
661: N=(long)(0.4*(bit_accuracy(prec) + 7));
662: d=gpuigs(e1,N); d=shiftr(addrr(d,divsr(1,d)),-1);
663: az=negi(gun); c=d; s=gzero;
664: for (k=0; ; k++)
665: {
666: x = lisexpr(ch); if (did_break()) err(breaker,"sumalt");
667: c = gadd(az,c); s = gadd(s,gmul(x,c));
668: az = divii(mulii(mulss(N-k,N+k),shifti(az,1)),mulss(k+1,k+k+1));
669: if (k==N-1) break;
670: a = addsi(1,a); ep->value = (void*) a;
671: }
672: tetpil=avma; pop_val(ep);
673: return gerepile(av,tetpil,gdiv(s,d));
674: }
675:
676: GEN
677: sumalt2(entree *ep, GEN a, char *ch, long prec)
678: {
679: long av = avma, tetpil,k,N;
680: GEN x,s,dn,pol;
681:
682: if (typ(a) != t_INT) err(talker,"non integral index in sumalt");
683:
684: push_val(ep, a);
685: N=(long)(0.31*(bit_accuracy(prec) + 5));
686: s=gzero; pol=polzagreel(N,N>>1,prec+1); dn=poleval(pol,gun);
687: pol[2]=lsub((GEN)pol[2],dn); pol=gdiv(pol,gsub(polx[0],gun));
688: N=lgef(pol)-2;
689: for (k=0; k<N; k++)
690: {
691: x = lisexpr(ch); if (did_break()) err(breaker,"sumalt2");
692: s=gadd(s,gmul(x,(GEN)pol[k+2]));
693: if (k==N-1) break;
694: a = addsi(1,a); ep->value = (void*) a;
695: }
696: tetpil=avma; pop_val(ep);
697: return gerepile(av,tetpil,gdiv(s,dn));
698: }
699:
700: GEN
701: sumpos(entree *ep, GEN a, char *ch, long prec)
702: {
703: long av = avma,tetpil,k,kk,N,G;
704: GEN p1,r,q1,reel,s,az,c,x,e1,d, *stock;
705:
706: if (typ(a) != t_INT) err(talker,"non integral index in sumpos");
707:
708: push_val(ep, NULL);
709: a=subis(a,1); reel=cgetr(prec);
710: e1=addsr(3,gsqrt(stoi(8),prec));
711: N=(long)(0.4*(bit_accuracy(prec) + 7));
712: d=gpuigs(e1,N); d=shiftr(addrr(d,divsr(1,d)),-1);
713: az=negi(gun); c=d; s=gzero;
714: G = -bit_accuracy(prec) - 5;
715: stock=(GEN*)new_chunk(N+1); for (k=1; k<=N; k++) stock[k]=NULL;
716: for (k=0; k<N; k++)
717: {
718: if (k&1 && stock[k]) x=stock[k];
719: else
720: {
721: x=gzero; r=stoi(2*k+2);
722: for(kk=0;;kk++)
723: {
724: long ex;
725: q1 = addii(r,a); ep->value = (void*) q1;
726: p1 = lisexpr(ch); if (did_break()) err(breaker,"sumpos");
727: gaffect(p1,reel); ex = expo(reel)+kk; setexpo(reel,ex);
728: x = gadd(x,reel); if (kk && ex < G) break;
729: r=shifti(r,1);
730: }
731: if (2*k<N) stock[2*k+1]=x;
732: q1 = addsi(k+1,a); ep->value = (void*) q1;
733: p1 = lisexpr(ch); if (did_break()) err(breaker,"sumpos");
734: gaffect(p1,reel); x = gadd(reel, gmul2n(x,1));
735: }
736: c=gadd(az,c); p1 = k&1? gneg_i(c): c;
737: s = gadd(s,gmul(x,p1));
738: az = divii(mulii(mulss(N-k,N+k),shifti(az,1)),mulss(k+1,k+k+1));
739: }
740: tetpil=avma; pop_val(ep);
741: return gerepile(av,tetpil,gdiv(s,d));
742: }
743:
744: GEN
745: sumpos2(entree *ep, GEN a, char *ch, long prec)
746: {
747: long av = avma,tetpil,k,kk,N,G;
748: GEN p1,r,q1,reel,s,pol,dn,x, *stock;
749:
750: if (typ(a) != t_INT) err(talker,"non integral index in sumpos2");
751:
752: push_val(ep, a);
753: a=subis(a,1); reel=cgetr(prec);
754: N=(long)(0.31*(bit_accuracy(prec) + 5));
755: G = -bit_accuracy(prec) - 5;
756: stock=(GEN*)new_chunk(N+1); for (k=1; k<=N; k++) stock[k]=NULL;
757: for (k=1; k<=N; k++)
758: {
759: if (odd(k) || !stock[k])
760: {
761: x=gzero; r=stoi(2*k);
762: for(kk=0;;kk++)
763: {
764: long ex;
765: q1 = addii(r,a); ep->value = (void*) q1;
766: p1 = lisexpr(ch); if (did_break()) err(breaker,"sumpos2");
767: gaffect(p1,reel); ex = expo(reel)+kk; setexpo(reel,ex);
768: x=gadd(x,reel); if (kk && ex < G) break;
769: r=shifti(r,1);
770: }
771: if (2*k-1<N) stock[2*k]=x;
772: q1 = addsi(k,a); ep->value = (void*) q1;
773: p1 = lisexpr(ch); if (did_break()) err(breaker,"sumpos2");
774: gaffect(p1,reel); stock[k] = gadd(reel, gmul2n(x,1));
775: }
776: }
777: pop_val(ep); s = gzero;
778: pol=polzagreel(N,N>>1,prec+1); dn=poleval(pol,gun);
779: pol[2]=lsub((GEN)pol[2],dn); pol=gdiv(pol,gsub(gun,polx[0]));
780: for (k=1; k<=lgef(pol)-2; k++)
781: {
782: p1 = gmul((GEN)pol[k+1],stock[k]);
783: if (k&1) p1 = gneg_i(p1);
784: s = gadd(s,p1);
785: }
786: tetpil=avma; return gerepile(av,tetpil,gdiv(s,dn));
787: }
788:
789: /********************************************************************/
790: /** **/
791: /** NUMERICAL INTEGRATION (Romberg) **/
792: /** **/
793: /********************************************************************/
794: GEN
795: intnum0(entree *ep, GEN a, GEN b, char *ch, long flag, long prec)
796: {
797: switch(flag)
798: {
799: case 0: return qromb (ep,a,b,ch,prec);
800: case 1: return rombint(ep,a,b,ch,prec);
801: case 2: return qromi (ep,a,b,ch,prec);
802: case 3: return qromo (ep,a,b,ch,prec);
803: default: err(flagerr);
804: }
805: return NULL; /* not reached */
806: }
807:
808: #define JMAX 25
809: #define JMAXP JMAX+3
810: #define KLOC 4
811:
812: /* we need to make a copy in any case, cf forpari */
813: static GEN
814: fix(GEN a, long prec)
815: {
816: GEN p = cgetr(prec);
817: gaffect(a,p); return p;
818: }
819:
820: GEN polint_i(GEN xa, GEN ya, GEN x, long n, GEN *ptdy);
821:
822: GEN
823: qromb(entree *ep, GEN a, GEN b, char *ch, long prec)
824: {
825: GEN ss,dss,s,h,p1,p2,qlint,del,x,sum;
826: long av = avma, av1, tetpil,j,j1,j2,lim,it,sig;
827:
828: a = fix(a,prec);
829: b = fix(b,prec);
830: qlint=subrr(b,a); sig=signe(qlint);
831: if (!sig) { avma=av; return gzero; }
832: if (sig<0) { setsigne(qlint,1); s=a; a=b; b=s; }
833:
834: s=new_chunk(JMAXP);
835: h=new_chunk(JMAXP);
836: h[0] = (long)realun(prec);
837:
838: push_val(ep, a);
839: p1=lisexpr(ch); if (p1 == a) p1=rcopy(p1);
840: ep->value = (void*)b;
841: p2=lisexpr(ch);
842: s[0]=lmul2n(gmul(qlint,gadd(p1,p2)),-1);
843: for (it=1,j=1; j<JMAX; j++,it=it<<1)
844: {
845: h[j]=lshiftr((GEN)h[j-1],-2);
846: av1=avma; del=divrs(qlint,it); x=addrr(a,shiftr(del,-1));
847: for (sum=gzero,j1=1; j1<=it; j1++,x=addrr(x,del))
848: {
849: ep->value = (void*) x; sum=gadd(sum,lisexpr(ch));
850: }
851: sum = gmul(sum,del); p1 = gadd((GEN)s[j-1],sum);
852: tetpil = avma;
853: s[j]=lpile(av1,tetpil,gmul2n(p1,-1));
854:
855: if (j>=KLOC)
856: {
857: tetpil=avma; ss=polint_i(h+j-KLOC,s+j-KLOC,gzero,KLOC+1,&dss);
858: j1=gexpo(ss); j2=gexpo(dss); lim=bit_accuracy(prec)-j-6;
859: if (j1-j2 > lim || (j1 < -lim && j2<j1-1))
860: {
861: pop_val(ep); if (gcmp0(gimag(ss))) ss=greal(ss);
862: tetpil=avma; return gerepile(av,tetpil,gmulsg(sig,ss));
863: }
864: avma=tetpil;
865: }
866: }
867: err(intger2); return NULL; /* not reached */
868: }
869:
870: GEN
871: qromo(entree *ep, GEN a, GEN b, char *ch, long prec)
872: {
873: GEN ss,dss,s,h,p1,qlint,del,ddel,x,sum;
874: long av = avma,av1,tetpil,j,j1,j2,lim,it,sig;
875:
876: a = fix(a, prec);
877: b = fix(b, prec);
878: qlint=subrr(b,a); sig=signe(qlint);
879: if (!sig) { avma=av; return gzero; }
880: if (sig<0) { setsigne(qlint,1); s=a; a=b; b=s; }
881:
882: s=new_chunk(JMAXP);
883: h=new_chunk(JMAXP);
884: h[0] = (long)realun(prec);
885:
886: p1 = shiftr(addrr(a,b),-1); push_val(ep, p1);
887: p1=lisexpr(ch); s[0]=lmul(qlint,p1);
888:
889: for (it=1, j=1; j<JMAX; j++, it*=3)
890: {
891: h[j] = ldivrs((GEN)h[j-1],9);
892: av1=avma; del=divrs(qlint,3*it); ddel=shiftr(del,1);
893: x=addrr(a,shiftr(del,-1)); sum=gzero;
894: for (j1=1; j1<=it; j1++)
895: {
896: ep->value = (void*)x; sum=gadd(sum,lisexpr(ch)); x=addrr(x,ddel);
897: ep->value = (void*)x; sum=gadd(sum,lisexpr(ch)); x=addrr(x,del);
898: }
899: sum = gmul(sum,del); p1 = gdivgs((GEN)s[j-1],3);
900: tetpil = avma; s[j]=lpile(av1,tetpil,gadd(p1,sum));
901:
902: if (j>=KLOC)
903: {
904: tetpil=avma; ss=polint_i(h+j-KLOC,s+j-KLOC,gzero,KLOC+1,&dss);
905: j1=gexpo(ss); j2=gexpo(dss); lim=bit_accuracy(prec)-(3*j/2)-6;
906: if ( j1-j2 > lim || (j1 < -lim && j2<j1-1))
907: {
908: pop_val(ep); if (gcmp0(gimag(ss))) ss=greal(ss);
909: tetpil=avma; return gerepile(av,tetpil,gmulsg(sig,ss));
910: }
911: avma=tetpil;
912: }
913: }
914: err(intger2); return NULL; /* not reached */
915: }
916:
917: #undef JMAX
918: #undef JMAXP
919: #define JMAX 16
920: #define JMAXP JMAX+3
921:
922: GEN
923: qromi(entree *ep, GEN a, GEN b, char *ch, long prec)
924: {
925: GEN ss,dss,s,h,q1,p1,p,qlint,del,ddel,x,sum;
926: long av = avma, av1,tetpil,j,j1,j2,lim,it,sig;
927:
928: p=cgetr(prec); gaffect(ginv(a),p); a=p;
929: p=cgetr(prec); gaffect(ginv(b),p); b=p;
930: qlint=subrr(b,a); sig= -signe(qlint);
931: if (!sig) { avma=av; return gzero; }
932: if (sig>0) { setsigne(qlint,1); s=a; a=b; b=s; }
933:
934: s=new_chunk(JMAXP);
935: h=new_chunk(JMAXP);
936: h[0] = (long)realun(prec);
937:
938: x=divsr(2,addrr(a,b)); push_val(ep, x);
939: p1=gmul(lisexpr(ch),mulrr(x,x));
940: s[0]=lmul(qlint,p1);
941: for (it=1,j=1; j<JMAX; j++, it*=3)
942: {
943: h[j] = ldivrs((GEN)h[j-1],9);
944: av1=avma; del=divrs(qlint,3*it); ddel=shiftr(del,1);
945: x=addrr(a,shiftr(del,-1)); sum=gzero;
946: for (j1=1; j1<=it; j1++)
947: {
948: q1 = ginv(x); ep->value = (void*)q1;
949: p1=gmul(lisexpr(ch), gsqr(q1));
950: sum=gadd(sum,p1); x=addrr(x,ddel);
951:
952: q1 = ginv(x); ep->value = (void*)q1;
953: p1=gmul(lisexpr(ch), gsqr(q1));
954: sum=gadd(sum,p1); x=addrr(x,del);
955: }
956: sum = gmul(sum,del); p1 = gdivgs((GEN)s[j-1],3);
957: tetpil=avma;
958: s[j]=lpile(av1,tetpil,gadd(p1,sum));
959:
960: if (j>=KLOC)
961: {
962: tetpil=avma; ss=polint_i(h+j-KLOC,s+j-KLOC,gzero,KLOC+1,&dss);
963: j1=gexpo(ss); j2=gexpo(dss); lim=bit_accuracy(prec)-(3*j/2)-6;
964: if (j1-j2 > lim || (j1 < -lim && j2 < j1-1))
965: {
966: pop_val(ep); if (gcmp0(gimag(ss))) ss=greal(ss);
967: tetpil=avma; return gerepile(av,tetpil,gmulsg(sig,ss));
968: }
969: }
970: }
971: err(intger2); return NULL; /* not reached */
972: }
973:
974: GEN
975: rombint(entree *ep, GEN a, GEN b, char *ch, long prec)
976: {
977: GEN aa,bb,mun,p1,p2,p3;
978: long l,av,tetpil;
979:
980: l=gcmp(b,a); if (!l) return gzero;
981: if (l<0) { bb=a; aa=b; } else { bb=b; aa=a; }
982: av=avma; mun = negi(gun);
983: if (gcmpgs(bb,100)>=0)
984: {
985: if (gcmpgs(aa,1)>=0) return qromi(ep,a,b,ch,prec);
986: p1=qromi(ep,gun,bb,ch,prec);
987: if (gcmpgs(aa,-100)>=0)
988: {
989: p2=qromo(ep,aa,gun,ch,prec); tetpil=avma;
990: return gerepile(av,tetpil,gmulsg(l,gadd(p1,p2)));
991: }
992: p2=qromo(ep,mun,gun,ch,prec); p3=gadd(p2,qromi(ep,aa,mun,ch,prec));
993: tetpil=avma; return gerepile(av,tetpil,gmulsg(l,gadd(p1,p3)));
994: }
995: if (gcmpgs(aa,-100)>=0) return qromo(ep,a,b,ch,prec);
996: if (gcmpgs(bb,-1)>=0)
997: {
998: p1=qromi(ep,aa,mun,ch,prec); p2=qromo(ep,mun,bb,ch,prec); tetpil=avma;
999: return gerepile(av,tetpil,gmulsg(l,gadd(p1,p2)));
1000: }
1001: return qromi(ep,a,b,ch,prec);
1002: }
1003:
1004: /********************************************************************/
1005: /** **/
1006: /** ZAGIER POLYNOMIALS (for sumiter) **/
1007: /** **/
1008: /********************************************************************/
1009:
1010: GEN
1011: polzag(long n, long m)
1012: {
1013: long d1,d,r,k,av=avma,tetpil;
1014: GEN p1,p2,pol1,g,s;
1015:
1016: if (m>=n || m<0)
1017: err(talker,"first index must be greater than second in polzag");
1018: d1=n-m; d=d1<<1; d1--; p1=gsub(gun,gmul2n(polx[0],1));
1019: pol1=gsub(gun,polx[0]); p2=gmul(polx[0],pol1); r=(m+1)>>1;
1020: g=gzero;
1021: for (k=0; k<=d1; k++)
1022: {
1023: s=binome(stoi(d),(k<<1)+1); if (k&1) s=negi(s);
1024: g=gadd(g,gmul(s,gmul(gpuigs(polx[0],k),gpuigs(pol1,d1-k))));
1025: }
1026: g=gmul(g,gpuigs(p2,r));
1027: if (!(m&1)) g=gadd(gmul(p1,g),gmul2n(gmul(p2,derivpol(g)),1));
1028: for (k=1; k<=r; k++)
1029: {
1030: g=derivpol(g); g=gadd(gmul(p1,g),gmul2n(gmul(p2,derivpol(g)),1));
1031: }
1032: g = m ? gmul2n(g,(m-1)>>1):gmul2n(g,-1);
1033: s=mulsi(n-m,mpfact(m+1));
1034: tetpil=avma; return gerepile(av,tetpil,gdiv(g,s));
1035: }
1036:
1037: #ifdef _MSC_VER /* Bill Daly: work around a MSVC bug */
1038: #pragma optimize("g",off)
1039: #endif
1040: GEN
1041: polzagreel(long n, long m, long prec)
1042: {
1043: long d1,d,r,j,k,k2,av=avma,tetpil;
1044: GEN p2,pol1,g,h,v,b,gend,s;
1045:
1046: if (m>=n || m<0)
1047: err(talker,"first index must be greater than second in polzag");
1048: d1=n-m; d=d1<<1; d1--; pol1=gadd(gun,polx[0]);
1049: p2=gmul(polx[0],pol1); r=(m+1)>>1; gend=stoi(d);
1050: v=cgetg(d1+2,t_VEC); g=cgetg(d1+2,t_VEC);
1051: v[d1+1]=un; b=mulri(realun(prec),gend); g[d1+1]=(long)b;
1052: for (k=1; k<=d1; k++)
1053: {
1054: v[d1-k+1]=un;
1055: for (j=1; j<k; j++)
1056: v[d1-k+j+1]=(long)addii((GEN)v[d1-k+j+1],(GEN)v[d1-k+j+2]);
1057: k2=k+k; b=divri(mulri(b,mulss(d-k2+1,d-k2)),mulss(k2,k2+1));
1058: for (j=1; j<=k; j++)
1059: g[d1-k+j+1]=(long)mpadd((GEN)g[d1-k+j+1],mulri(b,(GEN)v[d1-k+j+1]));
1060: g[d1-k+1]=(long)b;
1061: }
1062: h=cgetg(d1+3,t_POL); h[1]=evalsigne(1)+evallgef(d1+3);
1063: for (k=0; k<=d1; k++) h[k+2]=g[k+1];
1064: g=gmul(h,gpuigs(p2,r));
1065: for (j=0; j<=r; j++)
1066: {
1067: if (j) g=derivpol(g);
1068: if (j || !(m&1))
1069: {
1070: h=cgetg(n+3,t_POL); h[1]=evalsigne(1)+evallgef(n+3);
1071: h[2]=g[2];
1072: for (k=1; k<n; k++)
1073: h[k+2]=ladd(gmulsg(k+k+1,(GEN)g[k+2]),gmulsg(k<<1,(GEN)g[k+1]));
1074: h[n+2]=lmulsg(n<<1,(GEN)g[n+1]); g=h;
1075: }
1076: }
1077: g = m ?gmul2n(g,(m-1)>>1):gmul2n(g,-1);
1078: s=mulsi(n-m,mpfact(m+1));
1079: tetpil=avma; return gerepile(av,tetpil,gdiv(g,s));
1080: }
1081: #ifdef _MSC_VER
1082: #pragma optimize("g",on)
1083: #endif
1084:
1085: /********************************************************************/
1086: /** **/
1087: /** SEARCH FOR REAL ZEROS of an expression **/
1088: /** **/
1089: /********************************************************************/
1090:
1091: GEN
1092: zbrent(entree *ep, GEN a, GEN b, char *ch, long prec)
1093: {
1094: long av=avma,tetpil,sig,iter,itmax;
1095: GEN c,d,e,tol,toli,min1,min2,fa,fb,fc,p,q,r,s,xm;
1096:
1097: a = fix(a,prec);
1098: b = fix(b,prec); sig=cmprr(b,a);
1099: if (!sig) { avma = av; return gzero; }
1100: if (sig<0) { c=a; a=b; b=c; }
1101:
1102: push_val(ep, a); fa = lisexpr(ch);
1103: ep->value = (void*)b; fb = lisexpr(ch);
1104: if (gsigne(fa)*gsigne(fb) > 0)
1105: err(talker,"roots must be bracketed in solve");
1106: itmax = (prec<<(TWOPOTBITS_IN_LONG+1)) + 1;
1107: tol = realun(3); setexpo(tol, 5-bit_accuracy(prec));
1108: fc=fb;
1109: for (iter=1; iter<=itmax; iter++)
1110: {
1111: if (gsigne(fb)*gsigne(fc) > 0)
1112: {
1113: c=a; fc=fa; d=subrr(b,a); e=d;
1114: }
1115: if (gcmp(gabs(fc,0),gabs(fb,0)) < 0)
1116: {
1117: a=b; b=c; c=a; fa=fb; fb=fc; fc=fa;
1118: }
1119: toli=mulrr(tol,gmax(tol,absr(b)));
1120: xm=shiftr(subrr(c,b),-1);
1121: if (cmprr(absr(xm),toli) <= 0 || gcmp0(fb))
1122: {
1123: tetpil=avma; pop_val(ep);
1124: return gerepile(av,tetpil,rcopy(b));
1125: }
1126: if (cmprr(absr(e),toli) >= 0 && gcmp(gabs(fa,0),gabs(fb,0)) > 0)
1127: {
1128: s=gdiv(fb,fa);
1129: if (cmprr(a,b)==0)
1130: {
1131: p=gmul2n(gmul(xm,s),1); q=gsubsg(1,s);
1132: }
1133: else
1134: {
1135: q=gdiv(fa,fc); r=gdiv(fb,fc);
1136: p=gmul2n(gmul(gsub(q,r),gmul(xm,q)),1);
1137: p=gmul(s,gsub(p,gmul(gsub(b,a),gsubgs(r,1))));
1138: q=gmul(gmul(gsubgs(q,1),gsubgs(r,1)),gsubgs(s,1));
1139: }
1140: if (gsigne(p) > 0) q=gneg_i(q); else p=gneg_i(p);
1141: min1=gsub(gmulsg(3,gmul(xm,q)),gabs(gmul(q,toli),0));
1142: min2=gabs(gmul(e,q),0);
1143: if (gcmp(gmul2n(p,1),gmin(min1,min2)) < 0)
1144: { e=d; d=gdiv(p,q); }
1145: else
1146: { d=xm; e=d; }
1147: }
1148: else { d=xm; e=d; }
1149: a=b; fa=fb;
1150: if (gcmp(gabs(d,0),toli) > 0) b=gadd(b,d);
1151: else
1152: {
1153: if (gsigne(xm)>0)
1154: b = addrr(b,toli);
1155: else
1156: b = subrr(b,toli);
1157: }
1158: ep->value = (void*)b; fb=lisexpr(ch);
1159: }
1160: err(talker,"too many iterations in solve");
1161: return NULL; /* not reached */
1162: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to sumiter.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari / src / language