#python2 -m primefac -m=p+1 13600689950702678449784835151815974930081348852322430023033556990750407302536196017142902387097892642345597564002803049369859021622154057824025261445149722108978168478137715149349659619454621524807710623319033897949123127670966633488318687764955391182209483781455502811562493433978547603095334295954771596786145240556520883744576353821136738684148011064935790222495233081914179923414412186454172809843508100994530119340028548341289369372240024738698118638987563907

#p =  26095988705098082141146804081349098362896803568710141131962347542899307107886799504316790900281317552305916634959743976639058747315050655853768513454406251

#q= n/p = 521179331597643717826719455288904380283486818776725416395655915960819344127365068024226118263530149099094250120477785823751838585076764808018002955416767480620795032894457484184359000800607418687273707348464760198903221109476238119882003449193031638494925360314500101182222492365386087041855756138309223907657


import gmpy2

N = 521179331597643717826719455288904380283486818776725416395655915960819344127365068024226118263530149099094250120477785823751838585076764808018002955416767480620795032894457484184359000800607418687273707348464760198903221109476238119882003449193031638494925360314500101182222492365386087041855756138309223907657

def fermat_factor(n):
    assert n % 2 != 0

    a = gmpy2.isqrt(n)
    b2 = gmpy2.square(a) - n

    while not gmpy2.is_square(b2):
        a += 1
        b2 = gmpy2.square(a) - n

    p = a + gmpy2.isqrt(b2)
    q = a - gmpy2.isqrt(b2)

    return int(p), int(q)

if __name__ == "__main__":
    (p, q) = fermat_factor(N)

    print("p = {}".format(p))
    print("q = {}".format(q))