Ξένη Γήινος's questions

Presumo que todos vocês saibam o que são números primos e o que é o Crivo de Eratóstenes, então não vou perder tempo explicando-os.

Agora, todos os números primos, exceto 2, são números ímpares, então precisamos verificar apenas os números ímpares. Isso é muito óbvio, mas vale a pena mencionar porque essa otimização simples reduziu pela metade os candidatos que precisamos verificar, e então precisamos marcar apenas os múltiplos ímpares de números primos diferentes de 2.

Outra otimização simples é verificar apenas números até a raiz quadrada do limite, e outra otimização simples é marcar como números compostos o início no quadrado do número primo. Como os motivos para isso devem ser óbvios, não explicarei o porquê, embora eu não consiga pensar em outras otimizações relacionadas à marcação de números compostos.

Mas podemos restringir ainda mais o espaço de busca: todos os números primos, exceto 2, 3 e 5, devem ser congruentes a [1, 7, 11, 13, 17, 19, 23, 29]30%. Isso é evidente pela natureza do módulo. Existem apenas 30 possibilidades; se o módulo for par, o número é múltiplo de 2, e as outras possibilidades significam que o número é múltiplo de 3, múltiplo de 5 ou ambos. Em outras palavras, todos os números primos devem ser coprimos a 30, exceto 2, 3 e 5.

Em Python, isso é:

wheel3 = [i for i in range(1, 30) if math.gcd(i, 30) == 1]
# [1, 7, 11, 13, 17, 19, 23, 29]

Agora calculamos a diferença entre pares consecutivos de elementos, começando em 7, e 1 vem imediatamente depois de 29 devido à natureza da operação de módulo, por exemplo, 31 % 30 == 1e então a diferença entre eles é 2.

Obtemos o seguinte: [4, 2, 4, 2, 4, 6, 2, 6].

Agora, de cada 30 números, precisamos verificar apenas 8, pulando 22. Isso representa uma melhoria significativa em relação à otimização anterior, que aplicava força bruta apenas aos números ímpares. Precisávamos processar 15 números de cada 30 números se processássemos todos os números ímpares. Agora, temos 7 números a menos, e o espaço de busca é reduzido para 4/15, que é 0,2666...

Podemos otimizar isso ainda mais usando uma roda maior, usando 2 * 3 * 5 * 7 = 210 como base, todos os números primos começando em 11 devem ser coprimos com 210.

wheel4 = [i for i in range(1, 210) if math.gcd(i, 210) == 1]
wheel4 == [
    1  , 11 , 13 , 17 , 19 , 23 ,
    29 , 31 , 37 , 41 , 43 , 47 ,
    53 , 59 , 61 , 67 , 71 , 73 ,
    79 , 83 , 89 , 97 , 101, 103,
    107, 109, 113, 121, 127, 131,
    137, 139, 143, 149, 151, 157,
    163, 167, 169, 173, 179, 181,
    187, 191, 193, 197, 199, 209
]

E a lista de alterações de índice é:

[
    2 , 4 , 2 , 4 , 6 , 2 ,
    6 , 4 , 2 , 4 , 6 , 6 ,
    2 , 6 , 4 , 2 , 6 , 4 ,
    6 , 8 , 4 , 2 , 4 , 2 ,
    4 , 8 , 6 , 4 , 6 , 2 ,
    4 , 6 , 2 , 6 , 6 , 4 ,
    2 , 4 , 6 , 2 , 6 , 4 ,
    2 , 4 , 2 , 10, 2 , 10
]

Agora, de cada 210 números, precisamos processar apenas 48 números, em comparação com os 56 números anteriores, precisamos processar 8 números a menos, reduzindo o espaço de busca para 8/35, que é 0,22857142857142857... e menos de um quarto.

Então, espero que a versão que usa a roda baseada em 210 leve apenas 6/7 ou 85,71% do tempo que a versão baseada em 30 leva para ser executada, mas não é isso que acontece.

import math
import numpy as np
import numba as nb


@nb.njit(cache=True)
def prime_wheel_sieve(n: int) -> np.ndarray:
    wheel = [4, 2, 4, 2, 4, 6, 2, 6]
    primes = np.ones(n + 1, dtype=np.uint8)
    primes[:2] = False
    for square, step in ((4, 2), (9, 6), (25, 10)):
        primes[square::step] = False

    k = 7
    lim = int(math.sqrt(n) + 0.5)
    i = 0
    while k <= lim:
        if primes[k]:
            primes[k**2 :: 2 * k] = False

        k += wheel[i]
        i = (i + 1) & 7

    return np.nonzero(primes)[0]


# fmt: off
WHEEL4 = np.array([
    2 , 4 , 2 , 4 , 6 , 2 ,
    6 , 4 , 2 , 4 , 6 , 6 ,
    2 , 6 , 4 , 2 , 6 , 4 ,
    6 , 8 , 4 , 2 , 4 , 2 ,
    4 , 8 , 6 , 4 , 6 , 2 ,
    4 , 6 , 2 , 6 , 6 , 4 ,
    2 , 4 , 6 , 2 , 6 , 4 ,
    2 , 4 , 2 , 10, 2 , 10
], dtype=np.uint8)
# fmt: on


@nb.njit(cache=True)
def prime_wheel_sieve4(n: int) -> np.ndarray:
    primes = np.ones(n + 1, dtype=np.uint8)
    primes[:2] = False
    for square, step in ((4, 2), (9, 6), (25, 10), (49, 14)):
        primes[square::step] = False

    k = 11
    lim = int(math.sqrt(n) + 0.5)
    i = 0
    while k <= lim:
        if primes[k]:
            primes[k**2 :: 2 * k] = False

        k += WHEEL4[i]
        i = (i + 1) % 48

    return np.nonzero(primes)[0]

In [549]: np.array_equal(prime_wheel_sieve(65536), prime_wheel_sieve4(65536))
Out[549]: True

In [550]: %timeit prime_wheel_sieve(65536)
161 μs ± 1.13 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [551]: %timeit prime_wheel_sieve4(65536)
163 μs ± 1.79 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [552]: %timeit prime_wheel_sieve4(131072)
330 μs ± 10.6 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [553]: %timeit prime_wheel_sieve(131072)
328 μs ± 7.4 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [554]: %timeit prime_wheel_sieve4(262144)
680 μs ± 14.3 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [555]: %timeit prime_wheel_sieve(262144)
669 μs ± 7.79 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [556]: %timeit prime_wheel_sieve(524288)
1.44 ms ± 16.2 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [557]: %timeit prime_wheel_sieve4(524288)
1.48 ms ± 13.4 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [558]: %timeit prime_wheel_sieve4(1048576)
3.25 ms ± 81.3 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [559]: %timeit prime_wheel_sieve(1048576)
3.23 ms ± 115 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [560]: %timeit prime_wheel_sieve(2097152)
7.08 ms ± 80.9 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [561]: %timeit prime_wheel_sieve4(2097152)
7.1 ms ± 85.9 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [562]: %timeit prime_wheel_sieve4(4194304)
14.8 ms ± 120 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [563]: %timeit prime_wheel_sieve(4194304)
14.2 ms ± 145 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [564]: %timeit prime_wheel_sieve(8388608)
39.4 ms ± 1.44 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [565]: %timeit prime_wheel_sieve4(8388608)
41.7 ms ± 2.56 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

De acordo com meus testes, usar uma roda maior torna o código mais lento, e não mais rápido. Por que isso acontece? Teoricamente, usar uma roda maior estreita o espaço de busca, então não deveria aumentar o tempo de execução. Por que usar uma roda maior torna o código mais lento?

Ok, usando o método científico, controlei as diferenças entre as duas funções de modo que a única diferença entre elas, a única quantidade que pode afetar o desempenho, é a roda usada.

Fiz a primeira função usar o operador módulo em vez do AND bit a bit, e fiz a segunda função usar uma lista local, assim como a primeira. Eu queria separar código e dados, mas tudo bem.

@nb.njit(cache=True)
def prime_wheel_sieve3(n: int) -> np.ndarray:
    wheel = [4, 2, 4, 2, 4, 6, 2, 6]
    primes = np.ones(n + 1, dtype=np.uint8)
    primes[:2] = False
    for square, step in ((4, 2), (9, 6), (25, 10)):
        primes[square::step] = False

    k = 7
    lim = int(math.sqrt(n) + 0.5)
    i = 0
    while k <= lim:
        if primes[k]:
            primes[k**2 :: 2 * k] = False

        k += wheel[i]
        i = (i + 1) % 8

    return np.nonzero(primes)[0]


@nb.njit(cache=True)
def prime_wheel_sieve4_1(n: int) -> np.ndarray:
    # fmt: off
    wheel = [
        2 , 4 , 2 , 4 , 6 , 2 ,
        6 , 4 , 2 , 4 , 6 , 6 ,
        2 , 6 , 4 , 2 , 6 , 4 ,
        6 , 8 , 4 , 2 , 4 , 2 ,
        4 , 8 , 6 , 4 , 6 , 2 ,
        4 , 6 , 2 , 6 , 6 , 4 ,
        2 , 4 , 6 , 2 , 6 , 4 ,
        2 , 4 , 2 , 10, 2 , 10
    ]
    # fmt: on
    primes = np.ones(n + 1, dtype=np.uint8)
    primes[:2] = False
    for square, step in ((4, 2), (9, 6), (25, 10), (49, 14)):
        primes[square::step] = False

    k = 11
    lim = int(math.sqrt(n) + 0.5)
    i = 0
    while k <= lim:
        if primes[k]:
            primes[k**2 :: 2 * k] = False

        k += wheel[i]
        i = (i + 1) % 48

    return np.nonzero(primes)[0]

Tive que adicionar comentários de formatação para evitar que o Black Formatter bagunçasse minha tabela formatada no VS Code e, claro, isso não afeta o desempenho em nada.

As únicas diferenças entre as funções são o valor inicial de k, os primos que tiveram que ser processados ​​antes de girar a roda propriamente dita e, claro, o tamanho da roda (e a própria roda). Estes tiveram que ser diferentes (porque usam rodas diferentes), mas todo o resto permanece inalterado.

In [679]: np.array_equal(prime_wheel_sieve3(65536), prime_wheel_sieve4_1(65536))
Out[679]: True

In [680]: %timeit prime_wheel_sieve3(65536)
162 μs ± 2.27 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [681]: %timeit prime_wheel_sieve4_1(65536)
158 μs ± 1.83 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [682]: %timeit prime_wheel_sieve3(131072)
326 μs ± 7.91 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [683]: %timeit prime_wheel_sieve4_1(131072)
322 μs ± 8.89 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [684]: %timeit prime_wheel_sieve3(262144)
659 μs ± 7.74 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [685]: %timeit prime_wheel_sieve4_1(262144)
655 μs ± 12.2 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [686]: %timeit prime_wheel_sieve3(524288)
1.45 ms ± 14.1 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [687]: %timeit prime_wheel_sieve4_1(524288)
1.42 ms ± 8.13 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [688]: %timeit prime_wheel_sieve3(1048576)
3.2 ms ± 68.4 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [689]: %timeit prime_wheel_sieve4_1(1048576)
3.2 ms ± 199 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Curiosamente, prime_wheel_sieve4_1o desempenho é um pouco mais rápido que prime_wheel_sieve3, mas é só um pouquinho. A aceleração é insignificante, mas sei que o tempo de execução do código é estocástico por natureza, então prime_wheel_sieve4_1um desempenho consistentemente mais rápido que , prime_wheel_sieve3um pouquinho, é estatisticamente significativo. Embora eu não tenha testado muito, isso não exclui a possibilidade de coincidência.

Mas, de acordo com minha teoria, eu deveria ver uma redução de 14,29% no tempo de execução, e não basicamente nenhuma melhoria. Meus testes fortaleceram meu argumento.

Então por que usar uma roda maior não acelera o código?

Tenho visto muitas postagens sobre esse assunto, mas nenhuma é exatamente o que estou procurando.

Quero encontrar todas as maneiras pelas quais um número inteiro positivo N maior que 1 pode ser expresso como a soma de no máximo N números inteiros de 1 a N, como de costume.

Por exemplo, na notação padrão, estas são todas as partições de 6:

[(1, 1, 1, 1, 1, 1),
 (1, 1, 1, 1, 2),
 (1, 1, 1, 2, 1),
 (1, 1, 1, 3),
 (1, 1, 2, 1, 1),
 (1, 1, 2, 2),
 (1, 1, 3, 1),
 (1, 1, 4),
 (1, 2, 1, 1, 1),
 (1, 2, 1, 2),
 (1, 2, 2, 1),
 (1, 2, 3),
 (1, 3, 1, 1),
 (1, 3, 2),
 (1, 4, 1),
 (1, 5),
 (2, 1, 1, 1, 1),
 (2, 1, 1, 2),
 (2, 1, 2, 1),
 (2, 1, 3),
 (2, 2, 1, 1),
 (2, 2, 2),
 (2, 3, 1),
 (2, 4),
 (3, 1, 1, 1),
 (3, 1, 2),
 (3, 2, 1),
 (3, 3),
 (4, 1, 1),
 (4, 2),
 (5, 1),
 (6,)]

Agora, a notação é de entropia muito baixa: primeiro, cada ocorrência do número aumenta o tamanho de uma partição específica. Isso é ineficiente e é difícil contar as ocorrências dos números quando eles se repetem muitas vezes. Quero substituir todas as ocorrências de um número por uma tupla de dois elementos, na qual o primeiro elemento é o número e o segundo é a contagem. Por exemplo, (1, 1, 1, 1, 1, 1)é equivalente a (1, 6). Ambos contêm a mesma informação, mas um é claramente muito mais conciso.

E segundo, há muitas duplicatas na saída; por exemplo, há cinco partições que contêm quatro 1s e um 2, que são contadas como cinco elementos separados. Isso também é ineficiente, já que a adição é comutativa; mudar a ordem dos números não altera o resultado, então todos são equivalentes, são todos o mesmo elemento.

Entretanto, se substituirmos todos os cinco por apenas um elemento, perdemos informações.

Em vez disso, desejo substituí-lo pelo seguinte formato:

Counter({((1, 2), (2, 2)): 6,
         ((1, 1), (2, 1), (3, 1)): 6,
         ((1, 4), (2, 1)): 5,
         ((1, 3), (3, 1)): 4,
         ((1, 2), (4, 1)): 3,
         ((1, 1), (5, 1)): 2,
         ((2, 1), (4, 1)): 2,
         ((1, 6),): 1,
         ((2, 3),): 1,
         ((3, 2),): 1,
         ((6, 1),): 1})

Então, quero que o resultado seja um Counterem que as chaves sejam as partições exclusivas e os valores sejam quantas maneiras os números podem ser organizados.

E sim, eu escrevi uma função para isso, usando força bruta e memorização. Ela se mostrou bem eficiente.

Primeiro, esta é a implementação que gera o resultado no formato padrão. Posto aqui para comparação:

def partitions(number: int) -> list[tuple[int, ...]]:
    result = []
    stack = [(number, ())]

    while stack:
        remaining, path = stack.pop()
        if not remaining:
            result.append(path)
        else:
            stack.extend((remaining - i, path + (i,)) for i in range(remaining, 0, -1))

    return result

São necessários 582 milissegundos para encontrar todas as partições de 20 no CPython e 200 milissegundos no PyPy3:

CPython

In [22]: %timeit partitions(20)
582 ms ± 4.22 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

PyPy3

In [36]: %timeit partitions(20)
199 ms ± 3.17 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Agora, a abordagem de força bruta com memorização que gera o formato pretendido:

PARTITION_COUNTERS = {}


def partition_counter(number: int) -> Counter:
    if result := PARTITION_COUNTERS.get(number):
        return result
    
    result = Counter()
    for i in range(1, number):
        for run, count in partition_counter(number - i).items():
            new_run = []
            added = False
            for a, b in run:
                if a == i:
                    new_run.append((a, b + 1))
                    added = True
                else:
                    new_run.append((a, b))
            
            if not added:
                new_run.append((i, 1))
            
            result[tuple(sorted(new_run))] += count
    
    result[((number, 1),)] = 1
    PARTITION_COUNTERS[number] = result
    return result

CPython

In [23]: %timeit PARTITION_COUNTERS.clear(); partition_counter(20)
10.4 ms ± 72.1 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

PyPy3

In [37]: %timeit PARTITION_COUNTERS.clear(); partition_counter(20)
9.75 ms ± 58.3 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Leva apenas 10 milissegundos para encontrar todas as partições de 20, muito, muito mais rápido que a primeira função, e o PyPy3 não o torna mais rápido.

Mas como podemos fazer melhor? Afinal, estou apenas usando força bruta. Sei que existem muitos algoritmos inteligentes para partição de números inteiros, mas nenhum deles gera saídas no formato pretendido.

Eu escrevi quatro funções que modificam uma matriz quadrada 2D no local, ela reflete metade da matriz quadrada delimitada por dois lados que se encontram e a diagonal correspondente de 45 graus, para a outra metade separada pela mesma diagonal.

Eu escrevi uma função para cada um dos quatro casos possíveis, para product(('upper', 'lower'), ('left', 'right'))refletir product(('lower', 'upper'), ('right', 'left')).

Eles usam o Numba para compilar Just-In-Time e são paralelizados usando numba.prangee, portanto, são muito mais rápidos que os métodos fornecidos pelo NumPy:

In [2]: sqr = np.random.randint(0, 256, (1000, 1000), dtype=np.uint8)

In [3]: %timeit x, y = np.tril_indices(1000); sqr[x, y] = sqr[y, x]
9.16 ms ± 30.9 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Como você pode ver, o código acima leva muito tempo para ser executado.

import numpy as np
import numba as nb


@nb.njit(cache=True, parallel=True, nogil=True)
def triangle_flip_LL2UR(arr: np.ndarray) -> None:
    height, width = arr.shape[:2]
    if height != width:
        raise ValueError("argument arr must be a square")

    for i in nb.prange(height):
        arr[i, i:] = arr[i:, i]


@nb.njit(cache=True, parallel=True, nogil=True)
def triangle_flip_UR2LL(arr: np.ndarray) -> None:
    height, width = arr.shape[:2]
    if height != width:
        raise ValueError("argument arr must be a square")

    for i in nb.prange(height):
        arr[i:, i] = arr[i, i:]


@nb.njit(cache=True, parallel=True, nogil=True)
def triangle_flip_LR2UL(arr: np.ndarray) -> None:
    height, width = arr.shape[:2]
    if height != width:
        raise ValueError("argument arr must be a square")

    last = height - 1
    for i in nb.prange(height):
        arr[i, last - i :: -1] = arr[i:, last - i]


@nb.njit(cache=True, parallel=True, nogil=True)
def triangle_flip_UL2LR(arr: np.ndarray) -> None:
    height, width = arr.shape[:2]
    if height != width:
        raise ValueError("argument arr must be a square")

    last = height - 1
    for i in nb.prange(height):
        arr[i:, last - i] = arr[i, last - i :: -1]

In [4]: triangle_flip_LL2UR(sqr)

In [5]: triangle_flip_UR2LL(sqr)

In [6]: triangle_flip_LR2UL(sqr)

In [7]: triangle_flip_UL2LR(sqr)

In [8]: %timeit triangle_flip_LL2UR(sqr)
194 μs ± 634 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [9]: %timeit triangle_flip_UR2LL(sqr)
488 μs ± 3.26 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [10]: %timeit triangle_flip_LR2UL(sqr)
196 μs ± 501 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [11]: %timeit triangle_flip_UL2LR(sqr)
486 μs ± 855 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

Por que eles têm tempos de execução tão diferentes? Dois deles levam cerca de 200 microssegundos para serem executados, os outros dois, cerca de 500 microssegundos, apesar de serem quase idênticos.

Descobri uma coisa. triangle_flip_UR2LL(arr)é o mesmo que triangle_flip_LL2UR(sqr.T)e vice-versa.

Agora, se eu transpor o array antes de chamar as funções, a tendência de desempenho se inverte:

In [109]: %timeit triangle_flip_UR2LL(sqr.T)
196 μs ± 1.15 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [110]: %timeit triangle_flip_LL2UR(sqr.T)
490 μs ± 1.24 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

Por que isso está acontecendo?

Então eu tenho esta função:

import numpy as np
import numba as nb


@nb.njit(cache=True, parallel=True, nogil=True)
def triangle_half_UR_LL(size: int, swap: bool = False) -> tuple[np.ndarray, np.ndarray]:
    total = (size + 1) * size // 2
    x_coords = np.full(total, 0, dtype=np.uint16)
    y_coords = np.full(total, 0, dtype=np.uint16)
    offset = 0
    side = np.arange(size, dtype=np.uint16)
    for i in nb.prange(size):
        offset = i * size - (i - 1) * i // 2
        end = offset + size - i
        x_coords[offset:end] = i
        y_coords[offset:end] = side[i:]
    
    return (x_coords, y_coords) if not swap else (y_coords, x_coords)

O que ele faz não é importante, o importante é que ele é compilado com JIT e Numba e, portanto, muito rápido:

In [2]: triangle_half_UR_LL(10)
Out[2]:
(array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2,
        2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5,
        5, 6, 6, 6, 6, 7, 7, 7, 8, 8, 9], dtype=uint16),
 array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 3, 4, 5, 6, 7, 8, 9, 2, 3, 4,
        5, 6, 7, 8, 9, 3, 4, 5, 6, 7, 8, 9, 4, 5, 6, 7, 8, 9, 5, 6, 7, 8,
        9, 6, 7, 8, 9, 7, 8, 9, 8, 9, 9], dtype=uint16))

In [3]: %timeit triangle_half_UR_LL(1000)
166 μs ± 489 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [4]: %timeit triangle_half_UR_LL(1000)
166 μs ± 270 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [5]: %timeit triangle_half_UR_LL(1000)
166 μs ± 506 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

Agora, se eu definir outra função e compilá-la JIT com o Numba, o desempenho da função rápida cai inexplicavelmente:

In [6]: @nb.njit(cache=True)
   ...: def dummy():
   ...:     pass

In [7]: dummy()

In [8]: %timeit triangle_half_UR_LL(1000)
980 μs ± 20 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [9]: %timeit triangle_half_UR_LL(1000)
976 μs ± 9.9 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [10]: %timeit triangle_half_UR_LL(1000)
974 μs ± 3.11 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

Isso é real, já reproduzi esse problema com sucesso muitas vezes, sem falhas. Inicio uma nova sessão do interpretador, colo o código e ele roda rápido. Defino a função fictícia, chamo a função fictícia e a função rápida fica inexplicavelmente mais lenta.

Captura de tela como prova:

Estou usando o Windows 11 e não tenho a mínima ideia do que está acontecendo.

Existe alguma explicação para isso? E como posso evitar esse problema?

Curiosamente, se eu me livrar do nogilparâmetro e sem alterar mais nada, o problema desaparece magicamente:

In [1]: import numpy as np
   ...: import numba as nb
   ...:
   ...:
   ...: @nb.njit(cache=True, parallel=True)
   ...: def triangle_half_UR_LL(size: int, swap: bool = False) -> tuple[np.ndarray, np.ndarray]:
   ...:     total = (size + 1) * size // 2
   ...:     x_coords = np.full(total, 0, dtype=np.uint16)
   ...:     y_coords = np.full(total, 0, dtype=np.uint16)
   ...:     offset = 0
   ...:     side = np.arange(size, dtype=np.uint16)
   ...:     for i in nb.prange(size):
   ...:         offset = i * size - (i - 1) * i // 2
   ...:         end = offset + size - i
   ...:         x_coords[offset:end] = i
   ...:         y_coords[offset:end] = side[i:]
   ...:
   ...:     return (x_coords, y_coords) if not swap else (y_coords, x_coords)

In [2]: %timeit triangle_half_UR_LL(1000)
186 μs ± 47.9 μs per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [3]: %timeit triangle_half_UR_LL(1000)
167 μs ± 1.61 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [4]: %timeit triangle_half_UR_LL(1000)
166 μs ± 109 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [5]: @nb.njit(cache=True)
   ...: def dummy():
   ...:     pass

In [6]: dummy()

In [7]: %timeit triangle_half_UR_LL(1000)
167 μs ± 308 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [8]: %timeit triangle_half_UR_LL(1000)
166 μs ± 312 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [9]: %timeit triangle_half_UR_LL(1000)
167 μs ± 624 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

Por que isso acontece?

Mas não, se eu definir outras funções, de alguma forma a primeira função fica lenta novamente. A maneira mais simples de reproduzir o problema é simplesmente redefinindo-a:

In [7]: dummy()

In [8]: %timeit triangle_half_UR_LL(1000)
168 μs ± 750 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [9]: import numpy as np

In [10]: %timeit triangle_half_UR_LL(1000)
167 μs ± 958 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [11]: import numba as nb

In [12]: %timeit triangle_half_UR_LL(1000)
167 μs ± 311 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [13]: @nb.njit(cache=True, parallel=True)
    ...: def triangle_half_UR_LL(size: int, swap: bool = False) -> tuple[np.ndarray, np.ndarray]:
    ...:     total = (size + 1) * size // 2
    ...:     x_coords = np.full(total, 0, dtype=np.uint16)
    ...:     y_coords = np.full(total, 0, dtype=np.uint16)
    ...:     offset = 0
    ...:     side = np.arange(size, dtype=np.uint16)
    ...:     for i in nb.prange(size):
    ...:         offset = i * size - (i - 1) * i // 2
    ...:         end = offset + size - i
    ...:         x_coords[offset:end] = i
    ...:         y_coords[offset:end] = side[i:]
    ...:
    ...:     return (x_coords, y_coords) if not swap else (y_coords, x_coords)

In [14]: %timeit triangle_half_UR_LL(1000)
1.01 ms ± 94.3 μs per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [15]: %timeit triangle_half_UR_LL(1000)
964 μs ± 2.02 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

A lentidão também acontece se eu definir a seguinte função e chamá-la:

@nb.njit(cache=True)
def Farey_sequence(n: int) -> np.ndarray:
    a, b, c, d = 0, 1, 1, n
    result = [(a, b)]
    while 0 <= c <= n:
        k = (n + b) // d
        a, b, c, d = c, d, k * c - a, k * d - b
        result.append((a, b))

    return np.array(result, dtype=np.uint64)

In [6]: %timeit triangle_half_UR_LL(1000)
166 μs ± 296 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [7]: %timeit Farey_sequence(16)
The slowest run took 6.25 times longer than the fastest. This could mean that an intermediate result is being cached.
6.03 μs ± 5.72 μs per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [8]: %timeit Farey_sequence(16)
2.77 μs ± 50.8 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

In [9]: %timeit triangle_half_UR_LL(1000)
966 μs ± 6.48 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

Dado um inteiro positivo N, podemos rotular todos os pontos da grade no quadrado N x N, começando em 1, o número total de pontos da grade é N x N, e os pontos da grade são list(itertools.product(range(1, N + 1), repeat=2)).

Agora, quero encontrar todas as tuplas (x, y)que satisfazem a condição x/y ser uma fração não reduzida. A seguir, uma implementação de força bruta que é garantidamente correta, mas é muito ineficiente:

import math
from itertools import product


def find_complex_points(lim: int) -> list[tuple[int, int]]:
    return [
        (x, y)
        for x, y in product(range(1, lim + 1), repeat=2)
        if math.gcd(x, y) > 1
    ]

Agora, a próxima função é um pouco mais inteligente, mas gera duplicatas e, como resultado, é apenas visivelmente mais rápida, mas não muito:

def find_complex_points_1(lim: int) -> set[tuple[int, int]]:
    lim += 1
    return {
        (x, y)
        for mult in range(2, lim)
        for x, y in product(range(mult, lim, mult), repeat=2)
    }

In [255]: %timeit find_complex_points(1024)
233 ms ± 4.44 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [256]: %timeit find_complex_points_1(1024)
194 ms ± 1.9 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

Existe uma maneira melhor de fazer isso?

(Meu objetivo é simples, quero criar uma matriz NumPy 2D do tipo uint8 com a forma (N, N), preenchê-la com 255 e tornar todos os pixels (x, y) 0 se (x+1)/(y+1) for uma fração não reduzida)

Eu criei um método que é muito mais inteligente do que os meus dois anteriores, e também tremendamente mais rápido, mas ele ainda gera duplicatas. Optei por não usar setaqui para que você possa copiar e colar o código como está e executar alguns testes e ver a saída exata na ordem em que são gerados:

def find_complex_points_2(lim: int) -> set[tuple[int, int]]:
    stack = dict.fromkeys(range(lim, 1, -1))
    lim += 1
    points = []
    while stack:
        x, _ = stack.popitem()
        points.append((x, x))
        mults = []
        for y in range(x * 2, lim, x):
            stack.pop(y, None)
            mults.append(y)
            points.extend([(x, y), (y, x)])
        
        for i, x in enumerate(mults):
            points.append((x, x))
            for y in mults[i + 1:]:
                points.extend([(x, y), (y, x)])
    
    return points

In [292]: sorted(set(find_complex_points_2(1024))) == find_complex_points(1024)
Out[292]: True

In [293]: %timeit find_complex_points_2(1024)
58.9 ms ± 580 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [294]: %timeit find_complex_points(1024)
226 ms ± 3.24 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Para esclarecer, a saída de find_complex_points_2(10)é:

In [287]: find_complex_points_2(10)
Out[287]:
[(2, 2),
 (2, 4),
 (4, 2),
 (2, 6),
 (6, 2),
 (2, 8),
 (8, 2),
 (2, 10),
 (10, 2),
 (4, 4),
 (4, 6),
 (6, 4),
 (4, 8),
 (8, 4),
 (4, 10),
 (10, 4),
 (6, 6),
 (6, 8),
 (8, 6),
 (6, 10),
 (10, 6),
 (8, 8),
 (8, 10),
 (10, 8),
 (10, 10),
 (3, 3),
 (3, 6),
 (6, 3),
 (3, 9),
 (9, 3),
 (6, 6),
 (6, 9),
 (9, 6),
 (9, 9),
 (5, 5),
 (5, 10),
 (10, 5),
 (10, 10),
 (7, 7)]

Como você pode ver, (10, 10)aparece duas vezes. Quero evitar cálculos redundantes.

Isso também acontece em find_complex_points_1, se eu não usar um conjunto, muitas duplicatas serão incluídas, porque o método usado inevitavelmente as gerará repetidamente, ao usar um conjunto ainda há computação desnecessária, ele simplesmente não coleta as duplicatas.

E não, na verdade eu quero que as coordenadas sejam substituídas pela soma de todos os números anteriores, então N é substituído por (N ² + N) / 2.

Acabei de implementar a geração de imagens para ilustrar melhor o que quero:

import numpy as np
import numba as nb


@nb.njit(cache=True)
def resize_img(img: np.ndarray, h_scale: int, w_scale: int) -> np.ndarray:
    height, width = img.shape
    result = np.empty((height, h_scale, width, w_scale), np.uint8)
    result[...] = img[:, None, :, None]
    return result.reshape((height * h_scale, width * w_scale))


def find_composite_points(lim: int) -> set[tuple[int, int]]:
    stack = dict.fromkeys(range(lim, 1, -1))
    lim += 1
    points = set()
    while stack:
        x, _ = stack.popitem()
        points.add((x, x))
        mults = []
        for y in range(x * 2, lim, x):
            stack.pop(y, None)
            mults.append(y)
            points.update([(x, y), (y, x)])

        for i, x in enumerate(mults):
            points.add((x, x))
            for y in mults[i + 1 :]:
                points.update([(x, y), (y, x)])

    return points


def natural_sum(n: int) -> int:
    return (n + 1) * n // 2


def composite_image(lim: int, scale: int) -> np.ndarray:
    length = natural_sum(lim)
    img = np.full((length, length), 255, dtype=np.uint8)
    for x, y in find_composite_points(lim):
        x1, y1 = natural_sum(x - 1), natural_sum(y - 1)
        img[x1 : x1 + x, y1 : y1 + y] = 0

    return resize_img(img, scale, scale)

composite_image(12, 12)

Desculpe pelo título longo. Não sei se isso é mais um problema de matemática ou de programação, mas acho que minha matemática está extremamente enferrujada e sou melhor em programação.

Então eu tenho essa expansão de fração contínua do arco tangente:

Eu peguei da Wikipedia

Tentei encontrar um algoritmo simples para calculá-lo:

E eu consegui, escrevi uma implementação de precisão infinita da expansão de frações contínuas sem usar nenhuma biblioteca, usando apenas aritmética básica de inteiros:

import json
import math
import random
from decimal import Decimal, getcontext
from typing import Callable, List, Tuple


Fraction = Tuple[int, int]


def arctan_cf(y: int, x: int, lim: int) -> Fraction:
    y_sq = y**2
    a1, a2 = y, 3 * x * y
    b1, b2 = x, 3 * x**2 + y_sq
    odd = 5
    for i in range(2, 2 + lim):
        t1, t2 = odd * x, i**2 * y_sq
        a1, a2 = a2, t1 * a2 + t2 * a1
        b1, b2 = b2, t1 * b2 + t2 * b1
        odd += 2

    return a2, b2

E converge mais rápido que a série arco-tangente de Newton que usei anteriormente .

Agora acho que se eu combinar isso com a fórmula do meio-ângulo do arco tangente, a convergência deve ser mais rápida.

def half_arctan_cf(y: int, x: int, lim: int) -> Fraction:
    c = (x**2 + y**2) ** 0.5
    a, b = c.as_integer_ratio()
    a, b = arctan_cf(a - b * x, b * y, lim)
    return 2 * a, b

E, de fato, converge ainda mais rápido:

def test_accuracy(lim: int) -> dict:
    result = {}
    for _ in range(lim):
        x, y = random.sample(range(1024), 2)
        while not x or not y:
            x, y = random.sample(range(1024), 2)

        atan2 = math.atan2(y, x)
        entry = {"atan": atan2}
        for fname, func in zip(
            ("arctan_cf", "half_arctan_cf"), (arctan_cf, half_arctan_cf)
        ):
            i = 1
            while True:
                a, b = func(y, x, i)
                if math.isclose(deci := a / b, atan2):
                    break

                i += 1

            entry[fname] = (i, deci)

        result[f"{y} / {x}"] = entry

    return result


print(json.dumps(test_accuracy(8), indent=4))

for v in test_accuracy(128).values():
    assert v["half_arctan_cf"][0] <= v["arctan_cf"][0]

{
    "206 / 136": {
        "atan": 0.9872880750087898,
        "arctan_cf": [
            16,
            0.9872880746658675
        ],
        "half_arctan_cf": [
            6,
            0.9872880746018052
        ]
    },
    "537 / 308": {
        "atan": 1.0500473287277563,
        "arctan_cf": [
            18,
            1.0500473281360896
        ],
        "half_arctan_cf": [
            7,
            1.0500473288158192
        ]
    },
    "331 / 356": {
        "atan": 0.7490241118247137,
        "arctan_cf": [
            10,
            0.7490241115996227
        ],
        "half_arctan_cf": [
            5,
            0.749024111913438
        ]
    },
    "744 / 613": {
        "atan": 0.8816364228048325,
        "arctan_cf": [
            13,
            0.8816364230439662
        ],
        "half_arctan_cf": [
            6,
            0.8816364227495634
        ]
    },
    "960 / 419": {
        "atan": 1.1592605364805093,
        "arctan_cf": [
            24,
            1.1592605359263286
        ],
        "half_arctan_cf": [
            7,
            1.1592605371181872
        ]
    },
    "597 / 884": {
        "atan": 0.5939827714677137,
        "arctan_cf": [
            7,
            0.5939827719895824
        ],
        "half_arctan_cf": [
            4,
            0.59398277135389
        ]
    },
    "212 / 498": {
        "atan": 0.40246578425167584,
        "arctan_cf": [
            5,
            0.4024657843859885
        ],
        "half_arctan_cf": [
            3,
            0.40246578431841773
        ]
    },
    "837 / 212": {
        "atan": 1.322727785860997,
        "arctan_cf": [
            41,
            1.322727786922624
        ],
        "half_arctan_cf": [
            8,
            1.3227277847674388
        ]
    }
}

Esse bloco assert é bastante longo para um grande número de amostras, mas nunca gera exceções.

Então, acho que posso usar a expansão de fração contínua do arco tangente com séries do tipo Machin para calcular π. (Usei a última série na seção vinculada porque ela converge mais rápido)

def sum_fractions(fractions: List[Fraction]) -> Fraction:
    while (length := len(fractions)) > 1:
        stack = []
        for i in range(0, length - (odd := length & 1), 2):
            num1, den1 = fractions[i]
            num2, den2 = fractions[i + 1]
            stack.append((num1 * den2 + num2 * den1, den1 * den2))

        if odd:
            stack.append(fractions[-1])

        fractions = stack

    return fractions[0]


MACHIN_SERIES = ((44, 57), (7, 239), (-12, 682), (24, 12943))


def approximate_loop(lim: int, func: Callable) -> List[Fraction]:
    fractions = []
    for coef, denom in MACHIN_SERIES:
        dividend, divisor = func(1, denom, lim)
        fractions.append((coef * dividend, divisor))

    return fractions


def approximate_1(lim: int) -> List[Fraction]:
    return approximate_loop(lim, arctan_cf)


def approximate_2(lim: int) -> List[Fraction]:
    return approximate_loop(lim, half_arctan_cf)


approx_funcs = (approximate_1, approximate_2)


def calculate_pi(lim: int, approx: bool = 0) -> Fraction:
    dividend, divisor = sum_fractions(approx_funcs[approx](lim))
    dividend *= 4
    return dividend // (common := math.gcd(dividend, divisor)), divisor // common

getcontext().rounding = 'ROUND_DOWN'
def to_decimal(dividend: int, divisor: int, places: int) -> str:
    getcontext().prec = places + len(str(dividend // divisor))
    return str(Decimal(dividend) / Decimal(divisor))


def get_accuracy(lim: int, approx: bool = 0) -> Tuple[int, str]:
    length = 12
    fraction = calculate_pi(lim, approx)
    while True:
        decimal = to_decimal(*fraction, length)
        for i, e in enumerate(decimal):
            if Pillion[i] != e:
                return (max(0, i - 2), decimal[:i])

        length += 10


with open("D:/Pillion.txt", "r") as f:
    Pillion = f.read()

Pillion.txt contém os primeiros 1000001 dígitos de π, Pi + Milhão = Pillion.

E funciona, mas apenas parcialmente. A expansão básica de fração contínua funciona muito bem com fórmulas do tipo Machin, mas combinadas com fórmulas de meio ângulo, eu só consigo obter 9 casas decimais corretas não importa o que aconteça, e de fato, eu obtenho 9 dígitos corretos na primeira iteração, e então essa coisa toda não melhora nunca:

In [2]: get_accuracy(16)
Out[2]:
(73,
 '3.1415926535897932384626433832795028841971693993751058209749445923078164062')

In [3]: get_accuracy(32)
Out[3]:
(138,
 '3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231')

In [4]: get_accuracy(16, 1)
Out[4]: (9, '3.141592653')

In [5]: get_accuracy(32, 1)
Out[5]: (9, '3.141592653')

In [6]: get_accuracy(1, 1)
Out[6]: (9, '3.141592653')

Mas os dígitos de fato mudam:

In [7]: to_decimal(*calculate_pi(1, 1), 32)
Out[7]: '3.14159265360948500093515231500093'

In [8]: to_decimal(*calculate_pi(2, 1), 32)
Out[8]: '3.14159265360945286794831052938917'

In [9]: to_decimal(*calculate_pi(3, 1), 32)
Out[9]: '3.14159265360945286857612896472974'

In [10]: to_decimal(*calculate_pi(4, 1), 32)
Out[10]: '3.14159265360945286857611676794770'

In [11]: to_decimal(*calculate_pi(5, 1), 32)
Out[11]: '3.14159265360945286857611676818392'

Por que a fração contínua com fórmula de meio-ângulo não está funcionando com fórmula tipo Machin? E é possível fazê-la funcionar, e se puder funcionar, então como? Eu quero uma prova de que é impossível, ou um exemplo funcional que prove que é possível.

Apenas uma verificação de sanidade, usando π/4 = arctan(1) consegui gerar half_arctan_cfdígitos de π, mas ele converge muito mais lentamente:

def approximate_3(lim: int) -> List[Fraction]:
    return [half_arctan_cf(1, 1, lim)]


approx_funcs = (approximate_1, approximate_2, approximate_3)

In [28]: get_accuracy(16, 2)
Out[28]: (15, '3.141592653589793')

In [29]: get_accuracy(16, 0)
Out[29]:
(73,
 '3.1415926535897932384626433832795028841971693993751058209749445923078164062')

E o mesmo problema se repete, ele atinge a precisão máxima de 15 dígitos na 10ª iteração:

In [37]: get_accuracy(9, 2)
Out[37]: (14, '3.14159265358979')

In [38]: get_accuracy(10, 2)
Out[38]: (15, '3.141592653589793')

In [39]: get_accuracy(11, 2)
Out[39]: (15, '3.141592653589793')

In [40]: get_accuracy(32, 2)
Out[40]: (15, '3.141592653589793')

Acabei de reescrever minha implementação de fração contínua arco-tangente e fiz com que ela evitasse cálculos redundantes.

No meu código, em cada iteração, t1 aumenta em 2 * y_sq, então não há necessidade de multiplicar repetidamente y_sq pelo número ímpar; em vez disso, basta usar uma variável cumulativa e um passo de 2 * y_sq.

E a diferença entre cada par de números quadrados consecutivos são apenas os números ímpares, então posso usar uma variável cumulativa de uma variável cumulativa.

def arctan_cf_0(y: int, x: int, lim: int) -> Fraction:
    y_sq = y**2
    a1, a2 = y, 3 * x * y
    b1, b2 = x, 3 * x**2 + y_sq
    odd = 5
    for i in range(2, 2 + lim):
        t1, t2 = odd * x, i**2 * y_sq
        a1, a2 = a2, t1 * a2 + t2 * a1
        b1, b2 = b2, t1 * b2 + t2 * b1
        odd += 2

    return a2, b2


def arctan_cf(y: int, x: int, lim: int) -> Fraction:
    y_sq = y**2
    a1, a2 = y, 3 * x * y
    b1, b2 = x, 3 * x**2 + y_sq
    t1_step, t3_step = 2 * x, 2 * y_sq
    t1, t2 = 5 * x, 4 * y_sq
    t3 = t2 + y_sq
    for _ in range(lim):
        a1, a2 = a2, t1 * a2 + t2 * a1
        b1, b2 = b2, t1 * b2 + t2 * b1
        t1 += t1_step
        t2 += t3
        t3 += t3_step

    return a2, b2

In [301]: arctan_cf_0(4, 3, 100) == arctan_cf(4, 3, 100)
Out[301]: True

In [302]: %timeit arctan_cf_0(4, 3, 100)
58.6 μs ± 503 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [303]: %timeit arctan_cf(4, 3, 100)
54.3 μs ± 816 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

Embora isso não melhore muito a velocidade, é definitivamente uma melhoria.

Isso pode parecer trivial, mas não encontrei uma boa solução para o problema. Eu até encontrei isso: gere todos os números binários de n bits da maneira mais rápida possível . mas não encontrei uma duplicata exata.

O problema é muito simples: dado um limite N que seja um inteiro positivo, gere a representação binária de cada número natural até N em ordem (excluindo N, o primeiro número natural é 0, então N - 1 é o número máximo a ser representado), na tupleforma, com cada tupla preenchida com zeros à esquerda para que todas as representações tenham o mesmo comprimento.

Por exemplo, se Nfor 4, a saída deverá ser [(0, 0), (0, 1), (1, 0), (1, 1)].

Neste ponto, o problema é de fato trivial, mas há um porém: não são permitidos " bin(n)e f'{n:b}'" e "afins". O algoritmo deve operar inteiramente no domínio binário, porque, pelo que entendi, tudo (texto, fotos, música, vídeos...) em computadores são numerais binários, então converter representações de um lado para o outro está adicionando cálculos desnecessários. Esses cálculos (conversão de base) são completamente redundantes e devem ser eliminados para produzir o programa mais eficiente (trata-se de manter os problemas específicos para o menor número possível de domínios, para que possamos operar apenas nesses domínios).

Eu escrevi um programa simples que faz exatamente o que descrevo:

from typing import Generator, Tuple

def count_in_binary(n: int) -> Generator[Tuple[int, ...], None, None]:
    if not isinstance(n, int) or n < 1:
        raise ValueError("The argument n must be a positive integer")

    l = (n - 1).bit_length() if n > 1 else 1
    numeral = [0] * l
    maxi = l - 1
    for _ in range(n):
        yield tuple(numeral)
        i = maxi
        while True:
            if not (d := numeral[i]):
                numeral[i] = 1
                break
            else:
                numeral[i] = 0
                i -= 1

Mas não tenho certeza se essa é a maneira mais eficiente de fazer isso em Python. Não usei muitas operações de bits e os computadores já representam números como binários, então deve haver maneiras mais rápidas de fazer isso.

A questão é: qual é o caminho mais rápido?

Para efeito de comparação, aqui está uma maneira que usa f'{n:b}' e, portanto, é mais concisa, mas na verdade é muito mais lenta e estúpida:

def count_in_binary1(n: int) -> Generator[Tuple[int, ...], None, None]:
    if not isinstance(n, int) or n < 1:
        raise ValueError("The argument n must be a positive integer")

    l = len(f'{n-1:b}')
    for i in range(n):
        yield tuple(map(int, f'{i:0{l}b}'))

In [50]: %timeit list(count_in_binary(256))
59.9 μs ± 209 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [51]: %timeit list(count_in_binary1(256))
452 μs ± 3.68 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

Editar

Não fiz muitos testes com a função original, apenas pensei que funcionaria, agora está corrigido.

E não, o escopo do escopo é limitado ao Python puro, então o NumPy não é permitido.

Editar 2

Agora acho que não há exceções.

Votei positivamente e aceitei a segunda resposta, pois ela responde à pergunta original, embora a pergunta original não inclua todas as informações relevantes, então o problema que eu estava tentando resolver postando uma pergunta continua sem solução.

Publiquei uma nova pergunta sem todas as informações relevantes. Por favor, responda: Qual é a maneira mais rápida de encontrar todas as permutações de 0, 1 de largura n sem itertools em Python puro?

Estou tentando implementar algoritmos eficientes de fatoração de primos em Python. Isso não é tarefa de casa ou trabalho, é completamente por curiosidade.

Aprendi que a fatoração prima é muito difícil :

Quero implementar algoritmos eficientes para isso como um desafio autoimposto. Eu me propus a implementar o método de fatoração de Fermat primeiro, pois parece simples o suficiente.

Código Python traduzido diretamente do pseudocódigo:

def Fermat_Factor(n):
    a = int(n ** 0.5 + 0.5)
    b2 = abs(a**2 - n)
    while int(b2**0.5) ** 2 != b2:
        a += 1
        b2 = a**2 - n

    return a - b2**0.5, a + b2**0.5

(Tenho que usar, abscaso contrário b2, será facilmente negativo e into cast falhará TypeErrorporque a raiz é complex)

Como você pode ver, ele retorna dois inteiros cujo produto é igual à entrada, mas ele retorna apenas duas saídas e não garante a primalidade dos fatores. Não tenho ideia de quão eficiente esse algoritmo é, mas a fatoração de semiprimos usando esse método é muito mais eficiente do que o método de divisão por tentativa usado na minha pergunta anterior: Por que a fatoração de produtos de primos próximos é muito mais lenta do que produtos de primos diferentes .

In [20]: %timeit FermatFactor(3607*3803)
2.1 μs ± 28.2 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

In [21]: FermatFactor(3607*3803)
Out[21]: [3607, 3803]

In [22]: %timeit FermatFactor(3593 * 3671)
1.69 μs ± 31 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

In [23]: FermatFactor(3593 * 3671)
Out[23]: [3593, 3671]

In [24]: %timeit FermatFactor(7187 * 7829)
4.94 μs ± 47.4 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

In [25]: FermatFactor(7187 * 7829)
Out[25]: [7187, 7829]

In [26]: %timeit FermatFactor(8087 * 8089)
1.38 μs ± 12.9 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

In [27]: FermatFactor(8087 * 8089)
Out[27]: [8087, 8089]

Então eu quero usar esse algoritmo para gerar todos os fatores primos de um inteiro qualquer dado (claro que eu sei que isso só funciona com inteiros ímpares, mas isso não é um problema, já que potências de dois podem ser fatoradas trivialmente usando hacking de bits). A maneira mais fácil que eu consigo pensar é chamar recursivamente Fermat_Factoraté que nseja primo. Eu não sei como verificar se um número é primo nesse algoritmo, mas eu notei algo:

In [3]: Fermat_Factor(3)
Out[3]: (1.0, 3.0)

In [4]: Fermat_Factor(5)
Out[4]: (1.0, 3.0)

In [5]: Fermat_Factor(7)
Out[5]: (1.0, 7.0)

In [6]: Fermat_Factor(11)
Out[6]: (1.0, 11.0)

In [7]: Fermat_Factor(13)
Out[7]: (1.0, 13.0)

In [8]: Fermat_Factor(17)
Out[8]: (3.0, 5.0)

In [9]: Fermat_Factor(19)
Out[9]: (1.0, 19.0)

In [10]: Fermat_Factor(23)
Out[10]: (1.0, 23.0)

In [11]: Fermat_Factor(29)
Out[11]: (3.0, 7.0)

In [12]: Fermat_Factor(31)
Out[12]: (1.0, 31.0)

In [13]: Fermat_Factor(37)
Out[13]: (5.0, 7.0)

In [14]: Fermat_Factor(41)
Out[14]: (1.0, 41.0)

O primeiro número na saída deste algoritmo para muitos primos é 1, mas não para todos, como tal, não pode ser usado para determinar quando a recursão deve parar. Aprendi isso da maneira mais difícil.

Então, eu apenas decidi usar a verificação de associação de um conjunto pré-gerado de primos. Naturalmente, isso ocorrerá RecursionError: maximum recursion depth exceededquando a entrada for um primo maior que o máximo do conjunto. Como não tenho memória infinita, isso deve ser considerado um detalhe de implementação.

Então implementei uma versão funcional (para algumas entradas), mas para algumas entradas válidas (produtos de primos dentro do limite), de alguma forma o algoritmo não fornece a saída correta:

import numpy as np
from itertools import cycle

TRIPLE = ((4, 2), (9, 6), (25, 10))
WHEEL = ( 4, 2, 4, 2, 4, 6, 2, 6 )
def prime_sieve(n):
    primes = np.ones(n + 1, dtype=bool)
    primes[:2] = False
    for square, double in TRIPLE:
        primes[square::double] = False
    
    wheel = cycle(WHEEL)
    k = 7
    while (square := k**2) <= n:
        if primes[k]:
            primes[square::2*k] = False
        
        k += next(wheel)
    
    return np.flatnonzero(primes)
    
PRIMES = list(map(int, prime_sieve(1048576)))
PRIME_SET = set(PRIMES)
TEST_LIMIT = PRIMES[-1] ** 2

def FermatFactor(n):
    if n > TEST_LIMIT:
        raise ValueError('Number too large')
    
    if n in PRIME_SET:
        return [n]
    
    a = int(n ** 0.5 + 0.5)
    if a ** 2 == n:
        return FermatFactor(a) + FermatFactor(a)
    
    b2 = abs(a**2 - n)
    while int(b2**0.5) ** 2 != b2:
        a += 1
        b2 = a**2 - n
    
    return FermatFactor(factor := int(a - b2**0.5)) + FermatFactor(n // factor)

Funciona para muitas entradas:

In [18]: FermatFactor(255)
Out[18]: [3, 5, 17]

In [19]: FermatFactor(511)
Out[19]: [7, 73]

In [20]: FermatFactor(441)
Out[20]: [3, 7, 3, 7]

In [21]: FermatFactor(3*5*823)
Out[21]: [3, 5, 823]

In [22]: FermatFactor(37*333667)
Out[22]: [37, 333667]

In [23]: FermatFactor(13 * 37 * 151 * 727 * 3607)
Out[23]: [13, 37, 727, 151, 3607]

Mas não todos:

In [25]: FermatFactor(5 * 53 * 163)
Out[25]: [163, 13, 2, 2, 5]

In [26]: FermatFactor(3*5*73*283)
Out[26]: [17, 3, 7, 3, 283]

In [27]: FermatFactor(3 * 11 * 29 * 71 *  137)
Out[27]: [3, 11, 71, 61, 7, 3, 3]

Por que é esse o caso? Como posso consertar?

Esta é uma questão puramente acadêmica, sem nenhuma consideração prática. Isto não é dever de casa, eu abandonei o ensino médio há muito tempo. Estou apenas curioso, e não consigo dormir bem sem saber o porquê.

Eu estava brincando com Python. Decidi fatorar números inteiros grandes e medir o tempo de execução das chamadas para cada entrada.

Usei vários números e descobri que alguns demoram muito mais para serem fatorados do que outros.

Decidi então investigar mais a fundo, escrevi rapidamente uma função de peneira de primos para gerar primos para teste. Descobri que um produto de um par de primos moderadamente grandes (dois primos de quatro dígitos) leva muito mais tempo para ser fatorado do que um produto de um primo muito grande (seis dígitos+) e um primo pequeno (<=três dígitos).

A princípio, pensei que minha primeira função simples para teste era ineficiente, e esse é realmente o caso. Então, escrevi uma segunda função que extraía números primos diretamente de uma lista pré-gerada de números primos. A segunda função era de fato mais eficiente, mas estranhamente ela exibe o mesmo padrão.

Aqui estão alguns números que usei:

13717421 == 3607 * 3803
13189903 == 3593 * 3671
56267023 == 7187 * 7829
65415743 == 8087 * 8089

12345679 == 37 * 333667
38760793 == 37 * 1047589
158202851 == 151 * 1047701
762312571 == 727 * 1048573

Código:

import numpy as np
from itertools import cycle

def factorize(n):
    factors = []
    while not n % 2:
        factors.append(2)
        n //= 2

    i = 3
    while i**2 <= n:
        while not n % i:
            factors.append(i)
            n //= i
        i += 2
    
    return factors if n == 1 else factors + [n]

TRIPLE = ((4, 2), (9, 6), (25, 10))
WHEEL = ( 4, 2, 4, 2, 4, 6, 2, 6 )
def prime_sieve(n):
    primes = np.ones(n + 1, dtype=bool)
    primes[:2] = False
    for square, double in TRIPLE:
        primes[square::double] = False
    
    wheel = cycle(WHEEL)
    k = 7
    while (square := k**2) <= n:
        if primes[k]:
            primes[square::2*k] = False
        
        k += next(wheel)
    
    return np.flatnonzero(primes)
    
PRIMES = list(map(int, prime_sieve(1048576)))
TEST_LIMIT = PRIMES[-1] ** 2

def factorize_sieve(n):
    if n > TEST_LIMIT:
        raise ValueError('Number too large')

    factors = []
    for p in PRIMES:
        if p**2 > n:
            break
        while not n % p:
            factors.append(p)
            n //= p
    
    return factors if n == 1 else factors + [n]

Resultado do teste:

In [2]: %timeit factorize(13717421)
279 μs ± 4.29 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [3]: %timeit factorize(12345679)
39.6 μs ± 749 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [4]: %timeit factorize_sieve(13717421)
64.1 μs ± 688 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [5]: %timeit factorize_sieve(12345679)
12.6 μs ± 146 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

In [6]: %timeit factorize_sieve(13189903)
64.6 μs ± 964 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [7]: %timeit factorize_sieve(56267023)
117 μs ± 3.88 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [8]: %timeit factorize_sieve(65415743)
130 μs ± 1.38 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [9]: %timeit factorize_sieve(38760793)
21.1 μs ± 232 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [10]: %timeit factorize_sieve(158202851)
21.4 μs ± 385 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [11]: %timeit factorize_sieve(762312571)
22.1 μs ± 409 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

Como você pode ver claramente, a fatoração de dois primos médios leva, em média, muito mais tempo do que dois extremos. Por que é esse o caso?

Estou tentando converter dados integrais não assinados em sua representação binária na memória da maneira mais eficiente possível.

Escrevi quatro funções de modelo para converter os inteiros em little endian e big endian, duas delas usam operações de bits e as outras duas usam ponteiros para copiar dados.

Eles foram verificados como corretos e também muito eficientes, pois determinei que as funções little endian são tão rápidas quanto std::memcpy, mas as funções big endian de alguma forma demoram um pouco mais.

Essas funções são:

#include <vector>

using std::vector;
typedef vector<uint8_t> bytes;

template<class T>
inline bytes LittleEndian(const T& data) {
    size_t size = sizeof(T);
    bytes _bytes(size);
    uint8_t mask = 255;
    for (size_t i = 0, shift = 0; i < size; i++, shift += 8) {
        _bytes[i] = (data >> shift) & mask;
    }
    return _bytes;
}

template<class T>
inline bytes BigEndian(const T& data) {
    size_t size = sizeof(T);
    bytes _bytes(size);
    uint8_t mask = 255;
    for (size_t i = size, shift = 0; i-- > 0; shift += 8) {
        _bytes[i] = (data >> shift) & mask;
    }
    return _bytes;
}

template<class T>
inline bytes CPU_Endian(const T& data) {
    size_t size = sizeof(T);
    bytes _bytes(size);
    uint8_t* dst = (uint8_t *)_bytes.data(), * src = (uint8_t *) & data;
    for (size_t i = 0; i < size; i++) {
        *dst++ = *src++;
    }
    return _bytes;
}

template<class T>
inline bytes Flip_CPU_Endian(const T& data) {
    size_t size = sizeof(T);
    bytes _bytes(size);
    uint8_t* dst = (uint8_t *)_bytes.data(), * src = (uint8_t *)&data + size - 1;
    for (size_t i = 0; i < size; i++) {
        *dst++ = *src--;
    }
    return _bytes;
}

E eu quero desenrolar os loops for usando std::index_sequence, e como eles estão relacionados, coloquei-os em uma pergunta. Tratam-se de três coisas: fazer algo repetidamente N vezes, criar uma sequência de índice que diminui em vez de aumentar e usar o índice para definir valores.

Eu tentei fazer isso sozinho, mas não funciona:

template<class T>
inline bytes CPU_Endian2(const T& data) {
    size_t size = sizeof(T);
    bytes _bytes(size);
    uint8_t* dst = (uint8_t*)_bytes.data(), * src = (uint8_t*)&data;
    [&]<std::size_t...N>(std::index_sequence<N...>){
        ((*dst++ = *src++),...);
    }(std::make_index_sequence<size>{});
    return _bytes;
}

Não compila, log de erros:

Build started at 18:54...
1>------ Build started: Project: hexlify_test, Configuration: Release x64 ------
1>hexlify_test.cpp
1>C:\Users\Estranger\source\repos\hexlify_test\hexlify_test.cpp(98,3): error C7515: a fold expression must contain an unexpanded parameter pack
1>C:\Users\Estranger\source\repos\hexlify_test\hexlify_test.cpp(99,3): error C3878: syntax error: unexpected token '(' following 'expression'
1>C:\Users\Estranger\source\repos\hexlify_test\hexlify_test.cpp(99,3): message : error recovery skipped: '( identifier ::  . . . {'
1>C:\Users\Estranger\source\repos\hexlify_test\hexlify_test.cpp(99,35): error C2760: syntax error: '}' was unexpected here; expected ';'
1>Done building project "hexlify_test.vcxproj" -- FAILED.
========== Build: 0 succeeded, 1 failed, 0 up-to-date, 0 skipped ==========
========== Build completed at 18:54 and took 01.796 seconds ==========

Como posso converter essas funções para aquelas que usam std::index_sequenceem vez de loops for?

Adicionar constexpra size_t size = sizeof(T);falhou ao compilá-lo.

Escrevi várias funções C++ simples para converter sequências de bytes em representações de strings.

Foi bem simples, tenho certeza que minha lógica está certa, achei extremamente fácil, até que comecei a imprimir as strings e descobri que a saída era um lixo:

#include <iostream>
#include <string>
#include <vector>

using std::vector;
typedef vector<uint8_t> bytes;
using std::string;
using std::cout;
using namespace std::literals;

string DIGITS = "0123456789abcdef"s;

static inline string hexlify(bytes arr) {
    string repr = ""s;
    for (auto& chr : arr) {
        repr += " " + DIGITS[(chr & 240) >> 4] + DIGITS[chr & 15];
    }
    repr.erase(0, 1);
    return repr;
}

bytes text = {
    84, 111, 32, 98, 101, 32,
    111, 114, 32, 110, 111, 116,
    32, 116, 111, 32, 98, 101
}; // To be or not to be

int main() {
    cout << hexlify(text);
}

2♠
÷82♠
÷82♠
÷82♠
÷

Por que isso está acontecendo?

Eu sei que minha lógica está certa, o seguinte é a tradução direta para Python:

digits = "0123456789abcdef"
def bytes_string(data):
    s = ""
    for i in data:
        s += " " + digits[(i & 240) >> 4] + digits[i & 15]
    return s[1:]

E funciona:

>>> bytes_string(b"To be or not to be")
'54 6f 20 62 65 20 6f 72 20 6e 6f 74 20 74 6f 20 62 65'

Mas por que não funciona em C++?

Estou usando o Visual Studio 2022 V17.9.7, sinalizadores do compilador:

/permissive- /ifcOutput "hexlify_test\x64\Release\" /GS /GL /W3 /Gy /Zc:wchar_t /Zi /Gm- /O2 /sdl /Fd"hexlify_test\x64\Release\vc143.pdb" /Zc:inline /fp:precise /D "NDEBUG" /D "_CONSOLE" /D "_UNICODE" /D "UNICODE" /errorReport:prompt /WX- /Zc:forScope /std:c17 /Gd /Oi /MD /std:c++20 /FC /Fa"hexlify_test\x64\Release\" /EHsc /nologo /Fo"hexlify_test\x64\Release\" /Ot /Fp"hexlify_test\x64\Release\hexlify_test.pch" /diagnostics:column 

Acabei de descobrir que a saída de lixo só ocorre no modo de depuração após a implementação da correção, mirei o C++ 20 no modo de depuração, de alguma forma o código causa saída de lixo no modo de depuração, mudar para o modo de liberação corrige o problema. Antes da correção ser implementada eu compilei no modo de lançamento e houve esse problema.

Acho as funções matemáticas básicas como pow, e outras extremamente ineficientes, então decidi implementar minhas próprias versões delas, que são precisas e eficientes exp. E eu fiz.logcmath

Eu uso aproximações polinomiais de 9ª ordem no intervalo [1, 2) ou [0, 1) para aproximar as funções, usei numpy.polyfitpara obter os polinômios e apliquei os polinômios usando std::index_sequencewhich é incrivelmente rápido e, embora meus métodos não sejam tão precisos quanto os da biblioteca, posso obter pelo menos 13 casas decimais corretas com fp:precise, e meus métodos fornecem a mesma precisão que o código da biblioteca com fp:fast.

E agora quero uma função para medir o tempo de execução das funções e determinar a precisão delas.

Anteriormente, eu apenas copiava e colava as mesmas linhas de código para cada nova função que desejo avaliar, e isso não é muito organizado e está sujeito a erros.

Minhas funções são simples, todas retornam um número de ponto flutuante, um valor único, que é floatou double. E eles pegam um ou dois argumentos, o segundo argumento é um número, pode ser int, floatou doublee seu tipo não faz diferença. O tipo do primeiro argumento é importante, porém, pode ser floatou double, seu tipo é igual ao tipo de retorno da função e seu tipo deve corresponder exatamente ao que está declarado na função, porque eu uso std::bit_cast.

Como exemplo, o seguinte é o que eu criei para medir o tempo de execução de uma função que recebe um único floatparâmetro e retorna a float:

#include <chrono>
#include <functional>
#include <vector>

using std::chrono::steady_clock;
using std::chrono::duration;
using std::vector;
using std::function;

double timeit(const function<float(float)>& func, const vector<float>& values, int runs = 1048576){
    auto start = steady_clock::now();
    size_t len = values.size();
    for (int64_t i = 0; i < runs; i++) {
        func(values[i % len]);
    }
    auto end = steady_clock::now();
    duration<double, std::nano> time = end - start;
    return time.count() / runs;
}

Ainda não adicionei a verificação, mas adicionarei mais tarde.

Quero saber como posso reutilizar a mesma função para avaliar uma função que retorna a double? Eu sei da sobrecarga de funções, mas isso significaria uma função com o mesmo corpo da função, apenas com uma lista de parâmetros diferente. Acho que terei que usar outra função para os dois argumentos, mas não quero usar sobrecarga para funções com um argumento.

Eu sei, std::variantmas quero comparar os valores e, pelo que sei, não oferece suporte à comparação com tipos constituintes.

Como resolver isso?