intfuncinit

Function: intfuncinit
Section: sums
C-Name: intfuncinit0
Prototype: V=GGED0,L,p
Help: intfuncinit(t=a,b,f,{m=0}): initialize tables for integrations
 from a to b using a weight f(t). For integral transforms such
 as Fourier or Mellin transforms.
Wrapper: (,,G)
Description:
  (gen,gen,gen,?small):gen:prec intfuncinit(${3 cookie}, ${3 wrapper}, $1, $2, $4, $prec)
Doc: initialize tables for use with integral transforms (such as Fourier,
 Laplace or Mellin transforms) in order to compute
 $$ \int_a^b f(t) k(t,z) \, dt $$
 for some kernel $k(t,z)$.
 The endpoints $a$ and $b$ are coded as in \kbd{intnum}, $f$ is the
 function to which the integral transform is to be applied and the
 nonnegative integer $m$ is as in \kbd{intnum}: multiply the number of
 sampling points roughly by $2^m$, hopefully increasing the accuracy. This
 function is particularly useful when the function $f$ is hard to compute,
 such as a gamma product.

 \misctitle{Limitation} The endpoints $a$ and $b$ must be at infinity,
 with the same asymptotic behavior. Oscillating types are not supported.
 This is easily overcome by integrating vectors of functions, see example
 below.

 \misctitle{Examples}

 \item numerical Fourier transform
 $$F(z) = \int_{-\infty}^{+\infty} f(t)e^{-2i\pi z t}\, dt. $$
 First the easy case, assume that $f$ decrease exponentially:
 \bprog
    f(t) = exp(-t^2);
    A = [-oo,1];
    B = [+oo,1];
    \p200
    T = intfuncinit(t = A,B , f(t));
    F(z) =
    { my(a = -2*I*Pi*z);
      intnum(t = A,B, exp(a*t), T);
    }
    ? F(1) - sqrt(Pi)*exp(-Pi^2)
    %1 = -1.3... E-212
 @eprog\noindent
 Now the harder case, $f$ decrease slowly: we must specify the oscillating
 behavior. Thus, we cannot precompute usefully since everything depends on
 the point we evaluate at:
 \bprog
    f(t) = 1 / (1+ abs(t));
    \p200
    \\ Fourier cosine transform
    FC(z) =
    { my(a = 2*Pi*z);
      intnum(t = [-oo, a*I], [+oo, a*I], cos(a*t)*f(t));
    }
    FC(1)
 @eprog
 \item Fourier coefficients: we must integrate over a period, but
 \kbd{intfuncinit} does not support finite endpoints.
 The solution is to integrate a vector of functions !
 \bprog
 FourierSin(f, T, k) =  \\ first k sine Fourier coeffs
 {
   my (w = 2*Pi/T);
   my (v = vector(k+1));
   intnum(t = -T/2, T/2,
      my (z = exp(I*w*t));
      v[1] = z;
      for (j = 2, k, v[j] = v[j-1]*z);
      f(t) * imag(v)) * 2/T;
 }
 FourierSin(t->sin(2*t), 2*Pi, 10)
 @eprog\noindent The same technique can be used instead of \kbd{intfuncinit}
 to integrate $f(t) k(t,z)$ whenever the list of $z$-values is known
 beforehand.

 Note that the above code includes an unrelated optimization: the
 $\sin(j w t)$ are computed as imaginary parts of $\exp(i j w t)$ and the
 latter by successive multiplications.

 \item numerical Mellin inversion
 $$F(z) = (2i\pi)^{-1} \int_{c -i\infty}^{c+i\infty} f(s)z^{-s}\, ds
  = (2\pi)^{-1} \int_{-\infty}^{+\infty}
     f(c + i t)e^{-\log z(c + it)}\, dt. $$
 We take $c = 2$ in the program below:
 \bprog
    f(s) = gamma(s)^3;  \\ f(c+it) decrease as exp(-3Pi|t|/2)
    c = 2; \\ arbitrary
    A = [-oo,3*Pi/2];
    B = [+oo,3*Pi/2];
    T = intfuncinit(t=A,B, f(c + I*t));
    F(z) =
    { my (a = -log(z));
      intnum(t=A,B, exp(a*I*t), T)*exp(a*c) / (2*Pi);
    }
 @eprog

 \synt{intfuncinit}{void *E, GEN (*eval)(void*,GEN), GEN a,GEN b,long m, long prec}.