{"id":1868,"date":"2016-07-11T12:00:42","date_gmt":"2016-07-11T17:00:42","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=1868"},"modified":"2016-07-11T00:11:41","modified_gmt":"2016-07-11T05:11:41","slug":"the-graeffe-root-squaring-method-for-computing-the-zeros-of-a-polynomial","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2016\/07\/11\/the-graeffe-root-squaring-method-for-computing-the-zeros-of-a-polynomial\/","title":{"rendered":"The Graeffe Root-Squaring Method for Computing the Zeros of a Polynomial"},"content":{"rendered":"\r\n\r\n<div class=\"content\"><!--introduction--><p>At a <a href=\"http:\/\/meetings.siam.org\/sess\/dsp_programsess.cfm?SESSIONCODE=23492\">minisymposium<\/a> honoring Charlie Van Loan this week during the <a href=\"http:\/\/www.siam.org\/meetings\/an16\/\">SIAM Annual Meeting<\/a>, I will describe several dubious methods for computing the zeros of polynomials. One of the methods is the Graeffe Root-squaring method, which I will demonstrate using my favorite cubic, $x^3-2x-5$.<\/p><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#2c42e800-68e1-41e3-9409-1ad1da3017ff\">The Method<\/a><\/li><li><a href=\"#40cb36b1-b641-412f-98a0-808af98d36fa\">Graffe.m<\/a><\/li><li><a href=\"#c63f1707-e884-43b0-b461-d9704bfa288b\">A Historic Cubic<\/a><\/li><li><a href=\"#20b84ceb-bf0f-4b25-b508-a099f13aa645\">Basic Iteration<\/a><\/li><li><a href=\"#5cf8a751-ce36-4fcc-b24d-ce9acade3302\">One Example<\/a><\/li><li><a href=\"#c3ebf386-0d26-4f7a-88f2-addddbd39284\">Math Software?<\/a><\/li><li><a href=\"#135fe385-9f05-4430-b606-024d6e9860ee\">References<\/a><\/li><\/ul><\/div><h4>The Method<a name=\"2c42e800-68e1-41e3-9409-1ad1da3017ff\"><\/a><\/h4><p>What is today often called the Graeffe Root-Squaring method was discovered independently by Dandelin, Lobacevskii, and Graeffe in 1826, 1834 and 1837.  A 1959 article by Alston Householder referenced below straightens out the history.<\/p><p>The idea is to manipulate the coefficients of a polynomial to produce a second polynomial whose roots are the squares of the roots of the first.  If the original has a dominant real root, it will become even more dominant.  When the process is iterated you eventually reach a point where the dominant root can be read off as the ratio of the first two coefficients.  All that remains is to take an n-th root to undo the iteration.<\/p><h4>Graffe.m<a name=\"40cb36b1-b641-412f-98a0-808af98d36fa\"><\/a><\/h4><p>Here is an elegant bit of  code for producing a cubic whose roots are the squares of the roots of a given cubic.<\/p><pre class=\"codeinput\">   type <span class=\"string\">graeffe<\/span>\r\n<\/pre><pre class=\"codeoutput\">\r\nfunction b = graeffe(a)\r\n% a = a(1)*x^3 + ... + a(4) is a cubic.\r\n% b = graffe(a) is a cubic whose roots are the squares of the roots of a.\r\nb = [ a(1)^2\r\n     -a(2)^2 + 2*a(1)*a(3)  \r\n      a(3)^2 - 2*a(2)*a(4) \r\n     -a(4)^2 ];\r\n<\/pre><h4>A Historic Cubic<a name=\"c63f1707-e884-43b0-b461-d9704bfa288b\"><\/a><\/h4><p>I discussed my favorite cubic, $z^3-2z-5$, in a series of posts beginning with <a href=\"https:\/\/blogs.mathworks.com\/cleve\/2015\/12\/21\/a-historic-cubic\/\">a historic cubic<\/a> last December 21st.<\/p><p>A contour plot of the magnitude of this cubic on a square region in the plane shows the dominant real root at approximately $x=2.09$ and the pair of complex conjugate roots with smaller magnitude approximately $|z|=1.55$.<\/p><pre class=\"codeinput\">    graphic\r\n<\/pre><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/graeffe_blog_01.png\" alt=\"\"> <h4>Basic Iteration<a name=\"20b84ceb-bf0f-4b25-b508-a099f13aa645\"><\/a><\/h4><p>Repeated application of the transformation essentially squares the coefficients.  So the concern is overflow.  When I first ran this years ago as a student on the Burroughs B205, I had a limited floating point exponent range and overflow was a severe constraint.  Today with IEEE doubles we're a little better off.<\/p><h4>One Example<a name=\"5cf8a751-ce36-4fcc-b24d-ce9acade3302\"><\/a><\/h4><p>Here is a run on my cubic.  I'm just showing a few significant digits of the polynomial coefficients because the important thing is their exponents.<\/p><pre class=\"codeinput\">   a = [1 0 -2 -5];\r\n   k = 0;\r\n   f = <span class=\"string\">'%5.0f %4.0f    %5.0f   %8.3e  %8.3e  %8.3e   %20.15f \\n'<\/span>;\r\n   fprintf(f,k,1,a,Inf);\r\n   <span class=\"keyword\">while<\/span> all(isfinite(a))\r\n       k = k+1;\r\n       a = graeffe(a);\r\n       n = 2^k;\r\n       r = (-a(2))^(1\/n);\r\n       fprintf(f,k,n,a,r)\r\n   <span class=\"keyword\">end<\/span>\r\n<\/pre><pre class=\"codeoutput\">    0    1        1   0.000e+00  -2.000e+00  -5.000e+00                    Inf \r\n    1    2        1   -4.000e+00  4.000e+00  -2.500e+01      2.000000000000000 \r\n    2    4        1   -8.000e+00  -1.840e+02  -6.250e+02      1.681792830507429 \r\n    3    8        1   -4.320e+02  2.386e+04  -3.906e+05      2.135184796196703 \r\n    4   16        1   -1.389e+05  2.316e+08  -1.526e+11      2.096144583759898 \r\n    5   32        1   -1.883e+10  1.125e+16  -2.328e+22      2.094553557477412 \r\n    6   64        1   -3.547e+20  -7.504e+32  -5.421e+44      2.094551481347086 \r\n    7  128        1   -1.258e+41  1.786e+65  -2.939e+89      2.094551481542327 \r\n    8  256        1   -1.582e+82  -4.203e+130  -8.636e+178      2.094551481542327 \r\n    9  512        1   -2.504e+164  -9.665e+260      -Inf      2.094551481542327 \r\n<\/pre><p>So after seven steps we have computed the dominant root to double precision accuracy, but it takes the eighth step to confirm that. And after nine steps we run out of exponent range.<\/p><h4>Math Software?<a name=\"c3ebf386-0d26-4f7a-88f2-addddbd39284\"><\/a><\/h4><p>This is a really easy example.  To make a robust polynomial root finder we would have to confront under\/overflow, scaling, multiplicities, complex roots, and higher degree.  As far as I know, no one has ever made a serious piece of mathematical software out of the Graeffe Root-squaring method.  If you know of one, please contribute a comment.<\/p><h4>References<a name=\"135fe385-9f05-4430-b606-024d6e9860ee\"><\/a><\/h4><p>Alston S. Householder, \"Dandelin, Lobacevskii, or Graeffe?,\" <i>Amer. Math. Monthly<\/i>, 66, 1959, pp. 464&#8211;466. <a href=\"http:\/\/www.jstor.org\/stable\/2310626?seq=1#\">&lt;http:\/\/www.jstor.org\/stable\/2310626?seq=1#<\/a>&gt;<\/p><script language=\"JavaScript\"> <!-- \r\n    function grabCode_297b3ce694f840dfb352c453efd0505a() {\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='297b3ce694f840dfb352c453efd0505a ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 297b3ce694f840dfb352c453efd0505a';\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 2016 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_297b3ce694f840dfb352c453efd0505a()\"><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; R2016a<br><\/p><\/div><!--\r\n297b3ce694f840dfb352c453efd0505a ##### SOURCE BEGIN #####\r\n%% The Graeffe Root-Squaring Method for Computing the Zeros of a Polynomial\r\n% At a\r\n% <http:\/\/meetings.siam.org\/sess\/dsp_programsess.cfm?SESSIONCODE=23492\r\n% minisymposium> honoring Charlie Van Loan this week during the\r\n% <http:\/\/www.siam.org\/meetings\/an16\/ SIAM Annual Meeting>,\r\n% I will describe several dubious methods for computing the\r\n% zeros of polynomials.\r\n% One of the methods is the Graeffe Root-squaring method, which I will\r\n% demonstrate using my favorite cubic, $x^3-2x-5$.\r\n\r\n%% The Method\r\n% What is today often called the Graeffe Root-Squaring method was\r\n% discovered independently by Dandelin, Lobacevskii, and Graeffe in\r\n% 1826, 1834 and 1837.  A 1959 article by Alston Householder referenced\r\n% below straightens out the history.\r\n\r\n%%\r\n% The idea is to manipulate the coefficients of a polynomial to produce\r\n% a second polynomial whose roots are the squares of the roots of the\r\n% first.  If the original has a dominant real root, it will become even\r\n% more dominant.  When the process is iterated you eventually reach a point\r\n% where the dominant root can be read off as the ratio of the first two \r\n% coefficients.  All that remains is to take an n-th root to undo the\r\n% iteration.\r\n\r\n%% Graffe.m\r\n% Here is an elegant bit of  code for producing a cubic whose roots are\r\n% the squares of the roots of a given cubic.\r\n\r\n   type graeffe\r\n   \r\n%% A Historic Cubic\r\n% I discussed my favorite cubic, $z^3-2z-5$, in a series of posts\r\n% beginning with\r\n% <https:\/\/blogs.mathworks.com\/cleve\/2015\/12\/21\/a-historic-cubic\/\r\n% a historic cubic> last December 21st.\r\n\r\n%%\r\n% A contour plot of the magnitude of this cubic on a square region\r\n% in the plane shows the dominant real root at approximately\r\n% $x=2.09$ and the pair of complex conjugate roots with smaller magnitude\r\n% approximately $|z|=1.55$.\r\n\r\n    graphic                                \r\n\r\n%% Basic Iteration\r\n% Repeated application of the transformation essentially squares the \r\n% coefficients.  So the concern is overflow.  When I first ran this years\r\n% ago as a student on the Burroughs B205, I had a limited floating point \r\n% exponent range and overflow was a severe constraint.  Today with IEEE\r\n% doubles we're a little better off.\r\n\r\n%% One Example\r\n% Here is a run on my cubic.  I'm just showing a few significant digits of\r\n% the polynomial coefficients because the important thing is their\r\n% exponents.\r\n\r\n   a = [1 0 -2 -5];\r\n   k = 0;\r\n   f = '%5.0f %4.0f    %5.0f   %8.3e  %8.3e  %8.3e   %20.15f \\n';\r\n   fprintf(f,k,1,a,Inf);\r\n   while all(isfinite(a))\r\n       k = k+1;\r\n       a = graeffe(a);\r\n       n = 2^k;\r\n       r = (-a(2))^(1\/n);\r\n       fprintf(f,k,n,a,r)\r\n   end\r\n   \r\n%%\r\n% So after seven steps we have computed the dominant root to double\r\n% precision accuracy, but it takes the eighth step to confirm that.\r\n% And after nine steps we run out of exponent range.\r\n%                                                                     \r\n\r\n%% Math Software?   \r\n% This is a really easy example.  To make a robust polynomial root\r\n% finder we would have to confront under\/overflow, scaling, multiplicities,\r\n% complex roots, and higher degree.  As far as I know, no one has ever\r\n% made a serious piece of mathematical software out of the Graeffe\r\n% Root-squaring method.  If you know of one, please contribute a comment.  \r\n\r\n%% References\r\n% Alston S. Householder, \"Dandelin, Lobacevskii, or Graeffe?,\"\r\n% _Amer. Math. Monthly_, 66, 1959, pp. 464\u00e2\u20ac\u201c466.\r\n% <http:\/\/www.jstor.org\/stable\/2310626?seq=1#\r\n% http:\/\/www.jstor.org\/stable\/2310626?seq=1#>  \r\n##### SOURCE END ##### 297b3ce694f840dfb352c453efd0505a\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img decoding=\"async\"  class=\"img-responsive\" src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/graeffe_blog_01.png\" onError=\"this.style.display ='none';\" \/><\/div><!--introduction--><p>At a <a href=\"http:\/\/meetings.siam.org\/sess\/dsp_programsess.cfm?SESSIONCODE=23492\">minisymposium<\/a> honoring Charlie Van Loan this week during the <a href=\"http:\/\/www.siam.org\/meetings\/an16\/\">SIAM Annual Meeting<\/a>, I will describe several dubious methods for computing the zeros of polynomials. One of the methods is the Graeffe Root-squaring method, which I will demonstrate using my favorite cubic, $x^3-2x-5$.... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2016\/07\/11\/the-graeffe-root-squaring-method-for-computing-the-zeros-of-a-polynomial\/\">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,5,4,16,8],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/1868"}],"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=1868"}],"version-history":[{"count":2,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/1868\/revisions"}],"predecessor-version":[{"id":1970,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/1868\/revisions\/1970"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=1868"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=1868"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=1868"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}