RWCTF 3rd - Crypto - Old Curse solve script

sol.sage

from Crypto.Util.number import *
PR.<qr>=PolynomialRing(ZZ)

def calc_params(e,N):
    delta = 0.1
    gama = 0.05
    beta = log(e,N).n()
    alpha = 0.25
    return alpha,beta,delta,gama

def condition_check(alpha,beta,delta,gama):
    print delta
    print 7/6 + (1/3)*alpha-gama-(1/3)*sqrt((2*alpha+1)*(2*alpha+6*beta-6*gama+1))
    return delta < 7/6 + (1/3)*alpha-gama-(1/3)*sqrt((2*alpha+1)*(2*alpha+6*beta-6*gama+1))

def calc_tau(alpha,beta,delta,gama):
    return (1-2*delta-2*alpha-2*gama)/(2*(1+2*alpha))

def calc_XYZ(alpha,beta,delta,gama,N):
    X = next_prime(floor(4*N^(beta+delta-1)))
    Y = next_prime(floor(3*sqrt(2)*N^(1/2+alpha)))
    Z = next_prime(floor(N^gama))
    return X,Y,Z

def multivariate(pol,X,Y,Z,tau,mm,alpha,beta,delta,gama,e,N):
    P.<x,y,z> = PolynomialRing(ZZ)
    tt = floor(tau*mm)+1
    #tt = 1
    ff = x*y - (N)*x + z
    polys = []
    monomials = []
    for k in range(mm+1):
        for i1 in range(k,mm+1):
            i2 = k
            i3 = mm-i1
            poly = x^(i1-k)*(z)^i3*ff^k*e^(mm-k)
            polys.append(poly)
    for k in range(mm+1):
        i1 = k
        for i2 in range(k+1,i1+tt+1):
            i3 = mm-i1
            poly = y^(i2-k)*(z)^(i3)*ff^k*e^(mm-k)
            polys.append(poly)
    
    polys = sorted(polys)
    for poly in polys:
        monomials += poly.monomials()
    monomials = sorted(set(monomials))

    nn = len(monomials)
    BB = Matrix(ZZ,nn)
    print len(polys)
    print len(monomials)
    for ii in range(nn):
        #BB[ii, 0] = polys[ii](0, 0, 0)
        for jj in range(nn):
            if monomials[jj] in polys[ii].monomials():
                #BB[ii, jj] = polys[ii].monomial_coefficient(monomials[jj]) * monomials[jj](X,Y,Z)
                BB[ii,jj] = polys[ii](X*x,Y*y,Z*z).monomial_coefficient(monomials[jj])
    #print BB
    print monomials
     
    BB = BB.LLL(delta=0.75)
    print "LLL done"
    #print BB
    
    dims = nn
    H = [(i, P(0)) for i in xrange(dims)]
    H = dict(H)
    for j in xrange(dims):
        for i in xrange(dims):
            H[i] += P((monomials[j] * BB[i, j]) / monomials[j](X,Y,Z))

    poly = P(gcd(H[0],H[1]))
    print poly
    print -poly.monomial_coefficient(x)
    print poly.monomial_coefficient(x*y)
    return N - ceil((-poly.monomial_coefficient(x))/poly.monomial_coefficient(x*y))
    
    
def gen(p,q,Ep,Eq,delta,gama):
    N = p*q
    u = random_prime(floor(N^delta), proof=False, lbound=ceil(N^(delta-0.02)))
    print u
    print log(u,N).n()
    w = random_prime(floor(N^gama), proof=False, lbound=ceil(N^(gama-0.02)))
    print w
    print log(w,N).n()
    e = (w*inverse_mod(u,Ep*Eq))%(Ep*Eq)
    print e
    print log(e,N).n()
    v = (e*u - w)/(Ep*Eq)
    print v
    print log(v,N).n()
    return e,u,w


def main(e,N,mm):
    P.<x,y,z> = PolynomialRing(ZZ)
    alpha,beta,delta,gama = calc_params(e,N)
    assert condition_check(alpha,beta,delta,gama)
    pol = x*y- (N)*x + z
    X,Y,Z = calc_XYZ(alpha,beta,delta,gama,N)
    tau = calc_tau(alpha,beta,delta,gama)
    res = multivariate(pol,X,Y,Z,tau,mm,alpha,beta,delta,gama,e,N)
    return res

if __name__ == "__main__":
    ''' 
    N = 2445641520497188372893910329538675824331454921520163900426562393634418526897575682249916293416221269674459540700624274860236238684609738360751815410091617
    e = 207753540686843587408555602893982678168821441852165899252123932416370824148707563812033872059010473801740084336709522813588017197501164099322578137710783
    e = 1385541386350864526341998465216697876045250157636252495557496884777332366331828079007542764765138085082758016229583509718779067112153670966735233957418107
    p = 68592042559830614325177858341490526647123827945850285759827931818992553395171
    q = 35654886911465489386559478739125595731698572982484092580287175860557076482027
    assert N == p*q
    a = 0
    b = 9
    mm = 2
    delta = 0.085
    gama = 0.045
    Ep = EllipticCurve(GF(p),[a,b])
    Eq = EllipticCurve(GF(q),[a,b])
    p_order = Ep.order()
    q_order = Eq.order()
    main(e,N,p,q,p_order,q_order,delta,gama,mm)
    '''
    '''
    N = 43115265506687226436196728756959707266402158394261294759458139340520129183826747
    e = 4429109683378321635373164354359544019395496659337936831132897706681971178351139
    delta = 0.15
    gama = 0.05
    mm = 2
    p = 6672224014662340178579721474326734152749
    q = 6461903169309154483833797011785886506503
    p_order = 6672224014662340178579721474326728185600
    q_order = (43115265506687226436196728756959661642298658940914579533819094510831187730169600/p_order)
    main(e,N,p,q,p_order,q_order,delta,gama,mm)
    '''
    
    a = 0x4a2e6dec145a807c30098f59a73305f5e0cadb3d6a90da183ded22d4dbee0865bbdbed56698289afb60a534a0af0226bbb94245671f747e820f6a70c928dfa07
    b = 0x120a461de47b36eeacd572a40bea21b0ed1475503905d1c2a1bab2c812bbc8f548e638f126f84fcddc2188c9f2bd3ee282bc5cfe253424d4594ee72985b3fd8c
    
    N = 80330528881183983072964816732300543404856810562533626369319300810697262966387144944887576330528743612839739692299784591097332512948890518183519167192046959230085412831864255497489112175176914874596237618253755256608956517757030073479666104923402013469283716999320744856718736837534911809839541660207743594867
    mm = 2
    
    e = 78452652317506438607956636739779994986676384637399723342738736371812868831141251164966879331214017314432739387076791674001159059604426825547538902010774841189596518785149221523738464397224366361779781148300651051284198636694801404816891957209985325619623109930150535820404950711233032177848101830061155574970
    #e,u,w = gen(p,q,p_order,q_order,delta,gama)
    res = main(e,N,mm)
    L = ecm.factor(res)
    print L
    q_order = 131*1108897087*505386797752007*120659691081137900860528439558149439256036479214584879088476613192185895986414329679519081477454257879221194033908435726005914629
    p_order = res//q_order
    lbound = floor(sqrt(N)-N^(1/4))
    ubound = ceil(sqrt(2*N)+N^(1/4))
    assert lbound< p_order and p_order < ubound

    phi = lcm(p_order,q_order)
    E = EllipticCurve(Zmod(N),[a,b])
    Stone = E(5316297494616251967087180573684467112077977207314228196651011473838683480275875989908990738740861375687186766156200219641981169308660139151062711296717379891376294785675104640775506724244803337279235747630215478504380272738204733311972022712834357078381541224632797503360732934454187646031643331529389570159,73177062713968648963738410812785853174528721431172461113561340178691492280271903912043554814810920745154304747328073913103230849027745226637330284520633847773874342467137552022725301429074046921710660867115557994943332628756632246059800601063580017261698262663178072317324978782579376388601713100806653808812)
    d = inverse_mod(e,phi)
    Heart = d*Stone
    print long_to_bytes(Integer(Heart[0]))
    R = Zmod(N)
    PR.<xn> = PolynomialRing(R)
    bits = 250
    for i in range(8):
        print i
        f = p_order - i*2^250 - xn
        f = f.monic()
        roots = f.small_roots(X=2^bits, beta=0.4,epsilon = 1/150)
        if len(roots) != 0:
            r = Integer(roots[0]+i*2^250)
            print "[+]found root : ",r
            break
    p = p_order - Integer(r)
    q = N/p
    assert N%p == 0
    flag = long_to_bytes(Integer(Heart[0])^^(p*2^248))
    print flag