Re: Integrating functions

[Karim BELABAS <Karim.Belabas@math.u-psud.fr>]

[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/