#include<stdio.h>
#include<assert.h>
#include<gsl/gsl_multiroots.h>
#include<gsl/gsl_vector.h>
#include<gsl/gsl_errno.h>

int tozero (const gsl_vector * v, void * params, gsl_vector * f){
	double y = *(double*)params;
	double x = gsl_vector_get(v,0);
	gsl_vector_set(f,0,x*x-y);
	return GSL_SUCCESS;
}

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

	gsl_multiroot_function F;
	F.f = tozero;
	F.n = 1;
	F.params = (void*)&y;
	
	gsl_multiroot_fsolver* s = gsl_multiroot_fsolver_alloc (gsl_multiroot_fsolver_broyden, 1);

	gsl_vector* v = gsl_vector_alloc(1);
	gsl_vector_set(v, 0, y/2);
	gsl_multiroot_fsolver_set (s, &F, v);

	int iter=0, status;
	do {
      iter++;
      status = gsl_multiroot_fsolver_iterate (s);
      if (status)   /* check if solver is stuck */
        break;
      status = gsl_multiroot_test_residual (s->f, 1e-7);
    }
  while (status == GSL_CONTINUE && iter < 1000);
	double result = gsl_vector_get(s->x,0);
	gsl_multiroot_fsolver_free (s);
	return result;
}
	
int main(){
	for(double x=0; x<11; x+=0.1) printf("%g %g\n",x,pow12(x));
return 0;
}
