Prime Power Predicate
March 13, 2015
I discovered Cohen’s algorithm while updating some of my prime number code and translating it to Python, so I’ll give the solution in Python:
from random import randint
from fractions import gcd
We begin with a version of the Miller-Rabin primality checker that returns 0 if the input number n is prime or the witness a that proves its compositeness:
def findWitness(n, k=5): # miller-rabin
s, d = 0, n-1
while d % 2 == 0:
s, d = s+1, d/2
for i in range(k):
a = randint(2, n-1)
x = pow(a, d, n)
if x == 1 or x == n-1: continue
for r in range(1, s):
x = (x * x) % n
if x == 1: return a
if x == n-1: break
else: return a
return 0
Now we have Cohen’s algorithm. The key to the algorithm is the variable q, which is initially n. At each step we find a witness a to the compositeness of the current q, stopping immediately if q is prime, then compute the greatest common divisor d described by Cohen. If that gcd is 1 or q, then Fermat’s Little Theorem tells us that n is not a prime power (because it is not a non-trivial divisor), and we can quit. Otherwise, we set q to the gcd d and make another step of the algorithm. Here’s the code:
# returns p,k such that n=p**k, or 0,0
# assumes n is an integer greater than 1
def primePower(n):
def checkP(n, p):
k = 0
while n > 1 and n % p == 0:
n, k = n / p, k + 1
if n == 1: return p, k
else: return 0, 0
if n % 2 == 0: return checkP(n, 2)
q = n
while True:
a = findWitness(q)
if a == 0: return checkP(n, q)
d = gcd(pow(a,q,n)-a, q)
if d == 1 or d == q: return 0, 0
q = d
We could do a better job in the checkP function, which computes the logarithm by repeated division, but Cohen did it that way, and left it as an exercise to the reader (Chapter 1, Exercise 4) to compute the logarithm by binary search, and we do the same. Here are some examples:
>>> primePower(1940255363782247**37)
(1940255363782247L, 37)
>>> primePower(7)
(7, 1)
The answer comes back instantaneously, much quicker on such large numbers than our previous algorithm.
You can run the program at http://ideone.com/cNzQYr.
Programming Praxis 
Port of the given python code to Java:
However this solution has two drawbacks:
1) The input is of type “long”, so it cannot compute for big numbers like 1940255363782247**37
2) It does *not* compute logarithm by binary search but by repeated division only.
I am working on the proper solution with these drawbacks addressed. I will post it if and when i crack it :)
randLong method is copied shamelessly from http://stackoverflow.com/a/2546186
package com.experiment; import java.math.BigInteger; import java.util.Random; public class MillerRabinPrimTestFinal { public static void main(String[] args) { //System.out.println(findWitness(1940255363782247L)); primePower(7L); primePower(25L); primePower(271818611107L); primePower(505447028499293771L); } private static long findWitness(long n) { long d = n-1; int s=0; long a, x; BigInteger xb; Random rng = new Random(); int k = 5; while(d%2==0) { s = s+1; d = d/2; } witnessloop: for(int i=0;i<k;i++) { a = randLong(rng, 2L, n-2); x = new BigInteger(String.valueOf(a)).modPow(new BigInteger(String.valueOf(d)), new BigInteger(String.valueOf(n))).longValue(); if(x == 1 || x == (n-1)) continue; for(int r=1;r<=s;r++) { xb = new BigInteger(String.valueOf(x)); xb = xb.multiply(xb).mod(new BigInteger(String.valueOf(n))); x = xb.longValue(); if(x == 1) return a; if(x == (n-1)) continue witnessloop; } return a; } return 0L; } private static void checkP(long n, long p) { int k = 0; while(n>1 && n%p == 0) { n = n/p; k = k+1; } if(n==1) { System.out.println(p+","+k); } else { System.out.println("0,0"); } } private static void primePower(long n) { if(n%2 == 0) { checkP(n, 2L); return; } long q = n; while (true) { long a = findWitness(q); if (a == 0) { checkP(n, q); return; } BigInteger ab = new BigInteger(String.valueOf(a)); BigInteger qb = new BigInteger(String.valueOf(q)); long d = ab.modPow(qb, new BigInteger(String.valueOf(n))) .subtract(ab).gcd(qb).longValue(); if (d == 1 || d == q) { System.out.println("0,0"); } q = d; } } private static long randLong(Random rng, long min, long max) { long bits, val; do { bits = (rng.nextLong() << 1) >>> 1; val = bits % max; } while(bits-val+(max-1) < min); return val; } }