diff --git a/math/sqrt.c b/math/sqrt.c index 9d2b077..83277e5 100644 --- a/math/sqrt.c +++ b/math/sqrt.c @@ -9,10 +9,28 @@ typedef union { } parts; } _dbl_cast; +double _iexp2(int x) { + double result = 1.0; + if(x > 0) { + for(int i = 0; i < x; i++) { + result *= 2; + } + } + else if(x < 0) { + for(int i = 0; i > x; i--) { + result /= 2; + } + } + return result; +} + +/* + First gets an initial estimate then refines it using the Bakhshali method +*/ double sqrt(double x) { _dbl_cast S = { .d = x }; double estimate = (0.5 + 0.5 * S.parts.mantissa) * - exp2(S.parts.exponent - 1023); + _iexp2(S.parts.exponent - 1023); double a_n, b_n;