fork download
copy 
  1. # sum of squares of divisors is a square
  2.  
  3. def isqrt(n): # newton
  4.     x = n; y = (x + 1) // 2
  5.     while y < x:
  6.         x=y; y=(x+n//x)//2
  7.     return x
  8.  
  9. def isSquare(n):
  10.   # def q(n):
  11.   # from sets import Set
  12.   # s, sum = Set(), 0
  13.   # for x in xrange(0,n):
  14.   #   t = pow(x,2,n)
  15.   #   if t not in s:
  16.   #     s.add(t)
  17.   #     sum += pow(2,t)
  18.   # return sum
  19.   # q(32) => 33751571
  20.   # q(27) => 38348435
  21.   # q(25) => 19483219
  22.   # q(19) =>   199411
  23.   # q(17) =>   107287
  24.   # q(13) =>     5659
  25.   # q(11) =>      571
  26.   # q(7)  =>       23
  27.   # 99.82% of non-squares
  28.   # caught by filters before
  29.   # square root calculation
  30.   if 33751571>>(n%32)&1==0 or \
  31.      38348435>>(n%27)&1==0 or \
  32.      19483219>>(n%25)&1==0 or \
  33.        199411>>(n%19)&1==0 or \
  34.        107287>>(n%17)&1==0 or \
  35.          5659>>(n%13)&1==0 or \
  36.           571>>(n%11)&1==0 or \
  37.            23>>(n% 7)&1==0:
  38.     return False
  39.   s = isqrt(n)
  40.   if s * s == n: return s
  41.   return False
  42.  
  43. def factors(n): # 2,3,5-wheel
  44.   wheel = [1,2,2,4,2,4,2,4,6,2,6]
  45.   w, f, fs = 0, 2, []
  46.   while f * f <= n:
  47.     while n % f == 0:
  48.       fs.append(f); n /= f
  49.     f, w = f + wheel[w], w + 1
  50.     if w == 11: w = 3
  51.   if n == 1: return fs
  52.   else: return fs + [n]
  53.  
  54. def sigma(k, n):
  55.   # sum of k'th powers of divisors of n
  56.   # set k = 0 to get count of divisors
  57.   def add(s, p, m):
  58.     if k == 0: return s*(m+1)
  59.     s *= (p**(k*(m+1))-1)
  60.     return s / (p**k-1)
  61.   fs = factors(n)
  62.   p,m,s = fs.pop(0),1,1
  63.   while len(fs) > 0:
  64.     f = fs.pop(0)
  65.     if f <> p:
  66.       s,p,m = add(s,p,m),f,1
  67.     else: m += 1
  68.   return add(s, p, m)
  69.  
  70. def task(m,n):
  71.   result = []
  72.   for x in xrange(m,n):
  73.     if x == 1:
  74.       result.append([1,1])
  75.     else:
  76.       s = sigma(2,x)
  77.       if isSquare(s):
  78.         result.append([x,s])
  79.   return result
  80.  
  81. print task(1,250)
Success #stdin #stdout 0s 23296KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
[[1, 1], [42, 2500], [246, 84100]]
https://ideone.com/hwHpiV
language:
Python (cpython 2.7.16)
created:
4 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!