{"id":6633,"date":"2020-01-06T05:46:15","date_gmt":"2020-01-06T04:46:15","guid":{"rendered":"http:\/\/www.walkingrandomly.com\/?p=6633"},"modified":"2020-01-06T06:11:54","modified_gmt":"2020-01-06T05:11:54","slug":"hypot-a-story-of-a-simple-function","status":"publish","type":"post","link":"https:\/\/walkingrandomly.com\/?p=6633","title":{"rendered":"Hypot – A story of a ‘simple’ function"},"content":{"rendered":"

My stepchildren are pretty good at mathematics for their age and have recently learned about Pythagora’s theorem<\/a><\/p>\n

$c=\\sqrt{a^2+b^2}$<\/p>\n

The fact that they have learned about this so early in their mathematical lives is testament to its importance. Pythagoras is everywhere in computational science and it may well be the case that you’ll need to compute the hypotenuse to a triangle some day.<\/p>\n

Fortunately for you, this important computation is implemented in every computational environment I can think of!
\nIt’s almost always called hypot<\/a> so it will be easy to find.
\nHere it is in action using Python’s numpy module<\/p>\n

import numpy as np\r\na = 3\r\nb = 4\r\nnp.hypot(3,4)\r\n\r\n5\r\n<\/pre>\n
When I’m out and about giving talks and tutorials about Research Software Engineering<\/a>, High Performance Computing<\/a> and so on<\/a>, I often get the chance to mention the hypot function and it turns out that fewer people know about this routine than you might expect.<\/p>\n
Trivial Calculation? Do it Yourself!<\/h3>\n
Such a trivial calculation, so easy to code up yourself! Here’s a one-line implementation<\/p>\n
def mike_hypot(a,b):\r\n    return(np.sqrt(a*a+b*b))\r\n<\/pre>\n
In use it looks fine<\/p>\n
mike_hypot(3,4)\r\n\r\n5.0\r\n<\/pre>\n
Overflow and Underflow<\/h3>\n
I could probably work for quite some time before I found that my implementation was flawed in several places. Here’s one<\/p>\n
mike_hypot(1e154,1e154)\r\n\r\ninf\r\n<\/pre>\n
You would, of course, expect the result to be large but not infinity. Numpy doesn’t have this problem<\/p>\n
np.hypot(1e154,1e154)\r\n\r\n1.414213562373095e+154\r\n<\/pre>\n
My function also doesn’t do well when things are small.<\/p>\n
a = mike_hypot(1e-200,1e-200)\r\n\r\n0.0\r\n<\/pre>\n
but again, the more carefully implemented hypot function in numpy does fine.<\/p>\n
np.hypot(1e-200,1e-200)\r\n\r\n1.414213562373095e-200\r\n<\/pre>\n
Standards Compliance<\/h3>\n
Next up — standards compliance. It turns out that there is a an official standard for how hypot implementations should behave in certain edge cases. The IEEE-754 standard for floating point arithmetic has something to say about how any implementation of hypot handles NaNs (Not a Number) and inf (Infinity).<\/p>\n
It states that any implementation of hypot should behave as follows (Here’s a human readable summary https:\/\/www.agner.org\/optimize\/nan_propagation.pdf)<\/a><\/p>\n
hypot(nan,inf) = hypot(inf,nan) = inf\r\n<\/pre>\n
numpy behaves well!<\/p>\n
np.hypot(np.nan,np.inf)\r\n\r\ninf\r\n\r\nnp.hypot(np.inf,np.nan)\r\n\r\ninf\r\n<\/pre>\n
My implementation does not<\/p>\n
mike_hypot(np.inf,np.nan)\r\n\r\nnan\r\n<\/pre>\n
So in summary, my implementation is<\/p>\n