Complex matrix inversion with issues

Discussion around the SwissMicros DM42 calculator
Thomas Okken
Posts: 1100
Joined: Tue May 02, 2017 5:48 pm
Location: Netherlands
Contact:

Re: Complex matrix inversion with issues

Post by Thomas Okken »

It looks like hypot() looks for a special case where one of the arguments is sufficiently small relative to the other, returning the larger argument, and performs this check by comparing the arguments' exponents while ignoring their mantissas. In floating-point formats where the mantissa has a dynamic range of 1 digit of the base, that would work fine, but with BID, the range is 10^34...

You can see this with 1 ENTER 30 1/X 3 * →POL

Or even more obvious,

1 ENTER 0.01 ENTER 1e-36 + →POL

Note that 0.01 + 1e-36 returns 0.01; representing that number exactly would require 35 digits, and of course we have only 34. And yet the calculation above returns a different result than that calculation without the "ENTER 1e-36 +" part. The only difference being the denormal form of 0.01 being used.

I'll take a look at the Intel source code later and see if I can come up with a patch.
whuyse
Posts: 198
Joined: Thu Dec 21, 2017 1:23 pm

Re: Complex matrix inversion with issues

Post by whuyse »

Fascinating.
Shouldn't Intel come up with a patch?
I think you are not using the latest version of the lib, because the differences didn't apply to Free42 - perhaps the newer version does have a correct hypot()?
Werner
41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE, DM15L
Thomas Okken
Posts: 1100
Joined: Tue May 02, 2017 5:48 pm
Location: Netherlands
Contact:

Re: Complex matrix inversion with issues

Post by Thomas Okken »

I don't think Intel is doing anything with this anymore. I guess they accept patches, or at least the 2.0u2 update looks like it was the result of patches provided by someone outside Intel.

I compared 2.0u1 and 2.0u2 file by file, and the differences consisted exactly of what was mentioned in the release notes, nothing more, so the hypot() bug must be in the newer version as well.
Thomas Okken
Posts: 1100
Joined: Tue May 02, 2017 5:48 pm
Location: Netherlands
Contact:

Re: Complex matrix inversion with issues

Post by Thomas Okken »

There you go:

bid128_hypot.c, lines 129-134:

Code: Select all

    if(exponent_x - exponent_y >= 35) 
    {   
     res.w[BID_HIGH_128W] = x.w[BID_HIGH_128W];
     res.w[BID_LOW_128W] = x.w[BID_LOW_128W];
     BID_RETURN (res);
    }
...without normalizing the numbers first.

The complete diff between the 2011 (2.0u1) and 2018 (2.0u2) versions of this file:

Code: Select all

2c2
<   Copyright (c) 2007-2018, Intel Corp.
---
>   Copyright (c) 2007-2011, Intel Corp.
Thomas Okken
Posts: 1100
Joined: Tue May 02, 2017 5:48 pm
Location: Netherlands
Contact:

Re: Complex matrix inversion with issues

Post by Thomas Okken »

I'm adding a patch to the Intel library, to change the number 35 in the offending check to 68. That covers the worst-case scenario, where the effective difference between the two exponents can be 33 less than the actual difference, if the larger number has a mantissa less than 10 while the smaller number has a mantissa greater than or equal to 10^33. That fixes the problem and removes the need for the work-around in Phloat::hypot().

I'll alert Intel to the bug but I'm not submitting a patch. If the exponent magnitude check is desired, they may be able to make it robust without preventing the performance optimization in legitimate cases like my fix does, but I'll leave that headache to them.
jneven
Posts: 12
Joined: Sat Jan 23, 2021 6:32 pm
Location: Belgium

Re: Complex matrix inversion with issues

Post by jneven »

Great work, Thomas.

I am looking to upgrade the DM42 to this new version of Intel-BCD. The library for the DM42 in SwissMicros-git is called gcciiilibbid_hard.a. I was wondering if this can be replaced by the library created via your shell script build-intel-lib.sh for arm7 arch.

Greetings,

Johan.
Thomas Okken
Posts: 1100
Joined: Tue May 02, 2017 5:48 pm
Location: Netherlands
Contact:

Re: Complex matrix inversion with issues

Post by Thomas Okken »

I haven't tried building the DM42 firmware myself, and I don't know if the library that's built by their build script is compatible with the one built by mine for armv7 on Android. It might be, but I think it would be safer and easier to just modify their version of the patch applied to the Intel library before building it like I did with mine.

If you take a look at the various *.patch files in my repo, you'll see that my latest commit modified all of them in exactly the same way: by appending a 19-line patch to LIBRARY/src/bid128_hypot.c, which changes the number 35 on line 129 to 68. Nothing else should be needed.
jneven
Posts: 12
Joined: Sat Jan 23, 2021 6:32 pm
Location: Belgium

Re: Complex matrix inversion with issues

Post by jneven »

Thanks for the great info, Thomas.
I tried to build a gcc111bcd library with the arm-tools, but I'm not shure what are the options taken by Swiss Micros, so I took another method to minimize the chance for trouble. I took the lib from the SM git repo, took bid128_hypot.o out of it. After disassembly (arm...objdump) of this object and your unpatched and patched object, I could see where to change this "35." (34 in object) into "68." (67 in object) . With xxd I made a hexdump, changed this number and reversed this into a new object. After insertion into SM gcc111bcd-lib I did a rebuild of the SM free42- and DM42PGM tree and all was OK.
So on my DM42, this odd matrix is now inverted without the workaround and with your patch in the intel-libary.
Again many thanks for your great software and help.
Johan.
Thomas Okken
Posts: 1100
Joined: Tue May 02, 2017 5:48 pm
Location: Netherlands
Contact:

Re: Complex matrix inversion with issues

Post by Thomas Okken »

jneven wrote:
Sun Mar 07, 2021 4:17 pm
I took the lib from the SM git repo, took bid128_hypot.o out of it. After disassembly (arm...objdump) of this object and your unpatched and patched object, I could see where to change this "35." (34 in object) into "68." (67 in object) . With xxd I made a hexdump, changed this number and reversed this into a new object. After insertion into SM gcc111bcd-lib I did a rebuild of the SM free42- and DM42PGM tree and all was OK.
Whoa. Nice!!! 8-)
Thomas Okken
Posts: 1100
Joined: Tue May 02, 2017 5:48 pm
Location: Netherlands
Contact:

Re: Complex matrix inversion with issues

Post by Thomas Okken »

I reported the bug to Intel, and I got a response saying they'll patch the code much like I did, except they're not changing the 35 on line 129 to 68, but to 69. I'm making that change in my patches as well.
Post Reply