18.783 Lecture 11 Pollard p-1.sagews

def pollard_pm1(N,B=0):
    if not B: B=ceil(sqrt(N))
    a = Integers(N).random_element()
    b = a
    for ell in primes(B):
        q = 1
        while q < N:
            q *= ell
        b = b^q
        d = gcd(b.lift()-1,N)
        if d == N: return 0
        if d > 1: return d
    return 0
    
def random_unsafe_prime(bits):
    while true:
        a=randint(0,bits)
        b=randint(0,floor((bits-a)/log(3,2)))
        c=randint(0,floor((bits-a-floor(b*log(3,2)))/log(5,2)))
        d=floor((bits-a-floor(b*log(3,2))-floor(c*log(5,2)))/log(7,2))
        p = 2^a*3^b*5^c*7^d+1
        if is_prime(p): return p

factor(random_unsafe_prime(512)-1)

2^43 * 3^10 * 5^130 * 7^54

%time
p1=random_unsafe_prime(512)
p2=random_prime(2^512,2^511)
print p1,p2
print pollard_pm1(p1*p2)

4807479287630479943604466481712153408482548810541337478124481833945698978620075887164841165409365474230693820407553062417241073987983769600000000000000001 13006152475104545524103623453283381207648191608217087131830391172329570734616989058507218055127554465392323393459642419857444986550994965844875037794996683 4807479287630479943604466481712153408482548810541337478124481833945698978620075887164841165409365474230693820407553062417241073987983769600000000000000001 CPU time: 1.08 s, Wall time: 1.13 s