Mic92/heron_sqrt
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

correct fix-point arithmetic

Browse Source
This commit is contained in:
Jörg Thalheim 2013-12-18 10:01:06 +01:00
parent 588ee8b710
commit df3750364a
1 changed files with 28 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

40
sqrt.c
View File

@ -6,9 +6,12 @@
typedef uint32_t u32; typedef uint32_t u32;
#define COMMA 12
#define COMMA_FIX (COMMA >> 1)
void heron_sqrt(u32 *mem) { void heron_sqrt(u32 *mem) {
int i = mem[0]; int i = mem[0];
while (i > 0) { while (i >= 1) {
u32 s = mem[i]; u32 s = mem[i];
if (s > 1) { if (s > 1) {
// x_0 = (s + 1) / 2 // x_0 = (s + 1) / 2
@ -24,23 +27,29 @@ void heron_sqrt(u32 *mem) {
} }
old_x = x; old_x = x;
} }
mem[i] = x; mem[i] = x << COMMA_FIX;
} }
i--; i--;
} }
} }
int main() { int main() {
u32 data[] = {0, float data[] = {
0,
1, 1,
2, 2,
4, 2.5,
2.25,
2.125,
2.0625,
2.03125, // precision test
3<<2, 3<<2,
1<<3, 1<<3,
1<<4, 1<<4,
1<<5, 1<<5,
1<<6, 1<<6,
5<<2, 5<<2,
11<<7,
1<<7, 1<<7,
1<<8, 1<<8,
1<<16, 1<<16,
@ -50,21 +59,26 @@ int main() {
1<<20, // 32 - 14 = 20 1<<20, // 32 - 14 = 20
1<<24, // 32 - 8 = 24 1<<24, // 32 - 8 = 24
1<<30, 1<<30,
0xffffffff}; 0xffffffff
const int length = sizeof(data)/sizeof(u32); };
const int length = sizeof(data)/sizeof(float);
u32* mem = malloc(sizeof(data) + sizeof(u32)); int i;
memcpy(&mem[1], data, sizeof(data));
u32* mem = malloc((length + 1) * sizeof(u32));
mem[0] = length; mem[0] = length;
for (i = 1; i <= length; i++) {
mem[i] = (u32) (data[i] * (1 << COMMA));
}
heron_sqrt(mem); heron_sqrt(mem);
int i = 0;
for (i = 0; i < length; i++) { for (i = 0; i < length; i++) {
printf("s = %u\n", data[i]); if (data[i] < (0xfffffff >> (COMMA))) {
printf("s = %f\n", data[i]);
printf("sqrt(s) = %f\n", sqrt(data[i])); printf("sqrt(s) = %f\n", sqrt(data[i]));
printf("heron_sqrt(s) = %d\n", mem[i+1]); printf("heron_sqrt(s) = %f\n", ((float) mem[i]) / (1 << COMMA));
}
puts("\n"); puts("\n");
} }
return 0; return 0;