diff --git a/math/dw_math.c b/math/dw_math.c index 50b35344c11f078495b6235f6a8ed5e5690d3eca..288a99eeb3c930d052d78975f0c61eea5642eda6 100644 --- a/math/dw_math.c +++ b/math/dw_math.c @@ -34,14 +34,19 @@ PRECISION AddLogs(PRECISION a, PRECISION b) /* Returns ln(x*exp(a) + y*exp(b)) computed to avoid overflow. If a = ln(c) and - b = ln(d), as is usually the case, then the routine returns ln(x*c + y*d). - Assumes that x*exp(a) + y*exp(b) is positive, but no checking is done. This - condition will always be satisfied if both x and y are positive. - + b = ln(d), as is usually the case, then the routine returns ln(x*c + y*d). + Both x and y must be non-negative. If either x or y is negative, it is + treated as if it were zero. */ PRECISION AddScaledLogs(PRECISION x, PRECISION a, PRECISION y, PRECISION b) { - return (a > b) ? a + log(x + y*exp(b-a)) : b + log(x*exp(a-b) + y); + if (x > 0) + if (y > 0) + return (a > b) ? a + log(x + y*exp(b-a)) : b + log(x*exp(a-b) + y); + else + return log(x) + a; + else + return (y > 0) ? log(y) + b : MINUS_INFINITY; } /*