Slightly faster:
j = 400001;
Total[# Floor[j/#] &[Flatten[Partition[Range[j], 3, 4, 1, {}]]]] // Timing
Very much faster:
j = 1111160000;
Total[Range[Floor[Sqrt[j]]] (With[{u = Ceiling[#, 4], v = Floor[#2, 4]},
(#2 - # + 1) (# + #2)/2 + (u - v - 4) (u + v)/8] & @@@ #)] +
Total[# Floor[j/#] &[Flatten[Partition[Range[#[[-1, 1]] - 1], 3, 4, 1, {}]]]] &[
Thread[Floor[{j/(Range[Floor[Sqrt[j]]] + 1) + 1, j/Range[Floor[Sqrt[j]]]}]]]
In the last one I say for e.g. j=600
we will encounter one case of 301
to 600
two cases of j/3
to j/2
, three cases of j/4
to j/3
etc. until Sqrt[j]
where this approach is inefficient because of alot of 0
s. From that point the approach above is used.
Edit: There was a mistake in a Part
argument, which however didn't matter for the result. To limit the memory, you can run this using an appropiate value for sz
:
j=10^16
Dynamic[{max, cur}]
Total[#, {1, 2}] &[
With[{max1 = Floor[j^(1./8 + 3/8 Tanh[Log[j]/20])], sz = 4000000},
With[{max2 = Floor[j/(max1 + 1)]}, With[{pars = MapAt[{Most[#] + 1, Rest[#]}\[Transpose] &, {
Partition[Round[Range[1, max1 + 1, max1/Ceiling[max1/sz]]], 2, 1],
Append[Round[Range[0, max2 - 1, max2/Ceiling[max2/sz]], 4], max2]}, 2]},
{max = {1, Length[pars[[1]]]}; cur = 0; Map[(++cur;
Total[Most[#] With[{x = Quotient[j, #]},
2 Differences[#] (Most[#] + Rest[#] + 1) -
Differences[x] (Most[x] + Rest[x] + 1)/2 &[
Quotient[x, 4]]] &[Range @@ #]]) &, pars[[1]]],
max = {2, Length[pars[[2]]]}; cur = 0; Map[(++cur;
Total[# Quotient[j, #] &[Drop[Range @@ #, {4, -1, 4}]]]) &, pars[[2]]]}]]]]
61685027506808493591863220290146
where max1
is merely a tuning integer. But by experiments I found this Floor[j^(1./8 + 3/8 Tanh[Log[j]/20])]
to be significantly faster than just Floor[j^(1/2)]
.