Updated 2013-09-19 14:07:35 by dkf

Purpose: Discuss the math function hypot.

The call,
`    [expr { hypot( \$x, \$y ) }]`

returns the length of the hypotenuse of a right triangle when the length of the legs are \$x and \$y. It is equivalent to
`    [expr { sqrt( \$x * \$x + \$y * \$y ) }]`

except that it can handle cases where \$x or \$y are extremely large or extremely small without getting overflow or underflow.

In cases where there are no problems with being very close to the limits of the underlying double representation and you don't actually need the distance itself (e.g. in circle intersection testing and checking to see if one point is close enough to another) it is often quicker to work with the square of the length of the hypotenuse and avoid the expensive square root operation. I do not know how this ranks up with the additional overheads of math in Tcl though. DKF

The most common use of hypot is in Converting between rectangular and polar co-ordinates.

hypot is also available in Tclx.

GPS The hypotenuse is very useful for drawing lines pixel by pixel. For example:
``` #Copyright 2003 George Peter Staplin
#You may use this under the same terms as Tcl.
proc draw.line.on.image {img x1 y1 x2 y2 color} {
set xDiff [expr {\$x2 - \$x1}]
set yDiff [expr {\$y2 - \$y1}]

set numPixels [expr {hypot(\$xDiff,\$yDiff)}]
set xRatio    [expr {\$xDiff / \$numPixels}]
set yRatio    [expr {\$yDiff / \$numPixels}]

for {set p 0} {\$p < \$numPixels} {incr p} {
set x [expr {round(\$xRatio * \$p) + \$x1}]
set y [expr {round(\$yRatio * \$p) + \$y1}]
\$img put \$color -to \$x \$y [expr {\$x + 1}] [expr {\$y + 1}]
}
}

proc main {} {
set img [image create photo -width 300 -height 300]
draw.line.on.image \$img 10 10 100 100 green
draw.line.on.image \$img 50 20  50 200 blue
draw.line.on.image \$img 40 50 300  50 maroon

pack [label .l -image \$img]
}
main```

Basically the display is like a right triangle. The left edge of the window is the side opposite the hypotenuse, while the bottom edge is the base side. Therefore using the Pythagorean Theorem which says that the length of the hypotenuse is the square root of: (x * x) + (y * y) we can find the number of pixels between the two points. If you have another use for the hypot function (perhaps for 3D) then I (GPS) would like to see it.

AMG: I most often use it to measure distance in Cartesian coordinate systems.

Regarding your example, it's actually not correct to use hypot() to determine the number of pixels to draw. Yes, it will give you the distance in pixels between the two points, but (except for horizontal and vertical lines) that will exceed the number of pixels that need to be drawn. (Also I assume that both endpoint pixels should be drawn, so one extra pixel needs to be rendered.) Depending on how the true line intersects a pixel used to approximate it, that single pixel can appear to be up to 1.4142 pixels in length. Assuming that the center of the pixel corresponds to integer coordinates, a line from (0,0) to (5,0) is rendered using six pixels and is five pixels in length. A line from (0,0) to (5,5) is also rendered using six pixels but is a little over seven pixels in length. In the former case, each of the endpoint pixels is effectively half a pixel in length, and all four interior pixels are each one pixel long, making for a total length of five pixels. In the latter case, each endpoint pixel is 0.7071 pixels in length, and the four interior pixels are each 1.4142 pixels long, totaling 7.0711.

Here's code:
``` proc line {img x1 y1 x2 y2 color} {
set num_pixels [expr {max(abs(\$x2 - \$x1), abs(\$y2 - \$y1))}]
if {\$num_pixels == 0} {
set dx 0
set dy 0
} else {
set dx [expr {double(\$x2 - \$x1) / \$num_pixels}]
set dy [expr {double(\$y2 - \$y1) / \$num_pixels}]
}
for {set i 0} {\$i <= \$num_pixels} {incr i} {
set x [expr {int(\$x1 + \$i * \$dx)}]
set y [expr {int(\$y1 + \$i * \$dy)}]
\$img put \$color -to \$x \$y [expr {\$x + 1}] [expr {\$y + 1}]
}
}```

The idea is that one (or both, in a 45 degree angle line) of \$x and \$y will advance by one (or negative one) every iteration. Except in the case of a 45 degree angle line, the other variable will advance by less than one and more than negative one. In the magnificently elegant language of mathematics, max(abs(\$dx), abs(\$dy)) = 1. Oh, and \$num_pixels = 0 is a degenerate case in which I arbitrarily choose to draw one pixel.

Of course, there exist far better methods of drawing lines. In particular, proper line drawing algorithms are capable of rendering runs of pixels, which are rendered efficiently using single-pixel-wide rectangles. Also it's good to balance the line so that the error is equally distributed on either side of the true line. If anyone is interested, I will post some code. But this is getting off-topic, since I am not aware of any quality line drawing algorithms that use hypot() or even sqrt().

AMG: Feature request! Why not extend hypot() to take any number of arguments? In general, hypot() is the square root of the sum of the squares of the arguments. Right now it only supports two arguments, making it painful to use in 3D and above. Single-argument hypot() should work the same as abs(), and zero-argument hypot() should return 0.

Lars H: Reasonable idea, although I'm not sure I'd count
`  expr {hypot(hypot(\$x,\$y),\$z)}`

as particularly painful.

AMG: True, I had that in mind, but I still find it painful to think that this code is harder to write, harder to read, and takes longer to execute, in comparison to the way it "should" be. It's gross to contemplate the fact that internally numbers are being square root'ed and then immediately squared again.

Lars H, insertion: If that is how it is really implemented. I vaguely recall that hypot (at least in older times) was rather implemented as
``` proc hypot {x y} {
set x [expr {abs(\$x)}]
set y [expr {abs(\$y)}]
if {\$x>=\$y} then {
set a \$x
set b [expr {\$y/\$x}]
set t [expr {\$b*\$b}]
} else {
set a \$y
set b [expr {\$x/\$y}]
set t [expr {\$b*\$b}]
}
return [expr {\$a * f(\$t)}]
}```

where the helper function f(t) = sqrt(1+t) is implemented by a MacLaurin expansion as f(t) = 1 + t/2 – 3*t**2/8 + 5*t**3/16 – ... (Whether it is worth the trouble probably depends on whether sqrt is available in hardware.)

AMG, continued: Also it turns out that this is more efficient:
` expr {sqrt(\$x * \$x + \$y * \$y + \$z * \$z)}`

According to [time], the sqrt method takes 2.63027 microseconds on my machine, and the double hypot takes 3.64441 microseconds, which is 38.6% longer.

DKF: The problem is that the C library hypot() function only takes two arguments, so to take more would require reimplementation from scratch (and hence we'd need to worry about numeric stability, etc.)

Lars H, 2008-04-23: A look in the source of MetaFont turns up the following algorithm for hypot (attributed to Cleve Moler and Donald Morrison, IBM Journal of Research and Development 27 (1983), 577–581):
``` proc hypot {a b} {
set a [expr {double(abs(\$a))}]
set b [expr {double(abs(\$b))}]
if {\$a<\$b} then {set r \$a; set a \$b; set b \$r}
while 1 {
set r [expr {1/(1+4*\$a*\$a/(\$b*\$b))}]
if {\$r<1e-20} then {break}
set a [expr {\$a + 2*\$r*\$a}]
set b [expr {\$b*\$r}]
}
return \$a
}```

The idea in the loop is to reflect the vector (a,b) in the vector (a,b/2) until b can be neglected; reflection preserves the length of the vector, so what one ends up with will be (hypot(a,b),0). Since the angle between (a,0) and (a,b/2) is slightly larger than that between (a,b/2) and (a,b), one doesn't get all the way to the a-axis in one step, but one gets very close — typically three rounds through the loop is sufficient.
``` % list [hypot 1 1] [tcl::mathfunc::hypot 1 1]
1.4142135623730951 1.4142135623730951
% list [hypot 3 4] [tcl::mathfunc::hypot 3 4]
5.000000000000001 5.0```

The slight overshoot is due to accumulated rounding errors.