Re: Integrating functions

[Jean-Marc Sac-Epée <jmse@poncelet.univ-metz.fr>]

Karim BELABAS a écrit :

> [Jean-Marc Sac-Epée:]
> > Thanks to your help, I know now how to integrate a function  using
> > intnum0:
> >
> > > #include <math.h>
> > > #include <stdio.h>
> > > #include <pari.h>
> > >
> > > GEN integrate(char *f, char *var, GEN a, GEN b, long prec)
> > > {
> > >   entree *ep_var = is_entry(var);
> > >   return intnum0(ep_var, a,b, f, 0,prec);
> > > }
> > >
> > > int
> > > main()
> > > {
> > >   pari_init(1000000,2);
> > >   output(integrate("x^2+3", "x", gzero, gun, DEFAULTPREC));
> > >   return 0;
> > > }
> >
> > Now, I would like to know if it is also possible with more complicated
> > fonctions. For example, it is easy to prove that for each t>0, the
> > polynomial 10*x^5 + 0.1*x^3 - t*x^2 - t has a unique real root. So, we
> > can define the fonction f: t -> this real root.
> >
> > Under gp, it is easy to integrate this function;
> > ? f(t)=real(polroots(10*x^5+0.1*x^3-t*x^2-t)[1])
> > ? intnum(t=0,1,f(t))
> >
> > In a C program, I can define this fonction for example with
> >
> > GEN f(GEN t)
> > {
> >  const long var = 0;
> >  const long prec = 4;
> >  GEN pol;
> >  pol = cgetg(deg+3,t_POL);
> >  pol[1] = evalsigne(1) | evalvarn(var) | evallgef(deg+3);
> >  pol[deg+2]=dbltor(10); pol[deg]=dbltor(0.1);
> >  pol[deg+1]=pol[deg-2]=gzero;
> >  pol[deg-1]=pol[deg-3]=gneg(t);
> >  return(greal(roots(pol,prec)[1]));
> > }
>
> I'd rather use
>
> GEN f(GEN t)
> {
>   const long prec = DEFAULTPREC; /* 64bit portable */
>   const long var = 0;
>   const long deg = 5;
>   GEN *pol, x = cgetg(deg+3,t_POL);
>
>   pol[1] = evalsigne(1) | evalvarn(var) | evallgef(deg+3);
>   pol = (GEN*)(x+2); /* strip codewords, avoid cast problems */
>   pol[0] = gneg(t);
>   pol[1] = gzero;
>   pol[2] = pol[0];
>   pol[3] = dbltor(0.1);
>   pol[4] = gzero;
>   pol[5] = dbltor(10.);
>   return greal(roots(x,prec)[1]);
> }
>
> [possibly not using the (GEN*) cast, but definitely stripping the codewords
> so as to be able to identify the coefficient associated to a given degree
> without mentally translating...]
>
> > But now, what is the correct syntax to integrate f?
>
> Two different ways, both of them kludges.
>
> 1) [not generic but works for your example...]
>   f = "real(polroots(10*x^5+0.1*x^3-t*x^2-t)[1])";
>   ep_var = is_entry("t");
>
> 2) [not recommended: you will need to link gp.o to your binary for that one...]
>   gpinstall(lib_f, "G" , "f", mylib.so);
>   f = "f(t)";
>   ep_var = is_entry("t");
>
> \begin{digression}
>
> In my opinion, it would be a _much_ better idea to copy and slightly modify
> the underlying intnum code (src/language/sumiter.c:qromb() in the above
> example) so as to replace the 'entree *ep' and 'char *ch' parameter with a C
> pointer to function GEN (*f)(GEN). Then all groups of the form
>
>    ep->value = (void*)x;   /* or push_val(ep, x) */
>    y = lisexpr(ch);
>
> can be replaced by
>
>     y = f(x);
>
> all other references to ep can be junked.
>
> This is straightforward (40 lines, about 2 minutes work), yields a routine
> which is easier to use, and should be more efficient [no parsing of GP
> expressions]. (in fact, I've just done it, see P.S below).
>
> \begin{technical}
> A generic, hence better, approach would be to use two parameters 'f' and
> 'data' (to avoid using global variables to properly define an 'f' needing
> external data)
>
>     y = f(x, data);
>
> 'f' being assumed to know how to use 'data' (possibly ignoring it).
>
> In fact all PARI functions using 'entree' arguments should have been written
> in this way, so as to be easy to integrate to library code. The GP version
> would have then been trivial to code using as 'f'
>
> GEN
> eval_GP_fun(GEN x, void *data)
> {
>   Ep(data)->value = (void*)x;
>   return lisexpr(Ch(data));
> }
> [with suitably defined macros Ep and Ch accessing the relevant fields of
> data, after proper casts, etc.]
>
> Of course, if needed, data could include pre/postprocessing, etc. Such a
> scheme is actually already used by a number of internal library functions,
> e.g fincke_pohst [ finds small vectors satisfying "certain conditions" in
> lattices using a generic Fincke-Pohst algorithm. The polredabs and qfminim
> functions use it for instance ]
>
> Since the stable version 2.1 is forthcoming, I'll probably rewrite everything
> in this style for version 2.2.0 (the next alpha...)  [should be an hour work
> now: there are MANY such functions]
> \end{technical}
> \end{digression}
>
> Hope this helps,
>
>     Karim.
>
> ===========================================================================
>
> P.S: In fact, I just wasted the necessary 2 minutes (was about 40 seconds in
> fact) adapting qromb [it's not generic, but you'll adapt it...].
> Now
>
>   integrate(f, stoi(1), stoi(2), DEFAULTPREC)
>
> works without any kludge nor intervention from the parser code. It's a
> shameless copy/paste job, I've not re-indented, or cleaned up the code in
> any way:
>
> #define JMAX 25
> #define JMAXP JMAX+3
> #define KLOC 4
>
> /* we need to make a copy in any case, cf forpari */
> static GEN
> fix(GEN a, long prec)
> {
>   GEN p = cgetr(prec);
>   gaffect(a,p); return p;
> }
>
> GEN
> integrate(GEN (*f)(GEN), GEN a, GEN b, long prec)
> {
>   GEN ss,dss,s,h,p1,p2,qlint,del,x,sum;
>   long av = avma, av1, tetpil,j,j1,j2,lim,it,sig;
>
>   a = fix(a,prec);
>   b = fix(b,prec);
>   qlint=subrr(b,a); sig=signe(qlint);
>   if (!sig)  { avma=av; return gzero; }
>   if (sig<0) { setsigne(qlint,1); s=a; a=b; b=s; }
>
>   s=new_chunk(JMAXP);
>   h=new_chunk(JMAXP);
>   h[0] = (long)realun(prec);
>
>   p1=f(a); if (p1 == a) p1=rcopy(p1);
>   p2=f(b);
>   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);
>     av1=avma; del=divrs(qlint,it); x=addrr(a,shiftr(del,-1));
>     for (sum=gzero,j1=1; j1<=it; j1++,x=addrr(x,del))
>     {
>       sum=gadd(sum,f(x));
>     }
>     sum = gmul(sum,del); p1 = gadd((GEN)s[j-1],sum);
>     tetpil = avma;
>     s[j]=lpile(av1,tetpil,gmul2n(p1,-1));
>
>     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)-j-6;
>       if (j1-j2 > lim || j1 < -lim)
>       {
>         if (gcmp0(gimag(ss))) ss=greal(ss);
>         tetpil=avma; return gerepile(av,tetpil,gmulsg(sig,ss));
>       }
>       avma=tetpil;
>     }
>   }
>   err(intger2); return NULL; /* not reached */
> }
> --
> Karim Belabas                    email: Karim.Belabas@math.u-psud.fr
> Dep. de Mathematiques, Bat. 425
> Universite Paris-Sud             Tel: (00 33) 1 69 15 57 48
> F-91405 Orsay (France)           Fax: (00 33) 1 69 15 60 19
> --
> PARI/GP Home Page: http://www.parigp-home.de/

Thank you very much, Karim, for this very detailed answer.

J.-Marc

begin:vcard 
n:Sac-Epée;Jean-Marc 
tel;fax:03 87 31 52 73
tel;work:03 87 54 72 69
x-mozilla-html:FALSE
url:http://www.mmas.univ-metz.fr/~jmse
org:Université de Metz;Mathématiques
adr:;;;;;;
version:2.1
email;internet:jmse@poncelet.univ-metz.fr
title:Ingénieur de Recherches
x-mozilla-cpt:;0
fn:Jean-Marc Sac-Epée
end:vcard