# Conservation of random bits:
#
# These four routines return uniformly distributed random numbers >= 0 and < m
# When neccessary, they call a routine to obtain more random bits.
# Careful use of floating point numbers allows operation on 52 bit
# positive integers.
#
# On the initial call, Ranval is loaded with 48 random bits. On subsequent
# calls, Ranval contains a random number consisting of the unused residual
# from the previous call normalized and right-filled (routines A and B) 
# or left-filled (routines C and D) with new random bits, between 0 and Range.
#
# The core of the routines find the largest integer r' <= Range divisible
# by m. If Ranval >= r' then we must discard some bits and try again.
# Routines A and D return floor(Ranval/(r'/m)) and retain Ranval mod (r'/m).
# Routines B and C return Ranval mod m and retain floor(Ranval/m).
#
package Ran_Int;
use strict;

use constant RAN_BITS => 52;
use constant RAN_NORM => 2**(RAN_BITS-8);

# inputs:
#	ranbyte: routine returns random int between 0 and 256
sub new {
  my $class = shift;
  my ($ranbyte) = @_;
  my $self = {};
  $self->{Ranbyte} = $ranbyte;
  $self->{Range} = 1;
  $self->{Ranval} = 0;
  return bless($self,$class);
}

# return int between 0 and m
# right fill with new random bits, return floor(V/(r'/m))
sub ranint_a {
  my ($self,$m) = @_;

  while (1) {
    while ($self->{Range}<=RAN_NORM) {		# normalize Range and Ranval
      $self->{Range} *= 256;
      $self->{Ranval} = $self->{Ranval}*256 + $self->{Ranbyte}->();
    }
    my $d = int($self->{Range}/$m);		# r'/m
    my $res = int($self->{Ranval}/$d);		# random integer
    $self->{Ranval} -= $res*$d;			# retained random bits
    if ($res==$m) { $self->{Range} -= $m*$d }
     else { $self->{Range} = $d; return $res }	# return random integer
  }
}

# return int between 0 and m
# right fill with new random bits, return V mod m
sub ranint_b {
  my ($self,$m) = @_;

  while (1) {
    while ($self->{Range}<=RAN_NORM) {		# normalize Range and Ranval
      $self->{Range} *= 256;
      $self->{Ranval} = $self->{Ranval}*256 + $self->{Ranbyte}->();
    }
    my $d = int($self->{Range}/$m);		# r'/m
    my $q = int($self->{Ranval}/$m);		# random integer
    my $res = $self->{Ranval}-$q*$m;		# retained random bits
    if ($q>=$d*$m) { $self->{Range} -= $d*$m; $self->{Ranval} = $res }
     else { $self->{Range} = $d; $self->{Ranval} = $q; return $res }	# return random integer
  }
}

# return int between 0 and m
# left fill with new random bits, return V mod m
sub ranint_c {
  my ($self,$m) = @_;

  while (1) {
    while ($self->{Range}<=RAN_NORM) {		# normalize Range and Ranval
      $self->{Ranval} += $self->{Range}*$self->{Ranbyte}->();
      $self->{Range} *= 256;
    }
    my $d = int($self->{Range}/$m);		# r'/m
    my $q = int($self->{Ranval}/$m);		# random integer
    my $res = $self->{Ranval}-$q*$m;		# retained random bits
    if ($q>=$d*$m) { $self->{Range} -= $d*$m; $self->{Ranval} = $res }
     else { $self->{Range} = $d; $self->{Ranval} = $q; return $res }	# return random integer
  }
}

# return int between 0 and m
# left fill with new random bits, return floor(V(r'/m))
sub ranint_d {
  my ($self,$m) = @_;

  while (1) {
    while ($self->{Range}<=RAN_NORM) {		# normalize Range and Ranval
      $self->{Ranval} += $self->{Range}*$self->{Ranbyte}->();
      $self->{Range} *= 256;
    }
    my $d = int($self->{Range}/$m);		# r'/m
    my $res = int($self->{Ranval}/$d);		# random integer
    $self->{Ranval} -= $res*$d;			# retained random bits
    if ($res==$m) { $self->{Range} -= $m*$d }
     else { $self->{Range} = $d; return $res }	# return random integer
  }
}

1;