#!/usr/bin/perl

sub general_chinese {
    use integer;
    my ($bases, $remainders) = @_;
    my ($i, $j);
  BASE: for ($i = 0; $i < lcm(@$bases); $i++) {
          for ($j = 0; $j < @$bases; $j++) {
               next BASE unless $i % $bases->[$j] == $remainders->[$j];
          }
          return $i;
    }
    return;
}

print "There are ", general_chinese( [2..7], [1,1,1,1,1,0] ), " eggs.\n";

# $lcm = lcm( @numbers ) computes the least common multiple of @numbers.
#
sub lcm {
    use integer;
    my $lcm = shift;
    foreach (@_) { $lcm *= $_ / gcd($_, $lcm) }
    return $lcm;
}

sub gcd {
    use integer;
    my $gcd = shift || 1;
    while (@_) {
        my $next = shift;
        # ($gcd, $next) = ($next, $gcd % $next) while $next;
        while( $next ) {
            my $r = $gcd % $next;
            $r += $next if $r < 0;   # fix for % in integer mode
            $gcd = $next;
            $next = $r;
        }
    }
    return $gcd;
}

# Store the list of moduli.
sub set_chinese {
    use integer;
    @bases = @_;
    $prod  = 1;

    foreach (@_) {

        # Make sure that the bases are all relatively prime,
        # as required by the theorem:
        die 'not relatively prime' unless gcd($_, $prod) == 1;

        # all clear, add it to the product
        $prod *= $_;
    }
    @inverses = map {
                          my $k = $prod / $_;
                          $k * mod_inverse( $k, $_ );
                    } @_;
}

# Convert from a list of remainders into an integer.
sub from_chinese {
    use integer;
    my $v = shift;
    my $t = 0;
    for (0..$#bases) {
        $t += $inverses[$_] * $v->[$_];
    }
    return $t % $prod;
}

# Convert from an integer into a list of remainders.
sub to_chinese {
    use integer;
    my $v = shift;
    my @v = map { $v%$_ } @bases;
    return \@v;
}

set_chinese(3, 4, 5, 7);
print "There are ", from_chinese( [1,1,1,0] ), " eggs.\n";     # prints 301

# Add two Chinese remainder lists.
sub add_chinese {
    use integer;
    my ($v1, $v2) = @_;

    my @v = map
            { ($v1->[$_] + $v2->[$_]) % $bases[$_] }
        0 .. $#bases;

    return \@v;
}

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

