excelbi-sum-of-sq-divs-perf-sq.py

import math
from sympy import divisor_sigma, isprime
import time

t1 = time.time()
results = []
n = 1
while len(results) < 50:
    if n > 2 and isprime(n):
        n += 1
        continue
    sosf = divisor_sigma(n, 2)
    if math.isqrt(sosf) ** 2 == sosf:
        results.append(n)
    n += 1
print(time.time() - t1)

results