{"id":11813,"date":"2024-10-22T20:53:53","date_gmt":"2024-10-23T00:53:53","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=11813"},"modified":"2024-10-23T18:07:56","modified_gmt":"2024-10-23T22:07:56","slug":"mobius-mertens-and-redheffer","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2024\/10\/22\/mobius-mertens-and-redheffer\/","title":{"rendered":"M\u00f6bius, Mertens and Redheffer"},"content":{"rendered":"\n<div class=\"content\"><!--introduction-->\n<p>Recently, I have made <a href=\"https:\/\/blogs.mathworks.com\/cleve\/\">a series of blog posts<\/a> about Redheffer matrices and the Mertens conjecture. After each of the posts, readers and colleagues offered suggestions to speed up the calculations. Here is a summary of what I have learned.<\/p>\n<!--\/introduction-->\n<h3>Contents<\/h3>\n<div>\n<ul>\n<li>\n<a href=\"#32215cb4-c6e7-4c82-a09c-be2a0e403218\">M&ouml;bius Function<\/a>\n<\/li>\n<li>\n<a href=\"#2421ff74-677d-4dda-9728-8ad0f481f610\">Mertens Function<\/a>\n<\/li>\n<li>\n<a href=\"#fbbf8ec1-839c-4bd7-9416-f48ae24df402\">Redheffer Matrices<\/a>\n<\/li>\n<li>\n<a href=\"#2c3c2db6-e1d2-41c4-8a2b-d1efbabe86e7\">OEIS<\/a>\n<\/li>\n<li>\n<a href=\"#a6e6f284-1b69-42f0-a04f-f9e20891cdbc\">Ten Million<\/a>\n<\/li>\n<li>\n<a href=\"#888c9931-3060-4b47-8405-ed484be2f796\">Mertens Conjecture<\/a>\n<\/li>\n<li>\n<a href=\"#21f44074-f0af-4aa2-b6a3-adb1a2756ea2\">Conjecture is False<\/a>\n<\/li>\n<li>\n<a href=\"#e485344f-bda2-4d4d-9fd0-e7102df79117\">Matrices<\/a>\n<\/li>\n<li>\n<a href=\"#5ffa1b4e-3fe8-4908-a32c-b41a2d162cd6\">Sparsity<\/a>\n<\/li>\n<li>\n<a href=\"#39aa4540-d55c-47d1-882b-33b04368d495\">Computing Mertens<\/a>\n<\/li>\n<li>\n<a href=\"#25c55aa9-7b26-4eb8-ae66-fbfbefea38b1\">#1<\/a>\n<\/li>\n<li>\n<a href=\"#4c32719f-3698-4320-a518-3b9d554916f9\">#2<\/a>\n<\/li>\n<li>\n<a href=\"#1d535f36-08d5-4044-87d1-b3aa7094574a\">#3<\/a>\n<\/li>\n<li>\n<a href=\"#4b32160d-f00c-46f0-861a-4fe8cf231394\">#4<\/a>\n<\/li>\n<li>\n<a href=\"#e624ea53-1063-4c25-a75f-c186f877d82c\">#5<\/a>\n<\/li>\n<li>\n<a href=\"#75297d30-41f2-45ca-ad62-8e710261ee95\">Performance<\/a>\n<\/li>\n<\/ul>\n<\/div>\n<h4>M&ouml;bius Function<a name=\"32215cb4-c6e7-4c82-a09c-be2a0e403218\"><\/a>\n<\/h4>\n<p>\n<img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"http:\/\/blogs.mathworks.com\/cleve\/files\/Mobius.jpg\" alt=\"\"> <\/p>\n<p>The function named after the mid-nineteenth century German mathematician <a href=\"https:\/\/en.wikipedia.org\/wiki\/August_Ferdinand_M%C3%B6bius\" target=\"_blank\">August M&ouml;bius<\/a> is fundamental to the study of prime numbers. The M&ouml;bius function maps the positive integers onto -1,0 and +1.<\/p>\n<pre>  mu(n) = 0 if n has a repeated prime factor,\n        = (-1)^(number prime factors) if the factors of n are not repeated<\/pre>\n<p>Here is my code for the M&ouml;bius function. It relies on <tt>factor(n)<\/tt> which uses a sieve to find the prime factors of <tt>n<\/tt>.<\/p>\n<pre>  function mu = mobius(n)\n     % mu = mobius(n) returns mu(1:n).\n     mu = -eye(1,n);\n     for k = 2:n\n         f = factor(k);\n         d = diff([1 f]);\n         if any(d == 0)\n             mu(k) = 0;\n         else\n             mu(k) = (-1)^length(f);\n         end\n     end\n  end<\/pre>\n<p>Here is a graph of <tt>mu(n)<\/tt> for <tt>n<\/tt> = 1:40. For example, <tt>mu(29:31)<\/tt> are all -1 because 29 and 31 are both prime and 30 has an odd number of prime factors, 2, 3 and 5.<\/p>\n<pre class=\"language-matlab\">\n    mu = moebius(40);\n    plot(1:40,mu,<span class=\"string\">'.-'<\/span>,<span class=\"string\">'linewidth'<\/span>,1.5,<span class=\"string\">'markersize'<\/span>,16)\n    axis([-242-1.51.5])\n    title(<span class=\"string\">'M&ouml;bius function'<\/span>)\n<\/pre>\n<img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"http:\/\/blogs.mathworks.com\/cleve\/files\/Mertens_finale_blog_01.png\" alt=\"\"> <h4>Mertens Function<a name=\"2421ff74-677d-4dda-9728-8ad0f481f610\"><\/a>\n<\/h4>\n<p>\n<img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"http:\/\/blogs.mathworks.com\/cleve\/files\/Mertens.jpg\" alt=\"\"> <\/p>\n<p>The Mertens function is named after the late-nineteenth century Polish mathematician <a href=\"https:\/\/en.wikipedia.org\/wiki\/Franz_Mertens\" target=\"_blank\">Franz Mertens<\/a>. The function, which we denote by <tt>M(n)<\/tt>, is simply the partial sums of the M&ouml;bius function. The MATLAB code is very short.<\/p>\n<pre>  function M = mertens(n)\n     % M = mertens(n) returns Mertens(1:n).\n     mu = moebius(n);\n     M = cumsum([1 mu(2:n)]);\n  end<\/pre>\n<p>Here is a graph of <tt>M(n)<\/tt> for <tt>n<\/tt> = 1:40.<\/p>\n<pre class=\"language-matlab\">\n    M = mertens(40);\n    plot(1:40,M,<span class=\"string\">'.-'<\/span>,<span class=\"string\">'linewidth'<\/span>,1.5,<span class=\"string\">'markersize'<\/span>,16)\n    axis([-242-4.51.5])\n    title(<span class=\"string\">'Mertens function'<\/span>)\n<\/pre>\n<img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"http:\/\/blogs.mathworks.com\/cleve\/files\/Mertens_finale_blog_02.png\" alt=\"\"> <h4>Redheffer Matrices<a name=\"fbbf8ec1-839c-4bd7-9416-f48ae24df402\"><\/a>\n<\/h4>\n<p>\n<img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"http:\/\/blogs.mathworks.com\/cleve\/files\/Redheffer.jpg\" alt=\"\"> <\/p>\n<p>The twentieth-century American mathematician <a href=\"https:\/\/en.wikipedia.org\/wiki\/Raymond_Redheffer\" target=\"_blank\">Ray Redheffer<\/a>, who was a professor at UCLA for 55 years, had a wide range of interests. He was the author of several popular textbooks, the inventor of the electronic game <a href=\"https:\/\/en.wikipedia.org\/wiki\/Nim\" target=\"_blank\">Nim<\/a> and, with Ray and Charles Eames, the creator of a four-meter long poster <a href=\"https:\/\/www.worthpoint.com\/worthopedia\/1966-ray-charles-eames-office-men-1724813663\" target=\"_blank\">Men of Modern Mathematics<\/a>.<\/p>\n<p>The n-by-n Redheffer's matrix <tt>R(n)<\/tt> is a matrix of zeros and ones. For <tt>j &gt; 1<\/tt>, the element <tt>R(k,j)<\/tt> equals one if <tt>j<\/tt> is a divisor of <tt>k<\/tt>. And importantly, for <tt>j = 1<\/tt>, the first column <tt>R(:,1)<\/tt> is all ones.<\/p>\n<pre>  function R = redheffer(n)\n     k = 1:n;\n     j = k';\n     R = (mod(k,j) == 0);\n     R(:,1) = 1;\n  end<\/pre>\n<p>Nick Higham put Redheffer's matrix in the <a href=\"https:\/\/blogs.mathworks.com\/cleve\/2019\/06\/24\/bohemian-matrices-in-the-matlab-gallery\/\">MATLAB Gallery<\/a> several years ago. Here is a spy plot of the 200-by-200.<\/p>\n<pre class=\"language-matlab\">\n    R = gallery(<span class=\"string\">'redheff'<\/span>,200);\n    spy(R)\n    title(<span class=\"string\">'Redheffer'<\/span>)\n<\/pre>\n<img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"http:\/\/blogs.mathworks.com\/cleve\/files\/Mertens_finale_blog_03.png\" alt=\"\"> <h4>OEIS<a name=\"2c3c2db6-e1d2-41c4-8a2b-d1efbabe86e7\"><\/a>\n<\/h4>\n<p>In the Online Encyclopedia of Integer Sequences, the Mertens function is <a href=\"http:\/\/oeis.org\/A002321\" target=\"_blank\">OEIS\/A002321<\/a>. One of the many cool features of the OEIS is \"listen\", which generates music driven by the sequences. Take a look -- and a listen -- to my 63 second movie about A002321.<\/p>\n<p>\n<a href=\"https:\/\/blogs.mathworks.com\/cleve\/files\/mertens_plot.mp4\">https:\/\/blogs.mathworks.com\/cleve\/files\/mertens_plot.mp4<\/a>\n<\/p>\n<h4>Ten Million<a name=\"a6e6f284-1b69-42f0-a04f-f9e20891cdbc\"><\/a>\n<\/h4>\n<p>Here are plots of <tt>M(1:n)<\/tt> with <tt>n<\/tt> ranging from 100 to 10 million. Each plot after the first also shows the range of the previous plot. I will discuss the red lines in the next section<\/p>\n<pre class=\"language-matlab\">\n    load<span class=\"string\"> mertens<\/span><span class=\"string\"> M<\/span>\n    tiledlayout(2,3)\n    <span class=\"keyword\">for<\/span> n = 10.^(2:7)\n        nexttile\n        mertens_plot(M(1:n))\n    <span class=\"keyword\">end<\/span>\n<\/pre>\n<img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"http:\/\/blogs.mathworks.com\/cleve\/files\/Mertens_finale_blog_04.png\" alt=\"\"> <h4>Mertens Conjecture<a name=\"888c9931-3060-4b47-8405-ed484be2f796\"><\/a>\n<\/h4>\n<p>The red lines above are plots of <tt>&plusmn;sqrt(n)<\/tt>. They clearly bound <tt>M(n)<\/tt> for <tt>n<\/tt> up to 10 million. The Mertens Conjecture is that this holds for all <tt>n<\/tt>.<\/p>\n<pre>  |M(n)| &lt; sqrt(n) for all n &gt; 0<\/pre>\n<p>The conjecture appears in a 1885 letter from Thomas Stieltjes to Charles Hermite and in a 1897 paper by Mertens.<\/p>\n<p>\n<a href=\"https:\/\/en.wikipedia.org\/wiki\/Mertens_conjecture\" target=\"_blank\">Wikipedia<\/a> says<\/p>\n<pre>  It is a striking example of a mathematical conjecture proven false\n  despite a large amount of computational evidence in its favor.<\/pre>\n<h4>Conjecture is False<a name=\"21f44074-f0af-4aa2-b6a3-adb1a2756ea2\"><\/a>\n<\/h4>\n<p>In 1985, 100 years after the Stieltjes letter, Andrew Odlyzko and Herman te Riele proved that the conjecture is false. Their proof is indirect. They prove the existence of infinitely many <tt>n<\/tt> for which<\/p>\n<pre>  |M(n)|\/sqrt(n) &gt; 1.06<\/pre>\n<p>But they do not know the value of any particular <tt>n<\/tt>. They informally estimate that such an <tt>n<\/tt> would be greater than 10^30 and probably much larger.<\/p>\n<h4>Matrices<a name=\"e485344f-bda2-4d4d-9fd0-e7102df79117\"><\/a>\n<\/h4>\n<p>Let's get back to the world of matrices. It is not obvious, but is true, that the determinants of Redheffer matrices are equal to the Mertens function.<\/p>\n<pre>  det(R(n)) = M(n)<\/pre>\n<p>So, if I could have proved that<\/p>\n<pre>  |det(R(n))| &lt; sqrt(n) for all n &gt; 0<\/pre>\n<p>I would have had a proof of the Riemann Hypothesis. [Added 10\/23. See the comment and reply below.]<\/p>\n<p>It might appear that I am out of the clutches of number theory and safely back to matrix computation. But that illusion does not last for long.<\/p>\n<h4>Sparsity<a name=\"5ffa1b4e-3fe8-4908-a32c-b41a2d162cd6\"><\/a>\n<\/h4>\n<p>The <tt>0(n^3)<\/tt> memory required by the Redheffer matrix from the gallery runs out of steam quickly. We need to take advantage of sparsity. This function generates the sparse representation of a Redheffer matrix directly.<\/p>\n<pre>  function R = spredheffer(n)\n      j = (1:n)';\n      k = ones(n,1);\n      m = n;\n      for i = 2:n\n          t = [1 i:i:n];\n          p = length(t);\n          j(m+(1:p)) = t;\n          k(m+(1:p)) = i;\n          m = m+p;\n      end\n      R = sparse(k,j,1,n,n);\n  end<\/pre>\n<h4>Computing Mertens<a name=\"39aa4540-d55c-47d1-882b-33b04368d495\"><\/a>\n<\/h4>\n<p>Here are five methods for computing the Mertens function.<\/p>\n<h4>#1<a name=\"25c55aa9-7b26-4eb8-ae66-fbfbefea38b1\"><\/a>\n<\/h4>\n<p>The first simply takes the determinant of the Redheffer matrix in the <tt>gallery<\/tt>.<\/p>\n<pre>  M = det(gallery('redheff',n))<\/pre>\n<h4>#2<a name=\"4c32719f-3698-4320-a518-3b9d554916f9\"><\/a>\n<\/h4>\n<p>The sparse Gaussian elimination function <tt>lu<\/tt> with one input and four sparse outputs is designed for solving sparse linear systems. Written primarily by Tim Davis and included in his UMFPACK package, the function uses an unsymmetric pattern multifrontal pivoting strategy to reduce fill-in while maintaining numerical stability. The determinant of the input matrix is the product of the four determinants of the output matrices. Two them are triangular and two are permutations, so it is easy, and quick, to compute their determinants.<\/p>\n<pre>  R = spredheffer(n);\n  [L,U,P,Q] = lu(R)\n  M = det(L)*det(U)*det(P)*det(Q);<\/pre>\n<h4>#3<a name=\"1d535f36-08d5-4044-87d1-b3aa7094574a\"><\/a>\n<\/h4>\n<p>Pat Quillen realized that by interchanging the first and last columns, there would be little fill-in. We need only one determinant.<\/p>\n<pre>  R = spredheffer(n);\n  R(:,[1 n]) = R(:,[n 1]);\n  M = -det(R);<\/pre>\n<h4>#4<a name=\"4b32160d-f00c-46f0-861a-4fe8cf231394\"><\/a>\n<\/h4>\n<p>A thoughtful reader of the blog submitted a suggestion based on the Schur complement. Suppose <tt>E<\/tt> is a block matrix,<\/p>\n<pre>  E = [A B\n       C D]<\/pre>\n<p>Schur's formula for the determinant of <tt>E<\/tt> is<\/p>\n<pre>  det(E) = det(D)*det(A - B*(D\\C))<\/pre>\n<p>Apply this to <tt>R(n)<\/tt> with <tt>A<\/tt> the (1,1) entry, which is 1, and <tt>D<\/tt> the lower (n-1)-by-(n-1) submatrix, which is upper triangular with ones on the diagonal and determinant equal 1. This leads to method #4 which uses backslash with a sparse, unit, upper triangular matrix. There is a Redheffer matrix, but no determinant.<\/p>\n<pre>  S = spredheffer(n);\n  e = ones(n-1,1);\n  M = 1 - e'*(S(2:n,2:n)\\e);<\/pre>\n<h4>#5<a name=\"e624ea53-1063-4c25-a75f-c186f877d82c\"><\/a>\n<\/h4>\n<p>Once we have generated <tt>R(n)<\/tt> and computed <tt>M(n)<\/tt>, how do we get <tt>R(n+1)<\/tt> and <tt>M(n+1)<\/tt>? After several iterations of this approach, I found myself without any matrices and only the original definitions of M&ouml;bius and Mertens.<\/p>\n<pre>  mu = mobius(n);\n  M = cumsum([1 mu(2:n)]);<\/pre>\n<h4>Performance<a name=\"75297d30-41f2-45ca-ad62-8e710261ee95\"><\/a>\n<\/h4>\n<p>Let's summarize the five methods. The first four generate a single, isolated value of <tt>M(n)<\/tt>. Method #5 generates <tt>M(n)<\/tt> for all <tt>n<\/tt> from 1 to any specified maximum.<\/p>\n<pre>        matrix     function    dets    M<\/pre>\n<pre>  #1    full       gallery       1     1\n  #2    sparse     lu            4     1\n  #3    sparse     swap          1     1\n  #4    sparse     \\             0     1\n  #5    none       factor        0    many<\/pre>\n<p>Time in seconds to compute <tt>M(n)<\/tt> on my Lenovo T14 laptop running Windows.<\/p>\n<pre>    n    2e4      2e5      2e6     2e7<\/pre>\n<pre>  #1    28.652     -        -       -\n  #2     0.344   21.77      -       -\n  #3     0.086    1.29    16.4\n  #4     0.055    0.65     6.7    70.1\n  #5     0.075    0.80    10.7   204.5<\/pre>\n<script language=\"JavaScript\"> <!-- \n    function grabCode_a7956cd39e4e4b2da601dae253a9ad13() {\n        \/\/ Remember the title so we can use it in the new page\n        title = document.title;\n\n        \/\/ Break up these strings so that their presence\n        \/\/ in the Javascript doesn't mess up the search for\n        \/\/ the MATLAB code.\n        t1='a7956cd39e4e4b2da601dae253a9ad13 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\n        t2='##### ' + 'SOURCE END' + ' #####' + ' a7956cd39e4e4b2da601dae253a9ad13';\n    \n        b=document.getElementsByTagName('body')[0];\n        i1=b.innerHTML.indexOf(t1)+t1.length;\n        i2=b.innerHTML.indexOf(t2);\n \n        code_string = b.innerHTML.substring(i1, i2);\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\n\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \n        \/\/ in the XML parser.\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\n        \/\/ doesn't go ahead and substitute the less-than character. \n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\n\n        copyright = 'Copyright 2024 The MathWorks, Inc.';\n\n        w = window.open();\n        d = w.document;\n        d.write('<pre>\\n');\n        d.write(code_string);\n\n        \/\/ Add copyright line at the bottom if specified.\n        if (copyright.length > 0) {\n            d.writeln('');\n            d.writeln('%%');\n            if (copyright.length > 0) {\n                d.writeln('% _' + copyright + '_');\n            }\n        }\n\n        d.write('<\/pre>\\n');\n\n        d.title = title + ' (MATLAB code)';\n        d.close();\n    }   \n     --> <\/script>\n<p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\">\n<br>\n<a href=\"javascript:grabCode_a7956cd39e4e4b2da601dae253a9ad13()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \n      the MATLAB code <noscript>(requires JavaScript)<\/noscript>\n<\/span><\/a>\n<br>\n<br>\n      Published with MATLAB&reg; R2024a<br>\n<\/p>\n<\/div>\n<!--\na7956cd39e4e4b2da601dae253a9ad13 ##### SOURCE BEGIN #####\n%% M\u00f6bius, Mertens and Redheffer\n% I recently have made\n% <https:\/\/blogs.mathworks.com\/cleve\/ a series of blog posts>\n% about Redheffer matrices and the Mertens conjecture.\n% After each of the posts, readers and colleagues\n% offered suggestions to speed up the calculations. \n% Here is a summary of what I have learned.\n\n%% M\u00f6bius Function\n%\n% <<Mobius.jpg>>\n%\n% The function named after the mid-nineteenth century German\n% mathematician\n% <https:\/\/en.wikipedia.org\/wiki\/August_Ferdinand_M%C3%B6bius\n% August M\u00f6bius>\n% is fundamental to the study of prime numbers.\n% The M\u00f6bius function maps the positive integers onto -1,0 and +1.\n%\n%    mu(n) = 0 if n has a repeated prime factor,\n%          = (-1)^(number prime factors) if the factors of n are not repeated\n%\n% Here is my code for the M\u00f6bius function.  It relies on |factor(n)|\n% which uses a sieve to find the prime factors of |n|.\n%\n%    function mu = mobius(n)\n%       % mu = mobius(n) returns mu(1:n). \n%       mu = -eye(1,n);\n%       for k = 2:n\n%           f = factor(k);\n%           d = diff([1 f]);\n%           if any(d == 0)\n%               mu(k) = 0;\n%           else\n%               mu(k) = (-1)^length(f);\n%           end\n%       end\n%    end\n%\n% Here is a graph of |mu(n)| for |n| = 1:40.  For example,\n% |mu(29:31)| are all -1 because 29 and 31 are both prime\n% and 30 has an odd number of prime factors, 2, 3 and 5.\n%\n     mu = moebius(40);\n     plot(1:40,mu,'.-','linewidth',1.5,'markersize',16)\n     axis([-2 42 -1.5 1.5])\n     title('M\u00f6bius function')\n\n%% Mertens Function\n%\n% <<Mertens.jpg>>\n%\n% The Mertens function is named after the late-nineteenth century Polish\n% mathematician\n% <https:\/\/en.wikipedia.org\/wiki\/Franz_Mertens Franz Mertens>.\n% The function, which we denote by |M(n)|, is simply the partial sums\n% of the M\u00f6bius function.  The MATLAB code is very short.\n%\n%    function M = mertens(n)\n%       % M = mertens(n) returns Mertens(1:n).\n%       mu = moebius(n);\n%       M = cumsum([1 mu(2:n)]);\n%    end\n%\n% Here is a graph of |M(n)| for |n| = 1:40.  \n%\n    M = mertens(40);\n    plot(1:40,M,'.-','linewidth',1.5,'markersize',16)\n    axis([-2 42 -4.5 1.5])\n    title('Mertens function')\n\n%% Redheffer Matrices\n%\n% <<Redheffer.jpg>>\n%\n% The twentieth-century American mathematician\n% <https:\/\/en.wikipedia.org\/wiki\/Raymond_Redheffer Ray Redheffer>,\n% who was a professor at UCLA for 55 years,\n% had a wide range of interests.\n% He was the author of several popular textbooks, the inventor\n% of the electronic game <https:\/\/en.wikipedia.org\/wiki\/Nim Nim>   \n% and, with Ray and Charles Eames, the creator of a four-meter long poster\n% <https:\/\/www.worthpoint.com\/worthopedia\/1966-ray-charles-eames-office-men-1724813663\n% Men of Modern Mathematics>.\n%\n% The n-by-n Redheffer's matrix |R(n)| is a matrix of zeros and ones.\n% For |j > 1|, the element |R(k,j)| equals one if |j| is a divisor of |k|.\n% And importantly, for |j = 1|, the first column |R(:,1)| is all ones.\n% \n%    function R = redheffer(n)\n%       k = 1:n;\n%       j = k';\n%       R = (mod(k,j) == 0);\n%       R(:,1) = 1;\n%    end\n%\n% Nick Higham put Redheffer's matrix in the \n% <https:\/\/blogs.mathworks.com\/cleve\/2019\/06\/24\/bohemian-matrices-in-the-matlab-gallery\/\n% MATLAB Gallery> several years ago.  \n% Here is a spy plot of the 200-by-200.\n%  \n     R = gallery('redheff',200);\n     spy(R)\n     title('Redheffer')\n    \n%% OEIS\n% In the Online Encyclopedia of Integer Sequences, the Mertens function\n% is <http:\/\/oeis.org\/A002321 OEIS\/A002321>.  One of the many cool\n% features of the OEIS is \"listen\", which generates music driven by \n% the sequences.\n% Take a look REPLACE_WITH_DASH_DASH and a listen REPLACE_WITH_DASH_DASH to my 63 second movie about A002321.\n%\n% <https:\/\/blogs.mathworks.com\/cleve\/files\/mertens_plot.mp4>\n\n%% Ten Million\n% Here are plots of |M(1:n)| with |n| ranging from 100 to 10 million.\n% Each plot after the first also shows the range of the previous plot.\n% I will discuss the red lines in the next section\n%\n    load mertens M\n    tiledlayout(2,3)\n    for n = 10.^(2:7)\n        nexttile\n        mertens_plot(M(1:n))\n    end\n\n%% Mertens Conjecture\n% The red lines above are plots of |\u00b1sqrt(n)|.\n% They clearly bound |M(n)| for |n| up to 10 million.\n% The Mertens Conjecture is that this holds for all |n|.\n%\n%    |M(n)| < sqrt(n) for all n > 0\n%\n% The conjecture appears in a 1885 letter from Thomas Stieltjes to\n% Charles Hermite and in a 1897 paper by Mertens.\n%\n% <https:\/\/en.wikipedia.org\/wiki\/Mertens_conjecture Wikipedia> says\n% \n%    It is a striking example of a mathematical conjecture proven false\n%    despite a large amount of computational evidence in its favor.\n\n%% Conjecture is False\n% In 1985, 100 years after the Stieltjes letter, Andrew Odlyzko\n% and Herman te Riele proved that the conjecture is false.\n% Their proof is indirect.  They prove the existence of\n% infinitely many |n| for which\n%\n%    |M(n)|\/sqrt(n) > 1.06\n%\n% But they do not know the value of any particular |n|.  They informally\n% estimate that such an |n| would be greater than 10^30 and probably\n% much larger.\n\n%% Matrices\n% Let's get back to the world of matrices.  It is not obvious,\n% but is true, that the determinants of Redheffer matrices\n% are equal to the Mertens function.\n% \n%    det(R(n)) = M(n)\n%\n% So, if I could have proved that \n%\n%    |det(R(n))| < sqrt(n) for all n > 0\n%\n% I would have had a proof of the Riemann Hypothesis.\n% [Added 10\/23. See the comment and reply below.]\n%\n% It might appear that I am out of the clutches of number theory\n% and safely back to matrix computation.  But that illusion does\n% not last for long.\n\n%% Sparsity\n% The |0(n^3)| memory required by the Redheffer matrix \n% from the gallery runs out of steam quickly.\n% We need to take advantage of sparsity.\n% This function generates the sparse representation of \n% a Redheffer matrix directly.\n%\n%    function R = spredheffer(n)\n%        j = (1:n)';\n%        k = ones(n,1);\n%        m = n;\n%        for i = 2:n\n%            t = [1 i:i:n];\n%            p = length(t);\n%            j(m+(1:p)) = t;\n%            k(m+(1:p)) = i;\n%            m = m+p;\n%        end\n%        R = sparse(k,j,1,n,n);\n%    end\n%\n%% Computing Mertens\n% Here are five methods for computing the Mertens function.\n%\n\n%% #1\n% The first simply takes the determinant of the Redheffer matrix\n% in the |gallery|.\n%\n%    M = det(gallery('redheff',n))\n%\n\n%% #2\n% The sparse Gaussian elimination function |lu| with one\n% input and four sparse outputs is designed for solving sparse linear\n% systems.  Written primarily\n% by Tim Davis and included in his UMFPACK package, the function\n% uses an unsymmetric pattern multifrontal pivoting strategy \n% to reduce fill-in while maintaining numerical stability.\n% The determinant of the input matrix is the product of the four\n% determinants of the output matrices. Two them are triangular and two\n% are permutations, so it is easy, and quick, to compute their determinants.\n%\n%    R = spredheffer(n);\n%    [L,U,P,Q] = lu(R)\n%    M = det(L)*det(U)*det(P)*det(Q);\n%\n\n%% #3\n% Pat Quillen realized that by interchanging the first and last\n% columns, there would be little fill-in.  We need only one determinant.\n%\n%    R = spredheffer(n);\n%    R(:,[1 n]) = R(:,[n 1]);\n%    M = -det(R);\n%\n\n%% #4\n%\n% A thoughtful reader of the blog submitted a suggestion based on the\n% Schur complement.  Suppose |E| is a block matrix, \n%\n%    E = [A B\n%         C D]\n% \n% Schur's formula for the determinant of |E| is\n%\n%    det(E) = det(D)*det(A - B*(D\\C))\n%\n% Apply  this to |R(n)| with |A| the (1,1) entry, which is 1,\n% and |D| the lower (n-1)-by-(n-1) submatrix, which is upper triangular\n% with ones on the diagonal and determinant equal 1.  This leads to method\n% #4 which uses backslash with a sparse, unit, upper triangular matrix.\n% There is a Redheffer matrix, but no determinant.\n%\n%    S = spredheffer(n);\n%    e = ones(n-1,1);\n%    M = 1 - e'*(S(2:n,2:n)\\e);\n%\n \n%% #5\n% Once we have generated |R(n)| and computed |M(n)|, how do we get\n% |R(n+1)| and |M(n+1)|?  After several iterations of this approach,\n% I found myself without any matrices and only the original definitions\n% of M\u00f6bius and Mertens.\n%\n%    mu = mobius(n);\n%    M = cumsum([1 mu(2:n)]);\n%\n\n%% Performance\n% Let's summarize the five methods.  The first four generate\n% a single, isolated value of |M(n)|.  Method #5 generates\n% |M(n)| for all |n| from 1 to any specified maximum.\n%\n%          matrix     function    dets    M\n%\n%    #1    full       gallery       1     1\n%    #2    sparse     lu            4     1\n%    #3    sparse     swap          1     1\n%    #4    sparse     \\             0     1\n%    #5    none       factor        0    many\n%\n% Time in seconds to compute |M(n)| on my Lenovo T14 laptop\n% running Windows.\n%\n%      n    2e4      2e5      2e6     2e7\n%\n%    #1    28.652     -        -       - \n%    #2     0.344   21.77      -       - \n%    #3     0.086    1.29    16.4       \n%    #4     0.055    0.65     6.7    70.1\n%    #5     0.075    0.80    10.7   204.5\n\n\n##### SOURCE END ##### a7956cd39e4e4b2da601dae253a9ad13\n-->\n","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img src=\"https:\/\/blogs.mathworks.com\/cleve\/files\/mertens_16M.png\" class=\"img-responsive attachment-post-thumbnail size-post-thumbnail wp-post-image\" alt=\"\" decoding=\"async\" loading=\"lazy\" \/><\/div><!--introduction-->\n<p>Recently, I have made <a href=\"https:\/\/blogs.mathworks.com\/cleve\/\">a series of blog posts<\/a> about Redheffer matrices and the Mertens conjecture. After each of the posts, readers and colleagues offered suggestions to speed up the calculations. Here is a summary of what I have learned.... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2024\/10\/22\/mobius-mertens-and-redheffer\/\">read more >><\/a><\/p>","protected":false},"author":78,"featured_media":11771,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[23,4,16,8,14,36],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/11813"}],"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=11813"}],"version-history":[{"count":6,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/11813\/revisions"}],"predecessor-version":[{"id":11831,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/11813\/revisions\/11831"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media\/11771"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=11813"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=11813"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=11813"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}