def kempner(p): for m in range(p+1): if factorial(m) % p == 0: return m return 0 def keller_olson_eta_function(p, k): runningSum = 0 for j in range(3, k+1): runningSum = runningSum + kempner(p ** j) return runningSum def keller_olson_formula(p, k): if k == 1: return factorial(p) else: partial = (factorial(p) * ((p - 1) ** p)) * (p ** p) if k == 2: return partial else: return partial * (p ** keller_olson_eta_function(p, k)) def count_pp_modn(n): factors = list(factor(n)) runningProduct = 1 for f in factors: runningProduct = runningProduct * keller_olson_formula(f[0], f[1]) return runningProduct