算法背景与原理：

1、Fermat 小定理：给定素数 p，a∈Z，则有 a^(p-1)%p=1
2、Fermat 素性检测算法：奇整数 m，若任取一整数 2<=a<=m-2，gcd（a，m）=1，使得 a^(m-1)(mod m)=1，则 m 至少有 1/2 的概率为素数

算法步骤：

1、从键盘输入待检测的大整数 m
2、给出安全参数 k
3、随机选取整数 a，满足 a∈[2,m-2]
4、计算 g=gcd（a，m），如果 g=1 进行下一步，否则不是素数，跳出
5、计算 r=a^(m-1)(mod m)，如果 r=1 则 m 可能是素数，否则不是素数，跳出
6、重复上述过程 k 次，如果每次判定 m 都可能是素数，那么 m 是素数的概率是 1-1/2^k

import random
import math

'''
def modpow(a,m):
    s=m-1
    r=1
    while s!=0:
        if s%2==1:
            r=r*a%m
        a=(a*a)%m
        s=s//2
    return r
'''
m=int(input('m = '))
k=int(input('安全系数k = '))

s=k
while s>0:
    a=random.randint(2, (m-2))
    print('a=',a,end="")
    #print(',(',a,',',m,')=',math.gcd(a, m),end="")
    if math.gcd(a, m)==1:
        #print(a,'^(',m,'-1)(mod',m,')≡',modpow(a, m),',',end="")
        #or  print(a,'^(',m,'-1)(mod',m,')≡',pow(a, m-1, m),',',end="")
        if pow(a, m-1, m)==1:
        #if modpow(a, m)==1:
            s=s-1
            print(',故m可能为素数  ',100*(1.0-(1.0/(2**(k-s)))),'%')
    else:
        print('m为合数')
        break
if s==0:
    print('m',100*(1.0-(1.0/(2**k))),'%可能为素数')

仅供参考，其实这玩意网上一搜一大把
主要解决模指数问题就行，如果仅用**做大数幂运算的话会很慢