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

Question

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?

Answer 1

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