Flottante casuale in C utilizzando getrandom

Question

Sto provando a generare un numero a virgola mobile casuale compreso tra 0 e 1 (se è su [0,1] o [0,1) non dovrebbe importare per me). Ogni domanda in linea su questo sembra coinvolgere la chiamata rand(), seminata con time(NULL), ma voglio essere in grado di invocare il mio programma più di una volta al secondo e ottenere numeri casuali diversi ogni volta. Questo mi ha portato alla syscall getrandom in Linux, che richiama da/dev/urandom. Sono arrivato fino a questo:

#include <stdio.h> 
#include <sys/syscall.h> 
#include <unistd.h> 
#include <stdint.h> 

int main() { 
    uint32_t r = 0; 
    for (int i = 0; i < 20; i++) { 
    syscall(SYS_getrandom, &r, sizeof(uint32_t), 0); 
    printf("%f\n", ((double)r)/UINT32_MAX); 
    } 
    return 0; 
} 

La mia domanda è semplice anche se non sto facendo questo in modo corretto. Sembra funzionare, ma sono preoccupato che io stia abusando di qualcosa, e non ci sono quasi esempi che usano getrandom() online.

Answer 1

OP ha 2 problemi:

1. Come ha iniziato la sequenza molto casualmente.

2. Come generare un double nell'intervallo [0 ... 1).

Il metodo usuale è scattare una fonte molto casuale come /dev/urandom o il risultato dal syscall() o forse anche seed = time()^process_id; e sementi tramite srand(). Quindi chiama rand() secondo necessità.

In basso è incluso un metodo rapido per generare una uniforme [0.0 to 1.0) (distribuzione lineare). Ma come tutte le funzioni generatrici casuali, quelle veramente buone sono basate su studi approfonditi. Questo si chiama semplicemente rand() un paio di volte sulla base di DBL_MANT_DIG e RAND_MAX,

[modifica] Original double rand_01(void) ha una debolezza che genera solo 2^52 diversi double s invece di 2^53. È stato modificato Alternativa: una versione double di rand_01_ld(void) molto al di sotto.

#include <assert.h> 
#include <float.h> 
#include <limits.h> 
#include <stdint.h> 
#include <stdio.h> 
#include <stdlib.h> 
#include <unistd.h> 

double rand_01(void) { 
    assert(FLT_RADIX == 2); // needed for DBL_MANT_DIG 
    unsigned long long limit = (1ull << DBL_MANT_DIG) - 1; 
    double r = 0.0; 
    do { 
    r += rand(); 
    // Assume RAND_MAX is a power-of-2 - 1 
    r /= (RAND_MAX/2 + 1)*2.0; 
    limit = limit/(RAND_MAX/2 + 1)/2; 
    } while (limit); 

    // Use only DBL_MANT_DIG (53) bits of precision. 
    if (r < 0.5) { 
    volatile double sum = 0.5 + r; 
    r = sum - 0.5; 
    } 
    return r; 
} 

int main(void) { 
    FILE *istream = fopen("/dev/urandom", "rb"); 
    assert(istream); 
    unsigned long seed = 0; 
    for (unsigned i = 0; i < sizeof seed; i++) { 
    seed *= (UCHAR_MAX + 1); 
    int ch = fgetc(istream); 
    assert(ch != EOF); 
    seed += (unsigned) ch; 
    } 
    fclose(istream); 
    srand(seed); 

    for (int i=0; i<20; i++) { 
    printf("%f\n", rand_01()); 
    } 

    return 0; 
} 

---

Se si volesse estendere ad una gamma ancora più FP, senza segno tipi interi ampi possono essere insufficienti. Di seguito è riportato un metodo portatile che non ha questa limitazione.

long double rand_01_ld(void) { 
    // These should be calculated once rather than each function call 
    // Leave that as a separate implementation problem 
    // Assume RAND_MAX is power-of-2 - 1 
    assert((RAND_MAX & (RAND_MAX + 1U)) == 0); 
    double rand_max_p1 = (RAND_MAX/2 + 1)*2.0; 
    unsigned BitsPerRand = (unsigned) round(log2(rand_max_p1)); 
    assert(FLT_RADIX != 10); 
    unsigned BitsPerFP = (unsigned) round(log2(FLT_RADIX)*LDBL_MANT_DIG); 

    long double r = 0.0; 
    unsigned i; 
    for (i = BitsPerFP; i >= BitsPerRand; i -= BitsPerRand) { 
    r += rand(); 
    r /= rand_max_p1; 
    } 
    if (i) { 
    r += rand() % (1 << i); 
    r /= 1 << i; 
    } 
    return r; 
} 

Answer 2

Se è necessario generare raddoppia, il seguente algoritmo potrebbe essere utile:

CPython generates random numbers utilizzando il seguente algoritmo (ho cambiato il nome della funzione, typedef e valori di ritorno, ma l'algoritmo rimane la stessa):

double get_random_double() { 
    uint32_t a = get_random_uint32_t() >> 5; 
    uint32_t b = get_random_uint32_t() >> 6; 
    return (a * 67108864.0 + b) * (1.0/9007199254740992.0); 
} 

La fonte di tale algoritmo è una Mersenne Twister generatore 19937 numero casuale per Takuji Nishimura e Makoto Matsumoto. Sfortunatamente il link originale menzionato nella fonte non è più disponibile per il download.

Il commento su questa funzione nel CPython osserva quanto segue:

[questa funzione] è la funzione denominata genrand_res53 nel codice originale; genera un numero casuale su [0,1) con una risoluzione di 53 bit; si noti che 9007199254740992 == 2**53; Presumo che stiano scrivendo "/2**53" come moltiplicando per reciproco nella speranza (probabilmente vana) che il compilatore sarà ottimizzare la divisione in fase di compilazione. 67108864 è 2**26. Nell'effetto , a contiene 27 bit casuali spostati a sinistra 26 e b riempimenti nei 26 bit inferiori del numeratore a 53 bit.

Il codice originale accreditato Wada Isaku per questo algoritmo, 2002/01/09

---

Semplificare da quel codice, se si desidera creare un float veloce, si dovrebbe mascherare i bit di uint32_t con (1 << FLT_MANT_DIG) - 1 e dividere per (1 << FLT_MANT_DIG) per ottenere il corretto intervallo di [0, 1):

#include <stdio.h> 
#include <sys/syscall.h> 
#include <unistd.h> 
#include <stdint.h> 
#include <float.h> 

int main() { 
    uint32_t r = 0; 
    float result; 
    for (int i = 0; i < 20; i++) { 
     syscall(SYS_getrandom, &r, sizeof(uint32_t), 0); 
     result = (float)(r & ((1 << FLT_MANT_DIG) - 1))/(1 << FLT_MANT_DIG); 
     printf("%f\n", result); 
    } 
    return 0; 
} 

Poiché si può presumere che il vostro Linux ha un compilatore C99, possiamo usare ldexpf al posto di quella divisione:

#include <math.h> 

result = ldexpf(r & ((1 << FLT_MANT_DIG) - 1), -FLT_MANT_DIG); 

Per ottenere l'intervallo chiuso [0, 1], si può fare leggermente meno efficiente

result = ldexpf(r % (1 << FLT_MANT_DIG), -FLT_MANT_DIG); 

---

Per generare un sacco di buona qualità numeri casuali veloci, vorrei semplicemente usare la chiamata di sistema per recuperare dati sufficienti per seminare un PRNG o CPRNG e procedere da lì.