Найти наименьшее екатериноекатерининское число

Натуральное число, превышающее 1, назовём екатериноекатерининским, если оно делится как на число своих делителей, так и на обоих его соседей по натуральному ряду. [то есть, если число это N, а число его делителей d, то искомое число N должно делиться на d, d-1 и d+1]

Нетрудно доказать, что екатериноекатерининских чисел бесконечно много (эта задача имеет красивое решение в одну строчку, постарайтесь до него додуматься).

Но чему равно наименьшее екатериноекатерининское число?

У меня получился очень плохой и медленно работающий код на Python:

def numdiv(k): # возвращает количество делителей числа
    if k==1:
        return 1
    d = 2 # 1 и само число
    for i in range(2, int(k/2)+1): #
        if k % i == 0:
            d += 1
    return d
 
for i in range (6, 200000, 6): # Очевидно, любое екатериноекатерининское число делится на 6.
    if i%(numdiv(i))==0 and i%(numdiv(i)-1)==0 and i%(numdiv(i)+1)==0:
        print(i)

Пожалуйста, посоветуйте, как изменить этот код в лучшую сторону.


Ответы (2 шт):

Автор решения: Zhihar

более эффективно вот так определять кол-во делителей

def numdiv(k):  # возвращает количество делителей числа
    if k == 1:
        return 1
    d = 3 if int(k**0.5)**2 == k else 2
    for i in range(2, int(k**0.5) + 1):  #
        if k % i == 0:
            d += 2

    return d

мы рассматриваем делители только до корня, но каждый найденный делитель - это 2 делителя i и k / i

и полный квадрат надо обработать отдельно

кроме того не нужно 3 раза вычислять кол-во множителей для одного и того же числа:

num = 6
while True:
    count = numdiv(num)
    if num % count == 0 and num % (count - 1) == 0 and num % (count + 1) == 0:
        print(num)
        break
    num += 6

таким образом весь код будет таким:

import math

num = 6
while True:
    count = 1 if int(math.sqrt(num))**2 == num else 0
    count += sum(2 for i in range(1, int(math.sqrt(num)) + 1) if num % i == 0)

    if num % count == 0 and num % (count - 1) == 0 and num % (count + 1) == 0:
        break

    num += 6

print(num)
→ Ссылка
Автор решения: Stanislav Volodarskiy

Если задано число делителей можно отыскать все екатериноекатерининские (екекие) числа ему соответствующие.

Пусть c число делителей. Тогда

НОК(c - 1, c, c + 1) == (c - 1) * c * (c + 1)    , если c - чётное
НОК(c - 1, c, c + 1) == (c - 1) * c * (c + 1) / 2, если c - нечётное

Зададим n_max - верхняя граница екеких чисел которые мы ищем. Из уравнения (c - 1) * c * (c + 1) / 2 <= n_max отыщем ограничение на c.

Например для n_max = 10^9, c <= 1260. Для каждого значения c в диапазоне [2, 1260] построим все екекие числа.

Число делителей числа связано с его разложением на простые. Разложим n = p1^e1 * p2^e2 * ... * pk^ek, где pi - попарно различные простые, ei - натуральные. Тогда число делителей n равно (e1 + 1) * (e2 + 1) * ... * (ek + 1).

Представим c в виде произведения (не обязательно простых) множителей: c = f1 * f2 * ... * fk. Этому представлению соответствуют n = p1^(f1 - 1) * p2^(f2 - 1) * ... * pk^(fk - 1) у которых число делителей равно ровно c. Простые pi могут быть любыми не равными.

Перебрав все возможные разложения c на множители мы переберём все семейства n которые им соответствуют.

n которые нам подходят должны делится на НОК(c - 1, c, c + 1). Это дает ограничения на pi в разложении n. Например c = 72.

НОК(71, 72, 73) = 186588 = 2^3 * 3^2 * 71 * 73

Разложения 72 на множители и соответствующие представления n:

72 = 2 * 2 * 2 * 3 * 3   n = p1 * p2 * p3 * p4^2 * p5^2
72 = 2 * 2 * 2 * 9       n = p1 * p2 * p3 * p4^8
72 = 2 * 2 * 3 * 6       n = p1 * p2 * p3^2 * p4^5
72 = 2 * 2 * 18          n = p1 * p2 * p3^17
72 = 2 * 3 * 3 * 4       n = p1 * p2^2 * p3^2 * p4^3
72 = 2 * 3 * 12          n = p1 * p2^2 * p3^11
72 = 2 * 4 * 9           n = p1 * p2^3 * p3^8
72 = 2 * 6 * 6           n = p1 * p2^5 * p3^5
72 = 2 * 36              n = p1 * p2^35
72 = 3 * 3 * 8           n = p1^2 * p2^2 * p3^7
72 = 3 * 4 * 6           n = p1^2 * p2^3 * p3^5
72 = 3 * 24              n = p1^2 * p2^23
72 = 4 * 18              n = p1^3 * p2^17
72 = 6 * 12              n = p1^5 * p2^11
72 = 8 * 9               n = p1^7 * p2^8
72 = 72                  n = p1^71

Нас интересуют только представления n которые кратны 2^3 * 3^2 * 71 * 73 если правильно подставить простые числа. Таких представлений только два, а подстановок простых чисел три:

72 = 2 * 2 * 3 * 6       n = 71 * 73 * 3^2 * 2^5

72 = 2 * 3 * 3 * 4       n = 71 * 3^2 * 73^2 * 2^3
                         n = 73 * 3^2 * 71^2 * 2^3

Аналогичный анализ для c = 128 среди прочих показывает серию в которой одно простое не задано: n = 2^7 * 3 * 43 * 127 * p1. Вместо p1 может быть подставлено любое простое не равное 2, 3, 43, 127.

Серии для c = 448 впервые включают в себя несколько простых:

n = 2^6 * 3 * 7 * 149 * 449 * p1 * p2
n = 2^6 * 3 * 7 * 149 * 449 * p1^3
n = 2^6 * 3 * 7 * 149^3 * 449 * p1
n = 2^6 * 3 * 7^3 * 149 * 449 * p1
n = 2^6 * 3^3 * 7 * 149 * 449 * p1
n = 2^13 * 3 * 7 * 149 * 449 * p1

На этих идеях построена программа, которая порождает все екекие числа в заданном диапазоне:

$ python generate.py 100000000
27241848    (72) = 2^3 * 3^2 * 71 * 73^2
26495496    (72) = 2^3 * 3^2 * 71^2 * 73
1492704    (72) = 2^5 * 3^2 * 71 * 73
884640    (96) = 2^5 * 3 * 5 * 19 * 97
5619264    (112) = 2^6 * 3 * 7 * 37 * 113
10485120... (128) = 2^7 * 3 * 43 * 127 * p1
18873216    (128) = 2^7 * 3^3 * 43 * 127
37931712    (168) = 2^6 * 3 * 7 * 13^2 * 167
95550336    (288) = 2^7 * 3^2 * 7 * 17^2 * 41
71662752    (288) = 2^5 * 3^3 * 7 * 17^2 * 41

Параметр - верхняя граница для поиска екеких чисел. Программа печатает найденные числа в таблице:

<екекое число> (<число его делителей>) = <разложение>

Для серий выводится наименьший представитель из серии и печатается многоточие.

Поиск всех чисел до 10^12 требует три секунды, возвращает 207 чисел и 64 серии.

Поиск всех чисел до 10^18 требует тринадцать с половиной часов, возвращает 4210 чисел и 3250 серий. Собственно все числа программа отыскивает за пятнадцать секунд, остальные часы уходят на проверку что других чисел нет (то есть разложения есть, но они порождают слишком большие произведения).

Код:

import math
import operator
import sys


def binary_search(lo, hi, p):
    while lo < hi - 1:
        m = (lo + hi) // 2
        if p(m):
            assert m < hi
            hi = m
        else:
            assert lo < m
            lo = m
    return hi


def max_c_by_n(n):
    return binary_search(1, n, lambda c: (c - 1) * c * (c + 1) // 2 > n)


def factors(n):
    i = 2
    j = n
    j_sqrt = math.isqrt(j)
    while i <= j_sqrt:
        if j % i == 0:
            p = 0
            while j % i == 0:
                j //= i
                p += 1
            yield i, p
            j_sqrt = math.isqrt(j)
        i += 1 if i == 2 else 2
    if j > 1:
        yield j, 1


def divisors(n):
    f = tuple(factors(n))

    def search(n):
        if n == len(f):
            yield 1
            return
        for d in search(n + 1):
            c = 1
            for i in range(f[n][1] + 1):
                yield d * c
                c *= f[n][0]

    return search(0)


def representations(c, n_max):
    lcm = math.lcm(c - 1, c, c + 1)
    f = sorted(factors(lcm), key=operator.itemgetter(-1), reverse=True)

    def search(n, p, i, tail):
        if p > n_max:
            return
        if i == len(f):
            yield n, tail
            return
        for d in divisors(n):
            if d > f[i][1]:
                b = f[i][0]
                e = d - 1
                yield from search(n // d, p * b ** e, i + 1, (b, e, tail))

    def retrieve(t):
        n, tail = t

        def traverse(t):
            while t is not None:
                b, e, t = t
                yield b, e

        return sorted(traverse(tail)), n

    return map(retrieve, search(c, 1, 0, None))


def is_prime(n):
    return all(n % d != 0 for d in range(2, math.isqrt(n) + 1))


def next_prime(n):
    n += 1
    while not is_prime(n):
        n += 1
    return n


def series(n, fixed_factor, n_max, excluded_primes):

    def search(n, max_d, factor, last_n, tail):
        if n == 1:
            yield factor, tail
            return

        p = next_prime(last_n)
        while p in excluded_primes:
            p = next_prime(p)

        for d in divisors(n):
            if 2 <= d <= max_d:
                e = d - 1
                next_factor = factor * p ** e
                if next_factor <= n_max:
                    yield from search(n // d, d, next_factor, p, (e, tail))

    def retrieve(t):

        def traverse(t):
            while t is not None:
                e, t = t
                yield e

        return tuple(traverse(t))[::-1]

    return ((f, retrieve(t)) for f, t in search(n, n, fixed_factor, 1, None))


def show_series(s):
    return ' * '.join(
        f'p{i}' if e == 1 else f'p{i}^{e}' for i, e in enumerate(s, start=1)
    )


def main():
    n = int(sys.argv[1])
    for c in range(2, max_c_by_n(n)):
        for f, e in representations(c, n):
            fixed_factor = math.prod(b ** e for b, e in f)
            fixed_expr = ' * '.join(str(b) if e == 1 else f'{b}^{e}' for b, e in f)
            if e == 1:
                print(f'{fixed_factor}    ({c}) = {fixed_expr}')
            else:
                for f, s in series(e, fixed_factor, n, set(b for b, _ in f)):
                    series_part = show_series(s)
                    print(f'{f}... ({c}) = {fixed_expr}', '*', series_part)


main()
→ Ссылка