# Toolkit for calculating with bignums and moduli
#
# Most of the routines in this package have been copied from two sources:
# "Mastering Algorithms":
#   Mastering Algorithms with Perl, by Jon Orwant, Jarkko Hietaniemi and
#   John Macdonald. Copyright 1999 O'Reilly & Associates, Inc.
#   ISBN: 1-56592-398-7
# "Applied Cryptography":
#   Applied Cryptograpy: Protocols, Algorithms and Source Code in C,
#   by Bruce Schneier. Copyright 1994 John Wiley & Sons, Inc.
#   ISBN: 0-471-59756-2 (page numbers refer to the first edition)

package ModToolkit;

use strict;
use FileHandle;
require Exporter;
use vars qw(@ISA @EXPORT_OK);

@ISA = qw(Exporter);
@EXPORT_OK = qw(
    $bignumlib
    exp_mod
    brand
    is_prime
    bprime
    gcd
    gcd_linear
    mod_inverse
    str2n
    n2str
    rsa_encode
    rsa_decode
    readkey
    blind
    npretty
);

# determine best bignum library. Preloaded lib if present, otherwise
# try Math::GMP first, but fall back to Math::BigInt
use vars qw($bignumlib);

BEGIN {
    if ( %Math::BigInt:: ) {
	$bignumlib = 'Math::BigInt';
    }
    elsif ( %Math::GMP:: ) {
	$bignumlib = 'Math::GMP';
    }
    else {
	eval {
	    require Math::GMP;
	    import Math::GMP;
	};
	if ( ! $@ ) {
	    $bignumlib = 'Math::GMP';
	}
	else {
	    require Math::BigInt;
	    import Math::BigInt;
	    $bignumlib = 'Math::BigInt';
	}
    }
}

# $i ** $j mod $n
# From Mastering Algorithms, chapter 12, page 519.
sub exp_mod {
    use integer;
    my($i, $j, $n) = @_;

    my $result = $i - $i + 1; # 1, but in the same type as $i
    return $result unless $j;

    my $pow2 = $i;

    while ( 1 ) {
	if ( $j % 2 ) {
	    $result = ($pow2 * $result) % $n;
	    return $result unless --$j;
	}
	$j /= 2;
	$pow2 = ($pow2 * $pow2) % $n;
    }
}

# points to os-specific rand routine, supposed to return 1 byte of
# randomness
my $osrand;

# devices or sockets where randomness can be found. /dev/urandom is
# used on freebsd and linux, /var/run/random.md5 is used on BSDi.
my @randdevs = qw(
    /dev/urandom
    /var/run/random.md5
);
my $randdev;
my $randdev_fh;

my($s1, $s2);

sub initrand {
    # prefer a random dev
    for my $dev ( @randdevs ) {
	if ( -r $dev && ! -f $dev ) {
	    $randdev = $dev;
	    $randdev_fh = new FileHandle $dev, "r" or
		die "Cannot open $dev: $!\n";
	    $osrand = \&dev_rand;
	    return;
	}
    }
    # no devices found, use double linear congruential generator
    # initialised with os-specific pseudo-random data.
    if ( $^O ne "MSWin32" ) {
        # unix specific pseudo-random data. Provide an extensive
        # path for finding the utils needed.
        local $ENV{PATH} = $ENV{PATH} . ":" . join(":",
            qw( /sbin /usr/sbin /bin /usr/bin /usr/local/bin /opt/bin
                /usr/local/gnu/bin /usr/contrib/bin /opt/gnu/bin
                /usr/libexec /usr/local/sbin )
        );
        $s1 = unpack("%32L*", `ps wwwaxle 2> /dev/null | gzip -c`);
        $s2 = unpack("%32L*", `netstat -na 2> /dev/null | gzip -c`);
    }
    else {
        # windows specific pseudo-random data.
        $s1 = unpack("%32L*", `netstat -na`);
        $s2 = unpack("%32L*", `netstat -se`);
    }
    # these values are available for all systems (but they are not
    # very random, unfortunately).
    $s1 ^= unpack( "%32L*", join("", %ENV, times, $^T) );
    $s2 ^= unpack( "%32L*", join("", $^O, $$, $^X, time) );
    $s1 = -$s1 if $s1 < 0;
    $s2 = -$s2 if $s2 < 0;
    $s1 %= 2_147_483_563;
    $s2 %= 2_147_483_399;
    $s1 ||= 1;
    $s2 ||= 1;
    $osrand = \&double_congruential;
}

sub dev_rand {
    my $byte;
    sysread($randdev_fh, $byte, 1) == 1 or
	die "Cannot read from $randdev: $!\n";
    unpack("C", $byte);
}

# the routines modmult and double_congruential are based on the C code in
# Applied Cryptography, chapter 15.1, page 349.

# calculates s = s*b mod m, provided m=a*b + c, where 0 <= c < b
sub modmult {
    use integer;

    my($sr, $a, $b, $c, $m) = @_;
    my $q = $$sr / $a;
    $$sr = $b * ( $$sr - $a * $q) - $c * $q;
    $$sr += $m if $$sr < 0;
}

sub double_congruential {
    use integer;

    modmult(\$s1, 53668, 40014, 12211, 2_147_483_563);
    modmult(\$s2, 52774, 40692,  3791, 2_147_483_399);
    my $z = $s1 - $s2;
    $z += 2_147_483_562 if $z < 0;
    $z % 256;
}

# rand($i), but return bignum if input is bignum
sub brand {
    initrand() unless $osrand;

    use integer;
    my $i = shift;

    my($r, $v, $max);
    do {
        $r = $i - $i;	# 0, but of the same type as $i
        $v = $r;
        my $x = 2 * $i; # be generous about number of bytes needed.
        while ( $x != 0 ) {
            my $byte = int rand 256;    # get perl's idea of a random byte
            $byte ^= $osrand->();       # add os-specific randomness to it
            $r = ($r * 256) + $byte;
            $v = ($v * 256) + 255;      # max possible random generated
            $x = int($x / 256);
        }
        ### max = greatest multiple of $i <= $v;
        $max = int(($v+1) / $i) * $i;
    }
    while ( $r >= $max );
    $r %= $i;
    $r += $i if $r < 0;
    $r;
}

# is_prime is from Mastering Algorithms, chapter 12, page 520.
# Slightly modified by JohnPC:
#  - added "quick" tests
#  - added auto-conversion to bignum for numbers > sqrt(2**31)
#  - added fix for core dump of Math::GMP
#  - added cluttery progress indicator on STDERR

# check whether $n is prime, by trying at least 5 random tests.
# This is the rabin-miller test.
sub is_prime {
    use integer;
    my $n = shift;
    return 1 if $n <= 2; # 0, 1 and 2 are considered prime
    # first test some "easy" cases
    return 0 unless $n % 2 && $n % 3 && $n % 5 && $n % 7 && $n % 11 &&
	$n % 13 && $n % 17 && $n % 19 && $n % 23 && $n % 29 && $n % 31;
    if ( !ref $n && $n > 46340 ) {
	# this is unsafe for numbers > sqrt(2**31). Convert to bignum
	$n = $bignumlib->new($n);
    }
    my $n1 = $n - 1;
    my $one = $n - $n1;		# 1, but of the same type as $n

    my $witness = $one * 100;

    # find the power of two for the top bit of $n1.
    my $p2 = $one;
    my $p2index = -1;
    ++$p2index, $p2 *= 2 while $p2 <= $n1;
    $p2 /= 2;

    # number of iterations: 5 for 260-bit number, go up to 25 for
    # much smaller numbers.
    my $last_witness = 5;
    $last_witness += (260 - $p2index)/13 if $p2index < 260;
    for my $witness_count ( 0..$last_witness ) {

	print STDERR $witness_count ? "+" : "";

	$witness *= 1024;
	# note: adding fractions to a Math::GMP makes it dump core
	$witness += int rand(1024);
	$witness = $witness % $n if $witness > $n;
	$witness = $one * 100, redo if $witness == 0;

	my $prod = $one;
	my $n1bits = $n1;
	my $p2next = $p2;

	# compute $witness ** ($n - 1), sortof...
	while ( 1 ) {
	    # is $prod, the product so far, a square root of 1? (-1 or +1)
	    my $rootone = $prod == 1 || $prod == $n1;

	    $prod = ($prod * $prod) % $n;

	    # an extra root of 1 disproves the primality.
	    if ( $prod == 1 && !$rootone ) {
		print STDERR "-\n";
		return 0;
	    }

	    if ( $n1bits >= $p2next ) {
		$prod = ($prod * $witness) % $n;
		$n1bits -= $p2next;
	    }
	    last if $p2next == 1;
	    $p2next /= 2;
	}
	if ( $prod != 1 ) {
	    print STDERR "-" if $witness_count;
	    return 0;
	}
    }
    print STDERR "!";
    return 1;
}

# return a prime < $x
sub bprime {
    use integer;
    my $x = shift;
RAND:
    {
	my $p = brand($x);
	$p++ unless $p % 2;
	redo RAND if $p >= $x;
	while ( !is_prime($p) ) {
	    print STDERR ".";
	    $p += 2;
	    redo RAND if $p >= $x;
	}
	print STDERR "\n";
	return $p;
    }
}

# calculate the GCD of $a and $b
# from Mastering Algorithms, chapter 12, page 501.
sub gcd {
    use integer;
    my($a, $b) = @_;

    while ( $b != 0 ) {
	my $r = $a % $b;
	$r += $b if $r < 0;

	$a = $b;
	$b = $r;
    }

    $a;
}

# ($gcd, $afactor, $bfactor) = gcd_linear($a, $b);
# computes GCD of $a and $b, as well as $afactor and $bfactor such that
#    $gcd == $a * $afactor + $b * $bfactor
# from Mastering Algorithms, chapter 12, page 503
sub gcd_linear {
    use integer;
    my($a, $b) = @_;

    # if either is zero, the other is trivially the GCD
    return ($a, 1, 0) if $b == 0;
    return ($b, 1, 0) if $a == 0;

    my($x1, $x2, $y1, $y2) = (0, 1, 1, 0);

    # if we remember the original value of $a and $b, calling them
    # A and B, then the following relations will be maintaind for
    # each iteration:
    #		$a == A * $x2 + B * $y2
    #		$b == A * $x1 + B * $y1

    while ( 1 ) {
	# we need the quotient as well as the remainder
	my ($q, $r);
	$r = $a % $b;
	# int % can do the wrong thing, fix it
	$r += $b if $r < 0;
	$q = ($a - $r) / $b;

	# when the remainder is zero, $b contains the GCD.
	# The relations that we have been maintaining tell us
	# that $b == A * $x1 + B * $y1.
	return ($b, $x1, $y1) if $r == 0;

	# when the remainder is not yet zero, progress to the
	# next iteration, but maintain the relations.
	($a, $b) = ($b, $r);
	($x1, $x2) = ($x2 - $q * $x1, $x1);
	($y1, $y2) = ($y2 - $q * $y1, $y1);
    }
}

# find $ans = mod_inverse($k, $n) such that:
#   ($k * $ans) is 1 (mod $n)
# from Mastering Algorithms, chapter 12, page 515
sub mod_inverse {
    use integer;
    my ($k, $n) = @_;

    my ($d, $kf, $nf) = gcd_linear($k, $n);

    # $d == $kf * $k + $nf * $n == ($kf * $k mod $n)
    return 0 unless $d == 1;
    $kf %= $n;
    $kf += $n if $kf < 0;

    return $kf;
}

# convert a string to a bignum, by treating the string as a very long binary
sub str2n {
    my $str = shift;
    my $n = $bignumlib->new(0);

    for my $c ( unpack("C*", $str) ) {
	$n = ($n * 256) + $c;
    }
    $n;
}

# convert a bignum to a string, by treating the int as a bytestream.
sub n2str {
    my $n = shift;
    my $str = '';

    while ( $n != 0 ) {
	my $r = $n % 256;
	$str = chr($r) . $str;
	$n = ($n - $r) / 256;
    }
    $str;
}

# encode $str using RSA, with key $e mod $n. Actually RSA encoding and
# decoding uses the same algorithm, but this also uses str2n to encode a
# string.
# returns an array of messages if the message contains more bits than $n.
# The RSA algorithm is described in Mastering Algorithms in chapter 13,
# around page 551, and in Applied Cryptography in chapter 12.4, around
# page 282.
sub rsa_encode {
    my($n, $e, $str) = @_;

    my $m = str2n($str);
    my @c;
    while ( $m != 0 ) {
	my $r = $m % $n;
	push @c, exp_mod($r, $e, $n);
	$m = ($m - $r) / $n;
    }
    @c;
}

# the other way around, decode with key $d mod $n all other message args,
# and return a single string.
sub rsa_decode {
    my $n = shift;
    my $d = shift;

    my $m = 0;
    while ( my $c = shift ) {
	$m = ($m * $n) + exp_mod($c, $d, $n);
    }
    n2str($m);
}

# reads a public and possibly private key from file. Returns the
# modulus, public exponent, and if available the private exponent
sub readkey {
    my $base = shift;

    my %k;
    my $err;
    for my $ext ( qw/pubkey seckey/ ) {
	if ( open KEY, "$base.$ext" ) {
	    while ( <KEY> ) {
		chomp;
		if ( /^([ned]): (\+?\d+)$/ ) {
		    $k{$1} = $bignumlib->new($2);
		}
		else {
		    warn "Bad entry in $base.$ext: $_\n";
		    return;
		}
	    }
	    close KEY;
	}
	else {
	    $err ||= "Cannot open $base.$ext: $!\n";
	}
    }
    unless ( $k{n} ) {
	warn $err;
	return;
    }
    ($k{n}, $k{e}, $k{d});
}

# blind a message m, that is supposed to be signed later by the private
# key belonging to pubkey e, n.
# returns the blinded message, signature blinding key $ks, and the
# message blinding key $km.
# to unblind the signature, multiply sig by $ks ** -1, and to unblind
# the message, multiply message by $km ** -1
# The algorithm for blind signatures based on RSA keys is described in
# Applied Cryptography, chapter 16.13, page 404.
sub blind {
    my($n, $e, $m) = @_;

    my $ks = brand($n);
    my $km = exp_mod($ks, $e, $n);
    return (($m * $km) % $n, $ks, $km);
}

# format a number in a somewhat more pretty way
sub npretty {
    my($prefix, $n, $EOL) = @_;

    $EOL ||= "\n";
    my $str = $prefix . $n;
    my $r = '';
    my $linelen = 68;
    while ( length($str) > $linelen ) {
	$r .= substr($str, 0, $linelen) . "\\$EOL" . " " x 4;
	substr($str, 0, $linelen) = '';
	$linelen = 64;
    }
    $r .= $str;
    $r;
}

1;