FENIX_libc/math/sqrt.c

27 lines
532 B
C

#include <math.h>
typedef union {
double d;
struct {
unsigned long int mantissa : 52;
unsigned int exponent : 11;
unsigned int sign_bit : 1;
} parts;
} _dbl_cast;
double sqrt(double x) {
_dbl_cast S = { .d = x };
double estimate = (0.5 + 0.5 * S.parts.mantissa) *
exp2(S.parts.exponent - 1023);
double a_n, b_n;
for(int i = 0; i < 15; i++) {
a_n = (x - (estimate * estimate)) / (2 * estimate);
b_n = estimate + a_n;
estimate = b_n - ((a_n * a_n) / (2 * b_n));
}
return estimate;
}