Hi all,
This puts together the reciprocal table interpolation with ieee
pack/unpack functions to create an ieee-compatible reciprocal
approximation suitable for use in the simulation. The input, output
and table widths are generalized. The ieee format is hard-wired to
float32 at the moment. Denormals, nans and overflow aren't handled.
The accuracy estimation code isn't very good in that it supplies test
points that can't actually be represented in the target precision,
which makes accuracy seem a little bit worse than it actually is.
Accuracy is about +-2**-17, using straight interpolation. Later, I'll
try incorporating some fancy rounding to improve accuracy by an
additional bit.
Regards,
Daniel
/*
* Floating point reciprocal by table interpolation
* With
* (c) 2005 Daniel Phillips
*/
#include <stdio.h>
#include <math.h>
#define table_size_bits 10
#define sample_bits 18
#define interp_bits 10
#define table_size (1 << table_size_bits)
#define interp_mask ((1 << interp_bits) - 1)
#define precision_bits (table_size_bits + interp_bits)
inline void unpack_ieee(float f, int *sign, int *exp, int *man)
{
unsigned *p = (int *)&f;
*sign = *p >> 31;
*exp = ((*p >> 23) & 0xff) - 127;
*man = ((*p & 0x7fffff))/*| (1 << 23)*/;
}
inline void pack_ieee(float *f, int sign, int exp, int man)
{
*(int *)f = (sign << 31) | ((exp + 127) << 23) | (man & 0x7fffff);
}
int table[table_size + 1];
float recip(float x)
{
int sign, exp, man;
unpack_ieee(x, &sign, &exp, &man);
man >>= (23 - precision_bits);
int i = man >> interp_bits, sample = table[i], delta = sample - table[i + 1];
pack_ieee(&x, sign, -exp - !!man,
((sample - (((man & interp_mask)*delta) >> interp_bits))
<< (23 - sample_bits)));
return x;
}
int main(void)
{
int i;
for (i = 0; i <= table_size; i++)
table[i] = (1 << (sample_bits + 1))/(1 + ((float)i/table_size));
float neg_err = 0, pos_err = 0;
for (i = 0; i < 10000; i++) {
float f = (float)i/1001 + 1, r = recip(f), r32 = 1/f, e = (r - r32)/r32;
printf("1/%f = %f (%f) err=%e\n", f, r, r32, e);
if (e < neg_err)
neg_err = e;
if (e > pos_err)
pos_err = e;
}
printf("max error = [%e, %e]\n", neg_err, pos_err);
return 0;
}
_______________________________________________
Open-graphics mailing list
[email protected]
http://lists.duskglow.com/mailman/listinfo/open-graphics
List service provided by Duskglow Consulting, LLC (www.duskglow.com)
