Core: pldexp could split 2^e into three factors instead of four

pldexp builds 2^e from four factors because an integer e/3 was once too slow to compute, but e arrives as a floating-point packet, where an approximation of e/3 is one multiplication. Splitting into three factors, b = \operatorname{trunc}(e \cdot \tfrac{89}{256}), drops one multiplication from the dependent chain: on x86 (SSE2, AVX2+FMA) pldexp latency falls 15–22% and throughput improves 1–25%. A prototype exists locally for pldexp_generic and the SSE/AVX/AVX-512 double specializations, and one formula covers half, bfloat16, float and double. This issue records the method, the alternatives for computing b, the costs, a pre-existing double-rounding bug that the prototype reduces but does not fix, and what applying it everywhere would take.

Current method

pldexp_generic (arch/Default/GenericPacketMathFrexpLdexp.h, 4cb563a0, reordered in !2706 (merged)) clamps e to \pm(2^{E} + M - 1), which is 278 for float and 2099 for double: beyond that every finite input over- or underflows. It then computes

b = \lfloor e/4 \rfloor = e \gg 2, \qquad a \cdot 2^e = \big(((a \cdot 2^b) \cdot 2^b) \cdot 2^b\big) \cdot 2^{e-3b}

where each factor is built by shifting b + \text{bias} into the exponent field. 2^{278} needs at least three normal floats. The 2020 comment explains the choice of four: "it turns out to be faster to split 'exponent' into four factors, since [exponent>>2] is much faster to compute than [exponent/3]." Every partial product contains a, so no intermediate over- or underflows where the result does not. !2527 (merged) regrouped the chain, and !2706 (merged) had to restore the order.

The SSE Packet2d, AVX Packet4d, AVX-512 Packet8d, AltiVec Packet2d and ZVector Packet2d specializations each carry their own copy of the same four-factor split.

Proposed method

e is already a floating-point packet, so an approximation of e/3 costs one multiplication, and the conversion to integer truncates it:

b = \operatorname{trunc}\big(e \cdot \tfrac{89}{256}\big), \qquad a \cdot 2^e = \big((a \cdot 2^b) \cdot 2^b\big) \cdot 2^{e-2b}
const Packet clamped = pmin(pmax(exponent, neg_max_exponent), max_exponent);
const PacketI e = pcast<Packet, PacketI>(clamped);
const PacketI b = pcast<Packet, PacketI>(pmul(clamped, pset1<Packet>(Scalar(0.34765625))));  // trunc(e*89/256)
const PacketI b_remainder = psub(e, padd(b, b));                                              // e - 2b
const Packet c1 = preinterpret<Packet>(plogical_shift_left<MantissaBits>(padd(b, bias)));      // 2^b
const Packet c2 = preinterpret<Packet>(plogical_shift_left<MantissaBits>(padd(b_remainder, bias)));  // 2^(e-2b)
return pmul(pmul(pmul(a, c1), c1), c2);

What b must satisfy. Below, c denotes the multiplier in b = \operatorname{trunc}(e \cdot c) or \lfloor e \cdot c \rfloor, an approximation of 1/3. b need not be e/3. It must satisfy two conditions for every e in the clamp range:

  1. Order: b and e - 2b have the sign of e, or are zero, so the partial products move monotonically from a to the result (equivalently, b lies between 0 and e/2).
  2. Range: 2^b and 2^{e-2b} are normal.

Condition 1 is what trips up floor-based formulas at small negative e: at e = -1, \lfloor e \cdot c \rfloor = -1 gives b = -1 and e - 2b = +1.

Why 89/256 and not 1/3. Half has no margin: at e = -41, b must be exactly -14 (then -14 - 14 - 13 = -41), and \operatorname{trunc}(-41/3) = -13 puts 2^{-15}, a subnormal, in the last factor. An exhaustive check over every e shows that b = \operatorname{trunc}(e \cdot c) meets both conditions for half, bfloat16, float and double exactly when c \in [0.3415, 0.3535], whether the conversion truncates or rounds to nearest. 89/256 = 0.34765625 is near the middle of that window and exact in all four formats.

type e_{\max} \max \lvert b\rvert \max \lvert e-2b\rvert normal exponents
half 41 14 14 [−14, 15]
float, bfloat16 278, 262 96 86 [−126, 127]
double 2099 729 641 [−1022, 1023]

Every backend's float-to-int conversion truncates: x86 cvtt*, NEON vcvtq_s32_f32, AltiVec vec_cts, SVE FCVTZS, LSX vftintrz, RVV vfcvt_rtz. Round to nearest would also do. The x86 double specializations use _mm*_cvttpd_epi32, so the MXCSR rounding mode doesn't matter.

SSE2 Packet2d. A direct port needs a second int32→int64 swizzle, which cost 3% on small exponents. Interleaving instead, unpacklo_epi32(b, e - 2b) plus one bias add, puts each pair in one 64-bit lane. slli_epi64(v, 52) then yields 2^b, since the high half shifts out, and slli(srli(v, 32), 52) yields 2^{e-2b}.

Ways to compute b

Six approaches were tried. All were checked exhaustively over the clamp range for float and double; the ones that pass are bit-identical to std::ldexp on the test bases.

  1. \lfloor e/4 \rfloor, four factors (master). One arithmetic shift, but the extra dependent multiplication.
  2. Floating point, \operatorname{trunc}(e \cdot c) (proposed). 2 µops: a multiply and a conversion. They start from the float exponent, in parallel with the conversion of e itself. Truncation toward zero gives condition 1 without any constant.
  3. Integer magic number, (21846\,e + 2^{15}) \gg 16. 21846 = \lceil 2^{16}/3 \rceil. The +2^{15} is required, because a plain multiply-high floors and fails condition 1 at e = -1. The cost is a 32-bit pmul<PacketI>: emulated with shuffles on SSE2, and pmulld (2 µops, 10-cycle latency) on SSE4.1/AVX2. It also has to wait for the conversion of e. On AVX2 its small-e throughput is worse than master's four factors, which may be why an earlier fast-division attempt did not pay off.
  4. Shift-add, \lfloor (3e + 4)/8 \rfloor = (e + (e + e) + 4) \gg 3. Condition 2 only needs a multiplier between about 0.27 and 0.45 (float) or 0.26 and 0.49 (double), so 3/8 works: \lvert b\rvert \le 3\lvert e\rvert/8 + 1/2 and \lvert e - 2b\rvert \le \lvert e\rvert/4 + 1. The +4 gives condition 1: for e \ge 4, e - 2b \ge (e-4)/4 \ge 0; for e \le -3, e - 2b \le (e+3)/4 \le 0; the rest were checked directly. It costs 4 integer ops and no multiply.
  5. Power-of-two multipliers are out. (2e + k) \gg s, or any (m e + k) \gg s with m a power of two, has multiplier 2^{-t}. 1/4 leaves e - 2b \approx e/2 (139 at e = 278, over 127); 1/2 makes \lvert b\rvert \approx e/2. A brute-force search over m \le 9, s \le 7 and every k finds no power of two. The cheapest valid forms are (3e + \{3, 4\}) \gg 3 and (5e + k) \gg 4 for k = 5 \dots 10, all 4 ops. The clamp can't be narrowed to widen the window: \text{denorm\_min} \cdot 2^e only overflows at e = 277.
  6. Difference of shifts, (e \gg 1) - (e \gg 3). 3 ops, multiplier 3/8 ((e \gg 1) - (e \gg 4), 7/16, also works). At e = -1 both floors are -1 and cancel, and in general 3e/8 - 1/2 \le b \le 3e/8 + 7/8. A sum such as (e \gg 2) + (e \gg 4) fails condition 1 with any constant from −4 to 4.

None of the integer forms covers half: 4 and 6 fail at e \in \{-41, -40, -39\}. The floating-point form covers it with c \in [0.3415, 0.3535].

Float pldexp_generic, ns per packet, best of 5, GCC 13.3 -O3, i7-13700KF under WSL2. "Small" is e \in [-20, 20]; "full" is the whole clamp range. The magic-number row comes from a separate run, in which the floating-point version measured 1.11 / 15.6 (SSE2) and 1.20 / 25.9 (AVX2).

method ops for b SSE2 small SSE2 full AVX2 small AVX2 full
1. \lfloor e/4 \rfloor, four factors (master) 1 shift 1.48 20.2 1.32 34.0
2. \operatorname{trunc}(e \cdot c) 1 mul + 1 cvt 1.12 15.6 1.21 26.0
3. (21846 e + 2^{15}) \gg 16 1 int mul + 2 1.58 16.2 1.43 26.6
4. (3e + 4) \gg 3 4 1.27 14.7 1.35 24.5
6. (e \gg 1) - (e \gg 3) 3 1.29 14.6 1.40 24.6

Fewer integer ops don't translate into speed on x86:

  • Shared units: vector shifts, multiplies and conversions contend for the same two of the three vector units on this core, and pldexp already spends two shifts building its factors.
  • Op totals: the loop has 4 integer adds/subs with method 2, 6 with method 6 and 7 with method 4. Method 2 instead adds one extra multiply and one conversion on the floating-point side.
  • Full-range gains: the integer forms' advantage there most likely comes from their larger b passing fewer intermediates through subnormals, which each trigger a microcode assist. That depends on the data.
  • Non-x86: on backends where integer shifts don't compete with multiplies (possibly NEON, SVE, RVV), method 6 may win. That needs measuring.

All benchmarks in this issue were taken with c = 1/3; 89/256 generates the same instructions.

Pre-existing: double rounding for subnormal results

Master's pldexp is not exact when the result is subnormal: an intermediate product becomes subnormal, is rounded, and is rounded again by the next multiplication. For example, in half, pldexp(1.25·2^-22, -3):

  • The exact result is 0.625 \cdot 2^{-24}, which rounds to 2^{-24}.
  • 1.25 \cdot 2^{-23} first rounds to 2 \cdot 2^{-24}, and later steps round the tie to 0.

Mismatches against std::ldexp, GCC 13.3 -O2, identical on SSE2 and AVX2+FMA:

case master three factors
float, random bases in [2^{-140}, 2^{-100}], e \in [-60, 0] 5,578 of 200,000 1,126
double, random bases in [2^{-1070}, 2^{-1000}], e \in [-200, 0] 2,583 of 200,000 594
half, every finite base, e \in [-45, 45] (via pldexp_generic<half>) 46,432 of 5,777,408 7,408
float, bases in [2^{100}, 2^{127}], e \in [-278, -200] 0 0

Each wrong result is off by one unit in the last place, and only in the subnormal range. ArrayBase::ldexp (9d5669f4) is documented as exact scaling, and packetmath does not catch this. Three factors mean one fewer intermediate rounding, which removes about 80% of the cases but not all.

An exact version must round only once: bring a to a normal m \in [\tfrac12, 1) with its exponent split off, which is exact. Scale m by the whole exponent t in one multiply when 2^t is normal. Otherwise use m \cdot 2^{t + 126}, which is normal and exact, and then one rounding multiply by 2^{-126} (float). Its cost against the three-factor split is unmeasured; it may belong in its own issue.

Pros and cons

Pros

  • One multiplication fewer on the dependent chain: latency −15% to −22%.
  • One fewer chance of a microcode assist when an intermediate is subnormal, which is visible as −12% to −23% on full-range exponents.
  • About 5× fewer double-rounded subnormal results than master (see above).
  • One formula for half, bfloat16, float and double.
  • Drops the integer pnmadd(3, b, e) of the generic path, a 10-cycle pmulld on SSE4.1/AVX2, and the three-op e - 3b of the double specializations.
  • The chain shrinks from three barriers to two once the fast-math fix (!3189 (merged)) is in.

Cons

  • One extra floating-point multiply and one extra float→int conversion per call, off the critical path. On a port-bound caller this can cost more than it saves: SSE2 double with small exponents went from −1% to +3% until the interleave above.
  • Every hand-written specialization has to be redone and measured separately. Each has its own int32→int64 widening.
  • Only measured on x86.

Benchmarks

benchmarks/Core/bench_packet_ops.cpp with two new cases:

  • Pldexp/range: throughput over 256 packets of bases in [1, 2) and integer exponents in [−range, range], checked against std::ldexp.
  • PldexpLatency: a dependent chain scaling by 2^{37} and back.

GCC 13.3 -O3, Google Benchmark, median of 5, 13th Gen Intel(R) Core(TM) i7-13700KF under WSL2, pinned to one core. Master is 53c33b64. Times are ns per iteration (256 packets, or 2 pldexp calls for latency), ± standard deviation.

ISA case master three factors change
SSE2 PldexpLatency_float 7.1 ± 0.0 5.6 ± 0.0 −22%
SSE2 PldexpLatency_double 7.1 ± 0.0 6.0 ± 0.0 −15%
SSE2 Pldexp_float/20 358.3 ± 1.9 270.1 ± 0.7 −25%
SSE2 Pldexp_float/278 5089 ± 34 4067 ± 16 −20%
SSE2 Pldexp_double/20 340.3 ± 0.7 336.8 ± 7.2 −1%
SSE2 Pldexp_double/2099 988.3 ± 3.9 862.2 ± 1.2 −13%
AVX2+FMA PldexpLatency_float 7.4 ± 0.0 5.8 ± 0.0 −21%
AVX2+FMA PldexpLatency_double 7.4 ± 0.1 5.8 ± 0.0 −22%
AVX2+FMA Pldexp_float/20 302.0 ± 1.9 281.9 ± 1.6 −7%
AVX2+FMA Pldexp_float/278 8201 ± 38 6288 ± 36 −23%
AVX2+FMA Pldexp_double/20 393.1 ± 2.5 366.9 ± 1.2 −7%
AVX2+FMA Pldexp_double/2099 1519 ± 10 1338 ± 3 −12%

A second run of the whole table agreed within 2 percentage points. Plain AVX (the Packet4d specialization without AVX2), measured with an ad-hoc harness on the same machine: latency −25% for both types, full-range throughput −21% (double) and −24% (float). pow(x, 2.5), which ends in pldexp, doesn't change: pldexp is a small part of it. exp uses pldexp_fast and is unaffected.

Correctness. Bit-identical to std::ldexp for every integer e \in [-300, 300] (float) and [-2200, 2200] (double) on ±0, ±denorm_min, 3·denorm_min, ±min, 0.75·min, 1, 1.5, 0.7, 1e-20, 3e15, ±max and ±inf. Checked on SSE2, AVX and AVX2+FMA, with and without -ffast-math. packetmath and packetmath_fastmath pass on those three. AVX-512DQ and AVX512-FP16 (the half path) compile without warnings but were not run (i7-13700KF under WSL2); the half path was exercised by instantiating pldexp_generic on the scalar Eigen::half. "Bit-identical" covers the listed bases; the double-rounding cases above are subnormal results from other bases.

Work to apply it throughout Eigen

  • pldexp_generic: three factors with b = \operatorname{trunc}(e \cdot 89/256) for every type. This covers float everywhere, double on NEON, SVE, RVV, LSX and the clang vector backend, native half on AVX512-FP16, and bfloat16 and half through their float wrappers. (prototyped)
  • SSE Packet2d with the interleaved factor construction; AVX Packet4d and AVX-512 Packet8d with cvttpd. (prototyped)
  • AltiVec Packet2d and ZVector Packet2d: rewrite their own four-factor copies. They also need !3189 (merged)'s barriers; that MR only reaches the generic path and the x86 copies.
  • Choose how to compute b per backend (see Ways to compute b): the floating-point form on x86, and possibly the difference of shifts on NEON, SVE and RVV for float and double. Measure the non-x86 backends before switching them: NEON (aarch64 and 32-bit arm), SVE, RVV, LSX, AltiVec. There the float→int conversion or the extra multiply may cost differently, and the double paths that go through pcast<Packet2d, Packet2l> may not pay off. CI labels: neon-tests, sve-tests, rvv-tests, lsx-tests, altivec-tests (ZVector has no CI job).
  • Run AVX-512 and AVX512-FP16 on hardware.
  • Fix the double rounding of subnormal results (above), in this change or its own, and test pldexp against std::ldexp on subnormal results.
  • Tests: packetmath samples exponents. Add an exhaustive sweep of every integer exponent across the clamp range on special, subnormal and normal bases, so a factor leaving the normal range fails for every type.
  • Benchmarks: add Pldexp and PldexpLatency to bench_packet_ops.cpp. (prototyped)
  • Order relative to !3189 (merged) (fast-math barriers in pldexp_apply_factors): the prototype is stacked on it, and the helper drops from three multiplications to two.
Edited by Charles Schlosser