72 lines
1.2 KiB
C
72 lines
1.2 KiB
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];
|
|
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) {
|
|
mem[i] = x;
|
|
break;
|
|
}
|
|
old_x = x;
|
|
}
|
|
mem[i] = x;
|
|
}
|
|
}
|
|
}
|
|
|
|
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;
|
|
}
|