We can program Euler's recursion for the sum of divisors, too. It only
makes sense, when we need to compute the sum of divisors for many integers
in 1:n, because then we can reuse the already computed values. It is the
case, when we search for amicable numbers: we compute the vector of
divisor-sums only once, and then search in that vector.
Euler's algorithm is remarkable, because it does not use any
multiplications or divisions as coded below, only additions/subtractions
and indexing. (Except the computation of the constant k, but that can be
approximated with shift operations - not done in the code, or use n in
place of 2k).
The auxiliary memory use is proportional to sqrt(n), which is small
compared to the result vector of size n. The running time is O(n^1.5),
which may or may not be competitive with the method in the previous post:
the sum of divisors is computed for each input by factoring and compute a
product over the prime factors.
~~~
function allDivSums(n) # returns an Int32 vector, v[i] = sum of divisors of
i, i <= n
k = iceil((sqrt(1+24n)+1)/6)
pn = zeros(Int32,2k)
pn[1:2:2k] = cumsum([1+3(0:k-1)])
pn[2:2:2k] = cumsum([2+3(0:k-1)])
sg = ones(Int32,2k)
sg[3:4:2k] = sg[4:4:2k] = -1
DS = ones(Int32,n)
for u = 2:n
s = 0; i = 1
while i <= 2k && pn[i] < u
s += flipsign(DS[u-pn[i]],sg[i])
i += 1
end
if i <= 2k && pn[i] == u
s += flipsign(u,sg[i])
end
DS[u] = s
end
return DS
end
~~~
