diff --git a/common/log1mexp.m b/common/log1mexp.m index 3892c70..51918e1 100644 --- a/common/log1mexp.m +++ b/common/log1mexp.m @@ -1,7 +1,7 @@ function y = log1mexp(x) -% Accurately compute y = log(1-exp(-x)) +% Accurately compute y = log(1-exp(x)) % reference: Accurately Computing log(1-exp(-|a|)) Martin Machler y = x; -i = x > log(2); -y(i) = log1p(-exp(-x(i))); -y(~i) = log(-expm1(-x(~i))); +i = x < -log(2); +y(i) = log1p(-exp(x(i))); +y(~i) = log(-expm1(x(~i))); diff --git a/common/log1pexp.m b/common/log1pexp.m index 7ad0b9d..10096e5 100644 --- a/common/log1pexp.m +++ b/common/log1pexp.m @@ -1,7 +1,8 @@ function y = log1pexp(x) % Accurately compute y = log(1+exp(x)) -% reference: Accurately Computing log(1-exp(|a|)) Martin Machler -seed = 33.3; +% reference: Accurately Computing log(1-exp(-|a|)) Martin Machler y = x; -idx = x 18; +j = i & (x <= 33.3); +y(~i) = log1p(exp(x(~i))); +y(j) = x(j)+exp(-x(j));