[PlaidCTF 2019] Project Eulernt

#!/usr/bin/env python3
import gmpy2

from decimal import Decimal

N = int(gmpy2.fac(333))
#N = 10334465434588059156093965538297516550622260041682062823432902469783188597914276568552700194849877929894375950252570477080418352732597658745665925604704669227133726477243854317836635130694123893711638533001980496229875665476598568821806170303765540489814402234159901540440432134155844542962445153646330595588291605924429211352279943471372817279938720974895260387784578239150931816946786416232516666251965421919651838044618050991294403546958930745419743836966520198735201123255884089263272829846640538826979843642885775791641575109178753509580001660392092396798648924375401024147883702298145910046889402880394195369984000000000000000000000000000000000000000000000000000000000000000000000000000000000
sN = int(gmpy2.isqrt(N))
#sN = 3214726338988757399964463840205273148284735043324463579894976679845078166928105412104944973948893914339037572694382785661727648297539107767478128297633669341356440278480314502443731079340424764653103468238563073341496690901434197268615240607985890327844073738551115260849983966971570699838147501655616953786428037017304945538845583678438817092853062

k = int(input("Enter number: "))

goodness = Decimal(abs(k - sN)) / sN 

if k and N % k == 0 and goodness < 1e-8:
    print(open('/home/eulernt/flag.txt').read())
elif k and N % k == 0 and goodness < 1e-4:
    print("Good work! You're getting there.")
else:
    print("Nope!")

You must find a divisor for $N$ closer to $N^{0.5}$($N = 333!$) with relative error 1e-8. We can easily factorize $N$. After that, I just apporached with simple climbing method such that

1. Randomly distribute the factors of $N$ to A or B
2. Select 3 factors of A and B, and swap it if diff is decreased.
3. If there is no such factors, terminate

It is not that elaborate apporach and requires a lot of time, but it works XD

from decimal import Decimal
import random
plist = [4,9,25,49,169,121,361,3,5,7,19,29,37,43,47,59,61,89,97,101,103,107,109,167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331,101*101,103*103,107*107,109*109,163*163, 157*157, 151*151, 149*149, 139*139, 137*137, 89*89, 97*97, 113*113,127*127,131*131,83*83, 79*79, 73*73]
sN = 3214726338988757399964463840205273148284735043324463579894976679845078166928105412104944973948893914339037572694382785661727648297539107767478128297633669341356440278480314502443731079340424764653103468238563073341496690901434197268615240607985890327844073738551115260849983966971570699838147501655616953786428037017304945538845583678438817092853062
N = 10334465434588059156093965538297516550622260041682062823432902469783188597914276568552700194849877929894375950252570477080418352732597658745665925604704669227133726477243854317836635130694123893711638533001980496229875665476598568821806170303765540489814402234159901540440432134155844542962445153646330595588291605924429211352279943471372817279938720974895260387784578239150931816946786416232516666251965421919651838044618050991294403546958930745419743836966520198735201123255884089263272829846640538826979843642885775791641575109178753509580001660392092396798648924375401024147883702298145910046889402880394195369984000000000000000000000000000000000000000000000000000000000000000000000000000000000
const = 2**163 * 3**81 * 5**39 * 7**25 * 11**15 * 13**12 * 17**10 * 19**7 * 23**7 * 29**5 * 31**5 * 37**4 * 41**4 * 43**3 * 47**3 * 53**3 * 59**2 * 61**2 * 67**2 * 71**2 * 73**1 * 79**1 * 83**1
#const = 23877790403724415895889742090969477905892583014070051248745074024151351725346141205015301537285705927194482434815550834956458585791415256383920842183350368172667073369935828137022480093631167679142456715988579872058431309674840741418723140111761796595358755717120000000000000000000000000000000000000000

def solve():
  for _ in range(10000):
    i = random.randint(0,len(plist)-1)
    j = random.randint(0,len(plist)-1)
    plist[i],plist[j]=plist[j],plist[i]
  A = []
  B = []
  totA = 1
  totB = 1
  for p in plist:
    if random.randint(0,1) == 0:
      A.append(p)
      totA *= p
    else:
      B.append(p)
      totB *= p

  while Decimal(abs(totA*const-sN)) / sN >= 1e-8:
    diff = abs(totA-totB)
    find = False
    for i1 in range(len(A)):
      if find: break
      for i2 in range(i1+1,len(A)):
        if find: break
        for i3 in range(i2+1,len(A)):
          if find: break
        for j1 in range(len(B)):
          if find: break
          for j2 in range(j1+1,len(B)):
            if find: break
            for j3 in range(j2+1,len(B)):
              if find: break
              totA_new = totA//A[i1]//A[i2]//A[i3]*B[j1]*B[j2]*B[j3]
              totB_new = totB//B[j1]//B[j2]//B[j3]*A[i1]*A[i2]*A[i3]
              diff2 = abs(totB_new-totA_new)
              if diff > diff2:
                totA = totA_new
                totB = totB_new
                A[i1],B[j1] = B[j1],A[i1]
                A[i2],B[j2] = B[j2],A[i2]
                A[i3],B[j3] = B[j3],A[i3]
                if Decimal(abs(totA*const-sN)) / sN < 1e-5:
                  print("---------------------------")
                  print(totA*const)
                  A.sort()               
                  print(A)
                  print("---------------------------")
                  print(totB*const)
                  B.sort()
                  print(B)
                  print("---------------------------")
                  print(Decimal(abs(totA*const-sN)) / sN)
                find = True
    if not find:
      print("nono...")
      return

  val = totA*const
  print("********DONE*********")
  print(val)

  print(Decimal(abs(val-sN))/sN)
  exit()

while True:
  solve()