Arnaud Liefooghe
email: arnaud.liefooghe@lifl.fr

Laetitia Jourdan
email: laetitia.jourdan@lifl.fr

Matthieu Basseur

El-Ghazali Talbi
email: el-ghazali.talbi@lifl.fr

Edmund K Burke

Metaheuristics for the Bi-objective Ring Star Problem

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

The purpose of the bi-objective ring star problem is to find a cycle (the ring) which visits a subset of nodes of a graph. The two objectives are the minimization of a cost associated to the ring itself and the minimization of a cost associated to the arcs directed from non-visited nodes to visited ones. Although this problem is clearly bi-objective, it has always been investigated in a single-objective way. It was introduced by Labbé et al. [START_REF] Labbé | The ring star problem: Polyhedral analysis and exact algorithm[END_REF], where the goal was to minimize the sum of both costs. Another mono-criterion formulation of the problem, where one of the objectives is regarded as a constraint, has been investigated, for instance, by Renaud et al. [START_REF] Renaud | Efficient heuristics for median cycle problems[END_REF]. These two formulations are commonly used to convert a multiobjective problem into a single-objective one by using scalar approaches [START_REF] Miettinen | Nonlinear Multiobjective Optimization[END_REF].

The ring star problem has a wide range of industrial applications, including telecommunication networks design, school bus routing, routing of essential medical care services, circular-shaped transportation, and post-box location. However, in spite of its real-world applications, this is the first time that such a problem is studied in a bi-objective way, perhaps because of its complexity. Indeed, it is particularly challenging because, once it is decided which nodes have to be visited or not, a classical travelling salesman problem still remains to be solved. Nevertheless, Current and Schilling [START_REF] Current | The median tour and maximal covering tour problems: Formulations and heuristics[END_REF] investigated two variants of a similar problem: the median tour problem and the maximal covering tour problem. In both, one criterion is the minimization of the total length of the tour, while another one is the maximization of the access to the tour for non-visited nodes.

To tackle these problems, the authors used a kind of lexicographic method, where one objective function is optimized after another. Furthermore, Dorner et al. [START_REF] Doerner | Multicriteria tour planning for mobile healthcare facilities in a developing country[END_REF] recently formulated a problem of tour planning for mobile healthcare facilities in Senegal. A mobile facility has to visit a subset of nodes. Non-visited nodes are then assigned to their closest tour stop or are regarded as unable to reach a tour stop. The obectives are the minimization of the ratio between medical working time and total working time, the minimization of the average distance to the nearest tour stops, and the maximization of a coverage criterion. The authors designed a Pareto ant colony optimization algorithm and two multi-objective genetic algorithms to solve real-world instances.

In this paper, we investigate metaheuristic solution methods for the problem under consideration. Three metaheuristics are designed to approximate the whole set of efficient solutions. A population-based local search and two evolutionary algorithms are compared on state-of-the-art instances involving up to 300 nodes. The reminder of the paper is organized as follows. In Section 2, we provide the basic definitions for multi-objective optimization and the formulation of the biobjective ring star problem. Section 3 presents three resolution methods designed to tackle the problem under consideration. Some experimental results and a comparative study are provided in Section 4, while the last section concludes the paper and discusses perspectives about this work.

Preliminaries

This section presents some basic concepts related to multi-objective optimization and provides the formulation of the bi-objective ring star problem.

Multi-objective Optimization

A Multi-objective Optimization Problem (MOP) aims to optimize a set of n ≥ 2 objective functions f 1 , f 2 , . . . , f n simultaneously. Without loss of generality, we assume that all n objective functions have to be minimized. Let X denote the set of feasible solutions in the decision space, and Z the set of feasible points in the objective space. To each decision vector x ∈ X is assigned exactly one objective vector z ∈ Z on the basis of a vector function f : X → Z with z = f (x) = (f 1 (x), f 2 (x), . . . , f n (x)). In the case of a Multi-objective Combinatorial Optimization Problem (MCOP), note that a decision vector x ∈ X has a finite set of possible values. Definition 1. An objective vector z ∈ Z weakly dominates another objective vector z ′ ∈ Z if and only if ∀i ∈ [1..n], z i ≤ z ′ i . Definition 2. An objective vector z ∈ Z dominates another objective vector

z ′ ∈ Z if and only if ∀i ∈ [1..n], z i ≤ z ′ i and ∃j ∈ [1.
.n] such as z j < z ′ j . Definition 3. An objective vector z ∈ Z is non-dominated if and only if there does not exist another objective vector z ′ ∈ Z such that z ′ dominates z.

A solution x ∈ X is said to be efficient (or Pareto optimal ) if its mapping in the objective space results in a non-dominated point. The set of all efficient solutions is the efficient (or Pareto optimal ) set, denoted by X E . The set of all nondominated vectors is the non-dominated front (or the trade-off surface), denoted by Z N . A common approach in solving MOPs is to find or to approximate the set of efficient solutions; or at least a solution x ∈ X E for each non-dominated vector z ∈ Z N such as f (x) = z. A reasonable basic introduction to multiobjective optimization can be found in [START_REF] Deb | Multi-objective optimization[END_REF].

Note that we will assume, throughout the paper, that objective values are normalized. To achieve this, the minimum and the maximum value of each objective function are used in order to adaptively replace each objective function by its corresponding normalized function, so that its values lie in the interval [0, 1].

The Bi-objective Ring Star Problem

The Ring Star Problem (RSP) can be described as follows. Let G = (V, E, A) be a complete mixed graph where

V = {v 1 , v 2 , . . . , v n } is a set of vertices, E = {[v i , v j ]|v i , v j ∈ V, i < j} is a set of edges, and A = {(v i , v j )|v i , v j ∈ V } is a set of arcs.
Vertex v 1 is the depot. To each edge [v i , v j ] we assign a non-negative ring cost c ij , and to each arc (v i , v j ) we assign a non-negative assignment cost d ij . The RSP consists of locating a simple cycle through a subset of nodes V ′ ⊂ V (with v 1 ∈ V ′ ) while (i) minimizing the sum of the ring costs related to all edges that belong to the cycle, and (ii) minimizing the sum of the assignment costs of arcs directed from every non-visited node to a visited one so that the associated cost is minimum. An example of solution is given in Figure 1, where solid lines represent edges that belong to the ring and dashed lines represent arcs of the assignments. The first objective is called the ring cost and is defined as:

[vi,vj ]∈E c ij x ij , (1)
where x ij is a binary variable equal to 1 if and only if the edge [v i , v j ] belongs to the cycle. The second objective, the assignment cost, can be computed as follows:

vi∈V \V ′ min vj ∈V ′ d ij . ( 2
)
The RSP is an NP-hard combinatorial problem since the particular case of visiting the whole set of nodes is equivalent to a traditional travelling salesman problem.

Fig. 1. An example of a solution for the ring star problem.

3 Metaheuristics for the Bi-objective Ring Star Problem

Three metaheuristics are proposed to tackle the bi-objective RSP: a variable neighbourhood iterative Local Search (LS) and two Evolutionary Algorithms (EAs). These algorithms are respectively steady-state variations of IBMOLS [START_REF] Basseur | Indicator-based multi-objective local search[END_REF], IBEA [START_REF] Zitzler | Indicator-based selection in multiobjective search[END_REF] and NSGA-II [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF]. IBMOLS and IBEA are both recent indicator-based metaheuristics, whereas NSGA-II is one of the most often used Pareto-based resolution methods. In this section, RSP-specific components are described after we have presented the main characteristics of the LS and of the EAs.

A Multi-objective Local Search

Since they are easily adaptable to the multi-objective context, many of the search algorithms proposed to tackle MOPs are EAs. However, LS algorithms are known to be effective metaheuristics for solving real-world applications [START_REF] Burke | Search Methodologies: Introductory Tutorials in Optimization and Decision Support Techniques[END_REF][START_REF] Glover | Handbook of Metaheuristics[END_REF]. Several multi-objective LS approaches have been proposed in the literature. In particular, the Indicator-Based Multi-Objective Local Search (IBMOLS for short) has recently been presented in [START_REF] Basseur | Indicator-based multi-objective local search[END_REF]. IBMOLS is a generic population-based multiobjective LS dealing with a fixed population size. This allows to obtain a set of efficient solutions in a single run without specifying any mechanism to control the number of solutions during the search process. Moreover, IBMOLS represents an alternative to aggregation-and Pareto-based multi-objective metaheuristics. Indeed, as proposed in [START_REF] Zitzler | Indicator-based selection in multiobjective search[END_REF], it is assumed that the optimization goal is given in terms of a binary quality indicator I [START_REF] Zitzler | Performance assessment of multiobjective optimizers: An analysis and review[END_REF] which can be regarded as an extension of the Pareto dominance relation. A value I(A, B) quantifies the difference in quality between two non-dominated sets A and B. So, if Z N denotes the optimal non-dominated front, the overall optimization goal can be formulated as:

arg min A∈Ω I(A, Z N ) , (3)
where Ω denotes the space of all non-dominated set approximations. As noted in [START_REF] Zitzler | Indicator-based selection in multiobjective search[END_REF], Z N does not have to be known, it is just required in the formalization of the optimization goal. Since Z N is fixed, I actually represents a unary function that assigns a real number reflecting the quality of each approximation set.

One of the main advantages of indicator-based optimization is that no specific diversity preservation mechanism is generally required, according to the indicator being used.

The IBMOLS algorithm maintains a population P . Then, it generates the neighbourhood of a solution contained in P until a good solution is found (i.e. one that is better than at least one solution of P in terms of the indicator being used). By iterating this simple principle to every solution of P , we obtain a local search step. The whole local search stops when the archive of potentially efficient solutions has not received any new item during a complete local search step. Moreover, as local search methods are usually performed in an iterative way, a population re-initialization scheme has to be designed after each local search. Several strategies can be used within an iterative IBMOLS [START_REF] Basseur | Indicator-based multi-objective local search[END_REF]. Solutions can be re-initialized randomly, and crossover or random noise can be applied to solutions of the efficient set approximation. The interested reader could refer to [START_REF] Basseur | Indicator-based multi-objective local search[END_REF] for more details about IBMOLS.

A benefical feature of this LS is the low number of parameters that are required. In addition to the population size, the binary quality indicator to be used and the population re-initialization strategy (between each local search) are the two only other problem-independent parameters. Indeed, several quality indicators can be used within IBMOLS. The binary additive ǫ-indicator [START_REF] Zitzler | Indicator-based selection in multiobjective search[END_REF] is particularly well-adapted to indicator-based search and seems to be efficient on different kinds of problems (see, for instance, [START_REF] Basseur | Indicator-based multi-objective local search[END_REF][START_REF] Zitzler | Indicator-based selection in multiobjective search[END_REF]). It is capable of obtaining both a well-converged and a well-diversified Pareto set approximation. This indicator computes the minimum value by which a solution x 1 ∈ X can be translated in the objective space to weakly dominate another solution x 2 ∈ X. For a minimization problem, it is defined as follows:

I ǫ+ (x 1 , x 2 ) = max i∈{1,...,n} (f i (x 1 ) -f i (x 2 )) . (4)
Furthermore, to evaluate the quality of a solution according to a whole population P and a binary quality indicator I, different approaches exist. As proposed in [START_REF] Zitzler | Indicator-based selection in multiobjective search[END_REF], we will here consider an additive technique that amplifies the influence of solutions mapping to dominating points over solutions mapping to dominated ones which can be outlined as follows:

I(P \ {x}, x) = x ⋆ ∈P \{x} -e -I(x ⋆ ,x)/κ , (5)
where κ > 0 is a scaling factor. However, the initial experiments were not satisfactory because the algorithm was not able to find the extreme points of the trade-off surface. This is known to be one of the drawbacks of the ǫ-indicator, apparently due to the high convexity of the front. To tackle this problem, we add a condition preventing the deletion of solutions corresponding to the extreme non-dominated vectors during the replacement step of IBMOLS. Additionally, the population re-initialization scheme used between each local search is based on random noise, such as in the basic simulated annealing algorithm [START_REF] Burke | Search Methodologies: Introductory Tutorials in Optimization and Decision Support Techniques[END_REF]. Random noise consists of multiple mutations applied to N different randomly chosen solutions contained in the archive of potentially efficient solutions. If the size of the archive is less than N , the population is filled with random solutions.

Multi-objective Evolutionary Algorithms

The multi-objective EAs designed for the RSP are variations of two state-ofthe-art methods, namely IBEA [START_REF] Zitzler | Indicator-based selection in multiobjective search[END_REF] and NSGA-II [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF]. Some minor modifications have been carried out to improve the algorithms for the particular case of the addressed problem, for which the set of non-dominated points is, in general, very large.

IBEA Introduced by Zitzler and Künzli [START_REF] Zitzler | Indicator-based selection in multiobjective search[END_REF], the Indicator-Based Evolutionary Algorithm (IBEA) is, like IBMOLS, an indicator-based metaheuristic. The fitness assignment scheme of this EA is based on a pairwise comparison of solutions contained in a population by using a binary quality indicator. As within IBMOLS, no diversity preservation technique is required, according to the indicator being used. The selection scheme for reproduction is a binary tournament between randomly chosen individuals. The replacement strategy is an environmental one that consists of deleting, one-by-one, the worst individuals, and in updating the fitness values of the remaining solutions each time there is a deletion; this is continued until the required population size is reached. Moreover, an archive stores solutions mapping to potentially non-dominated points, in order to prevent their loss during the stochastic search process. However, in our case, and in contrast to the IBEA defined in [START_REF] Zitzler | Indicator-based selection in multiobjective search[END_REF], this archive is updated at each generation since the beginning of the EA, so that the output size is not necessarily less than or equal to the population size. Just like for the IBMOLS algorithm, the indicator used within IBEA is the additive ǫ-indicator; and the same mechanism has been used to prevent the loss of the extreme points on the trade-off surface.

NSGA-II At each generation of NSGA-II (Non-dominated Sorting Genetic Algorithm II [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF]), the solutions contained in the population are ranked into several classes. Individuals mapping to vectors from the first front all belong to the best efficient set; individuals mapping to vectors from the second front all belong to the second best efficient set; and so on. Two values are then computed for every solution of the population. The first one corresponds to the rank the corresponding solution belongs to, and represents the quality of the solution in terms of convergence. The second one, the crowding distance, consists of estimating the density of solutions surrounding a particular point of the objective space, and represents the quality of the solution in terms of diversity. A solution is said to be better than another if it has the best rank, or in the case of a tie, if it has the best crowding distance. The selection strategy is a deterministic tournament between two random solutions. At the replacement step, only the best individuals survive, with respect to the population size. Likewise, an external population is added to the steady-state NSGA-II in order to store every potentially efficient solution found during the search.

Application to the Bi-objective Ring Star Problem

This section provides the problem-specific steps of the metaheuristics introduced earlier. Components designed for the particular case of the bi-objective RSP, such as the encoding mechanism, the population initialization as well as the neighbourhood, mutation and crossover operators, are described below.

Solution Encoding

The encoding mechanism used to represent a RSP solution, for both the LS and the EAs, is based on the random keys concept proposed by Bean [START_REF] Bean | Genetic algorithms and random keys for sequencing and optimization[END_REF]. This implementation has already been successfully applied for a single-objective version of the RSP in [START_REF] Renaud | Efficient heuristics for median cycle problems[END_REF]. To each node v i belonging to the ring we assign exactly one random key x i ∈ [0, 1[, where x 1 = 0. A special value is assigned to unvisited nodes. Thus, the ring route associated to a solution corresponds to the nodes ordered according to their random keys in an increasing way; i.e. if x i < x j , then v j comes after v i . As an example, a possible representation for the cycle (v 1 ,v 7 ,v 4 ,v 9 ,v 2 ,v 6 ) is given in Figure 2. Vertices v 3 ,v 5 ,v 8 and v 10 are assigned to a visited node in such a way that the associated assignment cost is minimum.

Vertex v1 v2 v3 v4 v5 v6 v7 v8 v9 v10

Random key 0 0.7 -0.3 -0.8 0.2 -0.5 -Fig. 2. A RSP solution represented by random keys.

Population Initialization For every optimization method, the initial population has been generated randomly. Each node has a probability p = 0.5 that it will be visited or not, and to each visited vertex we associate a key randomly generated between 0 and 1.

Neighbourhood and Mutation Operators

As the RSP is both a routing problem and an assignment problem, different neighbourhood and mutation operators have to be designed. Here, we consider the following:

insert operator : adds an unvisited node v i in the cycle, the position where to insert v i is chosen in order to minimize the ring cost remove operator : removes a vertex v j of the ring -2-opt operator : applies a 2-opt operator between two nodes of the cycle v i and v j , i.e. it reverses the sequence of visited nodes between v i and v j .

For the LS, the neighbours of a solution are randomly explored, without considering any order between these three operators; and each neighbour is at most visited once. Moreover, note that it is not necessary to completely re-evaluate a solution each time a neighbourhood operator is applied. Thus, after an insert neighbourhood operator, we just have to re-assign unvisited nodes in order to minimize the assignment cost. After a remove neighbourhood operator, we just have to re-assign the nodes that were previously assigned to the one that has been removed. And, after a 2-opt neighbourhood operator, we just have to recompute the ring cost, the assignment cost being unchanged. In the case of mutations, the operators are applied to randomly chosen vertices.

Crossover Operator The crossover operator is a quadratic crossover closely related to the one proposed in [START_REF] Renaud | Efficient heuristics for median cycle problems[END_REF]. Two randomly selected solutions s 1 and s 2 are divided in a particular position. Then, the first part of s 1 is combined with the second part of s 2 to build a first offspring, and the first part of s 2 is combined with the second part of s 1 to build a second offspring. Every node retains its random key so that it enables an easy reconstruction of the new individuals. Thanks to the random keys encoding mechanism, solutions having a different ring size can easily be recombined, even if the initial ring structures are generally broken in the offspring solutions.

Experiments

The metaheuristics described in the previous section have all been implemented using the ParadisEO-MOEO library3 [START_REF] Liefooghe | ParadisEO-MOEO: A framework for evolutionary multi-objective optimization[END_REF]. ParadisEO-MOEO is a C++ whitebox object-oriented framework dedicated to the reusable design of metaheuristics for multi-objective optimization. All the algorithms share the same base components for a fair comparison between them. Computational runs were performed on an Intel Core 2 Duo 6600 (2 × 2.40 GHz) machine, with 2 GB RAM.

Experimental Protocol

Benchmarks The performance of the metaheuristics has been tested on different instances taken from the TSPLIB 4 [START_REF] Reinelt | TSPLIB -A traveling salesman problem library[END_REF]. These instances involve between 50 and 300 nodes. The number at the end of an instance's name represents the number of nodes for the instance under consideration. Let l ij denote the distance between two nodes v i and v j of a TSPLIB file. Then, the ring cost c ij and the assignment cost d ij have both been set to l ij for every pair of nodes v i and v j .

Parameter Setting For each one of the metaheuristics proposed to tackle the bi-objective RSP, the search process stops after a certain ammount of run time.

As shown in Table 1, this run time is defined according to the size of the instance under consideration. Likewise, the population size depends on the number of vertices involved in the instance (see Table 1). For each instance, a small (S), a medium (M), a large (L) and an extra-large (XL) population size have been tested. The noise rate for the population re-initialization in the iterated version of IBMOLS is set to a fixed percentage of the instance's size. We investigate three different values for this noise rate: 5%, 10% and 20%. Then, 0.05 × n, 0.1 × n and 0.2 × n random mutations are applied respectively for a problem with n nodes. For both IBMOLS and IBEA, the scaling factor κ is set to 0.05. Finally, for the EAs, the crossover probability is set to 0.25, and the mutation probability to 1.00, with a probability of 0.25, 0.25 and 0.50 for the remove, the insert and the 2-opt operator, respectively. Performance Assessment For each TSPLIB instance and each metaheuristic proposed in Section 3, a set of 20 runs, with different initial populations, has been performed. In order to evaluate the quality of the non-dominated front approximations obtained for a specific test instance, we follow the protocol given in [START_REF] Knowles | A tutorial on the performance assessment of stochastic multiobjective optimizers[END_REF]. First, we compute a reference set Z ⋆ N of non-dominated points extracted from the union of all these fronts. Second, we define z max = (z max ) denotes the upper bound of the first (respectively second) objective in the whole non-dominated front approximations. Then, to measure the quality of an output set A in comparison to Z ⋆ N , we compute the difference between these two sets by using the unary hypervolume metric [START_REF] Zitzler | Multiobjective evolutionary algorithms: A comparative case study and the strength pareto approach[END_REF], (1.05 × z max 1 , 1.05 × z max 2 ) being the reference point. The hypervolume difference indicator (I - H ) computes the portion of the objective space that is weakly dominated by Z ⋆ N and not by A. Furthermore, we also consider the R2 indicator proposed in [START_REF] Hansen | Evaluating the quality of approximations of the non-dominated set[END_REF] with a Chebycheff utility function defined by z ⋆ = (1, 1), ρ = 0.01 and a set Λ of 500 uniformly distributed normalized weighted vectors. As a consequence, for each test instance, we obtain 20 hypervolume differences and 20 R2 measures, corresponding to the 20 runs, per algorithm. As suggested by Knowles et al. [START_REF] Knowles | A tutorial on the performance assessment of stochastic multiobjective optimizers[END_REF], once all these values are computed, we perform a statistical analysis on pairs of optimization methods for a comparison on a specific test instance. To this end, we use the Mann-Whitney statistical test as described in [START_REF] Knowles | A tutorial on the performance assessment of stochastic multiobjective optimizers[END_REF], with a p-value lower than 5%. Note that all the performance assessment procedures have been achieved using the performance assessment tool suite provided in PISA5 [START_REF] Bleuler | PISA -a platform and programming language independent interface for search algorithms[END_REF].

Computational Results and Discussion

Table 2 presents the results obtained by the metaheuristics on eight different test instances. Due to space limitation and in order to simplify the reading of the table, only the results obtained by a large population size and by a noise rate of 5% for IBMOLS and by an extra-large population size for NSGA-II and IBEA are reported in the paper. These parameters have respectively been chosen as they were globally more efficient for each one of the algorithms. Overall, with respect to the metrics we used, we can see that IBMOLS performs significantly better than IBEA and NSGA-II on most test instances. Nevertheless, it is not the case on large problems (pr264 and pr299), where IBMOLS is outperformed by both algorithms according to the R2 metric. Additionally, although IBEA is in general statistically outperformed by IBMOLS, it performs significantly better than NSGA-II on a large number of the tested instances, and never performs significantly worse on each one of them. Furthermore, we can see that the overall efficiency of NSGA-II is very poor since it is statistically outperformed on most problems, except occasionally where it performs better than the IBMOLS algorithm, as pointed out above. To summarise, IBMOLS performs well on small-size RSP instances, but seems to have more trouble in dealing with large ones. Moreover, we also compared our results to the ones given in [START_REF] Labbé | The ring star problem: Polyhedral analysis and exact algorithm[END_REF] for a mono-objective version of the problem. For each test instance, the error ratio between the point belonging to Z ⋆ N that minimizes the single objective function investigated in [START_REF] Labbé | The ring star problem: Polyhedral analysis and exact algorithm[END_REF] and the exact optimal value is is always under 2% and is averagely under 0.5%.

One of the main characteristics of the problem under consideration seems to be the high number of points located in the trade-off surface. Then, after a certain number of iterations, a large part of the population involved in all the algorithms might map to potentially non-dominated points. This could explain the low efficiency of NSGA-II. Since the same rank is assigned to the major part of the population, only the crowding distance is used to compare solutions. However, the indicator-based fitness assignment scheme is obviously much more suited to determine potentially efficient solutions than the single crowding distance. Moreover, the high performance of IBMOLS in comparison to IBEA might depends on how close are the solutions which map to non-dominated points in the decision space. If these solutions are close to each other according to the neighbourhood operators, a LS is known to be particularly well-suited to find additional interesting solutions by exploring the neighbourhood of a potentially efficient solution. On the contrary, an EA usually explores the decision space in a more random way. Thus, a landscape analysis could be interesting to study the bi-objective RSP in more depth.

Conclusion

In this paper, a multi-objective routing problem, the bi-objective ring star problem, has been investigated. It has already been studied in a single-objective form where either both objectives have been combined [START_REF] Labbé | The ring star problem: Polyhedral analysis and exact algorithm[END_REF] or one objective has been treated as a constraint [START_REF] Renaud | Efficient heuristics for median cycle problems[END_REF]. Here, for the first time, this problem is formulated H and the R2 metrics by using a Mann-Whitney statistical test with a p-value of 5%. According to the metric under consideration, either the results of the algorithm located at a specific row are significantly better than those of the algorithm located at a specific column (≻), either they are worse (≺), or there is no significant difference between both (≡).

I - H R2 IBMOLS IBEA NSGA-II IBMOLS IBEA NSGA-II eil51 IBMOLS - ≻ ≻ - ≻ ≻ IBEA ≺ - ≻ ≺ - ≻ NSGA-II ≺ ≺ - ≺ ≺ - st70 IBMOLS - ≡ ≻ - ≻ ≻ IBEA ≡ - ≻ ≺ - ≻ NSGA-II ≺ ≺ - ≺ ≺ - kroA100 IBMOLS - ≡ ≻ - ≡ ≻ IBEA ≡ - ≻ ≡ - ≻ NSGA-II ≺ ≺ - ≺ ≺ - bier127 IBMOLS - ≻ ≻ - ≻ ≻ IBEA ≺ - ≻ ≺ - ≻ NSGA-II ≺ ≺ - ≺ ≺ - kroA150 IBMOLS - ≻ ≻ - ≻ ≻ IBEA ≺ - ≻ ≺ - ≻ NSGA-II ≺ ≺ - ≺ ≺ - kroA200 IBMOLS - ≻ ≻ - ≻ ≻ IBEA ≺ - ≡ ≺ - ≡ NSGA-II ≺ ≡ - ≺ ≡ - pr264 IBMOLS - ≡ ≻ - ≺ ≺ IBEA ≡ - ≻ ≻ - ≻ NSGA-II ≺ ≺ - ≻ ≺ - pr299 IBMOLS - ≡ ≻ - ≺ ≺ IBEA ≡ - ≻ ≻ - ≻ NSGA-II ≺ ≺ - ≻ ≺ -
in such a way that multiple conflicting criteria have to be optimized simultaneously. Three metaheuristics have been proposed to approximate the minimal complete set of efficient solutions: a population-based local search with variable neighbourhood and two evolutionary algorithms. Experiments were conducted using various test instances. We concluded that the local search method was significantly more efficient than the evolutionary algorithms on a large majority of instances, with respect to the performance metrics we used. The only instances for which the local search was outperformed were large ones. As a next step, we will try to solve ring star problem instances involving an even bigger number of nodes, to verify if our observations are still valid. If it is the case, it could be interesting to design a cooperation scheme between two different methods (i.e. the local search procedure and an evolutionary algorithm) in order to benefit from the respective quality of each one of them. The resulting hybrid metaheuristic could be particularly effective for solving large size problems.

1 , z max 2 ), where z max 1 (respectively z max 2

 1212

Table 1 .

 1 Instance-dependant parameters setting.

	Instance		Population size		Running
		S	M	L	XL	time
	eil51	5	10	15	100	20"
	st70	5	10	15	100	1'
	kroA100	10	15	20	100	2'
	bier127	10	15	20	100	5'
	kroA150	15	20	30	100	10'
	kroA200	15	20	30	100	20'
	pr264	15	20	30	100	50'
	pr299	15	20	30	100	100'

Table 2 .

 2 Comparison of different metaheuristics for the I -

ParadisEO is available at http://paradiseo.gforge.inria.fr.

The package is available at http://www.tik.ee.ethz.ch/pisa/assessment.html.


Charles Bouillaguet
email: charles.bouillaguet@univ-lille.fr

Florette Martinez
email: florette.martinez@lip6.fr

Julia Sauvage
email: julia.sauvage@etu.upmc.fr

John Von Neumann

Predicting the PCG Pseudo-Random Number Generator In Practice

Keywords: Pseudo-random number generator, guess-and-determine attack, truncated congruential generator, euclidean lattices, closest vector problem, practical attack

The Permuted Congruential Generators (PCG) are popular conventional (non-cryptographic) pseudo-random generators designed in 2014. They are used by default in the NumPy scientific computing package. Even though they are not of cryptographic strength, their designer stated that predicting their output should be nevertheless be "challenging". In this article, we present a practical algorithm that recovers all the hidden parameters and reconstructs the successive internal states of the generator. This enables us to predict the next "random" numbers, and output the seeds of the generator. We have successfully executed the reconstruction algorithm using 512 bytes of challenge input; in the worst case, the process takes 20 000 CPU hours. This reconstruction algorithm makes use of cryptanalytic techniques, both symmetric and lattice-based. In particular, the most computationally expensive part is a guessand-determine procedure that solves about 2 52 instances of the Closest Vector Problem on a very small lattice.

Introduction

Pseudo-random generators (PRG) are well-studied primitives in symmetric cryptography. A PRG is an efficient deterministic algorithm that stretch a small random seed into a longer pseudo-random stream. To achieve cryptographic-grade pseudo-randomness, a PRG must ensure that the pseudo-random stream is computationally indistinguishable from a "truly" random sequence of bits by efficient adversaries. Alternatively, it is possible to define pseudo-randomness by asking that no efficient algorithm is capable of predicting the next pseudo-random bit with non-negligible accuracy. The two definitions are in fact equivalent. It is well-known that pseudo-random generators can be turned into symmetric encryption algorithm, by generating "random" masks to be used in the one-time pad. This is precisely what stream ciphers do.

Not all pseudo-random generators are of cryptographic strength. In some applications, it is simply not necessary: to be used in Monte-Carlo numerical simulations or generate random choices in games, a relaxed, non-cryptographic notion of pseudo-randomness may be sufficient. This allows for faster algorithms. For instance, python standard library's random module uses the Mersenne Twister [START_REF] Matsumoto | Mersenne twister: A 623dimensionally equidistributed uniform pseudo-random number generator[END_REF]. The C library that comes along gcc (the glibc) uses a (poor) truncated linear congruential generator by default to implement the rand function.

In the realm of non-cryptographic random generators, a PRG is deemed "good enough" it is passes some efficient statistical tests -whereas the cryptographic notion of pseudorandomness asks that it passes any efficient test. There are de facto statistical test suites; an initial battery of randomness tests for RNGs was suggested by Knuth in the 1969 first edition of The Art of Computer Programming. In 1996, Knuth's tests were then supplanted by Marsaglia's Diehard tests. In 2007, L'Ecuyer proposed the TestU01 [LS07] library, whose "BigCrush" test is considered state-of-the-art by the relevant community, to the best of our knowledge. In 2010, the NIST proposed its own statistical test suite (Special Publication 800-22), to which improvements were later suggested [ZML + 16]. We understand that the PractRand test suite also has a good reputation.

In any case, designers of conventional pseudo-random generators try to obtain the simplest and fastest algorithm that passes the day's favourite test suite. A few selected ones are shown in Fig. 1. lehmer64 is a truncated linear congruential generator, touted as "the fastest PRNG that passes BigCrush" [START_REF] Lemire | The fastest conventional random number generator that can pass Big Crush?[END_REF]. xorshift128 is a clever implementation of a 128-bit LFSR with period 2 128 -1 due to Marsaglia [START_REF] Marsaglia | Xorshift RNGs[END_REF], using only a few simple 32-bit operations. xorshift128+ is a improved version due to Vigna [Vig17] that returns the sum of two consecutive outputs of a Xorshift LFSR; it passes the "BigCrush" test suite and is the default PRNG in many Javascript implementations, including that in Google's V8 engine (Chrome), Firefox and Safari.

Failures in cryptographic pseudo-random generators have catastrophic security implications. Let us mention for instance the well-known problem in Debian Linux from 2008, where a bug in the OpenSSL package led to insufficient entropy gathering and to practical attacks on the SSH and SSL protocols (the only remaining source of entropy comes from the PID of the process, i.e. 16 bits or less of effective entropy) [YRS + 09].

However, problems in non-cryptographic random number generators can also have dire consequences (barring the obvious case where they are used in lieu of their cryptographic counterparts). When they are used in scientific computing for Monte-Carlo methods, their defects have the potential to alter the results of numerical simulations. Ferrenberg et al. [START_REF] Ferrenberg | Monte carlo simulations: Hidden errors from "good" random number generators[END_REF] ran a classical Ferromagnetism Ising model Monte-Carlo simulation, in a special case where the exact results could be computed analytically, and compared the results of the simulation with the "true" answer. They used several conventional pseudo-random generators: a 32-bit linear congruential generator, two LFSRs, various combinations thereof, etc. They observed that changing the source of random numbers significantly altered the outcome of the numerical simulation. Different generators produced different biases: in particular a given LFSRs yielded energy levels that were systematically too low and critical temperatures that were always too high, while another kind of generator yielded the opposite (in many, repeated, trials).

The scientific computing community also realized that the need for fast parallel random number generation could be satisfied by the use of block ciphers in counter mode [START_REF] Salmon | Parallel random numbers: As easy as 1, 2, 3[END_REF]. The need for speed then leads to the use of weakened cryptographic primitives (roundreduced AES or custom and presumably weak block-ciphers)

In most cases, it is fairly easy to see that a given conventional PRG does not meet the cryptographic notion of pseudo-randomness, and there are few exceptions. Most are fairly easy to predict, meaning that after having observed a prefix of the output, it is easy to produce the next "pseudo-random" bits. This makes a good source of exercises for cryptology students.

In this paper, we study the PCG family of non-cryptographic pseudo-random generators proposed by O'Neil [O'N14b, O'N14a]. She did not claim that the algorithm has cryptographic strength, but that predicting its output ought to be "challenging". We therefore took up the challenge.

PCG stands for "Permuted Congruential Generator": it essentially consists in applying a non-linear filtering function on top of a linear congruential generator (in a way reminiscent to the venerable filtered LFSRs). The resulting combination is fast and passes current statistical test suites. The PCG family contains many members, but we focus on the strongest one, named either PCG64 or PCG-XSL-RR. It has a 128-bit internal state and produces 64 bits when clocked. It is the default pseudo-random number generator in the popular NumPy [START_REF] Van Der Walt | The numpy array: A structure for efficient numerical computation[END_REF] scientific computing package for Python.

The internal state of the PCG64 generator is made of a 128-bit "state" and a 128-bit "increment", whose intended use is to provide several pseudo-random streams with the same seed (just as the initialisation vectors do in stream ciphers). A default increment is provided in case the end-user just want one pseudo-random stream with a single 128-bit seed.

Contribution.

We describe an algorithm that reconstructs the full internal state of the strongest member of the PCG family. This allows to predict the pseudo-random stream deterministically and clock the generator backwards. The original seeds can also easily be reconstructed. The state reconstruction algorithm is practical and we have executed it in practice. It follows that predicting the output of the PCG should be considered practically feasible.

While the PCG pseudo-random generator is not meant as a cryptographic primitive, obtaining an actual prediction algorithm requires the use of cryptanalytic techniques. Making it practical requires in addition a non-trivial implementation effort.

Our algorithm reconstruct the internal state using a "guess-and-determine" approach: some bits of the internal state are guessed ; assuming the guesses are correct, some other information is computed ; a consistency check discards bad guesses early on ; then candidate internal states are computed and fully tested. The problem actually come in two distinct flavors.

When the increment is known (for instance when it is the default value), a simplified prediction algorithm recovers the internal state from 192 bits of pseudo-random stream. The process runs in 20 CPU minutes. It guesses 37 bits of the internal state, then solves an instance of the Closest Vector Problem (CVP) in a 3-dimensional euclidean lattice. This requires about 50 arithmetic operations in total and reveals the entire internal state if the guesses are correct.

When the increment is unknown, things are a bit more complicated. This is the default situation in NumPy, where both the state and the increment are initialised using an external source of entropy. In this case, our prediction algorithm requires 4096 bits of pseudo-random stream ; it guesses between 51 and 55 bits, then for each guess it solves an instance of CVP in dimension 4 (using about 75 arithmetic operations). This recovers 64 more bits of information about the difference between two successive states, and this is enough to filter the bad guesses. This information can then be used in a subsequent and comparably inexpensive phase to recover the entire internal state. On average, the whole process requires a bit less than 20 000 CPU hours to complete.

We implemented our algorithms, then asked the designer of the PCG family to send us "challenge" pseudo-random streams ; we ran our code and emailed back the (correct) seeds used to generate the challenge streams the next day.

Related Work. Deterministic pseudo-random generators can be traced back to the work of Von Neumann and Metropolis on the ENIAC computer [vN51]; they suggested around 1946 to use the "middle-square" method: if u n is a k-digit number, form u n+1 by taking the square of the k 2 middle digits of u n . This is a venerable precursor of the Blum-Blum-Shub "provably secure" PRNG. The main problem of this method is that it produces sequences that quickly enters short cycles.

Lehmer later proposed linear congruential generators in 1949, also for use on the ENIAC computer [START_REF] Lehmer | Mathematical methods in large-scale computing units[END_REF]. He gave the sequence defined by u 0 = 47594118, u n+1 = 23u n mod 10 8 + 1 and proved that it had period 5882353, a clear improvement compared to the middle-square approach. More details on early pseudo-random generators can be found in [START_REF] Knuth | of The Art of Computer Programming[END_REF].

Knuth discussed whether truncated linear congruential generators could be good stream ciphers; he therefore studied the problem of recovering the internal state of a truncated linear congruential generator x i+1 = ax i + c mod 2 k when a and c are unknown [START_REF] Knuth | Deciphering a linear congruential encryption[END_REF]; he gave an algorithm exponential in the number of truncated bits.

Boyar studied further the problem [START_REF] Boyar | Inferring sequences produced by a linear congruential generator missing low-order bits[END_REF] and presented an algorithm which could predict a linear congruential generator when all the parameters (multiplier, increment, modulus and initial state) are unknown; she extended her idea to the case of truncated linear congruential generators, under the condition that the number of bit unrevealed is really small in comparison to the size of the modulus.

Frieze et al [FHK + 88] improved the efficiency of reconstruction algorithms in simpler cases. They supposed that the multiplier a and the modulus 2 k were known and used lattice-based techniques to recover a truncated linear congruential generator with more truncated bits.

Later on, Joux and Stern extended this result to the case where the multiplier a and the modulus 2 k are unknown, also using lattice techniques [START_REF] Joux | Lattice reduction: A toolbox for the cryptanalyst[END_REF].

The PCG Pseudo-Random Number Generator Family

This section introduces some notations and describes the PCG64 non-cryptographic pseudorandom number generator (a.k.a. PCG-XSL-RR in the designer's terminology).

If x ∈ {0, 1} n is an n-bit string, then x[i:j] denotes the bit string x i x i+1 . . . x j-2 x j-1 , where x = x 0 . . . x n-1 (this is the "slice notation" used in Python). The set Z 2 k of integers modulo 2 k is seen as the set of k-bit strings. If x is a floating point number, then x denotes the nearest integer (using the "rounding half to even" tie-breaking rule -this is the default in IEEE754 arithmetic). If U is a vector or a sequence, then U i is the i-th element (we use capital letters for these). If U is such a sequence, we denote by U mod M the sequence (U 0 mod M, U 1 mod M, . . . ). The XOR operation is denoted ⊕, left and right rotations are denote ≪ and ≫ respectively. Modular addition is denoted + (or to make it even more explicit).

PCG64 has an internal state of 128-bit, which operate as a linear congruential generator modulo 2 128 . More precisely:

S i+1 = aS i + c mod 2 128 ,
Where the "multiplier" a is a fixed 126-bit constant. The first initial state S 0 is the seed of the generator. The increment c may be specified by the user of the PRNG to produce different output streams with the same seed (just as the IV acts in a stream cipher). If no value of c is specified, then a default increment is provided. Note that c must be odd. The default values are:

a = 47026247687942121848144207491837523525 (fixed) c = 117397592171526113268558934119004209487 (default value, user-definable)
Each time the PRNG is clocked, 64 output bits are extracted from the internal state using a non-linear function that makes use of data-dependent rotations, in a way reminiscent of the RC5 block cipher [START_REF] Ronald | The RC5 encryption algorithm[END_REF]. The six most significant bits of the internal state encode a number 0 ≤ r < 64. The two 64-bit halves of the internal state are XORed together, and this 64-bit result is rotated right by r positions.

The successive 64-bit outputs of the generator are X 0 , X 1 , . . . where:

X i = (S i [0:64] ⊕ S i [64:128] Yi ) ≫ S i [122:128] ri . ( 1
)
For the sake of convenience, we denote by Y i the XOR of the two halves of the state (before the rotation) and by r i the number of shifts of the "i-th rotation". Fig. 2 summaries the process. The overall design strategy is similar to that of a filtered LFSR: the successive states of a weak internal generator with a strong algebraic structure are "filtered" by a non-linear function.

Updating the internal state requires a 128 × 128 → 128 multiplication operation. In fact, this can be done with three 64 × 64 → 128 multiplication and two 64-bit additions. High-end desktop CPUs all implement these operations in hardware, so the generator is quite fast on these platforms.

Tools

In the rest of this paper, we often perform arithmetic operations on integers where only some bits are known. This leads to generation of unknown carries. If a, b are integers modulo 2 128 and 0 ≤ i < j < 128, then there is a carry 0 ≤ γ ≤ 1 (resp. a borrow 0 ≤ β ≤ 1) such that:

(a b)[i:j] = a[i:j] b[i:j] γ, ( 2
) (a b)[i:j] = a[i:j] b[i:j] β. ( 3

Linear Congruential Generators and Lattices

Given an integer k, a fixed multiplier a, an increment c and a "seed" x, define the sequence:

U 0 = x, U i+1 = aU i + c.
When reduced modulo 2 k , the sequence U form the successive states of a linear congruential generator (LCG). Let LCG k (x, c) denote the vector (u 0 , u 1 , u 2 , . . . ) of integers modulo 2 k . It is easy to check that:

LCG k (x + y, c + d) = LCG k (x, c) + LCG k (y, d), ( 4
) LCG k (λx, λc) = λ LCG k (x, c). ( 5
)
Let L denote the euclidean lattice spawned by the rows of the following n × n matrix:

G n,k =       1 a a 2 . . . a n-1 0 2 k 0 . . . 0 0 0 2 k . . . 0 . . . . . . . . . . . . . . . 0 0 0 . . . 2 k      
This lattice contains all n-terms geometric progressions of common ratio a modulo 2 k ; therefore the first n terms of the sequence LCG k (x, 0) give the coordinates of a vector in this lattice.

Reconstructing the state of a truncated linear congruential generator can generally be seen as the problem of finding a point of this lattice given only an approximation thereof, for instance when the least-significant bits of each components have been dropped. The "lattice approach" to truncated linear congruential generators is due to Frieze et al.

[FHK + 88].
Let U be a vector of n integers modulo 2 k such that

U i+1 = aU i mod 2 k ; it is a geometric progression of common ratio a, modulo 2 k , therefore (U 0 , . . . , U n-1 ) ∈ L. Let T i = U i [k -:k]
denote the top bits of U i , and let N denote an arbitrary "noise vector" such that N i ∈ {-1, 0, 1}. Finally, set T i = T i + N i mod 2 . We will be facing the following problem ("reconstructing noisy truncated geometric series") several times:

INPUT T = ( T 0 , . . . , T n-1 ) ∈ (Z 2 )
n , a (noisy) version of U truncated to the top bits.

OUTPUT U 0 mod 2 k , the first term of the (non-truncated) geometric sequence.

We will be facing a "high-dimension" instance in section 5.3 and many "low-dimension" in sections 4 and 5. The rest of this section discusses algorithmic tools to solve these problems. We first claim that 2 k-T is "close" to a point of the lattice L n,k .

Lemma 1. There exists U ∈ L such that U i ≡ U i mod 2 k and U -2 k-T ≤ 2 √ n2 k-.
Proof. We first observe that U belongs to the lattice L. We start by setting U ← U , and we examine all coordinates of U :

• If T i = T i + N i (without modulo
), then we have:

2 k-T i -U i = 2 k-(T i + N i ) -U i = 2 k-N i -U i [0, k -] ≤ 2 k-N i | + |U i [0, k -] ≤ 2 k-+1 .
• Otherwise, there are two possible "wraparound" cases:

-Either T i = 0 and N i = -1, which leads to T i = 2 -1. In this case, we have

U i = U i [0,
k -] and we set U i = U i + 2 k (note that this amount to adding a lattice vector to U , so U stays in the lattice). We have:

2 k-T i -(U i + 2 k ) = 2 k-+ U i [0, k -] ≤ 2 k-+1 .
-Or T i = 2 -1 and N i = +1, which leads to T i = 0. This implies that

U i = 2 k -2 k-+ U i [0, k -]; we set U i = U i -2 k (
again, with this modification U stays in the lattice), and we find:

2 k-T i -(U i -2 k ) = 2 k -2 k-+ U i [0, k -] -2 k ≤ 2 k-+1 .
In the end, we have U -2 k-T ≤ 2 √ n2 k-, and U ≡ U mod 2 k .

Reconstruction in "High" Dimension Using an Exact CVP Solver

Lemma 1 tells us that the approximation of a geometric sequence obtained by dropping least-significant bits cannot be arbitrarily far from a lattice point which reveals U 0 mod 2 k . Therefore, we may possibly reconstruct truncated geometric series by finding the lattice vector closest to the approximation we have. This means solving instances of the wellknown Closest Vector Problem (CVP), a fundamental hard problem on lattices. It is NP-hard, and all known algorithms are exponential in the dimension of the lattice, yet they can be fairly practical up to dimension ≈ 70. Let CVP(L, x) denote the vector of L closest to the input vector x. Using the same notations as above, we want to know if CVP L, 2 k-T is indeed U . This will necessarily be the case when 2 k-T -U is smaller than the length of the shortest non-zero vector of L -this quantity, the first minimum of the lattice, is denoted by λ 1 (L). By the triangular inequality, we have:

| CVP(L, 2 k-T ) -U | ≤ CVP L, 2 k-T -2 k-T + 2 k-T -U .
But, as the vector U belongs to the lattice, by definition of the closest vector and by lemma 1:

| CVP(L, 2 k-T ) -U | ≤ 2 2 k-T -U ≤ 4 √ n2 k-.
If we can prove that the right side of this inequality is smaller than the first minimum of the lattice λ 1 (L), then we would have proved that CVP(L, 2 k-T ) indeed reveals U 0 mod 2 k .

In section 5.3, we will be facing the problem of reconstructing a geometric sequence modulo 2 128 given arbitrarily many (noisy versions of the) most-significant 6 bits of successive elements of the sequence. Therefore we have k = 128 and = 6, and we wish to determine the required number of samples, i.e. the value of n. This means finding the values of n such that √ n2 124 ≤ λ 1 (L). Starting from n = 122/6 , we computed the length of the shortest vector of the lattice spanned by G n,128 for each successive n until the condition holds. The Shortest Vector Problem (SVP) is another well-known lattice NP-hard problem; we used the (almost) off-the-shelf G6K library [ADH + 19], which gave results very quickly by sieving. fplll [dt16] was too slow above dimension 50, in the default settings.

After this computation, we found that the minimal possible n is 63: with n = 63, the shortest vector of L has length greater than 2 127.02 , which is high enough. This vector can be obtained by bootstrapping the geometric sequence with U 0 = 12144252875850345479015002205241987363 then reducing the terms modulo 2 128 in zero-centered representation (subtracting 2 128 to U i if U i > 2 127 ). It follows that when n ≥ 63, k = 128 and = 6, any CVP oracle will return a vector congruent to the original U when given T .

Reconstruction in Low Dimension Using Babai's Rounding

In sections 4 and 5.1 we will need to reconstruct billions of noisy truncated geometric series modulo 2 64 with very few terms, of which a large fraction of most-significant bits are known. In this setting, the CVP problem becomes much easier. This enables us to use faster and more ad hoc methods, such as Babai's rounding algorithm [START_REF] Babai | On lovász' lattice reduction and the nearest lattice point problem[END_REF].

If M is a square matrix, we denote by ~M ~the induced matrix norm :

~M ~= sup x∈R n xM x
In the case of the • 2 norm used throughout this paper, M is the largest singular value of G; equivalently, it is the square root of the absolute value of the largest eigenvalue of M t M . Denote again by L the n-dimensional lattice spanned by the rows of G n,64 , and let H denote the LLL-reduction of G n,64 . The same lattice is also spanned by the rows of H. For instance, with n = 3:

H =   -1241281756092
3827459685972 -728312298332 -5001120657083 -2117155768935 5479732607037 8655886039732 3303731088004 6319848582548

 

Set S = 2 k-T H -1 ; as 2 k-T is not a priori an element of the lattice, S need not be an integer vector. Let then R denote the rounding of S, i.e. R i = S i . Then RH is an element of the lattice. Under the right conditions, it will be the vector of the lattice closest to 2 k-T . Indeed:

U -RH = U -R -2 k-T H -1 + 2 k-T H -1 H = U -2 k-T -R -2 k-T H -1 H ≤ U -2 k-T + R -2 k-T H -1 × ~H~.
By definition, R is the closest integer vector to 2 k-T H -1 . But as U is an element of the lattice, U H -1 is an integer vector. Thus R-

2 k-T G -1 is shorter than U G -1 -2 k-T G -1 .
Hence :

U -RH ≤ U -2 k-T + R -2 k-T H -1 × ~H≤ U -2 k-T + U H -1 -2 k-T H -1 × ~H≤ U -2 k-T + U -2 k-T × H -1 × ~H≤ U -2 k-T × 1 + ~H-1 ~× ~H~ .
Note that ~H-1 ~× ~H~is the condition number of the matrix H. Lattice reduction has the side effect of reducing the condition number, therefore it makes sense to use an LLL-reduced basis. If we can prove that the right side of the inequality is smaller than the first minimum of the lattice, then we would have proved that RH is indeed the closest vector we were searching for. Because we have fixed k = 64, by lemma 1 we have

2 k-T -U ≤ 2 √ n2 64--1.
So, if we fix n we can search for the minimum number of known most-significant bits such that: When is greater than the values given in table 1, then Babai's rounding technique will always return the closest vector, and allow us to reconstruct a truncated geometric series.

1 + ~H~× H -1 2 √ n2 64-≤ λ 1 (L)

Application to the lehmer64 generator

Adapting the previous reasoning enables an efficient state reconstruction algorithm for the lehmer64() generator shown in Fig. 1. When clocked, it outputs the top 64 bits of a geometric sequence (k = 128 and = 64). Three successive outputs are sufficient to reconstruct the internal state using Babai's rounding technique. This yields the following reconstruction algorithm: def reconstruct(X):

""" Produce the internal state of the generator given three consecutive outputs of lehmer64(). 16 multiplications, 1 division, 11 additions and 3 roundings only.

""" We first consider the easier case where the "increment" (the c term in the definition of the underlying linear congruential generator) is known -recall that a default value is specified in case the user of the pseudo-random generator does not want to provide one. In this case, reconstructing the 128-bit internal state S i of the generator is sufficient to produce the pseudo-random flow with 100% accuracy (the generator can also be clocked backwards if necessary, so that the seed can be easily reconstructed). We therefore focus on reconstructing S 0 (the seed) from X 0 , X 1 , X 2 , . . . . A very simple strategy could be the following:

a = 0xda942042e4dd58b5 r = round(2.64929081169728e-7 * X[0] + 3.51729342107376e-7 * X[1] + 3.89110109147656e-8 * X[2]) s = round(3.12752538137199e-7 * X[0] -1.00664345453760e-7 * X[1] -2.16685184476959e-7 * X[2]) t = round(3.54263598631140e-8 * X[0] -2.05535734808162e-7 * X[1] + 2.73269247090513e-7 * X[2]) u =
1. Guess the 64 upper bits of S 0 (this includes the rotation).

2. Compute the missing 64 lower bits using (1), with:

S 0 [0:64] = S 0 [64:128] ⊕ (X 0 ≪ S[122:128]). 3. Compute S 1 then extract X 1 ; if X 1 is correct, then output S 0 .
This "baseline" procedure requires 2 64 iterations of a loop that does a dozen arithmetic operations; it always output the correct value of S 0 , and may output a few other ones (they can be easily discarded by checking X 2 ). An improved "guess-and-determine" state reconstruction algorithm is possible, which essentially amounts to expose a truncated version of the underlying linear congruential generator, and attack it using the tools exposed in section 3. This is possible by combining the following ingredients:

• The underlying linear congruential generator uses a power-of-two modulus, therefore the low-order bits of S i+1 are entirely determined by the low-order bits of S i . More precisely, we have:

S i+1 = aS i + c mod 2 , for all 0 ≤ ≤ 128 (6)
Therefore, guessing the least-significant bits of S 0 yields a "long-term advantage" that holds for all subsequent states.

• Guessing a 6-bit rotation r i gives access to Y i (the XOR of the two halves of the internal state). Thus, if a part of the state is known, then this transfers existing knowledge to the other half.

In figure 3, we see that guessing S 0 [0: ] and a few 6-bit rotations r i give access to S i [58:64 + ] for the corresponding states. Therefore, looking at S i [ :64 + ], we are facing a truncated linear congruential generator on 64 bits, where we have access to the most 6 + bits of each state (denoted by T ), for a few consecutive states. This is sufficient to reconstruct entirely the successive states of this truncated linear congruential generator. This reveals S 0 [ :64 + ], and using (1) the entire S 0 can be reconstructed. The precise details follow.

We consider the sequence of internal states S = (S 0 , S 1 , . . . ) = LCG 128 (S 0 , c). We will guess the least-significant bits of S 0 , therefore let us assume that their value is known and denote it by w. We define S = LCG 128 (S 0 -w, 0) and K = LCG 128 (w, c) -this r 0 w S 0 64 bits 64 bits ??????????????????? T 0 r 1 w 1 S 1 ??????????????????? T 1 r 2 w 2 S 2 ??????????????????? T 2 bits 6 64 bits bits Figure 3: A guess-and-determine algorithm to reconstruct the first internal state S 0 . Magenta bits are guessed; cyan bits are obtained using the linear congruence relation (6) modulo 2 ; yellow bits are obtained from the output and the guessed rotations using (1).

is Known. By (4), we have S = S -K. The point is that the elements of S follow a geometric progression of common ratio a; in addition, the least significant bits of each components are equal to zero. It follows that S [ :64 + ] also follows a geometric progression of common ratio a, this time modulo 2 64 . The crux of the reconstruction algorithm is to find S [ :64 + ].

We know K i [58:64 + ] for all i, and for each guessed rotation r i we have access to

T i def = S i [58:64 + ]. We want T i def = (S i K i )[58:64 + ],
which is the truncation of S i . Thanks to (3), we know that there is an unknown vector B of borrows, whose components are either 0 or 1, such that T = S[58:64 + ] K[58:64 + ] B. Because the borrows are unknown, we in fact compute T = S[58:64 + ] K[58:64 + ], and clearly T = T B. We are thus in the context of the problem discussed in section 3, namely reconstructing a geometric sequence given 6 + (noisy) most-significant bits. The "noise" is the unknown vector B of borrows.

We will guess n rotations and least-significant bits of the state, for a total of 2 6n+ guessed bits. Table 1 gives a lower-bound on given n, and we see that the total number of guessed bits reaches a minimum of 38 when n = 3 and = 20. Therefore, success is guaranteed if we guess = 20 low-order bits of the state and three consecutive rotations.

The algorithm that reconstructs the internal state of the PCG64 generator with known increment proceeds as shown in algorithm 1. The point is that when the guesses are correct, then from the truncated geometric series T , the solution of the CVP instance reveals U j = S j [ :64 + ]. From there, the correction of the algorithm is easily established.

The procedure is completely practical. More details are given in section 6. Let us just mention that the procedure often works (twice faster) with = 19 or even four times faster with = 18 (with a reduced success probablity).

Algorithm 1 State reconstruction Algorithm (case where c is known)

1: procedure ReconstructState (X 0 , X 1 , X 2 ) 2:
// Statement involving j must be repeated for j = 0, 1, 2. for 0 ≤ w < 2 do Guess least-significant bits of S 0 6:

K j ← a j w + c(a j -1)(a -1) -1 mod 2 128 Known part 7:
for 0 ≤ r 0 , r 1 , r 2 < 64 do Guess rotations 8:

Y j ← X j ≪ r j Undo rotations 9:

T j ← (r j ⊕ Y j [58:64]) + 64 • (K j ⊕ Y j ) [0: ] Truncated LCG output 10:
T j ← T j K j [58:64 + ] Truncated geometric series on 6 + bits 11:

(U 0 , U 1 , U 2 ) ← 2 58-• ( T 0 , T 1 , T 2 ) • H -1 • H CVP (Babai rounding) 12: S 0 [0:64] ← K 0 [0:64] + 2 • U 0 [0:64 -] Reconstruct S 0 13: S 0 [64:128] ← S 0 [0:64] ⊕ Y 0 14: S 1 ← aS 0 + c Recompute X 1 15: Y 1 = S 1 [0:64] ⊕ S 1 [64:128] 16: if Y 1 = Y 1 then Check consistency 17:
output S 0 as a candidate internal state.

State Reconstruction for PCG64 With Secret Increment

The algorithm of section 4 does not apply directly to the general case where the value of c is unknown. A "baseline" procedure would consist in guessing S 0 [64:128] and S 1 [64:128]; using eq. ( 1), this would reveal S 0 and S 1 ; from there, the increment c is easily to obtain, and every secret information has been reconstructed. This would take 2 128 iterations of a very simple procedure, which is completely infeasible.

Set ∆S i = S i+1 S i ; it is easily checked that ∆S i is a geometric progression of common ratio a. Therefore, reconstructing both S 0 and ∆S 0 is sufficient to compute all subsequent states (and recover the unknown increment c). The global "guess-and-determine" strategy is essentially the same as before: gaining access to a truncated version of ∆S i , solving a small SVP instance, reconstructing ∆S 0 , then checking consistency.

Let us set:

∇S i def = S i -S 0 ≡ i-1 j=0 ∆S j ≡ ∆S 0 • i-1 j=0 a j ≡ ∆S 0 a i -1 a -1 mod 2 128 (7)
Note that ∇S 0 = 0 and ∇S 1 = ∆S 0 . Therefore, knowledge of ∆ 0 entails that of the whole sequence of ∇S i . The prediction algorithm we propose proceeds in three phases:

1. Reconstruct ∆S 0 [0:64 + ] from X 0 , . . . , X 4 , check consistency with X 5 , . . . , X 63 .

2. Reconstruct all rotations r i from this partial knowledge.

3. Fully reconstruct ∆S 0 from the rotations.

4. Reconstruct S 0 from ∆S 0 and the rotations.

Only the first phase is computationally intensive. The four steps are discussed in the next four subsections.

Partial Difference Reconstruction

In order to access to a part of ∆S i , we use the same "guess-and-determine" strategy as in section 4: we guess the least significant bits of S 0 and some rotations, then check consistency. The difference is that, since c is unknown, we must in addition guess the least significant bits of c to obtain the same "long-term advantage" (c is always odd; this makes one less bit to guess). We must also guess k + 1 successive rotation to get information on k successive differences ∆S i .

Confirming that the guesses are correct is less immediate. When c was known, we could reconstruct the internal state; from there, filtering out the bad guesses was easy. When c is unknown, the same strategy does not work, but a very strong consistency check can still be implemented.

We consider again the sequence of internal states S = (S 0 , S 1 , . . . ) = LCG 128 (S 0 , c). We will guess the least-significant bits of S 0 and of c, therefore let us assume that their value is known and denote it by w 0 and c 0 . We define S = LCG 128 (S 0 -w, c -c 0 ) and K = LCG 128 (w 0 , c 0 ) -again, K is Known and S = S -K. This time, the components of S do not follow a geometric progression; but we still have that the least significant bits of each S i are zero. Set ∆S i def = S i+1 -S i ; ∆S [ :64 + ] follows a geometric progression of common ratio a modulo 2 64 (again). This time, we have to find ∆S 0 [ :64 + ].

As in section 4, we have access to T i def = S i [58:64 + ]. We want to subtract the known part to obtain T i def = (S i K i )[58:64 + ], which is the truncation of S i . This again introduces an unknown vector B of borrows, and in fact we can only compute T = S[58:64 + ] K[58:64 + ], with T = T B. As explained above, to access a geometric sequence, we would like to obtain ∆T i def = T i+1 -T i , but we can only compute:

∆ T i def = T i+1 -T i = (T i+1 T i ) (B i+1 B i )
We are thus still in the context of the problem discussed in section 3, but this time the "noise" caused by the carries is given by B i+1 -B i . When the guesses are correct, then Babai's rounding will reconstruct ∆ S [ :64 + ] from ∆ T . This in turn yields ∆S 0 [0:64 + ].

Once we have found ∆S 0 [0:64 + ], we can compute ∇S i [0:64 + ] for any i because eq. ( 7) holds modulo 2 64+ ; because we have guessed the first rotation and the least significant bits of the state, using (1) we gain access to S 0 [58:64 + ]; combined with the "differences" ∇S i , this reveals S i [58:64 + ] for any i (and we already had S i [0: ]). This allows us to compute Y i [0: ] = S i [0: ] ⊕ S i [64:64 + ] for any i. Given a "fresh" output X i , and assuming that the guesses are correct, then we should have:

S i [0: ] ⊕ S i [64:64 + ] = (X i ≪ r i )[0: ]. (8)
In particular, if the guesses were correct, then we should have for any i:

S i [0: ] ⊕ S i [64:64 + ] ∈ (X i ≪ r)[0: ] | 0 ≤ r < 64 . ( 9
)
If none of the 64 possible rotations yields a match, then the guesses made beforehand have to be wrong. As a consequence, bad guesses can be filtered with an arbitrarily low probability of false positives, by trying several indices i.

A few details still need to be fleshed out. To be precise, let us assume that we have guessed the least-significant bits of S 0 (we denote them by w 0 ) and the first rotation r 0 . Set Y 0 = X 0 ≪ r 0 . We obtain the i-th state by S i ≡ ∇S i S 0 ; however, because the "middle" of S 0 is unknown, then an unknown carry may cross the 64-th bit during the addition and perturb S i [64:64 + ]. As a result, there is an unknown vector C, whose components are either 0 or 1, such that such that:

S i [64:64 + ] = C i ∇S i [64:64 + ] (w 0 ⊕ Y 0 [0: ] S0[64:64+ ]
) Algorithm 2 Partial difference reconstruction algorithm (when c is unknown).

1: procedure ConsistencyCheck(∆S 0 , w 0 , Y 0 , X 5 , . . . , X k ) 2: v 0 = w 0 ⊕ Y 0 [0: ] v 0 = S 0 [64:64 + ]
3:

for i = 5, . . . , k do 4:

u i ← ∆S 0 (a i -1)(a -1) -1 mod 2 64+ u i = ∇S i [0:64 + ]
5:

w i = w 0 u i [0: ] w i = S i [0: ] 6: v i = v 0 u i [64:64 + ] S i [64:64 + ] ∈ {v i , v i } 7: v i = v i 1 8: C i ← {w i ⊕ (X i ≪ r i )[0: ] | 0 ≤ r i < 64} Check eq. (9) 9: if {v i , v i } ∩ C i = ∅ then 10:
return False Bad Guesses 11:

return True No inconsistency 12:

13: procedure ReconstructPartialDifference(X 0 , . . . , X k )

// Statement involving j must be repeated for j = 0, 1,

H ← LLL reduction of G 4,64 16:

← 14 17:

for 0 ≤ w 0 < 2 and 0 ≤ c 0 < 2 -1 do Guess least-significant bits 18:

K j ← a j w 0 + (2c 0 + 1)(a j -1)(a -1) -1 mod 2 128 Known part 19:

for 0 ≤ r 0 , r 1 , r 2 , r 3 , r 4 < 64 do Guess rotations 20:

Y j ← X j ≪ r j Undo rotations 21:

T j ← (r j ⊕ Y j [58:64]) + 64 • (K j ⊕ Y j ) [0: ] Truncated LCG 22:
T j ← T j K i [58:64 + ] Cancel known part 23:

∆ T j = T j+1 T j Difference (truncated geom. seq.)

24:

(∆U 0 , . . . , ∆U 3 ) ← (∆ T 0 , . . . ,

∆ T 3 ) • 2 58-• H -1 • H CVP 25: ∆S 0 [0:64 + ] ← (K 1 K 0 ) [0: ] + 2 • ∆U 0 [0:64] Check 26:
if ConsistencyCheck(∆ 0 , w 0 , Y 0 , X 5 , . . . , X k ) then 27:

return (w 0 , c 0 , r 0 , . . . , r 4 , ∆S 0 ).

In algorithm 2, ConsistencyCheck uses eq. ( 9) combined with this observation to discard bad guesses.

The heart of the algorithm is again the reconstruction of a truncated geometric progression. Looking at table 1, we see that the best choice consists in guessing 5 consecutive rotations and = 14 least-significant bits. Therefore, ReconstructPartialDifference does 2 57 iterations of the inner loop, and succeeds deterministically.

Predicting all the Rotations

Knowing the values of ∆S 0 [0:64 + ] as well as the least-significant bits of S 0 and c is sufficient to get rid the nastier feature PCG64: armed with this knowledge, we can determine all the subsequent rotations deterministically, at negligible cost, using eq (8). For each index i, it suffices to try the 64 possible values of r i ; only one should satisfy eq (8). The complete pseudo-code is shown in algorithm 3.

It is unlikely that several possible values of r i match: each value is "checked" on bits, so an accidental match happens with probability 2 -6 . The total number of lists returned by ReconstructRotations then follows a binomial distribution of parameters 2 -6 , k. With = 14 and k = 64, then only one rotation vector should pass the test for 0 ≤ i < 64 on average. T ← ReconstructRotations(∆S 0 , v 0 , i + 1, k) Find all the (r i+1 , . . . , r k )

Algorithm 3 Rotations and full difference reconstruction algorithm

6:

H ← [] List of possible r i 's 7:

u i ← ∆S 0 (a i -1)(a -1) -1 mod 2 64+ u i = ∇S i [0:64 + ]
8:

w i = w 0 + u i [0: ] mod 2 w i = S i [0: ] 9: v i = v 0 + u

Full Difference Reconstruction

Using X 0 , X 1 , . . . , X 63 , we recover all rotations and thus we recover the 6 most-significant bits of S 0 , S 1 , . . . , S 63 . This allows us to compute the 6 most significant bits of the differences ∆S i between consecutive states (up to missing carries), and we are faced with the problem of reconstructing a 128-bit geometric progression using 63 consecutive outputs truncated to their 6 most-significant bits. There is again an unknown vector of borrows B such that ∆S i [122:128] C i = r i+1 r i .

Reconstructing ∆S 0 from the r i is exactly the problem discussed in section 3.2. This can be done by solving an instance of CVP in dimension 63. We use the off-the-shelf CVP solver embedded in fplll: it runs in no measurable time.

Complete State Reconstruction

Once all the rotations have be recovered and ∆S 0 has been found entirely, the only thing that remain is to actually find the entire S 0 . For this, we use again eq. (1), coupled with the "differences":

S i = S 0 ∇S i Y i = S i [0:64] ⊕ S i [64:128].
The Y i and ∇S i are known, ∇S 0 = 0, and the problem consists in recovering S 0 . We could probably encode it as an instance of SAT, feed it to a SAT-solver and be done with it.

Nevertheless, here is a detailed recovery procedure which obtain all bits of S 0 , from right to left, by exploiting the non-linearity of modular addition. It takes negligible time. Let C i the vector of (incoming) carries generated during the addition of S 0 and ∇S i :

S i [j] = S 0 [j] ⊕ ∇S i [j] ⊕ C i [j] C i [j] = 0 if j = 0 MAJ(S 0 [j -1], ∇S i [j -1], C i [j -1]) if j > 0
Combining all the above, we have:

Y i [j] = Y 0 [j] ⊕ ∇S i [j] ⊕ ∇S i [64 + j] ⊕ C i [j] ⊕ C i [64 + j] ( 10
)
This useful equation enables an induction process.

• When j = 0, the 0-th carries are zero, and therefore eq. ( 10) reveals the 64-th carries:

C i [64 + j] = Y 0 [j] ⊕ Y i [j] ⊕ ∇S i [j] ⊕ ∇S i [64 + j] .
• Next, suppose that C i [0:j], S 0 [0:j -1], C i [64:64 + j] and S 0 [64:64 + j -1] are known, for all i. We can compute C i [j] ⊕ C i [64 + j] for any i using eq. (10). We then look a a specific index i > 0 such that

∇S i [j -1] = C i [j -1] and ∇S i [64 + j -1] = C i [64 + j -1].
The point is that, thanks to the majority functions,

C i [j] = S 0 [j -1] and C i [64 + j] = ∇S i [64 + j -1]. It follows that: S 0 [j -1] = Y 0 [j -1] ⊕ Y i [j -1] ⊕ ∇S i [j -1] ⊕ ∇S i [64 + j -1] ⊕ ∇S i [64 + j -1]
From there, we also have

S 0 [64 + j -1] = Y 0 [64 + j -1] ⊕ S 0 [j -1]
, and the j-th carry bits can be computed normally.

The whole procedure is shown in algorithm 4. Note that once S 0 has been found, then all subsequent states can be computed with error using S i = S 0 ∇S i . In particular, computing S 1 gives c by c ← S 1 aS 0 . This complete the reconstruction procedure for PCG64.

Implementation and Practical Results

We have implemented the state reconstruction algorithms described above using a mixture of C (for the computationally expensive parts) and Python (for the rest). We used the fplll library [dt16] to solve CVP instances exactly in dimension 63.

In this section, we briefly outline important aspects of our implementations and present practical results. Our codes are available in the supplementary material as well as online at:

Algorithm 4 Full state reconstruction algorithm 1: function ReconstructState(∆S 0 , r 0 , . . . , r k , X 0 , . . . , X k ) 2:

for i = 0, 1, . . . , k do Setup 3:

∇S i ← ∆S 0 (a i -1)(a -1) -1 mod 2 128 4: Y i ← X i ≪ r i Undo rotations 5: C i [0] ← 0 Bootstrap induction 6: C i [64] ← Y i [0] ⊕ Y i [j] ⊕ ∇S i [j] ⊕ ∇S i [64 + j] 7:
for j = 1, 2, . . . , 64 do Induction 8:

i ← ⊥ Find good index 9:

for k = 1, 2, . . . , k do 10: if ∇S k [j -1] = C k [j -1] ∧ ∇S k [64 + j -1] = C k [64 + j -1] then 11: i ← k 12: if i = ⊥ then
No suitable indice found?

13:

Abort with Failure 14:

Compute next state bit 15:

S 0 [j -1] ← Y 0 [j -1] ⊕ Y i [j -1] ⊕ ∇S i [j -1] ⊕ ∇S i [64 + j -1] ⊕ ∇S i [64 + j -1] 16: S 0 [64 + j -1] ← Y 0 [64 + j -1] ⊕ S 0 [j -1] 17:
for i = 0, 1, . . . , k do Compute next carries 18: The designer of PCG was kind enough to send us two sets challenge inputs: one with the default (known) increment and one with a random secret increment. She generated random seeds and provided us with the first outputs of the pseudo-random generator. We were able to reconstruct the seed with an extremely high confidence level, because they re-generate the same outputs. We emailed back the seeds and received confirmation that they were indeed correct.

C i [j] ← MAJ(S 0 [j -1], ∇S i [j -1], C i [j -1]) 19: C i [64 + j] ← MAJ(S 0 [64 + j -1], ∇S i [64 + j -1], C i [64 + j -1]) 20: return S 0
We have therefore successfully taken the challenge of predicting the output of the PCG64 generator.

The analysis of section 3 yields parameters that guarantee that the reconstruction procedure always succeeds. In most cases, these parameters are pessimistic. We ran a series of experiments to determine more practical choices: using smaller-than-guaranteed values of (the number of guessed least-significant bits), we measured the success probability of the state reconstruction procedure. The results are shown in table 2.

Known Increment

When the increment c is known, algorithm 1 is all it takes to reconstruct the internal state of the generator and predict it (or output the seed). We implemented it in C, using return magic.l; } This hack exploits the IEEE754 representation of double-precision floats: the mantissa lies in bits [0:52] while the sign bit and the exponents take the 12 most significant bits. Adding 2 52 + 2 51 forces the mantissa to shift to the correct position and inserts an extra 1 bit at position 51. The two shifts clear the extra bit and the exponent, while correctly expanding the sign bit.

Unknown Increment

When the increment c is known, the internal state of PCG64 can be practically reconstructed from X 0 , . . . , X 63 using the algorithms shown in section 5. Only algorithm 2 is computationally expensive; we implemented it in C, while we implemented algorithms 3 and 4 in Python.

We have shown that algorithm 2 is correct when = 14. The procedure does 2 29+2 iterations of the inner loop, so decreasing would really be interesting. Looking at table 2, we settle for = 13 in the worst case; let T denotes the running time when = 13.

It seems that the most promising strategy consists in choosing = 11; if the reconstruction procedure fails, then we try again with different inputs. The expected running time of this approach number of trials is T /(16 × 0.64) ≈ T /10.25. In our implementation, T = 200, 000 CPU hours, so the expected running time of the reconstruction procedure is about 20, 000 CPU hours. In fact we were lucky: on the challenge input, the first attempt with = 11 succeeded, so the whole process took only 12, 500 CPU hours.

It actually ran in 35 wall-clock minutes using 512 cluster nodes, each equipped with two 20-cores Intel Xeon Gold 6248 @ 2.5Ghz ("Cascade Lake"). The actual machine is the jean-zay computer located at the IDRIS national computation center. Note that on this particular parallel computer, running the algorithm with = 13 would take 10 hours using the same amount of resources, so the whole procedure is practical, even in the absolute worst case.

The outer loop of algorithm 2 makes 2 2 -1 iterations while the inner loop makes 2 30 iterations. Using a single hardware execution context, we measured that one of the outer loop takes between 41.5s and 44s (apparently not all nodes of the cluster are running at exactly the same speed, potentially because of "turbo boost" adjustments and thermal constraints). Because of this variability, we implemented a master-slave work distribution pattern, in which a master process dispatches iterations of the outer loop to slave processes. This also made checkpointing very easy. We used MPI for inter-process communication.

With = 11, the whole process took 2 56.74 CPU cycles, which makes less than 6 cycles per iteration of the inner loop. We used essentially the same implementation tricks discussed above. However, this time we had to additionally implement the ConsistencyCheck procedure, which is called in the inner loop. We observed that the set of possible candidate values C only depends on w 0 (the variable of the outer loop). Therefore, before entering the inner loop, we precompute a bit field of size 2 describing C i . To simplify the implementation, we flatten them by computing C = ∪ i C i . This slightly increase the probability of false positives, but makes our code slightly simpler.

Conclusion

We have presented a practical state reconstruction algorithm for the PCG64 conventional pseudo-random number generator, the default in the NumPy library. In the worst case, we recovers all the secret information using 512 consecutive output bytes, using 2.3 CPU years of computation. We have executed the algorithm in practice using a large parallel computer. The PCG64 generator is fast and not intended for cryptographic purposes; we have shown that, in practice, this comes at the price of strong pseudo-randomness. It should absolutely not be used when unpredictability of the random numbers is required, for fear of practical attacks.

On the other hand, our results do not mean that PCG64 should be deprecated for scientific computing. But they do mean that its output has detectable properties. Whether these properties may affect the results of Monte-Carlo numerical simulations is another matter entirely.

Figure 1 :

 1 Figure 1: Some conventional pseudo-random generators, designed for speed and simplicity.

1:

  function ReconstructRotations(∆S 0 , v 0 , i, k) 2:// Return a list of potential [r i , r i+1 , . . . , r k ]; assume that v 0 = S 0 [64:64 + ]

Table 1

 1

		: minimal needed for a given n
	n ~H~× H -1	λ 1 (L)	minimum
	3	2.87	4.09e 12 2 41.9	26
	4	2.06	2.44e 14 2 47.8	20
	5	3.77	1.72e 15 2 50.6	18
	6	2.69	1.03e 16 2 53.2	15

Table 2 :

 2 Empirical success probabilities with smaller parameters.

	n = 3 (section 4)	n = 4 (section 5.1)
		Success proba.		Success proba.
	16	≈ 0.125	10	≈ 0.12
	17	≈ 0.25	11	≈ 0.64
	18	≈ 0.5	12	≈ 0.995
	19	1	13	1
	20	1 (proved)	14	1 (proved)
		https://github.com/cbouilla/pcg/

Acknowledgements. The authors thank Thomas Espitau, Damien Vergnaud and Vincent Zucca for helpful discussions. We also thank Léo Ducas for helping us with the G6K library and Melissa O'Neil for kindly indulging our request for challenge inputs.

This work was granted access to the HPC resources of IDRIS under the allocation 2019-A0060610749 made by GENCI.

Predicting the PCG Pseudo-Random Number Generator In Practice OpenMP to parallelize the outer loop that guess the least-significant bits of the state. This yields a simple multi-core implementation. We used the gcc 8.3.0 compiler.

From section 3.3, we know that guessing = 20 least-significant bits ensures deterministic success. However, we observed empirically that = 19 works with probability ≈ 1, and runs twice as fast. = 18 and = 17 run with probability ≈ 1/2 and ≈ 1/4 respectively, therefore are much less useful. In practice, we used = 19.

We ran it on a server equipped with two 16-core Intel Xeon Gold 6130 CPU @ 2.10GHz ("Skylake") CPUs. The inner loop does 2 37 iterations and terminates in 42.3s, which makes 23 core minutes.

These processors operate at a different frequency depending on the number of cores used and the type of instructions executed. Our code uses only scalar instructions, so the CPUs runs at the highest frequency tier when executing it. Using a single software thread per physical core (each core presents two hardware execution contexts, commercially called HyperThreads) allows the CPU to run at 2.8Ghz, the maximum "Turbo" frequency on all cores. Using one software thread per hardware thread reduces the frequency to ≈ 2.6Ghz, but allows to better saturate the execution units of the CPU and yields a nearly 20% speedup overall.

Therefore the algorithm requires 2 41.67 CPU cycles in total; this makes less than 5 cycles per iteration of the inner loop. We used several implementation tricks to reach this level of efficiency:

• We used the __uint128_t type provided by most C compilers to do 128-bit arithmetic when computing S 1 from S 0 . Apart from that, the algorithm has been designed to do mostly 64-bit arithmetic, for the sake of efficiency.

• Looking at the algorithm, it is clear that U 1 and U 2 are actually not needed, so we just don't compute them.

• T J is a function of w, j and r J (with j = 0, 1, 2). therefore, for each new value of w, we precompute once and for all an array indexed by (J, R J ) of the 192 possible values of T J .

• Pushing the same idea a bit further, we precompute parts of the matrix-vector product inside the rounding: this computes a linear combination of the rows of G -1 3 , in which T j is the coefficient of the j-th row. So we precompute the 576 possible products

• We enumerate the possible rotations in lexicographic order. This means that T 0 changes in each iteration while T 1 (resp T 2 ) changes every 64 (resp 4096) iterations. Therefore, in 98% of the iterations, two-thirds of the matrix-vector product inside the rounding are the same as from the previous iteration. Therefore, we fully compute the matrix-vector product only when r 1 changes and only update it when r 0 changes.

• The rounding operation, when done naively by writing llround(x), is actually a bottleneck: it calls a library function that accounted for about 20% of the total running time. We instead used the following technique, which correctly returns x whenever |x| < 2 51 : long long fast_round(double x) { union { double d; long long l; } magic; magic.d = x + 6755399441055744.0; magic.l <<= 13; magic.l >>= 13;