use strict;
use warnings;
use Math::BigInt lib => "Pari,GMP";
use Math::BigRat;

# Random number generator parameters
#     $x = ($a*$x + $b) % $m
# These need to match the RNG your perl binary was compiled with.
# Run this one-liner and check to see if the output matches one
# of the signatures below.
#     srand 1;@x=map int(256*rand),1..10;print "@x\n"

# drand48() -- the default on Linux
# Signature: 10 116 213 86 144 0 48 253 192 93
my $a = Math::BigInt->new("0x5DEECE66D");
my $b = Math::BigInt->new("0xB");
my $m = Math::BigInt->new(1) << 48;
my $n = 7;  # number of input bytes needed

# ActiveState perl for windows
# Signature: 0 144 49 207 149 122 89 229 210 191
#my $a = Math::BigInt->new("0x343FD");
#my $b = Math::BigInt->new("0x269EC3");
#my $m = Math::BigInt->new(1) << 31;
#my $n = 4;  # number of bytes needed

my $one_half = Math::BigRat->new("1/2");
my $three_fourths = Math::BigRat->new("3/4");

my @y = map int(256*rand), 1..$n;
print "rand = @y\n";

my @mat = map [(0) x ($n+1)], 0..$n;
my $k = $m >> 8;
$mat[0][0] = $k;
$mat[1][1] = 1;
foreach (2..$n) {
  $mat[0][$_] = ($mat[0][$_-1] * $a + $b) % $m;
  $mat[1][$_] = ($mat[1][$_-1] * $a) % $m;
  $mat[$_][$_] = $m;
}
$mat[0][$_] -= $y[$_-1] * $k for 1..$n;

l3_reduce(\@mat);

CHECK: foreach my $row (@mat) {
  if ($row->[0] < 0) { $_ = -$_ for @$row }
  next CHECK unless $row->[0] == $k;
  foreach (1..$n) {
    next CHECK unless $row->[$_] >= 0 && $row->[$_] < $k;
  }
  my $x = $y[$n-1]*$k + $row->[$n];
  my @p;
  foreach (1..10) {
    $x = ($a*$x + $b) % $m;
    push @p, $x / $k;
  }
  print "prediction = @p\n";
}

my @z = map int(256*rand), 1..10;
print "rand = @z\n";

# l3 lattice basis reduction algorithm
# @$m is an array of lattice basis vectors
# $m is reduced in-place (the original is destroyed)
sub l3_reduce {
  my ($m) = @_;
  my $n = @$m;
  my (@m2, @mag, @coeff);

  foreach my $i (0..$n-1) {
    my @v = map Math::BigRat->new($_), @{$m->[$i]};
    foreach my $j (0..$i-1) {
      my $w = $m2[$j];
      my $c = dot_product(\@v, $w) / $mag[$j];
      $coeff[$i][$j] = $c;
      foreach my $k (0..$#v) {
        $v[$k] -= $c * $w->[$k];
      }
    }
    $m2[$i] = \@v;
    $mag[$i] = dot_product(\@v, \@v);
  }
  undef @m2;

  my $i = 1;
  while ($i < $n) {
    reduce_one($m, \@coeff, $i, $i-1);
    my $u = $coeff[$i][$i-1];
    my $u2 = $u*$u;
    if (rcmp($mag[$i], ($three_fourths - $u2) * $mag[$i-1]) < 0) {
      my $t = $mag[$i] + $u2*$mag[$i-1];
      my $s = $mag[$i-1] / $t;
      $mag[$i-1] = $t;
      $mag[$i] *= $s;
      $coeff[$i][$i-1] *= $s;

      @$m[$i, $i-1] = @$m[$i-1, $i];
      foreach my $j (0..$i-2) {
        my $t = $coeff[$i][$j];
        $coeff[$i][$j] = $coeff[$i-1][$j];
        $coeff[$i-1][$j] = $t;
      }

      foreach my $j ($i+1..$n-1) {
        my $t = $coeff[$j][$i];
        $coeff[$j][$i] = $coeff[$j][$i-1] - $u*$t;
        $coeff[$j][$i-1] = $t + $coeff[$i][$i-1]*$coeff[$j][$i];
      }

      $i-- if $i > 1;
    }
    else {
      foreach my $j (reverse 0..$i-2) {
        reduce_one($m, \@coeff, $i, $j);
      }
      $i++;
    }
  }

  return $m;
} # l3_reduce

sub reduce_one {
  my ($m, $coeff, $i, $j) = @_;
  if (rcmp(abs($coeff->[$i][$j]), $one_half) > 0) {
    my $c = $coeff->[$i][$j] + $one_half;
    my $d = $c->numerator() / $c->denominator();
    my $r = Math::BigRat->new($d);
    if (rcmp($r, $c) > 0) {
      --$d;
      $r = Math::BigRat->new($d);
    }
    my ($u, $v) = @$m[$i, $j];
    foreach my $k (0..$#$u) {
      $u->[$k] -= $d * $v->[$k];
    }
    foreach my $k (0..$j-1) {
      $coeff->[$i][$k] -= $r * $coeff->[$j][$k];
    }
    $coeff->[$i][$j] -= $r;
  }
} # reduce_one

sub dot_product {
  my ($x, $y) = @_;
  my $sum = $x->[0] * $y->[0];
  foreach my $i (1..$#$x) {
    $sum += $x->[$i] * $y->[$i];
  }
  return $sum;
} # dot_product

# comparison not implemented in current Math::BigRat
sub rcmp {
  my ($r, $s) = @_;
  my $t = $r->numerator() * $s->denominator();
  my $u = $s->numerator() * $r->denominator();
  return $t->bcmp($u);
} # rcmp