¿Es 1.0 una salida válida de std :: generate_canonical?

Siempre pensé que los números aleatorios estarían entre cero y uno, sin ellos1 , es decir, son números del intervalo medio abierto [0,1]. La documentación en cppreference.com de std::generate_canonicalconfirma esto.

Sin embargo, cuando ejecuto el siguiente 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;
}

Me da el siguiente resultado:

Bug!

es decir, me genera un perfecto 1, lo que causa problemas en mi integración MC. ¿Es ese comportamiento válido o hay un error de mi parte? Esto da el mismo resultado con G ++ 4.7.3

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

y clang 3.3

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

Si este es el comportamiento correcto, ¿cómo puedo evitarlo 1?

Edición 1 : G ++ de git parece sufrir el mismo problema. Estoy en

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

y compilando ~/temp/prefix/bin/c++ -std=c++11 -Wl,-rpath,/home/cschwan/temp/prefix/lib64 test.c && ./a.outda la misma salida, lddrendimientos

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)

Edición 2 : informé el comportamiento aquí: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63176

Edición 3 : El equipo clang parece estar al tanto del problema: http://llvm.org/bugs/show_bug.cgi?id=18767

El problema está en el mapeo del codominio de std::mt19937( std::uint_fast32_t) a float; el algoritmo descrito por el estándar da resultados incorrectos (inconsistentes con su descripción de la salida del algoritmo) cuando se produce una pérdida de precisión si el modo de redondeo IEEE754 actual no es redondo a infinito negativo (tenga en cuenta que el valor predeterminado es redondo -al más cercano).

La salida 7549723a de mt19937 con su semilla es 4294967257 ( 0xffffffd9u), que cuando se redondea a flotación de 32 bits 0x1p+32, que es igual al valor máximo de mt19937, 4294967295 ( 0xffffffffu) cuando eso también se redondea a flotación de 32 bits.

El estándar podría garantizar un comportamiento correcto si especificara que al convertir de la salida de la URNG a la RealTypede generate_canonical, el redondeo debe realizarse hacia el infinito negativo; Esto daría un resultado correcto en este caso. Como QOI, sería bueno para libstdc ++ hacer este cambio.

Con este cambio, 1.0ya no se generará; en cambio, los valores límite 0x1.fffffep-Npara 0 < N <= 8se generarán con más frecuencia (aproximadamente 2^(8 - N - 32)por N, dependiendo de la distribución real de MT19937).

Yo recomendaría no usar floatcon std::generate_canonicaldirectamente; en lugar de generar el número en doubley luego redondear hacia el 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());
    }

Este problema también puede ocurrir con std::uniform_real_distribution<float>; la solución es la misma, para especializar la distribución doubley redondear el resultado hacia infinito negativo en float.

Según el estándar, 1.0no es válido.

C ++ 11 §26.5.7.2 Plantilla de función generate_canonical

Cada función instanciada a partir de la plantilla descrita en esta sección 26.5.7.2 asigna el resultado de una o más invocaciones de un generador de números aleatorios uniforme suministrado ga un miembro del RealType especificado de modo que, si los valores producidos por g _ig están distribuidos uniformemente, el Los resultados de instanciación t _j , 0 ≤ t _j <1 , se distribuyen de la manera más uniforme posible como se especifica a continuación.

Me encontré con una pregunta similar con uniform_real_distribution, y así es como interpreto la redacción parsimoniosa de la Norma sobre el tema:

El estándar siempre define las funciones matemáticas en términos de matemáticas , nunca en términos de punto flotante IEEE (porque el estándar todavía finge que punto flotante podría no significar punto flotante IEEE). Entonces, cada vez que ve una redacción matemática en el Estándar, se trata de matemáticas reales , no de IEEE.

La norma dice que tanto uniform_real_distribution<T>(0,1)(g)y generate_canonical<T,1000>(g)debe devolver valores en el rango medio abierta [0,1). Pero estos son valores matemáticos . Cuando toma un número real en el rango medio abierto [0,1) y lo representa como punto flotante IEEE, bueno, una fracción significativa del tiempo se redondeará T(1.0).

Cuando Tes float(24 bits de mantisa), esperamos ver uniform_real_distribution<float>(0,1)(g) == 1.0faproximadamente 1 de cada 2 ^ 25 veces. Mi experimentación de fuerza bruta con libc ++ confirma esta 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); });
}

Salida de ejemplo:

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

Cuando Tes double(53 bits de mantisa), esperamos ver uniform_real_distribution<double>(0,1)(g) == 1.0aproximadamente 1 de cada 2 ^ 54 veces. No tengo la paciencia para probar esta expectativa. :)

Tengo entendido que este comportamiento está bien. Puede ofender nuestro sentido de "rango medio abierto" de que una distribución que dice devolver números "menores que 1.0" puede de hecho devolver números que son iguales a 1.0; pero esos son dos significados diferentes de "1.0", ¿ves? El primero es el matemático 1.0; el segundo es el número de coma flotante de precisión simple IEEE 1.0. Y nos han enseñado durante décadas a no comparar números de punto flotante para la igualdad exacta.

Cualquier algoritmo en el que alimente los números aleatorios no le importará si a veces es exacto 1.0. No hay nada que pueda hacer con un número de punto flotante, excepto las operaciones matemáticas, y tan pronto como realice alguna operación matemática, su código tendrá que lidiar con el redondeo. Incluso si pudieras asumir eso legítimamente generate_canonical<float,1000>(g) != 1.0f, aún no podrías asumir eso generate_canonical<float,1000>(g) + 1.0f != 2.0f, debido al redondeo. Simplemente no puedes alejarte de eso; Entonces, ¿por qué fingiríamos en esta única instancia que puedes?