Why classical division ("/") for large integers is much slower than integer division ("//")?

Question

I ran into a problem: The code was very slow for 512 bit odd integers if you use classical division for (p-1)/2. But with floor division it works instantly. Is it caused by float conversion?

def solovayStrassen(p, iterations):    
    for i in range(iterations):         
        a = random.randint(2, p - 1)         
        if gcd(a, p) > 1:
            return False        
        first = pow(a, int((p - 1) / 2), p)    
        j = (Jacobian(a, p) + p) % p                            
        if first != j: 
            return False 
    return True

---

The full code

import random
from math import gcd

#Jacobian symbol
def Jacobian(a, n): 
    if (a == 0): 
        return 0 
    ans = 1 
    if (a < 0):         
        a = -a
        if (n % 4 == 3): 
            ans = -ans 
    if (a == 1): 
        return ans
    while (a): 
        if (a < 0): 
            a = -a
            if (n % 4 == 3):
                ans = -ans 
        while (a % 2 == 0): 
            a = a // 2 
            if (n % 8 == 3 or n % 8 == 5): 
                ans = -ans  
        a, n = n, a
        if (a % 4 == 3 and n % 4 == 3): 
            ans = -ans
        a = a % n
        if (a > n // 2): 
            a = a - n 
    if (n == 1): 
        return ans
    return 0 

def solovayStrassen(p, iterations):    
    for i in range(iterations):         
        a = random.randint(2, p - 1)         
        if gcd(a, p) > 1:
            return False        
        first = pow(a, int((p - 1) / 2), p)    
        j = (Jacobian(a, p) + p) % p                            
        if first != j: 
            return False 
    return True

def findFirstPrime(n, k): 
    while True:              
        if solovayStrassen(n,k):            
            return n
        n+=2    

a = random.getrandbits(512)
if a%2==0:
    a+=1    
print(findFirstPrime(a,100))

Answer 1

As noted in comments, int((p - 1) / 2) can produce garbage if p is an integer with more than 53 bits. Only the first 53 bits of p-1 are retained when converting to float for the division.

>>> p = 123456789123456789123456789
>>> (p-1) // 2
61728394561728394561728394
>>> hex(_)
'0x330f7ef971d8cfbe022f8a'
>>> int((p-1) / 2)
61728394561728395668881408
>>> hex(_) # lots of trailing zeroes
'0x330f7ef971d8d000000000'

Of course the theory underlying the primality test relies on using exactly the infinitely precise value of (p-1)/2, not some approximation more-or-less good to only the first 53 most-significant bits.

As also noted in a comment, using garbage is likely to make this part return earlier, not later:

        if first != j: 
            return False 

So why is it much slower over all? Because findFirstPrime() has to call solovayStrassen() many more times to find garbage that passes by sheer blind luck.

To see this, change the code to show how often the loop is trying:

def findFirstPrime(n, k):
    count = 0
    while True:
        count += 1
        if count % 1000 == 0:
            print(f"at count {count:,}")
        if solovayStrassen(n,k):            
            return n, count
        n+=2    

Then add, e.g.,

random.seed(12)

at the start of the main program so you can get reproducible results.

Using floor (//) division, it runs fairly quickly, displaying

(6170518232878265099306454685234429219657996228748920426206889067017854517343512513954857500421232718472897893847571955479036221948870073830638539006377457, 906)

So it found a probable prime on the 906th try.

But with float (/) division, I never saw it succeed by blind luck:

at count 1,000
at count 2,000
at count 3,000
...
at count 1,000,000

Gave up then - "garbage in, garbage out".

One other thing to note, in passing: the + p in:

        j = (Jacobian(a, p) + p) % p                            

has no effect on the value of j. Right? p % p is 0.