diff --git a/README b/README
index f616f3f..248d5eb 100644
--- a/README
+++ b/README
@@ -1 +1,4 @@
 The "imbaud python library" (imp lib), or just imp for short!
+
+TODO:
+- define a getPrime function like PyCryptodome's
diff --git a/imp/math/crt.py b/imp/math/crt.py
new file mode 100644
index 0000000..a50679c
--- /dev/null
+++ b/imp/math/crt.py
@@ -0,0 +1,124 @@
+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)])
+