From: j@bitminer.ca
Subject: Re: exact floating point calculations in roff(7)
To: Ingo Schwarze <schwarze@usta.de>
Cc: tech@openbsd.org
Date: Thu, 03 Apr 2025 15:10:44 -0400



On 2025-04-03 08:09, Ingo Schwarze wrote:
> Hello John,
> 
> On Wed, Apr 02, 2025 at 02:52:39PM -0400, j@bitminer.ca wrote:
> 
>> Hello, I think you left out the best part:
>> 
>>         switch (unit) {
>>         case 'f':
>>                 myres *= 65536.0;
>>                 break;
>>     ...
>>     <snip>
>>     ...
>>         case 'M':
>>                 myres *= 24.0 / 100.0;
>>                 break;
>>         default:
>>                 break;
>>         }
>>         if (res != NULL)
>>                 *res = myres;    <--- ***here***
>> 
>> The whole point seems to be to calculate an integer,
>> returned in *res.
> 
> That is completely correct.
> The point is that the roff(7) language allows users to specify
> arbitrary floating point numbers and these can be scaled with
> so-called scaling units - but the end result is indeed an integer.
> 
>> Some programmer
> 
> That was me.  :-)

Ahh, not knowing the history here, I tried to be neutral.

No offense intended.

> 
>> who might have been less inclined to write scaled
>> integer arithmetic code
> 
> Before i wrote the double arithmetic code, that code did use integer
> arithmetic.  I had written that integer arithmetic code several years
> earlier.
> 
>> decided, instead, to adopt doubles.
> 
> I deemed the change necessary because the roff(7) input syntax
> supports decimal fractions.
> 
> Admittedly, it might be possible to parse 23.999 as *two* integers
> (23 and 999), scale both parts separately, do some kind of carry,
> and finally add both results, all using integer arithmetics.
> I admit i did not consider that possibility when switching this
> code to floating point, but it does not feel like a simple solution.
> 
>> OK, fine, but I think an implied goal here is to generate an exact
>> *rounded* result.  Which is necessary because as often seen doubles
>> don't calculate integer values very reliably.
>> 
>> While your fix is obviously valid I think the fix could also have
>> been
>> 
>>     *res = trunc(myres + 0.5);
>> 
>> which eliminates the 23.999999 problem in the usual floating-point
>> way.
> 
> No, that would be incorrect behaviour.
> If myres is 23.9 or 23.99 or even 23.999999, then roff(7) semantics
> requires that *res must become 23, but your rounding would make it 24.

So now that makes me confused, as your original example was

> 100M = 100 * 24/100 = 23   # sic, 23 instead of 24

which result was correct/incorrect depending on compiler options.  Is
it supposed to be (100*24)/100 or the likely broken 100*(24/100)?  It
seems roff(7) is way outside my comfort zone after all.

Your test case attempt:

> 	double x, y;
> 	x = 100.0;
> 	y = 24.0/100.0;
> 	x *= y;

will not work -O2 as the compiler reduces the arithmetic to a
load/store of a compile-time computed constant.

I suggest inserting

void
altc(double *a, double *b) {
return;
}

as a separate compilation unit, and adding

void
altc( double *a, double *b);

near the top, and adding

         altc(&x, &y);

before
         x *= y;

which outsmarts the compiler and forces arithmetic even -O2.

As a suggestion, perhaps there should be specific comments
to explain the ordered and simple arithmetic statements.

cheers

J