API tools faq
paste
bling_yu

Predict truncated Galois LFSRs

bling_yu
Jan 18th, 2024
229
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 21.32 KB | Source Code | 0 0
raw download clone embed print report
  1. '''
  2. The first experiment!
  3. This code is used to predict truncated Galois LFSRs with known parameters.
  4. Suppose that the modulus and the state transition matrix are known, given the first Ni (i=1,2,...,n) high-order truncated digits of the i-th sequence output by the i-th register, we need to recover the initial state.
  5. The library used is the Sagemath.
  6. '''
  7.  
  8. # Choose n1 registers from n registers and set the length of each sequence output by register. The total number of output digits is equal to d.
  9. def distribution(n, n1, d):
  10.     N = vector(ZZ, n1)
  11.     D = vector(ZZ, n)
  12.     C = Combinations(range(n), n1)
  13.     N = C.list()[randint(0, C.cardinality() - 1)]
  14.     for i in N:
  15.         D[i] = 1
  16.     for i in range(n1, d):
  17.         k = randint(0, n1 - 1)
  18.         D[N[k]] = D[N[k]] + 1
  19.     return N, D
  20.  
  21. # Recover the initial state by the high-order truncated digits
  22. def recover_the_initial_state_by_high_order_truncated_digits(m, n, d, n1, beta):
  23.     C = random_matrix(ZZ, n, n, x = 0, y = m)
  24.     A = random_matrix(ZZ, d, n, x = 0, y = m)
  25.     B = Matrix(ZZ, d, n)
  26.     Y = Matrix(ZZ, d, n)
  27.     N, D = distribution(n, n1, d)
  28.     n2 = n - n1
  29.     N2 = vector(ZZ, n2)
  30.     a = Matrix(ZZ, d + 2 * n2 + 1, n)
  31.    
  32.     i2 = 0
  33.     for i in range(0, n):
  34.         if i in N:
  35.             0
  36.         else:
  37.             N2[i2] = i
  38.             i2 = i2 + 1
  39.            
  40.     for i in range(1, d):
  41.         A[i] = (A[i - 1] * C) % m
  42.     for i in range(0, n):
  43.         for j in range(0, D[i]):
  44.             Y[j, i] = A[j, i] >> beta
  45.    
  46.     L = Matrix(ZZ, d + 2 * n2 + 1, d + 2 * n2 + 1)
  47.     L[0, 0] = 2^(beta - 1)
  48.     for i in range(1, n + 1):
  49.         L[0, i] = 2^(beta - 1)
  50.         L[i, i] = 1
  51.     for i in range(n + 1, n + n2 + 1):
  52.         L[0, i] = 2^(beta - 1)
  53.     col = n + n2
  54.    
  55.     for i in N:
  56.         CC = identity_matrix(ZZ, n)
  57.         for j in range(1, D[i]):
  58.             CC = (CC * C) % m
  59.             B[j] = (2^beta * (Y[0] * CC - Y[j])) % m
  60.             L[0, col + j] = -B[j, i] + 2^(beta - 1)
  61.             for k in range(0, n2):
  62.                 L[n + k + 1, col + j] = 2^beta * CC[N2[k], i]
  63.             L[col + j, col + j] = m
  64.             for k in range(0, n):
  65.                 L[k + 1, col + j] = CC[k, i]
  66.         col = col + D[i] - 1
  67.        
  68.     if beta >= alpha:
  69.         for i in range(n + 1, n + n2 + 1):
  70.             L[i, i] = 2^(beta - alpha)
  71.     else:
  72.         for i in range(0, d + 2 * n2 + 1):
  73.             for j in range(0, d + 2 * n2 + 1):
  74.                 L[i, j] = L[i, j] * 2^(alpha - beta)
  75.         for i in range(n + 1, n + n2 + 1):
  76.             L[i, i] = 1
  77.            
  78.     reduced_L = L.BKZ(block_size = 20)
  79.  
  80.     flag = 0
  81.     if beta >= alpha:
  82.         for i in range(0, d + 2 * n2 + 1):
  83.             if reduced_L[i, 0] == 2^(beta - 1):
  84.                 for j in range(1, n + 1):
  85.                     reduced_L[i, j] = -reduced_L[i, j] + 2^(beta - 1)
  86.                 for j in range(n + 1, n + n2 + 1):
  87.                     reduced_L[i, j] = (-reduced_L[i, j] + 2^(beta - 1)) * 2^alpha
  88.                 for j in range(0, n2):
  89.                     a[flag, N2[j]] = reduced_L[i, N2[j] + 1] + reduced_L[i, n + j + 1]
  90.                 for j in range(0, n1):
  91.                     a[flag, N[j]] = reduced_L[i, N[j] + 1] + Y[0, N[j]] * 2^beta
  92.                 flag = flag + 1
  93.             if reduced_L[i, 0] == -2^(beta - 1):
  94.                 for j in range(1, n + 1):
  95.                     reduced_L[i, j] = reduced_L[i, j] + 2^(beta - 1)
  96.                 for j in range(n + 1, n + n2 + 1):
  97.                     reduced_L[i, j] = (reduced_L[i, j] + 2^(beta - 1)) * 2^alpha
  98.                 for j in range(0, n2):
  99.                     a[flag, N2[j]] = reduced_L[i, N2[j] + 1] + reduced_L[i, n + j + 1]
  100.                 for j in range(0, n1):
  101.                     a[flag, N[j]] = reduced_L[i, N[j] + 1] + Y[0, N[j]] * 2^beta
  102.                 flag = flag + 1
  103.  
  104.     if beta < alpha:
  105.         for i in range(0, d + 2 * n2 + 1):
  106.             if reduced_L[i, 0] == 2^(alpha - 1):
  107.                 for j in range(1, n + 1):
  108.                     reduced_L[i, j] = (-reduced_L[i, j] + 2^(alpha - 1)) / 2^(alpha - beta)
  109.                 for j in range(n + 1, n + n2 + 1):
  110.                     reduced_L[i, j] = -reduced_L[i, j] + 2^(alpha - 1)
  111.                 for j in range(0, n2):
  112.                     a[flag, N2[j]] = reduced_L[i, N2[j] + 1] + reduced_L[i, n + j + 1] * 2^beta
  113.                 for j in range(0, n1):
  114.                     a[flag, N[j]] = reduced_L[i, N[j] + 1] + Y[0, N[j]] * 2^beta
  115.                 flag = flag + 1
  116.             if reduced_L[i, 0] == -2^(alpha - 1):
  117.                 for j in range(1, n + 1):
  118.                     reduced_L[i, j] = (reduced_L[i, j] + 2^(alpha - 1)) / 2^(alpha - beta)
  119.                 for j in range(n + 1, n + n2 + 1):
  120.                     reduced_L[i, j] = reduced_L[i, j] + 2^(alpha - 1)
  121.                 for j in range(0, n2):
  122.                     a[flag, N2[j]] = reduced_L[i, N2[j] + 1] + reduced_L[i, n + j + 1] * 2^beta
  123.                 for j in range(0, n1):
  124.                     a[flag, N[j]] = reduced_L[i, N[j] + 1] + Y[0, N[j]] * 2^beta
  125.                 flag = flag + 1
  126.                
  127.     return A[0], a
  128.  
  129. # Determine whether the initial state is recovered
  130. def judge(m, n, d, n1, beta):
  131.     A, a = recover_the_initial_state_by_high_order_truncated_digits(m, n, d, n1, beta)
  132.     if A == a[0]:
  133.         return 1
  134.     else:
  135.         return 0
  136.  
  137.  
  138. n = 16
  139. m = 2^31 - 1
  140.  
  141. alpha = 17
  142. n1 = 4
  143.  
  144. n2 = n - n1
  145. beta = 31 - alpha
  146. s = 100
  147. for d in range(31, 20, -1):
  148.     tol = 0
  149.     for i in range(0, s):
  150.         jud = judge(m, n, d, n1, beta)
  151.         if jud:
  152.             tol = tol + 1
  153.     print(d, tol)
  154.  
  155.  
  156.  
  157. '''
  158. The second experiment!
  159. This code is used to predict truncated Galois LFSRs with known parameters.
  160. Suppose that the modulus and the state transition matrix are known, given the first Ni (i=1,2,...,n) low-order truncated digits of the i-th sequence output by the i-th register, we need to recover the initial state.
  161. The library used is the Sagemath.
  162. '''
  163.  
  164. # Choose n1 registers from n registers and set the length of each sequence output by register. The total number of output digits is equal to d.
  165. def distribution(n, n1, d):
  166.     N = vector(ZZ, n1)
  167.     D = vector(ZZ, n)
  168.     C = Combinations(range(n), n1)
  169.     N = C.list()[randint(0, C.cardinality() - 1)]
  170.     for i in N:
  171.         D[i] = 1
  172.     for i in range(n1, d):
  173.         k = randint(0, n1 - 1)
  174.         D[N[k]] = D[N[k]] + 1
  175.     return N, D
  176.  
  177. # Recover the initial state by the low-order truncated digits
  178. def recover_the_initial_state_by_low_order_truncated_digits(m, n, d, n1, alpha, beta):
  179.     C = random_matrix(ZZ, n, n, x = 0, y = m)
  180.     A = random_matrix(ZZ, d, n, x = 0, y = m)
  181.     B = Matrix(ZZ, d, n)
  182.     Z = Matrix(ZZ, d, n)
  183.     N, D = distribution(n, n1, d)
  184.     n2 = n - n1
  185.     N2 = vector(ZZ, n2)
  186.     dim = d + 2 * n2 + 1
  187.     a = Matrix(ZZ, 2 * dim, n)
  188.    
  189.     i2 = 0
  190.     for i in range(0, n):
  191.         if i in N:
  192.             0
  193.         else:
  194.             N2[i2] = i
  195.             i2 = i2 + 1
  196.            
  197.     for i in range(1, d):
  198.         A[i] = (A[i - 1] * C) % m
  199.     for i in range(0, n):
  200.         for j in range(0, D[i]):
  201.             Z[j, i] = A[j, i] % (2^beta)
  202.    
  203.     gamma = gcd(m, 2^beta)
  204.     m_new = m / gamma
  205.      
  206.     L = Matrix(ZZ, dim, dim)
  207.     L[0, 0] = 2^(alpha - 1)
  208.     for i in range(1, n + 1):
  209.         L[0, i] = 2^(alpha - 1)
  210.         L[i, i] = 1
  211.     for i in range(n + 1, n + n2 + 1):
  212.         L[0, i] = 2^(alpha - 1)
  213.     col = n + n2
  214.    
  215.     for i in N:
  216.         CC = identity_matrix(ZZ, n)
  217.         for j in range(1, D[i]):
  218.             CC = (CC * C) % m
  219.             B[j] = (2^(-beta) * (Z[0] * CC - Z[j])) % m_new
  220.             L[0, col + j] = -B[j, i] + 2^(alpha - 1)
  221.             for k in range(0, n2):
  222.                 L[n + k + 1, col + j] = (2^(-beta) * CC[N2[k], i]) % m_new
  223.             L[col + j, col + j] = m_new
  224.             for k in range(0, n):
  225.                 L[k + 1, col + j] = CC[k, i] % m_new
  226.         col = col + D[i] - 1
  227.        
  228.     if alpha >= beta:
  229.         for i in range(n + 1, n + n2 + 1):
  230.             L[i, i] = 2^(alpha - beta)
  231.     else:
  232.         for i in range(0, dim):
  233.             for j in range(0, dim):
  234.                 L[i, j] = L[i, j] * 2^(beta - alpha)
  235.         for i in range(n + 1, n + n2 + 1):
  236.             L[i, i] = 1
  237.  
  238.     reduced_L = L.BKZ(block_size = 20)
  239.  
  240.     flag = 0
  241.     if alpha >= beta:
  242.         for i in range(0, dim):
  243.             if reduced_L[i, 0] == 2^(alpha - 1):                
  244.                 for j in range(1, n + 1):
  245.                     reduced_L[i, j] = -reduced_L[i, j] + 2^(alpha - 1)
  246.                 for j in range(n + 1, n + n2 + 1):
  247.                     reduced_L[i, j] = (-reduced_L[i, j] + 2^(alpha - 1)) / 2^(alpha - beta)
  248.                 for j in range(0, n2):
  249.                     a[flag, N2[j]] = reduced_L[i, N2[j] + 1] * 2^beta + reduced_L[i, n + j + 1]
  250.                 for j in range(0, n1):
  251.                     a[flag, N[j]] = reduced_L[i, N[j] + 1] * 2^beta + Z[0, N[j]]
  252.                 flag = flag + 1
  253.                 for j in range(0, n2):
  254.                     if reduced_L[i, N2[j] + 1] < 2^(alpha - 1):
  255.                         a[flag, N2[j]] = (reduced_L[i, N2[j] + 1] + 2^(alpha - 1)) * 2^beta + reduced_L[i, n + j + 1] - 1
  256.                     if reduced_L[i, N2[j] + 1] > 2^(alpha - 1):
  257.                         a[flag, N2[j]] = (reduced_L[i, N2[j] + 1] - 2^(alpha - 1)) * 2^beta + reduced_L[i, n + j + 1] + 1
  258.                 flag = flag + 1
  259.                              
  260.             if reduced_L[i, 0] == -2^(alpha - 1):
  261.                 for j in range(1, n + 1):
  262.                     reduced_L[i, j] = reduced_L[i, j] + 2^(alpha - 1)
  263.                 for j in range(n + 1, n + n2 + 1):
  264.                     reduced_L[i, j] = (reduced_L[i, j] + 2^(alpha - 1)) / 2^(alpha - beta)
  265.                 for j in range(0, n2):
  266.                     a[flag, N2[j]] = reduced_L[i, N2[j] + 1] * 2^beta + reduced_L[i, n + j + 1]
  267.                 for j in range(0, n1):
  268.                     a[flag, N[j]] = reduced_L[i, N[j] + 1] * 2^beta + Z[0, N[j]]
  269.                 flag = flag + 1
  270.                 for j in range(0, n2):
  271.                     if reduced_L[i, N2[j] + 1] < 2^(alpha - 1):
  272.                         a[flag, N2[j]] = (reduced_L[i, N2[j] + 1] + 2^(alpha - 1)) * 2^beta + reduced_L[i, n + j + 1] - 1
  273.                     if reduced_L[i, N2[j] + 1] > 2^(alpha - 1):
  274.                         a[flag, N2[j]] = (reduced_L[i, N2[j] + 1] - 2^(alpha - 1)) * 2^beta + reduced_L[i, n + j + 1] + 1
  275.                 flag = flag + 1
  276.                
  277.  
  278.     if alpha < beta:
  279.         for i in range(0, dim):
  280.             if reduced_L[i, 0] == 2^(beta - 1):
  281.                 for j in range(1, n + 1):
  282.                     reduced_L[i, j] = (-reduced_L[i, j] + 2^(beta - 1)) / 2^(beta - alpha)
  283.                 for j in range(n + 1, n + n2 + 1):
  284.                     reduced_L[i, j] = -reduced_L[i, j] + 2^(beta - 1)
  285.                 for j in range(0, n2):
  286.                     a[flag, N2[j]] = reduced_L[i, N2[j] + 1] * 2^beta + reduced_L[i, n + j + 1]
  287.                 for j in range(0, n1):
  288.                     a[flag, N[j]] = reduced_L[i, N[j] + 1] * 2^beta + Z[0, N[j]]
  289.                 flag = flag + 1
  290.                 for j in range(0, n2):
  291.                     if reduced_L[i, N2[j] + 1] < 2^(alpha - 1):
  292.                         a[flag, N2[j]] = (reduced_L[i, N2[j] + 1] + 2^(alpha - 1)) * 2^beta + reduced_L[i, n + j + 1] - 1
  293.                     if reduced_L[i, N2[j] + 1] > 2^(alpha - 1):
  294.                         a[flag, N2[j]] = (reduced_L[i, N2[j] + 1] - 2^(alpha - 1)) * 2^beta + reduced_L[i, n + j + 1] + 1
  295.                 flag = flag + 1
  296.                
  297.             if reduced_L[i, 0] == -2^(beta - 1):
  298.                 for j in range(1, n + 1):
  299.                     reduced_L[i, j] = (reduced_L[i, j] + 2^(beta - 1)) / 2^(beta - alpha)
  300.                 for j in range(n + 1, n + n2 + 1):
  301.                     reduced_L[i, j] = reduced_L[i, j] + 2^(beta - 1)
  302.                 for j in range(0, n2):
  303.                     a[flag, N2[j]] = reduced_L[i, N2[j] + 1] * 2^beta + reduced_L[i, n + j + 1]
  304.                 for j in range(0, n1):
  305.                     a[flag, N[j]] = reduced_L[i, N[j] + 1] * 2^beta + Z[0, N[j]]
  306.                 flag = flag + 1
  307.                 for j in range(0, n2):
  308.                     if reduced_L[i, N2[j] + 1] < 2^(alpha - 1):
  309.                         a[flag, N2[j]] = (reduced_L[i, N2[j] + 1] + 2^(alpha - 1)) * 2^beta + reduced_L[i, n + j + 1] - 1
  310.                     if reduced_L[i, N2[j] + 1] > 2^(alpha - 1):
  311.                         a[flag, N2[j]] = (reduced_L[i, N2[j] + 1] - 2^(alpha - 1)) * 2^beta + reduced_L[i, n + j + 1] + 1
  312.                 flag = flag + 1
  313.                
  314.     return A[0], a
  315.  
  316. # Determine whether the initial state is recovered
  317. def judge(m, n, d, n1, alpha, beta):
  318.     A, a = recover_the_initial_state_by_low_order_truncated_digits(m, n, d, n1, alpha, beta)
  319.     count = 0
  320.     for j in range(0, n):
  321.         if A[j] == a[0, j] or A[j] == a[1, j]:
  322.             count = count + 1
  323.     if count == n:
  324.         return 1
  325.     else:
  326.         return 0
  327.    
  328.  
  329. n = 16
  330. m = 2^31 - 2
  331.  
  332. beta = 17
  333. n1 = 4
  334.  
  335. n2 = n - n1
  336. alpha = 31 - beta
  337. s = 100
  338.  
  339. for d in range(32, 20, -1):
  340.     tol = 0
  341.     for i in range(0, s):
  342.         result = judge(m, n, d, n1, alpha, beta)
  343.         if result:
  344.             tol = tol + 1
  345.     print(d, tol)
  346.  
  347.  
  348.  
  349. '''
  350. The third experiment!
  351. This code is used to predict truncated Galois LFSRs with unknown parameters.
  352. Suppose that the modulus and the state transition matrix are unknown, given the first Ni (Ni>r, i=1,2,...,n) low-order truncated digits of the i-th sequence output by the i-th register, we need to recover the modulus, the characteristic polynomial, the initial state, and the state transition matrix.
  353. The library used is the Sagemath.
  354. '''
  355. import time
  356.  
  357. # Recover the modulus by the resultant    
  358. def recover_the_modulus(m, n, r, t, D, A, Y):
  359.     time_start = time.clock()
  360.     L = Matrix(ZZ, r, t + r)    
  361.     for i in range(0, r):
  362.         L[i, t + i] = 1
  363.         p = 0
  364.         for j in range(0, n):
  365.             for k in range(0, D[j] - r + 1):
  366.                 L[i, p + k] = Y[i + k, j]
  367.             p = p + D[j] - r + 1      
  368.     reduced_L = L.BKZ(block_size = 20)
  369.    
  370.     count = vector(ZZ, r)
  371.     for i in range(0, r):
  372.         yy = 0
  373.         for j in range(0, n):
  374.             for k in range(0, D[j] - r + 1):
  375.                 z = 0
  376.                 for zz in range(0, r):
  377.                     z = z + reduced_L[i, t + zz] * A[k + zz, j]
  378.                 if z == 0:
  379.                     yy = yy + 1
  380.                 if z != 0:
  381.                     yy = 0
  382.                     break
  383.             if yy == 0:
  384.                 break
  385.         if yy == t:
  386.             count[i] = count[i] + 1
  387.     time_end = time.clock()
  388.     print('Step 1: time of searching polynomials that annihilate the sequence: ', time_end - time_start, 's\n')
  389.    
  390.     time_start = time.clock()
  391.     R.<x> = ZZ[]
  392.     f = []
  393.     fi = reduced_L[0, t] * x^0
  394.     for i in range(1, r):
  395.         fi = fi + x^i * reduced_L[0, i + t]
  396.     f.append(fi)
  397.     fi = reduced_L[1, t] * x^0
  398.     for i in range(1, r):
  399.         fi = fi + x^i * reduced_L[1, i + t]
  400.     f.append(fi)
  401.     recovered_m_n = f[0].resultant(f[1])
  402.    
  403.     flag = 2
  404.     while recovered_m_n != m^n:
  405.         if flag == r - t:
  406.             flag = -1
  407.             break
  408.         fi = reduced_L[flag, t] * x^0
  409.         for i in range(1, r):
  410.             fi = fi + x^i * reduced_L[flag, i + t]
  411.         f.append(fi)
  412.         for i in range(0, flag):
  413.             res = f[flag].resultant(f[i])
  414.             recovered_m_n = gcd(recovered_m_n, res)
  415.         flag = flag + 1
  416.    
  417.     if flag > 0:
  418.         recovered_m = recovered_m_n^(1/n)
  419.         print("Success to recover m by ", flag, " polynomials!")
  420.         print("m == ", recovered_m, "\n")
  421.     else:
  422.         print("Fail to recover m!")
  423.         print("The wrong m == ", recovered_m_n^(1/n), "\n")
  424.    
  425.     time_end = time.clock()
  426.     print('Step 2: time of recovering the modulus: ', time_end - time_start, 's\n')
  427.     return count, reduced_L, flag
  428.  
  429. # Recover the coefficients by the Chinese Remainder Theorem
  430. def recover_the_coefficients(L, m, F, r, t):
  431.     time_start = time.clock()
  432.     R.<x> = PolynomialRing(Integers(m), implementation='NTL')
  433.    
  434.     for i in range(0, r):
  435.         mark = r - 1
  436.         while L[i, mark + t] == 0:
  437.             mark = mark - 1
  438.         for j in range(0, t + r):
  439.             L[i, j] = L[i, j] * ((L[i, mark + t]^(-1)) % m)
  440.    
  441.     f = []
  442.     fi = L[0, t]
  443.     for i in range(1, r):
  444.         fi = fi + x^i * L[0, i + t]
  445.     f.append(fi)
  446.     fi = L[1, t]
  447.     for i in range(1, r):
  448.         fi = fi + x^i * L[1, i + t]
  449.     f.append(fi)
  450.     g = f[0].gcd(f[1])
  451.    
  452.     flag = 2
  453.     while g != F:
  454.         if flag == r - t:
  455.             flag = -1
  456.             break
  457.         fi = L[flag, t]
  458.         for i in range(1, r):
  459.             fi = fi + x^i * L[flag, i + t]
  460.         f.append(fi)
  461.         for i in range(0, flag):
  462.             g = g.gcd(f[flag])
  463.         flag = flag + 1
  464.    
  465.     if flag > 0:
  466.         print("Success to recover f by ", flag, " polynomials!")
  467.         print("f == ", g, "\n")
  468.     else:
  469.         print("Fail to recover f!")
  470.         print("The correct f == ", F, "\n")
  471.         print("The wrong f == ", g, "\n")
  472.        
  473.     time_end = time.clock()
  474.     print('Step 3: time of recovering the coefficients: ', time_end - time_start, 's\n')
  475.     return flag
  476.  
  477. # Recover the initial state by Kannan's embedding technique
  478. def recover_the_initial_state(m, n, beta, F, Y, A):
  479.     time_start = time.clock()
  480.     flag = vector(ZZ, n)
  481.     for ni in range(0, n):
  482.         Q = Matrix(ZZ, n, n)
  483.         Q[0, n - 1] = (-F.list()[0]) % m
  484.         for i in range(1, n):
  485.             Q[i, i - 1] = 1
  486.             Q[i, n - 1] = (-F.list()[i]) % m
  487.         for d in range(25, 40):
  488.             Q_power = identity_matrix(ZZ, n)
  489.             recovered_a = vector(ZZ, n)
  490.             recovered_z = vector(ZZ, n)
  491.             L = Matrix(ZZ, d + 1, d + 1)
  492.             for i in range(1, n):
  493.                 Q_power = (Q_power * Q) % m
  494.             for i in range(n, d):
  495.                 Q_power = (Q_power * Q) % m
  496.                 b = 0
  497.                 for j in range(0, n):
  498.                     L[j + 1, i + 1] = Q_power[j, 0]
  499.                     b = b + Q_power[j, 0] * Y[j, ni]
  500.                 L[0, i + 1] = (2^beta * (Y[i, ni] - b)) % m + 2^(beta - 1)
  501.            
  502.             L[0, 0] = 2^(beta - 1)
  503.             for i in range(1, n + 1):
  504.                 L[0, i] = 2^(beta - 1)
  505.                 L[i, i] = 1
  506.             for i in range(n + 1, d + 1):
  507.                 L[i, i] = m
  508.            
  509.             reduced_L = L.BKZ(block_size = 20)
  510.             if reduced_L[0, 0] == -2^(beta - 1):
  511.                 for i in range(0, n):
  512.                     recovered_z[i] = reduced_L[0, i + 1] + 2^(beta - 1)
  513.                     recovered_a[i] = Y[i, ni] * 2^(beta) + recovered_z[i]
  514.             if reduced_L[0, 0] == 2^(beta - 1):
  515.                 for i in range(0, n):
  516.                     recovered_z[i] = -reduced_L[0, i + 1] + 2^(beta - 1)
  517.                     recovered_a[i] = Y[i, ni] * 2^(beta) + recovered_z[i]
  518.            
  519.             tol = 0
  520.             for i in range(0, n):
  521.                 if recovered_a[i] == A[i, ni]:
  522.                     tol = tol + 1
  523.             if tol == n:
  524.                 break
  525.                
  526.         flag[ni] = d
  527.     print("Success to recover the initial states by ", flag, " truncated digits!")
  528.     time_end = time.clock()
  529.     print('Step 4: time of recovering the initial states: ', time_end - time_start, 's\n')
  530.     return flag
  531.  
  532. # Recover the state transition matrix
  533. def recover_the_state_transition_matrix(m, n, A, C):
  534.     A0 = A[0:n]
  535.     A1 = A[1:n + 1]
  536.     recovered_C = (A0^(-1) * A1) % m
  537.     if recovered_C == C:
  538.         print("Success to recover the state transition matrix!")
  539.     else:
  540.         print("Fail to recover the state transition matrix!")
  541.     return recovered_C == C
  542.    
  543.    
  544. m = 2^31 - 1
  545. n = 16
  546. I = identity_matrix(ZZ, n)
  547. R.<x> = PolynomialRing(Integers(m))
  548. c = vector(ZZ, n)
  549. tol = 0
  550. r = 170
  551. t = 65
  552. beta = 14
  553. d = n * r + t - n
  554.  
  555. count = vector(ZZ, r)
  556. sum_count = vector(ZZ, r)
  557.  
  558. for i in range(0, 100):
  559.     time_start = time.clock()
  560.     print(i)
  561.     find_polynomial = 0
  562.     while find_polynomial == 0:
  563.         C = random_matrix(ZZ, n, n, x = 0, y = m)
  564.         F = x * I - C
  565.         find_polynomial = F.determinant().is_irreducible()
  566.    
  567.     di = r + t - n
  568.     A = random_matrix(ZZ, di, n, x = 0, y = m)
  569.     Y = Matrix(ZZ, di, n)
  570.     Z = Matrix(ZZ, di, n)
  571.     D = vector(ZZ, n)
  572.     for i in range(0, n):
  573.         D[i] = r
  574.     for i in range(n * r, d):
  575.         p = randint(1, n - 1)
  576.         D[p] = D[p] + 1
  577.        
  578.     for i in range(1, di):
  579.         A[i] = (A[i - 1] * C) % m
  580.     for i in range(0, n):
  581.         for j in range(0, D[i]):
  582.             Y[j, i] = A[j, i] >> beta
  583.             Z[j, i] = A[j, i] % (2^beta)
  584.    
  585.     count, reduced_L, flag_m = recover_the_modulus(m, n, r, t, D, A, Y)
  586.     flag_f = recover_the_coefficients(reduced_L, m, F.determinant(), r, t)
  587.     sum_count = count + sum_count
  588.     print("sum_count == ", sum_count, "\n")
  589.     if flag_m > 0 and flag_f > 0:
  590.         flag_a = recover_the_initial_state(m, n, beta, F.determinant(), Y, A)
  591.         flag_C = recover_the_state_transition_matrix(m, n, A, C)
  592.     time_end = time.clock()
  593.     print('Total time: ', time_end - time_start, 's\n\n')
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!