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 2 + 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)
O algoritmo do Weeble é executado em
O(n² m²)ondemé o tamanho dos inteiros em bits (usando uma multiplicação ingênua). Como podemos assumir que a multiplicação de números é feita em tempo constante (devido aos inteiros nativos limitados usados pelo Numpy), isso significa,O(n²)mas com uma constante oculta significativa que não deve ser negligenciada. Em termos de desempenho, o algoritmo é limitado por operações de falha de página ineficientes e pelo preenchimento de grandes matrizes temporárias. Está longe de ser ótimo.O algoritmo de Pete Kirkham deve ser executado
O(n²)(embora seja difícil de provar formalmente) com uma constante oculta relativamente pequena. Esta é uma boa abordagem. No entanto, é muito lento devido à ineficiência das operações Numpy escalares, em vez das vetorizadas. Felizmente, ele pode ser facilmente vetorizado:Observe que corrigi a implementação para retornar resultados corretos (com valores 0 e 255).
Uma solução alternativa muito simples é usar o Numba para acelerar o código de Pete Kirkham. Dito isso, o código não é eficiente porque itera em itens de linhas diferentes no loop mais interno. Podemos corrigir isso facilmente trocando as variáveis:
Abordagem alternativa mais rápida
Observe que a matriz de saída é simétrica, então nem precisamos calcular a parte inferior esquerda. De fato,
gcd(a, b) == gcd(b, a). Infelizmente, não acho que possamos usar essa propriedade para tornar o código Numpy vetorizado, mas provavelmente podemos tornar o código Numba mais rápido.Além disso, a diagonal pode ser trivialmente definida como 0 (exceto o primeiro item), pois se
gcd(a, a) == a. Tecnicamente, também podemos trivialmente definir a vizinhança direta da diagonal (ie ) como 255, já que . deve ser preenchido com valores alternados (ie ) já que . Uma estratégia semelhante pode ser aplicada para . Para outras diagonais, a situação começa a ficar mais complexa, pois certamente precisamos fatorar números. Fatorar números é sabidamente caro, mas esse custo é amortizado aqui, já que precisamos fazer isso várias vezes.gcd(a, a) > 1a > 1array.diagonal(1)gcd(a, a-1) = 1array.diagonal(1)[255, 0, 255, 0, ...]gcd(a, a-2) = gcd(a, 2) = 2 - (a % 2)array.diagonal(2)O(n)Outra simetria do mdc é
gcd(a, b) = gcd(a, b-a) = gcd(a-b, b).Podemos aproveitar toda essa simetria do MDC para escrever uma implementação significativamente mais rápida usando programação dinâmica . Uma implementação ingênua (que combina todas as simetrias de forma bastante eficiente) é a seguinte:
Infelizmente, a transposição é bastante ineficiente e leva quase todo o tempo... Otimizar é possível, mas não é fácil. Podemos dividir a computação em blocos dependentes (semelhante ao funcionamento do algoritmo de decomposição de LU em bloco ). Isso torna o código muito mais complexo e rápido, graças a um padrão de acesso mais eficiente:
A transposição ainda é a parte mais lenta deste código. Ela pode ser otimizada ainda mais, mas às custas de um código significativamente maior. Prefiro manter isso razoavelmente simples, mas observe que uma maneira de tornar a transposição muito mais rápida é usar intrínsecos SIMD (certamente pelo menos >2 vezes mais rápido).
Referência
Aqui estão os resultados para
N=1024minha máquina (CPU i5-9600KF):O código Numba vetorizado é muito mais rápido que a outra implementação e os códigos Numba otimizados superam todas as implementações por uma grande margem (exceto o novo
find_complex_points_3)!É possível paralelizar alguns códigos do Numba para torná-los ainda mais rápidos, mas isso não é trivial e certamente é rápido o suficiente de qualquer maneira, sem mencionar que não escala bem porque o código é bastante limitado pela memória para grandes
N. Na verdade, uma cópia básica do Numpy leva cerca de 0,3 ms, o que pode ser considerado um tempo de execução de limite inferior.Como o que você quer é um array numpy, uma abordagem diferente seria começar com o array e usar algo como o crivo de Eratóstenes para marcar aqueles que não são reduzidos:
Até mesmo o algoritmo mais simples funciona muito mais rápido se você usar o numpy para todos os cálculos. Os loops e cálculos do Numpy são muito mais rápidos que os do CPython. Você pode usar numpy.gcd para os cálculos do MDC. Algo assim:
Para mim, isso é cerca de 7 vezes mais rápido que o original
find_complex_pointse cerca de 50% mais rápido que ofind_complex_points_2. Tenho certeza de que um algoritmo inteligente pode fazer ainda melhor, mas isso parece um bom ponto de complexidade/desempenho.Eu consegui, vetorizei meu algoritmo mais inteligente usando o Numba. Agora, como o resultado deveria ser um array NumPy, eu deveria gerar um array diretamente.
Agora, usar matrizes bidimensionais é ineficiente, usei uma matriz plana para garantir que o acesso à memória seja contínuo.
A segunda otimização que fiz foi, em vez de usar um dicionário e remover repetidamente números compostos, posso simplesmente usar uma matriz para rastrear quais números são marcados como compostos, o que garante que a variável de estado não mude de tamanho.
E então, uma grande fonte de duplicatas são as tuplas (x, x), que têm a garantia de serem geradas múltiplas vezes se o número composto tiver múltiplos divisores, então, em vez de gerar as tuplas repetidamente, simplesmente não as gere e gere todas essas tuplas no final para garantir que cada uma delas seja gerada exatamente uma vez.
Entretanto, isso não garante que nenhuma duplicata seja gerada. Não tenho ideia do porquê ainda pode haver duplicatas, mas como estou usando uma matriz NumPy, isso não importa.
Código
É mais rápido do que meu método anterior por uma margem muito ampla, porém ainda é menos eficiente do que o método mais rápido de Jérôme Richard :
Mas acho que posso otimizá-lo ainda mais.
E um exemplo de imagem de saída:
Eu me superei mais uma vez. Em vez de usar um array unidimensional, usei um array bidimensional aqui. Porque confirmei experimentalmente que usar um array unidimensional é de alguma forma menos eficiente do que usar um array bidimensional, possivelmente por causa de
.reshape.Agora, em vez de usar um método inteligente e evitar gerar duplicatas, basta fazer força bruta e atribuir os múltiplos a 0 um por um. Os valores atribuídos permanecem atribuídos e, no NumPy, a atribuição pode ser feita usando fatiamento, que é vetorizado e muito rápido.
Mas a principal otimização que economiza muito tempo é aplicar força bruta somente a números primos, porque não há necessidade alguma de aplicar força bruta a números compostos, já que eles são garantidos para serem processados ao processar pelo menos um número primo, e há muito mais números compostos do que primos para qualquer intervalo finito dado, essa otimização simples economiza muito tempo.
Isso pode, é claro, ser combinado com uma matriz de primos pré-gerados, basta fazer uma peneira de roda, mas em vez de obter os índices de 1, basta usar a matriz booleana como está, então o primeiro loop na segunda função pode ser usado
nb.prange(2, lim)para paralelizar, o que seria ainda mais rápido, mas já há tanto código que não postarei a última otimização.