#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;
}