include "gsl.idp"
// FFCS - avoid using P2 and P3 which already define something very precise
real [int] p2=[ 2,-3,1];
real [int] p3=[ 8,-10,1,1];
real [int] X(3);
complex[int] Z(4);
real [int] Q4=[ -1,0,0,0,1];
gslpolysolvequadratic(p2,X);
cout << X[0] << " " << X[1] << endl; 
gslpolysolvecubic(p3,X);
cout << X[0] << " " << X[1] << " " << X[2] << endl; 
gslpolycomplexsolve(Q4,Z); 
cout << Z[0] << " " << Z[1] << " " << Z[2] << " " << Z[3] <<endl;

cout << "\n\n gsl random :\n\n";
 cout << " list of generator " << ngslrng << endl;
 for(int i=0; i <ngslrng ; ++i)
   cout << i << " name  " <<  gslname(gslrngtype(i)) << endl; 
gslrng r=gslrngtype(1);// chose a generator 
// cout << r << endl;  no print .. 
for (int i = 0; i < 10; i++) 
   cout << gslrnguniform(r) << endl;;  
 
 real x0=0.2;
 int n=10;
 real[int,int] dspline(2,n+6);
 for(int i=0; i<= n+5; ++i)
  {
    real xx = square(real(i)/n)+x0; 
    real yy = sin(xx*pi*2); 
    dspline(0,i)= xx;
    dspline(1,i)= yy;    
  }
  /*
 the 7  type of spline of gsl are 
  0: gslinterpcspline (default .
  1: gslinterpakima
  2: gslinterpsteffen
  3: gslinterplinear
  4:  gslinterppolynomial
  5: gslinterpcsplineperiodic
  6:  gslinterpakimaperiodic 
  */
  gslspline spline(gslinterplinear,dspline);
  int mm=100;
 real[int] xx(mm+1),yy(mm+1),zz(mm+1),dyy(mm+1),ddyy(mm+1),dzz(mm+1),ddzz(mm+1);
 real err=0; 
 for (int i=0; i<=mm; ++i)
 { xx[i]= real(i)/mm+x0;
   yy[i]= spline(xx[i]);
   dyy[i]= spline.d(xx[i]);
   ddyy[i]= spline.dd(xx[i]);
   zz[i]= sin(pi*xx[i]*2);
   dzz[i]= 2*pi*cos(pi*xx[i]*2);
   ddzz[i]= -2*pi*2*pi*sin(pi*xx[i]*2);
   err = max( err , abs(yy[i]-zz[i]));
   cout << xx[i] << " " << yy[i] << "  ~ " << zz[i] << ",  " << dyy[i] << " ~ " << dzz [i] <<  " , " << ddyy[i] << " " << ddzz[i] <<  endl; 
 }
 //  test out of  bound v 3.50
 // OUT => VALUE AT BOUND ... FH OCT 2015 
 cout << dspline(1,0) << " ==  " << spline(-1) << " " << spline(100) << " == "<< dspline(1,dspline.m-1) << endl; 
 cout << " err = " << err << endl; 
 plot([xx,yy],[xx,zz],cmm= "err = " + err);
 assert(err < 0.3);
