Sometimes I am asked questions like "why did you introduce function xyz into MATLAB? And sometimes I have a compelling answer, even for something that looks simple on the surface. Consider the function hypot.
Contents
What Does hypot Compute?
hypot essentially computes the square root of the sum of the squares of 2 inputs. So what's the big deal, right? For "reasonable" input values, there is no big issue.
fhypot = @(a,b) sqrt(abs(a).^2+abs(b).^2);
Other Ways to Compute hypot
Let's create a set of values and compute hypot 3 ways for these: with hypot, fhypot, and sqrt(2)*x. These should all give the same answers, provide a is real and positive.
a = logspace(0,5,6)
format long
aHypot = hypot(a,a)
aFhypot = fhypot(a,a)
aSqrt2 = sqrt(2)*aa =
1 10 100 1000 10000 100000
aHypot =
1.0e+005 *
Columns 1 through 3
0.000014142135624 0.000141421356237 0.001414213562373
Columns 4 through 6
0.014142135623731 0.141421356237310 1.414213562373095
aFhypot =
1.0e+005 *
Columns 1 through 3
0.000014142135624 0.000141421356237 0.001414213562373
Columns 4 through 6
0.014142135623731 0.141421356237310 1.414213562373095
aSqrt2 =
1.0e+005 *
Columns 1 through 3
0.000014142135624 0.000141421356237 0.001414213562373
Columns 4 through 6
0.014142135623731 0.141421356237310 1.414213562373095
Results
The results are the same, to within round-off, for the 3 methods here. But what happens if the magnitude of a is larger, near realmax perhaps?
realmax
ans =
1.797693134862316e+308
For my computer, a 32-bit Windows machine, realmax for doubles is on the order of 10^308. Let's see what happens if a value nearly that large is used with the different versions of hypot.
a = 1e308 aHypot = hypot(a,a) aFhypot = fhypot(a,a) aSqrt2 = sqrt(2)*a
a =
1.000000000000000e+308
aHypot =
1.414213562373095e+308
aFhypot =
Inf
aSqrt2 =
1.414213562373095e+308
What you can see is that the straight-forward method returns Inf instead of a finite answer. And that's why hypot was added to MATLAB, to compute the hypotenuse robustly, avoiding both underflow and overflow.
Get
the MATLAB code
Published with MATLAB® 7.5
That was a pretty interesting post Loren. I’m not sure if it will affect the way I do things around here, but it’s something that could come in useful in the future. Thanks for the tip!
Before hypot existed, I used a=abs(complex(a,b)); In doing so, I was able to tap the robustness of the abs() function.
Duane-
Thanks for pointing out that abs is also implemented robustly.
–Loren
Loren,
Does hypot(a,a) any different from norm([a a])?
Hal
Hal-
norm is also careful to scale results, like abs and hypot.
–Loren