import libnum
def getFullP(low_p, n):
    R.<x> = PolynomialRing(Zmod(n), implementation='NTL')
    p = x*2^512 + low_p
    root = (p-n).monic().small_roots(X = 2^128, beta = 0.4)
    if root:
        return p(root[0])
    return None
    
def phase4(low_d, n, c):
    maybe_p = []
    for k in range(1, 4):
        p = var('p')
        p0 = solve_mod([3*p*low_d  == p + k*(n*p - p^2 - n + p)], 2^512)
        maybe_p += [int(x[0]) for x in p0]
    #print(maybe_p)
    
    for x in maybe_p:
        P = getFullP(x, n)
        if P: break
    
    P = int(P)
    Q = n // P
   
    assert P*Q == n
    print(P)
    print(Q)
    
    d = inverse_mod(3, (P-1)*(Q-1))
    print(d)
    print(libnum.n2s(int(power_mod(c, d, n))))
    


n = 92896523979616431783569762645945918751162321185159790302085768095763248357146198882641160678623069857011832929179987623492267852304178894461486295864091871341339490870689110279720283415976342208476126414933914026436666789270209690168581379143120688241413470569887426810705898518783625903350928784794371176183
c = 56164378185049402404287763972280630295410174183649054805947329504892979921131852321281317326306506444145699012788547718091371389698969718830761120076359634262880912417797038049510647237337251037070369278596191506725812511682495575589039521646062521091457438869068866365907962691742604895495670783101319608530
low_d = 787673996295376297668171075170955852109814939442242049800811601753001897317556022653997651874897208487913321031340711138331360350633965420642045383644955

phase4(low_d, n, c)