#include<stdio.h>
#include<math.h>
#include<assert.h>
#define M_PI 3.1415926535897932384626433
#include<gsl/gsl_odeiv2.h>
#include<gsl/gsl_errno.h>

int n=0;
int rhs (double t, const double y[], double dydt[], void * params) {
	n++;
	dydt[0]=1.0/2/y[0];
	return GSL_SUCCESS;
}

double mysqrt(double x) {
	assert(x>=0);
	if(x==0) return 0;
	if(x<1) return 1/mysqrt(1/x);
	if(x>9) return 3*mysqrt(x/9);

	gsl_odeiv2_system s = { .function = rhs, .jacobian = NULL, .dimension = 1, .params = NULL };

	gsl_odeiv2_driver* d =
	gsl_odeiv2_driver_alloc_y_new( &s, gsl_odeiv2_step_rk8pd, 1e-3, 1e-15, 1e-15);

	double t=1, y[]={1};
	int status = gsl_odeiv2_driver_apply( d, &t, x, y );
	if(status != GSL_SUCCESS) fprintf(stderr,"my_sin: status = %i\n",status);
	return y[0];
}

int main(void) {
	n=0;
	fprintf(stderr,"%20.15g\n", mysqrt(9));
	fprintf(stderr,"n=%i\n", n);
	for(double x=0; x<=10; x+=0.1) printf("%g %g\n", x, mysqrt(x));
	printf("\n\n");
	printf("%g %g\n", 1./4, 1./2);
	printf("%i %i\n", 1, 1);
	printf("%i %i\n", 4, 2);
	printf("%i %i\n", 9, 3);
return 0;
}