new implementation of logaddexp (#41)

* new implementation of logaddexp

* Add threshold for zero results in `log1pexp`

* changed type promotion to pass the new tests

* clarified logaddexp signature

* relaxed signature type to Real

* removed inaccurate threshold for logaddexp

* removed internal function _logaddexp

* Update src/basicfuns.jl

Looks good to me as well.

Co-authored-by: David Widmann <[email protected]>

Co-authored-by: David Widmann <[email protected]>
Co-authored-by: David Widmann <[email protected]>

Author: lucasondel

src/basicfuns.jl

@@ -312,9 +312,18 @@ Return `log(exp(x) + exp(y))`, avoiding intermediate overflow/undeflow, and hand
 non-finite values.
 """
 function logaddexp(x::Real, y::Real)
-    # ensure Δ = 0 if x = y = ± Inf
-    Δ = x == y ? zero(x - y) : abs(x - y)
-    max(x, y) + log1pexp(-Δ)
+    # Compute max = Base.max(x, y) and diff = x == y ? zero(x - y) : -abs(x - y)
+    # in a faster type-stable way
+    a, b = promote(x, y)
+    if a < b
+        diff = a - b
+        max = b
+    else
+        # ensure diff = 0 if a = b = ± Inf
+        diff = a == b ? zero(a - b) : b - a
+        max = !isnan(b) ? a : b
+    end
+    return max + log1pexp(diff)
 end
 
 Base.@deprecate logsumexp(x::Real, y::Real) logaddexp(x, y)