Implementazioni di algoritmi/Mersenne Twister
Mersenne Twister è un algoritmo per la generazione di numeri pseudocasuali di tipo lineare congruenziale sviluppato nel 1997 da Makoto Matsumoto (松本 眞) e Takuji Nishimura (西村 拓士). È un algoritmo che genera un ottimo insieme di numeri pseudocasuali, supplendo a varie mancanze presenti negli altri algoritmi per generare numeri pseudocasuali oggi diffusi e usati (come il generatore LCG presente nel nucleo di base del C, la funzione rand()).
Ci sono almeno due varianti conosciute di queste algoritmo, che differiscono solo nel valore del numero primo di Mersenne usato. Il più nuovo ed usato è il Mersenne Twister MT 19937.
Esiste anche una nuova versione del Mersenne twister, molto più veloce, chiamata SIMD-oriented Fast Mersenne Twister (SFMT).
Questo PRNG non è da considerarsi sicuro per la crittografia, come si può leggere sul sito del creatore, è ottimo invece grazie alla sua equidistribuzione per simulazioni come il metodo montecarlo
Link
modifica#include <stdio.h>
#include <radc++.h>
#define int64 long long
//
//Varie Routine di conversione
//
//Calcola la potenza con complessità O(log ex)
inline unsigned long int pow_n(unsigned long int a, unsigned long int ex){
unsigned long int s = 1;
while(ex > 0){
if((ex & 1) != 0)
s *= a;
a *= a;
ex >>= 1;
}
return s;
}
inline unsigned long int rec(char bin) {
if (bin == 49) {
return 1 ; } ;
if (bin == 48) {
return 0 ; } ; } ;
inline void dec2bin(long int num, char *str) {//str deve essere un array di 33 elementi
int i;
for(i=1;num!=0;i++){
if(num%2==0){
str[32-i]='0';}
else{
str[32-i]='1';};
num=num/2;};
for(i;i<=32;i++){
str[32-i]='0';};
str[32]=NULL;
};
inline void double2bin(int64 num, char *str) {//str deve essere un array di 65 elementi
int i;
for(i=1;num!=0;i++){
if(num%2==0){
str[64-i]='0';}
else{
str[64-i]='1';};
num=num/2;};
for(i;i<=64;i++){
str[64-i]='0';};
str[64]=NULL;
};
inline long int bin2dec(const char bin[]) {
unsigned long int len = strlen(bin) ;
unsigned long int i, sum = 0 ;
for(i=1;i<=(len);i++) {
char a = bin[len-i] ;
sum = sum + ((rec(a)) * (pow_n(2, (i-1)))) ; }
return sum ; } ;
inline void c_first32bitof(int64 num, char *str) {//str deve essere un array di 33 elementi
char ret[65];
double2bin(num, ret);
int i;
for(i=0;i<32;i++){
str[i]=ret[i];};
str[32]=NULL;
};
inline void c_last32bitof(int64 num, char *str) {//str deve essere un array di 33 elementi
char ret[65];
double2bin(num, ret);
int i;
for(i=32;i<64;i++){
str[i-32]=ret[i];};
str[32]=NULL;
};
inline long int first32bitof(int64 n) {
char data[33] ;
(c_first32bitof(n,data)) ;
return (bin2dec(data)) ;
};
inline long int last32bitof(int64 n) {
char data[33] ;
(c_last32bitof(n,data)) ;
return (bin2dec(data)) ;
};
inline long int mergebit(long int a,long int n) {
char bina[33] ;
dec2bin(a, bina) ;
char first = bina[0] ;
char bin[33] ;
dec2bin(n, bin);
int i, len=strlen(bin) ;
if(len<31) {
char bin2[33];
bin2[0] = first;
for(i=1;i<len;i++) {
bin2[i] = '0' ; } ;
for(i;i<31;i++) {
bin2[i] = bin[i] ; } ;
return (bin2dec(bin2)) ; }
else {
char bin2[33];
bin2[0] = first;
for(i=1;i<31;i++) {
bin2[i] = bin[i] ; } ;
return (bin2dec(bin2)) ; } ;
} ;
//
//Mersenne Twister
//
long int mt[624] ;
inline void m_twister_init(long int seed) {//L'algoritmo vero e proprio
mt[0]=seed ;
int i;
for(i=1;i<624;i++) {
int64 tmp = (((0x10DCD)*(mt[i-1]))+1) ;
mt[i] = last32bitof(tmp) ; } ;
for(i=0;i<623;i++) {
long int y = mergebit(mt[i],mt[i+1]);
if(y%2==0) {
mt[i] = (mt[((i+397)%624)])^(y >> 1) ; }
else {
mt[i] = ((mt[((i+397)%624)])^(y >> 1)^(0x9908B0DF)) ; } ; } ;
long int y = mergebit(mt[i],mt[i+1]);
if(y%2==0) {
mt[i] = ((mt[((i+397)%624)])^(y >> 1)) ; }
else {
mt[i] = ((mt[((i+397)%624)])^(y >> 1)^(0x9908B0DF)) ; } ;
} ;
inline long int extract_n(long int i) {
long int y = mt[i] ;
y = y ^ (y >> 11) ;
y = y ^ ((y << 7)&&(0x9D2C5680)) ;
y = y ^ ((y << 15)&&(0xEFC60000)) ;
y = y ^ (y >> 18) ;
return y ;
} ;
inline long int quick_m_t(long int seed) {//Funzione di esempio
srand(time(NULL)+clock()) ;
srand((rand())^seed) ;
long int s = ((rand())%624) ;
m_twister_init(seed) ;
return extract_n(s%623) ;
} ;
Altri progetti
modifica
Wikipedia contiene una voce su questo algoritmo
