Algorithm for finding partitionings of hard variants of boolean satisfiability problem with application to inversion of some cryptographic functions

In this paper we propose an approach for constructing partitionings of hard variants of the Boolean satisfiability problem (SAT). Such partitionings can be used for solving corresponding SAT instances in parallel. For the same SAT instance one can construct different partitionings, each of them is a set of simplified versions of the original SAT instance. The effectiveness of an arbitrary partitioning is determined by the total time of solving of all SAT instances from it. We suggest the approach, based on the Monte Carlo method, for estimating time of processing of an arbitrary partitioning. With each partitioning we associate a point in the special finite search space. The estimation of effectiveness of the particular partitioning is the value of predictive function in the corresponding point of this space. The problem of search for an effective partitioning can be formulated as a problem of optimization of the predictive function. We use metaheuristic algorithms (simulated annealing and tabu search) to move from point to point in the search space. In our computational experiments we found partitionings for SAT instances encoding problems of inversion of some cryptographic functions. Several of these SAT instances with realistic predicted solving time were successfully solved on a computing cluster and in the volunteer computing project SAT@home. The solving time agrees well with estimations obtained by the proposed method.

such that for any i, j : i � = j formula C ∧ G i ∧ G j is unsatisfiable and (hereinafter by "≡" we denote logical equivalence). Obviously when one has a partitioning of the original SAT instance, SAT for formulas C ∧ G j , j ∈ {1, . . . , s} can be solved independently in parallel.
There exist various partitioning techniques. For example one can construct {G j } s j=1 using a scattering procedure, a guiding path solver, lookahead solver and a number of other techniques described in Hyvärinen (2011). Unfortunately, for these partitioning methods it is hard in general case to estimate the time required to solve an original problem. From the other hand in a number of papers about SAT-based cryptanalysis of several keystream ciphers there was used a partitioning method that makes it possible to construct such estimations in quite a natural way. In particular, in Eibach et al. (2008), Soos et al. (2009), Soos (2010), Zaikin and Semenov (2008) for this purpose the information about the time to solve small number of subproblems randomly chosen from the partitioning of an original problem was used. In our paper we give strict formal description of this idea within the borders of the Monte Carlo method in its classical form (Metropolis and Ulam 1949). Also we focus our attention on some important details of the method that were not considered in previous works. Consider SAT for an arbitrary CNF C over a set of Boolean variables X = {x 1 , . . . , x n } . To an arbitrary set X = x i 1 , . . . , x i d , X ⊆ X we refer as a decomposition set. Consider a partitioning of C that consists of a set of 2 d formulas where G j , j ∈ {1, . . . , 2 d } are all possible minterms over X . Note that an arbitrary formula G j takes a value of true on a single truth assignment α we will refer as a decomposition family produced by X . It is easy to see that the decomposition family is the partitioning of the SAT instance C.
Let A be some SAT solving algorithm. Hereinafter we presume that A is complete, i.e. it halts on every input. We also presume that A is a non-randomized deterministic algorithm. We denote the total runtime of A on all the SAT instances from C X as t C,A X . Below we suggest a method for estimating t C,A X .
Define the uniform distribution on the set {0, 1} d . With each randomly chosen truth assignment (α 1 , . . . , α d ) from {0, 1} d we associate a value ξ C,A (α 1 , . . . , α d ) that is equal to the runtime of A on CNF C X /(α 1 , . . . , α d ) . Let ξ 1 , . . . , ξ Q be all the different values that ξ C,A (α 1 , . . . , α d ) takes on all the possible (α 1 , . . . , α d ) ∈ {0, 1} d . Below we use the following notation Denote the number of (α 1 , . . . , α d ), such that ξ C,A (α 1 , . . . , α d ) = ξ j , as ♯ξ j . Associate with (1) the following set We say that the random variable ξ C,A X has distribution P ξ C,A X . Note that the following equality holds Therefore, To estimate the expected value E ξ C,A X we will use the Monte Carlo method (Metropolis and Ulam 1949). According to this method, a probabilistic experiment that consists of N independent observations of values of an arbitrary random variable ξ is used to approximately calculate E[ξ ]. Let ζ 1 , . . . , ζ N be results of the corresponding observations. They can be considered as a single observation of N independent random variables with the same distribution as ξ. If E[ξ ] and Var(ξ ) are both finite then from the Central Limit Theorem (Feller 1971) we have the main formula of the Monte Carlo method Here σ = √ Var(ξ ) stands for a standard deviation, γ -for a confidence level, γ = � δ γ , where �(·) is the normal cumulative distribution function. It means that under the considered assumptions the value is a good approximation of E[ξ ], when the number of observations N is large enough.
Due to completeness of A the expected value and variance of random variable ξ C,A (X) are finite. Since A is deterministic (i.e. it does not use randomization) the observed values will have the same distribution. One can use the preprocessing stage to estimate the effectiveness of the considered partitioning because N can be significantly less than 2 d .
So the process of estimating the value (2) for a given X is as follows. We randomly choose N truth assignments of variables from X Below we refer to (4) as random sample. Then consider values ζ j = ξ C,A α j , j = 1, . . . , N and calculate the value If N is large enough then the value of F C,A X can be considered as a good approximation of (2). That is why one can search for a decomposition set with minimal value of F C,A (·) instead of finding a decomposition set with minimal value (2). Below we refer to function F C,A (·) as predictive function.

Algorithms for minimization of predictive function
Below we will describe the algorithm for finding good partitionings. This algorithm is based on the procedure minimizing the predictive function in the special search space.
Let C be an arbitrary CNF over the set of Boolean variables X = {x 1 , . . . , x n } . Let X ⊆ X be an arbitrary decomposition set. We can represent X by binary vector For an arbitrary χ ∈ {0, 1} n we compute the value of function F (χ) in the following way. For vector χ we construct the corresponding set X (it is formed by variables from X that correspond to 1 positions in χ). Then we construct a random sample α 1 , . . . , α N , α j ∈ {0, 1} |X| [see (4)] and solve SAT for CNFs C X /α j . For each of these SAT instances we measure ζ j -the runtime of algorithm A on the input C X /α j . After this we calculate the value of F C,A X according to (5). As a result we have the value of F (χ) in the considered point of the search space. Now we will solve the problem F (χ) → min over the set {0, 1} n . Of course, the problem of search for the exact minimum of function F (χ) is extraordinarily complex. Therefore our main goal is to find in affordable time the points in {0, 1} n with relatively good values of function F (·). Note that the function F (·) is not specified by some formula and therefore we do not know any of its analytical properties. That is why to minimize this function we use metaheuristic algorithms: simulated annealing and tabu search.
First, we need to introduce the notation. By R we denote the search space, for example, R = {0, 1} n , however, as we will see later, for the problems considered one can use the search spaces of much less power. During the minimization of function F (·) we iteratively move from one point of the search space to another: By N ρ (χ ) we denote the neighborhood of point χ of radius ρ in the search space R. The point from which the search starts we denote as χ start . We will refer to the decomposition set specified by this point as X start . The current Best Known Value of F (·) is denoted by F best . The point in which the F best was achieved we denote as χ best . By χ center we denote the point the neighborhood of which is processed at the current moment. We call the point, in which we computed the value F (·), a checked point. The neighborhood N ρ (χ ) in which all the points are checked is called checked neighborhood. Otherwise the neighborhood is called unchecked.
According to the scheme of the simulated annealing (Kirkpatrick et al. 1983), the transition from χ i to χ i+1 is performed in two stages. First we choose a point χ i from N ρ χ i . The point χ i becomes the point χ i+1 with the probability denoted as Pr χ i → χ i+1 |χ i . This probability is defined in the following way: In the pseudocode of the algorithm demonstrated below, the function that tests if the point χ i becomes χ i+1 , is called PointAccepted (this function returns the value of true if the transition occurs and false otherwise). The change of parameter T i corresponds to decreasing the "temperature of the environment" (Kirkpatrick et al. 1983) (in the pseudocode by decreaseTemperature() we denote the function which implements this procedure). Usually it is assumed that The process starts at some initial value T 0 and continues until the temperature drops below some threshold value T inf (in the pseudocode the function that checks this condition is called temperatureLimitReached()).

Algorithm 1: Simulated annealing algorithm for minimization of the predictive function
Input: CNF C, initial point χstart Output: Pair χ best , F best , where F best is a prediction for C, χ best is a corresponding decomposition set Another metaheuristic scheme that we used for minimization of F (·) is the tabu search algorithm (Glover and Laguna 1997). According to this algorithm we store the points from the search space, in which we already calculated the values of function F (·), in special tabu lists. When we try to improve the current Best Known Value of F (·) in the neighborhood of some point χ center then for an arbitrary point χ from the neighborhood we first check if we haven't computed F (χ) earlier. If we haven't and, therefore, the point χ is not contained in tabu lists, then we compute F (χ). This strategy is justified in the case of the minimization of predictive function F (·) because the computing of values of the function in some points of the search space can be very expensive. The use of tabu lists makes it possible to significantly increase the number of points of the search space processed per time unit. Let us describe the tabu search algorithm for minimization F (·) in more detail. To store the information about points, in which we already computed the value of F (·) we use two tabu lists L 1 and L 2 . The L 1 list contains only points with checked neighborhoods. The L 2 list contains checked points with unchecked neighborhoods. Below we present the pseudocode of the tabu search algorithm for F (·) minimization.
Algorithm 2: Tabu search altorithm for minimization of the predictive function In this algorithm the function markPointInTabuLists(χ , L 1 , L 2 ) adds the point χ to L 2 and then marks χ as checked in all neighborhoods of points from L 2 that contain χ. If as a result the neighborhood of some point χ ′ becomes checked, the point χ ′ is removed from L 2 and is added to L 1 . If we have processed all the points in the neighborhood of χ center but could not improve the F best then as the new point χ center we choose some point from L 2 . It is done via the function getNewCenter(L 2 ). To choose the new point in this case one can use various heuristics. In our current implementation the tabu search algorithm chooses the point for which the total conflict activity (Marques-Silva et al. 2009) of Boolean variables, contained in the corresponding decomposition set, is the largest.
As we already mentioned above, taking into account the features of the considered SAT problems makes it possible to significantly decrease the size of the search space. For example, knowing the so called Backdoor Sets (Williams et al. 2003) can help in that matter. Let us consider the SAT instance that encodes the inversion problem of the function of the kind f : {0, 1} k → {0, 1} l . Let S(f) be the Boolean circuit implementing f. Then the set X in , formed by the variables encoding the inputs of the Boolean circuit S(f), is the so called Strong Unit Propagation Backdoor Set (Järvisalo and Junttila 2009). It means that if we use X in as the decomposition set, then the CDCL (Conflict-Driven Clause Learning Marques-Silva et al. 2009) solver will solve SAT for any CNF of the kind C X in /α , α ∈ {0, 1} |X in | on the preprocessing stage, i.e. very fast. Therefore the set X in can be used as the set X start in the predictive function minimization procedure.
Moreover, in this case it is possible to use the set 2X in in the role of the search space R. In all our computational experiments we followed this path.

Computational experiments
We implemented the algorithms from the previous section in the form of PDSAT MPIprogram (https://github.com/Nauchnik/pdsat). One process of PDSAT is the leader process, all the other are computing processes (each process corresponds to 1 CPU core).
The leader process selects points of the search space (we use neighborhoods of radius ρ = 1). For every new point χ = χ X it generates a random sample (4) of size N. Each assignment from (4) combined with the original CNF C defines the SAT instance from the decomposition family C X . These instances are solved by computing processes. When computing the value of the predictive function we assume that the decomposition family will be processed by 1 CPU core. We can extrapolate the estimation obtained to an arbitrary parallel (or distributed) computing system because the processing of C X consists in solving independent subproblems. In the computing processes Mini-Sat solver (Eén and Sörensson 2003) is used. This solver was modified to be able to stop computations upon receiving corresponding messages from the leader process.
Below we present the estimations produced by PDSAT for SAT-based cryptanalysis of the A5/1 (Biryukov et al. 2000), Bivium (Cannière 2006) and Grain (Hell et al. 2007) keystream generators. We used the Transalg system (Otpuschennikov et al. 2015) to construct SAT instances for these problems.

Time estimations for SAT-based cryptanalysis of A5/1
For the first time the SAT-based cryptanalysis of the A5/1 keystream generator was considered in Semenov et al. (2011). Further we study this problem in the following form: to find the secret key of length 64 bits based on the given 114-bit keystream fragment. The PDSAT program was used to find partitionings with good time estimations for CNFs encoding this problem. The computational experiments were performed on the computing cluster "Academician V.M. Matrosov" of ISDCT SB RAS (http://hpc.icc.ru/index. php). One computing node of this cluster consists of 2 AMD Opteron 6276 CPUs (32 CPU cores in total). In each experiment PDSAT was launched for 1 day using 2 computing nodes (i.e. 64 CPU cores). We used random samples of size N = 10 4 .
On Figs. 1, 2a, b three decomposition sets are shown. We described the first decomposition set (further referred to as S 1 ) in the paper Semenov et al. (2011). This set (consisting of 31 variables) was constructed "manually" based on the analysis of algorithmic features of the A5/1 generator. The second one (S 2 ), consisting of 31 variables, was found as a result of the minimization of F (·) by the simulated annealing algorithm (see "Algorithms for minimization of predictive function" section). The third decomposition set (S 3 ), consisting of 32 variables, was found as a result of minimization of F (·) by the tabu search algorithm. In the Table 1 the values of F (·) (in seconds) for all three decomposition sets are shown. Note that each of decomposition sets S 2 and S 3 was found for one 114 bit fragment of keystream that was generated according to the A5/1 algorithm for a randomly chosen 64-bit secret key.

Solving cryptanalysis instances for A5/1 in the volunteer computing project SAT@home
The values of predictive function presented in Table 1 show that the SAT-based cryptanalysis of the A5/1 generator requires quite significant computing power. Specifically for the purpose of solving hard SAT instances we developed the volunteer computing project SAT@home (Posypkin et al. 2012   SAT@home was added to the official list of BOINC projects (http://boinc.berkeley.edu/ projects.php).
The experiment aimed at solving 10 cryptanalysis instances for the A5/1 keystream generator was held in SAT@home from December 2011 to May 2012. We used the rainbow-tables (http://opensource.srlabs.de/projects/a51-decrypt) to construct the corresponding instances. When analyzing 8 bursts of keystream (i.e. 912 bits) these tables allow to find the secret key with probability about 88%. We randomly generated 1000 instances and applied the rainbow-tables technique to analyze 8 bursts of keystream, generated by A5/1. Among these 1000 instances the rainbow-tables could not find the secret key for 125 problems. From these 125 instances we randomly chose 10 and in the computational experiments applied the SAT approach to the analysis of the first burst of each corresponding keystream fragment (114 bits). For each SAT instance we constructed the partitioning generated by the S 1 decomposition set (see Fig. 1) and processed it in the SAT@home project. All 10 instances constructed this way were successfully solved in SAT@home (i.e. we managed to find the corresponding secret keys) in about 5 months (the average performance of the project at that time was about 2 teraflops). The second experiment on the cryptanalysis of A5/1 was launched in SAT@home in May 2014. It was done with the purpose of testing the decomposition set found by tabu search algorithm. In particular we took the decomposition set S 3 (see Fig. 2b). On September 26, 2014 we successfully solved in SAT@home all 10 instances from the considered series.
It should be noted that in all the experiments the time required to solve the problem agrees with the predictive function value computed for the desomposition sets S 1 and S 3 . Our computational experiments clearly demonstrate that the proposed method of automatic search for decomposition sets makes it possible to construct SAT partitionings with the properties close to that of "reference" partitionings, i.e. partitionings constructed based on the analysis of algorithmic features of the considered cryptographic functions.

Time estimations for SAT-based cryptanalysis of Bivium and Grain
The Bivium keystream generator (Cannière 2006) is constructed from two shift registers of a special kind. The first one contains 93 cells and the second one contains 84 cells. The Grain keystream generator (Hell et al. 2007) also uses 2 shift registers: first is 80-bit nonlinear feedback shift register (NFSR), second is 80-bit linear feedback shift register (LFSR). To mix registers outputs the generator uses a special filter function h(x). In accordance with Maximov et al. (2007), Soos (2010) we considered cryptanalysis problems for Bivium and Grain in the following formulation. Based on the known fragment of keystream we search for the values of all registers cells at the end of the initialization phase. It means that we need to find 177 bits in case of Bivium and 160 bits in case of Grain.
Usually it is sufficient to consider keystream fragment of length comparable to the total length of shift registers to uniquely identify the secret key. Here we followed Eibach et al. (2008), Soos (2010) and set the keystream fragment length for Bivium cryptanalysis to 200 bits and for Grain cryptanalysis to 160 bits.
In our computational experiments we applied PDSAT to SAT instances that encode the cryptanalysis of Bivium and Grain according to the formulations described above.
In these experiments to minimize the predictive functions we used only the tabu search algorithm, since compared to the simulated annealing it traverses more points of the search space per time unit. Also we noticed that the decomposition set for the A5/1 cryptanalysis, constructed by the tabu search algorithm, is closer to the "reference" set than that constructed with the help of simulated annealing.
During the cryptanalysis of Bivium and Grain in the role of X start we used the set formed by the variables encoding the cells of registers of the generator considered at the end of the initialization phase. Further we refer to these variables as starting variables. Thus X start = 177 in case of Bivium, and X start = 160 in case of Grain. When computing predictive function values PDSAT used random samples of size N = 10 5 . It was launched for 1 day using 5 computing nodes (160 CPU cores in total) within the computing cluster "Academician V.M.Matrosov". So there was 1 leader process and 159 computing processes. Time estimations obtained are F best = 3.769 × 10 10 for Bivium and F best = 4.368 × 10 20 seconds for Grain. Corresponding decomposition set X best for Bivium is marked with gray on Fig. 3 (50 variables) and the decomposition set for Grain is marked with gray on Fig. 4 (69 variables). Interesting fact is that X best for Grain contains only variables corresponding to the LFSR cells.
In Eibach et al. (2008), Soos et al. (2009), Soos (2010) a number of time estimations for SAT-based cryptanalysis of Bivium were proposed. In particular, in Eibach et al. (2008) several fixed types of decomposition sets (strategies in the notation of Eibach et al.  (2008)) were analyzed. The best decomposition set from Eibach et al. (2008) consists of 45 variables encoding the last 45 cells of the second shift register. Note that in Eibach et al. (2008) the corresponding estimation of time equal to 1.637 × 10 13 was calculated using random samples of size 10 2 . In Soos et al. (2009), Soos (2010 the estimations of runtime for CryptoMiniSat solver, working with SAT instances encoding Bivium cryptanalysis, were presented. From the description of experiments in these papers it can be seen that authors used probabilistic experiment to estimate the sets of variables chosen by CryptoMiniSat during the solving process and extrapolated the estimations obtained to time points of the solving process that lay in the distant future. Note that in Soos et al. (2009), Soos (2010 the problem of estimating the effectiveness of a particular partitioning is not considered as the problem of estimating the expected value of some random variable (that is necessary for it to correspond to the Monte Carlo method in its classical sense). Apparently, as it is described in Soos et al. (2009), Soos (2010, the random samples of size 10 2 and 10 3 were used. In the Table 2 all three estimations mentioned above are shown. The performance of one CPU core we used in our experiments is comparable with that of one CPU core used in Soos et al. (2009), Soos (2010.

Solving cryptanalysis instances for Bivium and Grain
Since the values of predictive functions for Bivium and Grain cryptanalysis turned out to be quite large, in our computational experiments we studied "weakened" variants of the corresponding instances. For this purpose we used the sets of the so called "guessing bits" (Bard 2009). The instances obtained were solved on a computing cluster (with the help of PDSAT) and in the SAT@home project.
In the solving mode of PDSAT for X best found during predictive function minimization all 2 X best assignments of variables from X best are generated. PDSAT solves all corresponding SAT instances. To compare obtained time estimations with real solving time we used PDSAT to solve several cryptanalysis problems for Bivium and Grain with several known guessing bits. Below we use the notation BiviumK (GrainK) to denote the cryptanalysis of Bivium (Grain) with known K guessing bits. In the role of guessing bits in all cases we chose known values of K starting variables encoding the last K cells of the second shift register. We solved 3 instances for each of the following problems: Bivium16, Bivium14, Bivium12, Grain44, Grain42 and Grain40. In the following experiments for each BiviumK (GrainK) problem we computed the estimation for the first instance from the corresponding series and used the obtained decomposition set for all 3 instances from the series. To get more statistical data we did not stop the solving process after the satisfying assignment was found, thus processing the whole decomposition family. In the Table 3 for each problem we show the time required to solve it using 15 computing nodes (480 CPU cores total) of "Academician We also solved the Bivium9 problem in SAT@home. With the help of PDSAT the decomposition set formed of 43 variables was found. Using this decomposition set 5 instances of Bivium9 were solved in SAT@home in about 4 months from September 2014 to December 2014. During this experiment the average performance of the project was about 4 teraflops.
It should be noted that for all considered BiviumK and GrainK problems the time required to solve the corresponding instances on the computing cluster and in SAT@ home agrees well with values of the predictive function found by our approach.

Related work
Apparently, the paper Cook and Mitchell (1997) was the first work in which it was proposed to use SAT encodings of inversion problems of cryptographic functions as justified hard SAT instances. One of the first examples of SAT encodings for a widely known ciphering algorithm was proposed in Massacci and Marraro (2000): in particular, in that paper the process of constructing SAT encoding for the DES algorithm was described.
To the best of our knowledge, the first example of successful application of SAT solvers to cryptanalysis of real-world cryptographic functions was given in Mironov and Zhang (2006). It used the SAT solvers to construct collisions for the hash functions from the MD family.
The monograph (Bard 2009) contains systematic research of various questions regarding algebraic cryptanalysis. A substantial part of this book studies the possibilities of the use of SAT solvers for solving cryptanalysis equations represented in the form of algebraic systems over finite fields.
The A5/1 algorithm is still used in many countries to cipher GSM traffic. During the long lifetime of this algorithm a lot of attacks on it have been created. However, the first attacks that allowed to find the secret key in manageable time were presented by the A5/1 Cracking Project Group in 2009 [27]. These attacks were in fact developed from the Rainbow method (Oechslin 2003  one uses 8 bursts of keystream. The success rate of the Rainbow method if one has only 1 burst of keystream is about 24 %. In all our computational experiments we analyzed the keystream fragment of size 114 bits, i.e. one burst, and considered only instances for which the solution could not be found using the Rainbow method. We successfully solved in SAT@home several dozens of such instances. In Semenov et al. (2011) we described our first experience on the application of the SAT approach to A5/1 cryptanalysis in the specially constructed grid system BNB-Grid. In that paper we found the set S 1 (see "Time estimations for SAT-based cryptanalysis of A5/1" section) manually, based on the peculiarities of the A5/1 algorithm. The Bivium generator is a weakened variant of the Trivium generator (Cannière 2006) developed within the context of the eSTREAM project. The detailed analysis of its vulnerabilities was performed in Maximov et al. (2007). As far as we know, the cryptanalysis estimations from that paper were not verified with the exception of the distinguishing attack. Later the Bivium generator became quite a popular object of the SAT-based cryptanalysis. The paper Mcdonald et al. (2007) was the first research in that direction. In Eibach et al. (2008) there was described the SAT-based attack on Bivium, which used specially constructed sets of guessing bits. One of the advantages of Eibach et al. (2008) consists in the fact that their computational experiments are easy to reproduce. In Soos et al. (2009), Soos (2010 there was constructed a time estimation for the SAT-based cryptanalysis of Bivium, that was much better than all previous estimations. Essentially, to construct it the Monte Carlo method was used (in Soos 2010 the author even uses the term "Monte Carlo algorithm"). However that paper does not really contain any references to theoretical basics of the method: there is no formal definition of the random variable, the expected value of which is estimated. The main novelty of our approach consists in strict justification of the applicability of the Monte Carlo method to estimating the effectiveness of SAT partitionings, and in using metaheuristic algorithms (simulated annealing and tabu search) for finding partitionings with good estimations of total time required to process them.
The Monte Carlo method for estimating the expected value of a random variable was first proposed in Metropolis and Ulam (1949). There are a lot of modern guides and handbooks containing the description and the results of application of this method, for example, Kalos and Whitlock (1986).
Simulated annealling was first described in Kirkpatrick et al. (1983). It is used to solve optimization problems from various areas. Tabu search is another widely used metaheuristic method originated from Glover and Laguna (1997).
The questions regarding solving SAT in parallel and distributed environments were considered in a number of papers. In particular, in Hyvärinen (2011) a systematic review of methods for constructing SAT partitionings is presented.
The grid systems aimed at solving SAT are relatively rare. In Schulz and Blochinger (2010) a desktop grid for solving SAT which used conflict clauses exchange via a peer-topeer protocol was described. Apparently, Black and Bard (2011) became the first paper about the use of a desktop grid based on the BOINC platform for solving SAT. Unfortunately, it did not evolve into a full-fledged volunteer computing project. The predecessor of the SAT@home was the BNB-Grid system (Evtushenko et al. 2009;Semenov et al. 2011), that was used to solve first large scale SAT-based cryptanalysis problems in 2009. At the present moment there are several common principles that lie in the basis of modern SAT solvers. From many years of our experience we believe that in application to cryptanalysis instances the best solvers are the ones based on the CDCL concept (Marques-Silva et al. 2009). It might seem surprising that CDCL solvers show good results even when we solve inversion problems for functions with large number of preimages (for example, when we search for collisions of cryptographic hash functions). Nowadays there are many CDCL-solvers that have a common basic architecture but differ in details and heuristics.

Conclusion
In the present paper we propose the method for constructing SAT partitionings for solving hard SAT instances in parallel. This approach is based on the Monte Carlo method (in its classical form) for estimating expected value of random variable. From our point of view the proposed method and the corresponding algorithms can be used in SATbased cryptanalysis, that is an actively developing direction in cryptography. We tested our method in application to cryptanalysis of several keystream generators (A5/1, Bivium, Grain). In the nearest future we are going to expand the list of metaheuristics used for minimization of predictive functions. Also we plan to investigate the question of accuracy of the estimations obtained by the Monte Carlo method for the considered class of problems in more detail.