[LCP]Questions about precision

[Chuck Martin]

Thanks for the reply.  I had begun to think no one wanted to tackle
this one at all.  :)

On Tue, Oct 15, 2002 at 03:34:46PM +0100, David Spencer wrote:
> It looks like it rounds according to the first undisplayed digit (FUD - 

This much I had figured out.

> Dealing with exact floating point numbers in C is a bit of a pain.  One 
> approach is to store an integer and exponent/scale, so you store 9075 
> and 10^3.  When it comes to displaying the represented value, you can 
> display I/S (9), a '.', and I-S*floor(I/S), not forgetting to use the 
> appropriate format for after the dot, e.g: (I'm storing E and S 
> separately for convenience.)

This is all well and good if you already know the size of numbers you're
working with.  When you're working with arbitrary input from the user
that could be any size, and have to also deal with intermediate results
that may not fit so nicely in this scheme, it becomes more difficult.

Just for the record, the program I've been working on as I've brought
up these questions on rounding and precision is the public domain
spreadsheet program, sc, which I've taken upon myself to maintain for
the last three years, and which is why I was looking at how Lotus and
Excel dealt with this.

> If Excel works to 30 significant digits then it almost certainly doesn't 
> use the standard C library.  It might use an alternative FP library, or 
> perhaps it stores the numbers differently, like text or BCD.

When I hadn't seen any responses to my previous message, I decided to
experiment with Excel again a couple of days ago, and I discovered that
Excel lies about what is actually being stored and used internally.
It seems that when displaying a number, it first rounds it to fifteen
significant digits, starting from the leftmost non-zero digit (the 9 to
the left of the decimal point in 9.075), which in this case makes it
appear to be exact.  It then pads the right with zeroes, which I think
is a bit deceptive, because it still uses the non-exact value internally.
This can be seen by entering "=9.075-9" into a cell, which effectively
removes two significant non-zero digits from the left, and forces Excel
to include two more digits to the right before rounding, resulting in
"0.074999999999999300000000000000" being displayed if the cell is
formatted to display 30 decimal places.  Multiply this result by 100
and subtract 7, and Excel will display another significant digit, and
repeating this several times shows that the decimal part of the number
is actually the same as that which you and I both got in our tests.

The fact that Excel rounds up instead of down means that the rounding
is either being done on the string itself instead of the double that
is obviously being used for other calculations, or they're using some
other trick.

In any case, I'd rather not emulate that behavior.  I think it would
be better to let the user know exactly what is being stored and used
internally, and now that I know what Excel is doing, if anyone e-mails
me asking why sc isn't as exact as Excel, I'll tell them that it is,
and that Excel is really lying to them about what value is really being
used.

> One way of rounding to an integer is to add a half and round down. 
> Adapting this, and adding .005 could work:
> 
> dave at feathers:~/prog/c > !ca
> cat float.c
> main()
> {
>  printf("%.20f %.20f %.20f\n",9.075+0.005,9.075,0.005);
> }
> 
> dave at feathers:~/prog/c > !./
> ./float
> 9.08000000000000007105 9.07499999999999928946 0.00500000000000000010

This is an interesting result.  If you add these values by hand (the
second and third numbers), you actually get 9.0799* because the .005 is
off by less than the 9.075, but I just confirmed your result in sc.

> Even if 0.005 were interpreted as 0.0049999999 this would still work to 
> 2dp as the value would be 9.079*.

I'm not sure I understand this.  If you "add a half and round down" as
you said, and if the result of adding a half was 9.079*, then rounding
down would give 9.07.

Anyway, I've come up with what I think is a nice solution to rounding
these inexact numbers in the proper direction, and I'd like to hear any
comments on this, especially if anyone sees any potential problems with
it.

What I've done is basically to scale the number being rounded, just as
it was being done before, so that the rounding is being done to the
nearest integer, but then, instead of comparing the scaled number to
the floor of the scaled number + .5 to determine the direction to round,
I take the floor of the scaled number, add .5, and scale it back to the
original scale, and then compare this unscaled midpoint to the original
number before scaling to determine the rounding direction.  If that seems
confusing or unclear, here is the actual code (e->e.o.left is the number
being rounded, and e->e.o.right is the number of decimal places to round
it to, with negative numbers meaning to round to the left of the decimal
point):

	case ROUND:
	    {	int precision = (int) eval(e->e.o.right);
		double temp, val, midpoint, fl, cl;
		double scale = 1;

		temp = val = eval(e->e.o.left);
		if (0 < precision)
		    do scale *= 10; while (0 < --precision);
		else if (precision < 0)
		    do scale /= 10; while (++precision < 0);
		temp *= scale;
		fl = floor(temp);
		cl = ceil(temp);
		midpoint = (fl + 0.5) / scale;

		if (rndtoeven)
		    if (val == midpoint)
			return (fl == floor(temp/2)*2 ? fl : cl) / scale;
		    else
			return (rint(temp) / scale);
		else
		    return (val < midpoint ? fl : cl) / scale;
	    }

Any comments or suggestions are welcome, but my preliminary tests seem to
give the desired results.

Chuck