#!/usr/bin/perl

use strict;
use warnings;

use bignum;

use constant CARD_POS  => 2020;

#use constant DECK_SIZE => 10007;
use constant DECK_SIZE => 119315717514047;
use constant SHUFFLES  => 101741582076661;

# Linear function representing shuffle: f(x) = $func[0] * x + $func[1]
my @func = (1, 0);

while (<>) {
    if (m#^deal into new stack#) {
        # compose (-x - 1) with Ax+B, mod N
        # -(Ax+B) - 1 => -Ax + (-B - 1)
        @func = (-$func[0] % DECK_SIZE, (-1 - $func[1]) % DECK_SIZE);

    } elsif (m#^cut (-?\d+)#) {
        # compose (x - c) with Ax+B, mod N
        # (Ax+B) - c => Ax + (B - c)
        @func = ($func[0], ($func[1] - $1) % DECK_SIZE);

    } elsif (m#^deal with increment (\d+)#) {
        # compose mx with Ax+B, mod N
        # m(Ax+B) => (Am)x + Bm
        @func = (($func[0] * $1) % DECK_SIZE, ($func[1] * $1) % DECK_SIZE);

    } else {
        die "Unrecognized line: $_\n";
    }
}

print "f(x) = $func[0]x + $func[1]\n";

my $fx = ($func[0] * 2019 + $func[1]) % DECK_SIZE;
print "Card 2019 is at position $fx\n";

#
#   f(x) = y mod N, need to find f_inverse to get x = f_inv(y) mod N
# Ax + B = y mod N
#     Ax = (y - B) mod N
#      x = A_inv * (y - B) mod N, A_inv is inverse of A in mod N
#
#my $y = (invmod( $func[0], DECK_SIZE ) * ($fx - $func[1])) % DECK_SIZE;
my $y = (&prime_mod_inverse( $func[0], DECK_SIZE ) * ($fx - $func[1])) % DECK_SIZE;
print "Card at position $fx is $y\n\n";

#
# Recursive function to compose linear function n times.
#
sub apply_n_times {
    my ($func, $n) = @_;

    return ($func->[0], $func->[1])  if ($n == 1);

    # divide and conqueur, get the function for half, h(x) = Ax + B
    my ($A, $B) = &apply_n_times( $func, int( $n / 2 ) );

    #
    # h(h(x)) = (A(Ax + B) + B
    #         = A^2x + AB + B
    #         = A'x + B'
    #
    ($A, $B) = (($A * $A) % DECK_SIZE, ($B * ($A + 1)) % DECK_SIZE);

    #
    # Let f(x) = Cx + D
    #
    # h(h(x)) o f(x) = A'(Cx + D) + B'
    #                = (A'C)x + A'D + B'
    #
    if ($n % 2 == 1) {
        ($A, $B) = (($A * $func->[0]) % DECK_SIZE, $A * $func->[1] + $B);
    }

    return ($A, $B);
}

#
# Calculate the inverse of a in mod n, where a and n are relatively prime,
# so that this holds:
#             a^n = a mod n
#     a * a^(n-1) = a mod n
#         a^(n-1) = 1 mod n
#     a * a^(n-2) = 1 mod n, thus a^(n-2) is the inverse mod n of a
#
# We're using the above linear composer to do the exponentiation, as
# composing (Ax + 0) with itself n times is (A^n)x.
#
sub prime_mod_inverse {
    my ($a, $n) = @_;

    my @pow = &apply_n_times( [$a, 0], $n - 2 );
    return ($pow[0]);
}


# 2423 is 101741582076661 % 10006, if testing with 10007 card deck, this
# should be the position in the first cycle that's congruent, so it
# should be equal to the big shuffle.
foreach my $i (1 .. 8, 2423) {
    my ($A, $B) = &apply_n_times( \@func, $i );

    #$y = (invmod( $A, DECK_SIZE ) * (CARD_POS - $B)) % DECK_SIZE;
    $y = (&prime_mod_inverse( $A, DECK_SIZE ) * (CARD_POS - $B)) % DECK_SIZE;
    print "Card at position ", CARD_POS, " after $i shuffles: $y\n";
}

# Caluculating part 2:
#
#       Ax + B = 2020 mod n
#           Ax = (2020 - B) mod n
#            x = A^-1 * (2020 - B) mod n
my ($A, $B) = &apply_n_times( \@func, SHUFFLES );

#$y = (invmod( $A, DECK_SIZE ) * (CARD_POS - $B)) % DECK_SIZE;
$y = (&prime_mod_inverse( $A, DECK_SIZE ) * (CARD_POS - $B)) % DECK_SIZE;
print "Card at position ", CARD_POS, " after ", SHUFFLES, " shuffles: $y\n";


# Extra Challenge: Find the card that doesn't move.
#       f(x) = x mod n
#     Ax + B = x mod n
#     x - Ax = B mod n
#   x(1 - A) = B mod n
#          x = B * mod_inv(1 - A) mod n

my $x_stable = ($B * &prime_mod_inverse( (1 - $A), DECK_SIZE )) % DECK_SIZE;
$fx = ($func[0] * $x_stable + $func[1]) % DECK_SIZE;

print "Card at position $x_stable is always $fx.\n";