#!/usr/bin/perl

# $ans = mod_divide( $n, $d, $p ) uses modular division to
# solve ($ans * $d) == $n (mod $p), where $p is prime.
sub mod_divide {
    use integer;
    my ( $n, $d, $p ) = @_;
    my @gcd = gcd_linear( $d, $p );
    my $inverse = ($n * $gcd[1]) % $p;
    $inverse += $p if $inverse < 0;  # fix if negative
    return $inverse;
}

print mod_divide( 1, 3, 5 ), "\n";  # prints 2

my @bases;
my @inverses;
my $prod;

# Find $ans = mod_inverse( $k, $n ) such that:
#    ($k * $ans) is 1 (mod $n)
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;
}

# ( $gcd, $afactor, $bfactor ) = gcd_linear( $a, $b )
# This computes the greatest common divisor of $a and $b,
# as well as $afactor and $bfactor such that
#           $gcd == $a * $afactor + $b * $bfactor
sub gcd_linear {
    use integer;

    my ( $a, $b ) = @_;

    # If either is zero, the other is trivially the GCD.
    return ( $a, 1, 0 ) unless $b;
    return ( $b, 0, 1 ) unless $a;

    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 maintained 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 ) unless $r;

        # 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);
    }
}