1.0 é uma saída válida de std :: generate_canonical?

Eu sempre pensei que números aleatórios ficariam entre zero e um, sem1 , ou seja, são números do intervalo semiaberto [0,1). A documentação em cppreference.com de std::generate_canonicalconfirma isso.

No entanto, quando executo o seguinte programa:

#include <iostream>
#include <limits>
#include <random>

int main()
{
    std::mt19937 rng;

    std::seed_seq sequence{0, 1, 2, 3, 4, 5, 6, 7, 8, 9};
    rng.seed(sequence);
    rng.discard(12 * 629143 + 6);

    float random = std::generate_canonical<float,
                   std::numeric_limits<float>::digits>(rng);

    if (random == 1.0f)
    {
        std::cout << "Bug!\n";
    }

    return 0;
}

Dá-me a seguinte saída:

Bug!

ou seja, ele me gera um perfeito 1, o que causa problemas na minha integração com o MC. Esse comportamento é válido ou existe um erro do meu lado? Isso fornece a mesma saída com o G ++ 4.7.3

g++ -std=c++11 test.c && ./a.out

e clang 3.3

clang++ -stdlib=libc++ -std=c++11 test.c && ./a.out

Se esse é o comportamento correto, como posso evitar 1?

Edit 1 : G ++ do git parece sofrer do mesmo problema. Estou em

commit baf369d7a57fb4d0d5897b02549c3517bb8800fd
Date:   Mon Sep 1 08:26:51 2014 +0000

e compilar com ~/temp/prefix/bin/c++ -std=c++11 -Wl,-rpath,/home/cschwan/temp/prefix/lib64 test.c && ./a.outfornece a mesma saída, lddgera

linux-vdso.so.1 (0x00007fff39d0d000)
libstdc++.so.6 => /home/cschwan/temp/prefix/lib64/libstdc++.so.6 (0x00007f123d785000)
libm.so.6 => /lib64/libm.so.6 (0x000000317ea00000)
libgcc_s.so.1 => /home/cschwan/temp/prefix/lib64/libgcc_s.so.1 (0x00007f123d54e000)
libc.so.6 => /lib64/libc.so.6 (0x000000317e600000)
/lib64/ld-linux-x86-64.so.2 (0x000000317e200000)

Edição 2 : relatei o comportamento aqui: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63176

Edit 3 : A equipe do clang parece estar ciente do problema: http://llvm.org/bugs/show_bug.cgi?id=18767

O problema está no mapeamento do codomain de std::mt19937( std::uint_fast32_t) para float; o algoritmo descrito pelo padrão fornece resultados incorretos (inconsistente com a descrição da saída do algoritmo) quando ocorre perda de precisão, se o modo de arredondamento IEEE754 atual for diferente de infinito arredondado para negativo (observe que o padrão é redondo -para-mais próximo).

A saída 7549723rd do mt19937 com sua semente é 4294967257 ( 0xffffffd9u), que quando arredondada para float de 32 bits fornece 0x1p+32, que é igual ao valor máximo de mt19937, 4294967295 ( 0xffffffffu) quando também é arredondada para float de 32 bits.

O padrão poderia garantir o comportamento correto se especificasse que, ao converter da saída do URNG para o RealTypede generate_canonical, o arredondamento deve ser realizado em direção ao infinito negativo; isso daria um resultado correto neste caso. Como QOI, seria bom para o libstdc ++ fazer essa alteração.

Com essa alteração, 1.0não será mais gerado; em vez disso, os valores limite 0x1.fffffep-Npara 0 < N <= 8serão gerados com mais frequência (aproximadamente 2^(8 - N - 32)por N, dependendo da distribuição real do MT19937).

Eu recomendaria não usar floatcom std::generate_canonicaldiretamente; em vez disso, gere o número doublee depois arredonde para o infinito negativo:

    double rd = std::generate_canonical<double,
        std::numeric_limits<float>::digits>(rng);
    float rf = rd;
    if (rf > rd) {
      rf = std::nextafter(rf, -std::numeric_limits<float>::infinity());
    }

Esse problema também pode ocorrer com std::uniform_real_distribution<float>; a solução é a mesma, para especializar a distribuição doublee arredondar o resultado para o infinito negativo em float.

De acordo com o padrão, 1.0não é válido.

C ++ 11 §26.5.7.2 Modelo de função generate_canonical

Cada função instanciado a partir do modelo descrito no presente ponto 26.5.7.2 mapeia o resultado de um ou mais invocações de um gerador de número aleatório uniforme fornecido ga um membro da RealType especificada de tal modo que, se os valores g _i produzido pela gsão uniformemente distribuídas, o os resultados da instanciação t _j , 0 ≤ t _j <1 , são distribuídos da maneira mais uniforme possível, conforme especificado abaixo.

Acabei de ter uma pergunta semelhante com uniform_real_distribution, e aqui está como eu interpreto a redação parcimoniosa do Padrão sobre o assunto:

O Padrão sempre define funções matemáticas em termos de matemática , nunca em termos de ponto flutuante IEEE (porque o Padrão ainda finge que ponto flutuante pode não significar ponto flutuante IEEE). Portanto, sempre que você vê palavras matemáticas no Padrão, trata-se de matemática real , não de IEEE.

A Norma diz que ambos uniform_real_distribution<T>(0,1)(g)e generate_canonical<T,1000>(g)devem retornar valores no intervalo semi-aberto [0,1). Mas esses são valores matemáticos . Quando você pega um número real no intervalo semi-aberto [0,1) e o representa como ponto flutuante IEEE, bem, uma fração significativa do tempo que ele arredondará para T(1.0).

Quando Té float(24 bits de mantissa), esperamos ver uniform_real_distribution<float>(0,1)(g) == 1.0fcerca de 1 em 2 ^ 25 vezes. Minha experimentação de força bruta com libc ++ confirma essa expectativa.

template<class F>
void test(long long N, const F& get_a_float) {
    int count = 0;
    for (long long i = 0; i < N; ++i) {
        float f = get_a_float();
        if (f == 1.0f) {
            ++count;
        }
    }
    printf("Expected %d '1.0' results; got %d in practice\n", (int)(N >> 25), count);
}

int main() {
    std::mt19937 g(std::random_device{}());
    auto N = (1uLL << 29);
    test(N, [&g]() { return std::uniform_real_distribution<float>(0,1)(g); });
    test(N, [&g]() { return std::generate_canonical<float, 32>(g); });
}

Exemplo de saída:

Expected 16 '1.0' results; got 19 in practice
Expected 16 '1.0' results; got 11 in practice

Quando Té double(53 mantissa bits), esperamos ver uniform_real_distribution<double>(0,1)(g) == 1.0cerca de 1 em 2 ^ 54 vezes. Não tenho paciência para testar essa expectativa. :)

Meu entendimento é que esse comportamento é bom. Isso pode ofender nosso senso de "meio intervalo aberto" de que uma distribuição que afirma retornar números "inferiores a 1,0" pode de fato retornar números iguais a 1.0; mas esses são dois significados diferentes de "1.0", entende? O primeiro é o 1.0 matemático ; o segundo é o número de ponto flutuante de precisão única IEEE 1.0. E somos ensinados há décadas a não comparar números de ponto flutuante para obter a igualdade exata.

Qualquer que seja o algoritmo no qual você alimenta os números aleatórios, não se importará se às vezes isso acontecer exatamente 1.0. Não há nada que você possa fazer com um número de ponto flutuante, exceto operações matemáticas, e assim que você fizer alguma operação matemática, seu código terá que lidar com arredondamentos. Mesmo se você pudesse legitimamente assumir isso generate_canonical<float,1000>(g) != 1.0f, ainda não seria capaz de assumir isso generate_canonical<float,1000>(g) + 1.0f != 2.0f- por causa do arredondamento. Você simplesmente não pode fugir disso; Então, por que fingiríamos neste caso único que você pode?