dsqrt.c

#include "mathutils.h"
#include <assert.h>

#define TOL 0.0001

/*****************************************************************
	Utilizes Newton's Method for computing square roots with a 
	*rough* estimate obtained first by first finding a suitable n
	for sqrt(x) = sqrt(a) * 2^n such that .5 <= a < 4.
	x must be a non-zero real; function terminates on a certain 
	degree of convergence
  ***************************************************************/
double newton_dsqrt(double x){
	assert(x >= 0.0);
	long int_x = (long)x;
	int n = 1;
	while(!(1 <= 1 + int_x && 1 + int_x < 5)){
		int_x = int_x >> 2;
		n++;
	}
	double seed = (double) (1 << n);
	double y_nx,r,y_pv;
	y_nx = r = 0.0;
	y_pv = seed;
	do{
		y_nx = 0.5*(y_pv + (x/y_pv));
		r = (ABS(y_nx - y_pv)/y_pv);
		y_pv = y_nx;
	}while(r > TOL);
	return y_nx;
}