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