Annotation of OpenXM_contrib/pari-2.2/src/language/sumiter.c, Revision 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