diff --git a/imp/extern/__init__.py b/imp/extern/__init__.py
deleted file mode 100644
index e69de29..0000000
diff --git a/imp/math/numbers/__init__.py b/imp/math/numbers/__init__.py
index 22a0d21..f4193b8 100644
--- a/imp/math/numbers/__init__.py
+++ b/imp/math/numbers/__init__.py
@@ -6,7 +6,7 @@ Terminology:
     the "prime proper divisors of n". 
 '''
 
-from imp.extern.primefac import primefac
+from imp.math.primefac import primefac
 
 def factors(n: int) -> int:
     pfactors: list[tuple[int, int]] = []
diff --git a/imp/math/numbers/functions.py b/imp/math/numbers/functions.py
index 8addd51..6cafa7f 100644
--- a/imp/math/numbers/functions.py
+++ b/imp/math/numbers/functions.py
@@ -1,3 +1,5 @@
 def factorial(n: int) -> int:
     if n == 0: return 1
     return n * factorial(n-1)
+
+def
diff --git a/imp/extern/primefac.py b/imp/math/primefac.py
similarity index 100%
rename from imp/extern/primefac.py
rename to imp/math/primefac.py
diff --git a/imp/math/primes.py b/imp/math/primes.py
index b130a0c..eaeebcb 100644
--- a/imp/math/primes.py
+++ b/imp/math/primes.py
@@ -1,10 +1,5 @@
-from math import gcd, inf
-
-from imp.math.numbers import bigomega, factors
-from imp.extern.primefac import (
-    isprime,
-    primegen as Primes,
-)
+from math import gcd
+from imp.math.numbers import bigomega
 
 def coprime(n: int, m: int) -> bool:
     return gcd(n, m) == 1
@@ -23,25 +18,68 @@ def semiprime(n: int) -> bool:
     '''
     return almostprime(n, 2)
     
-def eulertotient(x: int | list) -> int:
-    '''
-    Evaluates Euler's Totient function.
-    Input: `x: int` is prime factorised by Lucas A. Brown's primefac.py
-           else `x: list` is assumed to the prime factorisation of `x: int`
-    '''
-    pfactors = x if isinstance(x, list) else factors(n)
-    return prod((p-1)*(p**(e-1)) for (p, e) in pfactors)
-# def eulertotient(n: int) -> int:
-#     '''
-#     Uses trial division to compute 
-#     Euler's Totient (Phi) Function.
-#     '''
-#     phi = int(n > 1 and n)
-#     for p in range(2, int(n ** .5) + 1):
-#         if not n % p:
-#             phi -= phi // p
-#             while not n % p:
-#                 n //= p
-#     #if n is > 1 it means it is prime
-#     if n > 1: phi -= phi // n 
-#     return phi
+'''
+Euler's Totient (Phi) Function
+'''
+def totient(n: int) -> int:
+    phi = int(n > 1 and n)
+    for p in range(2, int(n ** .5) + 1):
+        if not n % p:
+            phi -= phi // p
+            while not n % p:
+                n //= p
+    #if n is > 1 it means it is prime
+    if n > 1: phi -= phi // n 
+    return phi
+
+'''
+Tests the primality of an integer using its totient.
+NOTE: If totient(n) has already been calculated 
+      then pass it as the optional phi parameter. 
+'''
+def is_prime(n: int, phi: int = None) -> bool:
+    return n - 1 == (phi if phi is not None else totient(n))
+
+'''
+Prime number generator function.
+Returns the tuple (p, phi(p)) where p is prime
+    and phi is Euler's totient function.
+'''
+def prime_gen(yield_phi: bool = False) -> int | tuple[int, int]:
+    n = 1
+    while True:
+        n += 1
+        phi = totient(n)
+        if is_prime(n, phi=phi):
+            if yield_phi:
+                yield (n, phi)
+            else:
+                yield n
+
+'''
+Returns the prime factorisation of a number.
+Returns a list of tuples (p, m) where p is
+a prime factor and m is its multiplicity.
+NOTE: uses a trial division algorithm
+'''
+def prime_factors(n: int) -> list[tuple[int, int]]:
+    phi = totient(n)
+    if is_prime(n, phi=phi):
+        return [(n, 1)]
+    factors = []
+    for p in prime_gen(yield_phi=False):
+        if p >= n:
+            break
+        # check if divisor
+        multiplicity = 0
+        while n % p == 0:
+            n //= p
+            multiplicity += 1
+        if multiplicity:
+            factors.append((p, multiplicity))
+            if is_prime(n):
+                break
+    if n != 1:
+        factors.append((n, 1))
+    return factors
+            
diff --git a/imp/math/util.py b/imp/math/util.py
index 4dacda2..69631c6 100644
--- a/imp/math/util.py
+++ b/imp/math/util.py
@@ -1,9 +1,8 @@
-from collections.abc import Iterable
 from itertools import chain, combinations
 
 def digits(n: int) -> int:
     return len(str(n))
 
-def powerset(iterable: Iterable) -> Iterable:
+def powerset(iterable):
     s = list(iterable)
     return chain.from_iterable(combinations(s, r) for r in range(len(s)+1))