#include "libscl.h"
using namespace std;
using namespace scl;

REAL f(REAL x, INTEGER k); 

const INTEGER selection = 5;
const REAL tol = 1.0e-10;

int main()
{
  #include "golub.cpp"

  INTEGER n=N[selection];

  INTEGER i,fail;

  realmat w,x;

  fail=hquad(n,x,w);
  cout << " fail = " << fail << " n = " << n << '\n';
  // cout << x << w << '\n';

  REAL r[N[selection]];
  INTEGER idx[N[selection]];
  for (i=0;i<N[selection];i++) r[i]=x[i+1];
  heapsort(r,n,idx);
  //for (i=0;i<N[selection];i++) cout << i << ' ' << idx[i] << '\n';

  realmat test(N[selection],6);

  for (i=1;i<=N[selection];i++) {
    test(i,1)=x[idx[i-1]+1];
    test(i,2)=X[selection][i-1];
    test(i,3)=x[idx[i-1]+1]-X[selection][i-1]; 
    if (abs(test(i,3))<fabs(test(i,2))*tol) test(i,3)=0.;
    test(i,4)=w[idx[i-1]+1];
    test(i,5)=W[selection][i-1]; 
    test.check(i,6)=w.check(idx[i-1]+1)-W[selection][i-1]; 
    if (abs(test(i,6))<fabs(test(i,5))*tol) test(i,6)=0.;
  }

  cout << test << '\n';

  REAL sum; 
  INTEGER k;
  fail=hquad(n,x,w); 
  for (k=1;k<=2*n-1;k++) {
    sum=0.0;
    for(i=1;i<=n;i++) sum+=f(x[i],k)*w[i]; 
    cout << "f(x)=pow(z*sqrt(2/pi)," << k << "), inegral = " << sum  << '\n';
  }
  return 0;
}

/* 
Moments of a normal, to test gauss-hermite
*/

REAL f(REAL x, INTEGER k) {
  REAL pi=3.141592653589793;
  REAL z,fz,fx;
  z=sqrt(2.0)*x;
  fz=pow(z,k);
  fx=fz/sqrt(pi);
  return fx;
}











