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

6

1

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?

J42161217

Posted 2017-03-27T13:29:43.967

Reputation: 6 355

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 – 2017-03-27T14:05:07.150

I just know that many people here can come up with a better algorithm than mine...including you! – J42161217 – 2017-03-27T14:09:20.023

Answers

7

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

Coolwater

Posted 2017-03-27T13:29:43.967

Reputation: 17 733

wow!this is it! can anybody give me the result for 10^16? I ran out off memory.. – J42161217 – 2017-03-27T16:49:51.610

3@Jenny_mathy FYI, the answer for 10^16 is 61685027506808493591863220290146. It took about 24GB of RAM to calculate – vapor – 2017-03-27T18:42:48.253

@happy fish thank you very much happy fish. I would never be able to compute this! thank you all – J42161217 – 2017-03-27T20:07:01.480