Re: [hatari-devel] FPU update

[Andreas Grabher <andreas.grabher@xxxxxxxxxxxx>]

Hello all,

while working to implement FGETEXP using SoftFloat routines i discovered an issue which seems to affect denormalized extended (80-bit) precision floats throughout all functions (add, mul, div, ...).

Background (description simplified, NaN, Infinity and zero not considered):
Every float is made up of a sign (s), an exponent (e) and a fraction (f). 32-bit single precision and 64-bit double precision floats have an implicit integer part of usually 1. The bias of the exponent is 0x7f (single) or 0x3ff (double). 

(-1) ^ s x 2 ^ (e-1023) x 1.f      (double)

This is true as long as the exponent is greater than 0. If the exponent is 0, the implicit integer part also becomes 0 and the bias changes to 0x7e (single) or 0x3fe (double). That is called a denormalized number.

(-1) ^ s x 2 ^ (e-1022) x 0.f      (double)

This means when crossing the border between normalized and denormalized numbers in fact the exponent would make a jump from 1 to -1, if the bias stayed the same. In SoftFloat this jump seems to be handled in normalizeFloatXXSubnormal by returning the exponent as 1 - shiftcount. That adds one to the exponent of all denormalized numbers and closes the gap that is caused by the change of bias. This works for single and double precision BUT

the 80-bit (extended) format does have an explicit integer part (j). Normalized and denormalized numbers are documented to have the same bias of 0x3fff in the Motorola data sheets. Numbers are denormalized, if the integer part is 0. A number with an exponent of 0 can still be normalized.

(-1) ^ s x 2 ^ (0-16383) x j.f      (extended)

Obviously this is not handled that way by SoftFloat, which has the same 1 - shiftcount in normalizeFloatx80Subnormal. This leads to following issue:

Experiment:
floatx80 f, two;
f.high = 0x0000;
f.low = 0x1000000000000000;
two.high = 0x4000;
two.low = 0x8000000000000000;

for (int i = 0; i < 6; i++) {
     f = floatx80_mul(f, two);
     printf("%s (%04x, %16llx)\n",fp_print(&f),f.high,f.low);
}

Result:
+4.20262892889011688e-4933D (0000, 2000000000000000)
+8.40525785778023377e-4933D (0000, 4000000000000000)
+3.36210314311209351e-4932 (0001, 8000000000000000)
+6.72420628622418701e-4932 (0002, 8000000000000000)
+1.34484125724483740e-4931 (0003, 8000000000000000)
+2.68968251448967481e-4931 (0004, 8000000000000000)


You can see the obvious jump in the exponent between the second and third iteration of the loop where the number becomes normalized. It causes the result to be off by factor 2. It seems the behavior of the 6888x is different from the x87 and SoftFloat and also different from what IEEE specifies for denormals: (-1) ^ s x 2 ^ (1 - b) x 0.f    (b = bias)

Fixing this would require some changes to SoftFloat and means, that the native float implementation with casts fails for denormals on all machines that handle it the x87 way.

I'm a little bit surprised about this discovery. I first thought there was something wrong with SoftFloat. Did you ever hear about this difference between x87/IEEE and 6888x? It seems this behavior does not comply with the IEEE standard. Or might this even be a mistake in the data sheets (i checked different versions of it, they all tell the same)?

Best wishes,

Andreas


Am 11.01.2017 um 19:18 schrieb Andreas Grabher <andreas.grabher@xxxxxxxxxxxx>:

> As Toni said, the merge is still work in progress. Most changes are not yet in WinUAE. That is also the cause for your compiler warnings and the unused headers.
> 
> In the meantime i have worked on FMOD and FREM. After some testing i think we have now a quite good SoftFloat implementation for them. Comparing it to math.h freml, remainderl and remquol shows identical results with in range numbers, even with weird inputs like 1e32 % 5. With invalid inputs remquol returns 0, same with SoftFloat. So that should be OK too.
> The native implementation now gives correct remainders, but bad quotients for invalid inputs and for input values with big differences in magnitude (difference of Exponent > 63). That would be hard to fix. Another reason to keep working on SoftFloat. Anyway the new native float implementation is much better than the old one, which gave wrong remainders for numbers with big differences in magnitude and returned no quotient at all.
> 
> Now all arithmetic functions are implemented, except FSINH, FETOXM1, FTANH, FATAN, FASIN, FATANH, FETOX, FTWOTOX, FTENTOX, FCOSH and FACOS.
> 
> 
> Am 11.01.2017 um 18:48 schrieb Nicolas Pomarède <npomarede@xxxxxxxxxxxx>:
> 
>> Le 11/01/2017 à 18:05, Toni Wilen a écrit :
>>>> Hi
>>>> 
>>>> I committed Andreas's changes that Toni merged into WinUAE to Hatari
>>>> devel version ; note that in Hatari there's no setting yet to switch
>>>> between "normal" and softfloat version, but at least everything compiles
>>>> fine (ie softfload is not used yet)
>>> 
>>> That was not complete merge. It does work, mostly but it also lacks some
>>> important fixes. (I forgot to mention it in commit message)
>>> 
>>> The rest will be merged later this week.
>>> 
>> 
>> Thanks, I will wait for the rest of the changes.
>> 
>> 
>> 
> 
> 
>