class Rng
   private

   # cambiando i due numeri seguenti si può avere un diverso generatore
   # Fibonacci lagged
   @@Gensize = 250
   @@Genlag = 103
      # queste costanti servono per il generatore lineare congruente
      A=16807
      M=2147483647
      Q=127773
      R=2836
      Duepi = 6.283185307

#   pointer # punta all'utimo elemento dello stack circolare

   public

   # costruttore e distruttore
   # Rng(unsigned long int seed=1)

   def initialize(seed)
      # qui sono memorizzati gli ultimi 'gensize' numeri random interi
      @buffer =  []
      # il seme e' il numero random iniziale
      k = seed

      for i in 1...@@Gensize
         # moltiplicazione modulo m che evita l'overflow
         k = A * (k % Q) - R * (k / Q)
         # con risultato negativo aggiungo m
         k += M if k < 0

         @buffer[i] = k

      end

      # inizializzazione di r250 presa dalle GSL

      msb  = 0x80000000
      mask = 0xffffffff

      for i in 0...32
         k = 7 * i + 3    # Select a word to operate on
         @buffer[k] &= mask    # Turn off bits left of the diagonal
         @buffer[k] |= msb     # Turn on the diagonal bit           */
         mask >>= 1
         msb >>= 1
      end

      # fine parte GSL

      @pointer = @@Gensize-1 # punta all ultimo posto riempito

   end

   def lfg_irand()

      # due dei 'gensize' numeri con xor (cioe' bit a bit modulo 2)
      # pointer punta all'ultimo elemento trovato
      # rispetto al nuovo elemento n punta a n-'gensize'
      # il secondo indice e' n-103, ovvero (n+147)%'gensize'
      # ma rispetto a pointer va avanzato di uno

      newrand = @buffer[@pointer] ^ @buffer[( @pointer + @@Gensize - @@Genlag + 1) % @@Gensize]
      @pointer = ( @pointer + 1 ) % @@Gensize
      @buffer[@pointer] = newrand
      # il nuovo numero e' inserito nello stack circolare

      return newrand
   end

   def lfg_uniform
      # genero un numero in doppia precisione compreso tra 0 e 1
      # a partire da un intero
      return lfg_irand / 2147483647.0
   end

   def lfg_boxmul
       # algoritmo di Box Mueller
       # genero un numero random con distribuzione gaussiana
       # a partire da due distribuiti uniformemente
      x1=lfg_uniform
      x2=lfg_uniform

      return Math.sqrt( - 2*Math.log(x1))*Math.cos(Duepi*x2)

   end

   def lfg_esponenziale
      return -Math.log(lfg_uniform.abs)
   end

   public

   def uniform
      return lfg_uniform
   end

   def gaussian
      return lfg_boxmul
   end

   def esponenziale
      return lfg_esponenziale
   end

end

# programma principale

r = Rng.new(23892384)
N = 100000

x = x2 = x4 = 0.0

# miscela iniziale
for j in 1...100000
   r.uniform
end

for j in 1..N
   t = r.uniform()
   x += t
   x2 += t*t
   x4 += t*t*t*t
end
puts "Uniforme\n"
puts"<x>"+(x / N).to_s+" <x^2> "+(x2 / N).to_s+" <x^4> "+(x4 / N).to_s+"\n"

x = x2 = x4 = 0.0
for j in 1..N
   t = r.gaussian()
   x += t
   x2 += t*t
   x4 += t*t*t*t
end
puts "Gaussiana\n"
puts"<x>"+(x / N).to_s+" <x^2> "+(x2 / N).to_s+" <x^4> "+(x4 / N).to_s+"\n"

x = x2 = x4 = 0.0
for j in 1..N
   t = r.esponenziale
   x += t
   x2 += t*t
   x4 += t*t*t*t
end
puts "esponenziale\n"
puts"<x>"+(x / N).to_s+" <x^2> "+(x2 / N).to_s+" <x^4> "+(x4 / N).to_s+"\n"