Поиск чисел с пятью нечетными делителями

"""
Найдите все натуральные числа, принадлежащие отрезку [35000000;40000000], у которых ровно пять различных нечётных делителей (количество чётных делителей может быть любым). В ответе перечислите найденные числа в порядке возрастания.

Для решения задачи воспользуемся тем, что число имеет ровно 5 различных нечётных делителей, только если оно представимо в виде 2^k * p^4, где k - степень двойки, p - простое число.
"""
import math
import time
from typing import Any, Callable, Generator
import typing


class Number:
    """
    Число n вместе с кешированными значениями n^2 и n^4

    - n (number) - само число
    - _n2 (number) - n^2, используется для ускорения проверки на простоту
    - _n4 (number) - n^4, используется для ускорения вычислений
    - _log2_n4 (number) - логарифм n^4 по основанию 2, используется для ускорения вычислений
    """

    def __init__(self, number: int):
        self.n = number
        self._n2 = number * number
        self._n4 = self._n2 * self._n2

        # В решателе используются логарифмы, поэтому кешируем их
        self._log2_n4 = math.log2(self._n4)

    def __repr__(self):
        return f"Number({self.n})"


class Solver:
    """
    Решатель задачи. Находит все числа вида 2^k * p^4, где k - степень двойки, p - простое число в заданном диапазоне.


    - left (int) - левая граница диапазона (включительно)
    - right (int) - правая граница диапазона (включительно)
    """

    PRIMES_100 = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37,
                  41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

    def __init__(self, left: int, right: int):
        self.left = min(left, right)
        self.right = max(left, right)
        if self.left < 2:
            raise ValueError("Left bound must be at least 2.")
        # В решателе используются логарифмы, поэтому кешируем их
        self._log2_left = math.log2(self.left)
        self._log2_right = math.log2(self.right)

        # Список простых чисел. Нужно обеспечить, что в нём есть все простые числа до \sqrt[4](right)
        # начинаем с простых чисел до 100
        self.primes = list(map(Number, Solver.PRIMES_100))
        # Добавляем простые числа до \sqrt[4](right)
        while self.primes[-1]._n4 < self.right:
            self.add_next_prime()

    def add_next_prime(self) -> Number:
        """
        Находит следующее простое число и добавляет его в список простых чисел.
        """
        p = self.primes[-1].n + 2
        while not self.is_prime(p):
            p += 2
        self.primes.append(Number(p))
        return self.primes[-1]

    def is_prime(self, n: int) -> bool:
        """
        Проверяет, является ли число простым.

        Более точно, проверяет, что среди self.primes нет делителей числа n.
        """
        assert (n >= 2)
        for prime in self.primes:
            if prime._n2 > n:
                break
            if n % prime.n == 0:
                return False
        return True

    class Result:
        """
        Результат вычислений.

        - prime (Number) - простое число.
        - degree (int) - степень двойки
        - value (int) - число вида 2^degree * prime^4
        """

        def __init__(self, prime: Number, degree: int):
            self.prime = prime
            self.degree = degree
            self.value = (1 << degree)*self.prime._n4

        def __repr__(self):
            return f"Result({self.prime}, {self.degree}, value={self.value})"

    def is_proper_num(self, p: Number) -> list[Result]:
        """
        Проверяет, является ли число p подходящим под условия задачи.

        Если число не подходит, возвращает пустой список.
        Если число подходит, возвращает список Result(p, k), где k - степень двойки, при которой 2^k * p^4 попадает в заданный диапазон.
        """
        if p._n4 > self.right:
            return []
        if p._n4 > self.left:
            return [Solver.Result(p, 0)]
        # Минимальная степень двойки, при которой 2^k * p^4 >= left
        deg = math.ceil(self._log2_left - p._log2_n4)
        # Максимальная степень двойки, при которой 2^k * p^4 <= right
        lim = math.floor(self._log2_right - p._log2_n4)

        if deg > lim:
            # В заданном диапазоне нет чисел вида 2^k * p^4
            return []
        return [Solver.Result(p, d) for d in range(deg, lim + 1)]

    def gen_proper_odd_primes(self) -> Generator[Result, None, None]:
        """
        Генератор, который перечисляет Result(p,k), то есть числа вида 2^k * p^4, где k - степень двойки, p - простое число, и 2^k * P^4 находится в заданном диапазоне.
        """
        for p in self.primes[1:]:
            yield from self.is_proper_num(p)
            # Если p^4 > right, то дальше нет смысла проверять
            if p._n4 > self.right:
                break

    def solve(self) -> list[Result]:
        """
        Решает задачу, находя все числа вида 2^k * p^4, где k - степень двойки, p - простое число в заданном диапазоне.

        Возвращает список объектов Result, соответствующих найденным простым числам.
        """
        return list(self.gen_proper_odd_primes())


if __name__ == "__main__":
    left = 35_000_000_000
    right = 40_000_000_000

if __name__ == "__main__":
    start_time = time.time()
    solver = Solver(left, right)
    results = solver.solve()
    end_time = time.time()
    print(f"Time taken: {end_time - start_time:.6f} seconds")

    for result in results:
        print(result)
    numbers = sorted([r.value for r in results])
    for n in numbers:
        print(n)
    print(f"Found {len(numbers)} results")

# Сравнение результатов


def fox(left, right):
    # Функция: проверка, является ли число простым (достаточно до sqrt)
    def is_prime(v_n):
        if v_n < 2:
            return False
        for v_i in range(2, int(v_n ** 0.5) + 1):
            if v_n % v_i == 0:
                return False
        return True

    # Список для хранения результатов
    lst_results = []

    # Перебираем все нечётные простые числа p, такие что p^4 * 2^k попадает в диапазон
    for v_p in range(3, math.ceil(math.pow(right, 0.25)), 2):  # 1000 — с запасом, т.к. 999^4 > 1e12
        if is_prime(v_p):
            v_p4 = v_p ** 4  # Вычисляем p^4

            v_n = v_p4
            while v_n <= right:
                if v_n >= left:
                    # Получаем и сортируем нечётные делители (всегда 5 штук)
                    lst_divs = [1, v_p, v_p**2, v_p**3, v_p**4]
                    lst_results.append((v_n, lst_divs))
                v_n *= 2  # Умножаем на степень двойки
    return lst_results


if __name__ == "__main__":
    fox_start_time = time.time()
    fox_results = fox(left, right)
    fox_end_time = time.time()
    print(f"Fox time taken: {fox_end_time - fox_start_time:.6f} seconds")
    for n, divs in sorted(fox_results, key=lambda x: x[0]):
        print(n, divs)
    print(f"Fox found {len(fox_results)} results")


def volodarski(left, right):

    def is_prime(n):
        return all(n % i != 0 for i in range(2, n))

    def numbers(m, n):
        i = 3
        while True:
            if is_prime(i):
                j = i ** 4
                if n < j:
                    break
                while j <= n:
                    if m <= j:
                        yield j
                    j *= 2
            i += 2

    return list(numbers(left, right))


def volodarski2(left, right):
    # из решения ru.stackoverflow.com/questions/1276027
    def odd_primes(n):
        sieve = bytearray(n // 2)
        if sieve:
            sieve[0] = 1
            for p in range(3, math.isqrt(n - 1) + 1, 2):
                for i in range(p * p // 2, n // 2, p):
                    sieve[i] = 1
        return (2 * i + 1 for i, k in enumerate(sieve) if k == 0)

    def numbers(m, n):
        e = 1 << m.bit_length()
        for p in odd_primes(math.isqrt(math.isqrt(n)) + 1):
            j = e * p ** 4
            # assert j >= m
            while j >= 2 * m:
                j //= 2
                e //= 2
            # assert m <= j <= 2 * m
            while j <= n:
                yield j
                j *= 2
    return sorted(numbers(left, right))


if __name__ == "__main__":
    volodarski_start_time = time.time()
    volodarski_results = volodarski2(left, right)
    volodarski_end_time = time.time()
    print(
        f"Volodarski time taken: {volodarski_end_time - volodarski_start_time:.6f} seconds")
    for n in sorted(volodarski_results):
        print(n)
    print(f"Volodarski found {len(volodarski_results)} results")
    # Сравнение результатов

# benchmark


def dt_n(f, n) -> float:
    """
    Функция для измерения времени выполнения функции f n раз.
    """
    t0 = time.time()
    for _ in range(n):
        f()
    return (time.time() - t0)


def guess_n(f: Callable, T: float) -> tuple[int, float, int]:
    """
    Функция для оценки количества итераций n, необходимых для достижения времени T.

    Возвращает кортеж из трех значений:
    - N (int) - оценка количества итераций,
    - dt (float) - время выполнения функции f во время подбора числа итераций,
    - n (int) - количество итераций при последнем измерении времени.
    """
    # Начинаем с n=1 и увеличиваем его, пока время выполнения не станет больше некоторого порога
    # (например, 100 миллисекунд)
    n = 1
    dt = dt_n(f, n)

    while dt < 1e-4:
        n *= 2
        dt = dt_n(f, n)

    N = math.floor((T/dt) * n)
    return max(1, N), dt, n


def benchmark(f: Callable, T: float = 1.0) -> tuple[float, int]:
    """
    Функция для бенчмарка функции f.
    """
    N, dT, n = guess_n(f, T)
    if dT < T/2:
        # Подбор числа итераций занял немного времени, поэтому
        # повторяем измерение времени с большим количеством итераций
        dT = dt_n(f, N)
        # Иногда время оказывается отрицательным, поэтому
        # повторяем измерение времени, пока оно не станет положительным
        while dT < 0:
            dT = dt_n(f, N)
    else:
        # Подбор числа итераций занял много времени, поэтому
        # используем его для оценки времени выполнения функции f
        N = n
    return dT, N


def benchmark_msg(prefix: str, f: Callable, T: float = 1):
    """
    Функция для бенчмарка функции f с выводом сообщения.
    """
    dT, N = benchmark(f, T)
    print(f"{prefix:16s}: {dT:.6f} seconds for {N} iterations", end="\t")
    dt = dT/N
    # if dt < 1e-4:
    #     print(f"average: {dt*1000:.6f}ms")
    # elif dt < 1:
    #     print(f"average: {dt:.6f}s")
    # else:
    #     print(f"average: {dt/60:.6f} min")
    print(f"average: {dt*1000:.6f}ms")
    return dT, N


if __name__ == "__main__":
    left_0 = 35_000_000
    right_0 = 40_000_000
    for i in range(0, 20, 3):
        left = left_0 * (10**i)
        right = right_0 * (10**i)

        print(f"# Testing with left={left:_}, right={right:_}")
        benchmark_msg("Solver", lambda: Solver(left, right).solve())
        benchmark_msg("Fox", lambda: fox(left, right))
        benchmark_msg("Volodarski2", lambda: volodarski2(left, right))