/* This code is hereby released to the public domain.
 * ~aaaaaa123456789, 2020-01-01                       */

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
#include <time.h>

#define FRACTIONAL_BITS 20

int32_t integer_inverse (unsigned value) {
  if (!value) return INT32_MAX;
  return ((2 << FRACTIONAL_BITS) / value + 1) >> 1;
}

int32_t multiply (int32_t first, int32_t second) {
  int64_t result = ((int64_t) first * (int64_t) second + ((int64_t) 1 << (FRACTIONAL_BITS - 1))) >> FRACTIONAL_BITS;
  if (result > INT32_MAX) return INT32_MAX;
  if (result < INT32_MIN) return INT32_MIN;
  return result;
}

int32_t approximate_square_root (int32_t operand, unsigned iterations) {
  unsigned iteration;
  int32_t result = (operand + (1 << FRACTIONAL_BITS)) >> 1;
  int32_t factor = ((1 << FRACTIONAL_BITS) - operand) >> 1;
  int32_t power = factor;
  int32_t coefficient = 1 << FRACTIONAL_BITS;
  for (iteration = 2; iteration <= iterations; iteration ++) {
    coefficient = multiply(coefficient, (2 * iteration - 3) << FRACTIONAL_BITS);
    coefficient = multiply(coefficient, integer_inverse(iteration));
    power = multiply(power, factor);
    result -= multiply(power, coefficient);
  }
  return result;
}

int32_t true_square_root (int32_t operand) {
  #if FRACTIONAL_BITS & 1
    #error This function requires an even number of fractional bits!
  #endif
  if (operand < 0) return 0;
  return ldexp(sqrt(operand), FRACTIONAL_BITS / 2) + 0.5;
}

unsigned calculate_error (int32_t approximate, int32_t real) {
  double error = (double) approximate / (double) real - 1.0;
  error = fabs(error);
  return error * 10000 + 0.5;
}

int main (void) {
  unsigned tests = 100, limit = 6, current;
  srand(clock() ^ time(NULL));
  fputs("  value    sqrt        1 term    ", stdout);
  for (current = 2; current <= limit; current ++) printf("     %2u terms   ", current);
  fputs("\n-------- -------- ---------------", stdout);
  for (current = 2; current <= limit; current ++) fputs(" ---------------", stdout);
  putchar('\n');
  while (tests --) {
    int32_t value = ((rand() ^ (rand() << 12)) >> 4) & ((1 << (FRACTIONAL_BITS + 1)) - 1); // random between 0.0 and 2.0
    int32_t true_root = true_square_root(value);
    printf("%.6f %.6f", ldexp(value, -FRACTIONAL_BITS), ldexp(true_root, -FRACTIONAL_BITS));
    for (current = 1; current <= limit; current ++) {
      int32_t root = approximate_square_root(value, current);
      printf(" %.6f", ldexp(root, -FRACTIONAL_BITS));
      unsigned error = calculate_error(root, true_root);
      if (error >= 10000)
        printf(" %5u%%", (error + 50) / 100);
      else
        printf(" %2u.%02u%%", error / 100, error % 100);
    }
    putchar('\n');
  }
  return 0;
}