This code finds the primitive roots of a relatively small number. It is not optimized for very large numbers.

Primitive_root.py

from fractions import gcd
import copy



def get_relative_primes(number):
    #   will hold all numbers that are relatively prime to number
    rel_primes = set()  
    # Dictionary:   key is primitive root 
    #               value is the set of powers mod number
    prim_root_dict = {}
    # Dictionary:   key is a relative prime to number
    #               value is the set of powers mod number
    all_powers = {}
    #   will hold the number or relative primes: phi(n)
    phi = 0
    #   get relative primes to the number
    for i in range(1, number):
        if gcd(i, number) is 1:
            rel_primes.add(i)
            phi = phi + 1
    #   iterate through all relative primes to find primitive
    #   roots       
    for i in rel_primes:
        temp_relprimes = copy.deepcopy(rel_primes)
        prim_set = set()
        #   get all rel_primes [0..phi(n)]th power mod number: i^[0..phi(i)] mod number
        for n in range(phi):
            rp = 1;
            #   get to the nth power the efficient way
            for p in range(n + 1):
                rp = (rp * i) % number
            '''
                determine if rp is in the set of relative primes.if so, 
                add n to the primitive set of this relative prime and 
                remove it from the temporary set of relative primes.
            '''
            if rp in temp_relprimes: 
                prim_set.add(rp)
                temp_relprimes.remove(rp)

        all_powers[i] = prim_set
        if len(temp_relprimes) is 0:
            prim_root_dict[i] = prim_set

    
    print("all numbers relatively prime to " + str(number) + ": ")
    print("\t" + str(list(rel_primes)))
    print("phi(" + str(number) + "): " + str(phi))
    print("="*100)
    print("bd"*50)
    print("pq"*50)
    print("="*100)

    print("all powers Z(" + str(number) + "):")
    for i in all_powers:
        print(str(i) + ":\t" + str(list(all_powers[i])))
    print("="*100)
    print("bd"*50)
    print("pq"*50)
    print("="*100)
    count = len(prim_root_dict)
    print("all primitive roots:")
    if count is not 0:
        print("Size:\t" + str(count))
        print("PR:\t" + str(prim_root_dict.keys()) + "\n" + "-"*10)
        for i in prim_root_dict:
            print(str(i) + ":\t" + str(list(prim_root_dict[i])))

    else:
        print("\tNo primitive root")
get_relative_primes(25)