Compare commits

...

2 commits

2 changed files with 20 additions and 2 deletions

View file

@ -15,7 +15,7 @@ double _agm(double x, double y) {
} }
double log(double x) { double log(double x) {
double s = x * exp2(24); double s = x * 16777216;
double nat_log = (M_PI / (2 * _agm(1, 4/s))) - (24 * M_LN2); double nat_log = (M_PI / (2 * _agm(1, 4/s))) - (24 * M_LN2);
return nat_log; return nat_log;
} }

View file

@ -9,10 +9,28 @@ typedef union {
} parts; } parts;
} _dbl_cast; } _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) { double sqrt(double x) {
_dbl_cast S = { .d = x }; _dbl_cast S = { .d = x };
double estimate = (0.5 + 0.5 * S.parts.mantissa) * double estimate = (0.5 + 0.5 * S.parts.mantissa) *
exp2(S.parts.exponent - 1023); _iexp2(S.parts.exponent - 1023);
double a_n, b_n; double a_n, b_n;