On Tuesday 15 February 2005 11:17, Lourens Veen wrote:
> Be sure to try every possible input for the mantissa. It gets
> interesting...
Today's code takes advantage of the ieee packing code to iterate through
252 * 2^16 distinct inputs, giving max error = [-3.76e-06, 3.68e-06],
or +-1/2^18, using the simple minded approach. I'll try improving that
later this week.
This code uses two tables as we discussed, a slight simplification.
There is an exponent problem in the range 2^-127, which I expected and
2^126, which I didn't expect and need to track down. If somebody beats
me to it, I'd like that.
The overflow in the zeroth table entry is handled subtly wrong. It
works for up to 9 interpolated bits before my crude hack breaks down.
This version is probably already good enough for our hardware, but I'll
keep working to knock the fuzz off. What's wrong with it is likely the
same sort of thing that will bite you when you graft the exponent
handling onto yours. Note, we also need to normalize the input because
the interpolants aren't normalized.
Regards,
Daniel
/*
* Floating point reciprocal by table interpolation
* (c) 2005 Daniel Phillips
* License: GPL version 2
*/
#include <stdio.h>
#include <math.h>
#define table_size_bits 10
#define sample_bits 18
#define interp_bits 6
#define table_size (1 << table_size_bits)
#define interp_mask ((1 << interp_bits) - 1)
#define sample_mask ((1 << sample_bits) - 1)
#define precision_bits (table_size_bits + interp_bits)
#define precision_mask ((1 << precision_bits) - 1)
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 value[table_size];
int delta[table_size];
float recip(float x)
{
int sign, exp, man, i;
unpack_ieee(x, &sign, &exp, &man);
man >>= (23 - precision_bits);
i = man >> interp_bits;
pack_ieee(&x, sign, -exp - !!man,
((value[i] - (((man & interp_mask)*delta[i]) >> interp_bits))
<< (23 - sample_bits)));
return x;
}
int calc_sample(int i)
{
return (1 << (sample_bits + 1)) / (1 + ((float)i/table_size));
}
int main(void)
{
int i;
for (i = 0; i < table_size; i++) {
value[i] = calc_sample(i) & sample_mask;
delta[i] = (value[i] - calc_sample(i + 1)) & sample_mask;
}
float x, neg_err = 0, pos_err = 0;
int start = 0, count = 252 << 16;
int sign = 0, exp = -126, man = 0;
for (i = start; i < start + count; i++) {
pack_ieee(&x, sign, exp, man);
float r = recip(x), r32 = 1/x, e = (r - r32)/r32;
#if 0
if (!(i % 99))
printf("[%x] 1/%f = %f (%f) err=%e\n", *(int *)&x, x, r, r32, e);
#endif
if (e < neg_err)
neg_err = e;
if (e > pos_err)
pos_err = e;
if (!(man = (man + (1 << (23 - precision_bits))) & 0x7fffff)) {
++exp;
printf("exp = %i\n", exp);
}
}
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)
