1850075156 = 43012 · 43013
sqrt_n = 43012.499997093866
int(sqrt_n) = 43012
range(2, int(sqrt_n)) = [2, 43011]
В цикле не проверяется делитель 43012. Ошибка.
18 минут
Дифы читать умеете? Вот так можно исправить ваш код. Работать будет 18 минут:
$ diff original.py fixed.py
7,8c7,8
< counter += 1
< for i in range(2, int(sqrt_n)):
---
> counter -= 1
> for i in range(2, int(sqrt_n) + 1):
30 секунд
Через разложение на простые задача решается за полминуты. Это лучше чем прямая проверка делителей. И теория за этим красивая, я её повторять не буду, всё есть в ответе Harry. Этот код хорош тем что factors
- процедура общего назначения, возвращает простые множители числа со степенями, а часть оптимизации спрятана внутри main
- цикл прерывается как только становится ясно, что число чётных делителей будет чётным:
import math
def factors(n):
i = 2
j = n
j_sqrt = math.isqrt(j)
while i <= j_sqrt:
if j % i == 0:
e = 0
while j % i == 0:
j //= i
e += 1
yield i, e
j_sqrt = math.isqrt(j)
i += 1 if i == 2 else 2
if j > 1:
yield j, 1
def main():
m = 0
k = 0
while True:
it = factors(1_850_000_000 + k)
p, e = next(it)
assert p == 2
if e % 2 != 0:
d = e
for _, e in it:
if e % 2 == 1:
break
d *= e + 1
else:
print(k, d)
m += 1
if m >= 5:
break
k += 2
main()
$ time python factors.py
22792 81
144450 81
266112 27
387778 9
509448 27
real 0m27.966s
user 0m27.748s
sys 0m0.008s
1.3 секунды
Сегментированное решето Эратосфена решает задачу быстрее чем за две секунды. Звучит громко, а идея простая - группу соседних чисел можно разложить на множители быстрее, чем если разлагать каждое число в отдельности. Функция factors(n1, n2)
возвращает простые множители для всех чисел в диапазоне [n1, n2). В сравнении с разложением одного числа часть оптимизаций теряется. Зато значительно уменьшается число бесполезных проверок делимости:
import math
def factors(n1, n2):
ns = list(range(n1, n2))
i = 2
n2_sqrt = math.isqrt(n2)
while i <= n2_sqrt:
for k in range((i - n1) % i, n2 - n1, i):
if ns[k] % i == 0:
e = 1
n = ns[k] // i
while n % i == 0:
n //= i
e += 1
ns[k] = n
yield n1 + k, i, e
i += 1 if i == 2 else 2
for n, j in enumerate(ns, start=n1):
if j > 1:
yield n, j, 1
def search(n1, n2):
n_even_divs = [0] * (n2 - n1)
for n, p, e in factors(n1, n2):
if p == 2:
n_even_divs[n - n1] = e
else:
n_even_divs[n - n1] *= e + 1
for n, d in enumerate(n_even_divs, start=n1):
if d % 2 == 1:
yield n, d
def main():
base = 1_850_000_000
m = 0
k = base
step = 25000
while True:
for n, d in search(k, k + step):
print(n - base, d)
m += 1
if m >= 5:
return
k += step
main()
$ time python sieve.py
22792 81
144450 81
266112 27
387778 9
509448 27
real 0m1.254s
user 0m1.244s
sys 0m0.004s
0.6 секунды
Каждый чётный делитель числа 2n находится во взаимнооднозначном соответствии с делителем числа n. Вместо того чтобы считать чётные делители больших чисел можно считать любые делители в два раза меньших чисел. Это даст двоякий выигрыш: чисел в два раза меньше и сами числа в два раза меньше. Процедура factors
не меняется. search
становится проще, main
- чуть сложнее. Время - 0.6c:
def search(n1, n2):
n_even_divs = [1] * (n2 - n1)
for n, p, e in factors(n1, n2):
n_even_divs[n - n1] *= e + 1
for n, d in enumerate(n_even_divs, start=n1):
if d % 2 == 1:
yield n, d
def main():
base = 1_850_000_000
m = 0
k = base // 2
step = 25000
while True:
for n, d in search(k, k + step):
print(2 * (n - base // 2), d)
m += 1
if m >= 5:
return
k += step
main()
$ time python sieve2.py
22792 81
144450 81
266112 27
387778 9
509448 27
real 0m0.539s
user 0m0.524s
sys 0m0.012s
⅓ секунды
factors
и search
можно объединить. Это позволит прекратить вычисления для конкретного числа если стало известно что число делителей будет чётным. Меньше трети секунды:
import math
def search(n1, n2):
sieve = [(1, n) for n in range(n1, n2)]
i = 2
n2_sqrt = math.isqrt(n2)
while i <= n2_sqrt:
for k in range((i - n1) % i, n2 - n1, i):
d, n = sieve[k]
if d != 0 and n % i == 0:
n //= i
e = 1
while n % i == 0:
n //= i
e += 1
if e % 2 == 0:
sieve[k] = d * (e + 1), n
else:
sieve[k] = 0, 0
i += 1 if i == 2 else 2
for n, (d, j) in enumerate(sieve, start=n1):
if d != 0 and j == 1:
yield n, d
def main():
base = 1_850_000_000
m = 0
k = base // 2
step = 25000
while True:
for n, d in search(k, k + step):
print(2 * (n - base // 2), d)
m += 1
if m >= 5:
return
k += step
main()
$ time python sieve3.py
22792 81
144450 81
266112 27
387778 9
509448 27
real 0m0.312s
user 0m0.308s
sys 0m0.000s
¼ секунды
До сих пор решето прореживалось двойкой и нечётными начиная с трёх. Так как прореживание повторяется вновь и вновь для каждого сегмента, оказывается выгодно построить список простых и использовать его в решете. Со всеми оптимизациями решето надо прореживать делителями до 30418. Текущая схема сделает 15209 итераций. Схема с простыми делителями - 3285 итераций. Выигрыша в пять раз не будет - удаляются итерации в которых исполняется мало кода. Но какой-то выигрыш получится.
Простые получаются классическим решетом Эратосфена. Выбрать его размер заранее сложно, будем увеличивать решето примерно в два раза, если простых не хватило. Обидно переделывать работу заново, но даже в таком виде решение эффективно. С другой стороны, получится проще чем настоящее сегментированное решето:
import math
def primes(n1, n2):
sieve = bytearray(n2)
for i in range(2, math.isqrt(n2)):
if sieve[i] == 0:
for j in range(i * i, n2, i):
sieve[j] = 1
return (i for i in range(max(2, n1), n2) if sieve[i] == 0)
class Primes:
def __init__(self, start_n=1):
self._n1 = 0
self._n2 = max(start_n, 1)
self._primes = list(primes(self._n1, self._n2))
def __iter__(self):
i1 = 0
while True:
if i1 >= len(self._primes):
self._n1 = self._n2
self._n2 = 2 * self._n2
self._primes.extend(primes(self._n1, self._n2))
i2 = len(self._primes)
yield from (self._primes[j] for j in range(i1, i2))
i1 = i2
def search(n1, n2, ps):
sieve = [(1, n) for n in range(n1, n2)]
n2_sqrt = math.isqrt(n2)
for i in ps:
if i >= n2_sqrt:
break
for k in range((i - n1) % i, n2 - n1, i):
d, n = sieve[k]
if d != 0 and n % i == 0:
n //= i
e = 1
while n % i == 0:
n //= i
e += 1
if e % 2 == 0:
sieve[k] = d * (e + 1), n
else:
sieve[k] = 0, 0
for n, (d, j) in enumerate(sieve, start=n1):
if d != 0 and j == 1:
yield n, d
def main():
base = 1_850_000_000
m = 0
k = base // 2
step = 25000
ps = Primes(25000)
while True:
for n, d in search(k, k + step, ps):
print(2 * (n - base // 2), d)
m += 1
if m >= 5:
return
k += step
main()
$ time python sieve4.py
22792 81
144450 81
266112 27
387778 9
509448 27
real 0m0.234s
user 0m0.232s
sys 0m0.000s
if sqrt_n.is_integer(): counter += 1
проверьте что он чётный...