## Summary The SSE2 `fpc_exp_real` in the `{$else FPC_HAS_TYPE_EXTENDED}` branch of `rtl/x86_64/math.inc` computes the integer part of its argument reduction with `addsd 0.5` followed by `cvttsd2si`. Because `cvttsd2si` truncates **toward zero**, that only implements round-to-nearest for `z >= -0.5`. For negative arguments `k` comes out one step too high and the reduced fraction `|f|` lands **outside `[-0.5, 0.5]`**, which is the interval the Padé approximant is valid on. On win64 this makes `exp(-1)` wrong by **13,445,713 ulp** (2.0e-9 relative) — incorrect from the 10th significant digit. Positive arguments are unaffected. This is a regression: on these targets 3.2.2 used the fdlibm implementation in `rtl/inc/genmath.inc` (≤1 ulp). The new routine defines `FPC_SYSTEM_HAS_EXP`, and `genmath.inc` guards its generic version with `{$ifndef FPC_SYSTEM_HAS_EXP}`, so it is no longer compiled in. **Replacing `addsd`+`cvttsd2si`+`subsd` with a single `cvtsd2si` takes the worst case from 13,445,713 ulp to 103 ulp and removes two instructions**, so it should not cost any of the speed the routine was introduced for. ## System Information - **Operating system:** Windows 10 (19045) — **affected**. Debian 13 (glibc 2.41, kernel 6.12) — **not affected**. - **Processor architecture:** x86-64 - **Compiler version:** `3.3.1-20608-gc5192fd1c5` [2026/07/28] (trunk). Also tested: `3.2.2` (not affected). - **Device:** Computer The affected code path is selected by `Extended = Double`, i.e. `FPC_HAS_TYPE_EXTENDED` not being defined — it is not a Windows-specific condition. I have not been able to test AArch64 or other targets that take the same branch. ## Steps to reproduce Compile and run the self-contained program below (no dependencies beyond `SysUtils`): ``` fpc -O2 exp_check.pas ./exp_check ``` The program prints the compiler version, target and `SizeOf(Extended)` itself, so its output is self-describing. ## Example Project ```pascal program exp_check; { Minimal accuracy check for System.Exp on Double. The reference values are the correctly rounded Double of exp(x), embedded as IEEE-754 bit patterns so nothing depends on how a literal is parsed. } {$mode objfpc}{$H+} uses SysUtils; type TCase = record x : Double; ref : QWord; // correctly rounded exp(x) tag : string[10]; end; const CASES: array[1..6] of TCase = ( (x: -1.0; ref: QWord($3FD78B56362CEF38); tag: 'exp(-1)'), (x: -0.5; ref: QWord($3FE368B2FC6F960A); tag: 'exp(-0.5)'), (x: -2.0; ref: QWord($3FC152AAA3BF81CC); tag: 'exp(-2)'), (x: -5.0; ref: QWord($3F7B993FE00D5376); tag: 'exp(-5)'), (x: -10.0; ref: QWord($3F07CD79B5647C9B); tag: 'exp(-10)'), (x: 1.0; ref: QWord($4005BF0A8B145769); tag: 'exp(1)') ); function Bits(const d: Double): QWord; begin Result := PQWord(@d)^; end; var FS: TFormatSettings; // '.' as decimal separator, so output is comparable // between users regardless of their locale function Show(const d: Double): string; begin Result := FloatToStrF(d, ffExponent, 17, 3, FS); end; var i: Integer; got: Double; gb, d, worst: QWord; begin FS := DefaultFormatSettings; FS.DecimalSeparator := '.'; WriteLn('FPC ', {$I %FPCVERSION%}, ' target ', {$I %FPCTARGETCPU%}, '-', {$I %FPCTARGETOS%}, ' SizeOf(Extended)=', SizeOf(Extended)); WriteLn; WriteLn('case Exp() returned bits (hex) expected bits ulp'); worst := 0; for i := Low(CASES) to High(CASES) do begin got := Exp(CASES[i].x); gb := Bits(got); if gb > CASES[i].ref then d := gb - CASES[i].ref else d := CASES[i].ref - gb; if d > worst then worst := d; WriteLn(Format('%-9s %s %.16x %.16x %d', [CASES[i].tag, Show(got), gb, CASES[i].ref, d])); end; WriteLn; WriteLn('worst ulp error: ', worst); if worst > 1000 then WriteLn('>>> AFFECTED (exp() of negative arguments is inaccurate here)') else WriteLn('>>> ok'); end. ``` ## What is the current bug behavior? FPC 3.3.1 (trunk), win64: ``` FPC 3.3.1 target x86_64-Win64 SizeOf(Extended)=8 case Exp() returned bits (hex) expected bits ulp exp(-1) 3.67879441917829E-001 3FD78B5636FA1989 3FD78B56362CEF38 13445713 exp(-0.5) 6.06530659713873E-001 3FE368B2FC6FC1A8 3FE368B2FC6F960A 11166 exp(-2) 1.35335283238742E-001 3FC152AAA3C0AD71 3FC152AAA3BF81CC 76709 exp(-5) 6.73794700152901E-003 3F7B993FE038502E 3F7B993FE00D5376 2817208 exp(-10) 4.53999298450679E-005 3F07CD79B61E727A 3F07CD79B5647C9B 12187103 exp(1) 2.71828182845909E+000 4005BF0A8B1457D0 4005BF0A8B145769 103 worst ulp error: 13445713 >>> AFFECTED (exp() of negative arguments is inaccurate here) ``` ## What is the expected (correct) behavior? FPC 3.2.2, same source, same flags, same machine: ``` FPC 3.2.2 target x86_64-Win64 SizeOf(Extended)=8 case Exp() returned bits (hex) expected bits ulp exp(-1) 3.67879441171442E-001 3FD78B56362CEF38 3FD78B56362CEF38 0 exp(-0.5) 6.06530659712633E-001 3FE368B2FC6F960A 3FE368B2FC6F960A 0 exp(-2) 1.35335283236613E-001 3FC152AAA3BF81CC 3FC152AAA3BF81CC 0 exp(-5) 6.73794699908547E-003 3F7B993FE00D5376 3F7B993FE00D5376 0 exp(-10) 4.53999297624848E-005 3F07CD79B5647C9A 3F07CD79B5647C9B 1 exp(1) 2.71828182845905E+000 4005BF0A8B14576A 4005BF0A8B145769 1 worst ulp error: 1 >>> ok ``` ## Relevant logs and/or screenshots ### Scope — measured on three configurations | | 3.2.2 · win64 | 3.3.1 · win64 | 3.3.1 · x86_64-linux | |---|---|---|---| | `SizeOf(Extended)` | 8 | 8 | 10 | | routine selected | fdlibm (`genmath.inc`) | **new SSE2 (`x86_64/math.inc`)** | x87 (`x86_64/math.inc`) | | ulp on `exp(-1)` | 0 | **13,445,713** | 0 | | wrong out of 12 | 2 (1 ulp each) | **11** | **0** | **x86_64-linux is not affected**, verified by running there (Debian 13), cross-compiled from the same source with the same compiler: 12/12 bit-exact. There `FPC_HAS_TYPE_EXTENDED` is defined, so the x87 routine is selected instead. That was confirmed directly rather than inferred: a test program carrying a **verbatim copy** of the SSE2 routine as inline assembler reproduces the RTL output **12/12 on win64**, but does **not** reproduce it on Linux (1/12, and that one is a coincidence) — i.e. a different routine is compiled there. ### The error tracks the reduced fraction exactly | input | abs(f) produced by the code | error | |---|---|---| | `exp(-0.001)` | 0.0014 (in range) | 1 ulp | | `exp(-0.5)` | 0.72 | 11,166 ulp | | `exp(-2)` | 0.89 | 76,709 ulp | | `exp(-5)` | 1.21 | 2,817,208 ulp | | `exp(-10)` | 1.43 | 12,187,103 ulp | | `exp(-1)` | 1.44 | 13,445,713 ulp | Re-implementing the routine in double-precision Python using the `trunc(z+0.5)` rule reproduces FPC's current output **bit-for-bit on all twelve inputs**, which confirms the diagnosis independently of the assembler test. ### Things that are ruled out - **Not the optimizer**: identical results with `-O3` and with `-O1 -g -Criot`. - **Not the process FP environment**: `System.Ln` is bit-exact (5/5) in every one of those binaries. - **Not a local patch**: the sequence is present in the official `main` branch of `rtl/x86_64/math.inc`. (My compiler is an `fpcup` build and reports as `-dirty`, so I checked the upstream source directly.) ### Downstream impact `Math.Power(base, exponent)` with a non-integer exponent is `exp(exponent*ln(base))`, so every non-integer power inherits the error — measured 2.99e-9 relative worst case over 20,000 random pairs. Anything computing `exp` of negative numbers is hit hardest: error functions, normal distributions, log-likelihoods, softmax, continuous discounting. In my case an `erf` implemented from CPython's published algorithm produced a normal CDF with 1.4e-9 absolute error where ~1e-16 was expected. What makes it easy to miss is that the same program is correct on Linux and wrong on Windows, with the same compiler and the same source. ## Possible fixes The offending sequence, in `rtl/x86_64/math.inc`: ```asm .LNotZero: mulsd %xmm4, %xmm0 addsd %xmm2, %xmm0 { Add 0.5 to make sure the fractional part falls between -0.5 and 0.5 } { %xmm0 won't be out of range due to the earlier checks } cvttsd2si %xmm0, %rax subsd %xmm2, %xmm0 cvtsi2sd %rax, %xmm1 ``` The comment states the intended invariant; the code does not establish it for negative `z`. `cvtsd2si` rounds using the current MXCSR rounding mode (round-to-nearest by default), which is exactly the reduction the algorithm needs: ```diff .LNotZero: mulsd %xmm4, %xmm0 - addsd %xmm2, %xmm0 { Add 0.5 to make sure the fractional part falls between -0.5 and 0.5 } - { %xmm0 won't be out of range due to the earlier checks } - cvttsd2si %xmm0, %rax - subsd %xmm2, %xmm0 + { cvtsd2si rounds with the current MXCSR rounding mode (round-to-nearest by + default), which is the reduction this algorithm needs. "addsd 0.5 + + cvttsd2si" only rounds to nearest for z >= -0.5, because cvttsd2si + truncates toward zero; for z < -0.5 it left abs(f) > 0.5, outside the + interval the Pade approximant is valid on. } + cvtsd2si %xmm0, %rax cvtsi2sd %rax, %xmm1 ``` (`FPC_HALF_DOUBLE` is then only needed for the `%xmm2` load at `.LNoUnderflow`, which can also be dropped.) ### Verification of the fix I did not rebuild the RTL. Instead the routine was copied verbatim into a test program together with a patched copy, both as inline assembler, and run natively. **The verbatim copy reproduces the RTL output bit-for-bit in 12/12 cases**, so it is measuring the same code: ``` case RTL ulp verbatim copy patched exp(-1) 13445713 13445713 103 exp(-0.5) 11166 11166 2 exp(-0.75) 985832 985832 2 exp(-2) 76709 76709 0 exp(-5) 2817208 2817208 3 exp(-10) 12187103 12187103 71 exp(-20) 54660 54660 5 exp(-0.001) 1 1 1 exp(1) 103 103 103 exp(0.5) 0 0 0 exp(2) 1 1 1 exp(10) 73 73 73 worst: 13,445,713 ulp -> 103 ulp ``` ### Where it came from, and one thing the patch does not fix The routine was added in `baa6f2a1` ("x86-64: New optimised 'exp' routine", MR !989), whose goal was speed — a 30–45% improvement. The merge request notes about accuracy: > "The least-significant bit of the mantissa is sometimes different compared to manually-calculated > values, but it is acceptable for this single bit to be off since it is at the limit of > floating-point precision and correctness cannot be guaranteed." That criterion seems entirely reasonable to me. My point is only that the routine does not appear to meet it for negative arguments: one bit would be 1 ulp, and what I measure is 1.3e7. It looks like the accuracy check simply did not cover negative inputs. With the reduction fixed, the worst case is **103 ulp, not ≤1 ulp**. That residue is the accuracy of the [5/5] Padé itself near the edge of the interval — it appears identically on `exp(1)`, which was never affected by the truncation. It is worth noting that on this same architecture both other code paths do better: the x87 routine gives 0 ulp on all 12 inputs on Linux, and the `genmath.inc` fdlibm code that this routine displaced on win64 gave ≤1 ulp. Whether to raise the approximant's degree, add a correction step, or not shadow `genmath.inc` here is a design call for whoever wrote it — I am reporting it as an observation, not proposing a change.