xref: /plan9/sys/src/cmd/map/libmap/cubrt.c (revision 59cc4ca53493a3c6d2349fe2b7f7c40f7dce7294)
1*59cc4ca5SDavid du Colombier #include <u.h>
2*59cc4ca5SDavid du Colombier #include <libc.h>
33e12c5d1SDavid du Colombier #include "map.h"
43e12c5d1SDavid du Colombier 
53e12c5d1SDavid du Colombier double
cubrt(double a)63e12c5d1SDavid du Colombier cubrt(double a)
73e12c5d1SDavid du Colombier {
83e12c5d1SDavid du Colombier 	double x,y,x1;
93e12c5d1SDavid du Colombier 	if(a==0)
103e12c5d1SDavid du Colombier 		return(0.);
113e12c5d1SDavid du Colombier 	y = 1;
123e12c5d1SDavid du Colombier 	if(a<0) {
133e12c5d1SDavid du Colombier 		y = -y;
143e12c5d1SDavid du Colombier 		a = -a;
153e12c5d1SDavid du Colombier 	}
163e12c5d1SDavid du Colombier 	while(a<1) {
173e12c5d1SDavid du Colombier 		a *= 8;
183e12c5d1SDavid du Colombier 		y /= 2;
193e12c5d1SDavid du Colombier 	}
203e12c5d1SDavid du Colombier 	while(a>1) {
213e12c5d1SDavid du Colombier 		a /= 8;
223e12c5d1SDavid du Colombier 		y *= 2;
233e12c5d1SDavid du Colombier 	}
243e12c5d1SDavid du Colombier 	x = 1;
253e12c5d1SDavid du Colombier 	do {
263e12c5d1SDavid du Colombier 		x1 = x;
273e12c5d1SDavid du Colombier 		x = (2*x1+a/(x1*x1))/3;
283e12c5d1SDavid du Colombier 	} while(fabs(x-x1)>10.e-15);
293e12c5d1SDavid du Colombier 	return(x*y);
303e12c5d1SDavid du Colombier }
31