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, |
|
Line 21 Foundation, Inc., 59 Temple Place - Suite 330, Boston, |
|
#include "pari.h" |
#include "pari.h" |
#include "anal.h" |
#include "anal.h" |
extern void changevalue_p(entree *ep, GEN x); |
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 |
void |
forpari(entree *ep, GEN a, GEN b, char *ch) |
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); |
b = gcopy(b); av=avma; lim = stack_lim(av,1); |
/* gcopy nedeed in case b gets overwritten in ch, as in |
/* 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); |
push_val(ep, a); |
while (gcmp(a,b) <= 0) |
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; |
if (loop_break()) break; |
a = (GEN) ep->value; a = gadd(a,gun); |
a = (GEN) ep->value; a = gadd(a,gun); |
if (low_stack(lim, stack_lim(av,1))) |
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 |
void |
forstep(entree *ep, GEN a, GEN b, GEN s, char *ch) |
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; |
GEN v = NULL; |
int (*cmp)(GEN,GEN); |
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; |
i = 0; |
while (cmp(a,b) <= 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 (loop_break()) break; |
if (v) |
if (v) |
{ |
{ |
Line 103 sinitp(long a, long c, byteptr *ptr) |
|
Line 105 sinitp(long a, long c, byteptr *ptr) |
|
byteptr p = *ptr; |
byteptr p = *ptr; |
if (a <= 0) a = 2; |
if (a <= 0) a = 2; |
if (maxprime() < (ulong)a) err(primer1); |
if (maxprime() < (ulong)a) err(primer1); |
while (a > c) c += *p++; |
while (a > c) |
|
NEXT_PRIME_VIADIFF(c,p); |
*ptr = p; return c; |
*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; |
*ptr = diffptr; |
prime[2] = sinitp(a, 0, ptr); |
prime[2] = sinitp(a, 0, ptr); |
} |
} |
changevalue_p(ep, prime); |
changevalue_p(ep, prime); |
} |
} |
|
|
static byteptr |
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 |
void |
forprime(entree *ep, GEN ga, GEN gb, char *ch) |
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; |
byteptr p; |
|
|
p = prime_loop_init(ga,gb, &a,&b, prime); |
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 |
void |
fordiv(GEN a, entree *ep, char *ch) |
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); |
GEN t = divisors(a); |
|
|
push_val(ep, NULL); l=lg(t); av2 = avma; |
push_val(ep, NULL); l=lg(t); av2 = avma; |
for (i=1; i<l; i++) |
for (i=1; i<l; i++) |
{ |
{ |
ep->value = (void*) t[i]; |
ep->value = (void*) t[i]; |
(void)lisseq(ch); if (loop_break()) break; |
(void)lisseq(ch); if (loop_break()) break; |
avma = av2; |
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 = 1: impose a1 <= ... <= an |
* fl = 2: 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 */ |
/* general case */ |
static void |
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) |
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 |
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 */ |
/* we run through integers */ |
static void |
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) |
if (c < 0) |
{ |
{ |
fv_a[i] = setloop(fv_a[i-1]); |
d->a[i] = setloop(d->a[i-1]); |
c = 0; |
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 |
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 |
void |
forvec(entree *ep, GEN x, char *c, long flag) |
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 (!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 |
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"); |
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 |
GEN |
somme(entree *ep, GEN a, GEN b, char *ch, GEN x) |
somme(entree *ep, GEN a, GEN b, char *ch, GEN x) |
{ |
{ |
long av,av0 = avma, lim; |
gpmem_t av, av0 = avma, lim; |
GEN p1; |
GEN p1; |
|
|
if (typ(a) != t_INT) err(talker,"non integral index in sum"); |
if (typ(a) != t_INT) err(talker,"non integral index in sum"); |
if (!x) x = gzero; |
if (!x) x = gzero; |
if (gcmp(b,a) < 0) return gcopy(x); |
if (gcmp(b,a) < 0) return gcopy(x); |
|
|
b = gfloor(b); |
b = gfloor(b); |
a = setloop(a); |
a = setloop(a); |
av=avma; lim = stack_lim(av,1); |
av=avma; lim = stack_lim(av,1); |
push_val(ep, a); |
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 |
GEN |
suminf(entree *ep, GEN a, char *ch, long prec) |
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); |
GEN p1,x = realun(prec); |
|
|
if (typ(a) != t_INT) err(talker,"non integral index in suminf"); |
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 |
GEN |
divsum(GEN num, entree *ep, char *ch) |
divsum(GEN num, entree *ep, char *ch) |
{ |
{ |
ulong av = avma; |
gpmem_t av = avma; |
GEN z, y = gzero, t = divisors(num); |
GEN z, y = gzero, t = divisors(num); |
long i, l = lg(t); |
long i, l = lg(t); |
|
|
push_val(ep, NULL); |
push_val(ep, NULL); |
for (i=1; i<l; i++) |
for (i=1; i<l; i++) |
{ |
{ |
ep->value = (void*) t[i]; |
ep->value = (void*) t[i]; |
z = lisseq(ch); |
z = lisseq(ch); |
if (did_break()) err(breaker,"divsum"); |
if (did_break()) err(breaker,"divsum"); |
y = gadd(y, z); |
y = gadd(y, z); |
Line 397 divsum(GEN num, entree *ep, char *ch) |
|
Line 399 divsum(GEN num, entree *ep, char *ch) |
|
GEN |
GEN |
produit(entree *ep, GEN a, GEN b, char *ch, GEN x) |
produit(entree *ep, GEN a, GEN b, char *ch, GEN x) |
{ |
{ |
long av,av0 = avma, lim; |
gpmem_t av, av0 = avma, lim; |
GEN p1; |
GEN p1; |
|
|
if (typ(a) != t_INT) err(talker,"non integral index in sum"); |
if (typ(a) != t_INT) err(talker,"non integral index in sum"); |
if (!x) x = gun; |
if (!x) x = gun; |
if (gcmp(b,a) < 0) return gcopy(x); |
if (gcmp(b,a) < 0) return gcopy(x); |
|
|
b = gfloor(b); |
b = gfloor(b); |
a = setloop(a); |
a = setloop(a); |
av=avma; lim = stack_lim(av,1); |
av=avma; lim = stack_lim(av,1); |
push_val(ep, a); |
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 |
GEN |
prodinf(entree *ep, GEN a, char *ch, long prec) |
prodinf(entree *ep, GEN a, char *ch, long prec) |
{ |
{ |
ulong av0 = avma, av,lim; |
gpmem_t av0 = avma, av, lim; |
long fl,G; |
long fl,G; |
GEN p1,x = realun(prec); |
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 |
GEN |
prodinf1(entree *ep, GEN a, char *ch, long prec) |
prodinf1(entree *ep, GEN a, char *ch, long prec) |
{ |
{ |
ulong av0 = avma, av,lim; |
gpmem_t av0 = avma, av, lim; |
long fl,G; |
long fl,G; |
GEN p1,p2,x = realun(prec); |
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 |
GEN |
prodeuler(entree *ep, GEN ga, GEN gb, char *ch, long prec) |
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; |
long a,b; |
ulong av,av0 = avma, lim; |
gpmem_t av, av0 = avma, lim; |
GEN p1,x = realun(prec); |
GEN p1,x = realun(prec); |
byteptr p; |
byteptr p; |
|
|
av = avma; |
av = avma; |
p = prime_loop_init(ga,gb, &a,&b, prime); |
p = prime_loop_init(ga,gb, &a,&b, prime); |
if (!p) { avma = av; return x; } |
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 |
GEN |
direulerall(entree *ep, GEN ga, GEN gb, char *ch, GEN c) |
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; |
long p,n,i,j,k,tx,lx,a,b; |
GEN x,y,s,polnum,polden; |
GEN x,y,s,polnum,polden; |
byteptr d; |
byteptr d; |
|
|
d = prime_loop_init(ga,gb, &a,&b, prime); |
d = prime_loop_init(ga,gb, &a,&b, prime); |
n = c? itos(c): b; |
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; |
if (n < b) b = n; |
push_val(ep, prime); |
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; |
GEN y,p1; |
long i,m; |
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) |
if (typ(nmax)!=t_INT || signe(nmax) < 0) |
err(talker,"bad number of components in vector"); |
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 |
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) |
vvecteur(GEN nmax, entree *ep, char *ch) |
{ |
{ |
GEN y = vecteur(nmax,ep,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; |
GEN y,z,p1; |
long i,j,m,n,s; |
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); |
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 (typ(ncol)!=t_INT || s<0) err(talker,"bad number of columns in matrix"); |
if (!s) return cgetg(1,t_MAT); |
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 |
GEN |
sumalt(entree *ep, GEN a, char *ch, long prec) |
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; |
GEN s,az,c,x,e1,d; |
|
|
if (typ(a) != t_INT) err(talker,"non integral index in sumalt"); |
if (typ(a) != t_INT) err(talker,"non integral index in sumalt"); |
|
|
push_val(ep, a); |
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++) |
for (k=0; ; k++) |
{ |
{ |
x = lisexpr(ch); if (did_break()) err(breaker,"sumalt"); |
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; |
if (k==N-1) break; |
a = addsi(1,a); ep->value = (void*) a; |
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)); |
return gerepile(av,tetpil,gdiv(s,d)); |
} |
} |
|
|
GEN |
GEN |
sumalt2(entree *ep, GEN a, char *ch, long prec) |
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; |
GEN x,s,dn,pol; |
|
|
if (typ(a) != t_INT) err(talker,"non integral index in sumalt"); |
if (typ(a) != t_INT) err(talker,"non integral index in sumalt"); |
|
|
push_val(ep, a); |
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"); |
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; |
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 |
GEN |
sumpos(entree *ep, GEN a, char *ch, long prec) |
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; |
GEN p1,r,q1,reel,s,az,c,x,e1,d, *stock; |
|
|
if (typ(a) != t_INT) err(talker,"non integral index in sumpos"); |
if (typ(a) != t_INT) err(talker,"non integral index in sumpos"); |
|
|
push_val(ep, NULL); |
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; |
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++) |
for (k=0; k<N; k++) |
{ |
{ |
if (k&1 && stock[k]) x=stock[k]; |
if (odd(k) && stock[k]) x = stock[k]; |
else |
else |
{ |
{ |
x=gzero; r=stoi(2*k+2); |
x = gzero; r = stoi(2*k+2); |
for(kk=0;;kk++) |
for(kk=0;;kk++) |
{ |
{ |
long ex; |
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"); |
p1 = lisexpr(ch); if (did_break()) err(breaker,"sumpos"); |
gaffect(p1,reel); ex = expo(reel)+kk; setexpo(reel,ex); |
gaffect(p1,reel); ex = expo(reel)+kk; setexpo(reel,ex); |
x = gadd(x,reel); if (kk && ex < G) break; |
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; |
q1 = addsi(k+1,a); ep->value = (void*) q1; |
p1 = lisexpr(ch); if (did_break()) err(breaker,"sumpos"); |
p1 = lisexpr(ch); if (did_break()) err(breaker,"sumpos"); |
gaffect(p1,reel); x = gadd(reel, gmul2n(x,1)); |
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)); |
s = gadd(s,gmul(x,p1)); |
az = divii(mulii(mulss(N-k,N+k),shifti(az,1)),mulss(k+1,k+k+1)); |
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 |
GEN |
sumpos2(entree *ep, GEN a, char *ch, long prec) |
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; |
GEN p1,r,q1,reel,s,pol,dn,x, *stock; |
|
|
if (typ(a) != t_INT) err(talker,"non integral index in sumpos2"); |
if (typ(a) != t_INT) err(talker,"non integral index in sumpos2"); |
|
|
push_val(ep, a); |
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; |
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++) |
for (k=1; k<=N; k++) |
{ |
{ |
if (odd(k) || !stock[k]) |
if (odd(k) || !stock[k]) |
{ |
{ |
x=gzero; r=stoi(2*k); |
x = gzero; r = stoi(2*k); |
for(kk=0;;kk++) |
for(kk=0;;kk++) |
{ |
{ |
long ex; |
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"); |
p1 = lisexpr(ch); if (did_break()) err(breaker,"sumpos2"); |
gaffect(p1,reel); ex = expo(reel)+kk; setexpo(reel,ex); |
gaffect(p1,reel); ex = expo(reel)+kk; setexpo(reel,ex); |
x=gadd(x,reel); if (kk && ex < G) break; |
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; |
q1 = addsi(k,a); ep->value = (void*) q1; |
p1 = lisexpr(ch); if (did_break()) err(breaker,"sumpos2"); |
p1 = lisexpr(ch); if (did_break()) err(breaker,"sumpos2"); |
gaffect(p1,reel); stock[k] = gadd(reel, gmul2n(x,1)); |
gaffect(p1,reel); stock[k] = gadd(reel, gmul2n(x,1)); |
} |
} |
} |
} |
pop_val(ep); s = gzero; |
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++) |
for (k=1; k<=lgef(pol)-2; k++) |
{ |
{ |
p1 = gmul((GEN)pol[k+1],stock[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); |
s = gadd(s,p1); |
} |
} |
tetpil=avma; return gerepile(av,tetpil,gdiv(s,dn)); |
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) **/ |
/** 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 */ |
/* we need to make a copy in any case, cf forpari */ |
static GEN |
static GEN |
Line 909 fix(GEN a, long prec) |
|
Line 942 fix(GEN a, long prec) |
|
gaffect(a,p); return p; |
gaffect(a,p); return p; |
} |
} |
|
|
extern GEN polint_i(GEN xa, GEN ya, GEN x, long n, GEN *ptdy); |
#if 0 |
|
|
GEN |
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); |
a = fix(a,prec); |
b = fix(b,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); |
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); |
a = fix(a, prec); |
b = fix(b, 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); |
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) |
for (it=1, j=1; j<JMAX; j++, it*=3) |
{ |
{ |
h[j] = ldivrs((GEN)h[j-1],9); |
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); |
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 |
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 |
GEN |
polzag(long n, long m) |
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"); |
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++) |
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++) |
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 */ |
#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 |
GEN |
polzagreel(long n, long m, long prec) |
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"); |
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++) |
for (k=1; k<=d1; k++) |
{ |
{ |
v[d1-k+1]=un; |
v[d1-k+1] = un; |
for (j=1; j<k; j++) |
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++) |
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++) |
for (j=0; j<=r; j++) |
{ |
{ |
if (j) g=derivpol(g); |
if (j) g = derivpol(g); |
if (j || !(m&1)) |
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++) |
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 |
#ifdef _MSC_VER |
#pragma optimize("g",on) |
#pragma optimize("g",on) |
Line 1183 polzagreel(long n, long m, long prec) |
|
Line 1242 polzagreel(long n, long m, long prec) |
|
GEN |
GEN |
zbrent(entree *ep, GEN a, GEN b, char *ch, long prec) |
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; |
GEN c,d,e,tol,tol1,min1,min2,fa,fb,fc,p,q,r,s,xm; |
|
|
a = fix(a,prec); |
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 */ |
else { d = xm; e = d; } /* bound decreasing too slowly, use bisection */ |
a = b; fa = fb; |
a = b; fa = fb; |
if (gcmp(gabs(d,0),tol1) > 0) b = gadd(b,d); |
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); |
ep->value = (void*)b; fb = lisexpr(ch); |
} |
} |
if (iter > itmax) err(talker,"too many iterations in solve"); |
if (iter > itmax) err(talker,"too many iterations in solve"); |