heron_sqrt/sqrt.c

#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;
  for (i = 1; i <= mem[0]; i++){
    u32 s = mem[i];
    // x_0 = (s + 1) / 2
    // without increment to avoid overflow for 0xffffffff
    u32 x = s >> 1;
    u32 old_x = x;

    if (x != 0) {
      while(1) {
        // x_{n + 1} = (x_n + (s / x_n)) / 2
        x = (x >> 1) + ((s/x) >> 1);
        printf("%d\n", x);
        if (old_x <= x) {
          mem[i] = x;
          break;
        }
        old_x = x;
      }
    }
    mem[i] = x;
  }
}

//void heron_sqrt2(u32 *mem) {
//  int i = 1;
//  while(i <= mem[0]) {
//    u32 s = mem[i];
//    // x_0 = (s + 1) / 2
//    // without increment to avoid overflow for 0xffffffff
//    u32 x = s >> 1;
//    u32 old_x = x;
//
//    if (x != 0) {
//      while(1) {
//        // x_{n + 1} = (x_n + (s / x_n)) / 2
//        //x = (x >> 1) + ((s/x) >> 1);
//        x = (x / 2) + ((s/x) / 2);
//        printf("%d\n", x);
//        if (old_x == x) {
//          mem[i] = x;
//          break;
//        }
//        old_x = x;
//      }
//    }
//    mem[i] = x;
//    i++;
//  }
//}

int main() {
  u32 data[] = {0,
    1,
    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};
  const int length = sizeof(data)/sizeof(u32);
  const int mem_size = sizeof(data) + sizeof(u32);
  u32* mem = malloc(mem_size);
  memcpy(mem, data, mem_size);
  u32* mem2 = malloc(mem_size);
  memcpy(mem2, data, mem_size);

  mem[0] = length;
  heron_sqrt(mem);
  //heron_sqrt2(mem2);

  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]);
    //printf("heron_sqrt2(s)   = %d\n", mem2[i]);
  }
  return 0;
}