{"id":11648,"date":"2024-09-30T10:28:40","date_gmt":"2024-09-30T14:28:40","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=11648"},"modified":"2024-10-22T18:54:17","modified_gmt":"2024-10-22T22:54:17","slug":"redheffer-and-mertens-accelerated","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2024\/09\/30\/redheffer-and-mertens-accelerated\/","title":{"rendered":"Redheffer and Mertens, Accelerated"},"content":{"rendered":"<div class=\"content\"><!--introduction-->\r\n<p>Shortly after I published the second post about the <a href=\"https:\/\/blogs.mathworks.com\/cleve\/2024\/09\/27\/redheffer-and-mertens-continued\/\">Mertens conjecture<\/a>, a reader's comment suggested a new approach to computing Redheffer determinants and the Mertens function. It is now possible to compute a half-million values of the Mertens function in about five hours.<\/p>\r\n<!--\/introduction-->\r\n<h3>Contents<\/h3>\r\n<div>\r\n<ul>\r\n<li>\r\n<a href=\"#7d362761-695e-4474-88df-0a49bd91d0a1\">Block matrices<\/a>\r\n<\/li>\r\n<li>\r\n<a href=\"#3237700b-50f8-4cce-a2a4-ce2a67875d0f\"><tt>redmert<\/tt><\/a>\r\n<\/li>\r\n<li>\r\n<a href=\"#c59e3dbe-0a74-4f7c-839b-a92e38cb63c3\">Inside <tt>redmert<\/tt><\/a>\r\n<\/li>\r\n<li>\r\n<a href=\"#6c005e64-ccd7-4f61-bab0-eec89c74cdd4\">mertens_plot<\/a>\r\n<\/li>\r\n<li>\r\n<a href=\"#dc3123c4-88e4-4acc-87ef-5dce4750f795\">Postscript<\/a>\r\n<\/li>\r\n<\/ul>\r\n<\/div>\r\n<h4>Block matrices<a name=\"7d362761-695e-4474-88df-0a49bd91d0a1\"><\/a>\r\n<\/h4>\r\n<p>The comment references the <a href=\"https:\/\/en.wikipedia.org\/wiki\/Determinant#Block_matrices\" target=\"_blank\" rel=\"noopener\">Wikipedia article<\/a> on block matrices.<\/p>\r\n<pre>  You could also consider the matrix as a 2x2 block matrix and\r\n  use the formula for the determinant of a block matrix [1].\r\n      A = redheffer(n);\r\n      M = full(A(1,1) - A(1, 2:end) * (A(2:end,2:end) \\ A(2:end, 1)));\r\n  Since the (n-1)x(n-1) block is upper triangular, the solve becomes\r\n  a back-substitution.<\/pre>\r\n<h4>\r\n<tt>redmert<\/tt><a name=\"3237700b-50f8-4cce-a2a4-ce2a67875d0f\"><\/a>\r\n<\/h4>\r\n<p>My new program is named <tt>redmert<\/tt>, an abbreviation of Redheffer-Mertens. It uses the fact that <tt>redheffer(n)<\/tt> is obtained from <tt>redheffer(n-1)<\/tt> by appending the last column.<\/p>\r\n<p>Let <tt>R(n)<\/tt> denote the upper or right-triangular part of <tt>redheffer(n)<\/tt>.<\/p>\r\n<pre>  R(n) = triu(redheffer(n))<\/pre>\r\n<p>\r\n<tt>R(n)<\/tt> is sparse, upper triangular and has ones on the diagonal. The indices of the nonzeros in the last column of <tt>R(n)<\/tt> are the factors of <tt>n<\/tt>. For example, here is <tt>R(8)<\/tt>.<\/p>\r\n<pre class=\"codeinput\">\r\nR8 = full(triu(redheffer(8)))\r\n<\/pre>\r\n<pre class=\"codeoutput\">R8 =\r\n     1     1     1     1     1     1     1     1\r\n     0     1     0     1     0     1     0     1\r\n     0     0     1     0     0     1     0     0\r\n     0     0     0     1     0     0     0     1\r\n     0     0     0     0     1     0     0     0\r\n     0     0     0     0     0     1     0     0\r\n     0     0     0     0     0     0     1     0\r\n     0     0     0     0     0     0     0     1\r\n<\/pre>\r\n<p>The idea behind <tt>redmert<\/tt> is to compute a sequence of Redheffer matrices, <tt>R<\/tt>, and associated values of the Mertens function, <tt>M<\/tt>.<\/p>\r\n<pre>  [M,R] = redmert(p,R)<\/pre>\r\n<p>The input is a scalar integer <tt>p<\/tt>, the desired sequence length, and a sparse matrix <tt>R<\/tt>, the upper triangle of a Redheffer matrix of some order, <tt>n<\/tt>. The output is an integer vector of values <tt>M(n+1:n+p)<\/tt> and the upper triangle of the Redheffer matrix of order <tt>n+p<\/tt>. This output <tt>R<\/tt> can then be used as the input <tt>R<\/tt> in another call to <tt>redmert<\/tt>.<\/p>\r\n<p>The sequence is started with an empty <tt>R<\/tt>.<\/p>\r\n<p>For example,<\/p>\r\n<pre class=\"codeinput\">\r\n[M,R] = redmert(8,[]);\r\n<\/pre>\r\n<p>The output is <tt>mertens(n), n = 1:8<\/tt>, and <tt>R8<\/tt> from the example above.<\/p>\r\n<pre class=\"codeinput\">\r\nMR8 = full(R)\r\n<\/pre>\r\n<pre class=\"codeoutput\">\r\nM =\r\n     1\r\n     0\r\n    -1\r\n    -1\r\n    -2\r\n    -1\r\n    -2\r\n    -2\r\nR8 =\r\n     1     1     1     1     1     1     1     1\r\n     0     1     0     1     0     1     0     1\r\n     0     0     1     0     0     1     0     0\r\n     0     0     0     1     0     0     0     1\r\n     0     0     0     0     1     0     0     0\r\n     0     0     0     0     0     1     0     0\r\n     0     0     0     0     0     0     1     0\r\n     0     0     0     0     0     0     0     1\r\n<\/pre>\r\n<h4>Inside <tt>redmert<\/tt><a name=\"c59e3dbe-0a74-4f7c-839b-a92e38cb63c3\"><\/a>\r\n<\/h4>\r\n<p>The entire code for <tt>redmert<\/tt> is twelve lines long. It manipulates sparse matrices and uses sparse backslash to solve a triangular system. Nothing else is required.<\/p>\r\n<p>Lines 7 and 8 generate the last column of <tt>R<\/tt> and lines 9 and 10 implement the new idea about block matrices.<\/p>\r\n<pre class=\"codeinput\">\r\ndbtype <span class=\"string\">redmert<\/span>\r\n<\/pre>\r\n<pre class=\"codeoutput\">\r\n1     function [M,R] = redmert(p,Rin)\r\n2         M = zeros(p,1);\r\n3         R = sparse(triu(Rin));\r\n4         n = size(R,1);\r\n5         for q = 1:p\r\n6             n = n+1;\r\n7             k = (mod(n,1:n) == 0);\r\n8             R(k,n) = 1;\r\n9             e = ones(n-1,1);\r\n10            M(q) = R(1,1) - e'*(R(2:n,2:n)\\e);\r\n11        end\r\n12    end\r\n<\/pre>\r\n<h4>mertens_plot<a name=\"6c005e64-ccd7-4f61-bab0-eec89c74cdd4\"><\/a>\r\n<\/h4>\r\n<p>It takes about five hours for <tt>redmert<\/tt> to compute half a million values on my laptop.<\/p>\r\n<pre>   n = 0.5e6;\r\n   p = 0.5e4;\r\n   R = sparse([]);\r\n   M = [];\r\n   for k = p:p:n\r\n       disp(k)\r\n       [Mout,R] = redmert(p,R);\r\n       M = [M; Mout];\r\n       mertens_plot(M)\r\n   end<\/pre>\r\n<pre class=\"codeinput\">\r\nmertens_plot\r\n<\/pre>\r\n<img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/redmert_blog_01.png\" alt=\"\"> <h4>Postscript<a name=\"dc3123c4-88e4-4acc-87ef-5dce4750f795\"><\/a>\r\n<\/h4>\r\n<p>I started this project by being surprised to find myself computing determinants. Now I am back to my long-time position disparaging determinants. They have been replaced by a good friend, backslash.<\/p>\r\n<script language=\"JavaScript\"> <!-- \r\n    function grabCode_3267a823a3ca4544ace75a24a1cc6636() {\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='3267a823a3ca4544ace75a24a1cc6636 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 3267a823a3ca4544ace75a24a1cc6636';\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 2024 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>\r\n<p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\">\r\n<br>\r\n<a href=\"javascript:grabCode_3267a823a3ca4544ace75a24a1cc6636()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n      the MATLAB code <noscript>(requires JavaScript)<\/noscript>\r\n<\/span><\/a>\r\n<br>\r\n<br>\r\n      Published with MATLAB&reg; R2024a<br>\r\n<\/p>\r\n<\/div>\r\n<!--\r\n3267a823a3ca4544ace75a24a1cc6636 ##### SOURCE BEGIN #####\r\n%% Redheffer and Mertens, Accelerated\r\n% Shortly after I published the second post about the\r\n% <https:\/\/blogs.mathworks.com\/cleve\/2024\/09\/27\/redheffer-and-mertens-continued\/\r\n% Mertens conjecture>,\r\n% a reader's comment suggested a new approach to\r\n% computing Redheffer determinants and the Mertens function.\r\n% It is now possible to compute a half-million values of the\r\n% Mertens function in about five hours.\r\n\r\n%% Block matrices\r\n% The comment references the\r\n% <https:\/\/en.wikipedia.org\/wiki\/Determinant#Block_matrices\r\n% Wikipedia article> on block matrices.\r\n%\r\n%    You could also consider the matrix as a 2x2 block matrix and\r\n%    use the formula for the determinant of a block matrix [1].\r\n%        A = redheffer(n);\r\n%        M = full(A(1,1) - A(1, 2:end) * (A(2:end,2:end) \\ A(2:end, 1)));\r\n%    Since the (n-1)x(n-1) block is upper triangular, the solve becomes\r\n%    a back-substitution.\r\n\r\n%% |redmert|\r\n% My new program is named |redmert|, an abbreviation of \r\n% Redheffer-Mertens.  It uses the fact that\r\n% |redheffer(n)| is obtained from |redheffer(n-1)| by appending\r\n% the last column.\r\n%\r\n% Let |R(n)| denote the upper or right-triangular\r\n% part of |redheffer(n)|.\r\n%\r\n%    R(n) = triu(redheffer(n))\r\n%   \r\n% |R(n)| is sparse, upper triangular and has ones on the\r\n% diagonal.  \r\n% The indices of the nonzeros in the last column of |R(n)|\r\n% are the factors of |n|.\r\n% For example, here is |R(8)|.\r\n\r\n    R8 = full(triu(redheffer(8)))\r\n\r\n%%\r\n% The idea behind |redmert| is to compute a sequence\r\n% of Redheffer matrices, |R|, and associated values of the\r\n% Mertens function, |M|.\r\n%\r\n%    [M,R] = redmert(p,R)\r\n%\r\n% The input is a scalar integer |p|, the desired sequence length,\r\n% and a sparse matrix |R|, the upper triangle of a Redheffer\r\n% matrix of some order, |n|.\r\n% The output is an integer vector of values |M(n+1:n+p)| and\r\n% the upper triangle of the Redheffer matrix of order |n+p|. \r\n% This output |R| can then be used as the input |R| in another\r\n% call to |redmert|.\r\n%\r\n% The sequence is started with an empty |R|.\r\n%\r\n% For example,\r\n%\r\n   [M,R] = redmert(8,[]);\r\n\r\n%%\r\n% The output is |mertens(n), n = 1:8|, and |R8| from the example\r\n% above.\r\n\r\n    M\r\n    R8 = full(R)\r\n\r\n%% Inside |redmert|\r\n% The entire code for |redmert| is twelve lines long.\r\n% It manipulates sparse matrices and uses sparse backslash\r\n% to solve a triangular system.  Nothing else is required.\r\n%\r\n% Lines 7 and 8 generate the last column of |R| and lines 9 and 10\r\n% implement the new idea about block matrices.\r\n\r\n      dbtype redmert\r\n\r\n%% mertens_plot\r\n% It takes about five hours for |redmert|\r\n% to compute half a million values on my laptop.\r\n% \r\n%     n = 0.5e6;\r\n%     p = 0.5e4;\r\n%     R = sparse([]);\r\n%     M = [];\r\n%     for k = p:p:n\r\n%         disp(k)\r\n%         [Mout,R] = redmert(p,R);\r\n%         M = [M; Mout];\r\n%         mertens_plot(M)\r\n%     end\r\n\r\n      mertens_plot\r\n\r\n%% Postscript\r\n% I started this project by being surprised to find\r\n% myself computing determinants.\r\n% Now I am back to my long-time position disparaging\r\n% determinants. They have been replaced by a good friend, backslash.\r\n\r\n##### SOURCE END ##### 3267a823a3ca4544ace75a24a1cc6636\r\n-->\r\n","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/redmert_plot_small.png\" class=\"img-responsive attachment-post-thumbnail size-post-thumbnail wp-post-image\" alt=\"\" decoding=\"async\" loading=\"lazy\" \/><\/div><!--introduction-->\r\n<p>Shortly after I published the second post about the <a href=\"https:\/\/blogs.mathworks.com\/cleve\/2024\/09\/27\/redheffer-and-mertens-continued\/\">Mertens conjecture<\/a>, a reader's comment suggested a new approach to computing Redheffer determinants and the Mertens function. It is now possible to compute a half-million values of the Mertens function in about five hours.... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2024\/09\/30\/redheffer-and-mertens-accelerated\/\">read more >><\/a><\/p>","protected":false},"author":78,"featured_media":11660,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[5,4,36,1],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/11648"}],"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=11648"}],"version-history":[{"count":5,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/11648\/revisions"}],"predecessor-version":[{"id":11804,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/11648\/revisions\/11804"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media\/11660"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=11648"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=11648"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=11648"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}