# Iterative calculation of a number-theoretical constant [closed]

Recently the decimal expansion of a number theoretic constant was searched for, which is the analog of the Landau-Ramanujan constant in a certain context. The constant starts 0.638909...

It can be computed as:

prec := 200;
A = N[ 2^(-1/2) 3^(-1/4), prec ];
For[ n = 1, n < 100000, n++, p = Prime[n];
If[ Mod[p, 3] == 2, A = A / Sqrt[(1 - 1/p) (1 + 1/p)] ]
]
Print[A]


I would like to know as many exact digits as possible, or say the first 100 digits.

What is the most efficient way to achieve this and at the same time ensure the validity of the digits?

## closed as off-topic by Daniel Lichtblau, Patrick Stevens, ciao, José Antonio Díaz Navas, J. M. will be back soon♦Mar 29 '18 at 5:24

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "This question cannot be answered without additional information. Questions on problems in code must describe the specific problem and include valid code to reproduce it. Any data used for programming examples should be embedded in the question or code to generate the (fake) data must be included." – J. M. will be back soon
If this question can be reworded to fit the rules in the help center, please edit the question.

• C is a system variable, so don't use it. Can you give a source for the algorithm you presented? Anyway, you don't really need to use For[]: Product[N[1/Sqrt[(1 - 1/p) (1 + 1/p)], 250], {p, Select[Prime[Range[1*^5]], Mod[#, 3] == 2 &]}]/(Sqrt[2] 3^(1/4)) – J. M. will be back soon Mar 24 '18 at 9:07
• For starters, you still haven't answered my question. Any additional background might help answer your question about "validity of the decimals". As for my code, maybe this version is more easier for you to "read and understand": Product[With[{p = Prime[k]}, If[Mod[p, 3] == 2, N[1/Sqrt[(1 - 1/p) (1 + 1/p)], 250], 1]], {k, 1, 10^5}]/(Sqrt[2] 3^(1/4))? – J. M. will be back soon Mar 24 '18 at 10:14
• @Apostrophos No need to be snippy. I am sure, if J.M. would consider his comment as a complete answer, he would have posted it as such. Really, using For is one of the worst things to use if it is about performance (which was also an aspect of your question). – Henrik Schumacher Mar 24 '18 at 10:33
• We neither know how the constant is called nor any paper where this algorithm is discussed. Usually, the one who proposes and advertises an algorithm is in charge to derive error estimates. – Henrik Schumacher Mar 24 '18 at 10:44
• If you truncate at prime p, the relative error is O(1/(p log p)). So this will need quite a lot of primes to get 100 digits correct. – Daniel Lichtblau Mar 24 '18 at 16:30

It's a good idea to switch as early as possible to finite precision. Currently, 1/ Sqrt[(1 - 1/p) (1 + 1/p)] is computed in exact arithmetic (and probably some time is wasted to decide that the Sqrt cannot be resolved exactly). As J.M. already pointed out, also For is not a performance booster. Better utilize that arithmetic operations such as Plus, Times, Mod etc. are Listable. Moreover, the rather expensive operations such as Sqrt and Divide commute with Times so that they may be applied less often. Moreover, using Subtract[1,x] and Divide[1,x] instead of 1-x and 1/x leads to a tiny speedup (this is because, e.g., y/x is interpreted by Mathematica as y * (1/x) which contains a superfluous multiplication).

ClearAll[A];
A[n_Integer?Positive, prec_: 200] := Module[{p, p2},
p = Prime[Range[n]];
p2 = N[Pick[p, Mod[p, 3], 2], prec];
N[2^(-1/2) 3^(-1/4), prec]/ Sqrt[Times @@ Subtract[1, Divide[1, p2^2]]]
]

A[100000] // AbsoluteTiming


{0.241356, 0.63890939726017130408440290791950881418165845460093007138302026669815 1893498896918805566151830415302075625388325192796690975902066185385403 03088163214273320600204154041174777934833253718149874713993272}

However, this is far away to be correct up to 100 digits: for example, A[100000] - A[1000000] is only of order 10^-8. But that's not my fault; it's the algorithm's.

The (truncated) sequence can be obtained (at a price of some performance loss) with the following function:

ClearAll[Alist];
Alist[n_Integer?Positive, prec_: 200] := Module[{p, p2, q},
p = Prime[Range[n]];
p2 = N[Pick[p, Mod[p, 3], 2], prec];
q = Sqrt[Divide[1, Subtract[1, Divide[1, N[p2^2, prec]]]]];
N[2^(-1/2) 3^(-1/4), prec] FoldList[Times, q]
]


One might try to apply extrapolation techniques (e.g. NumericalMathNSequenceLimit) to the result (the Differences of the truncated sequence seem to have a rather smooth tail), but I am not an expert in this matter. J.M. is. Maybe if you kindly ask him...

• I would use Prime[Range[n]] to exploit listability. I was going to try doing an analysis of the decay rate of the terms, but the uncooperativeness of the OP put me off. (+1 anyway.) – J. M. will be back soon Mar 24 '18 at 10:45

An approach using Fold; however, slower than Henrik's solution

ClearAll[A2];
A2[n_Integer?Positive, prec_: 200] := Module[{p, p2},
p = Prime[Range[n]];
p2 = N[Pick[p, Mod[p, 3], 2], prec];
Fold[#1/Sqrt[1 - 1/#2^2] &, N[2^(-1/2) 3^(-1/4), prec], p2]]

A2[100000] // AbsoluteTiming

(* {0.391856, \
0.6389093972601713040844029079195088141816584546009300713830202666981518934988\
969188055661518304153020756253883251927966909759020661853854030308816321427332\
0600204154041174777934833253718149874713993272} *)
`