heron_sqrt/c/sqrt.c

93 lines
1.6 KiB
C
Raw Normal View History

2013-11-11 09:26:27 +00:00
#include <stdio.h>
#include <math.h>
#include <stdint.h>
2013-12-04 20:32:13 +00:00
#include <stdlib.h>
2013-11-11 09:26:27 +00:00
typedef uint32_t u32;
2014-01-02 12:01:07 +00:00
// have to a multiply of 2
2014-01-13 08:08:21 +00:00
#define COMMA 8
2013-12-18 09:01:06 +00:00
#define COMMA_FIX (COMMA >> 1)
2013-12-04 20:32:13 +00:00
void heron_sqrt(u32 *mem) {
2013-12-13 09:36:19 +00:00
int i = mem[0];
2014-01-02 12:01:07 +00:00
while (i > 0) {
2013-12-04 20:32:13 +00:00
u32 s = mem[i];
2014-01-02 12:01:07 +00:00
u32 x = s;
2013-12-06 14:11:41 +00:00
if (s > 1) {
// x_0 = (s + 1) / 2
// without increment to avoid overflow for 0xffffffff
2014-01-02 12:01:07 +00:00
x = x >> 1;
2013-12-06 14:11:41 +00:00
u32 old_x = x;
2013-11-11 09:26:27 +00:00
2013-12-04 20:32:13 +00:00
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;
}
2013-11-11 09:26:27 +00:00
}
2014-01-02 12:01:07 +00:00
mem[i] = x << COMMA_FIX;
2013-12-13 09:36:19 +00:00
i--;
2013-11-11 09:26:27 +00:00
}
}
int main() {
2013-12-18 09:01:06 +00:00
float data[] = {
2014-01-13 08:08:21 +00:00
0,
2014-01-02 12:01:07 +00:00
//0.5,
1.0 / (1 << COMMA),
1.0 / (1 << (COMMA - 1)),
2013-12-04 20:32:13 +00:00
1,
2,
2013-12-18 09:01:06 +00:00
2.5,
2.25,
2.125,
2.0625,
2.03125, // precision test
2013-12-06 14:11:41 +00:00
3<<2,
2013-11-11 09:26:27 +00:00
1<<3,
1<<4,
1<<5,
1<<6,
2013-12-06 14:11:41 +00:00
5<<2,
2013-12-18 09:01:06 +00:00
11<<7,
2013-11-11 09:26:27 +00:00
1<<7,
1<<8,
1<<16,
1<<17,
1<<18,
1<<19,
1<<20, // 32 - 14 = 20
1<<24, // 32 - 8 = 24
1<<30,
2013-12-18 09:01:06 +00:00
0xffffffff
};
const int length = sizeof(data)/sizeof(float);
int i;
2013-11-11 09:26:27 +00:00
2013-12-18 09:01:06 +00:00
u32* mem = malloc((length + 1) * sizeof(u32));
2013-12-04 20:32:13 +00:00
mem[0] = length;
2013-12-18 09:01:06 +00:00
for (i = 1; i <= length; i++) {
2014-01-02 12:01:07 +00:00
// float to fixed point
2013-12-18 09:01:06 +00:00
mem[i] = (u32) (data[i] * (1 << COMMA));
}
2013-12-04 20:32:13 +00:00
heron_sqrt(mem);
2013-11-11 09:26:27 +00:00
2013-12-04 20:32:13 +00:00
for (i = 0; i < length; i++) {
2014-01-02 12:01:07 +00:00
// out-range-check
if (data[i] <= (0xffffffff >> COMMA)) {
printf("s = %10.10f\n", data[i]);
printf("sqrt(s) = %10.10f\n", sqrt(data[i]));
// fixed point to float
printf("heron_sqrt(s) = %10.10f\n", ((float) mem[i]) / (1 << COMMA));
2013-12-18 09:01:06 +00:00
}
2013-12-06 14:11:41 +00:00
puts("\n");
2013-11-11 09:26:27 +00:00
}
return 0;
}