import math,time,random
 
def miller(N,bases=(2,3,5,7,11,13)):
    if N<2:return False
    if N in bases:return True
    if N%2==0:return False
    d=N-1;s=0
    while d%2==0:d>>=1;s+=1
    for a in bases:
        if a>=N:continue
        x=pow(a,d,N)
        if x==1 or x==N-1:continue
        for _ in range(s-1):
            x=pow(x,2,N)
            if x==N-1:break
        else:return False
    return True
 
SP=[2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97]
 
print("=== pow(a,d,N) speed ===")
for bits in [64,128,256,512,1024]:
    N=(1<<bits)+643;it=max(1,min(500,50000//bits))
    t0=time.perf_counter()
    for _ in range(it):pow(2,N-1,N)
    t=(time.perf_counter()-t0)/it
    print(f"{bits:>5}b: {t*1e6:>8.0f} us")
 
print("\n=== Next prime ===")
for s in [10**6,10**9,10**12]:
    N=s+1
    if N%2==0:N+=1
    sk=0;mt=0;t0=time.perf_counter()
    while True:
        if any(N%p==0 for p in SP if p<N):sk+=1;N+=2;continue
        mt+=1
        if miller(N):break
        N+=2
    t=time.perf_counter()-t0
    print(f"dopo {s:.0e}: gap={N-s} {t*1000:.1f}ms sieve={sk} miller={mt}")
 
print("\n=== RSA prime gen ===")
random.seed(42)
for bits in [256,512]:
    times=[]
    for _ in range(3):
        t0=time.perf_counter()
        while True:
            N=random.getrandbits(bits)|(1<<(bits-1))|1
            if any(N%p==0 for p in SP):continue
            if miller(N):break
        times.append(time.perf_counter()-t0)
    print(f"{bits}b: {sum(times)/3*1000:.1f}ms avg")
 
print("\n=== Mersenne ===")
for p in [31,61,89,107,127]:
    M=(1<<p)-1;t0=time.perf_counter()
    r=miller(M);t=time.perf_counter()-t0
    print(f"M{p}: {'PRIMO' if r else 'COMP'} {t*1e6:.0f}us")
 
print("\nDone!")

aW1wb3J0IG1hdGgsdGltZSxyYW5kb20KCmRlZiBtaWxsZXIoTixiYXNlcz0oMiwzLDUsNywxMSwxMykpOgogICAgaWYgTjwyOnJldHVybiBGYWxzZQogICAgaWYgTiBpbiBiYXNlczpyZXR1cm4gVHJ1ZQogICAgaWYgTiUyPT0wOnJldHVybiBGYWxzZQogICAgZD1OLTE7cz0wCiAgICB3aGlsZSBkJTI9PTA6ZD4+PTE7cys9MQogICAgZm9yIGEgaW4gYmFzZXM6CiAgICAgICAgaWYgYT49Tjpjb250aW51ZQogICAgICAgIHg9cG93KGEsZCxOKQogICAgICAgIGlmIHg9PTEgb3IgeD09Ti0xOmNvbnRpbnVlCiAgICAgICAgZm9yIF8gaW4gcmFuZ2Uocy0xKToKICAgICAgICAgICAgeD1wb3coeCwyLE4pCiAgICAgICAgICAgIGlmIHg9PU4tMTpicmVhawogICAgICAgIGVsc2U6cmV0dXJuIEZhbHNlCiAgICByZXR1cm4gVHJ1ZQoKU1A9WzIsMyw1LDcsMTEsMTMsMTcsMTksMjMsMjksMzEsMzcsNDEsNDMsNDcsNTMsNTksNjEsNjcsNzEsNzMsNzksODMsODksOTddCgpwcmludCgiPT09IHBvdyhhLGQsTikgc3BlZWQgPT09IikKZm9yIGJpdHMgaW4gWzY0LDEyOCwyNTYsNTEyLDEwMjRdOgogICAgTj0oMTw8Yml0cykrNjQzO2l0PW1heCgxLG1pbig1MDAsNTAwMDAvL2JpdHMpKQogICAgdDA9dGltZS5wZXJmX2NvdW50ZXIoKQogICAgZm9yIF8gaW4gcmFuZ2UoaXQpOnBvdygyLE4tMSxOKQogICAgdD0odGltZS5wZXJmX2NvdW50ZXIoKS10MCkvaXQKICAgIHByaW50KGYie2JpdHM6PjV9Yjoge3QqMWU2Oj44LjBmfSB1cyIpCgpwcmludCgiXG49PT0gTmV4dCBwcmltZSA9PT0iKQpmb3IgcyBpbiBbMTAqKjYsMTAqKjksMTAqKjEyXToKICAgIE49cysxCiAgICBpZiBOJTI9PTA6Tis9MQogICAgc2s9MDttdD0wO3QwPXRpbWUucGVyZl9jb3VudGVyKCkKICAgIHdoaWxlIFRydWU6CiAgICAgICAgaWYgYW55KE4lcD09MCBmb3IgcCBpbiBTUCBpZiBwPE4pOnNrKz0xO04rPTI7Y29udGludWUKICAgICAgICBtdCs9MQogICAgICAgIGlmIG1pbGxlcihOKTpicmVhawogICAgICAgIE4rPTIKICAgIHQ9dGltZS5wZXJmX2NvdW50ZXIoKS10MAogICAgcHJpbnQoZiJkb3BvIHtzOi4wZX06IGdhcD17Ti1zfSB7dCoxMDAwOi4xZn1tcyBzaWV2ZT17c2t9IG1pbGxlcj17bXR9IikKCnByaW50KCJcbj09PSBSU0EgcHJpbWUgZ2VuID09PSIpCnJhbmRvbS5zZWVkKDQyKQpmb3IgYml0cyBpbiBbMjU2LDUxMl06CiAgICB0aW1lcz1bXQogICAgZm9yIF8gaW4gcmFuZ2UoMyk6CiAgICAgICAgdDA9dGltZS5wZXJmX2NvdW50ZXIoKQogICAgICAgIHdoaWxlIFRydWU6CiAgICAgICAgICAgIE49cmFuZG9tLmdldHJhbmRiaXRzKGJpdHMpfCgxPDwoYml0cy0xKSl8MQogICAgICAgICAgICBpZiBhbnkoTiVwPT0wIGZvciBwIGluIFNQKTpjb250aW51ZQogICAgICAgICAgICBpZiBtaWxsZXIoTik6YnJlYWsKICAgICAgICB0aW1lcy5hcHBlbmQodGltZS5wZXJmX2NvdW50ZXIoKS10MCkKICAgIHByaW50KGYie2JpdHN9Yjoge3N1bSh0aW1lcykvMyoxMDAwOi4xZn1tcyBhdmciKQoKcHJpbnQoIlxuPT09IE1lcnNlbm5lID09PSIpCmZvciBwIGluIFszMSw2MSw4OSwxMDcsMTI3XToKICAgIE09KDE8PHApLTE7dDA9dGltZS5wZXJmX2NvdW50ZXIoKQogICAgcj1taWxsZXIoTSk7dD10aW1lLnBlcmZfY291bnRlcigpLXQwCiAgICBwcmludChmIk17cH06IHsnUFJJTU8nIGlmIHIgZWxzZSAnQ09NUCd9IHt0KjFlNjouMGZ9dXMiKQoKcHJpbnQoIlxuRG9uZSEiKQ==