discrete_logarithm_pollard_rho.sf

#!/usr/bin/ruby

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

# Pollard's rho for groups of prime order p
func __prime_order_log__(g, h, p, n, max_tries = 100) {

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

    for attempt in (1..max_tries) {

        # Random starting point (a,b) with X = g^a * h^b
        var a1 = irand(0, p-1)
        var b1 = irand(0, p-1)
        var x1 = mulmod(powmod(g, a1, n), powmod(h, b1, n), n)

        var a2 = a1
        var b2 = b1
        var x2 = x1

        # Floyd's cycle detection
        var iter = { |a,b,x|
            var 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)
            }
            (a, b, x)
        }

        while (true) {
            # 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}
                var da = submod(a1, a2, p)
                var db = submod(b2, b1, p)

                if (db == 0) {
                    break      # degenerate case, restart
                }

                var x = mulmod(da, invmod(db, p), p)

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

    return nil      # 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.
func __prime_power_log__(a, g, n, p, e, full_order) {

    var L = full_order
    var r = idiv(L, ipow(p,e))          # co‑factor

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

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

    for i in (0 ..^ e) {

        # Create an element of order p by raising to p^{e-1-i}
        var exp = ipow(p, e - i - 1)
        var sub_g = powmod(cur_g, exp, n)   # generator of order p
        var sub_a = powmod(cur_a, exp, n)   # corresponding element

        # Solve the discrete log in the prime‑order subgroup
        var d = __prime_order_log__(sub_g, sub_a, p, n) \\ return (nil, false)

        x += (d * f)
        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, true)
}

func discrete_log(a, g, n, order = nil) {

    # Normalise inputs
    a %= n
    g %= n

    # g must be invertible modulo n
    if (!g.is_coprime(n)) {
        return nil
    }

    # Determine the order of g if not provided
    order := znorder(g, n) \\ return nil

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

    # Trivial cases
    if (order == 1) {
        return (a == 1 ? 0 : nil)
    }

    # Factor the order into prime powers
    var factors = factor_exp(order)

    # Solve for x modulo each prime power
    var residues = []

    for p,e in (factors) {
        var (x, ok) = __prime_power_log__(a, g, n, p, e, order)
        ok || return nil
        residues << [x, ipow(p, e)]
    }

    # Combine via CRT
    var x = Math.chinese(residues...)

    # Verify the result (should always hold if the algorithm succeeded)
    (powmod(g, x, n) == a) ? x : nil
}

assert_eq(discrete_log(5678, 5, 10007), 8620)
#assert_eq(discrete_log(232752345212475230211680, 23847293847923847239847098123812075234, 804842536444911030681947), 13)

[
 [[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];
].each_2d {|t, v|
    say "Testing: discrete_log(#{t.join(', ')}) == #{v}"
    assert_eq(discrete_log(t...), v)
}

assert_eq(
    36.of{|n| discrete_log(2**n - 5, 3, 2**(n+1)) }.grep.join(' '),
    %n[NaN, 0, NaN, 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].grep.join(' ')
)