Newsgroup: it.comp.lang.c
Titolo: FAQ 7 - Numeri pseudocasuali

Per oggi è l'ultima ;-) !
-------------------------------------------------------------
Come faccio a generare un numero (pseudo) casuale ?
Il C dispone della funzione rand() definita in stdlib.h per la generazione
di numeri pseudocasuali nell'intervallo compreso tra o e RAND_MAX (costante
definita in stdlib.h).
Essi sono valori ottenuti per via algoritmica, per cui in realtà casuali non
lo sono affatto, difatti se non si inizializza l'algoritmo con un opportuno
"seme" la sequenza di valori risulta essere sempre la stessa. La funzione
utilizzata allo scopo è la srand(seme) anch'essa definita in stdlib.h.
Perché sia efficace, il seme deve poter variare. Si usa perciò chiamarla con
il valore restituito dalla funzione time(NULL) definita in time.h. Essa
ritorna il numero di secondi trascorsi dal 1° gennaio 1970 ed è perciò
differente ogni volta che la si chiama. Vediamo ad esempio come generare un
numero compreso tra 0 e 9: NB: rand() ritorna un int
#include <stdlib.h>
#include <time.h>
int main(void) {
int numerocasuale ;
srand(time(NULL)) ;
numerocasuale=(rand() %10);
...}
--
Marco Coppo
University Of Sheffield

 ...}
Non sono sicuro che questo esempio sia portabile per sistemi dove un intero e'
rappresentato da 2 byte
Ivan Molella

 when you write.
******************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
/**
 | uint32 must be an unsigned integer type capable of holding at least
32
 | bits; exactly 32 should be fastest, but 64 is better on an Alpha with
 | GCC at -O3 optimization so try your options and see what's best for
you
**/
typedef unsigned long uint32;
#define N            (624)                /* length of state
vector          */
#define M            (397)                /* a period
parameter              */
#define K            (0x9908B0DFU)        /* a magic
constant                */
#define hiBit(u)     ((u) & 0x80000000U)  /* mask all but highest   bit
of u */
#define loBit(u)     ((u) & 0x00000001U)  /* mask all but lowest    bit
of u */
#define loBits(u)    ((u) & 0x7FFFFFFFU)  /* mask     the highest   bit
of u */
#define mixBits(u,v) (hiBit(u)|loBits(v)) /* move hi bit of u to hi bit
of v */
#define UINT32MAX    (0xFFFFFFFF)         /* the biggest
value               */
#define TWOM32 (2.3283064365386963E-10)   /* (double) 2^(-32) */
static uint32  state[N+1];   /* state vector + 1 extra to not violate
ANSI C */
static uint32 *next;         /* next random value is computed from
here      */
static int     left = -1;    /* can *next++ this many times before
reloading */
void seedMT(
  uint32 seed
) {
/**
 |
 |  We initialize state[0..(N-1)] via the generator
 |
 |    x_new = (69069 * x_old) mod 2^32
 |
 |  from Line 15 of Table 1, p. 106, Sec. 3.3.4 of Knuth's
 |  _The Art of Computer Programming_, Volume 2, 3rd ed.
 |
 |  Notes (SJC): I do not know what the initial state requirements
 |  of the Mersenne Twister are, but it seems this seeding generator
 |  could be better.  It achieves the maximum period for its modulus
 |  (2^30) iff x_initial is odd (p. 20-21, Sec. 3.2.1.2, Knuth); if
 |  x_initial can be even, you have sequences like 0, 0, 0, ...;
 |  2^31, 2^31, 2^31, ...; 2^30, 2^30, 2^30, ...; 2^29, 2^29 + 2^31,
 |  2^29, 2^29 + 2^31, ..., etc. so I force seed to be odd below.
 |
 |  Even if x_initial is odd, if x_initial is 1 mod 4 then
 |
 |    the          lowest bit of x is always 1,
 |    the  next-to-lowest bit of x is always 0,
 |    the 2nd-from-lowest bit of x alternates      ... 0 1 0 1 0 1 0 1
... ,
 |    the 3rd-from-lowest bit of x 4-cycles        ... 0 1 1 0 0 1 1 0
... ,
 |    the 4th-from-lowest bit of x has the 8-cycle ... 0 0 0 1 1 1 1 0
... ,
 |     ...
 |
 |  and if x_initial is 3 mod 4 then
 |
 |    the          lowest bit of x is always 1,
 |    the  next-to-lowest bit of x is always 1,
 |    the 2nd-from-lowest bit of x alternates      ... 0 1 0 1 0 1 0 1
... ,
 |    the 3rd-from-lowest bit of x 4-cycles        ... 0 0 1 1 0 0 1 1
... ,
 |    the 4th-from-lowest bit of x has the 8-cycle ... 0 0 1 1 1 1 0 0
... ,
 |     ...
 |
 |  The generator's potency (min. s>=0 with (69069-1)^s = 0 mod 2^32) is
 |  16, which seems to be alright by p. 25, Sec. 3.2.1.3 of Knuth.  It
 |  also does well in the dimension 2..5 spectral tests, but it could be
 |  better in dimension 6 (Line 15, Table 1, p. 106, Sec. 3.3.4, Knuth).
 |
 |  Note that the random number user does not see the values generated
 |  here directly since reloadMT() will always munge them first, so
maybe
 |  none of all of this matters.  In fact, the seed values made here
could
 |  even be extra-special desirable if the Mersenne Twister theory says
 |  so-- that's why the only change I made is to restrict to odd seeds.
 |
**/
  register uint32  x = (seed | 1U) & UINT32MAX;
  register uint32 *s = state;
  register int     j;
  for (left=0, *s++=x, j=N;   --j;   *s++ = (x *= 69069U) & UINT32MAX) {
}
}
#if defined(REAL) || defined(OPENREAL)
double reloadMT(void)
#else
uint32 reloadMT(void)
#endif
{
  register uint32 *p0 = state;
  register uint32 *p2 = state+2;
  register uint32 *pM = state+M;
  register uint32  s0;
  register uint32  s1;
  register int     j;
  if (left < -1) seedMT(4357U);
  left = N-1;
  next = state+1;
  for (s0=state[0], s1=state[1], j=N-M+1;   --j;   s0=s1, s1=*p2++) {
    *p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U);
  }
  for(pM=state, j=M;   --j;   s0=s1, s1=*p2++) {
    *p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U);
  }
  s1  = state[0];
  *p0 = *pM ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U);
  s1 ^= (s1 >> 11);
  s1 ^= (s1 <<  7) & 0x9D2C5680U;
  s1 ^= (s1 << 15) & 0xEFC60000U;
#if defined(REAL)
  return ((double) (s1 ^ (s1 >> 18)) / UINT32MAX);
#elif defined(OPENREAL)
  return ((double) (s1 ^ (s1 >> 18)) * TWOM32);
#else
  return (s1 ^ (s1 >> 18));
#endif
}
/* 'inline' the following procedure for maximum speed */
#if defined(REAL) || defined(OPENREAL)
double randomMT(void)
#else
uint32 randomMT(void)
#endif
{
  uint32 y;
  if (--left < 0) return(reloadMT());
  y  = *next++;
  y ^= (y >> 11);
  y ^= (y <<  7) & 0x9D2C5680U;
  y ^= (y << 15) & 0xEFC60000U;
#if defined(REAL)
  return ((double) (y ^ (y >> 18)) / UINT32MAX);
#elif defined(OPENREAL)
  return ((double) (y ^ (y >> 18)) * TWOM32);
#else
  return (y ^ (y >> 18));
#endif
}
int main()
{
  int j;
  /* you can seed with any uint32, but the best are odds in 0..(2^32 -
1) */
/*seedMT(4357U);*/
  seedMT(time(0));
  /* print the first 2,002 random numbers seven to a line as an example
*/
  for(j=0; j<2002; j++) {
#if defined(REAL) || defined(OPENREAL)
    printf("%10.8f%c", randomMT(), (j%7)==6 ? '\n' : ' ');
#else
    printf("%10lu%c", randomMT(), (j%7)==6 ? '\n' : ' ');
#endif
  }
  return EXIT_SUCCESS;
}
EMG

 Ivan Molella
Mah, l'ho tratto da un esempio del compilatore Turbo C della Borland, e
funziona sul mio computer dove un intero è rappresentato da 2 byte. L'ho
compilato anche sotto Linux ed era Ok. Perchè dici che non sei sicuro? dove
vedi una possibile non portabilità? Fammi sapere.
Ciao
Marco

Riporto di seguito l'intervento di Maurizio Loreti, giuntomi via email,
cosicchè tutti possano stare aggiornati sull'evoluzione della discussione.:
-------------------------------------------------------------------
aaggiungerei che prendendo il modulo si rovina pazzamente la
pseudocasualita' della sequenza.  tipicamente, prendendo il numero
pseudocasuale modulo due si ottiene una sequenza alternata di zeri e
uni...  metti un pointer al sito 
http://131.113.34.4/matumoto/emt.html
per chi vuole un generatore migliore di rand (tipicamente di bassa
qualita').
--
Maurizio Loreti
http://www.pd.infn.it/~loreti/mlo.html
Univ. of Padova, Dept. of Physics - Padova, Italy
lor...@pd.infn.it

 rappresentato da 2 byte
E quale sarebbe il problema? Più che altro quella tecnica è molto poco
soddisfacente dal punto di vista della "buona casualità" del risultato.
EMG

 Io ho questa funzione:
.........
Grazie mille, ma ho un problema: mi sembra un po' troppo lunga da inserire
nel file delle FAQ. non avresti un collegamento da cui scaricarla?
In alternativa io la metterei zippata in un file che poi linkiamo al file
delle FAQ.
fammi sapere
ciao


 vedi una possibile non portabilità? Fammi sapere.
int main(void) {
int numerocasuale ;
 srand(time(NULL)) ;
 numerocasuale=(rand() %10);
 ...}
Il problema e' che i secondi dal 1970 sono tanti;
Nel mio compilatore borland con int a 16 bit, time_t il tipo di dato restituito
da time() e' definito long, mentre srand accetta come argomento un unsigned int
L'esempio riportato dal mio borland 2.0 e' cosi':
#include <stdio.h>
#include <time.h>
int main(void)
{
   int i;
   time_t t;
   srand((unsigned) time(&t));
   printf("Ten random numbers from 0 to 99\n\n");
   for(i=0; i<10; i++)
       printf("%d\n", rand() % 100);
   return 0;
}

 delle FAQ.
Questa cosa dei file esterni andrebbe organizzata bene. Ma magari
aspettiamo se ci sarà un numero di files sufficienti a giustificarla.
EMG

 EMG
Yes, sir
Ciao
Marco

 Io metterei solo il link, se sei d'accordo.
Piů che d'accordo.
EMG