微信公众号：Julia语言
每周一三五更新Julia语言；
每周二四六更新Python进阶；
What’s so hard about finding a hypotenuse?
https://www.johndcook.com/blog/2010/06/02/whats-so-hard-about-finding-a-hypotenuse/
Numerical libraries usually have a function for finding the hypotenuse of a right triangle. Isn’t this trivial? If the sides of a triangle are x and y, just write this:
sqrt(x*x + y*y)
That works in theory, but in practice it may fail. If x is so large that x*x overflows, the code will produce an infinite result.
We’d like to be able to say to the computer “Now xx and yy might be too big, but just wait a minute. I’m going to take a square root, and so the numbers will get smaller. I just need to borrow some extra range for a minute.” You can do something like that in environments that have extended precision, but that would be inefficient and unnecessary. And of course it would fail completely if you’re already using the largest numeric type available.
Here’s how to compute sqrt(xx + yy) without risking overflow.
- max = maximum(|x|, |y|)
- min = minimum(|x|, |y|)
- r = min / max
- return max*sqrt(1 + r^2)
Notice that the square root argument is between 1 and 2 and so there is no danger of overflow in taking the square root. The only way the algorithm could overflow is if the result itself is too large to represent. In that case the fault lies with the problem, not the algorithm: the algorithm was asked to compute a quantity larger than the largest representable number, and so it does the correct thing by overflowing and returning an infinite value. However, the algorithm will not unnecessarily return an infinite value due to an overflow in an intermediate computation.
To see how the algorithm above may succeed when the naive implementation would fail, let M be the largest representable number. For IEEE double precision, M is on the order of 10^308. All that matters for this example is that M is greater than 4. Let x and y equal 0.5 M. The algorithm above would return 0.707 M. But the naive algorithm would compute xx and overflow since 0.25 M² > M. Then xx + y*y would evaluate to infinity and the square root would evaluate to infinity.
备注
Julia语言中，可以直接使用hypot(x, y)函数求取斜边。
julia> hypot(3,4)
5.0
julia> hypot(5,12)
13.0
欢迎关注微信公众账号Julia语言.jpg
点击阅读原文可查看历史文章
网友评论