tux_sol_dicectf.sage

import itertools
from sympy.solvers.diophantine.diophantine import diop_linear
from sympy import symbols

def getPolyInfo(poly):
    HM = poly.monomials()[0]            # HM: head monomial
    HC = poly.monomial_coefficient(HM)  # HC: head coefficient
    HT = HC*HM                          # HT: head term
    HI = poly.exponents()[0]            # HI: head index
    return {'HT':HT,'HC':HC,'HM':HM,'HI':HI}

# https://eprint.iacr.org/2012/675.pdf
def multivariate(N,E,beta,m,improve=False):

    # Setup multivariate simultaneous equations
    l = len(E)
    P,x = PolynomialRing(ZZ,l+1,'x').objgens() # y = x[-1]
    F = []
    for k in range(l):
        F.append((-1+x[k]*(x[-1]+N)))
    X = [floor(N^beta) for _ in range(l)] + [floor(3*N^.5)] # upper bounds

    # Construct polynomials: G
    G = []
    I = []
    for k in range(l):
        g = {}
        i_k = []
        for i in range(m+1):
            for j in range(i+1):
                g[i,j] = x[k]^(i-j) * F[k]^j * E[k]^(m-j)

                # check HM,HI
                polyinfo = getPolyInfo(g[i,j])
                assert polyinfo['HM'] == x[k]^i * x[-1]^j
                assert tuple(polyinfo['HI']) == tuple([i if (_ == k) else j if(_ == l) else 0 for _ in range(l+1)])
                i_k.append(tuple(polyinfo['HI']))
        G.append(g)
        I.append(i_k)

    # Get index set I+
    I_plus = set()
    for i in itertools.product(*I):
        I_plus.add(tuple([sum(i[k][j] for k in range(len(i))) for j in range(l+1)]))
    I_plus = sorted(I_plus)

    if improve:
        # Get index set I++
        I_pp = set()
        for idx in I_plus:
            if idx[-1] <= max(idx[:l]) or .5*idx[-1] + sum(idx[:l])*beta - sum(min(idx[k],idx[-1]) for k in range(l)) < 0:
                I_pp.add(idx)
        I_pp = sorted(I_pp)
        I_plus = I_pp

    # Construct polynomials: G+
    G_plus = []
    for idx in I_plus:
        J = []
        for j_tuples in itertools.product(range(idx[-1]+1),repeat=l):
            if sum(j_tuples) == idx[-1] and not any([idx[k] < j_tuples[k] for k in range(l)]):
                J.append(j_tuples)

        # Find GCD(HC)
        gcdn = -1
        for j in J:
            polyinfo = getPolyInfo(prod(G[k][idx[k],j[k]] for k in range(l)))
            gcdn = polyinfo['HC'] if (gcdn == -1) else gcd(gcdn,polyinfo['HC'])

        if len(J) > 1:
            arr = [getPolyInfo(prod(G[k][idx[k],J[j][k]] for k in range(l)))['HC'] for j in range(len(J))]
            sols = diop_linear(sum(arr[t]*x0 for t,x0 in enumerate(symbols("x0:%d"%len(J),integer=True)))-gcdn)
            a = [int(expr.subs({t:0 for t in sols[-1].free_symbols})) for expr in sols]
            assert sum(arr[t]*a[t] for t in range(len(J))) == gcdn
            poly = sum(a[j] * prod(G[k][idx[k],J[j][k]] for k in range(l)) for j in range(len(J)))
        else:
            poly = sum(prod(G[k][idx[k],j[k]] for k in range(l)) for j in J)

        # check HM
        for j in J:
            polyinfo = getPolyInfo(prod(G[k][idx[k],j[k]] for k in range(l)))
            assert polyinfo['HM'] == prod(x[k]^idx[k] for k in range(l+1))

        g_poly = poly.subs({x[i]:x[i]*X[i] for i in range(l+1)})
        G_plus.append(g_poly)

    # Get monomials
    monomials = set()
    for poly in G_plus:
        for m1 in poly.monomials():
            monomials.add(m1)
    monomials = sorted(monomials) # sort monomials

    # Setup matrix for LLL
    assert len(G_plus) == len(monomials)
    print("Matrix dimensions: {}x{}".format(len(G_plus),len(monomials)))
    L = Matrix(ZZ,len(G_plus),len(monomials))
    for i,poly in enumerate(G_plus):
        for j,m1 in enumerate(monomials):
            L[i,j] = poly.monomial_coefficient(m1)

    L = L.LLL()

    # Get system of polynomials
    H = [0 for _ in range(L.nrows())]
    for i in range(L.nrows()):
        for j,m1 in enumerate(monomials):
            H[i] += P(m1 * L[i,j] / m1.subs({x[k]:X[k] for k in range(l+1)}))

    # Find roots using Groebner basis
    H = [eq.change_ring(QQ) for eq in H]
    I = Ideal(*H[0:l+2])
    g = I.groebner_basis('giac')
    if g == [1]:
        return (0,0)
    x0 = var('x',n=l+1)
    sols = solve([eq(x0) for eq in g],x0,solution_dict=True)
    for s in sols:
        if [s[i].is_integer() for i in x0] == [True]*(l+1):
            p_plus_q = 1-s[x0[-1]]
            x1 = PolynomialRing(ZZ,'x').gen()
            f = x1**2 - p_plus_q*x1 + N
            roots = f.roots()
            p,q = roots[0][0],roots[1][0]
            return (p,q)
    return (0,0)

N = 14337596625981852268356733199577363454740211423306675087406912203612106754740920535831186446286722686623181504391319837160175722628785779086065934404508907495375575604993274597567269696998989342554512169936041805398614285019559611603192791406141673196062487187176974406170426326985080300401992871878272750861304966084697316255932229517848590059063403143320068885919541934394285204068719613486701354616465507009243987608727662705537132041395688365562715140968368814327462123434153164039701357637150078015089287814859891259359522535418157646847992980552057214329745061692026530855072281235664669710276172313808317219377
E = [12976676506376070110040406570149996101394527411415881605860843983349363069127746407886818103004002117366801754683449911816105354951654828869006253634088640931647307722477393307231188093852730642851927610384508472618758673499740002362888463449907850898256709493427721244039346744775543365490896192997672156351539772375297375865701816236826858372370828737647663910668758559951734259138880399000934641709724370982987684871113999721232035923836779554293705654708805788634929466224221204317582909006400475748926200070330888314906329563110519417492641071166678480550894995905707586159556311266401315984011972916692724662007, 1863763332200533355187495032998617435554266559806936864499583761336209979224578994667995146126874905256288657356393344293005098598793270446199080399410824354463753330703855276639236442323809498935079594864731622806981650099779625802394492866772783164556122039084479891630004937148495825596230840040565205548593767270423848984375003632523408070963142589587241988428801951874519861970126281101472254516564569333579422627371073253518765484477877683153600580032505599620122184765699186189626337136723831818155748441783703710300026814764444157324740641952652611754903513683898296023658398333790727757433569553518827781081, 2266107108355028257265149613659513002904527665415922083354756255551387373798380715368636460529830761909106090742940804896135018713646665813103359939678904950388594536893820576541648929534374197197699236090310478668004566290689920782095321065479404918697144906086349325096417620451929567373057958837794405051832027393710209669487434964110954150175622748681581358286557194922971711892121496052503670245966738246893335673371911853174980221559458430193451035071192245185753464656500832472716390638460484671724823329824542466725661315185971990283547100285181536347175215757494437839674240264046716971285767053511953014531, 3012665180743215926044185696769375844027893756273454464899436474243670235496407935713055380725328632598984826511997157495169513834153639995899362029217480834870286205900296318828132585598137285805346657507373383217564586845718076975940327895945639149538433468219499891765910172868182009270172797753293191815210844128651385425158004075468732463275400326825588500413826705324683883847629057842100412354271506810539904381290499384122685388561018876503682490390411350939677572500607680288271353610724980926630513275655559323486108263101876350397264023226898671971480444539423455934221932088753929657082051926857943772249, 11654698355361642253478315394509778991884549939645189298025844618346681492919905598771246194066168408610846920719586231628616008528010499420499104967297209036853716237681916872940149475834811605883772954478356917401690932414133714754145102545282932055458620466374791600444965981474694964466374906501557735419193196210081985354200587280684723865152043847523132554069361865253297610018545961901224789197822102327845824463263202308386229169765022168048658546459296356873527223601643713792899944747399360948976054350376734975776305888045491739605171010206699460442078038128020076159311221789613452312845058576145149846029]
c = 3956820731464413675248976184330090386023400385276540416612412624089512558297526443036270899672623474843471114663110309683900437072074275592048007465581536399064139915757208597226420873829762312175028263779702542540471556091196109301321307501425411404463271934422940357892918540470499135323435236378549758184835624591770086862646876844851191370941059743180644939927085439332975325190261804739679687807083505313731296227816013191899297272702691347829966232657275540293776915188791505265719622572712788357792237668230216117983278509746622406227163590935715374994406994094600128249001135177407180401006889949317057904629
p,q = multivariate(N,E,.5,1,True)
print(p,q)
assert p*q == N
phi = (p-1)*(q-1)
d = pow(0x1337,-1,phi)
from Crypto.Util.number import long_to_bytes
print(long_to_bytes(pow(c,d,N)))