commit c35e2ca28eb847b2b556c5bc017723d88dc426bc
parent 4eddb4d8cf371b4fb2196183b3fca241021544bd
Author: Jared Tobin <jared@jtobin.io>
Date: Thu, 2 Jul 2026 16:35:05 -0230
perf: log1p for the per-step wealth factor, fast-path log_sum_exp
Two small changes to the per-observation work in the update
functions:
1. Replace 'log (1 + lambda * z)' with 'log1p (lambda * z)' in
all four update paths. Improves accuracy for small 'lambda *
z' (the standard log1p rationale) and avoids the '1 + ...'
add.
2. In the two-sided combinations (Bounded, BernTS), skip the
'log_sum_exp' call when the cheap upper bound
'max(logw_p, logw_n) + log 2' is already at or below the
running max-log-sum. Under H_0 the wealth process is a
supermartingale that stays near its early peak, so this
fast-path fires on the vast majority of steps once the peak
is set.
Impact from criterion on a 1000-sample cycle at p = p_0:
* Bernoulli.TwoSided fold newton: 25.62 us -> 15.65 us (-39%)
* Bernoulli.TwoSided fold fixed: 16.21 us -> 8.27 us (-49%)
* Bounded fold adaptive: 14.09 us -> 12.88 us (-9%)
* Bounded fold newton: 14.83 us -> 14.57 us (-2%)
The BernTS gap closes because its H_0 workload triggers the
fast-path aggressively. Bounded is more modest because the
input stream ('cycle [0.3, 0.7]') keeps both directional
wealths active at similar magnitudes, so the log_sum_exp fires
more often.
Also moves 'log 2' to a shared 'log2_dbl' constant in Common
(used both as the initial two-sided max-log-sum and as the
fast-path slack).
I checked GHC Core for the Newton fold specialisation
question: the current bench's shared 'run_b' helper causes GHC
to emit one worker taking 'Bettor' as an argument, with a
per-iteration case dispatch. Rewriting the bench to inline the
fold at each arm produces the specialised code, but the
measurable difference is in the noise (<1 us on a 15 us
fold). Left the library shape alone.
All 59 tests still pass.
Diffstat:
4 files changed, 38 insertions(+), 18 deletions(-)
diff --git a/lib/Numeric/Eproc/Bernoulli.hs b/lib/Numeric/Eproc/Bernoulli.hs
@@ -75,6 +75,7 @@ module Numeric.Eproc.Bernoulli (
, samples
) where
+import GHC.Float (log1p)
import Numeric.Eproc.Common (
Bettor(..), Verdict(..), ConfigError(..)
, BetState, init_bet, bet_lambda, step_bet
@@ -200,8 +201,7 @@ update Config{..} State{..} !x =
let !xd = if x then 1 else 0
!z = xd - cfg_p0
!lam = bet_lambda cfg_bettor cfg_lam_max st_bet
- !fac = 1 + lam * z
- !logw' = st_log_w + log fac
+ !logw' = st_log_w + log1p (lam * z)
!maxw' = max st_max_log_w logw'
!s' = step_bet cfg_bettor cfg_lam_max st_bet z
in State (st_n + 1) logw' maxw' s'
diff --git a/lib/Numeric/Eproc/Bernoulli/TwoSided.hs b/lib/Numeric/Eproc/Bernoulli/TwoSided.hs
@@ -65,10 +65,11 @@ module Numeric.Eproc.Bernoulli.TwoSided (
, samples
) where
+import GHC.Float (log1p)
import Numeric.Eproc.Common (
Bettor(..), Verdict(..), ConfigError(..)
, BetState, init_bet, bet_lambda, step_bet
- , finite, log_sum_exp
+ , finite, log_sum_exp, log2_dbl
)
-- types ----------------------------------------------------------------------
@@ -156,7 +157,7 @@ initial Config{..} =
st_n = 0
, st_log_w_pos = 0
, st_log_w_neg = 0
- , st_max_log_sum = log 2
+ , st_max_log_sum = log2_dbl
, st_bet_pos = s0
, st_bet_neg = s0
}
@@ -178,12 +179,15 @@ update Config{..} State{..} !x =
!z = xd - cfg_p0
!lam_p = bet_lambda cfg_bettor cfg_lam_max_pos st_bet_pos
!lam_n = bet_lambda cfg_bettor cfg_lam_max_neg st_bet_neg
- !fac_p = 1 + lam_p * z
- !fac_n = 1 - lam_n * z
- !logw_p = st_log_w_pos + log fac_p
- !logw_n = st_log_w_neg + log fac_n
- !log_sum = log_sum_exp logw_p logw_n
- !max_sum = max st_max_log_sum log_sum
+ !logw_p = st_log_w_pos + log1p (lam_p * z)
+ !logw_n = st_log_w_neg + log1p (negate lam_n * z)
+ -- see the twin comment in 'Numeric.Eproc.Bounded.update' for
+ -- why we can skip 'log_sum_exp' via a cheap upper bound.
+ !cheap_ub = max logw_p logw_n + log2_dbl
+ !max_sum
+ | cheap_ub <= st_max_log_sum = st_max_log_sum
+ | otherwise =
+ max st_max_log_sum (log_sum_exp logw_p logw_n)
!sp = step_bet cfg_bettor cfg_lam_max_pos st_bet_pos z
!sn = step_bet cfg_bettor cfg_lam_max_neg st_bet_neg (negate z)
in State (st_n + 1) logw_p logw_n max_sum sp sn
diff --git a/lib/Numeric/Eproc/Bounded.hs b/lib/Numeric/Eproc/Bounded.hs
@@ -86,10 +86,11 @@ module Numeric.Eproc.Bounded (
, samples
) where
+import GHC.Float (log1p)
import Numeric.Eproc.Common (
Bettor(..), Verdict(..), ConfigError(..)
, BetState, init_bet, bet_lambda, step_bet
- , finite, log_sum_exp
+ , finite, log_sum_exp, log2_dbl
)
-- types ----------------------------------------------------------------------
@@ -209,7 +210,7 @@ initial Config{..} =
st_n = 0
, st_log_w_pos = 0
, st_log_w_neg = 0
- , st_max_log_sum = log 2
+ , st_max_log_sum = log2_dbl
, st_bet_pos = s0
, st_bet_neg = s0
}
@@ -242,12 +243,17 @@ update Config{..} State{..} !x =
let !z = x - cfg_null_mean
!lam_p = bet_lambda cfg_bettor cfg_lam_max_pos st_bet_pos
!lam_n = bet_lambda cfg_bettor cfg_lam_max_neg st_bet_neg
- !fac_p = 1 + lam_p * z
- !fac_n = 1 - lam_n * z
- !logw_p = st_log_w_pos + log fac_p
- !logw_n = st_log_w_neg + log fac_n
- !log_sum = log_sum_exp logw_p logw_n
- !max_sum = max st_max_log_sum log_sum
+ !logw_p = st_log_w_pos + log1p (lam_p * z)
+ !logw_n = st_log_w_neg + log1p (negate lam_n * z)
+ -- Skip 'log_sum_exp' when the cheap upper bound
+ -- log_sum_exp a b <= max a b + log 2
+ -- already sits at or below the running max: no update can
+ -- move it. Under H_0 (calibration) this is the common case.
+ !cheap_ub = max logw_p logw_n + log2_dbl
+ !max_sum
+ | cheap_ub <= st_max_log_sum = st_max_log_sum
+ | otherwise =
+ max st_max_log_sum (log_sum_exp logw_p logw_n)
!sp = step_bet cfg_bettor cfg_lam_max_pos st_bet_pos z
!sn = step_bet cfg_bettor cfg_lam_max_neg st_bet_neg (negate z)
in State (st_n + 1) logw_p logw_n max_sum sp sn
diff --git a/lib/Numeric/Eproc/Common.hs b/lib/Numeric/Eproc/Common.hs
@@ -33,6 +33,7 @@ module Numeric.Eproc.Common (
-- * Internal: helpers
, finite
, log_sum_exp
+ , log2_dbl
) where
import GHC.Float (log1p)
@@ -130,6 +131,15 @@ log_sum_exp !a !b
| otherwise = b + log1p (exp (a - b))
{-# INLINE log_sum_exp #-}
+-- | @log 2@ as a shared constant. Used both as the initial value of
+-- the two-sided running max-log-sum (since @K^+_0 + K^-_0 = 2@) and
+-- as the tight upper-bound slack in the fast-path skip inside
+-- 'Numeric.Eproc.Bounded.update' /
+-- 'Numeric.Eproc.Bernoulli.TwoSided.update'.
+log2_dbl :: Double
+log2_dbl = log 2
+{-# INLINE log2_dbl #-}
+
-- | Per-bettor state. One constructor per 'Bettor' alternative; the
-- constructor used in any given state matches the 'Bettor' chosen
-- in the enclosing 'Config'.
