# by Charlie Neder, Sep 01 2018 from math import log2 from itertools import combinations limit = 10**10 bases = [2] powers = [2,7,53,151] #potential primes in the exponent #next prime would be p(128) = 719, for when limit gets pushed to 10**230+ primes = [] arr = [] for n in range(2,int(log2(limit))): #get primes <= log2(limit) ..if all(n%i for i in primes): primes.append(n) size = 1 while 1: #build array of products of distinct primes to ensure perfect powers < n are counted correctly ..res = [] ..for i in combinations(primes,size): ....prod = 1 ....for j in i: prod *= j ....if prod <= int(log2(limit)): res.append(prod) ..if len(res): ....arr.append(res) ....size += 1 ..else: break powerset = set(powers) #precompute valid exponents for nextpower(n) using powers for p in powers: ..while 1: #iteratively populate powerset with new products of the allowed primes ....ps = powerset.copy() ....for q in powerset: ......if not p*q in powerset and p*q <= int(log2(limit)): ........ps.add(p*q) ....if ps == powerset: ......break ....else: powerset = ps def powmath(n,direction): #either add or subtract number of perfect powers <= n, depending on direction ..diff = 1 ..for i in range(len(arr)): ....for j in arr[i]: ......diff += (int(n**(1/j)) - 1) * (-1)**i ..return n + diff * direction def rad(n): ..res = powmath(n,1) ..while not powmath(res,-1) == n: res += 1 ..return res #powmath(rad(n),-1) == n, by definition def nextval(n): #locate next valid e-number > n ..nxt = None ..for i in powerset: ....for b in bases: ......if nxt: ........if b**i < nxt[0] and b**i > n: nxt = (b**i,b,i) ......elif b**i > n: nxt = (b**i,b,i) ..return nxt n = 2 count = 1 while n < limit: ..n, base, exp = nextval(n) ..print(count,n) ..#print(n,rad(n),base,exp) ..bases.append(rad(n)) ..count += 1