math/log.c

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

math/sqrt.c

@@ -9,28 +9,10 @@ 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) *
-    _iexp2(S.parts.exponent - 1023);
+    exp2(S.parts.exponent - 1023);
 
   double a_n, b_n;