OpenXM_contrib/pari-2.2/src/language/sumiter.c - diff

Return to sumiter.c CVS log

Up to [local] / OpenXM_contrib / pari-2.2 / src / language

Diff for /OpenXM_contrib/pari-2.2/src/language/Attic/sumiter.c between version 1.1 and 1.2

version 1.1, 2001/10/02 11:17:10

version 1.2, 2002/09/11 07:27:04

Line 21 Foundation, Inc., 59 Temple Place - Suite 330, Boston,

#include "pari.h"

#include "anal.h"

extern void changevalue_p(entree *ep, GEN x);

extern GEN polint_i(GEN xa, GEN ya, GEN x, long n, GEN *ptdy);

/********************************************************************/

/** **/

Line 31 extern void changevalue_p(entree *ep, GEN x);

Line 32 extern void changevalue_p(entree *ep, GEN x);

void

forpari(entree *ep, GEN a, GEN b, char *ch)

{

ulong av,av0 = avma, lim;

gpmem_t av, av0 = avma, lim;

b = gcopy(b); av=avma; lim = stack_lim(av,1);

/* gcopy nedeed in case b gets overwritten in ch, as in

Line 40 forpari(entree *ep, GEN a, GEN b, char *ch)

Line 41 forpari(entree *ep, GEN a, GEN b, char *ch)

push_val(ep, a);

while (gcmp(a,b) <= 0)

{

ulong av1=avma; (void)lisseq(ch); avma=av1;

gpmem_t av1=avma; (void)lisseq(ch); avma=av1;

if (loop_break()) break;

a = (GEN) ep->value; a = gadd(a,gun);

if (low_stack(lim, stack_lim(av,1)))

Line 58 static int negcmp(GEN x, GEN y) { return gcmp(y,x); }

Line 59 static int negcmp(GEN x, GEN y) { return gcmp(y,x); }

void

forstep(entree *ep, GEN a, GEN b, GEN s, char *ch)

{

long ss, av,av0 = avma, lim, i;

long ss, i;

gpmem_t av, av0 = avma, lim;

GEN v = NULL;

int (*cmp)(GEN,GEN);

Line 75 forstep(entree *ep, GEN a, GEN b, GEN s, char *ch)

Line 77 forstep(entree *ep, GEN a, GEN b, GEN s, char *ch)

i = 0;

while (cmp(a,b) <= 0)

{

long av1=avma; (void)lisseq(ch); avma=av1;

gpmem_t av1=avma; (void)lisseq(ch); avma=av1;

if (loop_break()) break;

if (v)

{

Line 103 sinitp(long a, long c, byteptr *ptr)

Line 105 sinitp(long a, long c, byteptr *ptr)

byteptr p = *ptr;

if (a <= 0) a = 2;

if (maxprime() < (ulong)a) err(primer1);

while (a > c) c += *p++;

while (a > c)

NEXT_PRIME_VIADIFF(c,p);

*ptr = p; return c;

}

Line 125 update_p(entree *ep, byteptr *ptr, long prime[])

Line 128 update_p(entree *ep, byteptr *ptr, long prime[])

*ptr = diffptr;

prime[2] = sinitp(a, 0, ptr);

}

changevalue_p(ep, prime);

}

static byteptr

Line 153 prime_loop_init(GEN ga, GEN gb, long *a, long *b, long

Line 156 prime_loop_init(GEN ga, GEN gb, long *a, long *b, long

void

forprime(entree *ep, GEN ga, GEN gb, char *ch)

{

long prime[] = {evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3), 0};

long prime[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 0};

long av = avma, a,b;

long a, b;

gpmem_t av = avma;

byteptr p;

p = prime_loop_init(ga,gb, &a,&b, prime);

Line 179 forprime(entree *ep, GEN ga, GEN gb, char *ch)

Line 183 forprime(entree *ep, GEN ga, GEN gb, char *ch)

void

fordiv(GEN a, entree *ep, char *ch)

{

long i,av2,l, av = avma;

long i, l;

gpmem_t av2, av = avma;

GEN t = divisors(a);

push_val(ep, NULL); l=lg(t); av2 = avma;

for (i=1; i<l; i++)

{

ep->value = (void*) t[i];

(void)lisseq(ch); if (loop_break()) break;

avma = av2;

}

Line 198 fordiv(GEN a, entree *ep, char *ch)

Line 203 fordiv(GEN a, entree *ep, char *ch)

* fl = 1: impose a1 <= ... <= an

* fl = 2: a1 < ... < an

static GEN *fv_a, *fv_m, *fv_M;

typedef struct {

static long fv_n, fv_fl;

GEN *a, *m, *M;

static char *fv_ch;

long n,fl;

char *ch;

} fvdat;

/* general case */

static void

fvloop(long i)

fvloop(long i, fvdat *d)

{

fv_a[i] = fv_m[i];

d->a[i] = d->m[i];

if (fv_fl && i > 1)

if (d->fl && i > 1)

{

GEN p1 = gsub(fv_a[i], fv_a[i-1]);

GEN p1 = gsub(d->a[i], d->a[i-1]);

if (gsigne(p1) < 0)

fv_a[i] = gadd(fv_a[i], gceil(gneg_i(p1)));

d->a[i] = gadd(d->a[i], gceil(gneg_i(p1)));

if (fv_fl == 2 && gegal(fv_a[i], fv_a[i-1]))

if (d->fl == 2 && gegal(d->a[i], d->a[i-1]))

fv_a[i] = gadd(fv_a[i], gun);

d->a[i] = gadd(d->a[i], gun);

}

if (i+1 == fv_n)

if (i+1 == d->n)

while (gcmp(fv_a[i], fv_M[i]) <= 0)

while (gcmp(d->a[i], d->M[i]) <= 0)

{

long av = avma; (void)lisseq(fv_ch); avma = av;

gpmem_t av = avma; (void)lisseq(d->ch); avma = av;

if (loop_break()) { fv_n = 0; return; }

if (loop_break()) { d->n = 0; return; }

fv_a[i] = gadd(fv_a[i], gun);

d->a[i] = gadd(d->a[i], gun);

}

else

while (gcmp(fv_a[i], fv_M[i]) <= 0)

while (gcmp(d->a[i], d->M[i]) <= 0)

{

long av = avma; fvloop(i+1); avma = av;

gpmem_t av = avma; fvloop(i+1, d); avma = av;

if (!fv_n) return;

if (!d->n) return;

fv_a[i] = gadd(fv_a[i], gun);

d->a[i] = gadd(d->a[i], gun);

}

/* we run through integers */

static void

fvloop_i(long i)

fvloop_i(long i, fvdat *d)

{

fv_a[i] = setloop(fv_m[i]);

d->a[i] = setloop(d->m[i]);

if (fv_fl && i > 1)

if (d->fl && i > 1)

{

int c = cmpii(fv_a[i], fv_a[i-1]);

int c = cmpii(d->a[i], d->a[i-1]);

if (c < 0)

{

fv_a[i] = setloop(fv_a[i-1]);

d->a[i] = setloop(d->a[i-1]);

c = 0;

}

if (c == 0 && fv_fl == 2)

if (c == 0 && d->fl == 2)

fv_a[i] = incloop(fv_a[i]);

d->a[i] = incloop(d->a[i]);

}

if (i+1 == fv_n)

if (i+1 == d->n)

while (gcmp(fv_a[i], fv_M[i]) <= 0)

while (gcmp(d->a[i], d->M[i]) <= 0)

{

long av = avma; (void)lisseq(fv_ch); avma = av;

gpmem_t av = avma; (void)lisseq(d->ch); avma = av;

if (loop_break()) { fv_n = 0; return; }

if (loop_break()) { d->n = 0; return; }

fv_a[i] = incloop(fv_a[i]);

d->a[i] = incloop(d->a[i]);

}

else

while (gcmp(fv_a[i], fv_M[i]) <= 0)

while (gcmp(d->a[i], d->M[i]) <= 0)

{

long av = avma; fvloop_i(i+1); avma = av;

gpmem_t av = avma; fvloop_i(i+1, d); avma = av;

if (!fv_n) return;

if (!d->n) return;

fv_a[i] = incloop(fv_a[i]);

d->a[i] = incloop(d->a[i]);

}

void

forvec(entree *ep, GEN x, char *c, long flag)

{

long i, av = avma, tx = typ(x), n = fv_n, fl = fv_fl;

gpmem_t av = avma;

GEN *a = fv_a, *M = fv_M, *m = fv_m;

long tx = typ(x);

char *ch = fv_ch;

fvdat D, *d = &D;

void (*fv_fun)(long) = fvloop_i;

if (!is_vec_t(tx)) err(talker,"not a vector in forvec");

if (fl<0 || fl>2) err(flagerr);

if (flag<0 || flag>2) err(flagerr);

fv_n = lg(x);

d->n = lg(x);

fv_ch = c;

d->ch = c;

fv_fl = flag;

d->a = (GEN*)cgetg(d->n,t_VEC); push_val(ep, (GEN)d->a);

fv_a = (GEN*)cgetg(fv_n,t_VEC); push_val(ep, (GEN)fv_a);

if (d->n == 1) (void)lisseq(d->ch);

fv_m = (GEN*)cgetg(fv_n,t_VEC);

fv_M = (GEN*)cgetg(fv_n,t_VEC);

if (fv_n == 1) lisseq(fv_ch);

else

{

for (i=1; i<fv_n; i++)

long i, t = t_INT;

d->fl = flag;

d->m = (GEN*)cgetg(d->n,t_VEC);

d->M = (GEN*)cgetg(d->n,t_VEC);

for (i=1; i<d->n; i++)

{

GEN *c = (GEN*) x[i];

GEN *e = (GEN*) x[i];

tx = typ(c);

tx = typ(e);

if (! is_vec_t(tx) || lg(c)!=3)

if (! is_vec_t(tx) || lg(e)!=3)

err(talker,"not a vector of two-component vectors in forvec");

if (gcmp(c[1],c[2]) > 0) fv_n = 0;

if (gcmp(e[1],e[2]) > 0) d->n = 0;

/* in case x is an ep->value and c destroys it, we have to copy */

if (typ(e[1]) != t_INT) t = t_REAL;

if (typ(c[1]) != t_INT) fv_fun = fvloop;

/* in case x is an ep->value and lisexpr(d->ch) kills it, have to copy */

fv_m[i] = gcopy(c[1]);

d->m[i] = gcopy(e[1]);

fv_M[i] = gcopy(c[2]);

d->M[i] = gcopy(e[2]);

}

fv_fun(1);

if (t == t_INT) fvloop_i(1, d); else fvloop(1, d);

}

pop_val(ep);

pop_val(ep); avma = av;

fv_n = n;

fv_ch = ch;

fv_fl = fl;

fv_a = a;

fv_m = m;

fv_M = M; avma = av;

}

/********************************************************************/

Line 314 forvec(entree *ep, GEN x, char *c, long flag)

Line 315 forvec(entree *ep, GEN x, char *c, long flag)

GEN

somme(entree *ep, GEN a, GEN b, char *ch, GEN x)

{

long av,av0 = avma, lim;

gpmem_t av, av0 = avma, lim;

GEN p1;

if (typ(a) != t_INT) err(talker,"non integral index in sum");

if (!x) x = gzero;

if (gcmp(b,a) < 0) return gcopy(x);

b = gfloor(b);

a = setloop(a);

av=avma; lim = stack_lim(av,1);

push_val(ep, a);

Line 343 somme(entree *ep, GEN a, GEN b, char *ch, GEN x)

Line 344 somme(entree *ep, GEN a, GEN b, char *ch, GEN x)

GEN

suminf(entree *ep, GEN a, char *ch, long prec)

{

long fl,G,tetpil, av0 = avma, av,lim;

long fl, G;

gpmem_t tetpil, av0 = avma, av, lim;

GEN p1,x = realun(prec);

if (typ(a) != t_INT) err(talker,"non integral index in suminf");

Line 373 suminf(entree *ep, GEN a, char *ch, long prec)

Line 375 suminf(entree *ep, GEN a, char *ch, long prec)

GEN

divsum(GEN num, entree *ep, char *ch)

{

ulong av = avma;

gpmem_t av = avma;

GEN z, y = gzero, t = divisors(num);

long i, l = lg(t);

push_val(ep, NULL);

for (i=1; i<l; i++)

{

ep->value = (void*) t[i];

z = lisseq(ch);

if (did_break()) err(breaker,"divsum");

y = gadd(y, z);

Line 397 divsum(GEN num, entree *ep, char *ch)

Line 399 divsum(GEN num, entree *ep, char *ch)

GEN

produit(entree *ep, GEN a, GEN b, char *ch, GEN x)

{

long av,av0 = avma, lim;

gpmem_t av, av0 = avma, lim;

GEN p1;

if (typ(a) != t_INT) err(talker,"non integral index in sum");

if (!x) x = gun;

if (gcmp(b,a) < 0) return gcopy(x);

b = gfloor(b);

a = setloop(a);

av=avma; lim = stack_lim(av,1);

push_val(ep, a);

Line 438 prodinf0(entree *ep, GEN a, char *ch, long flag, long

Line 440 prodinf0(entree *ep, GEN a, char *ch, long flag, long

GEN

prodinf(entree *ep, GEN a, char *ch, long prec)

{

ulong av0 = avma, av,lim;

gpmem_t av0 = avma, av, lim;

long fl,G;

GEN p1,x = realun(prec);

Line 465 prodinf(entree *ep, GEN a, char *ch, long prec)

Line 467 prodinf(entree *ep, GEN a, char *ch, long prec)

GEN

prodinf1(entree *ep, GEN a, char *ch, long prec)

{

ulong av0 = avma, av,lim;

gpmem_t av0 = avma, av, lim;

long fl,G;

GEN p1,p2,x = realun(prec);

Line 492 prodinf1(entree *ep, GEN a, char *ch, long prec)

Line 494 prodinf1(entree *ep, GEN a, char *ch, long prec)

GEN

prodeuler(entree *ep, GEN ga, GEN gb, char *ch, long prec)

{

long prime[] = {evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3), 0};

long prime[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 0};

long a,b;

ulong av,av0 = avma, lim;

gpmem_t av, av0 = avma, lim;

GEN p1,x = realun(prec);

byteptr p;

av = avma;

p = prime_loop_init(ga,gb, &a,&b, prime);

if (!p) { avma = av; return x; }

Line 530 prodeuler(entree *ep, GEN ga, GEN gb, char *ch, long p

Line 532 prodeuler(entree *ep, GEN ga, GEN gb, char *ch, long p

GEN

direulerall(entree *ep, GEN ga, GEN gb, char *ch, GEN c)

{

long prime[] = {evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3), 0};

long prime[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 0};

ulong av0 = avma,av, lim = (av0+bot)>>1;

gpmem_t av0 = avma, av, lim = stack_lim(av0, 1);

long p,n,i,j,k,tx,lx,a,b;

GEN x,y,s,polnum,polden;

byteptr d;

d = prime_loop_init(ga,gb, &a,&b, prime);

n = c? itos(c): b;

if (!d || b < 2 || n <= 0) { x=cgetg(2,t_VEC); x[1]=un; return x; }

if (!d || b < 2 || n <= 0) return _vec(gun);

if (n < b) b = n;

push_val(ep, prime);

Line 636 vecteur(GEN nmax, entree *ep, char *ch)

Line 638 vecteur(GEN nmax, entree *ep, char *ch)

{

GEN y,p1;

long i,m;

long c[]={evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3), 0};

long c[]={evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 0};

if (typ(nmax)!=t_INT || signe(nmax) < 0)

err(talker,"bad number of components in vector");

Line 657 vecteur(GEN nmax, entree *ep, char *ch)

Line 659 vecteur(GEN nmax, entree *ep, char *ch)

}

GEN

vecteursmall(GEN nmax, entree *ep, char *ch)

{

GEN y,p1;

long i,m;

long c[]={evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 0};

if (typ(nmax)!=t_INT || signe(nmax) < 0)

err(talker,"bad number of components in vector");

m=itos(nmax); y=cgetg(m+1,t_VECSMALL);

if (!ep || !ch)

{

for (i=1; i<=m; i++) y[i] = 0;

return y;

}

push_val(ep, c);

for (i=1; i<=m; i++)

{

c[2]=i; p1 = lisseq(ch);

if (did_break()) err(breaker,"vector");

y[i] = itos(p1);

}

pop_val(ep); return y;

}

GEN

vvecteur(GEN nmax, entree *ep, char *ch)

{

GEN y = vecteur(nmax,ep,ch);

Line 668 matrice(GEN nlig, GEN ncol,entree *ep1, entree *ep2, c

Line 696 matrice(GEN nlig, GEN ncol,entree *ep1, entree *ep2, c

{

GEN y,z,p1;

long i,j,m,n,s;

long c1[]={evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3), 1};

long c1[]={evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 1};

long c2[]={evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3), 1};

long c2[]={evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3), 1};

s=signe(ncol);

if (ep1 == ep2) err(talker, "identical index variables in matrix");

if (ep1 == ep2 && ep1) err(talker, "identical index variables in matrix");

if (typ(ncol)!=t_INT || s<0) err(talker,"bad number of columns in matrix");

if (!s) return cgetg(1,t_MAT);

Line 712 matrice(GEN nlig, GEN ncol,entree *ep1, entree *ep2, c

Line 740 matrice(GEN nlig, GEN ncol,entree *ep1, entree *ep2, c

/********************************************************************/

/** **/

/** SOMMATION DE SERIES **/

/** SUMMING SERIES **/

/** **/

/********************************************************************/

Line 743 sumpos0(entree *ep, GEN a, char *ch, long flag, long p

Line 771 sumpos0(entree *ep, GEN a, char *ch, long flag, long p

GEN

sumalt(entree *ep, GEN a, char *ch, long prec)

{

long av = avma, tetpil,k,N;

long k, N;

gpmem_t av = avma, tetpil;

GEN s,az,c,x,e1,d;

if (typ(a) != t_INT) err(talker,"non integral index in sumalt");

push_val(ep, a);

e1=addsr(3,gsqrt(stoi(8),prec));

e1 = addsr(3,gsqrt(stoi(8),prec));

N=(long)(0.4*(bit_accuracy(prec) + 7));

N = (long)(0.4*(bit_accuracy(prec) + 7));

d=gpuigs(e1,N); d=shiftr(addrr(d,divsr(1,d)),-1);

d = gpowgs(e1,N);

az=negi(gun); c=d; s=gzero;

d = shiftr(addrr(d, ginv(d)),-1);

az = negi(gun); c = d;

s = gzero;

for (k=0; ; k++)

{

x = lisexpr(ch); if (did_break()) err(breaker,"sumalt");

Line 761 sumalt(entree *ep, GEN a, char *ch, long prec)

Line 792 sumalt(entree *ep, GEN a, char *ch, long prec)

if (k==N-1) break;

a = addsi(1,a); ep->value = (void*) a;

}

tetpil=avma; pop_val(ep);

tetpil = avma; pop_val(ep);

return gerepile(av,tetpil,gdiv(s,d));

}

GEN

sumalt2(entree *ep, GEN a, char *ch, long prec)

{

long av = avma, tetpil,k,N;

long k, N;

gpmem_t av = avma, tetpil;

GEN x,s,dn,pol;

if (typ(a) != t_INT) err(talker,"non integral index in sumalt");

push_val(ep, a);

N=(long)(0.31*(bit_accuracy(prec) + 5));

N = (long)(0.31*(bit_accuracy(prec) + 5));

s=gzero; pol=polzagreel(N,N>>1,prec+1); dn=poleval(pol,gun);

pol = polzagreel(N,N>>1,prec+1);

pol[2]=lsub((GEN)pol[2],dn); pol=gdiv(pol,gsub(polx[0],gun));

dn = poleval(pol,gun);

N=lgef(pol)-2;

pol[2] = lsub((GEN)pol[2],dn);

for (k=0; k<N; k++)

pol = gdiv(pol, gsub(polx[0],gun));

N = degpol(pol);

s = gzero;

for (k=0; k<=N; k++)

{

x = lisexpr(ch); if (did_break()) err(breaker,"sumalt2");

s=gadd(s,gmul(x,(GEN)pol[k+2]));

s = gadd(s, gmul(x,(GEN)pol[k+2]));

if (k==N-1) break;

if (k == N) break;

a = addsi(1,a); ep->value = (void*) a;

}

tetpil=avma; pop_val(ep);

tetpil = avma; pop_val(ep);

return gerepile(av,tetpil,gdiv(s,dn));

return gerepile(av,tetpil, gdiv(s,dn));

}

GEN

sumpos(entree *ep, GEN a, char *ch, long prec)

{

long av = avma,tetpil,k,kk,N,G;

long k, kk, N, G;

gpmem_t av = avma, tetpil;

GEN p1,r,q1,reel,s,az,c,x,e1,d, *stock;

if (typ(a) != t_INT) err(talker,"non integral index in sumpos");

push_val(ep, NULL);

a=subis(a,1); reel=cgetr(prec);

a = subis(a,1); reel = cgetr(prec);

e1=addsr(3,gsqrt(stoi(8),prec));

e1 = addsr(3,gsqrt(stoi(8),prec));

N=(long)(0.4*(bit_accuracy(prec) + 7));

N = (long)(0.4*(bit_accuracy(prec) + 7));

d=gpuigs(e1,N); d=shiftr(addrr(d,divsr(1,d)),-1);

d = gpowgs(e1,N); d = shiftr(addrr(d, ginv(d)),-1);

az=negi(gun); c=d; s=gzero;

az = negi(gun); c = d; s = gzero;

G = -bit_accuracy(prec) - 5;

stock=(GEN*)new_chunk(N+1); for (k=1; k<=N; k++) stock[k]=NULL;

stock = (GEN*)new_chunk(N+1); for (k=1; k<=N; k++) stock[k] = NULL;

for (k=0; k<N; k++)

{

if (k&1 && stock[k]) x=stock[k];

if (odd(k) && stock[k]) x = stock[k];

else

{

x=gzero; r=stoi(2*k+2);

x = gzero; r = stoi(2*k+2);

for(kk=0;;kk++)

{

long ex;

Line 818 sumpos(entree *ep, GEN a, char *ch, long prec)

Line 855 sumpos(entree *ep, GEN a, char *ch, long prec)

p1 = lisexpr(ch); if (did_break()) err(breaker,"sumpos");

gaffect(p1,reel); ex = expo(reel)+kk; setexpo(reel,ex);

x = gadd(x,reel); if (kk && ex < G) break;

r=shifti(r,1);

r = shifti(r,1);

}

if (2*k<N) stock[2*k+1]=x;

if (2*k < N) stock[2*k+1] = x;

q1 = addsi(k+1,a); ep->value = (void*) q1;

p1 = lisexpr(ch); if (did_break()) err(breaker,"sumpos");

gaffect(p1,reel); x = gadd(reel, gmul2n(x,1));

}

c=gadd(az,c); p1 = k&1? gneg_i(c): c;

c = gadd(az,c); p1 = k&1? gneg_i(c): c;

s = gadd(s,gmul(x,p1));

az = divii(mulii(mulss(N-k,N+k),shifti(az,1)),mulss(k+1,k+k+1));

}

Line 836 sumpos(entree *ep, GEN a, char *ch, long prec)

Line 873 sumpos(entree *ep, GEN a, char *ch, long prec)

GEN

sumpos2(entree *ep, GEN a, char *ch, long prec)

{

long av = avma,tetpil,k,kk,N,G;

long k, kk, N, G;

gpmem_t av = avma, tetpil;

GEN p1,r,q1,reel,s,pol,dn,x, *stock;

if (typ(a) != t_INT) err(talker,"non integral index in sumpos2");

push_val(ep, a);

a=subis(a,1); reel=cgetr(prec);

a = subis(a,1); reel = cgetr(prec);

N=(long)(0.31*(bit_accuracy(prec) + 5));

N = (long)(0.31*(bit_accuracy(prec) + 5));

G = -bit_accuracy(prec) - 5;

stock=(GEN*)new_chunk(N+1); for (k=1; k<=N; k++) stock[k]=NULL;

stock = (GEN*)new_chunk(N+1); for (k=1; k<=N; k++) stock[k] = NULL;

for (k=1; k<=N; k++)

{

if (odd(k) || !stock[k])

{

x=gzero; r=stoi(2*k);

x = gzero; r = stoi(2*k);

for(kk=0;;kk++)

{

long ex;

Line 858 sumpos2(entree *ep, GEN a, char *ch, long prec)

Line 897 sumpos2(entree *ep, GEN a, char *ch, long prec)

p1 = lisexpr(ch); if (did_break()) err(breaker,"sumpos2");

gaffect(p1,reel); ex = expo(reel)+kk; setexpo(reel,ex);

x=gadd(x,reel); if (kk && ex < G) break;

r=shifti(r,1);

r = shifti(r,1);

}

if (2*k-1<N) stock[2*k]=x;

if (2*k-1 < N) stock[2*k] = x;

q1 = addsi(k,a); ep->value = (void*) q1;

p1 = lisexpr(ch); if (did_break()) err(breaker,"sumpos2");

gaffect(p1,reel); stock[k] = gadd(reel, gmul2n(x,1));

}

pop_val(ep); s = gzero;

pol=polzagreel(N,N>>1,prec+1); dn=poleval(pol,gun);

pol = polzagreel(N,N>>1,prec+1);

pol[2]=lsub((GEN)pol[2],dn); pol=gdiv(pol,gsub(gun,polx[0]));

dn = poleval(pol,gun);

pol[2] = lsub((GEN)pol[2],dn);

pol = gdiv(pol, gsub(gun,polx[0]));

for (k=1; k<=lgef(pol)-2; k++)

{

p1 = gmul((GEN)pol[k+1],stock[k]);

if (k&1) p1 = gneg_i(p1);

if (odd(k)) p1 = gneg_i(p1);

s = gadd(s,p1);

}

tetpil=avma; return gerepile(av,tetpil,gdiv(s,dn));

Line 883 sumpos2(entree *ep, GEN a, char *ch, long prec)

Line 924 sumpos2(entree *ep, GEN a, char *ch, long prec)

/** NUMERICAL INTEGRATION (Romberg) **/

/** **/

/********************************************************************/

GEN

typedef struct {

intnum0(entree *ep, GEN a, GEN b, char *ch, long flag, long prec)

entree *ep;

{

char *ch;

switch(flag)

} exprdat;

{

case 0: return qromb (ep,a,b,ch,prec);

case 1: return rombint(ep,a,b,ch,prec);

case 2: return qromi (ep,a,b,ch,prec);

case 3: return qromo (ep,a,b,ch,prec);

default: err(flagerr);

}

return NULL; /* not reached */

}

#define JMAX 25

typedef struct {

#define JMAXP JMAX+3

GEN (*f)(void *, GEN);

#define KLOC 4

void *E;

} invfun;

/* we need to make a copy in any case, cf forpari */

static GEN

Line 909 fix(GEN a, long prec)

Line 942 fix(GEN a, long prec)

gaffect(a,p); return p;

}

extern GEN polint_i(GEN xa, GEN ya, GEN x, long n, GEN *ptdy);

#if 0

GEN

qromb(entree *ep, GEN a, GEN b, char *ch, long prec)

int_loop(entree *ep, char *ch)

{

GEN ss,dss,s,h,p1,p2,qlint,del,x,sum;

gpmem_t av = avma;

long av = avma, av1, tetpil,j,j1,j2,lim,it,sig;

intstruct T;

T.in = NULL;

for(;;)

{

GEN x = Next(&T);

if (!x) return gerepileupto(av, T.out);

ep->value = (void*)x;

T.in = lisexpr(ch);

}

#endif

/* f(x) */

static GEN

_gp_eval(void *dat, GEN x)

{

exprdat *E = (exprdat*)dat;

E->ep->value = x;

return lisexpr(E->ch);

}

/* 1/x^2 f(1/x) */

static GEN

_invf(void *dat, GEN x)

{

invfun *S = (invfun*)dat;

GEN y = ginv(x);

return gmul(S->f(S->E, y), gsqr(y));

}

#define swap(a,b) { GEN _x = a; a = b; b = _x; }

static GEN

interp(GEN h, GEN s, long j, long lim, long KLOC)

{

gpmem_t av = avma;

long e1,e2;

GEN dss, ss = polint_i(h+j-KLOC,s+j-KLOC,gzero,KLOC+1,&dss);

e1 = gexpo(ss);

e2 = gexpo(dss);

if (e1-e2 <= lim && e1 >= -lim) { avma = av; return NULL; }

if (gcmp0(gimag(ss))) ss = greal(ss);

return ss;

}

static GEN

qrom3(void *dat, GEN (*eval)(void *, GEN), GEN a, GEN b, long prec)

{

const long JMAX = 25, KLOC = 4;

GEN ss,s,h,p1,p2,qlint,del,x,sum;

long j, j1, it, sig;

gpmem_t av;

a = fix(a,prec);

b = fix(b,prec);

qlint=subrr(b,a); sig=signe(qlint);

qlint = subrr(b,a); sig = signe(qlint);

if (!sig) { avma=av; return gzero; }

if (!sig) return gzero;

if (sig<0) { setsigne(qlint,1); s=a; a=b; b=s; }

if (sig < 0) { setsigne(qlint,1); swap(a,b); }

s=new_chunk(JMAXP);

s = new_chunk(JMAX+KLOC-1);

h=new_chunk(JMAXP);

h = new_chunk(JMAX+KLOC-1);

h[0] = (long)realun(prec);

push_val(ep, a);

p1 = eval(dat, a); if (p1 == a) p1 = rcopy(p1);

p1=lisexpr(ch); if (p1 == a) p1=rcopy(p1);

p2 = eval(dat, b);

ep->value = (void*)b;

s[0] = lmul2n(gmul(qlint,gadd(p1,p2)),-1);

p2=lisexpr(ch);

for (it=1,j=1; j<JMAX; j++, it<<=1)

s[0]=lmul2n(gmul(qlint,gadd(p1,p2)),-1);

for (it=1,j=1; j<JMAX; j++,it=it<<1)

{

h[j]=lshiftr((GEN)h[j-1],-2);

h[j] = lshiftr((GEN)h[j-1],-2);

av1=avma; del=divrs(qlint,it); x=addrr(a,shiftr(del,-1));

av = avma; del = divrs(qlint,it);

for (sum=gzero,j1=1; j1<=it; j1++,x=addrr(x,del))

x = addrr(a, shiftr(del,-1));

{

for (sum = gzero, j1 = 1; j1 <= it; j1++, x = addrr(x,del))

ep->value = (void*) x; sum=gadd(sum,lisexpr(ch));

sum = gadd(sum, eval(dat, x));

}

sum = gmul(sum,del); p1 = gadd((GEN)s[j-1], sum);

sum = gmul(sum,del); p1 = gadd((GEN)s[j-1],sum);

s[j] = lpileupto(av, gmul2n(p1,-1));

tetpil = avma;

s[j]=lpile(av1,tetpil,gmul2n(p1,-1));

if (j>=KLOC)

if (j >= KLOC && (ss = interp(h, s, j, bit_accuracy(prec)-j-6, KLOC)))

{

return gmulsg(sig,ss);

tetpil=avma; ss=polint_i(h+j-KLOC,s+j-KLOC,gzero,KLOC+1,&dss);

j1=gexpo(ss); j2=gexpo(dss); lim=bit_accuracy(prec)-j-6;

if (j1-j2 > lim || j1 < -lim)

{

pop_val(ep); if (gcmp0(gimag(ss))) ss=greal(ss);

tetpil=avma; return gerepile(av,tetpil,gmulsg(sig,ss));

}

avma=tetpil;

}

err(intger2); return NULL; /* not reached */

return NULL;

}

GEN

static GEN

qromo(entree *ep, GEN a, GEN b, char *ch, long prec)

qrom2(void *dat, GEN (*eval)(void *, GEN), GEN a, GEN b, long prec)

{

GEN ss,dss,s,h,p1,qlint,del,ddel,x,sum;

const long JMAX = 16, KLOC = 4;

long av = avma,av1,tetpil,j,j1,j2,lim,it,sig;

GEN ss,s,h,p1,qlint,del,ddel,x,sum;

long j, j1, it, sig;

gpmem_t av;

a = fix(a, prec);

b = fix(b, prec);

qlint=subrr(b,a); sig=signe(qlint);

qlint = subrr(b,a); sig = signe(qlint);

if (!sig) { avma=av; return gzero; }

if (!sig) return gzero;

if (sig<0) { setsigne(qlint,1); s=a; a=b; b=s; }

if (sig < 0) { setsigne(qlint,1); swap(a,b); }

s=new_chunk(JMAXP);

s = new_chunk(JMAX+KLOC-1);

h=new_chunk(JMAXP);

h = new_chunk(JMAX+KLOC-1);

h[0] = (long)realun(prec);

p1 = shiftr(addrr(a,b),-1); push_val(ep, p1);

p1 = shiftr(addrr(a,b),-1);

p1=lisexpr(ch); s[0]=lmul(qlint,p1);

s[0] = lmul(qlint, eval(dat, p1));

for (it=1, j=1; j<JMAX; j++, it*=3)

{

h[j] = ldivrs((GEN)h[j-1],9);

av1=avma; del=divrs(qlint,3*it); ddel=shiftr(del,1);

av = avma; del = divrs(qlint,3*it); ddel = shiftr(del,1);

x=addrr(a,shiftr(del,-1)); sum=gzero;

x = addrr(a, shiftr(del,-1));

for (j1=1; j1<=it; j1++)

for (sum = gzero, j1 = 1; j1 <= it; j1++)

{

ep->value = (void*)x; sum=gadd(sum,lisexpr(ch)); x=addrr(x,ddel);

sum = gadd(sum, eval(dat, x)); x = addrr(x,ddel);

ep->value = (void*)x; sum=gadd(sum,lisexpr(ch)); x=addrr(x,del);

sum = gadd(sum, eval(dat, x)); x = addrr(x,del);

}

sum = gmul(sum,del); p1 = gdivgs((GEN)s[j-1],3);

tetpil = avma; s[j]=lpile(av1,tetpil,gadd(p1,sum));

s[j] = lpileupto(av, gadd(p1,sum));

if (j>=KLOC)

if (j >= KLOC && (ss = interp(h, s, j, bit_accuracy(prec)-(3*j/2)-6, KLOC)))

{

return gmulsg(sig, ss);

tetpil=avma; ss=polint_i(h+j-KLOC,s+j-KLOC,gzero,KLOC+1,&dss);

j1=gexpo(ss); j2=gexpo(dss); lim=bit_accuracy(prec)-(3*j/2)-6;

if ( j1-j2 > lim || j1 < -lim)

{

pop_val(ep); if (gcmp0(gimag(ss))) ss=greal(ss);

tetpil=avma; return gerepile(av,tetpil,gmulsg(sig,ss));

}

avma=tetpil;

}

err(intger2); return NULL; /* not reached */

return NULL;

}

#undef JMAX

/* integrate after change of variables x --> 1/x */

#undef JMAXP

static GEN

#define JMAX 16

qromi(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, long prec)

#define JMAXP JMAX+3

{

GEN A = ginv(b), B = ginv(a);

invfun S;

S.f = eval;

S.E = E; return qrom2(&S, &_invf, A, B, prec);

}

GEN

/* a < b, assume b "small" (< 100 say) */

qromi(entree *ep, GEN a, GEN b, char *ch, long prec)

static GEN

rom_bsmall(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, long prec)

{

GEN ss,dss,s,h,q1,p1,p,qlint,del,ddel,x,sum;

if (gcmpgs(a,-100) >= 0) return qrom2(E,eval,a,b,prec);

long av = avma, av1,tetpil,j,j1,j2,lim,it,sig;

if (b == gun || gcmpgs(b, -1) >= 0)

{ /* a < -100, b >= -1 */

GEN _1 = negi(gun); /* split at -1 */

return gadd(qromi(E,eval,a,_1,prec),

qrom2(E,eval,_1,b,prec));

}

/* a < -100, b < -1 */

return qromi(E,eval,a,b,prec);

}

p=cgetr(prec); gaffect(ginv(a),p); a=p;

static GEN

p=cgetr(prec); gaffect(ginv(b),p); b=p;

rombint(void *E, GEN (*eval)(void*, GEN), GEN a, GEN b, long prec)

qlint=subrr(b,a); sig= -signe(qlint);

{

if (!sig) { avma=av; return gzero; }

long l = gcmp(b,a);

if (sig>0) { setsigne(qlint,1); s=a; a=b; b=s; }

GEN z;

s=new_chunk(JMAXP);

if (!l) return gzero;

h=new_chunk(JMAXP);

if (l < 0) swap(a,b);

h[0] = (long)realun(prec);

if (gcmpgs(b,100) >= 0)

x=divsr(2,addrr(a,b)); push_val(ep, x);

p1=gmul(lisexpr(ch),mulrr(x,x));

s[0]=lmul(qlint,p1);

for (it=1,j=1; j<JMAX; j++, it*=3)

{

h[j] = ldivrs((GEN)h[j-1],9);

if (gcmpgs(a,1) >= 0)

av1=avma; del=divrs(qlint,3*it); ddel=shiftr(del,1);

z = qromi(E,eval,a,b,prec);

x=addrr(a,shiftr(del,-1)); sum=gzero;

else /* split at 1 */

for (j1=1; j1<=it; j1++)

z = gadd(rom_bsmall(E,eval,a,gun,prec), qromi(E,eval,gun,b,prec));

{

q1 = ginv(x); ep->value = (void*)q1;

p1=gmul(lisexpr(ch), gsqr(q1));

sum=gadd(sum,p1); x=addrr(x,ddel);

q1 = ginv(x); ep->value = (void*)q1;

p1=gmul(lisexpr(ch), gsqr(q1));

sum=gadd(sum,p1); x=addrr(x,del);

}

sum = gmul(sum,del); p1 = gdivgs((GEN)s[j-1],3);

tetpil=avma;

s[j]=lpile(av1,tetpil,gadd(p1,sum));

if (j>=KLOC)

{

tetpil=avma; ss=polint_i(h+j-KLOC,s+j-KLOC,gzero,KLOC+1,&dss);

j1=gexpo(ss); j2=gexpo(dss); lim=bit_accuracy(prec)-(3*j/2)-6;

if (j1-j2 > lim || j1 < -lim)

{

pop_val(ep); if (gcmp0(gimag(ss))) ss=greal(ss);

tetpil=avma; return gerepile(av,tetpil,gmulsg(sig,ss));

}

err(intger2); return NULL; /* not reached */

else

z = rom_bsmall(E,eval,a,b,prec);

if (l < 0) z = gneg(z);

return z;

}

GEN

rombint(entree *ep, GEN a, GEN b, char *ch, long prec)

intnum(void *E, GEN (*eval)(void*,GEN), GEN a, GEN b, long flag, long prec)

{

GEN aa,bb,mun,p1,p2,p3;

gpmem_t av = avma;

long l,av,tetpil;

GEN z;

switch(flag)

l=gcmp(b,a); if (!l) return gzero;

if (l<0) { bb=a; aa=b; } else { bb=b; aa=a; }

av=avma; mun = negi(gun);

if (gcmpgs(bb,100)>=0)

{

if (gcmpgs(aa,1)>=0) return qromi(ep,a,b,ch,prec);

case 0: z = qrom3 (E,eval,a,b,prec); break;

p1=qromi(ep,gun,bb,ch,prec);

case 1: z = rombint(E,eval,a,b,prec); break;

if (gcmpgs(aa,-100)>=0)

case 2: z = qromi (E,eval,a,b,prec); break;

{

case 3: z = qrom2 (E,eval,a,b,prec); break;

p2=qromo(ep,aa,gun,ch,prec); tetpil=avma;

default: err(flagerr); z = NULL;

return gerepile(av,tetpil,gmulsg(l,gadd(p1,p2)));

}

p2=qromo(ep,mun,gun,ch,prec); p3=gadd(p2,qromi(ep,aa,mun,ch,prec));

tetpil=avma; return gerepile(av,tetpil,gmulsg(l,gadd(p1,p3)));

}

if (gcmpgs(aa,-100)>=0) return qromo(ep,a,b,ch,prec);

if (!z) err(intger2);

if (gcmpgs(bb,-1)>=0)

return gerepileupto(av, z);

{

p1=qromi(ep,aa,mun,ch,prec); p2=qromo(ep,mun,bb,ch,prec); tetpil=avma;

return gerepile(av,tetpil,gmulsg(l,gadd(p1,p2)));

}

return qromi(ep,a,b,ch,prec);

}

GEN

intnum0(entree *ep, GEN a, GEN b, char *ch, long flag, long prec)

{

exprdat E;

GEN z;

E.ch = ch;

E.ep = ep; push_val(ep, NULL);

z = intnum(&E, &_gp_eval, a,b, flag, prec);

pop_val(ep); return z;

}

/********************************************************************/

/** **/

/** ZAGIER POLYNOMIALS (for sumiter) **/

/** ZAGIER POLYNOMIALS (for sumalt) **/

/** **/

/********************************************************************/

GEN

polzag(long n, long m)

{

long d1,d,r,k,av=avma,tetpil;

gpmem_t av = avma;

GEN p1,p2,pol1,g,s;

long d1, d, r, k;

GEN A, B, Bx, g, s;

if (m>=n || m<0)

if (m >= n || m < 0)

err(talker,"first index must be greater than second in polzag");

d1=n-m; d=d1<<1; d1--; p1=gsub(gun,gmul2n(polx[0],1));

d1 = n-m; d = d1<<1; d1--;

pol1=gsub(gun,polx[0]); p2=gmul(polx[0],pol1); r=(m+1)>>1;

A = gsub(gun, gmul2n(polx[0],1)); /* 1 - 2x */

g=gzero;

B = gsub(gun, polx[0]); /* 1 - x */

Bx = gmul(polx[0],B); /* x - x^2 */

r = (m+1)>>1;

g = gzero;

for (k=0; k<=d1; k++)

{

s=binome(stoi(d),(k<<1)+1); if (k&1) s=negi(s);

s = binome(stoi(d), (k<<1)+1); if (k&1) s = negi(s);

g=gadd(g,gmul(s,gmul(gpuigs(polx[0],k),gpuigs(pol1,d1-k))));

g = gadd(g, gmul(s, gmul(gpowgs(polx[0],k), gpowgs(B,d1-k))));

}

g=gmul(g,gpuigs(p2,r));

g = gmul(g, gpowgs(Bx,r));

if (!(m&1)) g=gadd(gmul(p1,g),gmul2n(gmul(p2,derivpol(g)),1));

if (!(m&1)) g = gadd(gmul(A,g), gmul2n(gmul(Bx,derivpol(g)),1));

for (k=1; k<=r; k++)

{

g=derivpol(g); g=gadd(gmul(p1,g),gmul2n(gmul(p2,derivpol(g)),1));

g = derivpol(g);

g = gadd(gmul(A,g), gmul2n(gmul(Bx,derivpol(g)),1));

}

g = m ? gmul2n(g,(m-1)>>1):gmul2n(g,-1);

g = m? gmul2n(g,(m-1)>>1): gmul2n(g,-1);

s=mulsi(n-m,mpfact(m+1));

s = mulsi(n-m, mpfact(m+1));

tetpil=avma; return gerepile(av,tetpil,gdiv(g,s));

return gerepileupto(av, gdiv(g,s));

}

#ifdef _MSC_VER /* Bill Daly: work around a MSVC bug */

Line 1132 polzag(long n, long m)

Line 1186 polzag(long n, long m)

GEN

polzagreel(long n, long m, long prec)

{

long d1,d,r,j,k,k2,av=avma,tetpil;

long d1, d, r, j, k, k2;

GEN p2,pol1,g,h,v,b,gend,s;

gpmem_t av = avma;

GEN Bx,B,g,h,v,b,s;

if (m>=n || m<0)

if (m >= n || m < 0)

err(talker,"first index must be greater than second in polzag");

d1=n-m; d=d1<<1; d1--; pol1=gadd(gun,polx[0]);

d1 = n-m; d = d1<<1; d1--;

p2=gmul(polx[0],pol1); r=(m+1)>>1; gend=stoi(d);

B = gadd(gun,polx[0]);

v=cgetg(d1+2,t_VEC); g=cgetg(d1+2,t_VEC);

Bx = gmul(polx[0],B);

v[d1+1]=un; b=mulri(realun(prec),gend); g[d1+1]=(long)b;

r = (m+1)>>1;

v = cgetg(d1+2,t_VEC);

g = cgetg(d1+2,t_VEC);

v[d1+1] = un; b = stor(d, prec);

g[d1+1] = (long)b;

for (k=1; k<=d1; k++)

{

v[d1-k+1]=un;

v[d1-k+1] = un;

for (j=1; j<k; j++)

v[d1-k+j+1]=(long)addii((GEN)v[d1-k+j+1],(GEN)v[d1-k+j+2]);

v[d1-k+j+1] = laddii((GEN)v[d1-k+j+1], (GEN)v[d1-k+j+2]);

k2=k+k; b=divri(mulri(b,mulss(d-k2+1,d-k2)),mulss(k2,k2+1));

k2 = k+k; b = divri(mulri(b,mulss(d-k2+1,d-k2)), mulss(k2,k2+1));

for (j=1; j<=k; j++)

g[d1-k+j+1]=(long)mpadd((GEN)g[d1-k+j+1],mulri(b,(GEN)v[d1-k+j+1]));

g[d1-k+j+1] = lmpadd((GEN)g[d1-k+j+1], mulri(b,(GEN)v[d1-k+j+1]));

g[d1-k+1]=(long)b;

g[d1-k+1] = (long)b;

}

h=cgetg(d1+3,t_POL); h[1]=evalsigne(1)+evallgef(d1+3);

g = gmul(vec_to_pol(g,0), gpowgs(Bx,r));

for (k=0; k<=d1; k++) h[k+2]=g[k+1];

g=gmul(h,gpuigs(p2,r));

for (j=0; j<=r; j++)

{

if (j) g=derivpol(g);

if (j) g = derivpol(g);

if (j || !(m&1))

{

h=cgetg(n+3,t_POL); h[1]=evalsigne(1)+evallgef(n+3);

h = cgetg(n+3,t_POL);

h[2]=g[2];

h[1] = evalsigne(1)|evallgef(n+3);

h[2] = g[2];

for (k=1; k<n; k++)

h[k+2]=ladd(gmulsg(k+k+1,(GEN)g[k+2]),gmulsg(k<<1,(GEN)g[k+1]));

h[k+2] = ladd(gmulsg(k+k+1,(GEN)g[k+2]), gmulsg(k<<1,(GEN)g[k+1]));

h[n+2]=lmulsg(n<<1,(GEN)g[n+1]); g=h;

h[n+2] = lmulsg(n<<1, (GEN)g[n+1]);

g = h;

}

g = m ?gmul2n(g,(m-1)>>1):gmul2n(g,-1);

g = m? gmul2n(g,(m-1)>>1): gmul2n(g,-1);

s=mulsi(n-m,mpfact(m+1));

s = mulsi(n-m, mpfact(m+1));

tetpil=avma; return gerepile(av,tetpil,gdiv(g,s));

return gerepileupto(av, gdiv(g,s));

}

#ifdef _MSC_VER

#pragma optimize("g",on)

Line 1183 polzagreel(long n, long m, long prec)

Line 1242 polzagreel(long n, long m, long prec)

GEN

zbrent(entree *ep, GEN a, GEN b, char *ch, long prec)

{

long av = avma,sig,iter,itmax;

long sig, iter, itmax;

gpmem_t av = avma;

GEN c,d,e,tol,tol1,min1,min2,fa,fb,fc,p,q,r,s,xm;

a = fix(a,prec);

Line 1238 zbrent(entree *ep, GEN a, GEN b, char *ch, long prec)

Line 1298 zbrent(entree *ep, GEN a, GEN b, char *ch, long prec)

else { d = xm; e = d; } /* bound decreasing too slowly, use bisection */

a = b; fa = fb;

if (gcmp(gabs(d,0),tol1) > 0) b = gadd(b,d);

else

else if (gsigne(xm) > 0) b = addrr(b,tol1);

{

else b = subrr(b,tol1);

if (gsigne(xm) > 0)

b = addrr(b,tol1);

else

b = subrr(b,tol1);

}

ep->value = (void*)b; fb = lisexpr(ch);

}

if (iter > itmax) err(talker,"too many iterations in solve");

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>