def rootmod(k, b, m, verbose=True):
     # Solve x^k==b mod m
     p = phi(m)
     if verbose: print("phi =", p)
     u = invert(k, p)
     if not u:
          if verbose: print("(k, phi) != 1")
          return
     else:
          if verbose: print("u =", u)
          out = pow(b, u, m)

          if gcd(b, m) > 1:
               if verbose: print("b^(ku-1) == {} mod {}".format(pow(b, k*u-1, m), m))
               if out == 0:
                    if verbose: print("(b,m) > 1, and the algorithm produced 0.")
                    return 0
               elif b % m == pow(out, k, m):
                    if verbose: print("Despite (b,m) > 1, the algorithm produced {}, which is correct.".format(out))
                    return out + m
               else:
                    if verbose: print("(b,m) > 1 and the algorithm produced an incorrect answer {}".format(out))
                    return -1
          else: # Guaranteed to be correct
               return out
     # Proof: Assume b^u is a solution. Then b^(ku) == b mod m.
     # Thus b^(ku-1) == 1 mod m --> phi(m) | ku-1 by Euler's Thm (and (b,m)=1).
     # But this is ku == 1 mod phi(m) by definition. So b^u is a solution
     # if u is k-inverse mod phi(m). (This also implies (k, phi(m)) = 1.)
     # Note: As some of the above if statements might suggest, the algorithm still sometimes works even if (b,m) > 1.

def huh(count=100, hi=100, verbose=False):
     # Perform a rootmod on "count" number of random triplets k,b,m
     from random import randint
     zero, correct, wrong = 0, 0, 0
     tot = count
     while count:
          if verbose: print()
          k, b, m = randint(1,hi), randint(1,hi), randint(1,hi)
          x = rootmod(k, b, m, verbose)
          if x is not None: # Ignore (k, phi) > 1
               count -= 1
               if x > m: # (b,m) > 1, but it worked anyways
                    correct += 1
               elif x == -1:
                    wrong += 1
               elif x == 0:
                    zero += 1
     print("There were {} zero-results, {} wrong results, {} correct results, and {} with (b, m) = 1".format(zero, wrong, correct, tot-zero-wrong-correct))