partial_d.sage

def partial_p(p0, kbits, n):
    PR.<x> = PolynomialRing(Zmod(n))
    nbits = n.nbits()

    f = 2^kbits*x + p0
    f = f.monic()
    roots = f.small_roots(X=2^(nbits//2-kbits), beta=0.3)  # find root < 2^(nbits//2-kbits) with factor >= n^0.3
    if roots:
        x0 = roots[0]
        p = gcd(2^kbits*x0 + p0, n)
        return ZZ(p)

def find_p(d0, kbits, e, n):
    P = var('P')
    for k in range(1, e+1):
        print(k)
        results = solve_mod([e*d0*P - k*(n-P+1)*P + k*n == P], 2^kbits)
        print(results)
        for x in results:
            p0 = ZZ(x[0])
            p = partial_p(p0, kbits, n)
            if p:
                return p


if __name__ == '__main__':
    print "start!"
    n = 22886501935336065182326527443310749425171647238626282413401383196346489429704102057721649923432485366546493527638108629871871958769468030807810941317738572491741321770829422490153900877578213060028200750443514516859983684513695362081632707707911149834332602372979434088912942912851938670551379551225494633726058082769820094808639196295355917935311972995448897110349031393317379557438029524114635337222762331676438910234927812848930884968824910365642189803231791506971024009754951423293404469711697860605662568095139192065818762238142932255427846376827236900905886912090980882932064525237833246818119279930638745034749
    e = 65537
    nbits = n.nbits()
    kbits = 1022
    d0 = 34853327299025755595075304279734536658859779367199829099864542202921236791274445811280627334793831109517046902319361193799141879670245490634859470527482992704459154837275611693065523725040305710310173322460987652380350573853183919952409167226508053603975272068685866320882516555872353263491319823526487495001
    print "lower %d bits (of %d bits) is given" % (kbits, nbits)

    p = find_p(d0, kbits, e, n)
    print "found p: %d" % p