Small amounts of 'modulo bias' still exists.

In the latest beta, I'm not detecting the older behavior of thinned-out output which had caused nearly a third of outputs in a 2^24 range to not exist within more than 2^27 random numbers. Each bit increasing this range size would double the testing time, so it's the kind of thing that's hard to detect in larger ranges without a huge number of tests. I'm not detecting the odd/even bias from the latest alias randbiasdemo3.

Modulo bias always exists when the output range's size is not a divisor of the input range's size, due to some outputs matching mor inputs than other outputs do. Increasing the size of the input range doesn't eliminate it, it just makes observing the bias difficult because the number of biased inputs becomes a smaller fraction of the total inputs. The biased inputs thown away is always less than the output range's size.

The modulo bias in the latest beta is hard to observe without either looking at a huge number of inputs or looking at a very large output range. It can be seen with the earlier randbiasdemo3 alias using either $rand or $rands, though dividing the inputs into smaller odd/even groups makes it a little harder to see.

/randbiasdemo3 0.75*2^52
/randbiasdemo3 0.75*2^52 s

The modulo bias is solved by discarding enough inputs to effectively shrinking the input range enough that it becomes a multiple of the output range. The number of biased inputs to discard can be calculated either as a %throwaway_above or as a %throwaway_below value. In the simplistic example of returning a range of 10 outputs from the 256 input values 0-255, there are 6 extra inputs that need to be discarded to bring the total inputs down to 250. The two ways are:

//var -s %in_max 255 , %in_min $calc(must_remain_zero) , %out_min 1 , %out_max 10 , %out_range $calc(%out_max + 1 - %out_min) , %throwaway_above $calc(%in_max - ((%in_max % %out_range) + 1) % %out_range)

//var -s %in_max 255 , %in_min $calc(must_remain_zero) , %out_min 1 , %out_max 10 , %out_range $calc(%out_max + 1 - %out_min) , %throwaway_below $calc((%in_max - (%out_range - 1)) % %out_range)

%throwaway_above $calc( %in_max - ((%in_max % %out_range) + 1) % %out_range)
%throwaway_below $calc((%in_max - (%out_range - 1) ) % %out_range)

These always results in a keep-range which is the highest multiple of out_range <= %in_range, and doesn't throw away any inputs when the out_range and in_max are both powers of 2. In the example of outputting a range-size 10 random number from an input of 0-255:

%throwaway_above $calc( 255 - (( 255 % 10 ) + 1) % 10)
= (255 - ((5 + 1) % 10)) = 255 - 6 = 249

%throwaway_below $calc(( 255 - ( 10 - 1) ) % 10)
= (255 - 9 ) % 10 = 246 % 10 = 6

... either throwing away the 6 values 250-255 or tossing the 6 values 0-5.

The $rand_withoutmodulobias alias below removes the bias by ignoring the biased inputs when reducing a 53-bit 'double' down to the requested output range. It assumes there wasn't any bias in how the double was created, and that each of the 2^53 outputs has an equal number of inputs. It would be difficult to detect small loss of precision without a debugging parameter which would force the first fetched 64bit value to be a specific value, to see what happens from various neighboring inputs and how they are used in various output range sizes.

A 53-bit int could be created without bias from a 64-bit int by taking as many bits as can be held in a double and optionally XOR'ing them with the remaining bits so that all 64 bits can have an effect on the outcome.

53bit_val_in_uint64_rand = (uint64_rand & 0x1ffffffffffff) ^ (uint64_rand >> 53)

The $rand_WithoutModuloBias alias calculates %throwaway_above after making sure that %in_max and the size of the output range both fit within a 53-bit double. It then handles reducing the 53-bit double down to the requested output range. If the $rand_WithoutModuloBias were substituted in place of $rand in the other aliases, it would eliminate the modulo bias as a source of any output bias, but wouldn't remove the hi-odd/lo-even bias in the prior beta since that exists in the 2^53-1 output range. It should always output the identical output that $rand would've otherwise returned except in cases where the double was originally filled with one of the biased outputs that it now throws away.

Because the modulo bias is independent of the random source, it also exists if the randbiasdemo4 alias is edited to use $rands. Modulo bias favoring smaller values:

/randbiasdemo4 0.75*2^53-1
/randbiasdemo4 0.75*2^51-1

Editing randbiasdemo to use $rand_withoutmodulobias instead of $rand (2 places) eliminates the modulo bias by throwing away the extra/biased inputs. This rand_WithoutModuloBias can also be edited in 2 places to use $rands instead of $rand.

Code:
alias rand_WithoutModuloBias {
  if (($1 == $null) || ($2 == $null)) goto syntax
  if (($1 !isnum) || ($2 !isnum)) var %lo $asc($1) , %hi $asc($2) , %num 0
  else var %lo $gettok($1,1,46) , %hi $gettok($2,1,46) , %num 1 | if (%lo > %hi) var %lo $v2 , %hi $v1
  var %diff %hi - %lo , %out_range %diff + 1 , %max.val $calc(2^53-1) , %min.val $calc(-(2^53))
  if ((%hi > %max.val) || (%diff > %max.val) || (%lo < %min.val)) goto syntax
  var %throwaway_above $calc(%max.val - (((%max.val % %out_range) + 1) % %out_range) )
  var %53bitint $rand(0,%max.val) | while (%53bitint > %throwaway_above) { var %53bitint $rand(0,%max.val) }
  if (%num) return      $calc(%lo + %53bitint % %out_range)
  return           $chr($calc(%lo + %53bitint % %out_range))
  :syntax
  echo -sc info2 *$rand_WithoutModuloBias( N1,N2 |char1,char2) MinN1:-2^53 MaxN2/range.max:+2^53-1
  halt
}


Code:
alias randbiasdemo4 {
  if ($calc($1) isnum 1- $+ $calc(2^53-1)) var -s %max $gettok($v1,1,46) | else var -s %max $int($calc(0.75*(2^53-1)))
  var %pockets 256 , %array $str(0 $chr(32),%pockets)
  var -s %pocket_size $calc( ((%max + 1 - 0) // %pockets)) , %out_range %pocket_size * %pockets
  var -s %throwaway_above $calc( %max - ((%max % %out_range) + 1) % %out_range)
  var %i 25600 , %div %max / %pockets
  while (%i) {
    var %rr $rand(0,%max) | while (%rr > %throwaway_above) var %rr $rand(0,%max)
    var %t $calc(1+ (%rr // %pocket_size)) , %a $gettok(%array,%t,32) + 1 , %array $puttok(%array,%a,%t,32)
    dec %i
  }
  echo 4 -a %array = $calc($replace(%array,$chr(32),+))
  echo -a $stddev(%array)
}
alias -l stddev {
  tokenize 32 $1- | if ($0 < 2) return | var %sum 0 , %sum.of.squares 0
  var %j $0 | while (%j) { inc %sum $gettok($1-,%j,32) | dec %j }
  var %mean %sum / $0 | var %j $0 | while (%j) {
  inc %sum.of.squares $calc( ($gettok($1-,%j,32) - %mean)^2 ) | dec %j }
  var %min $gettok($sorttok($1-,32,n),1,32) , %max $gettok($sorttok($1-,32,nr),1,32)
  var %pop.variance    $calc(%sum.of.squares / ($0   ))
  var %sample.variance $calc(%sum.of.squares / ($0 -1))
  var %pop.std.dev    $calc((%sum.of.squares / ($0   ))^.5)
  var %sample.std.dev $calc((%sum.of.squares / ($0 -1))^.5)
  var %std.err $calc(%sample.std.dev / ($0 ^ .5) )
  return count: $0 sum: %sum mean: %mean min: %min max: %max popvar: %pop.variance pop.stddev: %pop.std.dev sample.var: %sample.variance sample.std.dev: %sample.std.dev std.err: %std.err
}