Redheffer and Mertens, Accelerated - Cleve Moler on Mathematics and Computing
 

Shortly after I published the second post about the [URL="https://blogs.mathworks.com/cleve/2024/09/27/redheffer-and-mertens-continued/"]Mertens conjecture[/URL], 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.

[B]Contents[/B]
[LIST][*][URL="https://www.labfans.com/bbs/#7d362761-695e-4474-88df-0a49bd91d0a1"]Block matrices[/URL][*][URL="https://www.labfans.com/bbs/#3237700b-50f8-4cce-a2a4-ce2a67875d0f"]redmert[/URL][*][URL="https://www.labfans.com/bbs/#c59e3dbe-0a74-4f7c-839b-a92e38cb63c3"]Inside redmert[/URL][*][URL="https://www.labfans.com/bbs/#6c005e64-ccd7-4f61-bab0-eec89c74cdd4"]mertens_plot[/URL][*][URL="https://www.labfans.com/bbs/#dc3123c4-88e4-4acc-87ef-5dce4750f795"]Postscript[/URL][/LIST][B]Block matrices[/B]

The comment references the [URL="https://en.wikipedia.org/wiki/Determinant#Block_matrices"]Wikipedia article[/URL] on block matrices.

You could also consider the matrix as a 2x2 block matrix and use the formula for the determinant of a block matrix [1]. A = redheffer(n); M = full(A(1,1) - A(1, 2:end) * (A(2:end,2:end) \ A(2:end, 1))); Since the (n-1)x(n-1) block is upper triangular, the solve becomes a back-substitution.[B]redmert[/B]

My new program is named redmert, an abbreviation of Redheffer-Mertens. It uses the fact that redheffer(n) is obtained from redheffer(n-1) by appending the last column.

Let R(n) denote the upper or right-triangular part of redheffer(n).

R(n) = triu(redheffer(n))R(n) is sparse, upper triangular and has ones on the diagonal. The indices of the nonzeros in the last column of R(n) are the factors of n. For example, here is R(8).

R8 = full(triu(redheffer(8)))R8 = 1 1 1 1 1 1 1 1 0 1 0 1 0 1 0 1 0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1The idea behind redmert is to compute a sequence of Redheffer matrices, R, and associated values of the Mertens function, M.

[M,R] = redmert(p,R)The input is a scalar integer p, the desired sequence length, and a sparse matrix R, the upper triangle of a Redheffer matrix of some order, n. The output is an integer vector of values M(n+1:n+p) and the upper triangle of the Redheffer matrix of order n+p. This output R can then be used as the input R in another call to redmert.

The sequence is started with an empty R.

For example,

[M,R] = redmert(8,[]);The output is mertens(n), n = 1:8, and R8 from the example above.

MR8 = full(R)M = 1 0 -1 -1 -2 -1 -2 -2R8 = 1 1 1 1 1 1 1 1 0 1 0 1 0 1 0 1 0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1[B]Inside redmert[/B]

The entire code for redmert is twelve lines long. It manipulates sparse matrices and uses sparse backslash to solve a triangular system. Nothing else is required.

Lines 7 and 8 generate the last column of R and lines 9 and 10 implement the new idea about block matrices.

dbtype redmert1 function [M,R] = redmert(p,Rin)2 M = zeros(p,1);3 R = sparse(triu(Rin));4 n = size(R,1);5 for q = 1:p6 n = n+1;7 k = (mod(n,1:n) == 0);8 R(k,n) = 1;9 e = ones(n-1,1);10 M(q) = R(1,1) - e'*(R(2:n,2:n)\e);11 end12 end[B]mertens_plot[/B]

It takes about five hours for redmert to compute half a million values on my laptop.

n = 0.5e6; p = 0.5e4; R = sparse([]); M = []; for k = p:p:n disp(k) [Mout,R] = redmert(p,R); M = [M; Mout]; mertens_plot(M) endmertens_plot[IMG]https://blogs.mathworks.com/cleve/files/redmert_blog_01.png[/IMG] [B]Postscript[/B]

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.

[RIGHT][COLOR=gray][I]
[URL="javascript:grabCode_3267a823a3ca4544ace75a24a1cc6636()"][I]Get the MATLAB code (requires JavaScript)[/I][/URL]

Published with MATLAB® R2024a
[/I][/COLOR][/RIGHT]