// 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);
}