package Crypt::Curve25519;
use base Math::Modular;
use bigint;

# TODO: constant timing code

# See also: http://cr.yp.to/ecdh.html

our @exclude=(0,
	      1,
	      325606250916557431795983626356110631294008115727848805560023387167927233504,
	      39382357235489614581723060781553021112529911719440698176882885853963445705823,
	      (2**255) - 19 - 1,
	      2**255 - 19, 2**255 - 19 + 1,
	      2**255 - 19 + 325606250916557431795983626356110631294008115727848805560023387167927233504,
	      2**255 - 19 + 39382357235489614581723060781553021112529911719440698176882885853963445705823,
	      2*(2**255 - 19) - 1,
	      2*(2**255 - 19),
	      2*(2**255 - 19) + 1);

our $a = 486662;
our $a_ = ($a - 2) / 4;
our $p = (2**255)-19;

sub makemod {
  my $x = shift;
  if (UNIVERSAL::isa($x, Math::BigInt)) {
    $x = $x->copy;
  } else {
    $x = Math::BigInt->new($x);
  }
  return $x->modular($p);
}

sub make {
  makemod(shift);
}

sub cast {
  my $x = shift;
  return $x if (UNIVERSAL::isa($x, __PACKAGE__));
  return make($x);
}

# sub curve25519 {
#   my ($n, $x1) = @_;
#   return mul(cast($x1), cast($n))
#   #return cast($x1) * cast($n);
# }

sub splitbits {
  my $i = shift()->copy->unmodular; @_ = ();
  my $r;
  while ($i > '0') {
    ($i, $r) = $i->bdiv(2);
    unshift @_, ($r == 1)?'1':'0';
  }
  return (@_);
}

sub doublexz {
  my ($x, $z) = @_;
  my $k = $x + $z;
  my $l = $x - $z;
  my $m = $k**2;
  my $n = $l**2;
  my $o = $m - $n;
  return ($n * $m, $o * ($m + $a_ * $o));
}

sub sum {
  my ($x, $z, $xp, $zp, $x1) = @_;
  my $p = ($x - $z) * ($xp + $zp);
  my $q = ($x + $z) * ($xp - $zp);
  return ((($p + $q) ** 2),
	  (($p - $q) ** 2) * $x1);
}

sub curve25519 {
  my ($f, $c) = @_;

  ($c, $f) = (cast($c), cast($f));
  return make(0) unless $f > '0';

  $f = makemod($f);

  my ($x1, $z1) = (makemod($c), makemod(1));
  my ($xa, $za) = (doublexz($x1, $z1));
  my ($x, $z) = ($x1, $z1);
  my $g = 1;
  my @b = splitbits($f);
  shift @b;

  my $b;
  for $b (@b) {
    if ($b == 0) {
      ($x, $z, $xa, $za) = (doublexz($x, $z), sum($xa, $za, $x, $z, $x1));
      $g = $g + $g;
    } else {
      ($x, $z, $xa, $za) = (sum($xa, $za, $x, $z, $x1), doublexz($xa, $za));
      $g = $g + $g + 1;
    }
  }

  warn "internal computation error: n=$n, f=$f" unless $g == $f;

  $x = $x * ($z ** ($p - 2));
  return make($x);
}


!undef;
