45 lines
839 B
C
45 lines
839 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 _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) *
|
|
_iexp2(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;
|
|
}
|