Excuse me, I attached attached yesterday's code to today's email.
Here's the updated code.
Regards,
Daniel
/*
* Analyze error of reciprocal approximations
* (c) 2005 Daniel Phillips
*/
#include <stdio.h>
#include <math.h>
#define table_size 1024
#define sample_bits 18
#define frac_bits 2
int main(void)
{
double neg_error = 0, pos_error = 0;
int table[table_size + 1];
int precision = (1 << sample_bits);
int i, j;
for (i = 0; i <= table_size; i++)
table[i] = (1 << sample_bits)/(1 + ((float)i / table_size));
for (i = 0; i < table_size; i++) {
int sample = table[i], delta = sample - table[i + 1];
printf("%i: [0x%x] error: ", i, sample);
for (j = 0; j < 1 << frac_bits; j++) {
#if 1
// more precise: keep extra interpolation bits around
int interp = (sample << frac_bits) - j*delta + (1 << (frac_bits - 1));
double crude = ((double)interp)/ (precision << frac_bits);
#else
// less precise: truncate product to same length as sample
int interp = sample - ((j*delta) >> frac_bits);
double crude = ((double)interp + .0) / precision;
#endif
double exact = 1/((double)((i << frac_bits) + j) /
(double)(table_size << frac_bits) + 1);
double error = crude - exact;
printf("%.2e ", error);
if (error > pos_error)
pos_error = error;
if (error < neg_error)
neg_error = error;
if ((int)(error*precision))
goto quit;
}
printf("\n");
}
quit:
printf("\nmax error = [%e, %e]\n", neg_error, pos_error);
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)
