Как оптимизировать алгоритм для проверки чисел на гипотезу Коллатца?

Question

Я пишу итоговый школьный проект, пришло в голову несколько идей оптимизации:

Множество, в котором записываю все числа, которые попадались в последовательности
При попадении чётного числа в функцию, сразу выхожу из функции (т.к. любое чётное число в любом случае станет нечётным, при этом в диапазоне такое нечётное число явно было)
В случае нечётности числа, после умножении на 3 и добавлении 1, сразу делю на 2, т.к. оно в любом случае становится чётным

Есть ли ещё какие-то способы оптимизации алгоритма (именно на языке программирования python), либо мне нужно полностью переделать код?

import time

checkedNumbers = set() #Множество, где хранятся все проверенные числа

def ThreeNPlusOne(N):
    if N % 2 == 0:
        return
    
    while N != 1 and N not in checkedNumbers:
        if N % 2 == 0:
            N //= 2
        else:
            N = 3 * N + 1
            N //= 2

numbersCount = int(input())

startTime = time.monotonic()

for number in range(1, numbersCount + 1):
    ThreeNPlusOne(number)

timeResult = time.monotonic() - startTime

print(timeResult) #Вывод времени, затраченного на проверку чисел

Answer 1

Стартовая точка

Ваша программа за десять секунд на моей машине обрабатывает полтора миллиона чисел. Это хороший результат, будем от него отталкиваться.

$ echo 1_500_000 | python original.py
9.894165792968124

Использование

Подумайте что будет с программой если она отыщет число, нарушающее гипотезу Коллаца. Программа работает, а пользователь гадает: она ещё перебирает числа? Или нужное число найдено и строится бесконечно растущая цепочка чисел?

Лучше сделать по другому: программа перебирает числа без ограничения, каждую секунду печатает текущее проверяемое число. Пока числа растут, вы оцениваете как далеко зашёл поиск. Если вы видите повторяющееся число, вы понимаете, что это число, скорее всего, нарушитель гипотезы.

К предыдущей идее я добавлю стартовое n. Если программа будет прервана, то её можно будет запустить с точки, где она была последний раз.

import threading
import time

n0 = int(input())


def log():
    while True:
        time.sleep(1)
        print(f'{n0:_}')


t = threading.Thread(target=log, daemon=True)
t.start()

while True:
    n = n0
    while n > 1:
        if n % 2 == 0:
            n //= 2
        else:
            n = (3 * n + 1) // 2
    n0 += 1

$ echo 1 | python search-1.py 
91_727
173_081
250_747
^CTraceback (most recent call last):
  File "/home/sv/desk/stackoverflow/collatz/search-1.py", line 18, in <module>
    while n > 1:
KeyboardInterrupt


$ echo 250_747 | python search-1.py 
320_502
388_221
455_140
...

Эта программа за 10 секунд добирается до семисот тысяч. В два раза хуже вашей.

Множество проверенных чисел

Это хорошая идея, благодаря которой ваша программа такая быстрая. К сожалению, это множество быстро растёт: в нём все числа меньше текущего и ещё сколько-то сверху. Когда программа уйдёт действительно далеко, памяти не хватит.

Есть идея почти такая же по эффективности, но которая не стоит ничего. В момент проверки очередного числа, мы можем быть уверены, что все меньшие числа сходятся к единице - они уже были проверены раньше. Достаточно условие while n > 1: заменить на while n >= n0:. Как только мы опустились ниже стартовой точки, проверку можно прервать.

Другими словами, чтобы опровергнуть гипотезу достаточно найти одно минимальное число. А для такого числа вся последовательность не меньше его самого. Поэтому более строгое условие цикла работает.

После этого изменения программа будет зацикливаться на единице, но единица не интересна, а для остальных чисел всё работает как прежде, только быстрее. Программа с этим изменением за десять секунд проверяет четырнадцать миллионов чисел. Ускорение в девять раз.

Если пропускать чётные стартовые точки (а они точно не интересны), можно успеть обработать семнадцать миллионов чисел. В одиннадцать раз быстрее:

import threading
import time

n0 = int(input())


def log():
    while True:
        time.sleep(1)
        print(f'{n0:_}')


t = threading.Thread(target=log, daemon=True)
t.start()

n0 += 1 - n0 % 2
assert n0 % 2 == 1
if n0 <= 1:
    n0 = 3
while True:
    n = n0
    while n >= n0:
        if n % 2 == 0:
            n //= 2
        else:
            n = (3 * n + 1) // 2
    n0 += 2

$ echo 1 | python search-3.py 
1_720_815
3_406_169
5_065_131
6_693_379
8_455_963
10_258_201
11_861_951
13_649_891
15_444_755
17_251_301
...

Длинные шаги

Как будет вести себя последовательность Коллатца для числа вида 16k + 7?

n = ak + b	комментарий
16k + 7 →	Для любых k это число нечётное, применяем правило (3n+1)/2.
→ (3·16/2)k + (3·7+1)/2 = 24k + 11 →	Опять нечётное, опять применяем правило (3n+1)/2.
→ (3·24/2)k + (3·11+1)/2 = 36k + 17 →	Снова нечётное.
→ (3·36/2)k + (3·17+1)/2 = 54k + 26 →	Это число чётное для любых k. Правило n/2.
→ (54/2)k + 26/2 = 27k + 13	Чётность этого числа зависит от чётности k. Для разных k требуются разные правила. Останавливаемся тут.

Любое число вида 16k + 7 через четыре шага перейдёт в число 27k + 13. Мы можем перескочить на четыре позиции вперёд не вычисляя промежуточные значения. В общем случае действуют правила:

Если a - чётное, b - нечётное, то n = ak + b - нечётное и действует правило (3n+1)/2: ak + b → (3a/2)k + (3b+1)/2.
Если a - чётное, b - чётное, то n = ak + b - чётное и действует правило n/2: ak + b → (a/2)k + b/2.
Если a - нечётное, то для разных k чётность n = ak + b меняется, останавливаемся.

Таблица для всех остатков по модулю 16:

16k +  0 ->  8k +  0 ->  4k +  0 ->  2k +  0 ->  1k +  0
16k +  1 -> 24k +  2 -> 12k +  1 -> 18k +  2 ->  9k +  1
16k +  2 ->  8k +  1 -> 12k +  2 ->  6k +  1 ->  9k +  2
16k +  3 -> 24k +  5 -> 36k +  8 -> 18k +  4 ->  9k +  2
16k +  4 ->  8k +  2 ->  4k +  1 ->  6k +  2 ->  3k +  1
16k +  5 -> 24k +  8 -> 12k +  4 ->  6k +  2 ->  3k +  1
16k +  6 ->  8k +  3 -> 12k +  5 -> 18k +  8 ->  9k +  4
16k +  7 -> 24k + 11 -> 36k + 17 -> 54k + 26 -> 27k + 13
16k +  8 ->  8k +  4 ->  4k +  2 ->  2k +  1 ->  3k +  2
16k +  9 -> 24k + 14 -> 12k +  7 -> 18k + 11 -> 27k + 17
16k + 10 ->  8k +  5 -> 12k +  8 ->  6k +  4 ->  3k +  2
16k + 11 -> 24k + 17 -> 36k + 26 -> 18k + 13 -> 27k + 20
16k + 12 ->  8k +  6 ->  4k +  3 ->  6k +  5 ->  9k +  8
16k + 13 -> 24k + 20 -> 12k + 10 ->  6k +  5 ->  9k +  8
16k + 14 ->  8k +  7 -> 12k + 11 -> 18k + 17 -> 27k + 26
16k + 15 -> 24k + 23 -> 36k + 35 -> 54k + 53 -> 81k + 80

В программе таблица для 2¹³ строится в начале, затем используется при построении последовательностей Коллатца. Вместо одного шага делается сразу тринадцать, что ускоряет проверку гипотезы. 39 миллионов чисел за десять секунд, в 26 раз быстрее:

import threading
import time


def make_fp(p):
    aps = []
    a0 = 2 ** p
    for b0 in range(2 ** p):
        a, b = a0, b0
        while a % 2 == 0:
            if b % 2 == 0:
                a, b = a // 2, b // 2
            else:
                a, b = 3 * a // 2, (3 * b + 1) // 2
        aps.append((a, b))

    mask = 2 ** p - 1

    def fp(n):
        k = n >> p
        b0 = n & mask
        a, b = aps[b0]
        return a * k + b

    return fp


p, n0 = map(int, input().split())
fp = make_fp(p)


def log():
    while True:
        time.sleep(1)
        print(f'{n0:_}')


t = threading.Thread(target=log, daemon=True)
t.start()

n0 += 1 - n0 % 2
assert n0 % 2 == 1
if n0 <= 1:
    n0 = 3
while True:
    n = n0
    while n >= n0:
        n = fp(n)
    n0 += 2

$ echo 13 1 | python search-4.py
4_003_361
7_987_545
11_902_293
15_889_673
19_897_363
23_905_145
27_819_173
31_825_151
35_779_689
39_745_659
...

Фильтрация остатков

Вернёмся к таблице остатков и посмотрим на эту строку:

16k + 3 -> 24k + 5 -> 36k + 8 -> 18k + 4 -> 9k + 2

Для любого неотрицательного k верно что 9k + 2 < 16k + 3. Напомню что мы ищем самого маленького нарушителя гипотезы. Если 16k + 3 нарушает, то и 9k + 2 нарушает, а он меньше. Следовательно, 16k + 3 не может оказаться минимальным нарушителем. Следовательно числа этого вида можно не проверять. Эта та же идея что с пропуском чётных чисел, только мы прошли немного дальше.

Можно показать что если в строке есть коэффицент при k, который меньше 16, эту строку рассматривать при поиске не нужно. Иногда надо будет проверить несколько маленьких k. Например в строке 16k + 1 → 9k + 1 возможно равенство для k = 0, его нужно проверить отдельно, перед удалением строки.

Вычистим таблицу от таких строк:

16k +  7 -> 24k + 11 -> 36k + 17 -> 54k + 26 -> 27k + 13
16k +  9 -> 24k + 14 -> 12k +  7 -> 18k + 11 -> 27k + 17
16k + 11 -> 24k + 17 -> 36k + 26 -> 18k + 13 -> 27k + 20
16k + 14 ->  8k +  7 -> 12k + 11 -> 18k + 17 -> 27k + 26
16k + 15 -> 24k + 23 -> 36k + 35 -> 54k + 53 -> 81k + 80

Смотреть можно не только на последний столбец, а на все столбцы. Тогда останутся только эти строки:

16k +  7 -> 24k + 11 -> 36k + 17 -> 54k + 26 -> 27k + 13
16k + 11 -> 24k + 17 -> 36k + 26 -> 18k + 13 -> 27k + 20
16k + 15 -> 24k + 23 -> 36k + 35 -> 54k + 53 -> 81k + 80

Из шестнадцати строк остались только три. Что обещает ускорение больше чем в пять раз. В таблице m - размер исходной таблицы остатков, r - сколько строк остаются после чистки, m/r - ожидаемое ускорение.

                ускорение                        ускорение
  m          r     m/r            m           r     m/r

                                 2^20     27328    38.4
 2^1         1     2.0           2^21     46611    45.0
 2^2         1     4.0           2^22     93222    45.0
 2^3         2     4.0           2^23    168807    49.7
 2^4         3     5.3           2^24    286581    58.5
 2^5         4     8.0           2^25    573162    58.5
 2^6         8     8.0           2^26   1037374    64.7
 2^7        13     9.8           2^27   1762293    76.2
 2^8        19    13.5           2^28   3524586    76.2
 2^9        38    13.5           2^29   6385637    84.1
2^10        64    16.0           2^30  12771274    84.1
2^11       128    16.0           2^31  23642078    90.8
2^12       226    18.1           2^32  41347483   103.9
2^13       367    22.3           2^33  82694966   103.9
2^14       734    22.3           2^34 151917636   113.1
2^15      1295    25.3           2^35 263841377   130.2
2^16      2114    31.0
2^17      4228    31.0
2^18      7495    35.0
2^19     14990    35.0

Строка m = 2¹ - эта оптимизация у нас уже есть, это отрасывание чётных кандидатов. Строка m = 2³² обещает ускорение больше чем в сто раз за счёт таблицы из 41 миллиона строк. И надо понимать что на построение таблицы тоже уходит время.

Программа строит две таблицы: одна - фильтр для остатков (функция make_aps), вторая - таблица прыжков вдоль конкретной последовательности Коллатца (функция make_fp). Они разные, потому что их размеры по разному влияют на производительность и потребление памяти.

import threading
import time


def check_n(n):
    while n > 1:
        if n % 2 == 0:
            n //= 2
        else:
            n = (3 * n + 1) // 2


def make_aps(p):
    a0 = 1
    aps = [(0, 1, 0)]

    for _ in range(p):
        prev_a0 = a0
        prev_aps = aps

        a0 = 2 * prev_a0 
        aps = []
        for r in (0, 1):
            for prev_b0, prev_a, prev_b in prev_aps:
                b0 = prev_b0 + r * prev_a0
                a = 2 * prev_a
                b = prev_b + r * prev_a
                if b % 2 == 0:
                    a //= 2
                    b //= 2
                else:
                    a = 3 * a // 2
                    b = (3 * b + 1) // 2

                assert a % 2 == 1
                assert a != a0
                if a < a0:
                    for k in range((b - b0) // (a0 - a) + 1):
                        check_n(a0 * k + b0)
                    continue

                aps.append((b0, a, b))

    return a0, aps


def make_fp(p):
    aps = []
    a0 = 2 ** p
    for b0 in range(2 ** p):
        a, b = a0, b0
        while a % 2 == 0:
            if b % 2 == 0:
                a, b = a // 2, b // 2
            else:
                a, b = 3 * a // 2, (3 * b + 1) // 2
        aps.append((a, b))

    mask = 2 ** p - 1

    def fp(n):
        k = n >> p
        b0 = n & mask
        a, b = aps[b0]
        return a * k + b

    return fp


def main():
    p_aps, p_f, n0 = map(int, input().split())

    a0, aps = make_aps(p_aps)
    print(f'aps ready {len(aps)} {a0 / len(aps):.1f}')

    fp = make_fp(p_f)
    print('fp ready')

    def log():
        while True:
            time.sleep(1)
            print(f'{n0:_}')

    t = threading.Thread(target=log, daemon=True)
    t.start()

    k = n0 // a0
    if aps[0][0] == 1 and k == 0:
        k = 1
    while True:
        base = a0 * k
        for b0, a, b in aps:
            n0 = base + b0
            n = a * k + b
            while n >= n0:
                n = fp(n)
        k += 1


main()

Фильтр по остаткам желательно делать как можно больше. Его ограничивает память. Я остановился на 2²⁹. Программа с таким фильтром занимает около одного гигабайта памяти.

Размер таблицы переходов подбирается по скорости работы. У меня лучший результат для 2¹¹. Дальнее увеличение шага эффекта не имеет.

В конфигурации 2²⁹ × 2¹¹ программа тратит одиннадцать секунд на подготовку таблиц и за следующие десять секунд обрабатывает больше одного миллиарда чисел. Время подготовки я исключаю, оно не влияет на производительность программы, которая предназначена работать сутками. Ускорение около 670 раз:

$ echo 29 11 1 | python search-5.py
aps ready 6385637 84.1
fp ready
111_491_695
222_265_087
333_651_739
443_895_143
554_737_823
661_690_015
768_067_583
872_693_503
981_254_847
1_088_366_587
...

Производительность программы мало деградирует с ростом n. Для стартовой точки 10¹⁸ производительность составит около 800 миллионов чисел за десять секунд. Падение около 20%, но нужно учесть что 10¹⁸ это далеко. Чтобы проверить все числа до этой границы на моём компьютере понадобится 341 год.

Куда двигаться дальше

Верификация кода

Надо проверить код. Дело в том что никто никогда не видел чисел нарушающих гипотезу, следовательно код не может быть протестирован полностью. Он должен быть вычитан и проверен другими людьми. Лучше всего, если разные люди напишут несколько реализаций ключевых процедур, которые можно будет сравнить в тестах. Будет обидно если в коде обнаружится ошибка, которая сведёт ваши усилия на нет. Такое уже было и именно c гипотезой Коллатца, проект BOINC потратил десятки тысяч долларов запуская программу, которая могла пропускать претендентов в нарушители гипотезы.

Всё что связано с железом

Естественно, нужно делать распределённую проверку. Надо будет добавить в программу верхнюю границу проверяемых чисел. Потом разбить все числа на блоки, проверку блоков распределить между компьютерами. Хорошо, что задача позволяет работать над разными блоками независимо.

Перенести вычисления на компилируемый язык вроде C++. Это даст прирост в десять-двадцать раз. Но код будет сложный - вам придётся следить за переполнениями целого типа. В Питоне таких трудностей нет - целые не ограничены.

Перенести вычисления на GPU. Это может дать ускорение ещё раз в пятьдесят.

Продвинутые алгоритмы

Всё что тут изложено - школьная арифметика и алгебра. Идеи взяты из раздела статьи в Википедии Collatz_conjecture#Optimizations. По гипотезе Коллатца написано много работ, которые исследуют её на гораздо более серьёзном уровне и позволяют написать ещё более производительные программы для проверки гипотезы.

Прочитайте Convergence verification of the Collatz problem. Там есть идеи, которые не требуют глубоких знаний в теории чисел.

Некоторые числа для ориентировки. Код в этом ответе за 341 год обработает все числа до 10¹⁸. Если вы его распределите на тысячу процессоров, проверка займёт четыре месяца. Известно что все числа до примерно 3·10²⁰ уже были проверены. Это в 300 раз больше. Не много, если помнить что мы сравниваем школьную арифметику помноженную на пару дней работы и знания специалистов по теории чисел помноженные на месяцы или годы усилий.

Удачи с вашим проектом!

БЛОГ НА HUSL