import random
from math import gcd

def expMod(a,k,n):
    ''' Calcul de a^n mod n'''
    res = 1
    while k > 0:
        if k % 2 == 1:
            res = (res * a) % n
        k = k // 2
        a = (a * a) % n
    return res

def jacobi(a,b,vb=0):
    """ Calcul de jacobi (a,b) pour 0 <= a < b avec b impair """
    res = 1
    if vb:
        print(res,a,b)
    while a > 0:
        while a % 2 == 0:
            if (b % 8) != 1 and  (b % 8) != 7:
              res = -res
            a //= 2
            if vb:
                print(res,a,b, " par simplification par 2")
        if a % 4 == 3 and  b % 4 == 3:
            res = -res
        r = b % a
        a , b = r , a
        if vb:
            print(res, a, b, " par loi de réciprocité")
    if b == 1:
        return res
    return 0

def solovayStrassen(n,k):
    for i in range(k):
        a = random.randint(1,n-1)
        d = gcd(a, n)
        if d > 1: # n 'est pas premier "
            return False
        j = jacobi(a,n)
        #r = expMod(a,(n-1)//2,n)
        r = pow(a,(n-1)//2,n)
        if (r - j) % n  != 0:
            return False
    #print(" Est premier avec une probabilité d'erreur < ", 2.0 ** -k)
    return True

def next_pseudo_prime(x, k=150):
    """ x est impair. Renvoie plus petit pseudo prime >= x """
    if x % 2 == 0:
        x += 1
    while not solovayStrassen(x,k):
        x += 2
    return x