[LCP]Rounding questions

[Chuck Martin]

On Fri, Sep 06, 2002 at 07:21:46AM +1000, Greg Black wrote:
> Steve Baker wrote:
> | well).  If you look in /usr/include/features.h you'll find all the "features"
> | you can select from.
> 
> And only if your system has features.h and all the rest of
> panoply.  None of the BSD systems I have access to here has
> this.

If this is true, I probably don't want to use round() then.

> Given that the functionality of round() is trivial to implement
> with the C89 floor() and ceil() functions, it would be much
> simpler to implement it in terms of those.

That was the way this program worked originally, but then I found
a bug that I haven't been able to figure out, and I thought maybe
using round() would be a better solution.  Maybe you can shed some
light on this bug.  Here is the relevant part of the code as it was
originally:

 	case ROUND:
	    {	int prec = (int) eval(e->e.o.right);
		double	scal = 1;
		if (0 < prec)
		    do scal *= 10; while (0 < --prec);
		else if (prec < 0)
		    do scal /= 10; while (++prec < 0);

		if (rndtoeven)
		    return (rint(eval(e->e.o.left) * scal) / scal);
		else {
		    double temp = eval(e->e.o.left);
		    temp *= scal;
/* xxx */
		    temp = ((temp - floor(temp)) < 0.5 ?
			    floor(temp) : ceil(temp));
		    return (temp / scal);
		}
	    }

I only added the "/* xxx */" line, which I'll explain shortly.  In the
above snippet, rndtoeven is a flag that determines how to handle the
case where the fractional part of a number is .5.  When set, it rounds
to the nearest even integer.  When unset, it should round away from
zero, as the round() function does.  The eval() function evaluates a
mathematical expression which is stored in a tree structure of operators
and values, and returns a double.  e->e.o.right is the precision and
e->e.o.left evaluates to the number to be rounded.

Most of the time, this works correctly, however I came across a case
when it didn't.  In this case, e->e.o.right (the precision) is 2 and
e->e.o.left is an expression which multiplies 390 and 0.0765.  This
product is 29.835, so the rounding should result in a return value of
29.84.  Instead, it is returning 29.83.

At the point indicate by "/* xxx */", I inserted an fprintf() that
printed the value of temp to stderr, which I then redirected to another
virtual console.  When I saw that this looked correct, I inserted
another fprintf() that printed the value of floor(temp) to stderr.
I could have just modified the first fprintf(), but I didn't.  In any
case, it still looked correct, so I inserted yet another fprintf()
that printed the result of "(temp - floor(temp)) < 0.5".  At this
point, not only did everything look correct, but the return value
suddenly changed to 29.84.  Deleting any one of the fprintf()'s
caused the erroneous value to return.  It didn't matter which one I
deleted.  All three had to be present for the return value to be
correct.

I decided to find out if the calculations being performed in the
fprintf()'s were effecting the result by replacing all three fprintf()'s
with fflush(stdout)'s, and the return value was still correct, but
a little experimenting showed that now, only two of the lines were
necessary.  Deleting one of the fflush(stdout)'s didn't cause the
result to be incorrect, but deleting two of them (leaving only one)
did.  Furthermore, if I change the initial assignment of temp to
"double temp = 390 * 0.0765", the result is correct even if all of
the fflush(stdout)'s are removed, however, it isn't evaluating the
expression it's supposed to be evaluating, so that isn't a proper
fix.  Unfortunately, I don't know how to isolate the problem to
simplify this into a much smaller program you can experiment with.  
The only way I know that you can experiment with this bug is to
take the program as a whole, which you're welcome to do.  It's
in the public domain, although you might need the version I'm
currently working on to see the bug.

As I write this, I've tried a few more variations of the fprintf()'s,
trying various combinations of temp and floor(temp), and the results
definitely depend not just on how many of them are present, but on
just what is being printed.  I just haven't been able to identify a
pattern yet.  It seems that the difference between temp and floor(temp)
can change between exactly 0.5 and 0.499999999999944 without either of
them changing individually, and this change will effect all of the
fprintf()'s at the same time even if it's a later one that's modified.
Any suggestions are welcome.

Chuck