Sum Of Squares Of Divisors Is Square

May 24, 2019

We break down the task into two parts: computing the sum of the squares of the divisors, and recognizing squares.

For the first half of the task, the divisor function sigma(k,n) from the realm of number theory returns the sum of the k‘th powers of the divisors of n; if k = 0, it returns the count of the divisors:

def sigma(k, n):
  def add(s, p, m):
    if k == 0: return s*(m+1)
    s *= (p**(k*(m+1))-1)
    return s / (p**k-1)
  fs = factors(n)
  p,m,s = fs.pop(0),1,1
  while len(fs) > 0:
    f = fs.pop(0)
    if f  p:
      s,p,m = add(s,p,m),f,1
    else: m += 1
  return add(s, p, m)

If n is large, that’s much faster than a brute-force list comprehension to compute the divisors, followed by squaring and summing. If n isn’t too big (up to about 10**12), trial division with a 2,3,5 prime wheel is a simple algorithm for finding the factors of n:

def factors(n):
  wheel = [1,2,2,4,2,4,2,4,6,2,6]
  w, f, fs = 0, 2, []
  while f * f <= n:
    while n % f == 0:
      fs.append(f); n /= f
    f, w = f + wheel[w], w + 1
    if w == 11: w = 3
  if n == 1: return fs
  else: return fs + [n]

Thus, factors(42) returns [2,3,7] and sigma(2,42) returns 2500. If your n is larger you will need a better factoring algorithm.

The other half of the problem is recognizing squares. You could use the built-in sqrt function, but that works on floating-point numbers and is thus inexact, so you probably want something better. Isaac Newton gave an algorithm for computing integer square roots by a method of successive approximations of the derivative:

def isqrt(n):
    x = n; y = (x + 1) // 2
    while y < x:
        x=y; y=(x+n//x)//2
    return x

This function returns the maximal x for which x² ≤ n. To identify squares, we could calculate the square root and multiply to check if the number is a square, but calculating square roots is expensive, so it pays to filter obvious non-squares before performing the square root calculation:

def isSquare(n):
  if 33751571>>(n%32)&1==0 or \
     38348435>>(n%27)&1==0 or \
     19483219>>(n%25)&1==0 or \
       199411>>(n%19)&1==0 or \
       107287>>(n%17)&1==0 or \
         5659>>(n%13)&1==0 or \
          571>>(n%11)&1==0 or \
           23>>(n% 7)&1==0:
    return False
  s = isqrt(n)
  if s * s == n: return s
  return False

The magic numbers are a set of bloom filters that discard numbers that aren’t quadratic residues; collectively, they identify 99.82% of non-squares before making the expensive square root calculation. The bloom filters are calculated by the following program:

def q(n):
from sets import Set
s, sum = Set(), 0
for x in xrange(0,n):
  t = pow(x,2,n)
  if t not in s:
    sum += pow(2,t)
return sum

Thus, for instance, q(32) = 33751571, and the quadratic residues (the set s that is computed inside the function) are 0, 1, 4, 9, 16, 17 and 25. Only 7/32 = 21.875% of non-square numbers survive the test. Combined with the other tests, only 0.18% of non-square numbers are not caught by one of the filters, so most of the calls to the expensive square root function are simply confirming that the number is square.

With all of that, it is easy to perform the required task:

for x in xrange(m,n):
  if x == 1:
    print [1,1]
    s = sigma(2,x)
    if isSquare(s):
      print [x,s]

Notice that 1 requires special handling, because 1 has no prime factors (it is neither prime nor composite). For m = 1, n = 250 we get:

[1, 1]
[42, 2500]

What a fun task! It requires a little bit of number theory (the sigma function) and some clever coding (bloom filters) to create a sufficiently fast program. It appears as sequence A046655 in the OEIS and as Problem 211 at Project Euler. You can run the program at

By the way, if the difference between m and n is large, instead of factoring each number individually you could sieve to find the factors. I give a program to sieve for the least prime factor of a range of number at my blog; you could modify that program to collect all the factors instead of just the least factor.


Pages: 1 2

2 Responses to “Sum Of Squares Of Divisors Is Square”

  1. Zack said

    Nice one. Here is my take in Julia (1.0):

    function IsSquare(x::Int64)
    sr = sqrt(x)
    return round(sr)^2 == x

    function FindDivisors(x::Int64)
    d = falses(x)
    d[1] = d[end] = true
    z = div(x, 2)

    if 2*z == x
        d[2] = d[z] = true
    for i = 3:(z-1)
        z = div(x, i)
        if i*z == x; d[i] = d[z] = true; end
    temp = (1:x).*d
    return temp[temp.>0]


    function test(x::Int64)
    d = FindDivisors(x)
    n = length(d)

    for i = 1:n
        d[i] *= d[i]
    s = sum(d)
    return IsSquare(s), s


    function main(m::Int64, n::Int64)
    d = n – m + 1
    Q = Array{Array{Int64, 1}}(undef, d)
    c = 0

    for x = m:n
        y, z = test(x)
        if y
            c += 1
            Q[c] = [x, z]
    return Q[1:c]


    julia> main(1,1000)
    5-element Array{Array{Int64,1},1}:
    [1, 1]
    [42, 2500]
    [246, 84100]
    [287, 84100]
    [728, 722500]

    The functions could use some optimization, but for small numbers the performance is decent. Cheers

  2. Milbrae said

    Some Pari/GP

    sosodis(m, n) = for (x = m, n, s = sigma(x, 2); if (issquare(s), print ("(", x, ", ", s, ")")))
    print (sosodis(1, 1000))

    Output is just like Zack’s

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: