Est-ce que 1.0 est une sortie valide de std::generate_canonical ?

Question

J'ai toujours pensé que les nombres aléatoires se situaient entre zéro et un, sans 1 c'est-à-dire qu'il s'agit de nombres de l'intervalle semi-ouvert [0,1]. Les documentation sur cppreference.com de std::generate_canonical le confirme.

Cependant, lorsque j'exécute le programme suivant :

#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;
}

Le résultat est le suivant :

Bug!

c'est-à-dire qu'il me génère un 1 ce qui cause des problèmes dans l'intégration de mon MC. Ce comportement est-il valable ou y a-t-il une erreur de ma part ? Cela donne le même résultat avec G++ 4.7.3

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

et clang 3.3

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

S'il s'agit d'un comportement correct, comment puis-je éviter 1 ?

Édition 1 : G++ à partir de git semble souffrir du même problème. Je suis sur

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

et en compilant avec ~/temp/prefix/bin/c++ -std=c++11 -Wl,-rpath,/home/cschwan/temp/prefix/lib64 test.c && ./a.out donne le même résultat, ldd rendements

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)

Édition 2 : J'ai signalé ce comportement ici : https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63176

Édition 3 : L'équipe de clang semble être consciente du problème : http://llvm.org/bugs/show_bug.cgi?id=18767

Answer 1

Le problème est de passer du codomaine de std::mt19937 ( std::uint_fast32_t ) à float l'algorithme décrit par la norme donne des résultats incorrects (incompatibles avec sa description de la sortie de l'algorithme) lorsque la perte de précision se produit si le mode d'arrondi actuel de l'IEEE754 est autre que l'arrondi à l'infini négatif (notez que la valeur par défaut est l'arrondi à l'infini le plus proche).

La 7549723ème sortie de mt19937 avec votre semence est 4294967257 ( 0xffffffd9u ), ce qui, une fois arrondi à un flottant de 32 bits, donne 0x1p+32 qui est égale à la valeur maximale de mt19937, 4294967295 ( 0xffffffffu ) lorsqu'il est également arrondi à un flottant de 32 bits.

La norme pourrait garantir un comportement correct si elle spécifiait que lors de la conversion de la sortie de l'URNG vers le format RealType de generate_canonical Dans ce cas, l'arrondi doit être effectué vers l'infini négatif, ce qui donnerait un résultat correct dans ce cas. En tant que QOI, il serait bon que libstdc++ fasse ce changement.

Avec cette modification, 1.0 ne seront plus générées ; à la place, les valeurs limites 0x1.fffffep-N para 0 < N <= 8 sera généré plus souvent (environ 2^(8 - N - 32) par N en fonction de la distribution effective de MT19937).

Je recommande de ne pas utiliser float con std::generate_canonical directement ; il s'agit plutôt de générer le nombre dans double et s'arrondit ensuite vers l'infini négatif :

    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());
    }

Ce problème peut également se produire avec std::uniform_real_distribution<float> ; la solution est la même, il s'agit de spécialiser la distribution sur double et arrondir le résultat vers l'infini négatif en float .

Answer 2

Selon la norme, 1.0 n'est pas valide.

C++11 §26.5.7.2 Fonction template generate_canonical

Chaque fonction instanciée à partir du modèle décrit dans la présente section 26.5.7.2 mappe le résultat d'une ou plusieurs invocations d'un générateur de nombres aléatoires uniforme fourni g à un membre du type RealType spécifié de telle sorte que, si les valeurs g _i produit par g sont uniformément distribués, les résultats de l'instanciation t _j , 0 t _j < 1 sont réparties aussi uniformément que possible, comme indiqué ci-dessous.

Answer 3

Il semble que GCC et Clang se soient trompés en produisant la version 1.0 et c'est un bogue qui devrait être signalé.

Pour être complet, j'ai essayé le code posté par l'OP avec VS2013 et j'obtiens des nombres qui sont >> 1.

1.879446e+018
1.774417e+019
3.401991e+018
1.844674e+019
2.513807e+018
4.104100e+018
1.147449e+019
5.705059e+018
1.407497e+019
1.018316e+019
1.714579e+018
5.609772e+018

En cherchant un peu, j'ai trouvé qu'il s'agit d'un bug pour VS2013 qui est toujours actif ( rapport de bogue ).

Answer 4

Je viens de rencontrer une question similaire avec uniform_real_distribution Voici comment j'interprète la formulation parcimonieuse de la norme à ce sujet :

La norme définit toujours les fonctions mathématiques en termes de mathématiques mais jamais en termes de virgule flottante IEEE (parce que la norme prétend toujours que la virgule flottante pourrait ne pas moyenne en virgule flottante IEEE). Par conséquent, chaque fois que vous voyez une formulation mathématique dans la norme, il s'agit de les vraies mathématiques et non l'IEEE.

La norme stipule que les deux uniform_real_distribution<T>(0,1)(g) y generate_canonical<T,1000>(g) doit renvoyer des valeurs comprises dans l'intervalle semi-ouvert [0,1]. Mais ces valeurs sont mathématique valeurs. Lorsque vous prenez un nombre réel dans l'intervalle semi-ouvert [0,1] et que vous le représentez en virgule flottante IEEE, une grande partie du temps, il s'arrondira à T(1.0) .

Quand T es float (24 bits de mantisse), nous nous attendons à voir uniform_real_distribution<float>(0,1)(g) == 1.0f environ 1 fois sur 2^25. Mon expérience de force brute avec libc++ confirme cette attente.

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); });
}

Exemple de sortie :

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

Quand T es double (53 bits de mantisse), nous nous attendons à voir uniform_real_distribution<double>(0,1)(g) == 1.0 environ 1 fois sur 2^54. Je n'ai pas la patience de tester cette attente :)

Je pense que ce comportement est acceptable. Le fait qu'une distribution prétendant renvoyer des nombres "inférieurs à 1,0" puisse en fait renvoyer des nombres qui sont "inférieurs à 1,0" peut heurter notre sens de l'"étrangeté semi-ouverte". égal a 1.0 ; mais il s'agit de deux significations différentes de "1.0", vous voyez ? La première est la mathématique 1,0 ; le second est le nombre IEEE à virgule flottante en simple précision 1.0 . Et on nous a appris pendant des décennies à ne pas comparer des nombres en virgule flottante pour obtenir une égalité exacte.

Quel que soit l'algorithme dans lequel vous introduisez les nombres aléatoires, il ne se souciera pas de savoir s'il obtient parfois exactement 1.0 . Il n'y a rien que vous puissiez do avec un nombre à virgule flottante, à l'exception des opérations mathématiques, et dès que vous effectuez une opération mathématique, votre code doit tenir compte des arrondis. Même si vous pourrait supposer légitimement que generate_canonical<float,1000>(g) != 1.0f , vous encore ne serait pas en mesure de supposer que generate_canonical<float,1000>(g) + 1.0f != 2.0f - en raison des arrondis. On ne peut tout simplement pas y échapper ; alors pourquoi prétendre, dans ce cas précis, que c'est possible ?