{"id":221,"date":"2012-07-30T12:35:17","date_gmt":"2012-07-30T17:35:17","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=221"},"modified":"2013-05-02T10:04:38","modified_gmt":"2013-05-02T15:04:38","slug":"pythagorean-addition","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2012\/07\/30\/pythagorean-addition\/","title":{"rendered":"Pythagorean Addition"},"content":{"rendered":"\r\n<!DOCTYPE html\r\n  PUBLIC \"-\/\/W3C\/\/DTD HTML 4.01 Transitional\/\/EN\">\r\n<style type=\"text\/css\">\r\n\r\nh1 { font-size:18pt; }\r\nh2.titlebg { font-size:13pt; }\r\nh3 { color:#4A4F55; padding:0px; margin:5px 0px 5px; font-family:Arial, Helvetica, sans-serif; font-size:11pt; font-weight:bold; line-height:140%; border-bottom:1px solid #d6d4d4; display:block; }\r\nh4 { color:#4A4F55; padding:0px; margin:0px 0px 5px; font-family:Arial, Helvetica, sans-serif; font-size:10pt; font-weight:bold; line-height:140%; border-bottom:1px solid #d6d4d4; display:block; }\r\n   \r\np { padding:0px; margin:0px 0px 20px; }\r\nimg { padding:0px; margin:0px 0px 20px; border:none; }\r\np img, pre img, tt img, li img { margin-bottom:0px; } \r\n\r\nul { padding:0px; margin:0px 0px 20px 23px; list-style:square; }\r\nul li { padding:0px; margin:0px 0px 7px 0px; background:none; }\r\nul li ul { padding:5px 0px 0px; margin:0px 0px 7px 23px; }\r\nul li ol li { list-style:decimal; }\r\nol { padding:0px; margin:0px 0px 20px 0px; list-style:decimal; }\r\nol li { padding:0px; margin:0px 0px 7px 23px; list-style-type:decimal; }\r\nol li ol { padding:5px 0px 0px; margin:0px 0px 7px 0px; }\r\nol li ol li { list-style-type:lower-alpha; }\r\nol li ul { padding-top:7px; }\r\nol li ul li { list-style:square; }\r\n\r\npre, tt, code { font-size:12px; }\r\npre { margin:0px 0px 20px; }\r\npre.error { color:red; }\r\npre.codeinput { padding:10px; border:1px solid #d3d3d3; background:#f7f7f7; }\r\npre.codeoutput { padding:10px 11px; margin:0px 0px 20px; color:#4c4c4c; }\r\n\r\n@media print { pre.codeinput, pre.codeoutput { word-wrap:break-word; width:100%; } }\r\n\r\nspan.keyword { color:#0000FF }\r\nspan.comment { color:#228B22 }\r\nspan.string { color:#A020F0 }\r\nspan.untermstring { color:#B20000 }\r\nspan.syscmd { color:#B28C00 }\r\n\r\n.footer { width:auto; padding:10px 0px; margin:25px 0px 0px; border-top:1px dotted #878787; font-size:0.8em; line-height:140%; font-style:italic; color:#878787; text-align:left; float:none; }\r\n.footer p { margin:0px; }\r\n\r\n  <\/style><div class=\"content\"><!--introduction--><p>How do you compute the hypotenuse of a right triangle without squaring the lengths of the sides and without taking any square roots?<\/p><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#3ece0897-c012-436a-b1dd-874b4bd34881\">Some Important Operations<\/a><\/li><li><a href=\"#a4aa5d92-820f-48c6-94cb-514a95b8522d\">Pythagorean Addition<\/a><\/li><li><a href=\"#8ba99af0-e4fb-4393-bfa7-afd1cf305553\">Underflow and Overflow<\/a><\/li><li><a href=\"#3e33acef-8c94-4fac-b2c2-9d6aa1dcb95c\">Don Morrison<\/a><\/li><li><a href=\"#2ae711ee-d40c-4fbf-bd8e-ee8cb8aafe59\">Don's Diagram<\/a><\/li><li><a href=\"#8bef8ced-dd4f-4995-96f1-5b93b425ac4e\">Function Pythag<\/a><\/li><li><a href=\"#fb276887-9027-42d8-8029-38a01f058f33\">Examples<\/a><\/li><li><a href=\"#4aebde3c-d445-4317-a66c-0f9aec6400d2\">Augustin Dubrulle<\/a><\/li><li><a href=\"#40e93728-8e20-44f2-966a-834ca7ede848\">Pythag Today?<\/a><\/li><\/ul><\/div><h4>Some Important Operations<a name=\"3ece0897-c012-436a-b1dd-874b4bd34881\"><\/a><\/h4><p>These are all important operations.<\/p><div><ul><li>Compute the 2-norm of a vector, $||v||_2$.<\/li><li>Find complex magnitude, $|x + iy|$.<\/li><li>Convert from cartestesian to polar coordinates, $x + iy = r e^{i \\theta}$<\/li><li>Compute an arctangent, $\\theta = \\arctan{y\/x}$<\/li><li>Find a plane rotation that zeros one component of a two-vector.<\/li><li>Find an orthogonal reflection that zeros $n-1$ components of an $n$-vector.<\/li><\/ul><\/div><p>All of them involve computing<\/p><p>$$ \\sqrt{x^2 + y^2} $$<\/p><p>in some way or another.<\/p><h4>Pythagorean Addition<a name=\"a4aa5d92-820f-48c6-94cb-514a95b8522d\"><\/a><\/h4><p>Let's introduce the notation $\\oplus$ for what we call <i>Pythagorean addition<\/i>.<\/p><p>$$ x \\oplus y = \\sqrt{x^2 + y^2} $$<\/p><p>This has some of the properties of ordinary addition, at least on the nonnegative real numbers.<\/p><p>You can use Pythagorean addition repeatedly to compute the 2-norm of a vector $v$ with components $v_1, v_2, \\ldots, v_n$.<\/p><p>$$||v||_2 = v_1 \\oplus v_2 \\oplus \\ldots \\oplus v_n$$<\/p><p>It is easy to see how Pythagorean addition is involved in the other operations listed above.<\/p><h4>Underflow and Overflow<a name=\"8ba99af0-e4fb-4393-bfa7-afd1cf305553\"><\/a><\/h4><p>Computationally, it is essential to avoid unnecessary overflow and underflow of floating point numbers. IEEE double precision has the following range. Any values outside this range are too small or too large to be represented.<\/p><pre class=\"codeinput\">format <span class=\"string\">compact<\/span>\r\nformat <span class=\"string\">short<\/span> <span class=\"string\">e<\/span>\r\nrange = [eps*realmin realmax]\r\n<\/pre><pre class=\"codeoutput\">range =\r\n  4.9407e-324  1.7977e+308\r\n<\/pre><p>This crude attempt to implement Pythagorean addition is not satisfactory because the intermediate results underflow or overflow.<\/p><pre class=\"codeinput\">bad_pythag = @(x,y) sqrt(x^2 + y^2)\r\n<\/pre><pre class=\"codeoutput\">bad_pythag = \r\n    @(x,y)sqrt(x^2+y^2)\r\n<\/pre><p>If <tt>x<\/tt> and <tt>y<\/tt> are so small that their squares underflow, then <tt>bad_pythag(x,y)<\/tt> will be zero even though the true result can be represented.<\/p><pre class=\"codeinput\">x = 3e-200\r\ny = 4e-200\r\nz = 5e-200 <span class=\"comment\">% should be the result<\/span>\r\nz = bad_pythag(x,y)\r\n<\/pre><pre class=\"codeoutput\">x =\r\n  3.0000e-200\r\ny =\r\n  4.0000e-200\r\nz =\r\n  5.0000e-200\r\nz =\r\n     0\r\n<\/pre><p>If <tt>x<\/tt> and <tt>y<\/tt> are so large that their squares overflow, then <tt>bad_pythag(x,y)<\/tt> will be infinity even though the true result can be represented.<\/p><pre class=\"codeinput\">x = 3e200\r\ny = 4e200\r\nz = 5e200 <span class=\"comment\">% should be the result<\/span>\r\nz = bad_pythag(x,y)\r\n<\/pre><pre class=\"codeoutput\">x =\r\n  3.0000e+200\r\ny =\r\n  4.0000e+200\r\nz =\r\n  5.0000e+200\r\nz =\r\n   Inf\r\n<\/pre><h4>Don Morrison<a name=\"3e33acef-8c94-4fac-b2c2-9d6aa1dcb95c\"><\/a><\/h4><p>Don Morrison was a mathematician and computer pioneer who spent most of his career at Sandia National Laboratory in Albuquerque. He left Sandia in the late '60s, founded the Computer Science Department at the University of New Mexico, and recruited me to join the university a few years later.<\/p><p>Don had all kinds of fascinating mathematical interests. He was an expert on cryptography. He developed fair voting systems for multi-candidate elections. He designed an on-demand public transportation system for the city of Albuquerque. He constructed a kind of harmonica that played computer punched cards. He discovered the Fast Fourier Transform before Culley and Tuckey, and published the algorithm in the proceedings of a regional ACM conference.<\/p><p>This is the first of two or three blogs that I intend to write about things I learned from Don.<\/p><h4>Don's Diagram<a name=\"2ae711ee-d40c-4fbf-bd8e-ee8cb8aafe59\"><\/a><\/h4><p>Don was sitting in on a class I was teaching on mathematical software. I was talking about the importance of avoiding underflow and overflow while computing the 2-norm of a vector. (It was particularly important back then because the IBM mainframes of the day had especially limited floating point exponent range.) We tried to do the computation with just one pass over the data, to avoid repeated access to main memory. This involved messy dynamic rescaling. It was also relatively expensive to compute square roots. Before the end of the class Don had sketched something like the following.<\/p><pre class=\"codeinput\">pythag_pic(4,3)\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/pythag_blog_01.png\" alt=\"\"> <p>We are at the point $(x,y)$, with $|y| \\le |x|$. We want to find the radius of the black circle without squaring $x$ or $y$ and without computing any square roots. The green line leads from point $(x,y)$ to its projection $(x,0)$ on the $x$-axis. The blue line extends from the origin through the midpoint of the green line. The red line is perpendicular to the blue line. The red line intersects the circle in the point $(x+,y+)$. Don realized that $x+$ and $y+$ could be computed from $x$ and $y$ with a few safe rational operations, and that $y+$ would be much smaller than $y$, so that $x+$ would be a much better approximation to the radius than $x$. The process could then be repeated a few times to get an excellent approximation to the desired result.<\/p><h4>Function Pythag<a name=\"8bef8ced-dd4f-4995-96f1-5b93b425ac4e\"><\/a><\/h4><p>Here, in today's MATLAB, is the algorithm. It turns out that the iteration is cubically convergent, so at most three iterations produces double precision accuracy.  It is not worth the trouble to check for convergence in fewer than three iterations.<\/p><pre class=\"codeinput\">type <span class=\"string\">pythag<\/span>\r\n<\/pre><pre class=\"codeoutput\">\r\nfunction x = pythag(a,b)\r\n% PYTHAG  Pythagorean addition\r\n% pythag(a,b) = sqrt(a^2+b^2) without unnecessary\r\n% underflow or overflow and without any square roots.\r\n   if a==0 &amp;&amp; b==0\r\n      x = 0;\r\n   else\r\n      % Start with abs(x) &gt;= abs(y)\r\n      x = max(abs(a),abs(b));\r\n      y = min(abs(a),abs(b));\r\n      % Iterate three times\r\n      for k = 1:3\r\n         r = (y\/x)^2;\r\n         s = r\/(4+r);\r\n         x = x + 2*s*x;\r\n         y = s*y;\r\n      end\r\n   end\r\nend\r\n\r\n<\/pre><p>Computing <tt>r = (y\/x)^2<\/tt> is safe because the square will not overflow and, if it underflows, it is negligible. There are only half a dozen other floating point operations per iteration and they are all safe.<\/p><p>It is not obvious, but the quantity $x \\oplus y$ is a loop invariant.<\/p><p>Surprisingly, this algorithm cannot be used to compute square roots.<\/p><h4>Examples<a name=\"fb276887-9027-42d8-8029-38a01f058f33\"><\/a><\/h4><p>Starting with $x = y$ is the slowest to converge.<\/p><pre class=\"codeinput\">format <span class=\"string\">long<\/span> <span class=\"string\">e<\/span>\r\nformat <span class=\"string\">compact<\/span>\r\npythag_with_disp(1,1)\r\n\r\nsqrt(2)\r\n<\/pre><pre class=\"codeoutput\">     1     1\r\n     1.400000000000000e+00     2.000000000000000e-01\r\n     1.414213197969543e+00     1.015228426395939e-03\r\n     1.414213562373095e+00     1.307981162604408e-10\r\nans =\r\n     1.414213562373095e+00\r\nans =\r\n     1.414213562373095e+00\r\n<\/pre><p>It's fun to compute Pythagorean triples, which are pairs of integers whose Pythagorean sum is another integer.<\/p><pre class=\"codeinput\">pythag_with_disp(4e-300,3e-300)\r\n<\/pre><pre class=\"codeoutput\">    4.000000000000000e-300    3.000000000000000e-300\r\n    4.986301369863013e-300    3.698630136986302e-301\r\n    4.999999974188252e-300    5.080526329415360e-304\r\n    5.000000000000000e-300    1.311372652398298e-312\r\nans =\r\n    5.000000000000000e-300\r\n<\/pre><pre class=\"codeinput\">pythag_with_disp(12e300,5e300)\r\n<\/pre><pre class=\"codeoutput\">    1.200000000000000e+301    5.000000000000000e+300\r\n    1.299833610648919e+301    2.079866888519135e+299\r\n    1.299999999999319e+301    1.331199999999652e+295\r\n    1.300000000000000e+301    3.489660928000008e+282\r\nans =\r\n    1.300000000000000e+301\r\n<\/pre><h4>Augustin Dubrulle<a name=\"4aebde3c-d445-4317-a66c-0f9aec6400d2\"><\/a><\/h4><p>Augustin Dubrulle is a French-born numerical analyst who was working for IBM in Houston in the 1970s on SSP, their Scientific Subroutine Package. He is the only person I know of who ever improved on an algorithm of J. H. Wilkinson.  Wilkinson and Christian Reinsch had published, in <i>Numerische Mathematik<\/i>, two versions of the symmetric tridiagonal QR algorithm for matrix eigenvalues.  The explicit shift version required fewer operations, but the implicit shift version had better numerical properties.  Dubrulle published a half-page paper in <i>Numerische Mathematik<\/i> that said, in effect, \"make the following change to the inner loop of the algorithm of Wilkinson and Reinsch\" to combine the superior properties of both versions.  This is the algorithm we use today in MATLAB to compute eigenvalues of symmetric matrices.<\/p><p>Augustin came to New Mexico to get his Ph. D. and became interested in the <i>pythag<\/i> algorithm.  He analyzed its convergence properties, showed its connection to Halley's method for computing square roots, and produced higher order generalizations. The three of us published two papers back to back, <a href=\"https:\/\/blogs.mathworks.com\/images\/cleve\/moler_morrison.pdf\">Moler and Morrison<\/a>, and <a href=\"https:\/\/blogs.mathworks.com\/images\/cleve\/dubrulle.pdf\">Dubrulle<\/a>, in the <i>IBM Journal of Research and Development<\/i> in 1983.<\/p><h4>Pythag Today?<a name=\"40e93728-8e20-44f2-966a-834ca7ede848\"><\/a><\/h4><p>What has become of <i>pythag<\/i>? Its functionality lives on in <tt>hypot<\/tt>, which is part of <tt>libm<\/tt>, the fundamental math library support for C. There is a <tt>hypot<\/tt> function in MATLAB.<\/p><p>On Intel chips, and on chips that use Intel libraries, we rely upon the Intel Math Kernel Library to compute <tt>hypot<\/tt>. Modern Intel architectures have two features that we did not have in the old days.  First, square root is an acceptably fast machine instruction, so this implementation would be OK.<\/p><pre class=\"codeinput\">type <span class=\"string\">ok_hypot<\/span>\r\n<\/pre><pre class=\"codeoutput\">\r\nfunction r = ok_hypot(x,y)\r\n   if x==0 &amp;&amp; y==0\r\n      r = 0;\r\n   elseif abs(x) &gt;= abs(y)\r\n      r = abs(x)*sqrt(1+abs(y\/x)^2);\r\n   else\r\n      r = abs(y)*sqrt(1+abs(x\/y)^2);\r\n   end\r\nend\r\n\r\n<\/pre><p>But even that kind of precaution isn't necessary today because of the other relevant feature of the Intel architecture, the extended floating point registers.  These registers are accessible only to library developers working in machine language. They provide increased precision and, more important in this situation, increased exponent range. So, if you start with standard double precision numbers and do the entire computation in the extended registers, you can get away with <tt>bad_pythag<\/tt>. In this case, clever hardware obviates the need for clever software.<\/p><script language=\"JavaScript\"> <!-- \r\n    function grabCode_1c54f057d0e44730a57f2ce2519b3afb() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='1c54f057d0e44730a57f2ce2519b3afb ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 1c54f057d0e44730a57f2ce2519b3afb';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        copyright = 'Copyright 2012 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add copyright line at the bottom if specified.\r\n        if (copyright.length > 0) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\\n');\r\n\r\n        d.title = title + ' (MATLAB code)';\r\n        d.close();\r\n    }   \r\n     --> <\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_1c54f057d0e44730a57f2ce2519b3afb()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n      the MATLAB code <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; 7.14<br><\/p><p class=\"footer\"><br>\r\n      Published with MATLAB&reg; 7.14<br><\/p><\/div><!--\r\n1c54f057d0e44730a57f2ce2519b3afb ##### SOURCE BEGIN #####\r\n%% Pythagorean Addition \r\n% How do you compute the hypotenuse of a right triangle without\r\n% squaring the lengths of the sides and without taking any square roots?\r\n\r\n%% Some Important Operations\r\n% These are all important operations.\r\n%\r\n% * Compute the 2-norm of a vector, $||v||_2$.\r\n% * Find complex magnitude, $|x + iy|$.\r\n% * Convert from cartestesian to polar coordinates, $x + iy = r e^{i \\theta}$\r\n% * Compute an arctangent, $\\theta = \\arctan{y\/x}$\r\n% * Find a plane rotation that zeros one component of a two-vector.\r\n% * Find an orthogonal reflection that zeros $n-1$ components of an $n$-vector.\r\n%\r\n% All of them involve computing\r\n%\r\n% $$ \\sqrt{x^2 + y^2} $$\r\n%\r\n% in some way or another.\r\n\r\n%% Pythagorean Addition\r\n% Let's introduce the notation $\\oplus$ for what we call _Pythagorean addition_.\r\n%\r\n% $$ x \\oplus y = \\sqrt{x^2 + y^2} $$\r\n%\r\n% This has some of the properties of ordinary addition, at least on the\r\n% nonnegative real numbers.\r\n\r\n%%\r\n% You can use Pythagorean addition repeatedly to compute the 2-norm\r\n% of a vector $v$ with components $v_1, v_2, \\ldots, v_n$.\r\n% \r\n% $$||v||_2 = v_1 \\oplus v_2 \\oplus \\ldots \\oplus v_n$$\r\n%\r\n% It is easy to see how Pythagorean addition is involved in the\r\n% other operations listed above.\r\n\r\n%% Underflow and Overflow\r\n% Computationally, it is essential to avoid unnecessary overflow and underflow\r\n% of floating point numbers. IEEE double precision has the following range.\r\n% Any values outside this range are too small or too large to be represented.\r\n\r\nformat compact\r\nformat short e\r\nrange = [eps*realmin realmax]\r\n\r\n%%\r\n% This crude attempt to implement Pythagorean addition is not satisfactory\r\n% because the intermediate results underflow or overflow.\r\nbad_pythag = @(x,y) sqrt(x^2 + y^2)\r\n\r\n%%\r\n% If |x| and |y| are so small that their squares underflow, then \r\n% |bad_pythag(x,y)| will be zero even though the true result can be represented.\r\n\r\nx = 3e-200\r\ny = 4e-200\r\nz = 5e-200 % should be the result\r\nz = bad_pythag(x,y)\r\n\r\n%%\r\n% If |x| and |y| are so large that their squares overflow, then \r\n% |bad_pythag(x,y)| will be infinity even though the true result can be\r\n% represented.\r\n\r\nx = 3e200\r\ny = 4e200\r\nz = 5e200 % should be the result\r\nz = bad_pythag(x,y)\r\n\r\n%% Don Morrison\r\n% Don Morrison was a mathematician and computer pioneer who spent most of his\r\n% career at Sandia National Laboratory in Albuquerque.\r\n% He left Sandia in the late '60s, founded the Computer Science Department\r\n% at the University of New Mexico, and recruited me to join the university\r\n% a few years later.\r\n\r\n%%\r\n% Don had all kinds of fascinating mathematical interests.\r\n% He was an expert on cryptography.\r\n% He developed fair voting systems for multi-candidate elections.\r\n% He designed an on-demand public transportation system for the city of\r\n% Albuquerque.\r\n% He constructed a kind of harmonica that played computer punched cards.\r\n% He discovered the Fast Fourier Transform before Culley and Tuckey,\r\n% and published the algorithm in the proceedings of a regional ACM conference.\r\n\r\n%%\r\n% This is the first of two or three blogs that I intend to write about\r\n% things I learned from Don.\r\n\r\n%% Don's Diagram\r\n% Don was sitting in on a class I was teaching on mathematical software.\r\n% I was talking about the importance of avoiding underflow and overflow\r\n% while computing the 2-norm of a vector.\r\n% (It was particularly important back then because the IBM mainframes of the\r\n% day had especially limited floating point exponent range.)\r\n% We tried to do the computation with just one pass over the data,\r\n% to avoid repeated access to main memory.\r\n% This involved messy dynamic rescaling.\r\n% It was also relatively expensive to compute square roots.\r\n% Before the end of the class Don had sketched something like the following.\r\n\r\npythag_pic(4,3)\r\n\r\n%%\r\n% We are at the point $(x,y)$, with $|y| \\le |x|$.\r\n% We want to find the radius of the black circle without squaring $x$ or $y$\r\n% and without computing any square roots.\r\n% The green line leads from point $(x,y)$ to its projection $(x,0)$ on the\r\n% $x$-axis.\r\n% The blue line extends from the origin through the midpoint of the green line.\r\n% The red line is perpendicular to the blue line.\r\n% The red line intersects the circle in the point $(x+,y+)$.\r\n% Don realized that $x+$ and $y+$ could be computed from $x$ and $y$ with a\r\n% few safe rational operations, and that $y+$ would be much smaller than $y$,\r\n% so that $x+$ would be a much better approximation to the radius than $x$.\r\n% The process could then be repeated a few times to get an excellent\r\n% approximation to the desired result.\r\n\r\n%% Function Pythag\r\n% Here, in today's MATLAB, is the algorithm.\r\n% It turns out that the iteration is cubically convergent, so at most three\r\n% iterations produces double precision accuracy.  It is not worth the trouble\r\n% to check for convergence in fewer than three iterations.\r\n\r\ntype pythag\r\n\r\n%%\r\n% Computing |r = (y\/x)^2| is safe because the square will not overflow and,\r\n% if it underflows, it is negligible.\r\n% There are only half a dozen other floating point operations per iteration\r\n% and they are all safe. \r\n\r\n%%\r\n% It is not obvious, but the quantity $x \\oplus y$ is a loop invariant.\r\n\r\n%%\r\n% Surprisingly, this algorithm cannot be used to compute square roots.\r\n\r\n%% Examples\r\n% Starting with $x = y$ is the slowest to converge.\r\n\r\nformat long e\r\nformat compact\r\npythag_with_disp(1,1)\r\n\r\nsqrt(2)\r\n\r\n%%\r\n% It's fun to compute Pythagorean triples, which are pairs of integers whose\r\n% Pythagorean sum is another integer.\r\n\r\npythag_with_disp(4e-300,3e-300)\r\n\r\n%%\r\n\r\npythag_with_disp(12e300,5e300)\r\n\r\n%% Augustin Dubrulle\r\n% Augustin Dubrulle is a French-born numerical analyst who was working\r\n% for IBM in Houston in the 1970s on SSP, their Scientific Subroutine Package.\r\n% He is the only person I know of who ever improved on an algorithm\r\n% of J. H. Wilkinson.  Wilkinson and Christian Reinsch had published,\r\n% in _Numerische Mathematik_, two versions of the symmetric tridiagonal QR\r\n% algorithm for matrix eigenvalues.  The explicit shift version required\r\n% fewer operations, but the implicit shift version had better\r\n% numerical properties.  Dubrulle published a half-page paper in\r\n% _Numerische Mathematik_ that said, in effect, \"make the following change\r\n% to the inner loop of the algorithm of Wilkinson and Reinsch\" to combine\r\n% the superior properties of both versions.  This is the algorithm we\r\n% use today in MATLAB to compute eigenvalues of symmetric matrices.\r\n\r\n%%\r\n% Augustin came to New Mexico to get his Ph. D. and became interested in\r\n% the _pythag_ algorithm.  He analyzed its convergence properties,\r\n% showed its connection to Halley's method for computing square roots,\r\n% and produced higher order generalizations.\r\n% The three of us published two papers back to back, \r\n% <https:\/\/blogs.mathworks.com\/images\/cleve\/moler_morrison.pdf\r\n% Moler and Morrison>, and\r\n% <https:\/\/blogs.mathworks.com\/images\/cleve\/dubrulle.pdf Dubrulle>,\r\n% in the _IBM Journal of Research and Development_ in 1983.\r\n\r\n%% Pythag Today?\r\n% What has become of _pythag_?\r\n% Its functionality lives on in |hypot|, which is part of |libm|,\r\n% the fundamental math library support for C.\r\n% There is a |hypot| function in MATLAB.\r\n\r\n%%\r\n% On Intel chips, and on chips that use Intel libraries, we rely\r\n% upon the Intel Math Kernel Library to compute |hypot|.\r\n% Modern Intel architectures have two features that we did not have in\r\n% the old days.  First, square root is an acceptably fast machine instruction,\r\n% so this implementation would be OK.\r\n\r\ntype ok_hypot\r\n\r\n%%\r\n% But even that kind of precaution isn't necessary today because of the other\r\n% relevant feature of the Intel architecture, the extended floating point\r\n% registers.  These registers are accessible only to library developers\r\n% working in machine language.\r\n% They provide increased precision and, more important in this situation,\r\n% increased exponent range.\r\n% So, if you start with standard double precision numbers and do the entire\r\n% computation in the extended registers, you can get away with |bad_pythag|.\r\n% In this case, clever hardware obviates the need for clever software.\r\n\r\n##### SOURCE END ##### 1c54f057d0e44730a57f2ce2519b3afb\r\n-->","protected":false},"excerpt":{"rendered":"<!--introduction--><p>How do you compute the hypotenuse of a right triangle without squaring the lengths of the sides and without taking any square roots?... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2012\/07\/30\/pythagorean-addition\/\">read more >><\/a><\/p>","protected":false},"author":78,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[11],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/221"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/users\/78"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/comments?post=221"}],"version-history":[{"count":13,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/221\/revisions"}],"predecessor-version":[{"id":233,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/221\/revisions\/233"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=221"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=221"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=221"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}