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

typedef uint32_t u32;

void heron_sqrt(u32 *mem) {
  int i = mem[0];
  while (i > 0) {
    u32 s = mem[i];
    if (s > 1) {
      // x_0 = (s + 1) / 2
      // without increment to avoid overflow for 0xffffffff
      u32 x = s >> 1;
      u32 old_x = x;

      while(1) {
        // x_{n + 1} = (x_n + (s / x_n)) / 2
        x = (x >> 1) + ((s/x) >> 1);
        if (old_x <= x) {
          break;
        }
        old_x = x;
      }
      mem[i] = x;
    }
    i--;
  }
}

int main() {
  u32 data[] = {0,
    1,
    2,
    4,
    3<<2,
    1<<3,
    1<<4,
    1<<5,
    1<<6,
    5<<2,
    1<<7,
    1<<8,
    1<<16,
    1<<17,
    1<<18,
    1<<19,
    1<<20, // 32 - 14 = 20
    1<<24, // 32 -  8 = 24
    1<<30,
    0xffffffff};
  const int length = sizeof(data)/sizeof(u32);

  u32* mem = malloc(sizeof(data) + sizeof(u32));
  memcpy(&mem[1], data, sizeof(data));

  mem[0] = length;
  heron_sqrt(mem);

  int i = 0;
  for (i = 0; i < length; i++) {
    printf("s             = %u\n", data[i]);

    printf("sqrt(s)       = %f\n", sqrt(data[i]));
    printf("heron_sqrt(s) = %d\n", mem[i+1]);
    puts("\n");
  }
  return 0;
}