# Efficient way to sum all divisors of numbers below N excluding divisors 4/d

So, I want to Sum all divisors off all numbers below N (for a big N, i.e. N=10^16) which (divisors) are NOT a multiple of 4

I tried DivisorSigma but then the following method was x2 faster

j = 10^6;
Sum[n*Floor[j/n], {n, 1, j}] -
Sum[4*m*Floor[j/(4*m)], {m, 1, Floor[j/4]}] // Timing


10^6 takes a few seconds. Can we do better? 10^16 is huge.. Can anybody suggest an algorithm that goes up to N/2 (or better) excluding every forth number?

• Looks like a difficult number theory problem. Is there any evidence that you can get the result within a reasonable time on a personal computer? – vapor Mar 27 '17 at 14:05
• I just know that many people here can come up with a better algorithm than mine...including you! – ZaMoC Mar 27 '17 at 14:09

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 0s. 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)].

• wow!this is it! can anybody give me the result for 10^16? I ran out off memory.. – ZaMoC Mar 27 '17 at 16:49
• @Jenny_mathy FYI, the answer for 10^16 is 61685027506808493591863220290146. It took about 24GB of RAM to calculate – vapor Mar 27 '17 at 18:42
• @happy fish thank you very much happy fish. I would never be able to compute this! thank you all – ZaMoC Mar 27 '17 at 20:07