Annotation of OpenXM_contrib/pari-2.2/src/language/sumiter.c, Revision 1.2
1.2 ! noro 1: /* $Id: sumiter.c,v 1.36 2002/07/15 13:30:01 karim Exp $
1.1 noro 2:
3: Copyright (C) 2000 The PARI group.
4:
5: This file is part of the PARI/GP package.
6:
7: PARI/GP is free software; you can redistribute it and/or modify it under the
8: terms of the GNU General Public License as published by the Free Software
9: Foundation. It is distributed in the hope that it will be useful, but WITHOUT
10: ANY WARRANTY WHATSOEVER.
11:
12: Check the License for details. You should have received a copy of it, along
13: with the package; see the file 'COPYING'. If not, write to the Free Software
14: Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15:
16: /********************************************************************/
17: /** **/
18: /** SUMS, PRODUCTS, ITERATIONS **/
19: /** **/
20: /********************************************************************/
21: #include "pari.h"
22: #include "anal.h"
23: extern void changevalue_p(entree *ep, GEN x);
1.2 ! noro 24: extern GEN polint_i(GEN xa, GEN ya, GEN x, long n, GEN *ptdy);
1.1 noro 25:
26: /********************************************************************/
27: /** **/
28: /** ITERATIONS **/
29: /** **/
30: /********************************************************************/
31:
32: void
33: forpari(entree *ep, GEN a, GEN b, char *ch)
34: {
1.2 ! noro 35: gpmem_t av, av0 = avma, lim;
1.1 noro 36:
37: b = gcopy(b); av=avma; lim = stack_lim(av,1);
38: /* gcopy nedeed in case b gets overwritten in ch, as in
39: * b=10; for(a=1,b, print(a);b=1)
40: */
41: push_val(ep, a);
42: while (gcmp(a,b) <= 0)
43: {
1.2 ! noro 44: gpmem_t av1=avma; (void)lisseq(ch); avma=av1;
1.1 noro 45: if (loop_break()) break;
46: a = (GEN) ep->value; a = gadd(a,gun);
47: if (low_stack(lim, stack_lim(av,1)))
48: {
49: if (DEBUGMEM>1) err(warnmem,"forpari");
50: a = gerepileupto(av,a);
51: }
52: changevalue_p(ep,a);
53: }
54: pop_val(ep); avma = av0;
55: }
56:
57: static int negcmp(GEN x, GEN y) { return gcmp(y,x); }
58:
59: void
60: forstep(entree *ep, GEN a, GEN b, GEN s, char *ch)
61: {
1.2 ! noro 62: long ss, i;
! 63: gpmem_t av, av0 = avma, lim;
1.1 noro 64: GEN v = NULL;
65: int (*cmp)(GEN,GEN);
66:
67: b = gcopy(b); av=avma; lim = stack_lim(av,1);
68: push_val(ep, a);
69: if (is_vec_t(typ(s)))
70: {
71: v = s; s = gzero;
72: for (i=lg(v)-1; i; i--) s = gadd(s,(GEN)v[i]);
73: }
74: ss = gsigne(s);
75: if (!ss) err(talker, "step equal to zero in forstep");
76: cmp = (ss > 0)? gcmp: negcmp;
77: i = 0;
78: while (cmp(a,b) <= 0)
79: {
1.2 ! noro 80: gpmem_t av1=avma; (void)lisseq(ch); avma=av1;
1.1 noro 81: if (loop_break()) break;
82: if (v)
83: {
84: if (++i >= lg(v)) i = 1;
85: s = (GEN)v[i];
86: }
87: a = (GEN) ep->value; a = gadd(a,s);
88:
89: if (low_stack(lim, stack_lim(av,1)))
90: {
91: if (DEBUGMEM>1) err(warnmem,"forstep");
92: a = gerepileupto(av,a);
93: }
94: changevalue_p(ep,a);
95: }
96: pop_val(ep); avma = av0;
97: }
98:
99: /* assume ptr is the address of a diffptr containing the succesive
100: * differences between primes, and c = current prime (up to *p excluded)
101: * return smallest prime >= a, update ptr */
102: static long
103: sinitp(long a, long c, byteptr *ptr)
104: {
105: byteptr p = *ptr;
106: if (a <= 0) a = 2;
107: if (maxprime() < (ulong)a) err(primer1);
1.2 ! noro 108: while (a > c)
! 109: NEXT_PRIME_VIADIFF(c,p);
1.1 noro 110: *ptr = p; return c;
111: }
112:
113: /* value changed during the loop, replace by the first prime whose
114: value is strictly larger than new value */
115: static void
116: update_p(entree *ep, byteptr *ptr, long prime[])
117: {
118: GEN z = (GEN)ep->value;
119: long a, c;
120:
121: if (typ(z) == t_INT) a = 1; else { z = gceil(z); a = 0; }
122: if (is_bigint(z)) { prime[2] = -1; /* = infinity */ return; }
123: a += itos(z); c = prime[2];
124: if (c < a)
125: prime[2] = sinitp(a, c, ptr); /* increased */
126: else if (c > a)
127: { /* lowered, reset p */
128: *ptr = diffptr;
129: prime[2] = sinitp(a, 0, ptr);
130: }
1.2 ! noro 131: changevalue_p(ep, prime);
1.1 noro 132: }
133:
134: static byteptr
135: prime_loop_init(GEN ga, GEN gb, long *a, long *b, long prime[])
136: {
137: byteptr p = diffptr;
138: ulong P;
139:
140: ga = gceil(ga); gb = gfloor(gb);
141: if (typ(ga) != t_INT || typ(gb) != t_INT)
142: err(typeer,"prime_loop_init");
143: if (is_bigint(ga) || is_bigint(gb))
144: {
145: if (cmpii(ga, gb) > 0) return NULL;
146: err(primer1);
147: }
148: P = maxprime();
149: *a = itos(ga); if (*a <= 0) *a = 1;
150: *b = itos(gb); if (*a > *b) return NULL;
151: if ((ulong)*a <= P) prime[2] = sinitp(*a, 0, &p);
152: if ((ulong)*b > P) err(primer1);
153: return p;
154: }
155:
156: void
157: forprime(entree *ep, GEN ga, GEN gb, char *ch)
158: {
1.2 ! noro 159: long prime[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 0};
! 160: long a, b;
! 161: gpmem_t av = avma;
1.1 noro 162: byteptr p;
163:
164: p = prime_loop_init(ga,gb, &a,&b, prime);
165: if (!p) { avma = av; return; }
166:
167: avma = av; push_val(ep, prime);
168: while ((ulong)prime[2] < (ulong)b)
169: {
170: (void)lisseq(ch); if (loop_break()) break;
171: if (ep->value == prime)
172: prime[2] += *p++;
173: else
174: update_p(ep, &p, prime);
175: avma = av;
176: }
177: /* if b = P --> *p = 0 now and the loop wouldn't end if it read 'while
178: * (prime[2] <= b)' */
179: if (prime[2] == b) { (void)lisseq(ch); (void)loop_break(); avma = av; }
180: pop_val(ep);
181: }
182:
183: void
184: fordiv(GEN a, entree *ep, char *ch)
185: {
1.2 ! noro 186: long i, l;
! 187: gpmem_t av2, av = avma;
1.1 noro 188: GEN t = divisors(a);
189:
190: push_val(ep, NULL); l=lg(t); av2 = avma;
191: for (i=1; i<l; i++)
192: {
1.2 ! noro 193: ep->value = (void*) t[i];
1.1 noro 194: (void)lisseq(ch); if (loop_break()) break;
195: avma = av2;
196: }
197: pop_val(ep); avma=av;
198: }
199:
200: /* Embedded for loops:
201: * fl = 0: execute ch (a), where a = (ai) runs through all n-uplets in
202: * [m1,M1] x ... x [mn,Mn]
203: * fl = 1: impose a1 <= ... <= an
204: * fl = 2: a1 < ... < an
205: */
1.2 ! noro 206: typedef struct {
! 207: GEN *a, *m, *M;
! 208: long n,fl;
! 209: char *ch;
! 210: } fvdat;
1.1 noro 211:
212: /* general case */
213: static void
1.2 ! noro 214: fvloop(long i, fvdat *d)
1.1 noro 215: {
1.2 ! noro 216: d->a[i] = d->m[i];
! 217: if (d->fl && i > 1)
1.1 noro 218: {
1.2 ! noro 219: GEN p1 = gsub(d->a[i], d->a[i-1]);
1.1 noro 220: if (gsigne(p1) < 0)
1.2 ! noro 221: d->a[i] = gadd(d->a[i], gceil(gneg_i(p1)));
! 222: if (d->fl == 2 && gegal(d->a[i], d->a[i-1]))
! 223: d->a[i] = gadd(d->a[i], gun);
! 224: }
! 225: if (i+1 == d->n)
! 226: while (gcmp(d->a[i], d->M[i]) <= 0)
! 227: {
! 228: gpmem_t av = avma; (void)lisseq(d->ch); avma = av;
! 229: if (loop_break()) { d->n = 0; return; }
! 230: d->a[i] = gadd(d->a[i], gun);
1.1 noro 231: }
232: else
1.2 ! noro 233: while (gcmp(d->a[i], d->M[i]) <= 0)
1.1 noro 234: {
1.2 ! noro 235: gpmem_t av = avma; fvloop(i+1, d); avma = av;
! 236: if (!d->n) return;
! 237: d->a[i] = gadd(d->a[i], gun);
1.1 noro 238: }
239: }
240:
241: /* we run through integers */
242: static void
1.2 ! noro 243: fvloop_i(long i, fvdat *d)
1.1 noro 244: {
1.2 ! noro 245: d->a[i] = setloop(d->m[i]);
! 246: if (d->fl && i > 1)
1.1 noro 247: {
1.2 ! noro 248: int c = cmpii(d->a[i], d->a[i-1]);
1.1 noro 249: if (c < 0)
250: {
1.2 ! noro 251: d->a[i] = setloop(d->a[i-1]);
1.1 noro 252: c = 0;
253: }
1.2 ! noro 254: if (c == 0 && d->fl == 2)
! 255: d->a[i] = incloop(d->a[i]);
1.1 noro 256: }
1.2 ! noro 257: if (i+1 == d->n)
! 258: while (gcmp(d->a[i], d->M[i]) <= 0)
1.1 noro 259: {
1.2 ! noro 260: gpmem_t av = avma; (void)lisseq(d->ch); avma = av;
! 261: if (loop_break()) { d->n = 0; return; }
! 262: d->a[i] = incloop(d->a[i]);
1.1 noro 263: }
264: else
1.2 ! noro 265: while (gcmp(d->a[i], d->M[i]) <= 0)
1.1 noro 266: {
1.2 ! noro 267: gpmem_t av = avma; fvloop_i(i+1, d); avma = av;
! 268: if (!d->n) return;
! 269: d->a[i] = incloop(d->a[i]);
1.1 noro 270: }
271: }
272:
273: void
274: forvec(entree *ep, GEN x, char *c, long flag)
275: {
1.2 ! noro 276: gpmem_t av = avma;
! 277: long tx = typ(x);
! 278: fvdat D, *d = &D;
1.1 noro 279:
280: if (!is_vec_t(tx)) err(talker,"not a vector in forvec");
1.2 ! noro 281: if (flag<0 || flag>2) err(flagerr);
! 282: d->n = lg(x);
! 283: d->ch = c;
! 284: d->a = (GEN*)cgetg(d->n,t_VEC); push_val(ep, (GEN)d->a);
! 285: if (d->n == 1) (void)lisseq(d->ch);
1.1 noro 286: else
287: {
1.2 ! noro 288: long i, t = t_INT;
! 289: d->fl = flag;
! 290: d->m = (GEN*)cgetg(d->n,t_VEC);
! 291: d->M = (GEN*)cgetg(d->n,t_VEC);
! 292: for (i=1; i<d->n; i++)
! 293: {
! 294: GEN *e = (GEN*) x[i];
! 295: tx = typ(e);
! 296: if (! is_vec_t(tx) || lg(e)!=3)
1.1 noro 297: err(talker,"not a vector of two-component vectors in forvec");
1.2 ! noro 298: if (gcmp(e[1],e[2]) > 0) d->n = 0;
! 299: if (typ(e[1]) != t_INT) t = t_REAL;
! 300: /* in case x is an ep->value and lisexpr(d->ch) kills it, have to copy */
! 301: d->m[i] = gcopy(e[1]);
! 302: d->M[i] = gcopy(e[2]);
1.1 noro 303: }
1.2 ! noro 304: if (t == t_INT) fvloop_i(1, d); else fvloop(1, d);
1.1 noro 305: }
1.2 ! noro 306: pop_val(ep); avma = av;
1.1 noro 307: }
308:
309: /********************************************************************/
310: /** **/
311: /** SUMS **/
312: /** **/
313: /********************************************************************/
314:
315: GEN
316: somme(entree *ep, GEN a, GEN b, char *ch, GEN x)
317: {
1.2 ! noro 318: gpmem_t av, av0 = avma, lim;
1.1 noro 319: GEN p1;
320:
321: if (typ(a) != t_INT) err(talker,"non integral index in sum");
322: if (!x) x = gzero;
323: if (gcmp(b,a) < 0) return gcopy(x);
324:
1.2 ! noro 325: b = gfloor(b);
1.1 noro 326: a = setloop(a);
327: av=avma; lim = stack_lim(av,1);
328: push_val(ep, a);
329: for(;;)
330: {
331: p1 = lisexpr(ch); if (did_break()) err(breaker,"sum");
332: x=gadd(x,p1); if (cmpii(a,b) >= 0) break;
333: a = incloop(a);
334: if (low_stack(lim, stack_lim(av,1)))
335: {
336: if (DEBUGMEM>1) err(warnmem,"sum");
337: x = gerepileupto(av,x);
338: }
339: ep->value = (void*) a;
340: }
341: pop_val(ep); return gerepileupto(av0,x);
342: }
343:
344: GEN
345: suminf(entree *ep, GEN a, char *ch, long prec)
346: {
1.2 ! noro 347: long fl, G;
! 348: gpmem_t tetpil, av0 = avma, av, lim;
1.1 noro 349: GEN p1,x = realun(prec);
350:
351: if (typ(a) != t_INT) err(talker,"non integral index in suminf");
352: a = setloop(a);
353: av = avma; lim = stack_lim(av,1);
354: push_val(ep, a);
355: fl=0; G = bit_accuracy(prec) + 5;
356: for(;;)
357: {
358: p1 = lisexpr(ch); if (did_break()) err(breaker,"suminf");
359: x = gadd(x,p1); a = incloop(a);
360: if (gcmp0(p1) || gexpo(p1) <= gexpo(x)-G)
361: { if (++fl==3) break; }
362: else
363: fl=0;
364: if (low_stack(lim, stack_lim(av,1)))
365: {
366: if (DEBUGMEM>1) err(warnmem,"suminf");
367: x = gerepileupto(av,x);
368: }
369: ep->value = (void*)a;
370: }
371: pop_val(ep); tetpil=avma;
372: return gerepile(av0,tetpil,gsub(x,gun));
373: }
374:
375: GEN
376: divsum(GEN num, entree *ep, char *ch)
377: {
1.2 ! noro 378: gpmem_t av = avma;
1.1 noro 379: GEN z, y = gzero, t = divisors(num);
380: long i, l = lg(t);
381:
382: push_val(ep, NULL);
383: for (i=1; i<l; i++)
384: {
1.2 ! noro 385: ep->value = (void*) t[i];
1.1 noro 386: z = lisseq(ch);
387: if (did_break()) err(breaker,"divsum");
388: y = gadd(y, z);
389: }
390: pop_val(ep); return gerepileupto(av,y);
391: }
392:
393: /********************************************************************/
394: /** **/
395: /** PRODUCTS **/
396: /** **/
397: /********************************************************************/
398:
399: GEN
400: produit(entree *ep, GEN a, GEN b, char *ch, GEN x)
401: {
1.2 ! noro 402: gpmem_t av, av0 = avma, lim;
1.1 noro 403: GEN p1;
404:
405: if (typ(a) != t_INT) err(talker,"non integral index in sum");
406: if (!x) x = gun;
407: if (gcmp(b,a) < 0) return gcopy(x);
408:
1.2 ! noro 409: b = gfloor(b);
1.1 noro 410: a = setloop(a);
411: av=avma; lim = stack_lim(av,1);
412: push_val(ep, a);
413: for(;;)
414: {
415: p1 = lisexpr(ch); if (did_break()) err(breaker,"prod");
416: x = gmul(x,p1); if (cmpii(a,b) >= 0) break;
417: a = incloop(a);
418: if (low_stack(lim, stack_lim(av,1)))
419: {
420: if (DEBUGMEM>1) err(warnmem,"prod");
421: x = gerepileupto(av,x);
422: }
423: ep->value = (void*) a;
424: }
425: pop_val(ep); return gerepileupto(av0,x);
426: }
427:
428: GEN
429: prodinf0(entree *ep, GEN a, char *ch, long flag, long prec)
430: {
431: switch(flag)
432: {
433: case 0: return prodinf(ep,a,ch,prec);
434: case 1: return prodinf1(ep,a,ch,prec);
435: }
436: err(flagerr);
437: return NULL; /* not reached */
438: }
439:
440: GEN
441: prodinf(entree *ep, GEN a, char *ch, long prec)
442: {
1.2 ! noro 443: gpmem_t av0 = avma, av, lim;
1.1 noro 444: long fl,G;
445: GEN p1,x = realun(prec);
446:
447: if (typ(a) != t_INT) err(talker,"non integral index in prodinf");
448: a = setloop(a);
449: av = avma; lim = stack_lim(av,1);
450: push_val(ep, a);
451: fl=0; G = -bit_accuracy(prec)-5;
452: for(;;)
453: {
454: p1 = lisexpr(ch); if (did_break()) err(breaker,"prodinf");
455: x=gmul(x,p1); a = incloop(a);
456: if (gexpo(gsubgs(p1,1)) <= G) { if (++fl==3) break; } else fl=0;
457: if (low_stack(lim, stack_lim(av,1)))
458: {
459: if (DEBUGMEM>1) err(warnmem,"prodinf");
460: x = gerepileupto(av,x);
461: }
462: ep->value = (void*)a;
463: }
464: pop_val(ep); return gerepilecopy(av0,x);
465: }
466:
467: GEN
468: prodinf1(entree *ep, GEN a, char *ch, long prec)
469: {
1.2 ! noro 470: gpmem_t av0 = avma, av, lim;
1.1 noro 471: long fl,G;
472: GEN p1,p2,x = realun(prec);
473:
474: if (typ(a) != t_INT) err(talker,"non integral index in prodinf1");
475: a = setloop(a);
476: av = avma; lim = stack_lim(av,1);
477: push_val(ep, a);
478: fl=0; G = -bit_accuracy(prec)-5;
479: for(;;)
480: {
481: p2 = lisexpr(ch); if (did_break()) err(breaker,"prodinf1");
482: p1 = gadd(p2,gun); x=gmul(x,p1); a = incloop(a);
483: if (gcmp0(p1) || gexpo(p2) <= G) { if (++fl==3) break; } else fl=0;
484: if (low_stack(lim, stack_lim(av,1)))
485: {
486: if (DEBUGMEM>1) err(warnmem,"prodinf1");
487: x = gerepileupto(av,x);
488: }
489: ep->value = (void*)a;
490: }
491: pop_val(ep); return gerepilecopy(av0,x);
492: }
493:
494: GEN
495: prodeuler(entree *ep, GEN ga, GEN gb, char *ch, long prec)
496: {
1.2 ! noro 497: long prime[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 0};
1.1 noro 498: long a,b;
1.2 ! noro 499: gpmem_t av, av0 = avma, lim;
1.1 noro 500: GEN p1,x = realun(prec);
501: byteptr p;
1.2 ! noro 502:
1.1 noro 503: av = avma;
504: p = prime_loop_init(ga,gb, &a,&b, prime);
505: if (!p) { avma = av; return x; }
506:
507: av = avma; push_val(ep, prime);
508: lim = stack_lim(avma,1);
509: while ((ulong)prime[2] < (ulong)b)
510: {
511: p1 = lisexpr(ch); if (did_break()) err(breaker,"prodeuler");
512: x = gmul(x,p1);
513: if (low_stack(lim, stack_lim(av,1)))
514: {
515: if (DEBUGMEM>1) err(warnmem,"prodeuler");
516: x = gerepilecopy(av, x);
517: }
518: if (ep->value == prime)
519: prime[2] += *p++;
520: else
521: update_p(ep, &p, prime);
522: }
523: /* cf forprime */
524: if (prime[2] == b)
525: {
526: p1 = lisexpr(ch); if (did_break()) err(breaker,"prodeuler");
527: x = gmul(x,p1);
528: }
529: pop_val(ep); return gerepilecopy(av0,x);
530: }
531:
532: GEN
533: direulerall(entree *ep, GEN ga, GEN gb, char *ch, GEN c)
534: {
1.2 ! noro 535: long prime[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 0};
! 536: gpmem_t av0 = avma, av, lim = stack_lim(av0, 1);
1.1 noro 537: long p,n,i,j,k,tx,lx,a,b;
538: GEN x,y,s,polnum,polden;
539: byteptr d;
540:
541: d = prime_loop_init(ga,gb, &a,&b, prime);
542: n = c? itos(c): b;
1.2 ! noro 543: if (!d || b < 2 || n <= 0) return _vec(gun);
1.1 noro 544: if (n < b) b = n;
545: push_val(ep, prime);
546:
547: y = cgetg(n+1,t_VEC); av = avma;
548: x = cgetg(n+1,t_VEC); x[1]=un; for (i=2; i<=n; i++) x[i]=zero;
549: p = prime[2];
550: while (p <= b)
551: {
552: s = lisexpr(ch); if (did_break()) err(breaker,"direuler");
553: polnum = numer(s);
554: polden = denom(s);
555: tx = typ(polnum);
556: if (is_scalar_t(tx))
557: {
558: if (!gcmp1(polnum))
559: {
560: if (!gcmp_1(polnum))
561: err(talker,"constant term not equal to 1 in direuler");
562: polden = gneg(polden);
563: }
564: }
565: else
566: {
567: ulong k1,q, qlim;
568: if (tx != t_POL) err(typeer,"direuler");
569: c = truecoeff(polnum,0);
570: if (!gcmp1(c))
571: {
572: if (!gcmp_1(c))
573: err(talker,"constant term not equal to 1 in direuler");
574: polnum = gneg(polnum);
575: polden = gneg(polden);
576: }
577: for (i=1; i<=n; i++) y[i]=x[i];
578: lx=degpol(polnum);
579: q=p; qlim = n/p; j=1;
580: while (q<=(ulong)n && j<=lx)
581: {
582: c=(GEN)polnum[j+2];
583: if (!gcmp0(c))
584: for (k=1,k1=q; k1<=(ulong)n; k1+=q,k++)
585: x[k1] = ladd((GEN)x[k1], gmul(c,(GEN)y[k]));
586: if (q > qlim) break;
587: q*=p; j++;
588: }
589: }
590: tx = typ(polden);
591: if (is_scalar_t(tx))
592: {
593: if (!gcmp1(polden))
594: err(talker,"constant term not equal to 1 in direuler");
595: }
596: else
597: {
598: if (tx != t_POL) err(typeer,"direuler");
599: c = truecoeff(polden,0);
600: if (!gcmp1(c)) err(talker,"constant term not equal to 1 in direuler");
601: lx=degpol(polden);
602: for (i=p; i<=n; i+=p)
603: {
604: s=gzero; k=i; j=1;
605: while (!(k%p) && j<=lx)
606: {
607: c=(GEN)polden[j+2]; k/=p; j++;
608: if (!gcmp0(c)) s=gadd(s,gmul(c,(GEN)x[k]));
609: }
610: x[i]=lsub((GEN)x[i],s);
611: }
612: }
613: if (low_stack(lim, stack_lim(av,1)))
614: {
615: if (DEBUGMEM>1) err(warnmem,"direuler");
616: x = gerepilecopy(av, x);
617: }
618: p += *d++; if (!*d) err(primer1);
619: prime[2] = p;
620: }
621: pop_val(ep); return gerepilecopy(av0,x);
622: }
623:
624: GEN
625: direuler(entree *ep, GEN a, GEN b, char *ch)
626: {
627: return direulerall(ep,a,b,ch,NULL);
628: }
629:
630: /********************************************************************/
631: /** **/
632: /** VECTORS & MATRICES **/
633: /** **/
634: /********************************************************************/
635:
636: GEN
637: vecteur(GEN nmax, entree *ep, char *ch)
638: {
639: GEN y,p1;
640: long i,m;
1.2 ! noro 641: long c[]={evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 0};
1.1 noro 642:
643: if (typ(nmax)!=t_INT || signe(nmax) < 0)
644: err(talker,"bad number of components in vector");
645: m=itos(nmax); y=cgetg(m+1,t_VEC);
646: if (!ep || !ch)
647: {
648: for (i=1; i<=m; i++) y[i] = zero;
649: return y;
650: }
651: push_val(ep, c);
652: for (i=1; i<=m; i++)
653: {
654: c[2]=i; p1 = lisseq(ch);
655: if (did_break()) err(breaker,"vector");
656: y[i] = isonstack(p1)? (long)p1 : (long)forcecopy(p1);
657: }
658: pop_val(ep); return y;
659: }
660:
661: GEN
1.2 ! noro 662: vecteursmall(GEN nmax, entree *ep, char *ch)
! 663: {
! 664: GEN y,p1;
! 665: long i,m;
! 666: long c[]={evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 0};
! 667:
! 668: if (typ(nmax)!=t_INT || signe(nmax) < 0)
! 669: err(talker,"bad number of components in vector");
! 670: m=itos(nmax); y=cgetg(m+1,t_VECSMALL);
! 671: if (!ep || !ch)
! 672: {
! 673: for (i=1; i<=m; i++) y[i] = 0;
! 674: return y;
! 675: }
! 676: push_val(ep, c);
! 677: for (i=1; i<=m; i++)
! 678: {
! 679: c[2]=i; p1 = lisseq(ch);
! 680: if (did_break()) err(breaker,"vector");
! 681: y[i] = itos(p1);
! 682: }
! 683: pop_val(ep); return y;
! 684: }
! 685:
! 686:
! 687: GEN
1.1 noro 688: vvecteur(GEN nmax, entree *ep, char *ch)
689: {
690: GEN y = vecteur(nmax,ep,ch);
691: settyp(y,t_COL); return y;
692: }
693:
694: GEN
695: matrice(GEN nlig, GEN ncol,entree *ep1, entree *ep2, char *ch)
696: {
697: GEN y,z,p1;
698: long i,j,m,n,s;
1.2 ! noro 699: long c1[]={evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 1};
! 700: long c2[]={evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 1};
1.1 noro 701:
702: s=signe(ncol);
1.2 ! noro 703: if (ep1 == ep2 && ep1) err(talker, "identical index variables in matrix");
1.1 noro 704: if (typ(ncol)!=t_INT || s<0) err(talker,"bad number of columns in matrix");
705: if (!s) return cgetg(1,t_MAT);
706:
707: s=signe(nlig);
708: if (typ(nlig)!=t_INT || s<0) err(talker,"bad number of rows in matrix");
709: m=itos(ncol)+1;
710: n=itos(nlig)+1; y=cgetg(m,t_MAT);
711: if (!s)
712: {
713: for (i=1; i<m; i++) y[i]=lgetg(1,t_COL);
714: return y;
715: }
716: if (!ep1 || !ep2 || !ch)
717: {
718: for (i=1; i<m; i++)
719: {
720: z=cgetg(n,t_COL); y[i]=(long)z;
721: for (j=1; j<n; j++) z[j] = zero;
722: }
723: return y;
724: }
725: push_val(ep1, c1);
726: push_val(ep2, c2);
727: for (i=1; i<m; i++)
728: {
729: c2[2]=i; z=cgetg(n,t_COL); y[i]=(long)z;
730: for (j=1; j<n; j++)
731: {
732: c1[2]=j; p1=lisseq(ch);
733: if (did_break()) err(breaker,"matrix");
734: z[j] = isonstack(p1)? (long)p1 : (long)forcecopy(p1);
735: }
736: }
737: pop_val(ep2);
738: pop_val(ep1); return y;
739: }
740:
741: /********************************************************************/
742: /** **/
1.2 ! noro 743: /** SUMMING SERIES **/
1.1 noro 744: /** **/
745: /********************************************************************/
746:
747: GEN
748: sumalt0(entree *ep, GEN a, char *ch, long flag, long prec)
749: {
750: switch(flag)
751: {
752: case 0: return sumalt(ep,a,ch,prec);
753: case 1: return sumalt2(ep,a,ch,prec);
754: default: err(flagerr);
755: }
756: return NULL; /* not reached */
757: }
758:
759: GEN
760: sumpos0(entree *ep, GEN a, char *ch, long flag, long prec)
761: {
762: switch(flag)
763: {
764: case 0: return sumpos(ep,a,ch,prec);
765: case 1: return sumpos2(ep,a,ch,prec);
766: default: err(flagerr);
767: }
768: return NULL; /* not reached */
769: }
770:
771: GEN
772: sumalt(entree *ep, GEN a, char *ch, long prec)
773: {
1.2 ! noro 774: long k, N;
! 775: gpmem_t av = avma, tetpil;
1.1 noro 776: GEN s,az,c,x,e1,d;
777:
778: if (typ(a) != t_INT) err(talker,"non integral index in sumalt");
779:
780: push_val(ep, a);
1.2 ! noro 781: e1 = addsr(3,gsqrt(stoi(8),prec));
! 782: N = (long)(0.4*(bit_accuracy(prec) + 7));
! 783: d = gpowgs(e1,N);
! 784: d = shiftr(addrr(d, ginv(d)),-1);
! 785: az = negi(gun); c = d;
! 786: s = gzero;
1.1 noro 787: for (k=0; ; k++)
788: {
789: x = lisexpr(ch); if (did_break()) err(breaker,"sumalt");
790: c = gadd(az,c); s = gadd(s,gmul(x,c));
791: az = divii(mulii(mulss(N-k,N+k),shifti(az,1)),mulss(k+1,k+k+1));
792: if (k==N-1) break;
793: a = addsi(1,a); ep->value = (void*) a;
794: }
1.2 ! noro 795: tetpil = avma; pop_val(ep);
1.1 noro 796: return gerepile(av,tetpil,gdiv(s,d));
797: }
798:
799: GEN
800: sumalt2(entree *ep, GEN a, char *ch, long prec)
801: {
1.2 ! noro 802: long k, N;
! 803: gpmem_t av = avma, tetpil;
1.1 noro 804: GEN x,s,dn,pol;
805:
806: if (typ(a) != t_INT) err(talker,"non integral index in sumalt");
807:
808: push_val(ep, a);
1.2 ! noro 809: N = (long)(0.31*(bit_accuracy(prec) + 5));
! 810: pol = polzagreel(N,N>>1,prec+1);
! 811: dn = poleval(pol,gun);
! 812: pol[2] = lsub((GEN)pol[2],dn);
! 813: pol = gdiv(pol, gsub(polx[0],gun));
! 814: N = degpol(pol);
! 815: s = gzero;
! 816: for (k=0; k<=N; k++)
1.1 noro 817: {
818: x = lisexpr(ch); if (did_break()) err(breaker,"sumalt2");
1.2 ! noro 819: s = gadd(s, gmul(x,(GEN)pol[k+2]));
! 820: if (k == N) break;
1.1 noro 821: a = addsi(1,a); ep->value = (void*) a;
822: }
1.2 ! noro 823: tetpil = avma; pop_val(ep);
! 824: return gerepile(av,tetpil, gdiv(s,dn));
1.1 noro 825: }
826:
827: GEN
828: sumpos(entree *ep, GEN a, char *ch, long prec)
829: {
1.2 ! noro 830: long k, kk, N, G;
! 831: gpmem_t av = avma, tetpil;
1.1 noro 832: GEN p1,r,q1,reel,s,az,c,x,e1,d, *stock;
833:
834: if (typ(a) != t_INT) err(talker,"non integral index in sumpos");
835:
836: push_val(ep, NULL);
1.2 ! noro 837: a = subis(a,1); reel = cgetr(prec);
! 838: e1 = addsr(3,gsqrt(stoi(8),prec));
! 839: N = (long)(0.4*(bit_accuracy(prec) + 7));
! 840: d = gpowgs(e1,N); d = shiftr(addrr(d, ginv(d)),-1);
! 841: az = negi(gun); c = d; s = gzero;
! 842:
1.1 noro 843: G = -bit_accuracy(prec) - 5;
1.2 ! noro 844: stock = (GEN*)new_chunk(N+1); for (k=1; k<=N; k++) stock[k] = NULL;
1.1 noro 845: for (k=0; k<N; k++)
846: {
1.2 ! noro 847: if (odd(k) && stock[k]) x = stock[k];
1.1 noro 848: else
849: {
1.2 ! noro 850: x = gzero; r = stoi(2*k+2);
1.1 noro 851: for(kk=0;;kk++)
852: {
853: long ex;
854: q1 = addii(r,a); ep->value = (void*) q1;
855: p1 = lisexpr(ch); if (did_break()) err(breaker,"sumpos");
856: gaffect(p1,reel); ex = expo(reel)+kk; setexpo(reel,ex);
857: x = gadd(x,reel); if (kk && ex < G) break;
1.2 ! noro 858: r = shifti(r,1);
1.1 noro 859: }
1.2 ! noro 860: if (2*k < N) stock[2*k+1] = x;
1.1 noro 861: q1 = addsi(k+1,a); ep->value = (void*) q1;
862: p1 = lisexpr(ch); if (did_break()) err(breaker,"sumpos");
863: gaffect(p1,reel); x = gadd(reel, gmul2n(x,1));
864: }
1.2 ! noro 865: c = gadd(az,c); p1 = k&1? gneg_i(c): c;
1.1 noro 866: s = gadd(s,gmul(x,p1));
867: az = divii(mulii(mulss(N-k,N+k),shifti(az,1)),mulss(k+1,k+k+1));
868: }
869: tetpil=avma; pop_val(ep);
870: return gerepile(av,tetpil,gdiv(s,d));
871: }
872:
873: GEN
874: sumpos2(entree *ep, GEN a, char *ch, long prec)
875: {
1.2 ! noro 876: long k, kk, N, G;
! 877: gpmem_t av = avma, tetpil;
1.1 noro 878: GEN p1,r,q1,reel,s,pol,dn,x, *stock;
879:
880: if (typ(a) != t_INT) err(talker,"non integral index in sumpos2");
881:
882: push_val(ep, a);
1.2 ! noro 883: a = subis(a,1); reel = cgetr(prec);
! 884: N = (long)(0.31*(bit_accuracy(prec) + 5));
! 885:
1.1 noro 886: G = -bit_accuracy(prec) - 5;
1.2 ! noro 887: stock = (GEN*)new_chunk(N+1); for (k=1; k<=N; k++) stock[k] = NULL;
1.1 noro 888: for (k=1; k<=N; k++)
889: {
890: if (odd(k) || !stock[k])
891: {
1.2 ! noro 892: x = gzero; r = stoi(2*k);
1.1 noro 893: for(kk=0;;kk++)
894: {
895: long ex;
896: q1 = addii(r,a); ep->value = (void*) q1;
897: p1 = lisexpr(ch); if (did_break()) err(breaker,"sumpos2");
898: gaffect(p1,reel); ex = expo(reel)+kk; setexpo(reel,ex);
899: x=gadd(x,reel); if (kk && ex < G) break;
1.2 ! noro 900: r = shifti(r,1);
1.1 noro 901: }
1.2 ! noro 902: if (2*k-1 < N) stock[2*k] = x;
1.1 noro 903: q1 = addsi(k,a); ep->value = (void*) q1;
904: p1 = lisexpr(ch); if (did_break()) err(breaker,"sumpos2");
905: gaffect(p1,reel); stock[k] = gadd(reel, gmul2n(x,1));
906: }
907: }
908: pop_val(ep); s = gzero;
1.2 ! noro 909: pol = polzagreel(N,N>>1,prec+1);
! 910: dn = poleval(pol,gun);
! 911: pol[2] = lsub((GEN)pol[2],dn);
! 912: pol = gdiv(pol, gsub(gun,polx[0]));
1.1 noro 913: for (k=1; k<=lgef(pol)-2; k++)
914: {
915: p1 = gmul((GEN)pol[k+1],stock[k]);
1.2 ! noro 916: if (odd(k)) p1 = gneg_i(p1);
1.1 noro 917: s = gadd(s,p1);
918: }
919: tetpil=avma; return gerepile(av,tetpil,gdiv(s,dn));
920: }
921:
922: /********************************************************************/
923: /** **/
924: /** NUMERICAL INTEGRATION (Romberg) **/
925: /** **/
926: /********************************************************************/
1.2 ! noro 927: typedef struct {
! 928: entree *ep;
! 929: char *ch;
! 930: } exprdat;
! 931:
! 932: typedef struct {
! 933: GEN (*f)(void *, GEN);
! 934: void *E;
! 935: } invfun;
! 936:
! 937: /* we need to make a copy in any case, cf forpari */
! 938: static GEN
! 939: fix(GEN a, long prec)
! 940: {
! 941: GEN p = cgetr(prec);
! 942: gaffect(a,p); return p;
! 943: }
! 944:
! 945: #if 0
1.1 noro 946: GEN
1.2 ! noro 947: int_loop(entree *ep, char *ch)
1.1 noro 948: {
1.2 ! noro 949: gpmem_t av = avma;
! 950: intstruct T;
! 951:
! 952: T.in = NULL;
! 953: for(;;)
1.1 noro 954: {
1.2 ! noro 955: GEN x = Next(&T);
! 956: if (!x) return gerepileupto(av, T.out);
! 957: ep->value = (void*)x;
! 958: T.in = lisexpr(ch);
1.1 noro 959: }
960: }
1.2 ! noro 961: #endif
1.1 noro 962:
1.2 ! noro 963: /* f(x) */
! 964: static GEN
! 965: _gp_eval(void *dat, GEN x)
! 966: {
! 967: exprdat *E = (exprdat*)dat;
! 968: E->ep->value = x;
! 969: return lisexpr(E->ch);
! 970: }
! 971: /* 1/x^2 f(1/x) */
! 972: static GEN
! 973: _invf(void *dat, GEN x)
! 974: {
! 975: invfun *S = (invfun*)dat;
! 976: GEN y = ginv(x);
! 977: return gmul(S->f(S->E, y), gsqr(y));
! 978: }
1.1 noro 979:
1.2 ! noro 980: #define swap(a,b) { GEN _x = a; a = b; b = _x; }
1.1 noro 981: static GEN
1.2 ! noro 982: interp(GEN h, GEN s, long j, long lim, long KLOC)
1.1 noro 983: {
1.2 ! noro 984: gpmem_t av = avma;
! 985: long e1,e2;
! 986: GEN dss, ss = polint_i(h+j-KLOC,s+j-KLOC,gzero,KLOC+1,&dss);
! 987:
! 988: e1 = gexpo(ss);
! 989: e2 = gexpo(dss);
! 990: if (e1-e2 <= lim && e1 >= -lim) { avma = av; return NULL; }
! 991: if (gcmp0(gimag(ss))) ss = greal(ss);
! 992: return ss;
1.1 noro 993: }
994:
1.2 ! noro 995: static GEN
! 996: qrom3(void *dat, GEN (*eval)(void *, GEN), GEN a, GEN b, long prec)
1.1 noro 997: {
1.2 ! noro 998: const long JMAX = 25, KLOC = 4;
! 999: GEN ss,s,h,p1,p2,qlint,del,x,sum;
! 1000: long j, j1, it, sig;
! 1001: gpmem_t av;
1.1 noro 1002:
1003: a = fix(a,prec);
1004: b = fix(b,prec);
1.2 ! noro 1005: qlint = subrr(b,a); sig = signe(qlint);
! 1006: if (!sig) return gzero;
! 1007: if (sig < 0) { setsigne(qlint,1); swap(a,b); }
1.1 noro 1008:
1.2 ! noro 1009: s = new_chunk(JMAX+KLOC-1);
! 1010: h = new_chunk(JMAX+KLOC-1);
1.1 noro 1011: h[0] = (long)realun(prec);
1012:
1.2 ! noro 1013: p1 = eval(dat, a); if (p1 == a) p1 = rcopy(p1);
! 1014: p2 = eval(dat, b);
! 1015: s[0] = lmul2n(gmul(qlint,gadd(p1,p2)),-1);
! 1016: for (it=1,j=1; j<JMAX; j++, it<<=1)
! 1017: {
! 1018: h[j] = lshiftr((GEN)h[j-1],-2);
! 1019: av = avma; del = divrs(qlint,it);
! 1020: x = addrr(a, shiftr(del,-1));
! 1021: for (sum = gzero, j1 = 1; j1 <= it; j1++, x = addrr(x,del))
! 1022: sum = gadd(sum, eval(dat, x));
! 1023: sum = gmul(sum,del); p1 = gadd((GEN)s[j-1], sum);
! 1024: s[j] = lpileupto(av, gmul2n(p1,-1));
! 1025:
! 1026: if (j >= KLOC && (ss = interp(h, s, j, bit_accuracy(prec)-j-6, KLOC)))
! 1027: return gmulsg(sig,ss);
1.1 noro 1028: }
1.2 ! noro 1029: return NULL;
1.1 noro 1030: }
1031:
1.2 ! noro 1032: static GEN
! 1033: qrom2(void *dat, GEN (*eval)(void *, GEN), GEN a, GEN b, long prec)
1.1 noro 1034: {
1.2 ! noro 1035: const long JMAX = 16, KLOC = 4;
! 1036: GEN ss,s,h,p1,qlint,del,ddel,x,sum;
! 1037: long j, j1, it, sig;
! 1038: gpmem_t av;
1.1 noro 1039:
1040: a = fix(a, prec);
1041: b = fix(b, prec);
1.2 ! noro 1042: qlint = subrr(b,a); sig = signe(qlint);
! 1043: if (!sig) return gzero;
! 1044: if (sig < 0) { setsigne(qlint,1); swap(a,b); }
1.1 noro 1045:
1.2 ! noro 1046: s = new_chunk(JMAX+KLOC-1);
! 1047: h = new_chunk(JMAX+KLOC-1);
1.1 noro 1048: h[0] = (long)realun(prec);
1049:
1.2 ! noro 1050: p1 = shiftr(addrr(a,b),-1);
! 1051: s[0] = lmul(qlint, eval(dat, p1));
1.1 noro 1052: for (it=1, j=1; j<JMAX; j++, it*=3)
1053: {
1054: h[j] = ldivrs((GEN)h[j-1],9);
1.2 ! noro 1055: av = avma; del = divrs(qlint,3*it); ddel = shiftr(del,1);
! 1056: x = addrr(a, shiftr(del,-1));
! 1057: for (sum = gzero, j1 = 1; j1 <= it; j1++)
1.1 noro 1058: {
1.2 ! noro 1059: sum = gadd(sum, eval(dat, x)); x = addrr(x,ddel);
! 1060: sum = gadd(sum, eval(dat, x)); x = addrr(x,del);
1.1 noro 1061: }
1062: sum = gmul(sum,del); p1 = gdivgs((GEN)s[j-1],3);
1.2 ! noro 1063: s[j] = lpileupto(av, gadd(p1,sum));
1.1 noro 1064:
1.2 ! noro 1065: if (j >= KLOC && (ss = interp(h, s, j, bit_accuracy(prec)-(3*j/2)-6, KLOC)))
! 1066: return gmulsg(sig, ss);
1.1 noro 1067: }
1.2 ! noro 1068: return NULL;
1.1 noro 1069: }
1070:
1.2 ! noro 1071: /* integrate after change of variables x --> 1/x */
! 1072: static GEN
! 1073: qromi(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, long prec)
! 1074: {
! 1075: GEN A = ginv(b), B = ginv(a);
! 1076: invfun S;
! 1077: S.f = eval;
! 1078: S.E = E; return qrom2(&S, &_invf, A, B, prec);
! 1079: }
1.1 noro 1080:
1.2 ! noro 1081: /* a < b, assume b "small" (< 100 say) */
! 1082: static GEN
! 1083: rom_bsmall(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, long prec)
1.1 noro 1084: {
1.2 ! noro 1085: if (gcmpgs(a,-100) >= 0) return qrom2(E,eval,a,b,prec);
! 1086: if (b == gun || gcmpgs(b, -1) >= 0)
! 1087: { /* a < -100, b >= -1 */
! 1088: GEN _1 = negi(gun); /* split at -1 */
! 1089: return gadd(qromi(E,eval,a,_1,prec),
! 1090: qrom2(E,eval,_1,b,prec));
! 1091: }
! 1092: /* a < -100, b < -1 */
! 1093: return qromi(E,eval,a,b,prec);
! 1094: }
1.1 noro 1095:
1.2 ! noro 1096: static GEN
! 1097: rombint(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, long prec)
! 1098: {
! 1099: long l = gcmp(b,a);
! 1100: GEN z;
1.1 noro 1101:
1.2 ! noro 1102: if (!l) return gzero;
! 1103: if (l < 0) swap(a,b);
! 1104: if (gcmpgs(b,100) >= 0)
! 1105: {
! 1106: if (gcmpgs(a,1) >= 0)
! 1107: z = qromi(E,eval,a,b,prec);
! 1108: else /* split at 1 */
! 1109: z = gadd(rom_bsmall(E,eval,a,gun,prec), qromi(E,eval,gun,b,prec));
! 1110: }
! 1111: else
! 1112: z = rom_bsmall(E,eval,a,b,prec);
! 1113: if (l < 0) z = gneg(z);
! 1114: return z;
! 1115: }
1.1 noro 1116:
1.2 ! noro 1117: GEN
! 1118: intnum(void *E, GEN (*eval)(void*,GEN), GEN a, GEN b, long flag, long prec)
! 1119: {
! 1120: gpmem_t av = avma;
! 1121: GEN z;
! 1122: switch(flag)
1.1 noro 1123: {
1.2 ! noro 1124: case 0: z = qrom3 (E,eval,a,b,prec); break;
! 1125: case 1: z = rombint(E,eval,a,b,prec); break;
! 1126: case 2: z = qromi (E,eval,a,b,prec); break;
! 1127: case 3: z = qrom2 (E,eval,a,b,prec); break;
! 1128: default: err(flagerr); z = NULL;
1.1 noro 1129: }
1.2 ! noro 1130: if (!z) err(intger2);
! 1131: return gerepileupto(av, z);
1.1 noro 1132: }
1133:
1134: GEN
1.2 ! noro 1135: intnum0(entree *ep, GEN a, GEN b, char *ch, long flag, long prec)
1.1 noro 1136: {
1.2 ! noro 1137: exprdat E;
! 1138: GEN z;
1.1 noro 1139:
1.2 ! noro 1140: E.ch = ch;
! 1141: E.ep = ep; push_val(ep, NULL);
! 1142: z = intnum(&E, &_gp_eval, a,b, flag, prec);
! 1143: pop_val(ep); return z;
1.1 noro 1144: }
1145: /********************************************************************/
1146: /** **/
1.2 ! noro 1147: /** ZAGIER POLYNOMIALS (for sumalt) **/
1.1 noro 1148: /** **/
1149: /********************************************************************/
1150:
1151: GEN
1152: polzag(long n, long m)
1153: {
1.2 ! noro 1154: gpmem_t av = avma;
! 1155: long d1, d, r, k;
! 1156: GEN A, B, Bx, g, s;
1.1 noro 1157:
1.2 ! noro 1158: if (m >= n || m < 0)
1.1 noro 1159: err(talker,"first index must be greater than second in polzag");
1.2 ! noro 1160: d1 = n-m; d = d1<<1; d1--;
! 1161: A = gsub(gun, gmul2n(polx[0],1)); /* 1 - 2x */
! 1162: B = gsub(gun, polx[0]); /* 1 - x */
! 1163: Bx = gmul(polx[0],B); /* x - x^2 */
! 1164: r = (m+1)>>1;
! 1165: g = gzero;
1.1 noro 1166: for (k=0; k<=d1; k++)
1167: {
1.2 ! noro 1168: s = binome(stoi(d), (k<<1)+1); if (k&1) s = negi(s);
! 1169: g = gadd(g, gmul(s, gmul(gpowgs(polx[0],k), gpowgs(B,d1-k))));
1.1 noro 1170: }
1.2 ! noro 1171: g = gmul(g, gpowgs(Bx,r));
! 1172: if (!(m&1)) g = gadd(gmul(A,g), gmul2n(gmul(Bx,derivpol(g)),1));
1.1 noro 1173: for (k=1; k<=r; k++)
1174: {
1.2 ! noro 1175: g = derivpol(g);
! 1176: g = gadd(gmul(A,g), gmul2n(gmul(Bx,derivpol(g)),1));
1.1 noro 1177: }
1.2 ! noro 1178: g = m? gmul2n(g,(m-1)>>1): gmul2n(g,-1);
! 1179: s = mulsi(n-m, mpfact(m+1));
! 1180: return gerepileupto(av, gdiv(g,s));
1.1 noro 1181: }
1182:
1183: #ifdef _MSC_VER /* Bill Daly: work around a MSVC bug */
1184: #pragma optimize("g",off)
1185: #endif
1186: GEN
1187: polzagreel(long n, long m, long prec)
1188: {
1.2 ! noro 1189: long d1, d, r, j, k, k2;
! 1190: gpmem_t av = avma;
! 1191: GEN Bx,B,g,h,v,b,s;
1.1 noro 1192:
1.2 ! noro 1193: if (m >= n || m < 0)
1.1 noro 1194: err(talker,"first index must be greater than second in polzag");
1.2 ! noro 1195: d1 = n-m; d = d1<<1; d1--;
! 1196: B = gadd(gun,polx[0]);
! 1197: Bx = gmul(polx[0],B);
! 1198: r = (m+1)>>1;
! 1199: v = cgetg(d1+2,t_VEC);
! 1200: g = cgetg(d1+2,t_VEC);
! 1201: v[d1+1] = un; b = stor(d, prec);
! 1202: g[d1+1] = (long)b;
1.1 noro 1203: for (k=1; k<=d1; k++)
1204: {
1.2 ! noro 1205: v[d1-k+1] = un;
1.1 noro 1206: for (j=1; j<k; j++)
1.2 ! noro 1207: v[d1-k+j+1] = laddii((GEN)v[d1-k+j+1], (GEN)v[d1-k+j+2]);
! 1208: k2 = k+k; b = divri(mulri(b,mulss(d-k2+1,d-k2)), mulss(k2,k2+1));
1.1 noro 1209: for (j=1; j<=k; j++)
1.2 ! noro 1210: g[d1-k+j+1] = lmpadd((GEN)g[d1-k+j+1], mulri(b,(GEN)v[d1-k+j+1]));
! 1211: g[d1-k+1] = (long)b;
1.1 noro 1212: }
1.2 ! noro 1213: g = gmul(vec_to_pol(g,0), gpowgs(Bx,r));
1.1 noro 1214: for (j=0; j<=r; j++)
1215: {
1.2 ! noro 1216: if (j) g = derivpol(g);
1.1 noro 1217: if (j || !(m&1))
1218: {
1.2 ! noro 1219: h = cgetg(n+3,t_POL);
! 1220: h[1] = evalsigne(1)|evallgef(n+3);
! 1221: h[2] = g[2];
1.1 noro 1222: for (k=1; k<n; k++)
1.2 ! noro 1223: h[k+2] = ladd(gmulsg(k+k+1,(GEN)g[k+2]), gmulsg(k<<1,(GEN)g[k+1]));
! 1224: h[n+2] = lmulsg(n<<1, (GEN)g[n+1]);
! 1225: g = h;
1.1 noro 1226: }
1227: }
1.2 ! noro 1228: g = m? gmul2n(g,(m-1)>>1): gmul2n(g,-1);
! 1229: s = mulsi(n-m, mpfact(m+1));
! 1230: return gerepileupto(av, gdiv(g,s));
1.1 noro 1231: }
1232: #ifdef _MSC_VER
1233: #pragma optimize("g",on)
1234: #endif
1235:
1236: /********************************************************************/
1237: /** **/
1238: /** SEARCH FOR REAL ZEROS of an expression **/
1239: /** **/
1240: /********************************************************************/
1241: /* Brent's method, [a,b] bracketing interval */
1242: GEN
1243: zbrent(entree *ep, GEN a, GEN b, char *ch, long prec)
1244: {
1.2 ! noro 1245: long sig, iter, itmax;
! 1246: gpmem_t av = avma;
1.1 noro 1247: GEN c,d,e,tol,tol1,min1,min2,fa,fb,fc,p,q,r,s,xm;
1248:
1249: a = fix(a,prec);
1250: b = fix(b,prec); sig=cmprr(b,a);
1251: if (!sig) { avma = av; return gzero; }
1252: if (sig < 0) { c=a; a=b; b=c; } else c = b;
1253:
1254: push_val(ep, a); fa = lisexpr(ch);
1255: ep->value = (void*)b; fb = lisexpr(ch);
1256: if (gsigne(fa)*gsigne(fb) > 0)
1257: err(talker,"roots must be bracketed in solve");
1258: itmax = (prec<<(TWOPOTBITS_IN_LONG+1)) + 1;
1259: tol = realun(3); setexpo(tol, 5-bit_accuracy(prec));
1260: fc=fb;
1261: e = d = NULL; /* gcc -Wall */
1262: for (iter=1; iter<=itmax; iter++)
1263: {
1264: if (gsigne(fb)*gsigne(fc) > 0)
1265: {
1266: c = a; fc = fa; e = d = subrr(b,a);
1267: }
1268: if (gcmp(gabs(fc,0),gabs(fb,0)) < 0)
1269: {
1270: a = b; b = c; c = a; fa = fb; fb = fc; fc = fa;
1271: }
1272: tol1 = mulrr(tol, gmax(tol,absr(b)));
1273: xm = shiftr(subrr(c,b),-1);
1274: if (cmprr(absr(xm),tol1) <= 0 || gcmp0(fb)) break; /* SUCCESS */
1275:
1276: if (cmprr(absr(e),tol1) >= 0 && gcmp(gabs(fa,0),gabs(fb,0)) > 0)
1277: { /* attempt interpolation */
1278: s = gdiv(fb,fa);
1279: if (cmprr(a,c)==0)
1280: {
1281: p = gmul2n(gmul(xm,s),1); q = gsubsg(1,s);
1282: }
1283: else
1284: {
1285: q = gdiv(fa,fc); r = gdiv(fb,fc);
1286: p = gmul2n(gmul(gsub(q,r),gmul(xm,q)),1);
1287: p = gmul(s, gsub(p, gmul(gsub(b,a),gsubgs(r,1))));
1288: q = gmul(gmul(gsubgs(q,1),gsubgs(r,1)),gsubgs(s,1));
1289: }
1290: if (gsigne(p) > 0) q = gneg_i(q); else p = gneg_i(p);
1291: min1 = gsub(gmulsg(3,gmul(xm,q)), gabs(gmul(q,tol1),0));
1292: min2 = gabs(gmul(e,q),0);
1293: if (gcmp(gmul2n(p,1), gmin(min1,min2)) < 0)
1294: { e = d; d = gdiv(p,q); } /* interpolation OK */
1295: else
1296: { d = xm; e = d; } /* failed, use bisection */
1297: }
1298: else { d = xm; e = d; } /* bound decreasing too slowly, use bisection */
1299: a = b; fa = fb;
1300: if (gcmp(gabs(d,0),tol1) > 0) b = gadd(b,d);
1.2 ! noro 1301: else if (gsigne(xm) > 0) b = addrr(b,tol1);
! 1302: else b = subrr(b,tol1);
1.1 noro 1303: ep->value = (void*)b; fb = lisexpr(ch);
1304: }
1305: if (iter > itmax) err(talker,"too many iterations in solve");
1306: pop_val(ep);
1307: return gerepileuptoleaf(av, rcopy(b));
1308: }
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-2.2 / src / language