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