5 Mar 2019

RSA从入门到入土

RSA从入门到入土

  • 直接分解n

    p-q 很大或很小

    p-1光滑 Pollard’s p − 1 算法分解N

    from Crypto.Util.number import *
      
    def Pollard_p_1(N):
        a = 2
        while True:
            f = a
            # precompute
            for n in range(1, 80000):
                f = pow(f, n, N)
            for n in range(80000, 104729+1):
                f = pow(f, n, N)
                if n % 15 == 0:
                    d = GCD(f-1, N)
                    if 1 < d < N:
                        return d
            print(a)
            a += 1
    

    p+1光滑 Williams’s p + 1 算法分解N

    yafu

  • 小公钥指数攻击

    • 条件

      e特别小,如e = 3

    • 原理

      $c≡m^{e}\ (mod\ n)$

      $ m^{e} = c + kn $

      $ m = (c + kn)^{1/e} $

      枚举k,开e次根,直到开出整数

      def small_exponent(pubkey, cipher):
      	N = 
          e = 
      	c = 
      	i = 0
      	i = 118719488
      	while 1:
      		if(gmpy.root(c+i*N, 3)[1]==1):
      			plaintext = gmpy.root(c+i*N, 3)[0]
      			break
      		i += 1
      	plaintext = transform.int2bytes(plaintext)
      	print(plaintext)
      
  • 共模攻击

    • 条件

      用相同的n不同的e对明文m加密

    • 原理

      $c_1 ≡ m^{e_1} (mod\ n)$

      $c_2 ≡ m^{e_2} (mod\ n)$

      $gcd(e_1,e_2) = 1$

      存在$s_1,s_2$使$s_1e_1+s_2e_2 = 1$

      $c_1≡m^{e_1}\ (mod\ n)$

      $c_2≡m^{e_2}\ (mod\ n)$

      所以

      $c_1^{s_1}c_2^{s_2}≡(m^{e_1})^{s_1}(m^{e_2})^{s_2}\ (mod\ n)$

      $c_1^{s_1}c_2^{s_2}≡m^{e_1s_1+e_2s_2}\ (mod\ n)$

      $c_1^{s_1}c_2^{s_2}≡m\ (mod\ n)$

      def common_modulus_attack(cipher1,cipher2):
      	c1 = 
      	c2 = 
      	e1 = 
      	e2 = 
      	N = 
      	s = gmpy2.gcdext(e1,e2)
      	s1 = s[1]
      	s2 = -s[2]
      	c2r = gmpy2.invert(c2, N)
      	plaintext = (pow(c1,s1,N) * pow(c2r,s2,N)) % N
      	plaintext = transform.int2bytes(plaintext)
      	print(plaintext)
      
  • dp

    import gmpy2
    e = 
    n = 
    dp = 
    for x in range(1, e):
        if(e*dp%x==1):
            p=(e*dp-1)//x+1
            if(n%p!=0):
                continue
            q=n//p
            phin=(p-1)*(q-1)
            d=gmpy.invert(e, phin)
    
  • dp、dq

    import gmpy2
    p = 
    q = 
    dp = 
    dq = 
    n = p*q
    phin = (p-1)*(q-1)
    dd = gmpy2.gcd(p-1, q-1)
    d=(dp-dq)//dd * gmpy2.invert((q-1)//dd, (p-1)//dd) * (q-1) +dq
    
  • Rabin算法

    • 条件

      e = 2

    • 原理

      $c=m^{2}\ mod\ n$

      $m_p =\sqrt{c}\ mod\ p $

      $m_q =\sqrt{c}\ mod\ q $

      若$p≡q≡3\ (mod\ 4)$

      $m_p = c^{1/4(p+1)}\ mod\ p$

      $m_q = c^{1/4(q+1)}\ mod\ q$

      $y_pp+y_qq=1$

      $gcdext(p,q)$解出$y_p,y_q$

      解出四个明文

      $a = (y_ppm_q+y_qqm_p)\ mod\ n$

      $ b = n -a$

      $c = (y_ppm_q-y_qq_mp)\ mod\ n$

      $d = n - c$

      Jarvis-OJ-hard-rsa

       # coding=utf-8
       import gmpy2
       import string
       from rsa import transform
       n = 87924348264132406875276140514499937145050893665602592992418171647042491658461L
       e = 2
       p = 275127860351348928173285174381581152299
       q = 319576316814478949870590164193048041239
       c = open("flag.enc" ,'rb').read()
       c = transform.bytes2int(c)
       # 计算yp和yq
       inv_p = gmpy2.invert(p, q)
       inv_q = gmpy2.invert(q, p)
             
       # 计算mp和mq
       mp = pow(c, (p + 1) / 4, p)
       mq = pow(c, (q + 1) / 4, q)
             
       # 计算a,b,c,d
       a = (inv_p * p * mq + inv_q * q * mp) % n
       b = n - int(a)
       c = (inv_p * p * mq - inv_q * q * mp) % n
       d = n - int(c)
             
       for i in (a, b, c, d):
           print(transform.int2bytes(i))
      
  • Wiener’s attack

  • Factoring with High Bits Known

    • 条件

      已知N的一个因子的较高位

      coppersmith1

      p.bit_length() == 1024 ,p的高位需已知约576位
      p.bit_length() == 512 ,p的高位需已知约288位
      
    • 例子

      '''
      n = 110884890902749085253001083431222443088115610795940152564793628519927092107501946446399003764508722709710121804620193329162066855289179887539537634989483300155392790067446377224025966917227342075570751172611456818461296838516185655681858001119900898375522640670694604696426035782721144065487316221499661637517
      e = 65537
      c = 56d1b214082fd508567e0a4e101dcaa4f3edf262d7330cae4d75d94b874f53dfe3c9ba66d62a41b7e9331e67ae6907e3c028701e53555fea0832b2908471d04ceb98dbedf576a504902d50c3c32050fa036573de4f466f9c5de6b6bd4ad2f96bd6cd235a62c6c9555eb5ecf5b793b514f60d3e75a8307983c0f1aab746477a7b
      front = 754471047130831460574350468751127056146566410666010180184022324900851348720910487519
      backLength = 512 - front.bit_length()
      '''
          
      import binascii
      n=110884890902749085253001083431222443088115610795940152564793628519927092107501946446399003764508722709710121804620193329162066855289179887539537634989483300155392790067446377224025966917227342075570751172611456818461296838516185655681858001119900898375522640670694604696426035782721144065487316221499661637517
      cipher = 0x56c5afbc956157241f2d4ea90fd24ad58d788ca1fa2fddb9084197cfc526386d223f88be38ec2e1820c419cb3dad133c158d4b004ae0943b790f0719b40e58007ba730346943884ddc36467e876ca7a3afb0e5a10127d18e3080edc18f9fbe590457352dca398b61eff93eec745c0e49de20bba1dd77df6de86052ffff41247d
      e2 = 0x10001
      pbits = 512
      for i in range(0,4095):
        p4 = 0x636c1b2209b27268ad05ff5d64802c40d509cefccd92953227264dab0f27187dea4fdf000
        p4 = p4 + int(hex(i),16)
        kbits = pbits - p4.nbits() 
        p4 = p4 << kbits 
        PR.<x> = PolynomialRing(Zmod(n))
        f = x + p4
        roots = f.small_roots(X=2^kbits, beta=0.4) 
        if roots: 
          p = p4+int(roots[0])
          assert n % p == 0
          q = n/int(p)
          phin = (p-1)*(q-1)
          d = inverse_mod(e2,phin)
          print(d)
          cipher = 0x56d1b214082fd508567e0a4e101dcaa4f3edf262d7330cae4d75d94b874f53dfe3c9ba66d62a41b7e9331e67ae6907e3c028701e53555fea0832b2908471d04ceb98dbedf576a504902d50c3c32050fa036573de4f466f9c5de6b6bd4ad2f96bd6cd235a62c6c9555eb5ecf5b793b514f60d3e75a8307983c0f1aab746477a7b
          flag = pow(cipher,d,n)
          flag = hex(int(flag))[2:-1]
          print (binascii.unhexlify(flag))
      #https://sagecell.sagemath.org/
      
  • Known High Bits Message Attack

    • 例子
      /*
      [+]n=0x7c3139d3be9a691abdf3ff49c712fcb84ba39bbd2189bb98d04e04d2d7cc086c9d31b06fdf828aaeeb3765e1ab8ea41a3f1b8c73b80a498f1e2eaad42a1ac7b8e54e705cd1e3e4a39940f9bdcd16d4b42ab71a826955cc78450d6915663c82ae80fd2f64b7e3a70f2b188b85a738759eeb0688dfa22525bbbe92d7934763445L
      [+]e=3
      [+]m=random.getrandbits(512)
      [+]c=pow(m,e,n)=0x20084d9c4fa81d903437a9fabea4a2ad025a00ddc961e4fcd0f52ff9ec750702c109ce0188ae96e540a5c3dcf55013ced9ee37ad9547240fc8773f81fbb509b0b8ab24ed0288a6e1f997b5c0b196236bc8da2df9cce77c559492963eeafbbe4f5a9cb18098bfac87a1e179b26f60948fb72327acc0675890009a04697b76073L
      [+]((m>>72)<<72)=0xb90f972f73ebb3952b3a8e50233f783732478d874795b44c33f685caf7637f4cd0c90cf3a599e1a01e84a28459220b31a490fd1892df58000000000000000000L
      */
          
      import time
      def matrix_overview(BB, bound):
          for ii in range(BB.dimensions()[0]):
              a = ('%02d ' % ii)
              for jj in range(BB.dimensions()[1]):
                  a += '0' if BB[ii,jj] == 0 else 'X'
                  a += ' '
              if BB[ii, ii] >= bound:
                  a += '~'
              print a
      def coppersmith_howgrave_univariate(pol, modulus, beta, mm, tt, XX):
          
          dd = pol.degree()
          nn = dd * mm + tt
          
          if not 0 < beta <= 1:
              raise ValueError("beta should belongs in (0, 1]")
          
          if not pol.is_monic():
              raise ArithmeticError("Polynomial must be monic.")
             
          polZ = pol.change_ring(ZZ)
          x = polZ.parent().gen()
          
          # compute polynomials
          gg = []
          for ii in range(mm):
              for jj in range(dd):
                  gg.append((x * XX)**jj * modulus**(mm - ii) * polZ(x * XX)**ii)
          for ii in range(tt):
              gg.append((x * XX)**ii * polZ(x * XX)**mm)
          
          BB = Matrix(ZZ, nn)
          
          for ii in range(nn):
              for jj in range(ii+1):
                  BB[ii, jj] = gg[ii][jj]
          
          # display basis matrix
          if debug:
              matrix_overview(BB, modulus^mm)
          
          # LLL
          BB = BB.LLL()
          
          # transform shortest vector in polynomial    
          new_pol = 0
          for ii in range(nn):
              new_pol += x**ii * BB[0, ii] / XX**ii
          
          # factor polynomial
          potential_roots = new_pol.roots()
          print "potential roots:", potential_roots
          
          # test roots
          roots = []
          for root in potential_roots:
              if root[0].is_integer():
                  result = polZ(ZZ(root[0]))
                  if gcd(modulus, result) >= modulus^beta:
                      roots.append(ZZ(root[0]))
          return roots
          
          
      # RSA gen options (for the demo)
      length_N = 1024  # size of the modulus
      Kbits = 72      # size of the root
      e = 3
      N = 0x7c3139d3be9a691abdf3ff49c712fcb84ba39bbd2189bb98d04e04d2d7cc086c9d31b06fdf828aaeeb3765e1ab8ea41a3f1b8c73b80a498f1e2eaad42a1ac7b8e54e705cd1e3e4a39940f9bdcd16d4b42ab71a826955cc78450d6915663c82ae80fd2f64b7e3a70f2b188b85a738759eeb0688dfa22525bbbe92d7934763445L
      ZmodN = Zmod(N);
          
      C = 0x20084d9c4fa81d903437a9fabea4a2ad025a00ddc961e4fcd0f52ff9ec750702c109ce0188ae96e540a5c3dcf55013ced9ee37ad9547240fc8773f81fbb509b0b8ab24ed0288a6e1f997b5c0b196236bc8da2df9cce77c559492963eeafbbe4f5a9cb18098bfac87a1e179b26f60948fb72327acc0675890009a04697b76073L
      msg = 0xb90f972f73ebb3952b3a8e50233f783732478d874795b44c33f685caf7637f4cd0c90cf3a599e1a01e84a28459220b31a490fd1892df58000000000000000000
      P.<x> = PolynomialRing(ZmodN) #, implementation='NTL')
      pol = (msg + x)^e - C
      dd = pol.degree()
          
      # Tweak those
      beta = 1                                # b = N
      epsilon = beta / 7                      # <= beta / 7
      mm = ceil(beta**2 / (dd * epsilon))     # optimized value
      tt = floor(dd * mm * ((1/beta) - 1))    # optimized value
      XX = ceil(N**((beta**2/dd) - epsilon))  # optimized value
          
      # Coppersmith
      start_time = time.time()
      roots = coppersmith_howgrave_univariate(pol, N, beta, mm, tt, XX)
          
      # output
      print "\n# Solutions"
      print "we found:", str(roots)
      print("in: %s seconds " % (time.time() - start_time))
      print "\n"
          
      
  • 部分私钥泄露

    • 例子
      /*
      [+]n=0x291b24eae63660849a91b7122663814918ae91d62e3431163c4f47ecdbf92c59c9c430bbcc9443e4ff3dedbe60b1c06f383771bf628cdd36e649aa0c96db4addac4885071b651d2b1ae4e131ae3c115f1a59b828999ca7af8f235b75ad5b757680249eaa9b531ec1edbf9204417f17df08ec550893ed36523fcfef7fb4b2415dL
      [+]e=3
      [+]m=random.getrandbits(512)
      [+]c=pow(m,e,n)=0x623dc16f9047da92278d94fe3cabbd89db4f8c4c612ac55a439df31e368133d697cb08a571e2aad2a194800a433bc00940967441bb7e0d30bfc0599c55aeefc4af8be67ffaac307b65a2096863ca87c6aad615535814758212baae7328ac1ae9bce9f39a52456852c4c0b9779edbb19016872f516e2be9fab463f3b405e25beL
      [+]d=invmod(e,(p-1)*(q-1))
      [+]d&((1<<512)-1)=0x91d03d35338acebcf703991efd4b3f9c88e2f022568c31a410a33062d3e3e24571dc3537e21741e6b1c9eba127db0a768842d79a3197dca5b86e2cd509cd3b93L
      */
          
      known_bits = 512
      e = 3
      X = var('X')
      N=0x291b24eae63660849a91b7122663814918ae91d62e3431163c4f47ecdbf92c59c9c430bbcc9443e4ff3dedbe60b1c06f383771bf628cdd36e649aa0c96db4addac4885071b651d2b1ae4e131ae3c115f1a59b828999ca7af8f235b75ad5b757680249eaa9b531ec1edbf9204417f17df08ec550893ed36523fcfef7fb4b2415d
      d0 = 0x91d03d35338acebcf703991efd4b3f9c88e2f022568c31a410a33062d3e3e24571dc3537e21741e6b1c9eba127db0a768842d79a3197dca5b86e2cd509cd3b93
      P.<x> = PolynomialRing(Zmod(N))
      for k in xrange(1, e+1):
          results = solve_mod([e * d0 * X - k * X * (N - X + 1) + k * N == X], 2 ** 512)
          
          for m in results:
              f = x * 2 ** known_bits + ZZ(m[0])
              f = f.monic()
              roots = f.small_roots(X = 1, beta=0.3)
          
              if roots:
                  x0 = roots[0]
                  p = gcd(2 ** known_bits * x0 + ZZ(m[0]), N)
                  print '[+] Found factorization!'
                  print 'p =', ZZ(p)
                  print 'q =', N / ZZ(p)
                  break
      n=0x291b24eae63660849a91b7122663814918ae91d62e3431163c4f47ecdbf92c59c9c430bbcc9443e4ff3dedbe60b1c06f383771bf628cdd36e649aa0c96db4addac4885071b651d2b1ae4e131ae3c115f1a59b828999ca7af8f235b75ad5b757680249eaa9b531ec1edbf9204417f17df08ec550893ed36523fcfef7fb4b2415d
      p = 4369408607185874842987791687972458181281635894126489505104950532427588844072160412371716367300194382551703650763840526184308340422066926730909955032516327
      q = 6606303039610996668393006981258443682930645126368365403476569537191826726070101133886617447839914797821331788273937106747740172175833966059802001454391579
      e = 3
      phin = (p-1)*(q-1)
      d = inverse_mod(e,phin)
      cipher = 0x623dc16f9047da92278d94fe3cabbd89db4f8c4c612ac55a439df31e368133d697cb08a571e2aad2a194800a433bc00940967441bb7e0d30bfc0599c55aeefc4af8be67ffaac307b65a2096863ca87c6aad615535814758212baae7328ac1ae9bce9f39a52456852c4c0b9779edbb19016872f516e2be9fab463f3b405e25be
      print(d)
      flag = pow(cipher,d,n)
      flag = hex(int(flag))[2:-1]
      print(flag)
          
      
  • Coppersmith’s short-pad attack

    • 例子
      [+]n=0xc9f2c02d0ce22b192b5a046f8311b3eb470394ef228bbe8bc31f2939e3d7472a62eea2468c06b7d7de3a155a2e5a10c98143ede2fdf2f60fe5d65c9ba9fa26f5f7d05591201c76765599fb35f13e00a5b089fd4215c57b1453aaefc911a73c9f39003153af5e4a2e882a1c6c02d0024a6b0dede6c159a65b0bfe5c57b616127L
      [+]e=3
      [+]m=random.getrandbits(512)
      [+]c=pow(m,e,n)=0xb6046b56183fcc80d8a7c5dbc1f39176e736e2054255002abe1947a6e51fb7c37bdd689235613aec0e2a2651fade4837b968d4d6396b908a407f35e742065a773499f3bcd6111f2a1d8b65a3c79c9d3b20d681b9bf8cb2f26d2c528bca82e76d45ec734647cb13ca1a327e88173a64839bd4d8e576427600c86e7bc7224832cL
      [+]x=pow(m+1,e,n)=0xb590cc6da005f5bae916d26dca52f3f8e4c6c77d3d24df9f1f6e4e1ef1e58dc3b2bb0571810f5f27b019be2a768a392057c83006cbb12363b9661089d3fae650017c64d218ebe2b48b2ae91128d7613e6e51fabb94e7aaaba01d711d40ddac122683060ca5416ff0a00fa7f043f834d3989f8240b677a0cdda107832abe56c4L
          
      import gmpy2
      def getM2(a,b,c1,c2,n):
          a3 = pow(a,3,n)
          b3 = pow(b,3,n)
          first = c1-a3*c2+2*b3
          first = first % n
          second = 3*b*(a3*c2-b3)
          second = second % n
          third = second*gmpy2.invert(first,n)
          third = third % n
          fourth = (third+b)*gmpy2.invert(a,n)
          return fourth % n
      n=0xc9f2c02d0ce22b192b5a046f8311b3eb470394ef228bbe8bc31f2939e3d7472a62eea2468c06b7d7de3a155a2e5a10c98143ede2fdf2f60fe5d65c9ba9fa26f5f7d05591201c76765599fb35f13e00a5b089fd4215c57b1453aaefc911a73c9f39003153af5e4a2e882a1c6c02d0024a6b0dede6c159a65b0bfe5c57b616127L
      e=3
      c1=0xb6046b56183fcc80d8a7c5dbc1f39176e736e2054255002abe1947a6e51fb7c37bdd689235613aec0e2a2651fade4837b968d4d6396b908a407f35e742065a773499f3bcd6111f2a1d8b65a3c79c9d3b20d681b9bf8cb2f26d2c528bca82e76d45ec734647cb13ca1a327e88173a64839bd4d8e576427600c86e7bc7224832cL
      c2=0xb590cc6da005f5bae916d26dca52f3f8e4c6c77d3d24df9f1f6e4e1ef1e58dc3b2bb0571810f5f27b019be2a768a392057c83006cbb12363b9661089d3fae650017c64d218ebe2b48b2ae91128d7613e6e51fabb94e7aaaba01d711d40ddac122683060ca5416ff0a00fa7f043f834d3989f8240b677a0cdda107832abe56c4L
      padding1 = 0
      padding2 = 1
      m = getM2(1,padding1-padding2,c1,c2,n)-padding2
      print(hex(m))
          
      
  • Boneh and Durfee attack

    • 例子
      /*
      [+]n=0xbadd260d14ea665b62e7d2e634f20a6382ac369cd44017305b69cf3a2694667ee651acded7085e0757d169b090f29f3f86fec255746674ffa8a6a3e1c9e1861003eb39f82cf74d84cc18e345f60865f998b33fc182a1a4ffa71f5ae48a1b5cb4c5f154b0997dc9b001e441815ce59c6c825f064fdca678858758dc2cebbc4d27L
      [+]d=random.getrandbits(1024*0.270)
      [+]e=invmod(d,phin)
      [+]hex(e)=0x11722b54dd6f3ad9ce81da6f6ecb0acaf2cbc3885841d08b32abc0672d1a7293f9856db8f9407dc05f6f373a2d9246752a7cc7b1b6923f1827adfaeefc811e6e5989cce9f00897cfc1fc57987cce4862b5343bc8e91ddf2bd9e23aea9316a69f28f407cfe324d546a7dde13eb0bd052f694aefe8ec0f5298800277dbab4a33bbL
      [+]m=random.getrandbits(512)
      [+]c=pow(m,e,n)=0xe3505f41ec936cf6bd8ae344bfec85746dc7d87a5943b3a7136482dd7b980f68f52c887585d1c7ca099310c4da2f70d4d5345d3641428797030177da6cc0d41e7b28d0abce694157c611697df8d0add3d900c00f778ac3428f341f47ecc4d868c6c5de0724b0c3403296d84f26736aa66f7905d498fa1862ca59e97f8f866cL
      */
          
      import time
          
      strict = False
          
      helpful_only = True
      dimension_min = 7 # stop removing if lattice reaches that dimension
          
      def helpful_vectors(BB, modulus):
          nothelpful = 0
          for ii in range(BB.dimensions()[0]):
              if BB[ii,ii] >= modulus:
                  nothelpful += 1
          
          print nothelpful, "/", BB.dimensions()[0], " vectors are not helpful"
          
      # display matrix picture with 0 and X
      def matrix_overview(BB, bound):
          for ii in range(BB.dimensions()[0]):
              a = ('%02d ' % ii)
              for jj in range(BB.dimensions()[1]):
                  a += '0' if BB[ii,jj] == 0 else 'X'
                  if BB.dimensions()[0] < 60:
                      a += ' '
              if BB[ii, ii] >= bound:
                  a += '~'
              print a
          
      # tries to remove unhelpful vectors
      # we start at current = n-1 (last vector)
      def remove_unhelpful(BB, monomials, bound, current):
          # end of our recursive function
          if current == -1 or BB.dimensions()[0] <= dimension_min:
              return BB
          
          # we start by checking from the end
          for ii in range(current, -1, -1):
              # if it is unhelpful:
              if BB[ii, ii] >= bound:
                  affected_vectors = 0
                  affected_vector_index = 0
                  # let's check if it affects other vectors
                  for jj in range(ii + 1, BB.dimensions()[0]):
                      # if another vector is affected:
                      # we increase the count
                      if BB[jj, ii] != 0:
                          affected_vectors += 1
                          affected_vector_index = jj
          
                  # level:0
                  # if no other vectors end up affected
                  # we remove it
                  if affected_vectors == 0:
                      print "* removing unhelpful vector", ii
                      BB = BB.delete_columns([ii])
                      BB = BB.delete_rows([ii])
                      monomials.pop(ii)
                      BB = remove_unhelpful(BB, monomials, bound, ii-1)
                      return BB
          
                  # level:1
                  # if just one was affected we check
                  # if it is affecting someone else
                  elif affected_vectors == 1:
                      affected_deeper = True
                      for kk in range(affected_vector_index + 1, BB.dimensions()[0]):
                          # if it is affecting even one vector
                          # we give up on this one
                          if BB[kk, affected_vector_index] != 0:
                              affected_deeper = False
                      # remove both it if no other vector was affected and
                      # this helpful vector is not helpful enough
                      # compared to our unhelpful one
                      if affected_deeper and abs(bound - BB[affected_vector_index, affected_vector_index]) < abs(bound - BB[ii, ii]):
                          print "* removing unhelpful vectors", ii, "and", affected_vector_index
                          BB = BB.delete_columns([affected_vector_index, ii])
                          BB = BB.delete_rows([affected_vector_index, ii])
                          monomials.pop(affected_vector_index)
                          monomials.pop(ii)
                          BB = remove_unhelpful(BB, monomials, bound, ii-1)
                          return BB
          # nothing happened
          return BB
          
      """ 
      Returns:
      * 0,0   if it fails
      * -1,-1 if `strict=true`, and determinant doesn't bound
      * x0,y0 the solutions of `pol`
      """
      def boneh_durfee(pol, modulus, mm, tt, XX, YY):
          """
          Boneh and Durfee revisited by Herrmann and May
              
          finds a solution if:
          * d < N^delta
          * |x| < e^delta
          * |y| < e^0.5
          whenever delta < 1 - sqrt(2)/2 ~ 0.292
          """
          
          # substitution (Herrman and May)
          PR.<u, x, y> = PolynomialRing(ZZ)
          Q = PR.quotient(x*y + 1 - u) # u = xy + 1
          polZ = Q(pol).lift()
          
          UU = XX*YY + 1
          
          # x-shifts
          gg = []
          for kk in range(mm + 1):
              for ii in range(mm - kk + 1):
                  xshift = x^ii * modulus^(mm - kk) * polZ(u, x, y)^kk
                  gg.append(xshift)
          gg.sort()
          
          # x-shifts list of monomials
          monomials = []
          for polynomial in gg:
              for monomial in polynomial.monomials():
                  if monomial not in monomials:
                      monomials.append(monomial)
          monomials.sort()
              
          # y-shifts (selected by Herrman and May)
          for jj in range(1, tt + 1):
              for kk in range(floor(mm/tt) * jj, mm + 1):
                  yshift = y^jj * polZ(u, x, y)^kk * modulus^(mm - kk)
                  yshift = Q(yshift).lift()
                  gg.append(yshift) # substitution
              
          # y-shifts list of monomials
          for jj in range(1, tt + 1):
              for kk in range(floor(mm/tt) * jj, mm + 1):
                  monomials.append(u^kk * y^jj)
          
          # construct lattice B
          nn = len(monomials)
          BB = Matrix(ZZ, nn)
          for ii in range(nn):
              BB[ii, 0] = gg[ii](0, 0, 0)
              for jj in range(1, ii + 1):
                  if monomials[jj] in gg[ii].monomials():
                      BB[ii, jj] = gg[ii].monomial_coefficient(monomials[jj]) * monomials[jj](UU,XX,YY)
          
          # Prototype to reduce the lattice
          if helpful_only:
              # automatically remove
              BB = remove_unhelpful(BB, monomials, modulus^mm, nn-1)
              # reset dimension
              nn = BB.dimensions()[0]
              if nn == 0:
                  print "failure"
                  return 0,0
          
          # check if vectors are helpful
          if debug:
              helpful_vectors(BB, modulus^mm)
              
          # check if determinant is correctly bounded
          det = BB.det()
          bound = modulus^(mm*nn)
          if det >= bound:
              print "We do not have det < bound. Solutions might not be found."
              print "Try with highers m and t."
              if debug:
                  diff = (log(det) - log(bound)) / log(2)
                  print "size det(L) - size e^(m*n) = ", floor(diff)
              if strict:
                  return -1, -1
          else:
              print "det(L) < e^(m*n) (good! If a solution exists < N^delta, it will be found)"
          
          # display the lattice basis
          if debug:
              matrix_overview(BB, modulus^mm)
          
          # LLL
          if debug:
              print "optimizing basis of the lattice via LLL, this can take a long time"
          
          BB = BB.LLL()
          
          if debug:
              print "LLL is done!"
          
          # transform vector i & j -> polynomials 1 & 2
          if debug:
              print "looking for independent vectors in the lattice"
          found_polynomials = False
              
          for pol1_idx in range(nn - 1):
              for pol2_idx in range(pol1_idx + 1, nn):
                  # for i and j, create the two polynomials
                  PR.<w,z> = PolynomialRing(ZZ)
                  pol1 = pol2 = 0
                  for jj in range(nn):
                      pol1 += monomials[jj](w*z+1,w,z) * BB[pol1_idx, jj] / monomials[jj](UU,XX,YY)
                      pol2 += monomials[jj](w*z+1,w,z) * BB[pol2_idx, jj] / monomials[jj](UU,XX,YY)
          
                  # resultant
                  PR.<q> = PolynomialRing(ZZ)
                  rr = pol1.resultant(pol2)
          
                  # are these good polynomials?
                  if rr.is_zero() or rr.monomials() == [1]:
                      continue
                  else:
                      print "found them, using vectors", pol1_idx, "and", pol2_idx
                      found_polynomials = True
                      break
              if found_polynomials:
                  break
          
          if not found_polynomials:
              print "no independant vectors could be found. This should very rarely happen..."
              return 0, 0
              
          rr = rr(q, q)
          
          # solutions
          soly = rr.roots()
          
          if len(soly) == 0:
              print "Your prediction (delta) is too small"
              return 0, 0
          
          soly = soly[0][0]
          ss = pol1(q, soly)
          solx = ss.roots()[0][0]
          
          #
          return solx, soly
          
      def example():
          #
          # The problem to solve (edit the following values)
          #
          
          # the modulus
          N = 0xbadd260d14ea665b62e7d2e634f20a6382ac369cd44017305b69cf3a2694667ee651acded7085e0757d169b090f29f3f86fec255746674ffa8a6a3e1c9e1861003eb39f82cf74d84cc18e345f60865f998b33fc182a1a4ffa71f5ae48a1b5cb4c5f154b0997dc9b001e441815ce59c6c825f064fdca678858758dc2cebbc4d27
          # the public exponent
          e = 0x11722b54dd6f3ad9ce81da6f6ecb0acaf2cbc3885841d08b32abc0672d1a7293f9856db8f9407dc05f6f373a2d9246752a7cc7b1b6923f1827adfaeefc811e6e5989cce9f00897cfc1fc57987cce4862b5343bc8e91ddf2bd9e23aea9316a69f28f407cfe324d546a7dde13eb0bd052f694aefe8ec0f5298800277dbab4a33bb
          
          # the hypothesis on the private exponent (the theoretical maximum is 0.292)
          delta = .18 # this means that d < N^delta
          
          #
          # Lattice (tweak those values)
          #
          
          # you should tweak this (after a first run), (e.g. increment it until a solution is found)
          m = 4 # size of the lattice (bigger the better/slower)
          
          # you need to be a lattice master to tweak these
          t = int((1-2*delta) * m)  # optimization from Herrmann and May
          X = 2*floor(N^delta)  # this _might_ be too much
          Y = floor(N^(1/2))    # correct if p, q are ~ same size
          
          #
          # Don't touch anything below
          #
          
          # Problem put in equation
          P.<x,y> = PolynomialRing(ZZ)
          A = int((N+1)/2)
          pol = 1 + x * (A + y)
          
          #
          # Find the solutions!
          #
          
          # Checking bounds
          if debug:
              print "=== checking values ==="
              print "* delta:", delta
              print "* delta < 0.292", delta < 0.292
              print "* size of e:", int(log(e)/log(2))
              print "* size of N:", int(log(N)/log(2))
              print "* m:", m, ", t:", t
          
          # boneh_durfee
          if debug:
              print "=== running algorithm ==="
              start_time = time.time()
          
          solx, soly = boneh_durfee(pol, e, m, t, X, Y)
          
          # found a solution?
          if solx > 0:
              print "=== solution found ==="
              if False:
                  print "x:", solx
                  print "y:", soly
          
              d = int(pol(solx, soly) / e)
              print "private key found:", d
          else:
              print "=== no solution was found ==="
          
          if debug:
              print("=== %s seconds ===" % (time.time() - start_time))
          
      if __name__ == "__main__":
          example()
          
      

Tags:
0 comments



本作品采用知识共享署名-非商业性使用-禁止演绎 4.0 国际许可协议CC BY-NC-ND 4.0)进行许可。

This work is licensed under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License (CC BY-NC-ND 4.0).