Code wars. Integers: Recreation One. I have error "Execution Timed Out"

Question

I try to do task from codewars.

Divisors of 42 are : 1, 2, 3, 6, 7, 14, 21, 42. These divisors squared are: 1, 4, 9, 36, 49, 196, 441, 1764. The sum of the squared divisors is 2500 which is 50 * 50, a square!

Given two integers m, n (1 <= m <= n) we want to find all integers between m and n whose sum of squared divisors is itself a square. 42 is such a number.

The result will be an array of arrays, each subarray having two elements, first the number whose squared divisors is a square and then the sum of the squared divisors.

Examples:

list_squared(1, 250) --> [[1, 1], [42, 2500], [246, 84100]]

list_squared(42, 250) --> [[42, 2500], [246, 84100]]

above is the instructor of the question.

My code have passed all tests, but have a error Execution Timed Out, maybe my code is not optimized. Maybe someone help me with optimize in my code.

Here is my code

import math
def list_squared(m, n):
    number = 0
    result = []
    if m==1:
        result.append([1,1])
    else:
        number = m
        squared_list(number)
        result += squared_list(number)

    for i in range(m + 1, n + 1):
        number = i
        squared_list(number)
        result += squared_list(number)

    return result

def squared_list(number):
    array, arrays = [], []
    x = sum([i*i for i in range(1, number//2+1)if not number%i]) + number*number

    if math.sqrt(x) % 1 == 0:
        array.append(number)
        array.append(x)
        arrays.append(array)

    return arrays

Answer 1

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:
    s.add(t)
    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:

def task(m,n):
  result = []
  for x in xrange(m,n):
    if x == 1:
      result.append([1,1])
    else:
      s = sigma(2,x)
      if isSquare(s):
        result.append([x,s])
  return result

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:

>>> task(1,250)
[[1, 1], [42, 2500], [246,84100]]

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 Ideone.

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.