1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:

def rab_test(N, rounds = 40):
    if N in [2, 3]:
        return True
    if N % 2 == 0:
        return False
    n = N // 2
    t = 1
    while n % 2 == 0:
        n = n // 2
        t += 1
    a = random.randrange(2, N-1)
    u = gmpy2.powmod(a, n, N)
    if u == 1 or N-u == 1:
        return True
    for k in range(1, t):
        u = gmpy2.powmod(u, 2, N)
        if N-u == 1:
            return True
    return False