50 lines
1.0 KiB
Python
50 lines
1.0 KiB
Python
from primesieve.numpy import primes
|
|
from math import isqrt
|
|
import numpy as np
|
|
|
|
max_order = 9
|
|
blocks = [primes(10**n, 10**(n + 1)) for n in range(max_order)]
|
|
|
|
def smallest_brilliant(lb):
|
|
pos = 1
|
|
root = isqrt(lb)
|
|
|
|
for blk in blocks:
|
|
n = len(blk)
|
|
if blk[-1]*blk[-1] < lb:
|
|
pos += n*(n + 1)//2
|
|
continue
|
|
|
|
i = np.searchsorted(blk, root, 'left')
|
|
i += blk[i]*blk[i] < lb
|
|
|
|
if not i:
|
|
return blk[0]*blk[0], pos
|
|
|
|
p = blk[:i + 1]
|
|
q = (lb - 1)//p
|
|
idx = np.searchsorted(blk, q, 'right')
|
|
|
|
sel = idx < n
|
|
p, idx = p[sel], idx[sel]
|
|
q = blk[idx]
|
|
|
|
sel = q >= p
|
|
p, q, idx = p[sel], q[sel], idx[sel]
|
|
|
|
pos += np.sum(idx - np.arange(len(idx)))
|
|
return np.min(p*q), pos
|
|
|
|
res = []
|
|
p = 0
|
|
for i in range(100):
|
|
p, _ = smallest_brilliant(p + 1)
|
|
res.append(p)
|
|
|
|
print(f'first 100 are {res}')
|
|
|
|
for i in range(max_order*2):
|
|
thresh = 10**i
|
|
p, pos = smallest_brilliant(thresh)
|
|
print(f'Above 10^{i:2d}: {p:20d} at #{pos}')
|