2

I am trying to calculate prime numbers on a single machine, around the size of 2^30-2^100.
My algorithm is included below for anyone interested.
I've optimized this Python code to be O(sqrt(n/2)) (I believe) for each number: it only accepts odd numbers, and I ensure the number passed to it is odd in another method.

I used the Fermat primality test to try and speed up the process. However, the numbers are too large for the built-in math.pow() method so I used Exponentiation by Squaring.
However, this is taking a very long time for larger numbers - which would be faster using just brute force.

Is my implementation wrong?
The time comes from the squaring algorithm, its recurrence stack also eats up my memory, is there a faster algorithm for this that I should research?

To calculate if the number 35184372088967 is prime, it took .00100111 seconds using my brute force algorithm, but took .40608 seconds to run the prime test.

Brute force prime number check:

def isPrime(n):
    for i in range(3,int(math.sqrt(n)),2):
        if(n%i==0):
            return False
    return True

Implementation of Fermat's algorithm:

def couldBePrime(n):
        if(n>308):
            return power(2,n-1)%n==1
        else:
            return math.pow(2,n-1)%n==1

Exponentiation by squaring algorithm (The time consuming part):

def power(base,exp):
    if exp == 0:
        return 1
    elif exp == 1:
        return base
    elif (exp & 1) != 0:
        return base * power(base * base, exp // 2)
    else:
        return power(base * base, exp // 2)
boardrider
  • 5,882
  • 7
  • 49
  • 86
CS2016
  • 331
  • 1
  • 3
  • 15
  • There are several problems in your code, some of them related to the scientific modeling of the problem (best handled on CS.SE but also on-topic on SO), some of them related to the Python implementation (off-topic on CS.SE and best handled on SO). So I am migrating to Stack Overflow (please don't repost). – Gilles 'SO- stop being evil' Oct 13 '17 at 16:46
  • 1
    Instead of your power function, say `pow(base,exponent,modulus)` to perform modular exponentiation. It's a built-in function, you don't need to write it yourself. – user448810 Oct 13 '17 at 17:05

1 Answers1

5

Bug: math.pow calculates floating-point values. Floating point calculations are approximate, and will give you meaningless results here. You need integer calculations, such as you do (inefficiently) in your power function. Python's built-in ** operator and pow function (not math.pow, which is a different function) both operate on integers.

In Python, like in many programming languages, the library called math is specifically about floating point computations, not about other kinds of mathematical computations such as calculations done on integers.

Inefficiency: to calculate b^e mod n, it is a lot more efficient to perform arithmetic modulo n, rather than first calculate b^e and then divide the result by n. Calculating b^e requires building a very large number, and this will be slow because the numbers get large pretty quickly as the calculation goes through higher and higher powers of b. (The optimal way to calculate b^e is not easy to determine, but all the ways involve calculating intermediate powers of b, the only practical uncertainty is in what order.) When you want the result modulo n, do all the successive multiplications modulo n: calculate b^2 mod n, then square and reduce modulo n to get b^4 mod n, etc. Each time you perform a multiplication, take the remainder of the division by n before you do anything else.

In Python, the standard library function pow (remember, not math.pow) will do that for you. It's as simple as

def couldBePrime(n):
    return pow(2, n-1, n) == 1

If Python didn't have this function, then your power function a reasonable way to implement it, if you reduced each intermediate result modulo n.

def power(base, exp, mod):
    if exp == 0:
        return 1
    elif exp == 1:
        return base % mod
    elif (exp & 1) != 0:
        return (base * power((base * base) % mod, exp // 2, mod)) % mod
    else:
        return power((base * base) % mod, exp // 2)

Calling a built-in function is a lot faster, of course, both because this is a decent but not extremely good way to perform the operation, and because Python is more optimized for ease of writing than for speed so it's best to leave the as much as possible of the numerical heavy lifting to built-in functions.

An additional note: to calculate a power of two, there's a much faster way than multiplications — do bit shifting. But that doesn't help here because you want to do calculate 2^e mod n and not 2^e.

Gilles 'SO- stop being evil'
  • 104,111
  • 38
  • 209
  • 254
  • This is exactly what I was looking for, not sure how I overlooked the built in `pow` function, but thank you for your very detailed and helpful response! – CS2016 Oct 13 '17 at 17:09