2008-10-03 05:43:44.
// Show effect of optimisation of fp results.
// Code taken from ~/src/tcl/tcl8.6a2/libtommath/bn_mp_sqrt.c
// Shows bug in 'if (d != 0.0)', with d being inexact
// Compile and run with -
// make CC=/usr/local/gcc-4.2.1/bin/gcc CFLAGS=-O2 LDFLAGS=-lm isqrt && ./isqrt
//
#include <math.h>
#include <stdio.h>
int main()
{
double d = 0.0;
int i;
unsigned long dig;
for (i = 0; i < 2; i++) {
d = ldexp(d, 28) + (double)((1<<28)-1);
//printf("d=%g=%a\n", d, d);
}
d=sqrt(d);
dig = (unsigned long) ldexp(d, -28);
printf("d=%g=%a dig=%d\n", d, d, dig);
d -= ldexp((double) dig, 28);
printf("d=%g=%a dig=%d\n", d, d, dig);
}