6
$\begingroup$

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?

Improve this question
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – vapor Mar 27 '17 at 14:05
  • $\begingroup$ I just know that many people here can come up with a better algorithm than mine...including you! $\endgroup$ – J42161217 Mar 27 '17 at 14:09
7
$\begingroup$

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.