#!/usr/bin/perl

# $gcd = euclid( $a, $b ) computes the GCD of $a and $b.
sub euclid {
    use integer;
    my ( $a, $b ) = @_;

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

        # MIDPOINT

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

    return $a;
}

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

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

print euclid(27, 39), "\n";	# prints 3
print gcd(60, 480, 105), "\n";	# prints 15

($gcd, $afactor, $bfactor) = gcd_linear(60, 480, 105);
print "gcd = $gcd, afactor = $afactor, bfactor = $bfactor.\n";