Newsgroup: it.comp.lang.c++
Titolo: randomize() , srand() ... come usarle?

Salve a tutti voi del NG.
Ho letto che randomize(), sebbene non standard, ha un comportamento similare
ad srand(), cioè inizializzano "il cosiddetto" seme del generatore random.
Si è anche ampiamente discusso sul NG it.comp.lang.c++ dell'inaffidabilità
di rand(), e che esistono metodi più efficienti per generare distribuzioni
uniformi.
Tuttavia io ho notato un "effetto collaterale" che a me porta un certo
vantaggio, e sarebbe che la sequenza generata da un ciclo del genere:
[Ad esempio i<30000]
while(i<30000)
{
int n1,n2; float p;
n1=rand()%1000;
randomize();
n2=rand%1000;
p=1000*n1+n2;
printf("p= %f\n", p);
}
è caratterizzata da sottosequenze, ciascuna delle quali è costituita da
numeri "fra loro vicini".
Per capire cosa intendo con questa frase un pò approssimativa pensate di
stabilire una soglia prefissata
int soglia= ... quello che stabiliamo noi di mettere, ( ad es soglia=100;
200, ecc. ecc.)
e di "complicare" un pochetto il ciclo effettuando un controllo del tipo
if(p<soglia) printf("*");
Quello che otteniamo è "una sequenza di asterischi del tipo"
******      **                      *******          *****
****          *****               *******
[la barra spaziatrice inica che la condizione non si è verificata.
Immaginate che questo sua l'output sullo screen.]
cioè, se ci si ragiona un pò su, questo avviene perchè p, in iterazioni
successive e vicine fra loro, assume valori pressochè "vicini fra loro",
cioè inferiori alla soglia.
Un comportamento di questo tipo, nelle telecomunicazioni, si chiama bursty
 che significa esplosivo).
E' come se la sorgente estraesse un tot di asterischi e poi stà zitta per un
pò, per poi voutare i sacco di un altro po di asterischi.
Non si ottiene mai, o quasi mai, una sequenza del tipo:
*  **   *  * *    **  **     ****   * * *   **    * *   *  *     * * *
**   **   * * **    ***   * * * * *
cioè dove la prob. di "attività" degli asterischi è uniformemente
distribuita secondo una distibuz. di prob. fissata dalla soglia.
Io credo che tale comportamento sia dovuto a randomize(), o meglio all'uso,
probabilmente improprio che ne faccio.
E allora quello che voglio chiedere è questo.
Se inizializzassi il seme solo prima del ciclo, cioè una chiamata di
randomize() effettuata solo una volta,
l'errore dovuto al non poter ottenere numeri a 6 cifre del tipo 304304,
oppure 589589 (è difficile, davvero molto difficile, che per due "tiri"
successivi la rand()%1000 fornisca esattamente 304) è o meno più grande
dell'errore con cui ottengo i numeri col ciclo in cui chiamo randomize()  (-
oppure srand() -) ogni due colpi?
Cioè, la sequenza che ottengo con una chiamata iniziale di randomize() e poi
basta ha una maggiore parvenza di essere una sequenza pseudo-casuale
rispetto a quella ottenuta col metodo inizialmente esposto?
Vi ringrazio moltissimo per le vostre risposte.
Mi servirebbe per la tesi.
Bye,
Sebastiano

Il 19 scorso ti sei lamentato che "è da una settimana che *non* parlate di
C++ ma della problematica alla base della generazionedei numeri *veramente*
casuali." adesso stai chiedendo risposte proprio a tal riguardo?!
Sii coerente, che ca**o!
Comunque qui bisogna aprire un ng chiamato 
rnd.it
; non ho mai visto una cosa
che susciti un tale interesse!
Cristiano

Leggi quello che ti ho risposto su it.comp.lang.c
Ah. Grazie per *l'estrema* disponibilità.
Chi non sa che rispondere è meglio che si stia zitto, anzichè offendere.

Allora,
voglio precisare delle cose una volta per tutte.
E' vero che sono stato prolisso sulla discussione ancora in corso su
it.comp.lang.c++,
(ieri sera ho riletto, fra gli altri, i miei messaggi sul *rumore bianco*),
ma ho esordito all'inizio
(questo spero che lo rammendiate) che le mie non volevano assumere il tono
di polemica o, peggio, di insulto.
Erano opinioni. Basta cosi.
E poi non mi sono *lamantato*, semmai ho tentato di precisare un concetto
che ritengo sia di base.
Quella ::cosa:: sul rumore bianco.
Ok. Forse il tono è risultato un pò aspro, addirittura si è capito, in base
alle vostre reply, che io cercassi la precisione assoluta nella generazione
dei casuali.
Non è affatto cosi. Anzi.
Volevo farmi portabandiera di chi diceva che per talune applicazioni "rand()
migliorata" poteva andare anche bene.
Ripeto, non volevo insultare.
Quindi a questo punto non capisco cosa mi dovrebbero ritornare... insulti???
(Però ditemi voi se "e che ca**o!" non lo è. O sono un puritano, ma non mi
pare, o forse non è stata *legittima difesa* rispondere al tono, diciamo
*sarcastico*, per non inasprire oltre, di Cristiano)
Le critiche sono sempre accettatissime. Rimango sempre aperto al dialogo...
con qualsiasi tono esso sia... ma dialogo!
Io parto sempre da un principio. (al solito, forse sbagliato)
Il NG non è stato instituito per fare sfoggio di sapienza, ma per scambiare
opinioni (sul tema indicato).
Ho capito male?
Per favore, rispondetemi almeno a questo. Potrebbe significare per me
cambiare il modo di postare sui NG. Voglio rimanere sempre entro *la
netiquette*.
Se ho capito bene, perchè vi arrabbiate quando chi conosce poco del
linguaggo fa delle osservazioni elementari?
Come faccio ad imparare se mi tirate le frecce? Dai libri? Ma a scuola ci
siamo andato tutti, e talvolta leggendo dai libri non è detto che si capisca
bene tutto.... parlarne fa sempre bene ;-)
E quale posto migliore se non il NG? (Almeno io ho pensato che fosse cosi)
Il tono di Cristiano mi è sembrato fuori luogo.
Lui si professa di animo mite, ma non risparmia le frecciatine ;-)
Riguardo, in ultima analisi, la questione del cross-mailing, non mi pare che
sia vietato. O si?
La risposta già la sapete. Se è vietato non lo faccio più.
Sono anticonformista, ma non rompiscatole.
Spero di aver chiarito un pò la mia posizione.
Che poi non si è discusso di c++, intendevo dire che non ho potuto
trattenermi quando ho visto delle *cose* che non erano corrette sul *rumore
bianco*.
Che le cose funzionino è un conto, ma bisogna usare i termini corretti.
Quel "famoso" rumore di cui si è tanto chiaccherato non è ne bianco, ne
termico, ne rosa.. ne a pois !!
E' solo rumore. E c'è una bella differenza.
Ok. Avete ragione sul dire che a ribattere sono stato polemico.
E' vero. Ripeto, avete ragione.
Ma se qualcuno di noi chiamasse template una classe voi non vi
"incacchiereste" un pochino?
I termini giusti al posto giusto.
Sono stato bocciato ad un esame solo per questo. Per un tecnico
specializzato (un ing., un buon programmatore, un bravo idraulico,
ecc.,ecc., metteteci pure qualunque professione ) usare i termini giusto al
momento giusto è , a mio parere, vitale per evitare di essere fraintesi.
Ci sarebbero "problemi di comunicazione" per dirla con un'espressione di
"altri" :-)))
Troppo *rigore* il mio?
Beh! allora quel "calcio di punizione" che avete già battuto lo meritavo :-)
Ma vi ripeto lo spirito dei mie post su NG (c++) volevano solo dire:
rumore bianco != rumore_comunemente_usato;
rumore_bianco_esistente_solo_in_teoria=true;
rumore_bianco_esistente_in_pratica!=bianco;
Tutto qui.
Continueremo a discutere ?
Bye,
Sebastiano

 Continueremo a discutere ?
E chi ti ferma piu' ? :-)
Ciao,
Nicola

 di L. Carroll.
Non ho letto il libro ... personalmente come gia' detto
penso che Cristiano non ragioni con la logica. Oltre a non dare
un peso ai dettagli (n non e' la stessa cosa di n+1, 10 non e'
la stessa cosa di MAX, rng_32 non e' la stessa cosa di rng)
penso anche che prenda i commenti ai suoi post come offese
personali. Per questo ho smesso qualche giorno fa di
rispondere... e ora anche di leggere (il rapporto
segnale/rumore e' troppo basso ;) ).
Andrea
--
----------------------------------------------------------------
 <\___/>  Andrea        email: agriff at tin dot it
 / . . \  "6502"          irc: _6502_@EFNET#italia
(  (")  ) Griffini   homepage: 
http://www.netservers.com/~agriff
Sent via Deja.com 
http://www.deja.com/
Before you buy.

Scusa la domanda, magari banale, ma hai provato ad usare e misurare
l'algoritmo r250 (periodicità dell'algoritmo 2^249), che è un algoritmo
congruente ed ha delle performance molto superiori agli algoritmi
convolutivi come quelli del Numerical Recipies o implementati di default nei
linguaggi? Io lo usai in tesi e poi con le trasformazioni di Box-Muller
ottenevo distribuzioni gaussiane. Si trova scritto sia in C, C++, Fortran,
etc. e naturalmente è anche velocissimo.
Non ho letto tutti i vostri post precedenti che erano veramente troppi,
perciò potrei aver detto cose scontate o banali.
Bye, Francy

 Continueremo a discutere ?
Io continuerň a postare fin quando ne avrň voglia nonostante le tue
intimazioni.
Cristiano

 Chi non sa che rispondere è meglio che si stia zitto, anzichè offendere.
Che non sappia rispondere, nessuno può dirlo (eccetto tu e la tua
maleducazione). Te lo dimostra il fatto che nella recente discussione sui
numeri casuali forse ho scritto più io che tutti gli altri messi insieme
(questo ovviamente non significa che io sono il più "bravo" come sembra che
vuoi far credere tu).
Cristiano

 l'algoritmo r250 (periodicità dell'algoritmo 2^249),
Non sono pronto a giurarci, ma credo che il periodo non sia 2^249.
Lo dico perché il generatore (detto di Tausworthe) ottenuto dalla
congruenza:
x[n]=x[n-460] + x[n-607]
ha un periodo pari a 2^(s-1)*(2^607-1) dove s è il numero dei bits di x[].
Per cui l'r250 che è ottenuto dalla congruenza:
x[n]=x[n-250] ^ x[n-147]
dovrebbe avere un periodo superiore a quello che dici (ma ripeto che non
posso giurarci), dopotutto lo xor è sempre un'addizione modulo 2.
Comunque dalle molte prove che ho svolto su questo generatore (r250), ho
notato che non è statisticamente molto valido (presenta una discreta
autocorrelazione tra i bits) ed inoltre ho notato che la sua
inizializzazione è piuttosto critica, cioè influisce molto
sull'autocorrelazione dei bits. Tengo a precisare che quello che ho detto
non ha alcun supporto matematico, scaturisce solo dalle mie decine di test
statistici a cui sto sottoponendo un gran numero di generatori.
Ovviamente alla base di tutto c'è sempre l'utilizzo che se ne deve fare con
un generatore...
Cristiano

 Cristiano
Ciao, Francesco

Non ci sono parole... finisco qui per non andare sul tono polemico.
Ti faccio riflettere però su una cosa: e che cazzo è una frase volgare (ora
sono io che la scrivo). Se non lo fosse perchè quegli infantili **???
Accetto in consiglio di Andrea, ma non posso farmi tirare fango addosso
gratuitamente.
Imparala tu la netiquette...
Almeno io so che differenza c'è tra la mancanza di educazione e
l'arroganza... anche se posso essere stato "inavvertitamente" maleducato nei
tuoi riguardi...
Con questo passo e chiudo sul serio.

 tuoi riguardi...
Hai fatto bene a mettere tra virgolette "inavvertitamente" perché è
veramente un eufemismo. Per la tua risposta dell'altra volta "Non ci sono
parole...".
Un'ultima cosa. Prima di rispondere a qualcuno bisogna che leggi interamente
quello che ti viene detto. Altra lezione di buone maniere...
Cristiano

 frequenza di taglio del filtro.
M'interessa molto ciò che dici perché sono alla ricerca (da ormai sei mesi)
di un valido generatore.
Purtroppo non riesco a capire cosa intendi per "tempo di campionamento" e il
resto.
Un genereatore di numeri pseudocasuli non contempla tali parametri.
Io inizializzo le 250 long dello stato e poi ad ogni chiamata esso mi
restituisce una long.
Il codice del generatore che ho testato ha questa intestazione:
" Kirkpatrick, S., and E. Stoll, 1981; "A Very Fast
  Shift-Register Sequence Random Number Generator",
  Journal of Computational Physics, V.40"
Cristiano

Chiedo scusa al NG. Ma A Cristiano rispondo privatamente, per evitare che la
piazza si sporchi di cose che non ne sono degne. Lo è già stata troppo.
Ma non posso tenermi il maleducato ( o meglio, gran maleducato) quando chi
mi accusa non mi pare che abbia fatto le giuste mosse sin dall'inizio.
Se il mio tono vi sembra *maleducato* beh, che dire, eviterò di postare da
questo momento in poi per sempre in avanti. Mi pare che il NG vada oltre
queste assurdità.
Se sono un maleducato.
Se invece credete che abbia un minimo di giustificazione, anche se i rendo
conto di aver surriscaldato un pò l'aria, fatelo sapere...
Sebastiano

 Cristiano
Si, si. E' quello che ho anche io e viene usato in ambito scinetifico, poi come
ho detto non l'ho testato a fondo, dato che era già usato e mi sono fidato.
Oltretutto le misurazioni che dovevo effettuare non erano particolarmente
critiche.
Parlare di tempo di campionamento è improprio, però rende l'idea. Immagina di
avere una grandezza (ad. ex. una tensione elettrica) che varia nel tempo, te
puoi pensare di fotografarla ad una frequenza ben precisa (cioè la campioni) e
quella è la frequenza di campionamento. Il risoltato ottenuto è rappresentato
da una sequenza di numeri che li pensi distanziati di un certo intervallo
temporale. Per analogia con il caso reale appena descritto, puoi avere una
sequenza di numeri generati numericamente (possono essere generati da rand(), o
sin(), o quello che vuoi). Una volta stabilità la base dei tempi, cioè una
volta che hai piazzato le ascisse di questi punti puoi definire tempi di
autocorrelazione e tutto quello che vuoi (è una simulazione della realtà,
ovvero simulazione numerica).
Quello che facevo io nelle mie sumulazioni era generare una sfilza di numeri,
applicargli le trasformazioni di box-muller per trasformare la distribuzione
uniforme in gaussiana e applicare un filtro passabasso. La filtratura è
essenziale alla stregua della realisticità delle simulazioni, perché quello che
esce fuori dal generatore casuale è uno spettro bianco ad un certo punto
tagliato di netto in corrispondenza della metà della frequenza di
campionamento. Il filtro passa basso fa in modo di abbassare gradualmentelmente
la coda dello spettro.
Comunque il concetto sperò che sia chiaro, poi io è da alcuni anni che non
faccio queste cose e potrei avere detto delle imprecisioni, chiedo venia.
Se vuoi ti posso mandare gli algoritmi che qualcuno ne ho ancora, sono in C e
sono veramente pochi byte.
Ciao, Francesco

 sono veramente pochi byte.
Sě, credo sia la cosa migliore.
Nonostante la tua dettagliata descrizione non sono ancora riuscito a capire
se per autocorrelazione intendiamo la stessa cosa.
Credo che puoi postare i listati qui nel ng (speriamo che nessuno del C++
s'offenda per i listati in C!).
Cristiano

 conto di aver surriscaldato un pò l'aria, fatelo sapere...
Te lo hanno fatto sapere le risposte che hai ricevuto: solo una
dell'argomento richiesto (alla quale non hai nemmeno risposto!) e una
valanga che ti dicono che sei un gran maleducato rompiscatole!
Se ancora non l'hai capito...
Cristiano

 a capirlo la dice lunga...
Per favore, potete continuare in privato?
Grazie.
EMG

> Per favore, potete continuare in privato?
No. Io la pianto qui e mi scuso per tale conversazione.
Cristiano

 Cristiano
Allora i listati sono questi (chiedo venia per la non eleganza del mio codice,
ma erano le prime righe che scrissi in c), che spero che possano essere di
qualche aiuto:
_______________________________________________________________
nwa.c (un file che usava gli algoritmi in oggetto):
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "noise.h"
int main(int argc,char *argv[])
{
 FILE *pe,*pf;
 int fnum,j,i,estp,rstp,dd,rseed,elo,rlo,arm,fondo,quad;
 double
ftmb,tend,tmb,time,amp,ramp1,ramp2,ramp,pi,fper,jfreq,tau,pws,pwsm,pwsn,pwsf,
frqrc,soglia,eps1,eps2,epsilon,jff,jffm,jffn,den,pwsr,pwsi,pwsrm,pwsrn,pwsim,pwsin;
 char jx,nametxt[13],namedat[13];
 if(argc==1){
  printf("manca nome file configurazione");
  exit(1);
 }
 for(dd=1;dd<argc;dd++){
  strcpy(namedat,argv[dd]);
  strcat(namedat,".dat");
  if((pf=fopen(namedat,"w"))==NULL){
   printf("non posso scrivere il file %s",namedat);
   exit(1);
  }
  fclose(pf);
  strcpy(nametxt,argv[dd]);
  strcat(nametxt,".txt");
  if((pe=fopen(nametxt,"r"))==NULL){
   printf("non posso aprire il file %s",nametxt);
   exit(1);
  }
  fscanf(pe,"%d%*s %lg%*s %lg%*s %lg%*s",&fnum,&jfreq,&soglia,&amp);
  fscanf(pe,"%lg%*s %lg%*s %d%*s %d%*s",&ftmb,&frqrc,&quad,&rseed);
  fscanf(pe,"%lg%*s %lg%*s %d%*s %d%*s",&ramp1,&ramp2,&rstp,&rlo);
  fscanf(pe,"%lg%*s %lg%*s %d%*s %d%*s",&eps1,&eps2,&estp,&elo);
  fscanf(pe,"%d%*s %d%*s",&arm,&fondo);
  printf("-------** amp=%g soglia=%g jfreq=%g frqrc=%g
**-------\n",amp,soglia,jfreq,frqrc);
  fclose(pe);
  pi=atan(1)*4;
  tmb=.5/ftmb;
  tend=fnum/jfreq;
  tau=1/frqrc;
  jff=2*pi*jfreq*arm;
  jffm=2*pi*jfreq*(arm-1.0/fnum);
  jffn=2*pi*jfreq*(arm+1.0/fnum);
  den=1.0/(2.0*fnum*ftmb/jfreq*fnum*ftmb/jfreq);
  noise_seed(rseed);
  for(j=0;j<estp;j++){
   if(estp==1)epsilon=eps1;
   else{
    if(elo)epsilon=eps1+j*(eps2-eps1)/(estp-1.0);
    else epsilon=eps1*pow(pow(eps2/eps1,1.0/(estp-1.0)),(double)j);
   }
   for(i=0;i<rstp;i++){
    if(rstp==1)ramp=ramp1;
    else{
     if(rlo)ramp=ramp1+i*(ramp2-ramp1)/(rstp-1.0);
     else ramp=ramp1*pow(pow(ramp2/ramp1,1.0/(rstp-1.0)),(double)i);
    }
    time=pwsr=pwsi=pwsrm=pwsrn=pwsim=pwsin=0;
    jx=-1;
    noise_init(tmb,tau,ramp*soglia);
    do{
     if(quad==0)fper=amp*sin(2*pi*jfreq*time)+noise();
     else{
      double fmix,fint;
      fmix=jfreq*time;
      if((fint=fmix-floor(fmix))<.5)fper=noise()-amp;
      else if(fint>.5)fper=noise()+amp;
     }
     if(fper<(soglia*(epsilon-1)))jx=-1;
     if(fper>(soglia*(epsilon+1)))jx=1;
     pwsr=pwsr+jx*cos(jff*time);
     pwsi=pwsi+jx*sin(jff*time);
     if(fondo==1){
      pwsrm=pwsrm+jx*cos(jffm*time);
      pwsim=pwsim+jx*sin(jffm*time);
      pwsrn=pwsrn+jx*cos(jffn*time);
      pwsin=pwsin+jx*sin(jffn*time);
     }
     time+=tmb;
    }while(time<tend);
    pws=(pwsr*pwsr+pwsi*pwsi)*den;
    if(fondo==1){
     pwsm=(pwsrm*pwsrm+pwsim*pwsim)*den;
     pwsn=(pwsrn*pwsrn+pwsin*pwsin)*den;
     pwsf=(pwsm+pwsn)/2;
    }
    pf=fopen(namedat,"a");
    printf("e= %g sdev= %g --> [pws= %g]",epsilon,ramp,pws);
    fprintf(pf,"%g %g %g",epsilon,ramp,pws);
    if(fondo==1){
     printf(" fondo= %g\n",pwsf);
     fprintf(pf," %g\n",pwsf);
    }
    else{
     printf("\n");
     fprintf(pf,"\n");
             }
    fclose(pf);
   }
  }
 }
 printf("\n*** ciao! ***\n");
 return 0;
}
_________________________________________________
A017wa.txt (file da passare in argomento all'eseguibile):
3000         statistica
100         frequenza
200         soglia
201         ampiezza
200000      campionamento
500       filtro
0           0:sinu_1:quad
2511        seme_rumore
.001           rumore_iniziale
100          rumore_finale
41          passi
0           0:log_1:lin
0           epsilon_iniziale
.4          epsilon_finale
1           passi
1           0:log_1:lin
1           armonica
0           0:nofondo_1:fondo
200         scansione
_____________________________________________
noise.c (mia implementazione del filtro e box-muller):
#include <math.h>
#include "r250.h"
extern double noise_sa,noise_r1,noise_r2;
double noise(){
 double w,r1,r2;
 static double rsave=0;
 do{
  r1=2.0*dr250()-1;
  r2=2.0*dr250()-1;
  w=r1*r1+r2*r2;
 }while(w>=1);
 w=sqrt(-2*log(w)/w);
 r1=r1*w;
 r2=r2*w;
 rsave=rsave*noise_sa+r1*noise_r1+r2*noise_r2;
 return rsave;
}
double noise_sa,noise_r1,noise_r2;
void noise_init(double h,double tau,double sdev){
 noise_sa=1-h/tau+.5*h/tau*h/tau;
   noise_r1=sdev*sqrt(2*tau*h)/tau-.5*h*sdev*sqrt(2*tau*h)/tau/tau;
   noise_r2=-.5*h*sdev*sqrt(2*tau*h/3)/tau/tau;
}
void noise_seed(int seed){
 r250_init(seed);
}
______________________________________________________
noise.h (relativo header):
double noise();
void noise_init(double h,double tau,double sdev);
void noise_seed(int seed);
_____________________________________________________
r250.c
/* r250.c the r250 uniform random number algorithm
  Kirkpatrick, S., and E. Stoll, 1981; "A Very Fast
  Shift-Register Sequence Random Number Generator",
  Journal of Computational Physics, V.40
  also:
  see W.L. Maier, DDJ May 1991
*/
static char rcsid[] = "@(#)r250.c 1.2 15:50:31 11/21/94   EFC";
#include <limits.h>
#include "r250.h"
/* set the following if you trust rand(), otherwise the minimal standard
   generator is used
*/
/* #define TRUST_RAND */
#ifndef TRUST_RAND
#include "randlcg.h"
#endif
/* defines to allow for 16 or 32 bit integers */
#define BITS 31
#if WORD_BIT == 32
#ifndef BITS
#define BITS 32
#endif
#else
#ifndef BITS
#define BITS    16
#endif
#endif
#if BITS == 31
#define MSB          0x40000000L
#define ALL_BITS     0x7fffffffL
#define HALF_RANGE   0x20000000L
#define STEP         7
#endif
#if BITS == 32
#define MSB          0x80000000L
#define ALL_BITS     0xffffffffL
#define HALF_RANGE   0x40000000L
#define STEP         7
#endif
#if BITS == 16
#define MSB         0x8000
#define ALL_BITS    0xffff
#define HALF_RANGE  0x4000
#define STEP        11
#endif
static unsigned int r250_buffer[ 250 ];
static int r250_index;
#ifdef NO_PROTO
void r250_init(sd)
int seed;
#else
void r250_init(int sd)
#endif
{
 int j, k;
 unsigned int mask, msb;
#ifdef TRUST_RAND
#if BITS == 32 || BITS == 31
 srand48( sd );
#else
 srand( sd );
#endif
#else
 set_seed( sd );
#endif
 r250_index = 0;
 for (j = 0; j < 250; j++)      /* fill r250 buffer with BITS-1 bit values */
#ifdef TRUST_RAND
#if BITS == 32 || BITS == 31
  r250_buffer[j] = (unsigned int)lrand48();
#else
  r250_buffer[j] = rand();
#endif
#else
  r250_buffer[j] = randlcg();
#endif
 for (j = 0; j < 250; j++) /* set some MSBs to 1 */
#ifdef TRUST_RAND
  if ( rand() > HALF_RANGE )
   r250_buffer[j] |= MSB;
#else
  if ( randlcg() > HALF_RANGE )
   r250_buffer[j] |= MSB;
#endif
 msb = MSB;         /* turn on diagonal bit */
 mask = ALL_BITS; /* turn off the leftmost bits */
 for (j = 0; j < BITS; j++)
 {
  k = STEP * j + 3; /* select a word to operate on */
  r250_buffer[k] &= mask; /* turn off bits left of the diagonal */
  r250_buffer[k] |= msb; /* turn on the diagonal bit */
  mask >>= 1;
  msb  >>= 1;
 }
}
unsigned int r250()  /* returns a random unsigned integer */
{
 register int j;
 register unsigned int new_rand;
 if ( r250_index >= 147 )
  j = r250_index - 147; /* wrap pointer around */
 else
  j = r250_index + 103;
 new_rand = r250_buffer[ r250_index ] ^ r250_buffer[ j ];
 r250_buffer[ r250_index ] = new_rand;
 if ( r250_index >= 249 ) /* increment pointer for next time */
  r250_index = 0;
 else
  r250_index++;
 return new_rand;
}
double dr250()  /* returns a random double in range 0..1 */
{
 register int j;
 register unsigned int new_rand;
 if ( r250_index >= 147 )
  j = r250_index - 147; /* wrap pointer around */
 else
  j = r250_index + 103;
 new_rand = r250_buffer[ r250_index ] ^ r250_buffer[ j ];
 r250_buffer[ r250_index ] = new_rand;
 if ( r250_index >= 249 ) /* increment pointer for next time */
  r250_index = 0;
 else
  r250_index++;
 return (double)new_rand / ALL_BITS;
}
#ifdef MAIN
/* test driver prints out either NMR_RAND values or a histogram */
#include <stdio.h>
#define NMR_RAND 5000
#define MAX_BINS 500
#ifdef NO_PROTO
void main(argc, argv)
int argc;
char **argv;
#else
void main(int argc, char **argv)
#endif
{
 int j,k,nmr_bins,seed;
 int bins[MAX_BINS];
 double randm, bin_inc;
 double bin_limit[MAX_BINS];
 if ( argc != 3 )
 {
  printf("Usage -- %s nmr_bins seed\n", argv[0]);
  exit(1);
 }
 nmr_bins = atoi( argv[1] );
 if ( nmr_bins > MAX_BINS )
 {
  printf("ERROR -- maximum number of bins is %d\n", MAX_BINS);
  exit(1);
 }
 seed = atoi( argv[2] );
 r250_init( seed );
 if ( nmr_bins < 1 ) /* just print out the numbers */
 {
  for (j = 0; j < NMR_RAND; j++)
   printf("%f\n", dr250() );
  exit(0);
 }
 bin_inc = 1.0 / nmr_bins;
 for (j = 0; j < nmr_bins; j++) /* initialize bins to zero */
 {
  bins[j] = 0;
  bin_limit[j] = (j + 1) * bin_inc;
 }
 bin_limit[nmr_bins-1] = 1.0e7; /* make sure all others are in last bin */
 for (j = 0; j < NMR_RAND; j++)
 {
  randm = r250() / (double)ALL_BITS;
  for (k = 0; k < nmr_bins; k++)
   if ( randm < bin_limit[k] )
   {
    bins[k]++;
    break;
   }
 }
 for (j = 0; j < nmr_bins; j++)
  printf("%d\n", bins[j]);
}
#endif
______________________________________________________
r250.h (relativo header):
/* r250.h prototypes for r250 random number generator,
  Kirkpatrick, S., and E. Stoll, 1981; "A Very Fast
  Shift-Register Sequence Random Number Generator",
  Journal of Computational Physics, V.40
  also:
  see W.L. Maier, DDJ May 1991
*/
#ifndef _R250_H_
#define _R250_H_ 1.2
#ifdef __cplusplus
extern "C" {
#endif
#ifdef NO_PROTO
void         r250_init();
unsigned int r250();
double      dr250();
#else
void         r250_init(int seed);
unsigned int r250( void );
double       dr250( void );
#endif
#ifdef __cplusplus
}
#endif
#endif

 #define BITS 31
Avevi qualche motivo particolare per utilizzare soltanto 31 bits (anziché
32)?
Ho esaminato gl'interessanti listati, in particolare noise.c che contiene
dr250() che genera 2^31 numeri da 0 a 1 inclusi (in teoria saresti potuto
incappare in un errore di dominio per log(w) o in una divisione per 0! Ma
questo forse non c'interessa.).
Perň non capisco perché ritieni che le caratteristiche di autocorrelazione
siano fondamentali per l'utilizzo che fai di r250(). Io credo che se tu
utilizzassi (per assurdo) l'orripilante rand() al posto di r250() forse
nemmeno te ne accorgeresti per quanto viene filtrata e strafiltrata.
In ogni caso mi piacerebbe avere un chiarimento: t'interessa (o interessava)
un generatore con i bits molto correlati o scarsamente correlati?
Vorrei aggiungere che utilizzare soltanto pochi bits per inizializzare un
generatore che dispone di una vasta area d'inizializzazione (r250() dispone
di uno stato da ben 250 long) mi sembra un peccato! :-)
Conviene inventarsi un modo migliore per riempire tutto lo stato anche con
numeri non proprio casuali; č sempre meglio che utilizzare un semplice
generatore (inizializzato con time(0) o simili) che non ha nemmeno delle
grandi caratteristiche statistiche.
Cristiano

 Cristiano
Per quanto riguarda le questioni riguardo r250.c ho preso da un sito di ricerca
già inizializzato (ho preferito questo perché già testato). Effettivamente a
posteriori potrei darti ragione rigurdo l'inizializzazione, ma oramai poco
importa... Per l'errore di dominio non mi preoccupavo, dovrebbe essere
abbastanza raro :-) un problema più da informatico che da fisico.
L'autocorrelazione per me era fondamentale che ce ne fosse meno possibile,
perché per le misure fatte era giusto che ci fosse solo quella determinata da
noise.c, che come hai visto fa un bell'integrale e perciò funge da filtro
passabasso. Comunque caratteristica di r250.c è che oltre ad avere buone misure
spettrali, è veloce, dato che un singolo punto del grafico poteva essere
running anche per qualche giorno....
In realtà forse è più critica la realizzazione del filtro che non del
generatore di rumore, dato che da quello dipendevano tutte le prorietà del mio
sistema.
Tanto per renderti partecipe, il sistema serviva per far delle misure correlate
al fenomeno della risonanza stocastica e pur essendo molto semplice il tutto ha
contribuito alla pubblicazione di quattro articoli.
Ciao, Francesco

 wrote in message
39FE1ED3...@libero.it...
