A grande maioria dos meus posts sobre geração de números aleatórios lidava principalmente com as propriedades de vários esquemas de geração. Isso pode ser inesperado, mas o desempenho do algoritmo de randomização pode depender não do esquema de geração escolhido, mas de outros fatores. Neste post (inspirado em um excelente
artigo de Daniel Lemyr ), examinaremos as principais razões para o declínio no desempenho da geração de números aleatórios, que geralmente superam o desempenho do mecanismo PRN.
Imagine esta situação:
Como lição de casa, Juan e Sasha implementam o mesmo algoritmo aleatório em C ++, que será executado no mesmo computador da universidade e com um conjunto de dados. Seu código é quase idêntico e difere apenas na geração de números aleatórios. Juan está com pressa para suas aulas de música, então ele simplesmente escolheu o turbilhão de Mersenne. Sasha, por outro lado, passou algumas horas extras pesquisando. Sasha conduziu benchmarks de vários dos PRNGs mais rápidos, que ele aprendeu recentemente nas redes sociais, e escolheu o mais rápido. Na reunião, Sasha estava impaciente em se gabar e perguntou a Juan: "Que sistema PRNG você usou?"
"Pessoalmente, peguei o vórtice de Mersenne - ele está embutido na linguagem e parece funcionar muito bem."
"Ha!" Respondeu Sasha. “Eu usei o
jsf32 . É muito mais rápido que o velho e lento turbilhão de Mersenne! Meu programa é executado em 3 minutos e 15 segundos! ”
"Hmm, nada mal, mas o meu pode fazer isso em menos de um minuto", diz Juan e encolhe os ombros. “Bem, então eu tenho que ir ao show. Você vem comigo?
"Não", responde Sasha. "Eu ... uh ... preciso olhar meu código novamente."
Essa situação ficcional embaraçosa
não é particularmente ficcional; é baseado em resultados reais. Se o seu algoritmo aleatório não está funcionando tão rápido quanto gostaríamos, e o gargalo parece ser a geração de números aleatórios, por incrível que pareça, o problema pode não estar no gerador de números aleatórios!
Introdução: Números Aleatórios na Prática
A maioria dos geradores de números aleatórios de alta qualidade modernos cria palavras de máquina preenchidas com bits aleatórios, ou seja, geralmente geram números no intervalo [0..2
32 ) ou [0..2
64 ). Mas em muitos casos de uso, os usuários precisam de números em um determinado intervalo - por exemplo, para rolar um dado ou escolher uma carta aleatória, números são necessários em pequenos intervalos constantes. No entanto, muitos algoritmos, desde a
mistura e
amostragem de reservatórios até
as árvores de busca binária aleatória, exigem números extraídos de outros intervalos.
Métodos
Veremos muitos métodos diferentes. Para simplificar a discussão, em vez de gerar números no intervalo [
i ..
j ) ou [
i ..
j ], geraremos números no intervalo [0 ..
k ). Tendo esse esquema, podemos, por exemplo, gerar números no intervalo [
i ..
j ) configurando
k =
j -
i , gerando um número no intervalo [0 ..
k ) e adicionando
i a ele.
Ferramentas C ++ incorporadas
Muitos idiomas possuem ferramentas internas para obter um número aleatório em um intervalo especificado. Por exemplo, para remover um cartão de um baralho com 52 cartões em linguagens de script como Perl e Python, podemos escrever
int(rand(52)) e
random.randint(0,52) . [Nota Usuário do
CryptoPirate :
Parece-me um erro aqui, no Python, o randint (a, b) gera números de a para b, incluindo b. E como existem 52 cartas no baralho e a primeira é “0”, deve haver random.randint (0,51) .] No C ++, podemos usar
uniform_int_distribution mesma
uniform_int_distribution .
O código C ++ para implementar essa abordagem é simples:
uint32_t bounded_rand(rng_t& rng, uint32_t range) { std::uniform_int_distribution<uint32_t> dist(0, range-1); return dist(rng); }
Geralmente, uma das técnicas descritas abaixo é usada nas ferramentas internas, mas a maioria dos usuários simplesmente as utiliza, sem pensar no que está acontecendo "oculto", acreditando que essas ferramentas foram projetadas corretamente e são bastante eficazes. No C ++, as ferramentas internas são mais complexas porque devem poder funcionar com mecanismos de geração bastante arbitrários - um gerador que produz valores no intervalo de -3 a 17 pode ser bastante válido e pode ser usado com
std::uniform_int_distribution para criar números em qualquer intervalo, por exemplo [0..1000). Ou seja, as ferramentas C ++ internas são muito complicadas para a maioria dos casos em que são usadas.
O restante clássico da divisão (distorcido)
Vamos passar de uma abordagem simplificada para uma simplista demais.
Quando estudei programação, geramos números no intervalo (por exemplo, para selecionar uma carta em um baralho de 52 cartas) usando o operador restante. Para obter o número no intervalo [0..52), escrevemos
rand() % 52 .
No C ++, essa abordagem pode ser implementada da seguinte maneira:
uint32_t bounded_rand(rng_t& rng, uint32_t range) { return rng() % range; }
Apesar da simplicidade dessa abordagem, ela demonstra o motivo pelo qual obter números no intervalo certo geralmente é uma tarefa lenta - requer divisão (para calcular o restante obtido pelo operador
% ). A divisão é geralmente pelo menos uma ordem de magnitude mais lenta que outras operações aritméticas; portanto, uma única operação aritmética leva mais tempo do que todo o trabalho realizado por um PRNG rápido.
Mas, além da baixa velocidade, também é
distorcida . Para entender por que
rand() % 52 retorna números assimétricos, suponha que
rand() crie números no intervalo [0..2
32 ) e observe que 52 não divide 2
32 completamente, divide-o 82 595 524 vezes com o restante 48. Ou seja, se usarmos
rand() % 52 , teremos 82 595 525 maneiras de selecionar as primeiras 48 cartas do baralho e apenas 82 595 524 maneiras de selecionar as últimas quatro cartas. Em outras palavras, há uma inclinação de 0,00000121% contra essas últimas quatro cartas (talvez sejam reis!). Quando eu era estudante e escrevia a lição de casa sobre jogar dados ou desenhar cartas, ninguém costumava se incomodar com distorções tão pequenas, mas com um aumento no intervalo, a distorção aumenta linearmente. Para um PRNG de 32 bits, um intervalo limitado de menos de
24 tem uma inclinação inferior a 0,5%, mas acima de 2
31 uma inclinação de 50% - alguns números retornam duas vezes mais que outros.
Neste artigo, consideraremos principalmente técnicas que usam estratégias para eliminar um erro sistemático, mas provavelmente vale a pena dizer que, para um PRNG de 64 bits, o valor de inclinação em aplicativos normais provavelmente será desprezível.
Outro problema pode ser que alguns geradores tenham bits baixos fracos. Por exemplo, as famílias GPRS Xoroshiro + e Xoshiro + têm bits baixos que não passam nos testes estatísticos. Quando executamos
% 52 (porque 52 é par), passamos o bit de ordem baixa diretamente para a saída.
Multiplicar números de ponto flutuante (inclinado)
Outra técnica comum é o uso de um PRNG que gera números de ponto flutuante no intervalo [0..1) com a conversão subsequente desses números no intervalo desejado. Essa abordagem é usada no Perl; é
recomendável usar
int(rand(10)) para gerar um número inteiro no intervalo [0..10), gerando um número de ponto flutuante seguido de arredondamento para baixo.
Em C ++, essa abordagem é escrita assim:
static uint32_t bounded_rand(rng_t& rng, uint32_t range) { double zeroone = 0x1.0p-32 * rng(); return range * zeroone; }
(Observe que
0x1.0p-32 é uma constante de ponto flutuante binário para 2
-32 , que usamos para converter um número inteiro aleatório no intervalo [0..2
32 ) para dobrar no intervalo da unidade; em vez disso, podemos realizar essa conversão usando
ldexp(rng(), -32) , mas quando
ldexp(rng(), -32) essa abordagem, ela ficou muito mais lenta.)
Essa abordagem é tão distorcida quanto o restante clássico da divisão, mas a distorção aparece de maneira diferente. Por exemplo, se selecionássemos números no intervalo [0..52), os números 0, 13, 26 e 39 ocorreriam uma vez menos que os outros.
Essa versão, ao generalizar para 64 bits, é ainda mais desagradável, pois requer um tipo de ponto flutuante cuja mantissa é de pelo menos 64 bits. Em máquinas x86 com Linux e macOS, podemos usar
long double para aproveitar os números de ponto flutuante x86 de precisão aumentada que possuem uma mantissa de 64 bits, mas
long double não
long double universalmente portado para todos os sistemas - em alguns sistemas,
long double equivalente a
double .
Existe um lado bom - essa abordagem é mais rápida que as soluções residuais para PRNGs com bits baixos fracos.
Multiplicação de números inteiros (enviesada)
O método de multiplicação pode ser adaptado à aritmética de ponto fixo, e não de ponto flutuante. De fato, apenas multiplicamos constantemente por 2
32 ,
uint32_t bounded_rand(rng_t& rng, uint32_t range) { uint32_t x = rng(); uint64_t m = uint64_t(x) * uint64_t(range); return m >> 32; }
Pode parecer que esta versão requer aritmética de 64 bits. Nos processadores x86, um bom compilador compila esse código em uma instrução
mult 32 bits (que nos fornece dois valores de saída de 32 bits, um dos quais é o valor de retorno). Pode-se esperar que esta versão seja rápida, mas é distorcida exatamente como o método de multiplicação de números de ponto flutuante.
Divisão de queda (sem inclinação)
Podemos modificar o esquema de multiplicação de ponto flutuante em um esquema baseado em divisão. Em vez de multiplicar
x * range / 2**32 calculamos
x / (2**32 / range) . Como trabalhamos com aritmética inteira, o arredondamento nesta versão será realizado de maneira diferente e, às vezes, gera valores fora do intervalo desejado. Se descartamos esses valores (por exemplo, nos livramos deles e geramos novos valores), como resultado, obtemos uma técnica sem distorções.
Por exemplo, no caso de retirar um cartão usando um PRNG de 32 bits, podemos gerar um número de 32 bits e dividi-lo por 2 32/52 = 82 595 524 para selecionar um cartão. Essa técnica funciona se o valor aleatório do PRNG de 32 bits for menor que 52 × 82595524 = 2 32/32 - 48. Se o valor aleatório do PRNR for um dos últimos 48 valores da parte superior do intervalo do gerador, será necessário descartá-lo e procurar outro.
Nosso código para esta versão usa um truque para dividir 2
32 por
range sem usar matemática de 64 bits. Para o cálculo direto de
2**32 / range precisamos representar o número 2
32 , que é muito grande (por um!) Para representar como um número inteiro de 32 bits. Em vez disso, levamos em consideração que, para números inteiros não assinados, o
range operação de negação unária calcula um valor positivo de 2
32 -
range ; Ao dividir esse valor por
range , obtemos uma resposta menor que
2**32 / range .
Portanto, o código C ++ para gerar números usando divisão e soltar se parece com o seguinte:
uint32_t bounded_rand(rng_t& rng, uint32_t range) {
Obviamente, essa abordagem requer duas operações lentas baseadas na divisão, que geralmente são mais lentas que outras operações aritméticas, portanto, você não deve esperar que seja rápida.
O restante da divisão (dupla) sem distorções - técnica OpenBSD
Também podemos adotar a abordagem de queda para eliminar a inclinação no método restante da divisão clássica. No exemplo com cartas de baralho, precisamos novamente soltar 48 valores. Nesta versão, em vez de descartar os
últimos 48 valores, (equivalente) descartamos os
primeiros 48 valores.
Aqui está a implementação dessa abordagem em C ++:
uint32_t bounded_rand(rng_t& rng, uint32_t range) {
Essa técnica remove a inclinação, mas requer duas operações de divisão demoradas com o restante de cada valor de saída (e você pode precisar de um gerador interno para criar vários números). Portanto, deve-se esperar que o método seja aproximadamente duas vezes mais lento que a abordagem de inclinação clássica.
arc4random_uniform OpenBSD (que também é usado no OS X e iOS) usa essa estratégia.
Restante de divisão (única) sem inclinação - metodologia Java
Java usa uma abordagem diferente para gerar um número em um intervalo que usa apenas uma operação de divisão restante, com exceção de casos bastante raros de descartar o resultado. Código:
static uint32_t bounded_rand(rng_t& rng, uint32_t range) { uint32_t x, r; do { x = rng(); r = x % range; } while (x - r > (-range)); return r; }
Para entender por que essa opção funciona, você precisa pensar um pouco. Diferentemente da versão anterior baseada em resíduos, que elimina o viés ao remover parte dos valores mais baixos do mecanismo de geração interno, esta versão filtra os valores da parte superior do intervalo do mecanismo.
Multiplicação de números inteiros inclinados - método Lemira
Da mesma maneira que removemos o viés do método restante da divisão, podemos eliminar o viés da técnica de multiplicação de números inteiros. Esta técnica foi inventada por Lemyr.
uint32_t bounded_rand(rng_t& rng, uint32_t range) { uint32_t t = (-range) % range; do { uint32_t x = rng(); uint64_t m = uint64_t(x) * uint64_t(range); uint32_t l = uint32_t(m); } while (l < t); return m >> 32; }
Drop bitmask (no skew) - técnica da Apple
Em nossa última abordagem, as operações de divisão e restante são completamente eliminadas. Em vez disso, utiliza uma operação de mascaramento simples para obter um número aleatório no intervalo [0..2
k ), onde
k é o menor valor, de modo que 2
k é maior que o intervalo. Se o valor for muito grande para o nosso intervalo, nós o descartamos e tentamos obter outro. O código é mostrado abaixo:
uint32_t bounded_rand(rng_t& rng, uint32_t range) { uint32_t mask = ~uint32_t(0); --range; mask >>= __builtin_clz(range|1); uint32_t x; do { x = rng() & mask; } while (x > range); return x; }
Essa abordagem foi adotada pela Apple quando (na versão macOS Sierra) executou
sua própria revisão do código
arc4random_uniform .
Técnicas básicas de benchmarking
Agora, temos várias abordagens que podem ser avaliadas. Infelizmente, quando estamos preocupados com os custos de uma operação de divisão única, o benchmarking se torna algo não trivial. Nenhuma referência pode levar em consideração todos os fatores que afetam o campo de aplicação, e não há garantia de que a melhor opção para sua aplicação seja certamente a melhor para a minha.
Usamos três parâmetros de referência e testamos as técnicas com muitos PRNGs diferentes.
Aleatório grande de referência
Provavelmente a referência mais óbvia é a mistura. Nesta referência, simulamos a execução de mixagem em larga escala. Para classificar uma matriz de tamanho
N, devemos gerar números nos intervalos [0 ..
N ), [0 .. (
N -1)), ..., [0..1). Nesta referência, assumiremos que
N é o número máximo possível (para
uint32_t é 2
32 -1). Código:
for (uint32_t i = 0xffffffff; i > 0; --i) { uint32_t bval = bounded_rand(rng, i); assert(bval < i); sum += bval; }
Observe que “usamos” cada número adicionando-o à
sum (para que não seja descartado pela otimização), mas não realizamos nenhuma mistura para focar na geração de números.
Para testar a geração de 64 bits, temos um teste semelhante, mas será impraticável executar um teste correspondente à mistura de uma matriz de tamanho 2
64 - 1 (porque levará muitos milhares de anos para concluir essa referência maior). Em vez disso, cruzamos o intervalo inteiro de 64 bits, mas geramos o mesmo número de valores de saída que no teste de 32 bits. Código:
for (uint32_t i = 0xffffffff; i > 0; --i) { uint64_t bound = (uint64_t(i)<<32) | i; uint64_t bval = bounded_rand(rng, bound ); assert(bval < bound); sum += bval; }
Resultados Mersenne vortex
Os resultados mostrados abaixo demonstram o desempenho desse parâmetro de referência para cada um dos métodos que examinamos ao usar o vórtice Mersenne e testá-lo em 32 bits (usando
std::mt19937 da
libstdc++ ) e código semelhante de 64 bits (usando
std:mt19937_64 da
libstdc++ ) Os resultados são a média geométrica de 15 execuções com diferentes valores de sementes, que são normalizadas para que o método restante da divisão clássica tenha um único tempo de execução.
Pode parecer que temos respostas claras sobre desempenho - parece que você pode criar técnicas para a perfeição deles e se perguntar o que os desenvolvedores do
libstdc++ estavam pensando quando escreveram uma implementação tão terrível para números de 32 bits. Mas, como costuma ser o caso do benchmarking, a situação é mais complicada do que parece com esses resultados. Em primeiro lugar, existe o risco de os resultados serem específicos do vórtice de Mersenne; portanto, expandiremos os muitos PRNGs testados. Em segundo lugar, pode haver um problema sutil com o próprio benchmark. Vamos primeiro lidar com a primeira pergunta.
Resultados de diferentes PRNGs
chacha8r gjrand32 32 bits com
arc4_rand32 ,
chacha8r ,
gjrand32 ,
jsf32 ,
mt19937 ,
pcg32 ,
pcg32_fast ,
sfc32 ,
splitmix32 ,
xoroshiro64+ ,
xorshift*64/32 xoshiro128+ ,
xoshiro128+ e
xoshiro128** e
gjrand64 64 bits
jsf64 ,
mcg128 ,
mcg128_fast ,
mt19937_64 ,
pcg64 ,
pcg64_fast ,
sfc64 ,
splitmix64 ,
xoroshiro128+ ,
xorshift*128/64 xoshiro256+ ,
xoshiro256+ e
xoshiro256* . Esses kits nos fornecerão alguns PRNs lentos e muitos muito rápidos.
Aqui estão os resultados:
Podemos ver as principais diferenças dos resultados com o vórtice de Mersenne. PRNGs mais rápidos deslocam o equilíbrio em direção ao código delimitador e, portanto, a diferença entre as diferentes abordagens se torna mais acentuada, especialmente no caso de PRNRs de 64 bits. Com um conjunto mais amplo de
libstc++ implementação
libstc++ deixa de parecer tão terrível.
Conclusões
Nesse benchmark por uma margem significativa, a abordagem baseada na multiplicação com viés ganha em velocidade. Existem muitas situações em que os limites serão pequenos em relação ao tamanho do PRNG, e o desempenho é absolutamente crítico. Nessas situações, é improvável que um leve viés tenha um efeito perceptível, mas a velocidade do PRNG terá. Um exemplo é o Quicksort com um ponto de referência aleatório. Dos métodos distorcidos, a técnica de máscara de bit parece promissora.
Mas antes de tirar conclusões sérias, precisamos apontar o enorme problema desse benchmark - a maior parte do tempo é gasta em limites muito altos, o que provavelmente dá importância excessiva a grandes intervalos. Portanto, precisamos ir para o segundo benchmark.
Aleatório pequeno de referência
, « » (). :
for (uint32_t j = 0; j < 0xffff; ++j) { for (uint32_t i = 0xffff; i > 0; --i) { uint32_t bval = bounded_rand(rng, i); assert(bval < i); sum += bval; } }
Conclusões
, .
, ; , , .
for (uint32_t bit = 1; bit != 0; bit <<= 1) { for (uint32_t i = 0; i < 0x1000000; ++i) { uint32_t bound = bit | (i & (bit - 1)); uint32_t bval = bounded_rand(rng, bound); assert(bval < bound); sum += bval; } }
Conclusões
. , , , , .
, , .
, - , . , .
, :
uint32_t bounded_rand(rng_t& rng, uint32_t range) {
range ,
. , , .
:
uint32_t bounded_rand(rng_t& rng, uint32_t range) { uint32_t r = rng(); if (r < range) { uint32_t t = (-range) % range; while (r < t) r = rng(); } return r % range; }
« Mod » (. ), « ». , ( ).
Large-Shuffle
64- ( mod ), 32- . , .
Small-Shuffle
, small-shuffle , . . (OpenBSD) (Java).
.
, : .
Normalmente, um cálculo a % brequer divisão, mas em situações em que o a < bresultado é simples ae a divisão não é necessária. E quando a/2 < b, o resultado é simples a - b. Portanto, em vez de computar a %= b;
nós podemos cumprir if (a >= b) { a -= b; if (a >= b) a %= b; }
O custo da divisão é tão significativo que o aumento do custo desse código mais complexo pode justificar-se economizando tempo devido à falta de divisão.Resultados de referência de reprodução aleatória grande
A adição dessa otimização aprimora muito os resultados do benchmark de shuffle grande. Isso é mais perceptível no código de 64 bits, onde a operação de tirar o restante é mais cara. O método de duplo restante (estilo OpenBSD) mostra versões com otimizações para apenas uma operação restante e para ambas.- , .
Small-Shuffle
small-shuffle, , . , .
.
:
, . .
32-
32- , , 32- :
, ,
pcg32_fast — Xoroshiro ( ). , - — . , 5%, , «».
64-
O gráfico mostra o desempenho de vários esquemas de geração de 64 bits, em média, entre todas as técnicas e quinze execuções normalizadas para o desempenho do vórtice Mersenne de 32 bits. Pode parecer estranho que a normalização seja realizada usando o vórtice Mersenne de 32 bits, mas isso nos permite ver os custos adicionais do uso da geração de 64 bits nos casos em que a geração de 32 bits é suficiente.,
mcg128_fast , 5%, .
pcg64 pcg64_fast mcg128_fast , 128- (, LCG) 128- (, MCG). , ,
pcg64 20% 64- .
, , 64- , 64- , 32-.
Conclusões
, (, 32- ) 45%. 66%; , .
( ) — ( ). :
uint32_t bounded_rand(rng_t& rng, uint32_t range) { uint32_t x = rng(); uint64_t m = uint64_t(x) * uint64_t(range); uint32_t l = uint32_t(m); if (l < range) { uint32_t t = -range; if (t >= range) { t -= range; if (t >= range) t %= range; } while (l < t) { x = rng(); m = uint64_t(x) * uint64_t(range); l = uint32_t(m); } } return m >> 32; }
, .
:
GitHub . 23
bounded_rand 26 (13 32- 13 64-), (GCC 8 LLVM 6), 26 * 23 * 2 = 1196 , 15 seed, 1196 * 15 = 17 940 , . 48- Xeon E7-4830v3 2,1 . .
. ,
jsf32.STD-libc++ , —
mt19937.BIASED_FP_MULT_SCALE. No benchmark 3, este último leva 69,6% menos tempo. Ou seja, o tempo dessa situação ficcional é baseado em dados da realidade.