Computing expressions from a number theory paper in MATLAB

I’m reading the paper Explicit estimates for the Riemann zeta function close to the 1-line (arXiv:2312.09412), and on page 6, the authors define the following expressions:

.

I would like to compute these quantities in MATLAB for given values of and k. For example, if and k = 1.5, the approximate expected outputs (based on the paper) are:

• ,
• ,
• ,
• .

 % Inputs 
t0 = 3;
k = 1.5;
w1=vpa(8/exp(1));

% Compute w2
if w1 < 1
    w2 = (1 - 1/exp(1) + log(sym(w1)) / log(log(t0))) / (sym(w1) * log(2));
else
    w2 = (1 - 1/exp(1)) / (sym(w1 )* log(2));
end

% Compute zeta(s)
s = 1 +k;
zeta_val = zeta(s);

% Compute B
B = 1 + 8 / (3 *sym( w1));

% Compute A_k
Ak = 1.546 * zeta_val * (1 + (2 + k)/t0)^(1/6) ...
    * (1 + (2 + k) / (t0 * log(t0))) ...
    * (1 + 2 * sqrt(1 + k) / t0);

% Display results
fprintf('For t0 = %.2f and k = %d:\n', t0, k);

For t0 = 3.00 and k = 1.500000e+00:

fprintf('w1 = %.8f\n', w1);

w1 = 2.94303553

fprintf('w2 = %.8f\n', w2);

w2 = 0.30986958

fprintf('zeta(s) = %.8f at s = %.8f\n', zeta_val, s);

zeta(s) = 1.34148726 at s = 2.50000000

fprintf('B = %.8f\n', B);

B = 1.90609394

fprintf('A_k = %.8f\n', Ak);

A_k = 9.99214293

My questions :

Is this implementation correct and numerically stable?

[Walter Roberson]

Theoretically more robust is the following implementation.

... It doesn't make any difference to the number of decimal places you display.

 % Inputs 
t0 = sym(3);
k = sym(1.5);
e = exp(sym(1));
w1 = (8/e);
log2 = log(sym(2));

% Compute w2
if w1 < 1
    w2 = (1 - 1/e + log(w1) / log(log(t0))) / (w1 * log2);
else
    w2 = (1 - 1/e) / (w1 * log2);
end

% Compute zeta(s)
s = 1 + k;
zeta_val = zeta(s);

% Compute B
B = 1 + 8 / (3 * w1);

% Compute A_k
Ak = sym(1546)/sym(10)^3 * zeta_val * (1 + (2 + k)/t0)^(1/sym(6)) ...
    * (1 + (2 + k) / (t0 * log(t0))) ...
    * (1 + 2 * sqrt(1 + k) / t0);

% Display results
fprintf('For t0 = %.2f and k = %d:\n', t0, k);

For t0 = 3.00 and k = 2:

fprintf('w1 = %.8f\n', w1);

w1 = 2.94303553

fprintf('w2 = %.8f\n', w2);

w2 = 0.30986958

fprintf('zeta(s) = %.8f at s = %.8f\n', zeta_val, s);

zeta(s) = 1.34148726 at s = 2.50000000

fprintf('B = %.8f\n', B);

B = 1.90609394

fprintf('A_k = %.8f\n', Ak);

A_k = 9.99214293

[John D'Errico]

Note that doing things like this:

w1=vpa(8/exp(1))

w1 = 

2.9430355293715386721942195435986

is BAD. Why? Because MATLAB computes exp(1) as a DOUBLE. VPA cannot be expected to know that 8/exp(1) ia what you think it is. The issue then is exp(1) LOOKS like the number e. But it is only an approximation, and when you need those extra digits, the extra digits you end up getting are complete garbage.

If instead you do this:

w1 = 8/exp(sym(1))

w1 = 

now w1 is a parameter where you know the true value, not a 16 decimal digit approximation to w1. Do you see the difference?

vpa(w1)

ans = 

2.9430355293715385727641901612917

Similarly, you did this:

w2 = (1 - 1/exp(1) + log(sym(w1)) / log(log(t0))) / (sym(w1) * log(2));

do you see that 1/exp(1) is a DOUBLE precision number? Again, you should not assume that MATLAB will see into your mind, and know to make that 1/exp(sym(1))