def crt(eqs: list[tuple[int, int]]) -> tuple[int, int]:
    '''
    Solves simultaneous linear diophantine equations
    using the Chinese-Remainder Theorem.
    Example:
    >>> crt([2, 3), (3, 5), (2, 7)])
    (23, 105)
    # aka x is congruent to 23 modulo 105
    '''
    # calculate the quotient and remainder form of the unknown
    q = 1
    r = 0
    # store remainders and moduli for each linear equation
    R = [eq[0] for eq in eqs]
    M = [eq[1] for eq in eqs]
    N = len(eqs)
    for i in range(N):
        (R_i, M_i) = (R[i], M[i])
        # adjust quotient and remainder
        r += R_i * q
        q *= M_i
        # apply current eq to future eqs
        for j in range(i+1, N):
            R[j] -= R_i
            R[j] *= pow(M_i, -1, M[j])
            R[j] %= M[j]
    return r # NOTE: forall integers k: qk+r is also a solution


from math import isqrt, prod

def bsgs(g: int, h: int, n: int) -> int:
    '''
    The Baby-Step Giant-Step algorithm computes the
    discrete logarithm (or order) of an element in a 
    finite abelian group.

    Specifically, for a generator g of a finite abelian
    group G order n, and an element h in G, the BSGS
    algorithm returns the integer x such that g^x = h.
    For example, if G = Z_n the cyclic group order n then:
        bsgs solves g^x = h (mod n) for x.
    '''
    # ensure g and h are reduced modulo n
    g %= n
    h %= n
    # ignore trivial case
    # NOTE: for the Pohlig-Hellman algorithm to work properly
    # NOTE: BSGS must return 1 NOT 0 for bsgs(1, 1, n)
    if g == h: return 1
    m = isqrt(n) + 1 # ceil(sqrt(n))
    # store g^j values in a hash table
    H = {j: pow(g, j, n) for j in range(m)}
    I = pow(g, -m, n)
    y = h
    for i in range(m):
        for j in range(m):
            if H[j] == y:
                return i*m + j
        y = (y * I) % n
    return None # No Solutions

def factors_prod(pf: list[tuple[int, int]]) -> int:
    return prod(p**e for (p, e) in pf)
        
def sph(g: int, h: int, pf: list[tuple[int, int]]) -> int:
    '''
    The Silver-Pohlig-Hellman algorithm for computing
    discrete logarithms in finite abelian groups whose
    order is a smooth integer.

    NOTE: this works in the special case, 
    NOTE: but I can't get the general case to work :(
    '''
    # ignore the trivial case
    if g == h:
        return 1
    R = len(pf) # number of prime factors
    N = factors_prod(pf) # pf is the prime factorisation of N
    print('N', N)
    # Special Case: groups of prime-power order:
    if R == 1:
        (p, e) = pf[0]
        x = 0
        y = pow(g, p**(e-1), N)
        for k in range(e):
            # temporary variables for defining h_k
            w = pow(g, -x, N)
            # NOTE: by construction the order of h_k must divide p
            h_k = pow(w*h, p**(e-1-k), N)
            # apply BSGS to find d such that y^d = h_k
            d = bsgs(y, h_k, N)
            if d is None:
                return None # No Solutions
            x = x + (p**k)*d
        return x

    # General Case:
    eqs = []
    for i in range(R):
        (p_i, e_i) = pf[i]
        # phi = (p_i - 1)*(p_i**(e_i - 1)) # Euler totient
        # pe = p_i**e_i
        # P = (N//(p_i**e_i)) % phi # reduce mod phi by Euler's theorem
        g_i = pow(g, N//(p_i**e_i), N)
        h_i = pow(h, N//(p_i**e_i), N)
        print("g and h", g_i, h_i)
        print("p^e", p_i, e_i)
        x_i = sph(g_i, h_i, [(p_i,e_i)])
        print("x_i", x_i)
        if x_i is None:
            return None # No Solutions
        eqs.append((x_i, p_i**e_i))
        print()
    # ignore the quotient, solve CRT for 0 <= x <= n-1
    x = crt(eqs) # NOTE: forall integers k: Nk+x is also a solution
    print('solution:', x)
    print(eqs)
    return x
    

if __name__ == '__main__':
    result = sph(3, 11, [(2, 1),(7,1)])