In this paper we study the problem of estimating the time required to process decompositions of hard SAT instances encoding inversion problems of some cryptographic functions. In particular, we consider one type of SAT decompositions, usually referred to as SAT partitioning. The effectiveness of a specific SAT parititioning is the total time required to solve all SAT instances from this partitioning. In the paper we construct statistical estimations of effectiveness of SAT partitionings with the help of computational scheme of the Monte Carlo method. We discuss the problem of accuracy of such estimations that arises because of drastic difference between the sizes of random samples, that can be processed in reasonable time, and the size of statistical population. We propose the method for improving constructed statistical estimations by using sets of random samples of increasing size followed by the extrapolation of obtained relation to the size of statistical population.
Copyright c by the paper’s authors. Copying permitted for private and academic purposes. The technologies that make it possible to solve such tests using SAT approach may prove useful for solving other combinatorial problems that can be effectively reduced to SAT.
In the present paper we describe an approach in which SAT solving algorithms are applied to the following problem. Assume we have a discrete function f : {0, 1}n → {0, 1}m which is total in {0, 1}n and is specified by some algorithm (program) A(f ) (by {0, 1}n, n ∈ N we denote the set formed by all binary words of length n). For an arbitrary y ∈ Range f , Range f ⊆ {0, 1}m we need to find such x ∈ {0, 1}n that f (x) = y. Below we will refer to this problem as the inversion problem for function f . It has an evident cryptographic context. For example, suppose that f is an algorithm implementing some keystream generator, y is a keystream fragment produced by the generator and x is an unknown secret key. Then the inversion problem for f is the problem of cryptanalysis of keystream generator based on the known keystream fragment.
One of the ways commonly employed to improve the effectiveness of solving
inversion problems for cryptographic functions is the use of parallel and distributed
computing environments. The choice of parallelization strategy that yields high effectiveness of
corresponding solving process is a very nontrivial task. For a number of functions with
secret key of small length the use of parallel computing makes it possible to solve
inversion problems via exhaustive search (also known as brute-force search or brute-force
attacks in the context of cryptography) [
In the present paper we study the problem of decomposing an original hard SAT instance into a family of easier SAT instanses, that can be solved independently of each other. Such families are called SAT partitionings. Partitionings can be constructed in many different ways, thus the time required to process different partitionings can greatly vary.
The problem of estimating the time of processing of an arbitrary SAT partitioning
is highly non-trivial in general case. Earlier [
In the present paper we propose an approach aimed at improving the accuracy of estimations constructed according to the Monte Carlo method by using a set of random samples of increasing size and then by extrapolating the obtained relation to the size of statistical population. We show that this approach works well in application to constructing SAT partitionings for SAT instances encoding inversion of several cryptographic functions. We applied this approach to improve the time estimations for solving inversion problems for some cryptographic functions. 2
Boolean Satisfiability Problem (SAT) consists in the following: for an arbitrary Boolean
formula F to answer the question whether it is satisfiable and to construct the satisfying
assignment in case if the answer is positive. This problem for an arbitrarty formula F
can be effectively (in polynomial time on the length of F binary code in general case)
reduced to the satisfiability problem for Boolean formula F ′ in conjunctive normal
form (CNF). Hereinafter we consider SAT in this context. SAT is a classical NP-hard
problem [
To reduce to SAT the inversion problem of f : {0, 1}m → {0, 1}m in an arbitrary
point y ∈ Rangef one can use various approaches, including automated systems that
construct SAT instance encoding the inversion of f in a considered point using the
text of a program specifying f . Such software systems employ the concept of symbolic
execution [
Let Cf (y) be a CNF, that encodes the inversion problem for some function f : {0, 1}n → {0, 1}m in point y ∈ Rangef . It means that Cf (y) is satisfiable and from each of its satisfying assignments one can effectively extract such set of Boolean values x ∈ {0, 1}n that f (x) = y. For the functions used in cryptography SAT for CNFs Cf (y) is usually quite hard. Essentially for ciphers used in practice it is pointless to believe that the problems of such kind can be solved without employing high performance computing using parallel computing architectures. Below we describe one general approach to parallelization of SAT usually referred to as Partitioning approach and study the problem of estimating the effectiveness of an arbitrary SAT partitioning, and also the problem of search for partitionings with good effectiveness estimations.
Let us assume that we consider SAT for an arbitrary CNF C. The SAT partitioning
of CNF C [
C ≡ (C ∧ C1) ∨ . . . ∨ (C ∧ Cs) .
It is clear that C is satisfiable if and only if at least one CNF of the kind C ∧ Ci, i ∈ {1, . . . , s} is satisfiable. It is important to note that SAT for these CNFs can be solved in parallel.
There is a lot of different approaches to constructing SAT partitionings. A
comprehensive study of many of them can be found in [
So, let C be an arbitrary CNF over the set of Boolean variables X = {x1, . . . , xn}. Consider an arbitrary set X˜ = {xi1 , . . . , xid }, X˜ ⊆ X . Let us call X˜ a decomposition set. Let g be an arbitrary minterm over X˜ and GX˜ = {gi}i=1,...,2d be the set of all possible minterms over X˜ . It is clear that an arbitrary formula g ∈ GX˜ is satisfiable only on a single truth assignment of variables from X˜ . Thus, the set
Δ(C, X˜ ) = {C ∧ g}g∈GX˜ is a SAT partitioning of CNF C. Hereinafter we work with partitionings of the type (1).
Let C be an arbitrary CNF and Δ(C, X˜ ) be its arbitrary SAT partitioning of the
type (1). Let A be an arbitrary SAT solving algorithm. By tA(C, X˜ ) we denote the
total time required to process all SAT instances from Δ(C, X˜ ) by algorithm A. In [
Theorem 1. Let A be an arbitrary complete SAT solving algorithm, i.e. its runtime is finite for an arbitrary input. Then for an arbitrary CNF C and any decomposition set X˜ there is a random variable ξA(C, X˜ ) with finite expected value E hξA(C, X˜ )i such that the following equality holds
tA(C, X˜ ) = 2|X˜ | · E hξA(C, X˜ )i .
Random variable ξ = ξA(C, X˜ ) from Theorem 1 is defined as follows. Define a uniform distribution on the set {0, 1}d, d = |X˜ |. With each assignment α ∈ {0, 1}d we link minterm g(α) over X˜ that is satisfiable on α. With each α randomly chosen from {0, 1}d we associate the value of random variable ξ equal to the time required by A to solve SAT for CNF C ∧ g(α).
Thus, if we know E[ξA(C, X˜ )] we can obtain exact value of time required to process
SAT partitioning Δ(C, X˜ ). Unfortunately, in general case we are not able to effectively
compute E[ξA(C, X˜ )]. However, we can estimate it using the Monte Carlo method
in its classical sence [
E[ξ] = E[ξ1] = . . . = E[ξN ], V ar(ξ) = V ar[ξ1] = . . . = V ar[ξN ]
Let E[ξ] and V ar(ξ) be both finite, also let γ be any fixed confidence level and
γ = Φ(δγ ), where Φ (·) is the normal cumulative distribution function. Then from
the Central Limit Theorem [
N · X ξj − E [ξ] < √ j=1 N ≥ γ.
Here σ = pV ar (ξ) stands for a standard deviation. It means that under the considered
N assumptions the value N1 · P ξj is a good approximation of E [ξ], when the number of j=1 observations N is large enough.
Taking into account all of the above, the procedure of statistical estimation of the value of tA(C, X˜ ) for a particular X˜ looks as follows. Assuming that there is specified a uniform distribution on {0, 1}d, d = |X˜ | we choose assignments α1, . . . , αN . The obtained set of assignments forms a random sample. For each αj, j ∈ {1, . . . , N } the algorithm A solves SAT for CNF C ∧ g(αj). The runtime of A is the value ξj of the observed random variable ξ. After solving all SAT instances of the kind C ∧ g(αj), j ∈ {1, . . . , N } we compute the following value: 1 N FA,C (X˜ ) = 2d · N · X ξj j=1
From (3) it follows that if N is large enough then the value of FA,C (X˜ ) can be considered as a good approximation of (2). Therefore, instead of searching for a decomposition set with minimal value (2) one can search for a decomposition set with minimal value of FA,C (X˜ ). Below we refer to function FA,C (X˜ ) as to a predictive function.
In [
The tabu search scheme for miminization of predictive function FA,C (X˜ ), that we
described in [
In this section we discuss the accuracy of estimations of the effectiveness of SAT
partitionings produced by the Monte Carlo method. State-of-the-art SAT-solvers are quite
complex programs, and even little modifications of their of parameters can greatly
influence its performance. Authors of [
The results of the papers mentioned above make it possible to conclude that the complex and inhomogenious behavior of the SAT solver can negatively influence the consistency of the estimations of the effectiveness of SAT partitionings produced using the Monte Carlo method. The main danger here lies in the risk of obtaining overly optimistic estimations. Below we propose one method for solving this problem and justify its adequacy using computational experiments.
So, what information can we use to estimate how well the sample mean approximates the expected value of the random variable considered? From the formula (3) it follows that the accuracy increases with the increase of the sample size (if expected value and variance are finite). In application to estimating the value of (2) the severity of the issue regarding the consistency of the Monte-Carlo estimations becomes apparent when we compare the random samples sizes with the total number of subproblems in the SAT partitioning. For example if the decomposition set consists of 50 variables, then the size of the corresponding partitioning is 250. Therefore the sample of size 10000 is approximately 8 · 10−10% of the total size of the statistical population. Can we consider the estimation produced with the use of such a small sample to be consistent?
The suggested approach for improving the estimations of effectiveness of SAT partitionings, obtained via Monte Carlo, is based on the following simple idea. We consider the problem of estimating E [ξ], assuming that random variable ξ has finite expected value and variance. Also we assume that the statistical population has finite size represented by number N ∗. Consider the sequence of random samples sizes N1 < N2 < . . . < Nl ≈ N ∗ (5) Below we will use the following notation:
1 FNi = Ni · Xj=1 ξj
Ni Here ξj jN=i1 is a random sample of size Ni from the statistical population of size N ∗. Hereinafter we assume that for each Ni a different random sample of corresponding size is considered. Suppose that we can in affordable time calculate the values FN1 , FN2, . . . , FNk , k < l. Based on this data we need to predict the behavior of function FN when N ∈ {Nk+1, . . . , Nl}. Below we study this problem in detail.
Unfortunately we cannot make any reasonable assumptions about the type of a function to be used to extrapolate FN . That is why we introduce the auxiliary function, to which we refer as jump function. This function is denoted as ε(N ) and is defined as follows: ε(Ni) = max
FNi−1 , FNi FNi
FNi−1 , i ∈ {2, . . . , l} In essence, each value of ε(Ni) shows how values of FN for N = Ni−1 and N = Ni differ from each other. Below we describe the behavior of function ε(N ) as N goes to infinity. It will allow us to choose the correct type of the function to extrapolate ε(N ) with. Let us show that the following theorem holds.
Theorem 2. Consider a probabilistic experiment in which the outcome is a random variable that takes positive values. Assume that this variable has finite expected value and variance. Consider the set
{N1, N2, . . . , Nl, . . .} such that Ni+1 = ⌊λ · Ni⌋ for each i and some constant λ > 1. Then for any confidence level γ there exists a constant α(γ) such that for some k0 = k0(γ) and all k ≥ k0 the following holds: α(γ) Pr 1 ≤ ε(Nk) ≤ 1 + √Nk ≥ γ2 Proof. Let ξ be a random variable satisfying all conditions of the theorem. The finiteness of the expected value and of the variance of ξ allows us to apply the Central Limit Theorem. Let us fix some confidence level γ. Assume a = E[ξ], b = δγ · σ. Since ξ takes only positive values, then a > 0. Consider an arbitrary set {N1, . . . , Nl, . . .} where Ni+1 = ⌊λ · Ni⌋, λ > 1 is some constant. Note, that the events Ak =
b |FNk − a| < √Nk , Ak+1 = (
FNk+1 − a <
b pNk+1 ) are independent and therefore from (3) for some k0 = k0(γ) and all k ≥ k0 the following holds:
Pr {Ak ∧ Ak+1} ≥ γ2 From (6) and from Nk+1 = ⌊λ · Nk⌋ for some λ > 1 it follows that with probability ≥ γ2 the following inequalities hold: ( a − √Nbk < FNk < a + √Nbk
a − √cN1 k < FNk+1 < a + √cN2 k where c1, c2 are some positive constants. Because ξ takes only positive values, then FNk > 0a,nFdNak+−1 √>cN1 k0.arSeinpcoesitEiv[ξe]. I>t m0etahnesnthfaotr fsoormaen ka0rbaitnrdaryalkl k≥ ≥k0 kfr0o mthe(7v)awluitehs b a − √Nk probability ≥ γ2 the following inequalities hold
b a − √Nk < a + √cN2 k
b FNk < a + √Nk FNk+1
a − √cN1 k 1 − a√NCk+c2 < FFNNk+k1 < 1 + a√NDk−c1 for some positive constants C, D. It is easy to construct similar bounds for FFNNk+k1 . Since it is evident that ε(Nk+1) ≥ 1 then from (8) the assertion of the theorem follows.
The theorem proved allows us to use the function of the kind 1 + √αN to predict the behavior of ε(N ). Here the constant α > 0 is picked up individually for each random variable ξ based on the known first values of ε(N ).
Now let us return to the problem of miminization of predictive function FA,C (X˜ ) over Boolean hypercube {0, 1}n. As we already noted, in fact we minimize function F (χ), the value of which in an arbitrary χ ∈ {0, 1}n is defined as follows. Ones in a vector χ specify decomposition set X˜ . For this X˜ we construct random sample α1, . . . , αN , αj ∈ {0, 1}d, d = |X˜ |. After this we compute the value of predictive function according to (4). The obtained value is the value of F (χ). Consider the following set of values of the parameter N (random sample size):
N1, N2 = λ · N1, . . . , Nl = λl−1 · N1 where λ, λ > 1, λ ∈ N is some constant. Also we suppose that Nl ≈ Δ(C, X˜ ) , where X˜ is a decomposition set represented by the vector χ. Let us fix some k∗, k∗ < l and compute values
FN1(χ), . . . , FNk∗ (χ) We will refer to this set as a training set. Our goal is to use the training set to compute values F˜Nk∗+1 (χ), . . . , F˜Nl (χ) that predict the behavior of function FN (χ) when N ∈ {Nk∗+1, . . . , Nl}.
Using (9) we compute ε(Ni), i ∈ {2, . . . , k∗} (7) (8) (9) We assume that (10) are the values of function of the type 1 + √αNi for corresponding Ni. Using this information we pick the value of the constant α. Denote the function obtained as ε˜(Ni). Let us compute the values of this function for Ni, i ∈ { ∗ k + 1, . . . , l}.
The values F˜Ni , i ∈ {k∗ + 1, . . . , l} we determine using the following formula ε˜(Ni) = max ( F˜Ni−1 , F˜Ni )
F˜Ni F˜Ni−1
, i ∈ {k∗ + 1, . . . , l}, F˜Nk∗ (χ) = FNk∗ (χ) When computing values of particular F˜Ni the main problem consists in avoiding overly optimistic estimations. That is why we choose to follow the “pessimistic scenario”, so for each i ∈ {k∗ + 1, . . . , l} we assume that
FNi+1 (χ) ≥ F˜Ni (χ) ˜ It is easy to see that in this case for any r > k∗ the following holds:
F˜Nr (χ) = FNk∗ (χ) · r
Y s=k∗+1 ε˜(Ns) (11) We consider the value F ∗(χ) = F˜Nl (χ) to be an improved estimation of the time required to process the partitioning Δ(C, X˜ ), where X˜ is a decomposition set defined by vector χ ∈ {0, 1}n.
The results of the application of the proposed technique aimed at improving the consistency of estimations are shown in the next section. 4
In this section we show the applicability of the suggested technique for improving the accuracy of estimations of predictive function values to SAT-based cryptanalysis of the Bivium cipher. 4.1
Time Estimations for Problems of SAT-based Cryptanalysis of
Bivium
The Bivium keystream generator [
In accordance with [
Therefore in our experiments we needed to find initial values of 177 Bivium
registers bits that lead to generation of known keystream fragment of size 200 bits. We
encoded the corresponding problem to SAT using the Transalg software system [
Further on the example of one instance of Bivium cryptanalysis we will show a typical situation, in which the sample size can critically affect the accuracy of the estimations. In the process of computing the random samples of size 105 were used. As a result PDSAT found the point χbest for which F (χbest) = 2.085 × 1011 seconds. The corresponding X˜best contains 50 variables. After this we launched the computational experiment in which we used random samples of size 104. After one day of work PDSAT found the point χ′best with predictive function value F (χ′best) = 5.98 × 106 seconds. We decided to compute in this point the values of predictive function using random samples of sizes N2 = 2 × 104, N3 = 4 × 104 and N4 = 8 × 104. For these random samples we also computed the values of jump function ε(N ). The corresponding results are presented in Table 1.
Sample size N1 = 104 FNi 5.98 × 106 ε(Ni) −
Note that ε(N4) = max n FFNN34 , FFNN43 o ≈ 3 × 105. The explanation of this phenomena is that samples of size N1 = 104, N2 = 2 × 104 and N3 = 4 × 104 contained only simple subproblems, resulting in overly optimistic estimations. However the sample of size N4 = 8 × 104 contained about a hundred of subproblems that were much harder than others. Therefore we can conclude that the estimation FN1 for sample size N1 = 104 is overly optimistic.
After this we returned to the point χbest found using random samples of size N1 = 105. Note that it defines the set X˜best, |X˜best| = 50. We computed the predictive function values in this point for random sample sizes
Nj = 2j−1 × 105, (12) where j ∈ {2, . . . , 14}. Obtained values FNj (χbest) we used as the training set. Then based on ε (Nj ), j ∈ {2, . . . , 14} we selected the value of constant α in the function εˆ(Nj ) = 1 + √αNj , that extrapolates the behavior of ε (Nj ) for j ∈ {15, . . . , 35}. Here X˜best = 50, so N35 ≈ 250. With the help of the Matlab system we got α = 445.9. The graphs of function ε (Nj ) for j ∈ {1, . . . , 14} and function εˆ(Nj ) for j ∈ {1, . . . , 35} are shown on Figure 1. Note that function εˆ(·) doesn’t have big jumps. On Figure 2 the extrapolation of the behavior of the predictive function in the considered point χbest constructed according to formula (11) is shown. The improved estimation for the predictive function value is Fˆ∗ (χbest) = 7.37 × 1011 2 1.5 n ispoe 1 l 0.5 0 9e+011 8e+011 7e+011 6e+011 snd 5e+011 ceo 4e+011 s 3e+011 2e+011 1e+011 0 4.00E+051.60E+066.40E+062.56E+071.02E+084.10E+081.64E+096.55E+092.62E+101.05E+114.19E+111.68E+126.71E+122.68E+131.07E+144.29E+141.13E+15
sample size Fig. 1. Values of ε and extrapolated ε˜ for Bivium
Fig. 2. Values of estimation FNj and extrapolated estimation F˜Nj for Bivium 4.2 Solving Cryptanalysis Problems for Weakened Ciphers and Justification of Proposed Method To evaluate the consistency of the proposed technique for improving estimations we solved a number of weakened cryptanalysis problems for the Bivium cipher using the computing cluster. For this purpose, the PDSAT program was used. In this program in addition to the prediction mode there is a solving mode. In this mode for X˜best found during predictive function minimization it generates all 2|X˜best| assignments of variables from X˜best and solves all corresponding SAT instances. By BiviumK we denote the families of weakened cryptanalysis problems of the Bivium cipher. We call each individual problem from such family an instance. An arbitrary instance from BiviumK is a cryptanalysis problem for the corresponding cipher in the formulation described above with an addition of known values of K initial variables that encode values of the last K cells of the second generator register. In the computational experiments we considered series of 3 different instances for Bivium14. In all cases on the prediction stage we used random samples of size N = 105. We processed these series in the following manner: for the first instance in the series of three we used PDSAT to find the decomposition set X˜best, defined by vector χbest. Then we computed the improved estimation Fˆ (χbest) for the corresponding set using the extrapolation via the jump function. After this the decomposition set X˜best found for the first instance was used to solve all 3 instances from the series in parallel. All instances were successfully solved. To solve one instance it took about 17 hours on 480 cores of Opteron 6276.
Thus for each problem of the kind Bivium14 we found exact time required to process the SAT partitioning defined by X˜best. After this we compared this value with the value obtained using the proposed technique for improving the estimated values of predictive function. Below we describe this experiment in more detail.
For each cryptanalysis instance Bivium14 the PDSAT system found the set X˜best containing 35 variables. Here PDSAT used random samples of size N1 = 105. Similar to the case of non-weakened Bivium we used λ = 2 and considered first 14 values Nj:
N1 = 105, Nj = 2j−1 × 105, j ∈ {2, . . . , 14}
The obtained values FNj (χbest) we used in the role of training set. Since |X˜best| = 35 and 235 ≈ N20 we needed to extrapolate the behaviour of functions ε(Nj ) and α .
FNj (χbest) to j ∈ {15, . . . , 20}. We extrapolate ε (Nj ) as a function of the type 1+ √Nj With the help of Matlab we got α = 72.273. The results of extrapolation are showed on Figure 3.
Then we extrapolated the behaviour of FNj (χbest) for j ∈ {15, . . . , 20} using formula (11). The results are shown on Figure 4.
The last column of the diagram on Figure 4 shows the difference between Fˆ∗(χbest) and the real time spent to process the corresponding SAT partitioning (solid line). Note that it does not exceed 7 % and Fˆ∗(χbest) is not overly optimistic estimation. 5
A lot of papers studied the questions regarding estimating the runtime of algorithms for
solving SAT and Constraint Satisfaction problem on hard instances. Thus in [
The first work that used SAT solvers for cryptanalysis was [
Note that in the cited papers neither random variables nor their expected values
estimated via Monte Carlo method were not strictly defined. It was first done in [
The main novelty of the results of the present article lies in the technique for improving the effectiveness estimations for SAT partitionings by means of extrapolating the values of predictive function to random samples of gradually increasing size. We showed the applicability of the corresponding technique on SAT instances encoding the inversion of some cryptographic functions.
Acknowledgments The research was funded by Russian Science Foundation (project No. 16-11-10046).