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

/*
 *        s   1
 *         \ /
 *          +
 *          |   1
 *          \  /
 *           >>
 *           |
 *       s   x
 *        \ / \
 *        (/) |\
 *         |  | \
 *         \ /   \
 *          +     \
 *          |  1  |
 *          \  /  |
 *           >>   |
 *            |  /
 *            \ /
 *             -
 */
typedef uint32_t u32;

u32 heron_sqrt(u32 s, int* iterations) {
  // x_0 = (s + 1) / 2
  // without increment to avoid overflow for 0xffffffff
  u32 x = s >> 1;
  u32 old_x = x;

  *iterations = 0; // statistic

  if (x == 0) return 0;
  while(1) {
    (*iterations)++; // statistic
    // x_{n + 1} = (x_n + (s / x_n)) / 2
    x = (x >> 1) + ((s/x) >> 1);
    if (old_x  <= x) {
      return x;
    }
    old_x = x;
  }
}

int main() {
  u32 data[] = {0, 2, 4,
    1<<2,
    1<<3,
    1<<4,
    1<<5,
    1<<6,
    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};
  int iter = 0, i = 0;
  for (i = 0; i < sizeof(data)/sizeof(int); i++) {
    if (i != 0) printf("------\n");
    printf("s = %u\n", data[i]);

    printf("sqrt(s)      = %f\n", sqrt(data[i]));
    printf("heron_sqrt(s)   = %d (%d iterations)\n",
        heron_sqrt(data[i], &iter), iter);

    // 1/16 = 0.0625
    u32 s8 = data[i] << 8;
    printf("heron_sqrt8(s)  = %f (%d iterations)%s\n",
        ((double) heron_sqrt(s8, &iter)) / 16,
        iter, s8 == 0 ? " Out of range!" : "");

    // 1/128 = 0.0078125
    u32 s14 = data[i] << 14;
    printf("heron_sqrt14(s) = %f (%d iterations)%s\n",
        ((double) heron_sqrt(data[i]<<14, &iter)) / 128,
        iter, s14 == 0 ? " Out of range!" : "");
  }
  return 0;
}