#include "doubledouble.h"
#include "ddnrutil.h"
#include "ddbsstep.h"

#include <math.h>

/* These are mutually exclusive */
#define EXP_1 0
#define PI 1

void f(int n, doubledouble t, doubledouble *y, doubledouble *ydot)
{
    int i;
    for(i=0;i<n;i++)
#if EXP_1
	ydot[i] = y[i];
#elif PI	// y' = sqrt(1-t^2) from -1 to 1 = PI
	ydot[i] = sqrt(fabs(1.0-t*t));
#endif
}


doubledouble Doit(doubledouble eps)
{
#if EXP_1
    doubledouble y = 1., x1=0, x2=1;
#elif PI
    doubledouble y = 0., x1=-1, x2=1;
#endif
    doubledouble h0=1, hmin=1e-36;
    int n=1, nok,nbad, status, maxSteps = 10000;

    //kmax = 0;
    //dxsav=0.1;
    cout << setprecision(4);
    cout << "y(" << x1 << ") = " << y << '\n';
    while((status = ddodeint(&y, n, &x1, x2, eps, h0, hmin, maxSteps, &nok, &nbad, f, ddbsstep)) != 0)
    {
	cout << "status(0=OK, 1=stepsizeUnderflow, 2=tooManySteps) = " << status << '\n';
	cout << "y(" << x1 << ") = " << y << '\n';
	cout << "nok = " << nok << ", nbad = " << nbad << '\n';
	cout << "keep going...\n" ;
	if(status == 1)
	    h0 *= 100;
    }
    cout << "status(0=OK, 1=stepsizeUnderflow, 2=tooManySteps) = " << status << '\n';
    cout << "y(" << x1 << ") = " << setprecision(30) << y << setprecision(0) << '\n';
    cout << "nok = " << nok << ", nbad = " << nbad << '\n';
    return y;
}

int main(void)
{
    doubledouble eps = 1e-1;
#if EXP_1
    doubledouble answer = exp(doubledouble(1.0,0));
#elif PI
    doubledouble answer = doubledouble::Pi/2;
#endif
    for(eps=0.1; eps >= doubledouble(1e-30); eps*=0.1)
    {
	doubledouble val;
	cout << " EPS = " << eps <<'\n';
	val = Doit(eps);
	cout << "\nshould=" << setprecision(30) << answer << setprecision(0) << '\n';
	cout << setprecision(3);
	cout << "for eps = " << eps << " err is " << fabs(answer-val) << '\n';
    }
    return 0;
}
