[jmh530]

Fixes #5305.

core.math.ldexp currently constructs subnormal outputs by shifting the significand, discarding bits without rounding. For example, ldexp(0.609375, -1071) returns four minimum-double units even though the exact result is 4.875 units and should round to five.

This change rounds subnormal results while preserving the existing ordinary-result path.

For float and double formats, express the result as a normal floating-point number of minimum-subnormal units, round with nearbyint, then construct the output encoding. This respects the active rounding mode, preserves negative zero, and handles a carry into the smallest normal value. Using nearbyint avoids introducing an inexact flag where the previous integer-shift path raised none.

Keep the helper out of line so targets that implement nearbyint as a library call do not add call setup to ordinary results. Wider formats use an exact normal intermediate followed by a final scaling multiplication.

Add runtime regressions for both signs, exact and inexact subnormal results, ties-to-even, underflow to signed zero, and rounding into the smallest normal value. Exercise the generic real implementation as well as the public path.

Validation

• Full core.math source-module unittests pass with LDC 1.43 on Windows and Linux at -O0 and -O2, including Linux's generic 80-bit path.
• 5,440 independently calculated boundary cases pass on Windows and Linux, covering all four rounding modes.
• Approximately 199,000 finite float/double comparisons per rounding mode agree with Linux scalbn; overflowing outputs were excluded.
• x86-32 and AArch64 unittest compilation succeeds. These targets were not executed.
• The new regression fails against the previous implementation.

Validation used the complete source module under a temporary module name with installed runtime libraries, rather than a complete compiler/runtime rebuild.

Performance

In a paired Windows benchmark, subnormal-result calls took approximately 3 ns versus 2.3 ns for the previous truncating implementation. A multiplication-only prototype took approximately 40 ns.

The ordinary-result path retains its original instruction sequence. These measurements describe the tested workload and machine, not a general performance guarantee.

Dependency

The signed-zero fix from #5304 is merged. This branch is rebased onto current master and contains only the subnormal-rounding change; it preserves the final signed-zero implementation and tests from #5304.

AI-Assisted-By: Codex (GPT-6 Astra Medium)

[jmh530]

The macOS arm64 failure appears unrelated to this change. It occurs during mimalloc’s CMake compiler check: the linker’s TAPI reader rejects arm64e.x1 in the selected macOS 27 SDK, before the modified D code is built.

[thewilsonator]

Can you rebase it and see if it goes away?

[jmh530]

Replayed the commit and force-pushed to trigger fresh CI. The branch was already based on the latest master (a207d95), so the patch is unchanged. Let’s see whether the macOS arm64 failure recurs.

[jmh530]

@thewilsonator The fresh run still hits the macOS arm64 failure during mimalloc’s CMake compiler check: the linker rejects arm64e.x1 in the macOS 27 SDK.

Both Windows jobs also fail during prerequisite setup because the downloaded curl archive lacks curl-*/bin/curl-ca-bundle.crt. All three failures occur before the modified D code is built.

[kinke]

Does this follow any algorithm/implementation you've found, or is that just AI coming up with one?

After a superficial glance, it doesn't look like it'd support big endian targets (based on the various *cast(Bits*)&value).

[kinke]

Also, it doesn't look like it'd really consider all 112 mantissa bits for the quadruple format.

[jmh530]

Disclosure: Codex did all the work on this PR, including the implementation and testing. It did not explicitly follow a reference implementation. This response and the follow-up analysis below were generated by GPT-6 Astra Pro.

I cannot attribute this particular helper to an established implementation.

The follow-up investigation found a related approach in fdlibm’s scalbn: construct an exact intermediate in the normal range, then multiply by a power of two to produce the rounded subnormal result. That is close to the strategy used by this PR’s wider-format path. This reference was identified afterward.

For the binary32/binary64 helper, the intended calculation is based on subnormal values being integer multiples of the smallest positive subnormal. Let (q) be that smallest value and (y=x,2^n) the exact mathematical result we want. Then:

[
k=\operatorname{round}_{\text{current mode}}\left(\frac{y}{q}\right),
\qquad
\text{result}=kq.
]

For example, if the exact result is (4.875q), round-to-nearest should produce (5q). The helper changes the exponent to express the result in these units, rounds that signed value to an integer, and constructs the output encoding directly. This avoids arithmetic with a subnormal result and preserves the significand until the rounding step. For deeper underflow, the clamp preserves the sign and substitutes a nonzero magnitude below half a unit, which gives the same rounded answer in all four modes.

On endianness, this helper reinterprets the whole binary32/binary64 representation as a same-sized uint/ulong, enforced by the existing static assertion. With matching integer and floating-point byte order, the masks and shifts work on both big- and little-endian targets. The wider-format path retains the existing endian-dependent word indexing. The module also cross-compiled for big-endian AArch64 and s390x, although those binaries were not executed.

The follow-up analysis did find a separate issue with the active-rounding claim. LLVM documents ordinary llvm.nearbyint as assuming the default rounding environment. Using the unchanged PR implementation under a temporary module name, LDC 1.43.0 / LLVM 22.1.8 with -O2 -flto=full eliminated a repeated helper call across a change from nearest to upward rounding. The test used runtime input, and the upward result incorrectly reused zero instead of returning the smallest positive subnormal. pragma(inline, false) did not prevent that optimization.

That needs to be resolved, with a regression covering repeated calls across rounding-mode changes. I also need to document the derivation above and clarify the representation assumptions in the source comments.

[kinke]

AI only, too bad.

The wider-format path retains the existing endian-dependent word indexing.

Where? I only see Bits casts for the lower 64-bit part, the upper 64 bits seem totally ignored. No usage of MANTISSA_MSB etc.

[jmh530]

Also, it doesn't look like it'd really consider all 112 mantissa bits for the quadruple format.

This response was generated by GPT-6 Astra Pro.

The quoted helper is only called under static if (T.mant_dig <= 53). Quadruple precision has T.mant_dig == 113—112 stored fraction bits plus the implicit leading bit—so it takes the separate wider-format branch.

In that branch, the ulong access only updates the exponent-containing half of the value. It preserves the 48 fraction bits in that half and leaves the lower 64 fraction bits untouched. The final vf * (T.epsilon / 4) operates on the complete quadruple value, so all 112 stored fraction bits are retained until that multiplication performs the subnormal rounding. The inspected AArch64 LLVM IR also shows the final operation as fmul fp128.

I could make the helper’s format restriction explicit in its declaration or comments, and add a regression where the lowest fraction bits affect the rounded result.

[kinke]

The quoted helper is only called under static if (T.mant_dig <= 53).

Oh right, I was mislead by the tests and the generic bool, thinking it'd test this too, not just ldexpImpl() (instead of the x87 asm version). Okay, so no endianness and > 64-bit problems.

I don't feel too comfortable with a pure AI solution though, even though the test coverage should be good, as this primitive is used quite often. I hoped you (the human behind the agent) would be into numerics stuff and could verify, ideally after comparing with solutions in well-established libs. E.g., musl's double-precision version seems very compact: https://git.musl-libc.org/cgit/musl/tree/src/math/scalbn.c

[jmh530]

This response was generated by GPT-6 Astra Pro.

Thanks for clarifying, and for the musl reference. I understand the concern about introducing a custom algorithm into such a frequently used primitive. The testing and follow-up analysis described so far were performed by AI; they don’t provide the human verification you’re asking for.

The current helper was intended to avoid the cost of subnormal arithmetic while preserving the previous exception-flag behavior. Those goals need to justify the additional complexity. A direct comparison with musl, covering the numerical reasoning, behavior, and performance, seems like the right next step before deciding whether to retain this approach.

[jmh530]

Added a follow-up fixing an optimization-sensitive rounding issue: under full LTO, the earlier implementation could reuse a result across rounding-mode changes. The helpers now use strict floating-point semantics and constrained rounding/conversion operations. The exact runtime-input reproducer is preserved as a regression test.

Algorithm and byte order

The units-and-integer-encoding approach was developed for this PR; it is not copied from musl. Comments above and inside ldexpSubnormal now explain its scaling, conversion bound, and carry into minimum normal.

For float/double, reinterpreting the complete same-width integer uses the target’s byte order and does not inherently assume little endian. Wider layouts retain explicit word selection.

Comparison with arithmetic alternatives

We compared encoding with a musl-style D implementation and a narrower arithmetic helper that changes only subnormal handling. With matched LDC 1.43.0 / LLVM 22.1.8 builds:

• On the tested AMD EPYC 9V74, the narrower arithmetic helper used roughly 13–24% less time in generic builds, while encoding won with SSE4.1, AVX, and native targeting.
• On the tested Intel Core Ultra 7 258V, that helper took roughly 14–22 times as long as encoding across the tested builds.

Correctness checks passed before benchmarking. Within each build, candidates used the same compiler options, inputs, calling context, and floating-point environment.

These are synthetic measurements on specific virtualized machines, not vendor-wide guarantees. They support retaining encoding rather than selecting arithmetic as the generic fallback.

Exception behavior: maintainer preference

This implementation preserves the affected float/double path’s existing suppression of new underflow/inexact flags. The arithmetic alternatives set those flags for tested inexact tiny results, including some values that round up to minimum normal.

The documentation and existing tests do not establish suppression as a mandatory API requirement, so I’ve retained it for compatibility in this fix. Would you prefer arithmetic exception semantics instead? If so, the implementation and performance tradeoff would need reassessment.

Selecting between the current encoding and arithmetic helpers based on target features would also make exception behavior depend on the build target. That would need an explicit policy decision.

Validation

The new Intel comparison passed:

• 468,416 numerical checks.
• 9,984 focused boundary checks.
• The expected failure of the original implementation’s minimal full-LTO reproducer.

Production regressions cover rounding-mode changes, both signs, signed zero, preservation of existing exception flags, and neighboring values around the minimum-normal boundary.

All RUN commands in the two added regression files passed, including O0/O2, full/thin LTO, separate compilation, and core.math source-module unittests. The expanded production tests also passed native full LTO.

This was focused validation with installed libraries, not a full compiler/runtime rebuild.

Wider-format limitation

Binary128 dynamic rounding still depends on the linked multiplication runtime. Constrained IR preserves the operation’s rounding requirements through optimization, but cannot repair a runtime implementation that always rounds to nearest.

The identified GCC 13 libgcc provider passed 24 emulated AArch64 little-endian cases. A separately linked big-endian executable using LLVM 22.1.8 compiler-rt’s fixed-nearest __multf3 failed 8 of 24 cases.

This change therefore does not claim universal binary128 dynamic-rounding support. It requires a verified compatible runtime provider. No new physical AArch64 execution or exhaustive wider-format validation is claimed.

Attribution

GPT-6 Astra Pro supplied the investigation and AMD/hybrid comparison. Codex implemented the correction and regressions and performed the Intel follow-up.

AI-Assisted-By: Codex (GPT-6 Astra Medium) and Chat-GPT (GPT-6 Astra Pro)

[jmh530]

Following up on the implementation and testing above: all CI checks now pass, including the platforms that initially failed.

The new tests exposed a few existing platform issues: reversed rounding-mode constants on Windows, exception-raising operations being removed by optimization, and an undersized floating-point environment structure on FreeBSD. Those fixes are in separate commits. I also corrected the test guard for WASI.

The exception checks now cover the five standard floating-point exceptions explicitly. They still verify underflow/inexact suppression and preservation of existing flags, but exclude the extra x86 denormal-operand flag that the existing input-normalization code can raise. The minimal runtime-input regression is unchanged.

The macOS tests now pass too, although I haven’t conclusively identified the cause of the earlier flag mismatch. I added failure diagnostics to make any recurrence easier to investigate.

The selected encoding approach and proposed exception policy remain unchanged; the green CI results provide broader platform validation of that implementation.

AI-Assisted-By: Codex (GPT-6 Astra Medium)