phi: avoid intermediate overflow in BE scale

src/phi_almohy.jl

@@ -92,6 +92,21 @@ end
     return r
 end
 
+@inline function _phi_log_be_coeff(m::Integer, p::Integer)
+    r = 0.0
+    for j in 1:m
+        r -= log(float(m + p + j))
+    end
+    for j in 1:(m + p + 1)
+        r -= log(float(m + j))
+    end
+    return r
+end
+
+@inline function _phi_be_scale(m::Integer, p::Integer, normTs, delta, K::Integer)
+    return exp((_phi_log_be_coeff(m, p) - delta * log(normTs)) / K)
+end
+
 """
     PhiPadeCache(A, p)
 
@@ -302,7 +317,7 @@ function _phi_ell!(cache::PhiPadeCache, A, normT, m::Integer, p::Integer, phat::
     normTs = normT / scalefac
     delta = (p - 1) * (p - phat) / p + 1
     K = 2m + p + 1
-    c = (_phi_be_coeff(m, p) / normTs^delta)^(1 / K)
+    c = _phi_be_scale(m, p, normTs, delta, K)
     absA = cache.absA
     @. absA = c * abs(A) / scalefac
     alpha = _normpow_nonneg!(cache, absA, K)
@@ -498,7 +513,7 @@ function _phi_ell(A::AbstractMatrix, normT, m::Integer, p::Integer, phat::Intege
     normTs = normT / scalefac
     delta = (p - 1) * (p - phat) / p + 1
     K = 2m + p + 1
-    c = (_phi_be_coeff(m, p) / normTs^delta)^(1 / K)
+    c = _phi_be_scale(m, p, normTs, delta, K)
     scaledT = c .* abs.(A ./ scalefac)
     alpha = _normpow_nonneg(scaledT, K)
     t = log2(2alpha / eps(Float64)) / (K - delta) + t0

test/basictests.jl

@@ -377,6 +377,19 @@ end
         @test P30[j] ≈ ref30[j] rtol = 1.0e-8 atol = 1.0e-13
     end
 
+    for normTs in (1.0e-8, 1.0e-2, 0.5, 1.0, 1.9), p in (60, 200, 500)
+        m = 1
+        phat = 0
+        delta = (p - 1) * (p - phat) / p + 1
+        K = 2m + p + 1
+        logcoeff_big = -sum(log(big(m + p + j)) for j in 1:m) -
+            sum(log(big(m + j)) for j in 1:(m + p + 1))
+        ref = Float64(exp((logcoeff_big - big(delta) * log(big(normTs))) / big(K)))
+        c = ExponentialUtilities._phi_be_scale(m, p, normTs, delta, K)
+        @test isfinite(c)
+        @test c ≈ ref rtol = 5.0e-14
+    end
+
     # Preallocated in-place entry point returns phi_j in out[j+1].
     n = 7
     A = randn(n, n) ./ 3