discrete_logarithm_pollard_rho.pl

#!/usr/bin/perl

# Pohlig-Hellman with Pollard's rho for each prime-power factor.

# Pollard's rho algorithm for logarithms
# https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm_for_logarithms

use 5.036;
use ntheory qw(:all);

# Pollard's rho for discrete logarithm in a group of prime order
sub _pollard_rho_log($g, $h, $p, $n, $max_tries = 10) {

    # Trivial cases
    return 0 if ($h == 1);
    return 1 if ($g == $h);

    # For very small prime orders, brute force is simpler and reliable
    if ($p <= 100) {
        my $t = 1;
        for my $i (0 .. $p - 1) {
            return $i if $t == $h;
            $t = mulmod($t, $g, $n);
        }
        return undef;
    }

    foreach my $attempt (1 .. $max_tries) {

        # Random starting point (a,b) with X = g^a * h^b
        my $a1 = urandomm($p);
        my $b1 = urandomm($p);
        my $x1 = mulmod(powmod($g, $a1, $n), powmod($h, $b1, $n), $n);

        my $a2 = $a1;
        my $b2 = $b1;
        my $x2 = $x1;

        # Floyd's cycle detection
        my $iter = sub($a, $b, $x) {
            my $r = ($x % 3);
            if ($r == 0) {
                $a = addmod($a, 1, $p);
                $x = mulmod($x, $g, $n);
            }
            elsif ($r == 1) {
                $b = addmod($b, 1, $p);
                $x = mulmod($x, $h, $n);
            }
            else {
                $a = mulmod(2,  $a, $p);
                $b = mulmod(2,  $b, $p);
                $x = mulmod($x, $x, $n);
            }
            return ($a, $b, $x);
        };

        while (1) {

            # Tortoise step
            ($a1, $b1, $x1) = $iter->($a1, $b1, $x1);

            # Hare step (two iterations)
            ($a2, $b2, $x2) = $iter->($a2, $b2, $x2);
            ($a2, $b2, $x2) = $iter->($a2, $b2, $x2);

            if ($x1 == $x2) {

                # Collision: g^{a1} h^{b1} = g^{a2} h^{b2}
                my $da = submod($a1, $a2, $p);
                my $db = submod($b2, $b1, $p);

                if ($db == 0) {
                    last;    # degenerate case, restart
                }

                my $x = mulmod($da, invmod($db, $p), $p);

                if (powmod($g, $x, $n) == $h) {
                    return $x;
                }
                last;    # verification failed, restart
            }
        }
    }

    return undef;    # failed after max_tries
}

# Solve g^x = a (mod n) where g has order exactly p^e * r,
# and we want x modulo p^e.
sub _prime_power_log($a, $g, $n, $p, $e, $full_order) {

    my $L = $full_order;
    my $r = divint($L, powint($p, $e));    # co-factor

    # Move into the subgroup of order p^e
    my $g0 = powmod($g, $r, $n);
    my $a0 = powmod($a, $r, $n);

    my $x     = 0;
    my $cur_g = $g0;                       # current generator, order p^{e-i}
    my $cur_a = $a0;                       # current element
    my $f     = 1;                         # current digit multiplier

    my $sub_g = powmod($g0, powint($p, $e - 1), $n);    # generator of order p

    foreach my $i (0 .. $e - 1) {

        # Create an element of order p by raising to p^{e-1-i}
        my $exp   = powint($p, $e - $i - 1);
        my $sub_a = powmod($cur_a, $exp, $n);           # corresponding element

        # Solve the discrete log in the prime-order subgroup
        my $d = _pollard_rho_log($sub_g, $sub_a, $p, $n) // return undef;

        $x = addint($x, mulint($d, $f));
        $f = mulint($f, $p);

        # Remove the already found part
        $cur_a = mulmod($cur_a, powmod($cur_g, -$d, $n), $n);
        $cur_g = powmod($cur_g, $p, $n);                        # next generator, order p^{e-1-i}
    }

    return $x;
}

# Solve g^x = a (mod n) where gcd(g, n) = 1, using Pohlig-Hellman over order of g.
# Suitable when n is a prime power (or when called per prime-power factor of n).
sub _dlog_coprime_prime_power_mod($a, $g, $n) {

    my $order = znorder($g, $n) // return undef;

    # Quick necessary condition: a must lie in the subgroup generated by g
    if (powmod($a, $order, $n) != 1) {
        return undef;
    }

    # Trivial case
    if ($order == 1) {
        return ($a == 1 ? 0 : undef);
    }

    # Factor the order into prime powers and solve for each
    my @factors  = factor_exp($order);
    my @residues = ();

    foreach my $pp (@factors) {
        my ($p, $e) = @$pp;
        my $x = _prime_power_log($a, $g, $n, $p, $e, $order) // return undef;
        push @residues, [$x, powint($p, $e)];
    }

    # Combine via CRT
    my $x = chinese(@residues);

    # Verify
    (defined($x) && powmod($g, $x, $n) == $a) ? $x : undef;
}

sub discrete_log($a, $g, $n) {

    # Normalise inputs
    $a = modint($a, $n);
    $g = modint($g, $n);

    # Handle non-coprime case: gcd(g, n) > 1
    if (gcd($g, $n) != 1) {

        my $g_pow = 1;     # g^k mod n (original n), for direct equality check
        my $n_red = $n;    # modulus being reduced
        my $a_red = $a;    # target being reduced
        my $d_acc = 1;     # accumulated product: (g/D_1)*(g/D_2)*...*(g/D_k) mod n_red
        my $k     = 0;

        while (gcd($g, $n_red) != 1) {

            # Check if g^k already equals a (mod n)
            return $k if $g_pow == $a;

            my $D = gcd($g, $n_red);
            return undef if $a_red % $D != 0;

            $n_red = divint($n_red, $D);
            $a_red = divint($a_red, $D);
            $d_acc = mulmod($d_acc, divint($g, $D), $n_red);
            $k++;
            $g_pow = mulmod($g_pow, $g, $n);
        }

        # Final direct check after stripping
        return $k if $g_pow == $a;

        # Phase 2: gcd(g, n_red) = 1 now; solve g^y = a_red * inv(d_acc) (mod n_red)
        my $inv_d = invmod($d_acc, $n_red) // return undef;
        my $a_new = mulmod($a_red, $inv_d, $n_red);
        my $y     = discrete_log($a_new, $g, $n_red);
        return defined($y) ? $k + $y : undef;
    }

    # Coprime case: gcd(g, n) = 1

    # Factor n into prime powers
    my @n_factors = factor_exp($n);

    # Composite n: solve g^x = a (mod p^e) for each prime-power factor, then CRT
    my @residues = ();

    foreach my $pp (@n_factors) {
        my ($p, $e) = @$pp;
        my $pe  = powint($p, $e);
        my $g_i = modint($g, $pe);
        my $a_i = modint($a, $pe);

        my $r = _dlog_coprime_prime_power_mod($a_i, $g_i, $pe);
        return undef unless defined $r;

        my $ord_i = znorder($g_i, $pe) // return undef;
        push @residues, [$r, $ord_i];
    }

    # Combine via CRT
    my $x = chinese(@residues) // return undef;

    # Verify the result
    (powmod($g, $x, $n) == $a) ? $x : undef;
}

use Test::More tests => 1309;

is(discrete_log(5678, 5, 10007), 8620);

foreach my $test (
                  [[5675,              5,      10000019],          2003974],            # 5675 = 5^2003974 mod 10000019
                  [[18478760,          5,      314138927],         34034873],
                  [[553521,            459996, 557057],            15471],
                  [[7443282,           4,      13524947],          6762454],
                  [[32712908945642193, 5,      71245073933756341], 5945146967010377],
  ) {
    my ($t, $v) = @$test;
    say "Testing: discrete_log(", join(', ', @$t), ") = ", $v;
    is(discrete_log($t->[0], $t->[1], $t->[2]), $v);
}

is_deeply(
          [map { discrete_log(powint(2, $_) - 5, 3, powint(2, $_ + 1)) } 0 .. 35],
          [undef,  0,       undef,    1,        7,        3,         27,       43,        75,        139,        11,         779,
           267,    1291,    3339,     7435,     32011,    48395,     81163,    146699,    277771,    15627,      1588491,    2637067,
           539915, 4734219, 13122827, 63454475, 29900043, 231226635, 97008907, 902315275, 365444363, 1439186187, 3586669835, 7881637131
          ]
         );

is_deeply([map { discrete_log(-1, 3, powint(3, $_) - 2) // 0 } 2 .. 30],
          [3, 10, 39, 60, 121, 0, 117, 4920, 0, 0, 0, 28322, 0, 1434890, 0, 0, 0, 116226146, 0, 0, 15690529803, 0, 108443565, 66891206007, 0, 0, 0, 0, 0]);

# Non-coprime tests
is(discrete_log(36, 44, 50), 2);    # 44^2 = 1936 = 36 (mod 50), gcd(44,50)=2
is(discrete_log(0,  2,  4),  2);    # 2^2 = 4 = 0 (mod 4)
is(discrete_log(4,  6,  8),  2);    # 6^2 = 36 = 4 (mod 8)

# Composite modulus, coprime base
is(discrete_log(130, 85, 177), 15);    # 177 = 3*59, gcd(85,177)=1
is(discrete_log(100, 52, 209), 10);    # 209 = 11*19, 52^10 = 100 (mod 209)

# Verify no-solution cases still return undef
is(discrete_log(3, 4, 6), undef);      # no solution exists

is(discrete_log(1, 2, 7), 0);
is(discrete_log(2, 2, 7), 1);
is(discrete_log(4, 2, 7), 2);
is(discrete_log(1, 3, 7), 0);

is(discrete_log(3, 2, 5), 3);          # 2^3 mod 5 = 3
is(discrete_log(4, 2, 5), 2);

is(discrete_log(2,     4,     7),      2);
is(discrete_log(4,     5,     7),      2);
is(discrete_log(5,     3,     7),      5);
is(discrete_log(130,   85,    177),    15);
is(discrete_log(79,    92,    129),    2);
is(discrete_log(115,   116,   141),    26);
is(discrete_log(67741, 90737, 120309), 146);
is(discrete_log(12,    42,    122),    13);
is(discrete_log(36,    44,    50),     2);
is(discrete_log(34,    170,   187),    5);

# Small modulus cycles

is(discrete_log(8, 2, 11), 3);
is(discrete_log(5, 2, 11), 4);
is(discrete_log(9, 3, 11), 2);

# Edge cases

is(discrete_log(1, 1, 13), 0);
is(discrete_log(1, 5, 13), 0);

# g == a
is(discrete_log(7, 7, 19), 1);

# modulus 2
is(discrete_log(1, 1, 2), 0);

# Non-prime modulus

is(discrete_log(4, 2, 15), 2);    # 2^2 = 4 mod 15
is(discrete_log(1, 4, 9),  0);

# Cases where solution may not exist

is(discrete_log(3, 4, 7), undef);
is(discrete_log(3, 2, 4), undef);
is(discrete_log(6, 4, 8), undef);

# Verify correctness by recomputing power

for my $n (7, 11, 13, 17) {
    for my $g (2 .. $n - 1) {
        for my $k (0 .. $n - 1) {

            my $a = powmod($g, $k, $n);
            my $r = discrete_log($a, $g, $n);

            ok(defined($r), "discrete_log($a, $g, $n)");
            is(powmod($g, $r, $n), $a) if defined($r);
        }
    }
}

# Randomized tests

for (1 .. 100) {
    my $n = urandomm(200000 - 50000) + 50000;
    my $g = urandomm($n - 2) + 2;
    my $k = urandomm(50000);

    my $a = powmod($g, $k, $n);
    my $r = discrete_log($a, $g, $n);

    ok(defined($r), "discrete_log($a, $g, $n)");
    is(powmod($g, $r, $n), $a) if defined($r);
}

# Computationally intensive tests

my $p = 1000003;
my $g = 2;
my $k = 123456;

my $a = powmod($g, $k, $p);

is(powmod($g, discrete_log($a, $g, $p), $p), $a);

# Larger exponent

my $k2 = 654321;
my $a2 = powmod($g, $k2, $p);

is(powmod($g, discrete_log($a2, $g, $p), $p), $a2);

# Large prime modulus stress test

my $p2 = 10000019;
my $g2 = 2;
my $k3 = 777777;

my $a3 = powmod($g2, $k3, $p2);

is(powmod($g2, discrete_log($a3, $g2, $p2), $p2), $a3);