project-euler/211_divisor_square_sum.py

43 lines
678 B
Python
Raw Normal View History

#!/usr/bin/env python3
from math import ceil
def is_square(n):
return int(n**0.5)**2 == n
def getpw(n, p):
pw = 0
while not n % p:
pw += 1
n //= p
return pw
N = 64 * 10**6
sieve = [1] * N
for p in range(2, N, 2):
sieve[p] = 2
for p in range(3, int(ceil(N**0.5 - 1)), 2):
if sieve[p] > 1:
continue
sieve[p] = p
for i in range(p**2, N, p):
sieve[i] = p
print("Sieved")
s = 1
for n in range(2, N):
if sieve[n] == 1:
sieve[n] = n
p = sieve[n]
pw = p ** getpw(n, p)
sieve[n] = sieve[n // p] + \
pw**2 * sieve[n // pw]
if is_square(sieve[n]):
s += n
print(s)