#!/usr/bin/perl -w
use strict;
use PDL;
use PDL::NiceSlice;
use Inline qw(Pdlpp);

my $dx = 0.1;
my $x = $dx * pdl [0..100];
my $y = $x**2 + 2*$x**3 + exp($x/2);
my $g_num1 = $dx * PDL::cumuintover($y);
my $g_ana = $x**3/3 + 2*$x**4/4 + 2 * exp($x/2) - 2;
print "     x   |    num  |    ana  |   rel err\n";
wcols "%8.3g ", $x(0:-1:10), $g_num1(0:-1:10), $g_ana(0:-1:10), 
  ($g_num1(0:-1:10) - $g_ana(0:-1:10))/$g_ana(0:-1:10);


no PDL::NiceSlice;

__DATA__

__Pdlpp__


# from equation 4.1.14, Press et al. 2nd Ed chapter 4.1
pp_def('cumuintover',
       Pars => 'f(n); int+ [o]g(n);',
       Code => '
int i;
int ns = $SIZE(n);
$GENERIC(g) tmp = 0;
threadloop %{
tmp = (3.0/8.0) * $f(n=>0) +
      (7.0/6.0) * $f(n=>1) +
      (23.0/24.0) * $f(n=>2);
for (i=5;i<ns;i++) { $g(n=>i) = tmp +
                             (23./24.) * $f(n=>i-2) + 
                             (7./6.)   * $f(n=>i-1) + 
                             (3./8.)   * $f(n=>i); 
                     tmp += $f(n=>i-2);
                    }
$g(n=>0) = 0;
$g(n=>1) = ($f(n=>0) + $f(n=>1))/2.0;
$g(n=>2) = ($f(n=>0) + 4. * $f(n=>1) + $f(n=>2) ) / 3.0;
$g(n=>3) = (3/8) * ($f(n=>0) + $f(n=>3) + 3. * ( $f(n=>1) + $f(n=>2) ) );
$g(n=>4) = (1./45.) * (  14. * ( $f(n=>0) + $f(n=>4) )
		    + 64. * ( $f(n=>1) + $f(n=>3) )
		    + 24. * $f(n=>2) 
		    );
%}
'
);