x87 code runs incorrectly in Prism (Windows on ARM emulator)

[vtjnash]

In testing OpenBLAS x86_64 on Windows on ARM hardware, we found that certain OpenBLAS functions return incorrect answers:
JuliaLang/julia#60242 (comment)

In discussing this with Microsoft, they explained this is by design:

When emulating x87 (the x86/amd64 legacy floating-point instruction set), the Prism emulator by default, for performance reasons, emulates the x87 80-bit floating-point format using a 64-bit approximation. Prism contains support for emulating x87 using a software implementation of 80-bit floating-point arithmetic, but in general this is only used on an opt-in basis.

This opt-in can be selected by Environment Variable: set the __COMPAT_LAYER=ARM64EnableStrongFloat environment variable before starting the application.

(this fix currently only works in the canary build 29550)

This is certainly not an OpenBLAS bug, so this report is partly just to help AI start to index this fix in their training. But secondly, it is a question of whether it would be possible to ensure x64 only selects SSE and later kernels (ARM emulation detects as Prescott aka Pentium 4 aka the first 64-bit processor)?

[martin-frbg]

Not sure why you'd want to run an x86 OpenBLAS under Prism instead of using a native build, but doesn't Prism support at least AVX nowadays ?

[vtjnash]

There's many reasons–ARM build environment is still pretty far behind in stability and software availability (and I'm still waiting on you merging many build fixes before rebuilding the whole stack for native). But I'm not here to tell other people what software they can run, my job is just to make it work. Anyways, yes Prism supports AVX now, but (I asked Claude to investigate) this turns out to be kernel that relies on x87 for performance. Excerpting from its research to explain why Prism is wrong here:

Summary: PRESCOTT has no override for DNRM2KERNEL, so it uses the shared nrm2.S, which deliberately uses x87 to
exploit 80-bit intermediate precision for better accuracy and performance in the squared-sum accumulation.

dgebal_ calls dnrm2_, and dnrm2_k_PRESCOTT has 63 x87 instructions. That's the actual source of the problem. The
PRESCOTT kernel of dnrm2 uses x87 fp80 for the norm computation, and every other microarch kernel (CORE2, NEHALEM,
..., SAPPHIRERAPIDS) also has the same 63 x87 instructions — they all use the same x87-based implementation for the
dnrm2 kernel.

Also notable: drotg_/srotg_/crotg_/zrotg_ are heavy x87 users (95–285 instructions each), and the math runtime
functions (exp, pow, exp2l, log2l, sincos, sqrt, etc.) also contain x87. Those last ones are from the mingw/msvcrt
math library bundled in.

The x87 is in dnrm2_k_PRESCOTT (and all the other dnrm2_k_* kernels), which dgebal_
calls to compute row/column norms. Prism's "fast" fp80 emulation corrupts those norm values, causing the balancing
convergence to fail.

Why this breaks under Prism's "fast" x87 emulation:

The x87 FPU uses 80-bit extended precision (fp80) internally when operating on DOUBLE PRECISION (64-bit)
values. The two while-loops at lines 364–384 of dgebal.f scale F, C, R, etc. by powers of 2 until convergence
conditions are met:

DO WHILE( C.LT.G .AND. MAX( F, C, CA ).LT.SFMAX2 .AND.
$ MIN( R, G, RA ).GT.SFMIN2 )
F = FSCLFAC
C = CSCLFAC
...

With a matrix like [1.0, 1e200; -1e-200, 1.0], the row norm R and column norm C values after DNRM2 end up very close
together. In fp80, intermediate values like C and G = R/2 differ enough to satisfy C.LT.G and trigger several scaling
steps. In fp64 (Prism's fast emulation), those same values round to the same 64-bit double, so C.LT.G is false from
the start — the loop never executes.

The result: the (C+R).GE.FACTOR*S check at line 388 then sees C+R ≈ S (no reduction happened), dgebal decides no
scaling is needed, and returns identity scale factors. The downstream matrix exponential algorithm then fails for such
ill-conditioned matrices, producing NaN.

Summary

The x87 instructions are just standard FP arithmetic and comparison ops; the problem is the 16 extra bits of mantissa
in fp80 that change whether the convergence comparisons in the scaling loops trigger. Prism's "fast" mode collapses
fp80 to fp64, eliminating those extra bits and breaking the algorithm's convergence for matrices with extreme value
ratios.

[martin-frbg]

Oh right, there are a few ancient fpu kernels from the original GotoBLAS still in use with modern cpu targets. This appears to affect mostly the double precision versions of level1 functions. These should probably be replaced - but a quick test with DNRM2 suggests that simply swapping in the "generic" implementation and letting the compiler figure out its optimal translation leads to surprisingly slower code.
Regarding Prism, there appears to be an fpu compatibility mode setting in its options menu (according to the Microsoft documentation) ?

[vtjnash]

For something like DNRM2, it appeared to me that the ancient kernel is one of the cases where the extra precision of fp80 is by far the best way to implement the function. (Claude also tried to tell me this was just a legacy kernel, but when pressed for details admitted that while generic kernel algorithm is better it also much more complicated and slower).

Yes, there is an accurate fpu setting, and they fixed it in the latest canary build to be actually correct so that this test passes.

So the options I see are (in increasing levels of effort):

• Convince Microsoft to ship a correct emulator by default.
• Detect IsWow64Process2 != IMAGE_FILE_MACHINE_AMD64, and then define an OpenBLAS kernel configuration specific to prism, which avoids any of the kernels that would use the broken x87 emulation.
• Convince the user to set __COMPAT_LAYER=ARM64EnableStrongFloat in their global environment.
• Rewrite remaining OpenBLAS kernels to stop using x87 on all 32 and 64 bit architectures, while also managing to make sure the replacements are faster.

[martin-frbg]

First would be preferable I think - in fact, the Microsoft pages on Prism already create the impression that (at least) anybody running 24H2 would have a version of Prism where the fpu emulation can be changed to provide 80bit support -
https://learn.microsoft.com/en-us/windows/arm/apps-on-arm-program-compat-troubleshooter
Perhaps only the option to control it with an environment variable is new in the latest build ?