La 1.0 est-elle une sortie valide de std :: generate_canonical?

J'ai toujours pensé que les nombres aléatoires se situeraient entre zéro et un, sans1 , c'est-à-dire que ce sont des nombres de l'intervalle semi-ouvert [0,1). La documentation sur cppreference.com de le std::generate_canonicalconfirme.

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

Cela me donne le résultat suivant:

Bug!

c'est à dire qu'il me génère un parfait 1, ce qui pose des problèmes dans mon intégration MC. Est-ce un comportement valable ou y a-t-il une erreur de mon côté? Cela donne la même sortie 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

Si c'est un comportement correct, comment puis-je éviter 1?

Edit 1 : G ++ de git semble souffrir du même problème. Je suis sur

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

et compiler avec ~/temp/prefix/bin/c++ -std=c++11 -Wl,-rpath,/home/cschwan/temp/prefix/lib64 test.c && ./a.outdonne le même résultat, ldddonne

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)

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

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

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

La sortie 7549723e de mt19937 avec votre valeur de départ est 4294967257 ( 0xffffffd9u), qui, lorsqu'elle est arrondie à un flottant de 32 bits 0x1p+32, donne , ce qui est égal à la valeur maximale de mt19937, 4294967295 ( 0xffffffffu) lorsqu'elle est également arrondie à 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 RealTypede generate_canonical, l'arrondi doit être effectué vers l'infini négatif; cela donnerait un résultat correct dans ce cas. En tant que QOI, il serait bon que libstdc ++ fasse ce changement.

Avec ce changement, 1.0ne sera plus généré; au lieu de cela, les valeurs limites 0x1.fffffep-Npour 0 < N <= 8seront générées plus souvent (approximativement 2^(8 - N - 32)par N, selon la distribution réelle de MT19937).

Je recommanderais de ne pas utiliser directement floatavec std::generate_canonical; plutôt générer le nombre dans doublepuis arrondir 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, pour spécialiser la distribution sur doubleet arrondir le résultat vers l'infini négatif en float.

Selon la norme, 1.0n'est pas valide.

C ++ 11 §26.5.7.2 Modèle de fonction generate_canonical

Chaque fonction instanciée à partir du modèle décrit dans cette 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 RealType spécifié de telle sorte que, si les valeurs g _i produites par gsont uniformément distribuées, le Les résultats d'instanciation t _j , 0 ≤ t _j <1 , sont distribués aussi uniformément que possible comme spécifié ci-dessous.

Je viens de rencontrer une question similaire avec uniform_real_distribution, et voici comment j'interprète le libellé parcimonieux du Standard sur le sujet:

La norme définit toujours les fonctions mathématiques en termes de mathématiques , jamais en termes de virgule flottante IEEE (car la norme prétend toujours que virgule flottante peut ne pas signifier virgule flottante IEEE). Ainsi, chaque fois que vous voyez un libellé mathématique dans la norme, il s'agit de vraies mathématiques , pas d'IEEE.

La norme dit que les deux uniform_real_distribution<T>(0,1)(g)et generate_canonical<T,1000>(g)devraient renvoyer des valeurs dans la plage semi-ouverte [0,1). Mais ce sont des valeurs mathématiques . Lorsque vous prenez un nombre réel dans la plage semi-ouverte [0,1) et que vous le représentez en virgule flottante IEEE, eh bien, une fraction significative du temps à laquelle il sera arrondi T(1.0).

Quand Test float(24 bits de mantisse), nous nous attendons à voir uniform_real_distribution<float>(0,1)(g) == 1.0fenviron 1 fois sur 2 ^ 25. Mon expérimentation 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 Test double(53 bits de mantisse), nous nous attendons à voir uniform_real_distribution<double>(0,1)(g) == 1.0environ 1 fois sur 2 ^ 54. Je n'ai pas la patience de tester cette attente. :)

Je crois comprendre que ce comportement est bien. Cela peut heurter notre sentiment de "demi-gamme ouverte" qu'une distribution prétendant renvoyer des nombres "inférieurs à 1.0" puisse en fait renvoyer des nombres égaux à 1.0; mais ce sont deux significations différentes de "1.0", vous voyez? Le premier est le 1.0 mathématique ; le second est le nombre à virgule flottante simple précision IEEE 1.0. Et on nous a appris pendant des décennies à ne pas comparer les nombres à virgule flottante pour une égalité exacte.

Quel que soit l'algorithme dans lequel vous introduisez les nombres aléatoires, cela ne se souciera pas de savoir s'il obtient parfois exactement 1.0. Il n'y a rien que vous puissiez faire avec un nombre à virgule flottante sauf des opérations mathématiques, et dès que vous effectuez une opération mathématique, votre code devra gérer l'arrondi. Même si vous pouviez légitimement supposer cela generate_canonical<float,1000>(g) != 1.0f, vous ne seriez toujours pas en mesure de le supposer generate_canonical<float,1000>(g) + 1.0f != 2.0f- à cause de l'arrondissement. Vous ne pouvez tout simplement pas vous en éloigner; alors pourquoi ferions-nous semblant dans ce cas unique que vous le pouvez?