Michel Tsamados

Anne Tanguy

Fabien Léonforte

Jean-Louis Barrat

On the study of local stress rearrangements during quasistatic plastic shear of a model glass: do local stress components contain enough information?

Keywords:

We present a numerical study of the mechanical response of a 2D Lennard-Jones amorphous solid under steady quasistatic and athermal shear. We focus here on the evolution of local stress components. While the local stress is usually taken as an order parameter in the description of the rheological behaviour of complex fluids, and for plasticity in glasses, we show here that the knowledge of local stresses is not sufficient for a complete description of the plastic behaviour of our system. The distribution of local stresses can be approximately described as resulting from the sum of localized quadrupolar events with an exponential distribution of amplitudes. However, we show that the position of the center of the quadrupoles is not related to any special evolution of the local stress, but must be described by another variable.

Introduction

While the mechanical plastic response of crystalline solids is now well described in terms of elementary defects called dislocations, the plastic response of amorphous solids is not so well understood, and has been the object of numerous studies in the past two decades. Most theoretical works [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19], on the forced dynamics of glassy materials assume that the irreversible evolution of glasses occurs through localized, irreversible events. Such events have been well identified in simulations of 2D models [1] and are referred to as shear transformation zones (STZ) or quadrupolar events [18,19,20,21] and are somewhat similar to T1 events in the rheology of 2D foams [22].

These localized and dissipative events should correspond to the "soft mode" of the material [18,19,23,24]), which is the first mode to become unstable, when the system reaches a limit of mechanical instability. However, the local and intermittent nature of the plastic events complicates the identifica-tion of the corresponding order and control parameters. Experimental studies on disordered materials, far below the glass transition temperature [22,25,26,27,28,29,30,31,32], support the existence of a collective behaviour of localized rearrangements leading to a strongly heterogeneous mechanical response as shown in the existence of shear bands in the macroscopic plasticity of such systems. It has also been shown experimentally [32,33,34] as well as numerically [1,18,21,35,36,37] while studying the stress-strain response of amorphous materials that, in the so called plastic regime, these events lead to an intermittent release of the macroscopic stress in the material. The understanding of the role played by local stresses in the macroscopic mechanical behaviour of such systems is thus very important. Since these systems are known to have a wide distribution of local stresses, mesoscopic models [7,10,11,38] propose precisely to relate the existence of shear bands to the local rearrangements due to local thresholds and activated dynamics in the local stress itself. In these models, the collective effects would thus be due to long-range elastic couplings in solid systems, and the irreversible behaviour would be associated with local stress thresholds, and avalanche dynamics.

In the present paper, we provide an extensive study of the local stress dynamics in the quasi-static plastic response of a model glass at very low temperature, and in the plastic flow regime. The system we have studied is a 2D Lennard-Jones glass, already studied previously for understanding the heterogeneous mechanical response in the elastic regime, far below plastic flow [39,40,41,42,21]. In the first part of this paper, we discuss the evolution of the local stresses (local pressure, deviatoric stress) in the plastic flow regime, within the framework of dynamical heterogeneities [43]. We show the existence of a finite size of cooperativity appearing in the plastic flow regime, even in the limit of zero temperature and quasi-static shear. In the second part of this paper, we provide a detailed analysis of the full statistical properties of the stress components (pressure, shear stress). We show that the temporal correlations shown in the local particle rearrangements for small imposed strain [21,19] are absent in the average evolution of local stress components. We then show in the third part of this paper, that the statistical evolution of stress components can be described as resulting simply from a sum of localized quadrupolar rearrangements. This description allows us to propose a simple equation of evolution for the local stress depending unfortunately on a local yet unknown criterion for instability. In the last part, we show that this local criterion for instability is not related to any stress threshold but could be related to a different variable, such as the local stress increments.

Dynamical heterogeneity 2.1 Two points correlation function

The system we study numerically has already been described extensively in previous articles [40,41,42,21]). It is a 2D glass obtained by quenching very quickly a liquid sample obtained by molecular dynamics simulations, and made of polydisperse particles interacting through a Lennard-Jones potential. In the following, all quantities will be expressed in Lennard-Jones units (u LJ ), that is an energy ǫ LJ = 1(≈ 0.1eV ), a distance σ LJ = 1(≈ 2 Å) and mass m = 1(≈ 10g.mol -1 ). The number of particles per unit surface is ρ = 0.925, corresponding to very small pressure P ≈ 0.2 at rest. The system is made of at least N = 10000 particles, corresponding to a square with a linear size of 104 interatomic distances. Larger systems have also been considered, especially in this paper a system of N = 23125 particles, corresponding t a rectangular box of size 50 × 500. These large system sizes are necessary in order to avoid strong finite size effects. The quenching starts from the liquid state, at a temperature T = 2. Then the system is relaxed using MD to a given temperature of 1 -0.5 -0.1 -0.05 -0.01 -0.005 -0.001, where it is aged during 1000 unit times. The final step is a relaxation at zero temperature, using the Conjugate Gradient Method in order to get an instantaneous minimization of the potential energy. The system being at zero temperature, it is sheared quasi-statically, by imposing successive steps of global shear strain δǫ = 5.10 -5 on the walls. After each step, the system is relaxed into its nearest equilibrium position, thus ensuring that the applied mechanical deformation is quasi-static. The local stress components are computed on each particle i by using the usual Irving-Kirkwood formula ( [44])

σ αβ,i ≡ -1/V i .Σ j t ij d α .d β .r ij
V i being the volume of the Voronoi cell, t ij the amplitude of the interacting force between particles i and j, r ij the distance between particles and -→ d the unit vector of the bond. The macroscopic stresses are obtained by averaging the local stress components over the whole system, or alternatively by measuring the force per unit length acting on the walls along appropriate directions. In Fig. 1, we show the (averaged) macroscopic shear stress as a function of the total strain applied to the system. The global shear stress has an intermittent behaviour, with an alternance of small and large jumps, giving rise to negative slopes that are the signature of dissipative events. Note that the distribution (figure 2) of these incremental stresses is peaked around the value given by the elastic response with an elastic shear modulus µ ≈ 11.7 ( [40]) and shows around this value softer, as well as more rigid, steps.

In order to characterize the evolution of the stress in the flowing regime, we first compute the two points correlation function of the local stress components acting on each particle inside the sample. The correlation function is "Lagrangian" in the sense that the stress is followed on a given particle as the system is sheared. In the plastic flow regime, we expect the system reaches a stationary steady state with time-translational invariant (TTI) observables and correlations. Note that in quasistatic simulation, time does not appear explicitly and must be interpreted in the following as the number of steps n of elementary shear strain δǫ applied to the system. The total imposed strain is thus ǫ = n.δǫ. In order to study the stress fluctuations correlation function, we introduce the notation ρ( -→ r , n) ≡ σ( -→ r , n)-σ(n) where A denotes a spatial average and -→ r is the position of the particle; σ is a component of the stress tensor, namely σ ≡ σ dev where σ dev ≡ σ 1 -σ 2 , σ 1 and σ 2 being the eigenvalues of the local 2x2 stress tensor; or σ ≡ σ p where σ p = -(σ 1 + σ 2 )/2 is the local pressure; or σ ≡ σ xy is the local shear stress. The two-time autocorrelation function of the stress fluctuations writes

C(∆n, n) = ρ( -→ r , n + ∆n)ρ( -→ r , n), (1)
and

C(∆n) = C(∆n, n) = ρ( -→ r , n + ∆n)ρ( -→ r , n) ≈ C(∆n, n), (2)
where the last equality assumes TTI. Numerically we have computed the stress components (as mentioned before from the Irving-Kirkwood definition, denoted by σ i (n)) on each particle of the system. The discrete version of the above equation writes

C(∆n) = 1 N i ρ i (n + ∆n)ρ i (n) . (3)
where

ρ i (n) = σ i (n) -σ(n).
We have computed this quantity in the plastic flow regime (i.e. for external shear deformation ǫ greater than 1.5%), where the average total stress saturates and fluctuates around a given average value σ Y . We checked that the correlations are invariant under time translation and in Fig. 3 we plot the autocorrelation function for the shear component of the stress, on a lin-log scale and normalized to 1 for ∆n = 0. The behaviour can be decomposed into two successive relaxations. It starts with a first exponential decay at small deformation with a characteristic strain of 3.1%, this decrease to a plateau is analogous to the β -relaxation in supercooled liquids, and is followed by an apparent stretched exponential decay. A similar Kohlrausch-Williams-Watts [45,46] is usually observed in supercooled liquids and ageing glassy systems, with a correlation that evolves as C(t) = A exp(-(t/τ ) β ) with β ≤ 1. This behaviour is approximately recovered in our systems at large deformations (see inset of Fig. 3), with the exponent β = 0.7. For small strains, the relaxation is strictly exponential. Note that the large strains stretched exponential decay could also be seen as the best fit for an intermediate regime between two successive exponential decays (see the fit by an exponential decay with characteristic strain ǫ ≈ 10% at very large strains in Fig. 3). In any case, the presence of dissipative events leads to the relaxation of local stress at large deformations, above a few percents of strain.

Four points correlation function or cooperativity number

In the recent years, the quest for a dynamical cooperativity length associated with the slowing down of the dynamics of supercooled liquids [43,47,48] has led to the development of new statistical tools to characterize such "dynamical heterogeneities". In our previous paper [21] we showed that the motion of the particules in our system under shear is highly non trivial. It shows a background of heterogeneous motion even in the elastic regime at very low temperature; for larger applied strains, in the plastic regime, it shows zones of very high mobility located in the vicinity of elementary shear bands and of localized quadrupolar events, and similarly zones of low mobility far from these irreversible events (the local displacement field can vary by many orders of magnitudes). This type of behaviour in the mechanical response of glass, even at zero temperature, can be seen as some kind of dynamical heterogeneity (while the "dynamics" here is overdamped) and it is of great interest to quantify a cooperativity degree or a cooperative length scale. To do so, Berthier et al. [43] have proposed to look at the so called χ 4 four point correlation function. An other possible observable is the cooperativity number introduced by Doliwa and Heuer [47], that estimates the spatial fluctuations of the two-point correlation function introduced previously. Following reference [47] the definition we use here is

N coop ≡ V ar[ΣX i ] ΣV ar[X i ] = { i X i - i X i } 2 i { X 2 i -X i 2 } , (4)
where

X i (t ′ , t) = (σ i (t+t ′ )-σ i (t + t ′ ))(σ i (t)-σ i (t))
is a dynamical quantity associated with every particle i, and computed on various components σ of the stress tensor. In Eq.( 4), one shows easily that in the case of totally uncorrelated variables X i ( X i X j = 0), one gets

N coop = 1, but in the opposite if the cor- relation is complete (X i = X j ) one gets N coop = N .
Finally in the case of L independent groups each comprising M identical variables and of zero mean (N = LM ) one gets

N coop = 1 + i =j X i X j X 2 i = 1 + 1 X 2 i N i=1 M-1 j=1 X 2 i = M. ( 5
)
The cooperativity number is obtained by averaging over the origins t and remains a function of time t ′ (N coop (t ′ )). Note that this function does not contain any useful information on the heterogeneous motion in the linear elastic regime, since in this case, the amplitude of the stress on each particle increases linearly with time. The cooperativity number is thus constant in this regime. A maximum can be reached only in the non-linear regime, like in the plastic flow regime here, or when thermal activation starts playing a role. In Fig. 4 we have plotted this function computed on all particles located at a distance of a few (typically 10) particle diameters from the boundaries, and for different components of the stress tensor. All the stress components have approximately the same contribution to the cooperativity number, with a maximum of a few tens at an imposed strain of a few percents (≈ 4%). As described before, the maximum is a measure of the maximum number of particles on which the local stress evolves in a correlated manner. In our case, it corresponds to an ensemble of 6 × 6 to 7 × 7 particles where the stress components evolve in a cooperative way. The strain at which this maximum is reached is about 4 percents, that is beyond the plastic threshold, and close to the characteristic strain obtained from the two points correlation function C. This picture corresponds to the description of localized plastic events, that we have made already in a preceding paper [21]. In the elastic regime (as in elastic branches present in the plastic flow plateau [21]), the cooperative evolution of stress can be more important, but does not contribute significantly to the cooperativity number. Fig. 4 shows that there are strong boundary effects measurable on this quantity. Indeed, the cooperativity number is much higher in the vicinity of the walls (≈ 100) corresponding to a layer of one particle diameter along the entire system size (L x = 104). Finally, note that this maximum cooperativity number is finite in the quasi-static limit studied here, suggesting a saturation at very low shear rate.

Statistical analysis of the local shear stress

In the previous section we have studied the average evolution of the stress components by looking at averaged quantities like the number of cooperativity, or the two-point correlation function. We will now examine in details the full histogram of stress changes. The results are shown here for the shear stress component, but we have checked that all the following results are also valid for other components (pressure and deviatoric stress). In figure 5 we show the typical evolution of the shear stress component on a given particle up to a total strain of 25%1 . The variations of stresses are larger than expected if they would be Brownian. In order to analyze this evolution, we have plotted in figure 6 the distribution P (∆σ xy , ∆n) of the shear stress increments

∆σ xy (∆n, n) ≡ ( σ xy (n + ∆n) -σ xy (n + ∆n)) -(σ xy (n) -σ xy (n)) (6)
for various numbers ∆n of incremental shear steps (∆n = 1, 2, 4, 8, 16...), averaged over the origins n and over the entire system. The distributions P (∆σ xy , ∆n) are all symmetric, with zero mean. For ∆n → ∞ one recovers a gaussian distribution, which is consistent with the central limit theorem. However, at small imposed strains, these distributions are not gaussian as would be the case for a Brownian evolution. A finer analysis of the distribution of the elementary increments P (∆σ xy , ∆n = 1) (see fig 6b andc) shows three zones. At small incremental stress jumps ∆σ xy , we can see a plateau of approximately constant probability, whose width evolves inversely proportional to the volume V of the sample, followed by an apparent power-law decay, in a zone of approximately three decades for 100/V ≤ ∆σ xy ≤ 10 where P (∆σ xy , ∆n = 1) ∝ 1/∆σ α+1 xy , and concluded by an exponential cut-off independent on the system size (characteristic shear stress ∆σ xy ≈ 1.4). This upper cut-off allows for a finite variance of the local stress evolution. In the absence of any temporal correlations, the entire process can thus be described by the central limit theorem, with a scale invariant distribution on a finite interval. In this intermediate stress range, the probability density function (hereafter, pdf) P (∆σ xy , ∆n) is well reproduced by the scale invariant relation

P (∆σ xy , ∆n) = ∆n -H f ( |∆σ xy | ∆n H ) or P (∆σ xy , ∆ǫ) = ∆ǫ -H f ( |∆σ xy | ∆ǫ H ) (7)
with

f (u) ∝ u 0 , for u ≪ 1 u -α-1 , for 1 ≪ u ≪ cste.∆n -H (8)
In the intermediate stress range, the process can thus be considered as self-similar. Fig 7(a) illustrates this scaling with a good superposition of the distributions for α = 0.7 and H ≃ 1 α (only for not to large ∆n since for very large ∆n the upper cut-off discussed before contributes significantly to the resulting distribution). The exponent α describes the algebraic (slow) asymptotical decay of the distribution of the incremental jumps P (∆σ = s, ∆n = 1) ∝ s -α-1 (as shown in fig 6(a)). The exponent H is related to the evolution of the stress jumps as a function of the applied shear strain ∆n. It characterizes the ∆n-dependence of the crossover between a regime of approximately uniform probability (for |∆σ| ≪ ∆n H ) and the power law regime (for |∆σ| ≫ ∆n H ). The coefficient H also accounts for possible temporal statistical correlation between jumps. As stated in Taqqu et al. [START_REF] Samorodnitsky | Stable Non-Gaussian Random Processes: Stochastic Models with Infinite Variance[END_REF] the only non degenerate α-stable self similar processes with stationary increments that verify H = 1 α and where 0 < α < 1 are the α -stable Lévy motions. As described above, the evolution of the shear stress is thus of Lévy flight type, but only in the intermediate stress (and applied strain) range. The Lévy flight evolution implies by definition, first that there is no temporal correlation between local stress jumps during the shear of the sample, second that the variance of elementary changes is infinite (as long as the exponentiel cut-off is neglected) 2 . These results can be compared with the study of the pdf P (∆y, ∆n = 1) of the positional jumps in the transverse direction, that showed in contrary that these jumps were correlated in time for small imposed deformation [21].

According to the above, a plausible equation that would describe the evolution within a Lévy 2 In the case of Brownian motion the properties stated above on the pdf and on the process are verified (self similarity, α-stable process, stationary increments) but with a finite variance (α > 2) allowing to a unique value

H = 1/2.
flight process for the stress component σ, averaged over the whole system, is

d(σ(i, ǫ) -σ(ǫ)) dǫ = η(i, ǫ) (9)
where η(i, ǫ) is a stochastic process whose spatial average η(ǫ) is the process that is entirely characterized above (ie by the distribution of the elementary increment and by the absence of correlation between successive increments), and ǫ is the external imposed shear strain. For small imposed shear strain, it can be mentioned that such an equation with a noise corresponding to a Lévy motion cannot be reinterpreted in terms of the usual Fokker-Planck equation as the moments of the shear stress σ are non vanishing for all orders. However, for large strain intervals ∆ǫ, the stochastic process becomes gaussian, due to the existence of the upper cut-off in the distribution of η.

We will now propose a microscopic origin for the elementary process η, by comparing it with a spatial redistribution of quadrupolar stress.

But before, it can be interesting (as in Ref. [19,20]) to identify the contribution to the distribution P (∆σ, ∆n) of events with a release of the macroscopic stress (plastic events with ∆σ macro < 0). This is shown in figure 7(b). We see here that for large imposed strains (∆n > 64 it means ∆ǫ > 0.32% -smaller than the observed plastic threshold) the contribution to the distribution is mainly due to plastic events meaning that the weight of elastic events is negligible in this range of applied strain. This corresponds also to the strain interval ∆n at which the distribution becomes gaussian with a good approximation (non gaussian parameter less than 1), and where the departure from Lévy motion is of course noticeable. For smaller strains, the contribution of steps with an increase of the macroscopic stress (∆σ macro > 0) influences the distribution mainly for small values of the local stresses. These positive steps are for example responsible for the different power-law decay (∝ 1/∆σ 2.5 xy in place of ∝ 1/∆σ 1.7 xy ) in the beginning of the first crossover in the distribution (at the end of the initial plateau). We are unable to explain the precise value of these two exponents. However, we will now propose a model that is able to explain the generic features of this distribution of shear stress increments.

Minimal model: sum of quadrupoles

As shown in [21], the plastic flow regime consists in a succession of elastic branches and plastic events. We will now focus on the contribution of plastic events to the evolution equation of the stress components in the system. We identify here a plastic event by the fact that the associated stress is relaxed macroscopically (corresponding to an evident dissipation of energy shown by a negative slope in the stress-strain response of the system, see Fig. 1). Some of them correspond to an isolated quadrupolar redistribution of local stresses around a well identified center inside the sample; other plastic events consists in an alignment of rotational structures in the particle displacements giving rise to elementary shear bands going through the system [21]. The two kinds of events coexist in the plastic flow regime, and the transition from one kind of plastic event to another is still matter of debate [19,34,[START_REF] Marmottant | [END_REF][START_REF] Graner | Marmottant[END_REF]. However, the isolated quadrupolar event is more frequent than the elementary shear band. In the linear initial regime for example, only the isolated quadrupolar events take place, while there are periods with only quadrupolar events in the plastic flow regime as well.

In the following, we take the quadrupolar isolated event as the elementary building block for explaining the plastic deformation of the material, and we propose a simple model that describes the plastic deformation (as shown in the redistribution of stresses) as a sum of uncorrelated quadrupoles of random amplitude A.

We first identify the distribution P (A) in our system. As shown by Picard et al. [START_REF] Picard | [END_REF], a quadrupolar event involves in 2D a long-range redistribution of stresses, due to a local pure shear. The corresponding stress change is of the form.

∆σ xy (r, θ) ≡ A r 2 0 r 2 , if r 0 ≤ r ≤ r max A , if 0 ≤ r ≤ r 0 (10)
where r 0 and r max are respectively the typical size of the quadrupole and the size of the system, A is the amplitude of the quadrupoles, and ∆σ xy denotes the incremental shear stress. We have neglected in this expression the quadrupolar angular dependence of the stress field and only considered its 1 r 2 spatial decay3 . The distribution of P (A) is measured in our data (see figure 8). It corresponds to the incremental stress at the center of the quadrupole, that is at the place where the displacement is maximum. The best fit is exponential, with P (A) = 1 2σA . exp(-|A| σA ) and σ A = 2 appears to be the characteristic amplitude of the quadrupolar event. It is independent on the system size, unlike the distribution of macroscopic stress release that is ∝ 1/V (figure 2), and is symmetric.

Using the exponential distribution P (A) shown above, we can now reproduce the distribution P (∆σ = s, ∆n = 1) of the incremental shear stress ∆σ xy averaged over the whole system. According to the previous analysis (see part III), the distribution for ∆n > 1 should then be reproduced by successive convolutions, assuming the absence of temporal correlations. Neglecting the angular dependence of ∆σ xy as in Eq. 10, we show that, for a given amplitude A, there is a simple bijective relation between the shear stress ∆σ xy and the radial coordinate r. We can thus write P A (∆σ xy )d∆σ xy = P (r)dr, with P (r)dr = 2πr/V dr. Using Eq.10, we get

P A (∆σ xy ) = πAr 2 0 V ∆σ 2 (11) with |A| r 2 0 r 2 max ≤ |∆σ xy | ≤ |A| (12)
In the following, we will consider only the case A > 0, the opposite case giving the symmetric distribution corresponding to ∆σ xy < 0. Assuming that the successive events contribute independently to the total distribution P (∆σ xy ) (an assumption justified by the fact that this distribution is averaged over the time origins), the average probability of having an incremental stress ∆σ xy inside the system is thus obtained by summing over all quadrupolar events P (∆σ xy ) = c ∆σxy r 2 max r 0 ∆σxy P A (∆σ xy )P (A)dA (13) with an additional normalisation factor c due to the limited range of allowed amplitudes. This gives

P ( ∆σ xy ) = c πr 2 0 2V 1 ∆σ 2 xy ((σ A + ∆σ xy )exp(-∆σ xy /σ A ) -( σ A + σ xy r 2 max r 2 0 )exp(-∆σ xy .r 2 max /(σ A .r 2 0 ))) ( 14
)
We show here that the upper exponential cut-off is proportional to exp(-∆σ xy /σ A ) and that the low ∆σ behavior4 is dominated by finite size scaling proportional to exp(-∆σ xy .r 2 max /(σ A .r 2 0 )). We have plotted this function in figure 9, and compared it with the result obtained numerically. We find a good agreement with the numerical result, for σ A = 2, r max = L = 100 and r 0 = 4. The parameter r 0 is the unique free parameter of the fit since r max scales like the system size and the fitted value for σ A corresponds to the value obtained in the measured distribution of P (A) (see figure 8), thus confirming our simple model. The small value obtained for r 0 shows that the plastic quadrupolar events are localized.

This simple model gives us a physical intuition on the origin of the distribution for the increments η observed in the averaged evolution equation 9. Therefore we can now make more precise the evolution equation at a local level as

∂(σ( -→ r , ǫ) -σ(ǫ)) ∂ǫ = d 3-→ r Q A(ǫ)r 2 0 || -→ r - ---→ r Q (ǫ)|| 2 ρ( -→ r Q , ǫ) (15)
where -→ r Q is the position where the quadrupolar event has occurred,|| -→ r --→ r Q || ≥ r 0 , and ρ( -→ r Q , ǫ) is the number of plastic event per unit strain and unit volume, taking place at -→ r Q when the external strain is ǫ.

As it appears clearly here, this equation is not complete since -→ r Q (ǫ) is unknown. In many models [11,10]) the occurrence of a quadrupolar event at position -→ r Q is triggered by a local stress threshold. In the following section, we discuss the existence of such a criterion.

Looking for a local threshold

The center of the quadrupolar plastic rearrangements is identified as the place where the particle displacement is the largest [21]. In order to relate the position of the center of the quadrupoles to local stress criteria, we have compared the distribution of the stress components obtained for the whole system, and the value of the same stress component on the particle in the center of the quadrupole, one step before the quadrupole takes place.

The result for the deviatoric stress σ dev ≡ σ 1σ 2 is shown in figure 10(a). The distribution of the deviatoric stress in the whole sample is stationary: it is the same in the beginning (for ex. one step before the first plastic event) and in the end (e.g. one step before the last measured plastic event) of the plastic flow. It is remarkable to see in this figure, that the distribution of the deviatoric stress restricted to the center of the quadrupoles (one step before each plastic event) computed for all the plastic events is very close to the previously discussed stationary distribution obtained over the whole sample, and not only in the center of the quadrupoles. The position of the center of the quadrupoles appears to be unrelated with any threshold value in the local stress components. The saturation of the macroscopic stress, and the associated well defined yield stress σ Y in these systems can thus be related only to the alternance of increase and decrease of stress and to its intermittent behavior, rather than to any identified local stress threshold. The existence of an apparent macroscopic yield stress is not related to a local yield stress. Note that criterion based on Tresca or Mohr-Coulomb criterion, involving the most probable relation between deviatoric stress and pressure have also been tested in our numerical systems. It shows a general tendency for the pressure to be affinely related to the deviatoric stress in the plastic flow regime. However, this tendency is shared by all the particles in the system, and not only by the particles at the center of the quadrupoles. It thus appears, like in the case of the deviatoric stresses, that a global Mohr-Coulomb criterion, or equivalently -for 2D systemsa Tresca criterion, is valid in the plastic flow regime, but is does not provide a locally selective criterion for plasticity.

Figure 10(b) shows the same distributions as for the deviatoric stress, but for the incremental deviatoric stresses: it means σ dev (n+1)-σ dev (n). In this case, quenched stresses are not taken into account, and a large incremental stress is generally the signature of a large local deformation. The distribution of incremental stresses is also stationary. The distribution restricted to the center of the quadrupoles is displaced to larger increases in the deviatoric stresses. This point has also been mentioned by Robbins et al. [12] and underlines the role played by incremental stresses in the dissipative dynamics of the systems, in comparison with total stresses that are more or less unchanged. It suggests that the incremental stress ∆σ can be a more important parameter for plastic purposes than the stress σ itself. However, its distribution is very broad, and a criterion based only on the incremental stresses is not very selective: it does not even exclude the possibility for the center of the quadrupole to take place where the increase of the deviatoric stress is not maximum but minimum... The question of the existence of a local criterion for plasticity appears thus clearly beyond our detailed numerical study of local stress components, and deserves a further studies.

Conclusion

We have shown in this paper that the local stress components evolve in a heterogeneous way during quasi-static plastic deformation. It is related to a finite cooperativity number identifying fluctuating zones whose maximum size is approximately 6 × 6 particle diameters in the plastic flow regime. The strain at which these zones form (≈ 4%) is related to the typical strain at which the memory of the local stress seems to be lost in the 2-point correlation function.

It is thus natural to propose a picture of the heterogeneous stress evolution, as a sum of localized plastic events. We have shown in this paper that the averaged evolution of the full distribution of local stress components can be described by a simple model. Within this model, the evolution of local stress components results from the uncorrelated sum of localized quadrupolar rearrangements. The model allows to recover the 3 regimes shown in the distribution: its 1/V dependence at small stresses, in agreement with a local description of the plastic instability, the cross-over to an apparent power-low decay ∝ 1/∆σ α+1 with α close to 1 in the intermediate range, and for larger stresses a sizeindependent upper cut-off directly related in our model to the size-independent distribution of stress rearrangements in the center of the quadrupoles.

However, a complete description of local stress rearrangements is still lacking: for example, the observed exponential nature of the distribution of the amplitudes of stress rearrangements in the center of the elementary quadrupoles -a crucial ingredient of our model -is still not explained, and the position of the center of the quadrupoles should be related to a criterion of instability not identified in this paper. This last question especially deserves further studies, since is appears clearly here that the local stresses are not determinent.

Any attempt to describe local plasticity in terms of local stress rearrangements should thus take into account these numerical results before proposing a complete description of the heterogeneous mechanical response in terms of a single evolution equation.

Fig. 1 .Fig. 2 .Fig. 3 .

 123 Fig.1. stress-strain response and incremental shear stress,for one configuraiton of 10 000 particles.

Fig. 4 .Fig. 5 .

 45 Fig. 4. (a) Cooperativity number as a function of time (applied shear strain) for various stress components: local shear stress σxy, local pressure P , eigenvalues of the stress tensor σ1 and σ2, and the deviatoric stress σ dev = σ1 -σ2 (from bottom to top). The curves have been obtained by averaging on 20 000 time origins t0. (b) Cooperativity number of the local shear stress at the borders, in comparison with the cooperativity number in the center of the sample

Fig. 6 .= 1 ,

 61 Fig. 6. (a) distribution P (∆σxy, ∆n) of the variation of local shear stress ∆σxy for different strain intervals ∆n. The distribution has been obtained by averaging over the whole sample, and for different time origins. (b) lin-log scale for P (∆σxy, ∆n = 1) and different system sizes N = 10000(104 × 104) and N = 23125(50 × 500). It shows a size independent exponential upper cut-off with a characteristic ∆σ ≈ 1.4. (c) log-log scale for P (∆σxy, ∆n = 1) and different system sizes N = 10000(104 × 104) and N = 23125(50 × 500). It shows the scaling of the first cross-over ∝ 1/V where V is the volume of the system. The stress variation is mainly due to local variation that averages over the whole system (see text).

Fig. 7 .

 7 Fig. 7. (a) scaling of P (∆σxy, ∆n) ∝ φ(∆σxy/∆n 1/z ) with z = α and α = 0.7. The scaling is very good in the intermediate range of the distribution, and for small ∆n < 64 where the exponential cut-off does not play a role. Inset: the same with α = 1 = z. This exponent corresponds also to the one discussed in the text. (b) separate contribution of steps associated with a release of stress at a macroscopic level, to the total distribution of local stresses, for time intervals ∆n = 1 and ∆n = 64. As shown here, the distribution of ∆σ for ∆n ≥ 64 is dominated by the contribution of stress changes due to plastic events.

Fig. 8 .

 8 Fig. 8. Distribution of shear stress releases in the center of the quadrupoles, during a plastic event. The fit is exponential with a characteristic ∆σ = 2.

  with r 0 = 4; σ A = 2; r max = 100

Fig. 9 .

 9 Fig. 9. Comparison of the distribution P (∆σxy, ∆n = 1) with the theoretical fit discussed in the text. The parameters of the fit are indicated in the legend box.

Fig. 10 .

 10 Fig. 10. (a) Distribution of the deviatoric stress in a configuration of 10 000 particles for 2 different strains in the beginning and in the end of the plastic flow plateau (triangles).Comparison with the distribution of the deviatoric stress in the center of the future quadrupole just before the quadrupolar event takes place (squares), computed on the 2500 plastic events of the plastic flow regime. The 3 curves are very similar. (b) Same as in (a), but for the incremental deviatoric stress. In this case, the distribution obtained on the center of the quadrupoles indicates a more pronounced correspondance with the highest incremental stresses inside the sample.

In a previous paper we have analyzed in a similar way the motion of a particle in the system under shear (see Ref.[21])

The angular dependence will mainly contribute by a scaling factor, and affect the weight of small ∆σ in a logarithmic way in the final result. See M.Tsamados et al. preprint (2007) for a detailed calculation.

with a logarithmic divergence if the angular dependence is taken into account see M. Tsamados et al. preprint (2007).

Acknowledgments During the course of this work, we have benefited of discussions from L. Bocquet, M. Robbins, S. Roux, P. Sollich, and D. Vandembroucq. Computational support by IDRIS/France, CINES/France and CEA/France is acknowledged.


Florent Capelli
email: florent.capelli@univ-lille.fr

Yann Strozecki

Enumerating models of DNF faster: breaking the dependency on the formula size

In this article, we study the problem of enumerating the models of DNF formulas. The aim is to provide enumeration algorithms with a delay that depends polynomially on the size of each model and not on the size of the formula. We succeed for two subclasses of DNF formulas: we provide a constant delay algorithm for k-DNF with fixed k by an appropriate amortization method and we give a polynomial delay algorithm for monotone formulas. We then focus on the average delay of enumeration algorithms and show that we can bring down the dependency of the average delay to the square root of the formula size and even to a logarithmic dependency for monotone formulas.

1

Introduction

An enumeration problem is the task of listing a set of elements without redundancies, usually corresponding to the solutions of a search problem, such as enumerating the spanning trees of a graph or the satisfying assignments of a formula. One way of measuring the complexity of an enumeration algorithm is the total time needed to compute all solutions. When the total time depends both on the input and output, an algorithm is called output sensitive. It is considered tractable and said to be output polynomial when it can be solved in polynomial time in the size of the input and the output.

Output sensitivity is relevant when all elements of a set must be generated, for instance to build a library of interesting objects to be studied by experts, as it is done in biology, chemistry or network analytics [START_REF] Barth | Efficient generation of stable planar cages for chemistry[END_REF][START_REF] Andrade | Enumeration of minimal stoichiometric precursor sets in metabolic networks[END_REF][START_REF] Böhmová | Computing and listing st-paths in public transportation networks[END_REF]. However, when the output is large with regard to the input, output polynomiality is not enough to capture tractability. Indeed, if one wants only a good solution or some statistic on the set of solutions, it can be interesting to generate only a fraction of the solutions. A good algorithm for this purpose must guarantee that the more time it has, the more solutions it generates. To measure the efficiency of such an algorithm, we need to compute the delay between two consecutive solutions. A first guarantee that we can expect is to have a good average delay (sometimes referred to as amortized delay), that is, to measure the total time divided by the number of solutions. There are many enumeration algorithms which are in constant amortized time or CAT, usually for the generation of combinatorial objects, such as the unrooted trees of a given size [START_REF] Robert | Constant time generation of free trees[END_REF], the linear extensions of a partial order [START_REF] Pruesse | Generating linear extensions fast[END_REF] or the integers by Gray code [START_REF] Ervin | The art of computer programming[END_REF]. Uno also proposed [START_REF] Uno | Constant time enumeration by amortization[END_REF] a general method to obtain constant amortized time algorithms, which can be applied, for instance, to find the matchings or the spanning trees of a graph.

However, when one wants to process a set in a streaming fashion such as the answers of a database query, it may be interesting to guarantee the regularity of the delay, usually by bounding it by a polynomial in the input size. We refer to such algorithms as polynomial delay algorithms. Many problems admit such algorithms, e.g. enumerate the cycles of a graph [START_REF] Read | Bounds on backtrack algorithms for listing cycles, paths, and spanning trees[END_REF], the satisfying assignments of variants of SAT [START_REF] Creignou | On generating all solutions of generalized satisfiability problems[END_REF] or the spanning trees and connected induced subgraphs of a graph [START_REF] Avis | Reverse search for enumeration[END_REF]. They are all based on similar techniques (see [START_REF] Mary | Énumération des Dominants Minimaux d'un graphe[END_REF] for a survey).

For some problems, the size of the input may be much larger than the size of the generated solutions, which makes polynomial delay an unsatisfactory measure of efficiency. In that case, algorithms whose delay depends on the size of a single solution are naturally more interesting than polynomial delay or output polynomial algorithms. We say that an algorithm is in strong polynomial delay when the delay between two consecutive solutions is polynomial in the size of the last solution. To make this notion robust and more relevant, a precomputation step is allowed before the enumeration, in time polynomial in the input size, so that the algorithm can read the whole input and set up useful data structures. Observe that the notion of strong polynomial delay is also well suited for infinite enumeration where the size of the solutions grows arbitrarily [START_REF] Costa Florêncio | Naive infinite enumeration of context-free languages in incremental polynomial time[END_REF].

In this paper, we focus on the notion of strong polynomial delay. While strong polynomial delay is a very desirable property for an algorithm it is rarely considered, since polynomial delay is often seen as the be-all end-all answer in enumeration. We feel that understanding the difference between polynomial delay and strong polynomial delay is a challenge that complexity of enumeration must address to be more relevant for practitioners. We are especially interested in finding such algorithms for the problem of enumerating the models of DNF formulas. Indeed, one can easily solve it with delay linear in the size of the formula (see Theorem 4) but it appears that removing the dependency in the number of terms in the delay is not easy, even when allowing exponential space. DNF formulas have a simple structure as it is easy to find a large number of solutions or to approximate their number [START_REF] Luby | On deterministic approximation of dnf[END_REF]. Nevertheless, exactly counting the solutions of a DNF formula is a canonical #P-complete problem. Enumerating the models of a DNF formulas thus appears to be a reasonable candidate to separate the notions of strong polynomial delay and polynomial delay. Surprisingly, we manage to remove the dependency of the delay in the number of terms for k-DNF and monotone DNF and to improve it for the general case. Our results are summarized in Table 1.

There are few examples of strong polynomial delay in the literature. Constant delay algorithms naturally fall into this category and a whole line of research is dedicated to design such algorithm for enumerating models of first order queries for restricted classes of structures [START_REF] Schweikardt | Enumeration for fo queries over nowhere dense graphs[END_REF] (see also the survey [START_REF] Segoufin | Enumerating with constant delay the answers to a query[END_REF]). However, while these algorithms are called constant delay because their delay does not depend on the database size, it often depends more than exponentially on the size of the solutions. Other examples naturally arise from logic such as the enumeration of assignments of MSO queries over trees or bounded width graphs [START_REF] Bagan | Mso queries on tree decomposable structures are computable with linear delay[END_REF][START_REF] Courcelle | Linear delay enumeration and monadic second-order logic[END_REF]. Recently, Amarilli et al. [START_REF] Amarilli | A circuit-based approach to efficient enumeration[END_REF] presented an algorithm to enumerate the solutions of a restricted class of Boolean circuits known as structured d-DNNF used in knowledge compilation with delay linear in the number of variables and hence independent of the size of the circuit.

The paper is organized as follows. We first introduce basic notions on formulas and enumeration complexity then present tries and Gray code in Section 2. We show in Sec 3 how to adapt three classical methods to generate the models of a DNF, the best having a linear delay in the formula size. In Section 4, we give a backtrack search algorithm, using a good amortization, which is in constant delay for k-DNF. In Section 5, we give another backtrack search algorithm, whose average delay is better than linear. Finally, in Section 6, we provide a strong polynomial delay for monotone DNF formulas but with an exponential memory and we specialize and adapt the algorithm of the previous section to obtain a logarithmic delay for monotone DNF formulas. Proposition with omitted proofs due to space restrictions are labeled with ( ) and a proof is given in the appendix.

Definitions and notations

Terms and DNF-formulas. Let X be a set of variables and let n and be the size of X. We fix some arbitrary order on X and write X = {x 1 , . . . , x n }. A literal is either a variable x ∈ X or the negation of a variable ¬x for some x ∈ X. A term C is a finite set of literals such that every two literals in C have a different variable. A Disjunctive Normal Form formula, DNF for short, is a finite set of terms. Given a literal , we denote its underlying variable by var( ). We extend this notation to terms by denoting var(C) := ∈C var( ) for a term C and to DNF by denoting var(D) := C∈D var(C) for a DNF D.

Given an assignment α : X → {0, 1}, we naturally extend α to literals by defining α(¬x) = 1 -α(x). An assignment α satisfies a term C if for every ∈ C, α(C) = 1. A model α is an assignment that satisfies a DNF D, that is, there exists C ∈ D such that α satisfies C. We write α |= D if α is a model of D. It is easy to see that given a term C, there exists a unique assignment of variables in var(C) satisfying C. We denote this assignment by 1 C .

Given a DNF D on variables X, we denote by sat(D) = {α | α |= D} the set of models of D. Let Y ⊆ X, τ : Y → {0, 1} and σ : X → {0, 1}, we say that σ is compatible with τ , denoted by σ τ , if the restriction of σ to Y is equal to τ . We denote by sat(D, τ ) = {α | α |= D, σ τ } the set of models of D compatible with τ . We also denote by D[τ ] the DNF defined as follows: we remove every term C from D such that there exists a literal ∈ C such that τ ( ) = 0. For the remaining terms, we remove every literal whose variable is in Y . Observe that since we consider DNF to be sets of terms, by definition, D[τ ] has no duplicated terms. It is clear that:

sat(D, τ ) = {τ ∪ α | α |= D[τ ]}.
The size of a DNF D is denoted by D and is equal to C∈D |C|.

Enumeration complexity. Let Σ be a finite alphabet and Σ * be the set of finite words built on Σ. Let A ⊆ Σ * × Σ * be a binary predicate, we write A(x) for the set of y such that A(x, y) holds. The enumeration problem Enum•A is the function which associates A(x) to x.

The computational model is the random access machine model (RAM) with addition, subtraction and multiplication as its basic arithmetic operations and an operation Output(i, j) which outputs the concatenation of the values of registers R i , R i+1 , . . . , R j . We assume that all operations are in constant time except the arithmetic instructions which are in time linear in the size of their inputs. While this is arguably not a real-life scenario, this allows us to formally write algorithms with constant time delay while it is theoretically impossible if the time needed to write the output is taken into account. A RAM machine solves Enum•A if, on every input x ∈ Σ * , it produces a sequence y 1 , . . . , y n such that A(x) = {y 1 , . . . , y n } and for all i = j, y i = y j . The space used by the machine at a given step is the sum of the number of bits required to store the integers in its registers.

The delay of a RAM machine which outputs the sequence {y 1 , . . . , y n } is the maximum over all i ≤ n of the time the machine uses between the generation of y i and y i+1 . Note that we allow a precomputation phase before the machine starts to enumerate solutions, which can be in time polynomial in the size of the input.

Trie. A trie is a data structure used to represent a set of words on an alphabet M which supports efficient insertion and deletion. We refer the reader to [START_REF] Fredkin | Trie memory[END_REF][START_REF] Ervin | The art of computer programming[END_REF] for more details. We use them to either store a formula, that is the sets of its terms, or to represent a set of assignments of a formula.

The trie is a M -ary tree, each of its inner nodes is labeled by a letter of the alphabet. A leaf represents a word: the sequence of labels traversed when going from the root to the leaf. The trie represents the set of words represented by its leaves. The children of a node are given by a linked list.

A model α of the DNF is represented by the sequence of its values: α(x 1 )α(x 2 ) . . . α(x n ). The trie we use to store these models are binary trees, since the labels are in {0, 1}. We represent a term by the sequence in ascending orders of its literals, with for all i, x i < xi+1 < x i+1 . The list of children of a node in the trie is maintained in this order.

We need to search, insert and suppress elements in a trie. Those three operations are in O(n) for the two tries we have defined, where n is the number of variables.

Gray Code. Gray codes are an efficient way to enumerate the integers between 0 and 2 n -1 written in binary or equivalently the subsets of a set of size n. They enjoy two important properties: the Hamming distance of two elements in the Gray enumeration order is one and each new element is produced in constant time using only additional O(log(n) 2 ) space (see [START_REF] Donald | Combinatorial algorithms, part 1, volume 4a of the art of computer programming[END_REF]). In other words, this enumeration algorithm has constant delay. We can generate all models of a term in constant delay, using Gray code and an additional array which contains the indexes of the free variables of the term.

Proposition 1 (( )). The models of a term C on variables X can be enumerated in constant delay.

Classical enumeration algorithms

In this section we present three generic enumeration methods applied to the generation of the models of a DNF formula. The best one has a linear delay in the size of the instance, and we study in later sections several restrictions to obtain a delay polynomial in the size of a solution.

Union of terms

The first two methods are based on the fact that the models of a DNF formula are the union of the models of its terms. The only problem is to avoid repetitions of models. We first use a method to enumerate the union of sets of elements which preserves polynomial delay (see [START_REF] Strozecki | Enumeration complexity and matroid decomposition[END_REF]). It relies on a priority rule between the sets to avoid repetitions, as recalled in the proof of the following proposition.

Proposition 2 (( ) Adapted from proposition 2.38 and 2.40 in [START_REF] Strozecki | Enumeration complexity and matroid decomposition[END_REF]). The models of a DNF formula D with m terms can be enumerated with delay O(m D ).

By improving the way we test whether a model of a term is the model of any term of larger index, we can drop the delay to O(m 2 ). In fact, by generating the solutions of each term in the same order, we can avoid completely the redundancy test and get a better delay. The following algorithm merges several ordered arrays which are generated dynamically by enumeration procedures.

Proposition 3 (( ) Adapted from proposition 2.41 in [START_REF] Strozecki | Enumeration complexity and matroid decomposition[END_REF]). The models of a DNF formula D with m terms and n variables can be enumerated with delay O(mn).

The average delay of the two previous algorithms is lower than their delays. Let r be the average number of times a solution is produced during the algorithm, we can replace m by r in the average delay of the previous algorithm. It is possible to prove that r is smaller than m by studying how terms share models, but the complexity gains are very small and we give an algorithm with a much better average delay in Section 5.

Flashlight method

We present a classical enumeration method called the Backtrack Search or sometimes the Flashlight Method used in many previous articles [START_REF] Read | Bounds on backtrack algorithms for listing cycles, paths, and spanning trees[END_REF][START_REF] Strozecki | On enumerating monomials and other combinatorial structures by polynomial interpolation[END_REF] in particular to solve auto-reducible problems. We describe the method in the context of the generation of the models of D a DNF formula.

We define a tree T D whose nodes are the assignments τ over variables {x 1 , . . . x k } such that there is a model σ of D which is compatible with τ . The children of a node labeled by τ are the partial assignments τ defined over {x 1 , . . . x k+1 } which are compatible with τ .

The leaves of T D are the models we want to generate, therefore a depth first traversal visits all leaves and thus outputs all solutions. Since a path from the root of the tree is of size n, it is enough to be able to find the children of a node in polynomial time to obtain a polynomial delay. Hence the Flashlight Method has a polynomial delay if and only if the following extension problem is in P: given τ over {x 1 , . . . x k } is there σ a model of D compatible with τ ?

The extension problem for a DNF is very simple to solve: compute the formula D[τ ] and decide whether it is satisfiable in time O( D ). This yields an enumeration algorithm with delay O(n D ). The delay can be improved by using the fact that we solve the extension problem several times on very similar instance as it has been done for other problems [START_REF] Murakami | Efficient algorithms for dualizing large-scale hypergraphs[END_REF][START_REF] Mary | Efficient enumeration of solutions produced by closure operations[END_REF] such as the enumeration of the models of a monotone CNF. Proof. We use the previous algorithm which does a depth first search in T D . When it visits the node τ , we need to decide whether D[τ ] is satisfiable which is equivalent to testing whether it has a non falsified term. To speed-up the flashlight search, we use a data structure to decide quickly this problem and we need to guarantee that it can be updated fast enough when the tree is traversed.

For each term C, we store an integer f C which represents how many literals of the term are falsified by the current partial assignment. A term C is valid when f C = 0. For each literal l, we store the list of terms which contain the literal. We also store an integer vc for valid terms, which counts how many terms are not falsified by the current partial assignment.

At the beginning of the algorithm, all f C are set to 0 and vc = m. Then visiting a child τ of τ corresponds to choosing the value of some variable x k . If x k is set to 0 then for each term C containing x k , f C is incremented. The number of terms such that f C is changed from 0 to 1 is subtracted to vc. If x k is set to 1, then the same is done for xk . The fact that D[τ ] is satisfiable is equivalent to vc > 0. Remark that going up the tree when backtracking works exactly as going down, but the variables f c are decremented and vc is incremented instead.

As a consequence, the complexity of the algorithm over a path in the tree is O( D ) since for each term C, the variable f C will be modified |C| times. When the algorithm goes down the tree it may first set x k to 0 and fails, but then it sets x k to 1 and goes down. Hence the cost of going down to a leaf is at most twice the cost of following the path to the leaf. Since between two outputted solutions, the algorithm follows one path up and one down, the delay is in O( D ).

This last algorithm has a delay linear in the size of the formula, however the size of the formula can be extremely large with regards to the size of a model. In the next sections we will try reduce or eliminate the dependency of the delay in the size of the input either for particular DNF formulas or by relaxing the notion of delay.

Enumerating models of k-DNF

A term C is a k-term if and only if |C| ≤ k. A DNF is a k-DNF if all its terms are k-terms. In this section, we present an algorithm to enumerate the models of a k-DNF with a 2 O(k) delay. The idea is to select a k-term and use its 2 n-k models to amortize more costly operations. More precisely, we prove the following: Theorem 5. The models of a k-DNF with n variables can be enumerated with precomputation in O(n) and O(k 2 2 5k ) delay.

To explain our algorithm, we need first to introduce notations. Let H(p) = -p log(p) -(1p) log(1 -p) be the binary entropy. From the following classical inequality (see [START_REF] Flum | Parameterized complexity theory[END_REF]):

k i=0 n i ≤ 2 nH(k/n) ,
we get for k ≤ n/2: Lemma 6. The number of k-terms on n variables is at most 2 nH(k/n)+k .

Let D be a DNF-formula on variables X. Assume wlog that X is ordered with <. Given a term C ∈ D, we denote by

1 C : var(C) → {0, 1} the only model of C, that is, for every ∈ C, 1 C (x) = 1 if = x and 1 C (x) = 0 if = ¬x.
If y ∈ var(C), we denote by

0 y C : {z ∈ var(C) | z ≤ y} → {0, 1} the assignment defined by 0 y C (z) = 1 C (z) for z < y and z ∈ var(C) and 0 y C (y) = 1 -1 C (y). For example, if C = x 1 ∧ x 2 ∧ ¬x 3 , we have 1 C = {x 1 → 1, x 2 → 1, x 3 → 0} and 0 x 1 C = {x 1 → 0}, 0 x 2 C = {x 1 → 1, x 2 → 0} and 0 x 3 C = {x 1 → 1, x 2 → 1, x 3 → 1}.
These assignments naturally induce a partitioning of the model of a DNF: Lemma 7. Given a DNF D and a term C ∈ D, we have: We start by picking a k-term C ∈ D. It is easy to see that for every τ : X \ var(C) → {0, 1}, we have τ ∪ 1 C |= D. Thus, enumerating sat(D, 1 C ) can be done very efficiently with a O(1) delay by simply enumerating 2 X\var(C) using a Gray code. Observe that the precomputation time boils down to choosing a term and outputting its first solution which can obviously be done in time O(n) as stated in the theorem.

sat(D) = sat(D, 1 C ) y∈var(C) sat(D, 0 y C ).
Between the output of two solutions of sat(D, 1 C ), we spend some time to precompute D[0 y C ] for every y ∈ var(C). Since |var(C)| ≤ k, we have to compute at most k such DNF. Moreover, to compute each one of them, we have to inspect every term of D and every variable of these terms. Since each term has at most k variables, we need O(kM ) steps to precompute all terms in D[0 y C ] for one y ∈ var(C). We then remove all copies but one of the same terms. This can be done in time O(kM 2 ) by comparing the terms pairwise. In the end, we need O(k 2 M 2 ) steps to precompute D[0 y C ] for every y ∈ var(C). Let A be a constant such that this precomputation can be done in time A • k 2 M 2 steps at most.

Assume that between the output of two solutions of sat(D, 1 C ), we allow dA steps for this precomputation (the value of d will be fixed later depending on our needs). Since |sat(D, 1 C )| ≥ 2 N -k , this gives us a total amount of 2 N -k dA steps for this precomputation. Thus, if 2

N -k dA > Ak 2 M 2 , that is, if 2 N -k d > k 2 M 2 (1)
we have enough time to compute D[0 y C ] for every y ∈ var(C). If this is the case, then we do the precomputation, finish the enumeration of sat(D, 1 C ) and then recursively start the enumeration of 0 y C × D[0 y C ] for each y ∈ var(C). The number n of variables of D[0 y C ] has of course decreased but we can still allow for dA steps of extra computation between the output of two solutions.

It is clear that as long as ( 1) is true, we can output solutions of D with delay dA. We may however have a problem if, at some point in the recursion, the DNF D we are using to compute solutions of D may have too few variables making (1) false. We are going to pick d sufficiently large so it never happens.

Since we only consider k-DNF, we know by Lemma 6 that at each step of the recursion with n unassigned variables, the number m of terms verifies m ≤ M 0 := 2 nH(k/n)+k . Our goal is thus to find d such that for every n ≥ 2k:

2 n-k d ≥ k 2 M 2 0
By rearranging the terms so that k is in factor in the exponent and denoting x = k/n, we conclude that if the following holds then ( 1) is true when k ≤ n/2:

d ≥ k 2 2 k((2H(x)-1)/x+3) (2)
It is easy to compute the maximum of the function (2H(x) -1)/x + 3 over [0, 1/2]. We obtain a bound of 5 at x = 1/2. Hence for all k ≤ n/2, if d ≥ k 2 2 5k then ( 1) is true at each recursive call. When n < 2k, we enumerate the models using the algorithm of Theorem 4. It gives a delay of D ≤ k2 3k since by applying Lemma 6 for n = 2k, we have that that the number of k-terms with 2k variables is at most 2 2kH(1/2)+k = 2 3k . Thus, a k-DNF with less than 2k variables has size at most k2 3k .

The O(k 2 2 5k ) delay of the previous algorithm could be improved to O(k 2 2 3k ) by using a more clever algorithm to detect duplicated terms when computing D[0 y C ]. This operation can easily be done in time O(k 2 M log(M )) by sorting the terms before eliminating duplicates and even in time O(k 2 M ) by using a trie with relevant pointers allowing to test if a term is already in the trie in time O(k). Both solutions would make the presentation of the proof heavier so we chose to present this simplified version since the delay is also 2 O(k) .

Theorem 5 gives an algorithm to enumerate the models of a k-DNF with delay d = 2 O(k) and O(n) precomputation. A pseudo code for this algorithm is given in the appendix (Algorithm 1). It has constant delay for constant k and it is in strong polynomial delay (polynomial in n the size of a solution) for terms of size O(log(n)). We conjecture however that the models of a general DNF cannot be enumerated with strong polynomial delay.

We can use Algorithm 1 to significantly improve a result on the enumeration of models of first order formula with free second order variables. In [START_REF] Durand | Enumeration complexity of logical query problems with second-order variables[END_REF], Theorem 10, it is proved that the models of a Σ 1 formula (a single block of existential quantifiers followed by a quantifier free formula) with free second order variables can be enumerated in polynomial delay using a method similar to Proposition 2. Moreover, this problem is shown to be equivalent to the enumeration of models of a k-DNF. As a corollary of Theorem 5, the enumeration of models of a Σ 1 formula can be done in constant delay. It is very surprising since it is the same complexity as for Σ 0 (quantifier free formulas), the first level of the hierarchy.

Average delay of enumerating the models of DNF

In this section, we analyze the average delay of the Flashlight method, using appropriate data structures, to show it is better than the delay. The idea is to amortize the cost of maintaining the formula D[τ ] during the traversal of the tree over all models of D[τ ] in the spirit of [START_REF] Uno | Constant time enumeration by amortization[END_REF]. To do that we exhibit a relation between the number of models of a DNF and its number of terms. For any class of DNF for which it is possible to guarantee, for all partial assignments τ , that the number of models of D[τ ] is large enough with regard to |D[τ ]|, we obtain a good average delay. We then improve the delay of the Flashlight method by using a different branching (in the spirit of Algorithm 1) which guarantees that the number of terms decreases sufficiently but not too much in some branch which allows us to amortize the cost of branching.

The following lemma shows that DNF whose terms share a model have many models. It applies directly to any monotone formula, since the all one assignment is always a model of a (positive) monotone term. Proof. Let α be a model of the m terms of the DNF formula D and x be a variable of D. We claim that x appears only negatively or positively in D. Indeed, assume w.l.o.g that that α(x) = 1 and that x ∈ C 1 and ¬x ∈ C 2 . Thus α |= C 2 which contradicts the definition of α. We construct a monotone DNF D by replacing every variable that appears only negatively in D by positive occurrences. There is an obvious one-to-one correspondence between the models of D and the models of D. Now, for C ∈ D , let m C be the model of C defined as m C (x) = 1 if and only if x ∈ var(C). Since D is monotone, it is clear that m C = m C if and only if C = C , which gives m distinct models for D and thus for D.

Using the previous lemma, we prove that a DNF has at most a number of terms which is quadratic in its number of models. Lemma 9. A DNF formula with m non empty distinct terms has at least m 1/2 models.

Proof. Let α be a model of the DNF D, we denote by r α the number of terms which are compatible with α: |{C |= α | C ∈ D}|. It corresponds to the number of redundancies of α, that is the number of terms which has it as model (and the number times it is generated in the algorithm of Proposition 3). Now, let us denote by A the average of r α over all models. By definition, we have

|sat(D)| = C∈D |sat(C)| A .
Since A is an average, there is an α such that r α ≥ A. It means that at least A distinct terms are compatible with α and then by Lemma 8 these terms have at least A distinct models. In the previous proof, the inequalities are not tight and we may certainly improve the result. In particular, we can avoid to bound C∈D |sat(C)| by m. We then get a lower bound on the number of models of m 1/2+c/ log(n) for some constant c. It improves the delay of the next algorithm only for superpolynomial m and require some additional computations, hence we chose to state Lemma 9 with the simpler √ m bound. The quality of the algorithms presented in this section is directly connected to this relation between the number of terms and the number of solutions of a DNF. Any improvement would have direct consequences on the complexity of enumerating the models of a DNF. Question 10. What is the best k such that we can guarantee that a DNF with m terms has at least m k solutions? We are not aware of any results concerning this fact. We only know that 1/2 ≤ k ≤ log 3 [START_REF] Andrade | Enumeration of minimal stoichiometric precursor sets in metabolic networks[END_REF]. The lower bound is Lemma 9 and the upper bound easily follows from the DNF having all terms of size at most n on n variables has 2 n models and 3 n terms.

To improve the average delay of the flashlight method, we need to use an adapted data structure. In particular we need that, when considering some inner node D[τ ] of T D , the cost to process it is in O(|D[τ ]|) and not in O(m). In particular, we need to guarantee that there are no redundancy of terms in the structure representing D[τ ] and that we can maintain it efficiently, that is why we use a trie.

Theorem 11. The models of a DNF can be enumerated with average delay O(n 2 √ m) and polynomial space.

Proof. We maintain the formula D[τ ] when we traverse T D using the trie containing its terms as explained in Section 2. On a node τ we can decide quickly whether D[τ ] has a model and we can maintain D[τ ] efficiently without redundancy of terms. In the flashlight search, we will fix the variables following their order x 1 , . . . , x k . Hence, visiting a child τ of τ corresponds to setting the value of some variable x k with all variables x i with i < k already fixed. If we set x k to 0 then we need to remove the subtree under the root of the trie, with first node x k . Then we remove the subtree under the root of the trie, with first node xk , and insert back all elements in this subtree into the trie without the first node xk . The complexity of the latter is in O(|D[τ ]|n) since the number of terms is bounded by |D[τ ]| (no term appear several times in the trie). To set x k to 1, we do the same operation where we exchange the roles of x k and xk . To be able to go up in the tree during the flashlight search, we must restore the trie to its previous state. To do that in time O(|D[τ ]|n), it is enough to store the list of elements which have been removed or added in the trie when going down the same edge and to reverse the operations. The additional memory used during the algorithm is bounded by O(mn 2 ).

We now compute the average delay, that is the total time of the algorithm divided by the number of outputted models. To do that we distribute the time spent on each inner node A leaf receives a cost of at most n √ m for each of its n ancestors, hence the time per leaf (or the average delay) is bounded by O(n 2 √ m).

We now improve the way we choose the next variable on which we branch in order to improve the delay. The idea is to branch on variables so that the recursion tree kept as balanced as possible. We use the following lemma: Lemma 12 (( )). Let D be a DNF formula with m terms. Let m 0 and m 1 be the number of terms of D[x → 0] and D[x → 1] respectively. Then m 0 + m 1 ≥ m/2. Theorem 13. The models of a DNF can be enumerated with average delay O(n √ m) and polynomial space.

Proof. Instead of branching on a variable x 1 in the flashlight search, the branching will have the structure of a potentially large comb, similar to the one used in Algorithm 1. We define a sequence of partial assignments α i , α i with 0 ≤ i ≤ n. The assignments α i and α i are defined over the variables {x 1 , . . . , x i } and they are compatible with α i . The assignment α 0 is the empty assignment, α i is the extension of α i-1 by x i → (x i ) and α i is the extension by

x i → 1 -(x i ). We write m i = |D[α i ]| and m i = |D[α i ]|,
these numbers are the number of distinct terms in D[α i ] and D[α i ] respectively. Let us now define the assignment , given D[α i-1 ], (x i ) is chosen to maximize m i .

By Lemma 12, we have that m i ≥ m i-1 /4 that is the number of terms is divided by at most 4 at each step. Since m 0 = m and m n ≤ 1, then there is a i 0 such that

m i 0 ∈ [m/8, m/2].
We use the same method as in Theorem 11 to maintain the trie representing D[τ ]. During the successive construction of D[α i ] and D[α i ], each term of D will be removed or inserted twice at most in the trie. Hence the total cost of this branching and of the computation of the m i and m i is in O(mn), that is as fast as a binary branching. We assign this cost to the solutions in the branch α i 0 . By construction, D[α i 0 ] has more than m/4 terms, hence it has by Lemma 9 at least m/8 models. Therefore the cost charged to the leaves is n √ m. Now let us analyze how much a leaf is charged. Remark that in the path from the root to the a leaf, not all nodes charge the leaf. By construction the nodes which charge the leaf have less than half the terms of the previous node which has charged the leaf. Hence the complexity charged to a node is bounded by

i n m/2 i = n √ m i 1/2 i/2 . Since i 1/2 i/2
is bounded by a constant, the cost of a leaf and thus the average delay is O(n √ m).

Enumerating models of monotone DNF

When the underlying formula is monotone, that is it does not contain any negated literal, we can enumerate the models with a delay polynomial in the number of variables only. However, our current techniques need an exponential memory to work.

Theorem 14. There is an algorithm that given a monotone DNF with n variables and m terms, enumerate the models of D with preprocessing O(nm 2 ) and delay O(n 3 ). The space needed for this algorithm is linear in the number of solutions of D.

Proof. We start by removing from D every term C such that there exists C ∈ D with C ⊆ C . Observe that it does not change the models of C since C ⇒ C. This preprocessing phase takes O(nm 2 ) since comparing two terms may be done in time O(n). Let D be the resulting minimized monotone DNF. The algorithm then work as follows: arbitrarily order the terms of D = {C 1 , . . . , C p } and its variables X = {x 1 , . . . , x n }. We initialize a trie T that will contain the solutions that we have already enumerated, that is, in the following algorithm, each time we output a solution, we store it in T , so we can check in time O(n) if a solution has already been enumerated.

We now enumerate the solutions as follows: we start by enumerating all solutions of C 1 . Then, we proceed by induction: once we have enumerated all solutions of

D i = C 1 ∨ • • • ∨ C i ,
we enumerate all solutions of C i+1 that are not solutions of D i until i = p. Once we are done, we have, by induction, enumerated all models of D.

We claim that we can do this with delay O(n 3 ) using a classical reverse search method. Let Y = (y 1 , . . . , y m ) be the variables that are not in var(C i ), ordered following the natural order we have on X. The solutions of C i are in one-to-one correspondence with 2 Y as follows: given S ⊆ Y , we have a solution m S defined as m S (x) = 1 if x ∈ var(C i ) ∪ S and m(x) = 0 otherwise. We explore the solution of C i by following a tree A whose nodes are labeled with m S for every S ⊆ Y . The root of the tree is label-led by m ∅ and for every S, the unique predecessor of node m S is m S with S = S \ max(S). In other words, given S, the successors of m S are m S∪{x k } for k > max(S).

We enumerate the solution by following the structure of A. We start from the root of A, that is, we enumerate m ∅ . We claim that m ∅ has not yet been enumerated. Indeed, assume toward a contradiction that m ∅ |= C j for j < i. Then since D is monotone, it means that C j ⊆ C i which is absurd since D is minimized. Thus, the first model we enumerate is guaranteed to be fresh. Now we follow the structure of A by depth first search. When we visit a node of A label-led with m S , we start by checking in the trie if m S has already been enumerated. It can be done in O(n). If m S has not been enumerated, then we output it and keep on going down in the tree A. If m S has been enumerated then we can discard the whole subtree of A rooted in m S . Indeed, by definition, every solution m S in this subtree verifies S ⊇ S. Thus, if m S has been enumerated then m S |= D i-1 and since D i-1 is monotone, m S |= D i-1 , so m S has already been enumerated.

Thus, if m S has already been enumerated, we backtrack in the tree to the next node in A that has not yet been explored and that is not in the discarded subtree. This can be done in O(1) if we maintain in each visited node a pointer to the first ancestor that has an unexplored child.

We claim that we will visit at most n 2 nodes of A before finding a new solution or realizing that we have enumerated all solutions of D i and should proceed with C i+1 . Indeed, by the previous argument, we only have to check the maximal unexplored nodes of A since if a solution has already been enumerated, so have all its descendants. Now, each node of A has at most n children and A is of depth n. Hence, during the depth-first search, we will have at most n 2 such maximal nodes (at most n per level). Since checking if a solution is already in the trie can be done in O(n), our algorithm will have delay at most O(n 3 ). By using a better data structure to store the solutions and another organization of the models of a term we should bring the delay down to O(n 2 ) or better, but we leave that for further researchs. Moreover, observe that, using the reduction in the proof of Lemma 8 to the monotone case, the algorithm described in Theorem 14 also works if every variable appears only positively or negatively in all terms of the formula.

We now consider algorithms with a good average delay for enumerating the models of a monotone DNF formula. This allows to obtain a better delay that the previous theorem, while only using a polynomial space. First, recall that monotone DNF formulas have at least as many models as terms because of Lemma 8. Hence, using the algorithm of Theorem 13 on a monotone formula, we obtain the following theorem.

Theorem 15. The models of a monotone DNF can be enumerated with average delay O(n) and polynomial space.

To improve the bound on the average delay using a similar algorithm, we should either guarantee a better relationship between the number of terms and solution or we should reduce the complexity of maintaining the trie during the algorithm. Note that the formula with all positive terms has 2 n models but also 2 n terms. If we further assume that no term are redundant, that is there are no C 1 , C 2 such that C 1 ⊆ C 2 , then the formula with all terms of size n -1 has n + 1 models and n terms. Even when m is large, the relationship is almost linear: the formula with all terms of size n/2 has 2 n /2 models and O(2 n / √ n) terms. Hence, to improve the average delay, it seems hard to rely on a better bound on the number of solutions. However, the cost to deal with the trie can be reduced, as long as m is not too large as shown in the next theorem.

Theorem 16. The models of a monotone DNF can be enumerated with average delay O(log(n)(log(m)+ log(n)) and polynomial space.

Proof. First remark that when a term contains n -k variables, then it has 2 k models. Hence, in the algorithm of Theorem 13, when we compute the cost distributed over all solutions of D[τ ], we can make the following adjustment. We consider two cases. First say that D[τ ] has n τ variables and assume there is a term with n τ -k variables such that k > log(|D[τ ]|) + 2 log(n). Then D[τ ] has at least 2 k models which is larger than |D[τ ]|n 2 . As a consequence, the complexity which is charged to the leaves is only 1/n and the total cost charged to a leaf by branching steps of this king is bounded by 1. Now assume that for some partial assignment τ , all terms in D[τ ] contain more than n τ -k literals, with k = log(|D[τ ]|) + 2 log(n). Then we change the encoding of the formula: each term of D is represented by its complementary. Since the formula is monotone, we need only to store

The complementary terms are stored in a trie where the children of a node are organized in an AVL tree instead of an ordered list. All operation can be done on this trie in time O(log(n) * k), where k is the size of the object to insert. The algorithm works as the one in Theorem 13 on this new data structure with the following difference in branching. When x i is set to 1, we remove the subtree under the roof with first node x i and add back all elements without x i . When x i is set to 0, we keep only the subtree under the roof with first node x i and remove the rest.

The cost of a step is now bounded by O(|D[τ ]|k log(n)) and the cost charged to the leaf is bounded by k log(n), that is by (log(m) + 2 log(n)) log(n) which proves the proposition.

A better data structure can be used in the previous theorem: a trie where the children of a node are stored in an array of size n. Then the operations are in time O(k) if we accept that we have an arbitrary supply of initialized memory. The average delay then drop to O(log(n) + log(m)) and we conjecture that there is a good data structure allowing such delay without needing initialized memory.

The result on monotone DNFs can easily be transferred to a problem studied in [START_REF] Mary | Efficient enumeration of solutions produced by closure operations[END_REF]: the generation of all unions of given subsets. An instance {s 1 , . . . , s m } is a set of m subsets of [1, n] and we want to generate all distinct unions of these s i . The delay of the algorithm in [START_REF] Mary | Efficient enumeration of solutions produced by closure operations[END_REF] is O(nm), using a flashlight search algorithm similar to the algorithm of Theorem 4. Among the enumeration problems captured by the framework of saturation by set operators [START_REF] Mary | Efficient enumeration of solutions produced by closure operations[END_REF] it is the only one not proved to be in strong polynomial delay and it is also proved to be at least as hard to enumerate as the models of a monotone DNF.

If there are m distinct subsets in the instance, these subsets are also solutions then we have an equivalent of Lemma 8. Moreover, if we restrict a set of subsets by fixing an element, we have an inequality similar to Lemma 12. Hence we obtain an algorithm to enumerate union of sets with an average delay O(n) using the branching scheme of Theorem 13. However, the algorithm of Theorem 14 do not seem usable for generating the union of sets, leaving the question of obtaining an algorithm with strong polynomial delay open.

have the only allowed value by C. The delay between two solutions is constant since we only look the value of σ(i) in an array and switch the value of a register between two outputs.

B Proof of Proposition 2

Proof. As explained in Proposition 1, a term C ∈ D, sat(C) can be enumerated in constant delay using Gray code enumeration. The terms of D are indexed from C 1 to C m in an arbitrary order. During the algorithm, the state of the enumeration by Gray Code for each term is maintained so that we can query in constant time the next model of a given term. At each step, the algorithm does a loop from C 1 to C m . For a term C i the next model is generated and the algorithm tests whether it is a model of some C j with j > i. If not, it is outputted. If a term has no more models we skip it.

By this method, we guarantee that each model is outputted when generated by the term of largest index it satisfies, hence all models are generated and without repetitions. Moreover, at each step of the algorithm, the model given by the last term which has still models will be outputted, therefore the delay is bounded by the time to execute one step of the algorithm.

The cost of the generating new models at each step is bounded by O(m) since each solution is produced in O(1). The cost of testing whether a model satisfies some term of larger index is bounded by D and it is done at most m times before outputting a solution which implies that the delay is O(m D ).

C Proof of Proposition 3

Proof. As in the previous algorithm, we run a simple enumeration algorithm on each term and we maintain their states so that we can easily query the next model of a term. We chose to enumerate the models of the terms in lexicographic ascending order (for some arbitrary order of the variables), which can be done with delay O(n) for each term.

The first model of each term which has not yet been outputted is stored in a trie; if all models of a term have been outputted then nothing is stored for this term. Moreover, for each model in the trie, we maintain the list of terms from which it has been generated.

At each step of the algorithm, the smallest model α is found in the trie, then outputted and removed from the trie all in time O(n). Then we use the list of terms which had α as a model, to generate for each of them their next model and add it to the trie in time O(mn) since the insertion can be done in time O(n) and the number of new models is bounded by O(m). The delay of this algorithm is thus bounded by O(mn). By induction, we prove that at each step the smallest non outputted model is outputted, which implies that all models are outputted without repetitions.

D Pseudocode for Theorem 5 E Proof of Lemma 12

Proof. Denote by n 0 the number of terms of D which contain x, n 1 the number which contain x and n the number which do not contain variable x. The terms of D[x → 0] are the terms of D which contains the literal x or which does not contain the variable x. Remark that fixing the variable x

Proposition 4 .

 4 The models of a DNF formula D can be enumerated with delay O( D ).

Proof.

  The right-to-left inclusion is clear as sat(D, τ ) ⊆ sat(D) for any τ . Moreover, these unions are clearly disjoint since for every y, z ∈ var(C), y < z, we have 0 z C (y) = 0 y C (y) and 1 C (y) = 0 y C (y). For the left-to-right inclusion, let τ ∈ sat(D). If τ 1 C , then τ ∈ sat(D, 1 C ). Otherwise, let y be the smallest variable of var(C) such that τ (y) = 1 C (y). Then we have τ 0 y C . Proof (of Theorem 5). Given a k-DNF D on variables X, we use Lemma 7 to enumerate sat(D). We denote by X = var(D), N = |X| and M = |D|.

Lemma 8 .

 8 A DNF formula D with m distinct terms which have a common model has at least m models.

  Hence we have |sat(D)| ≥ A. Now, multiplying both equations and observing the fact that |sat(C)| ≥ 1 for every C ∈ D, we obtain |sat(D)| 2 ≥ C∈D |sat(C)| ≥ m.

  D[τ ] to the models in the leaves of the subtree rooted at D[τ ]. As we have explained, the time spent on the node D[τ ] (to do the branching on some variable x k ) is in O(|D[τ ]|n). By Lemma 9, we have that D[τ ] has at least |D[τ ]| models. The time spent in D[τ ] is distributed uniformly over all leaves of the tree rooted at D[τ ]. Hence each leaf receive at most n |D[τ ]| = n|D[τ ]|/ |D[τ ]|.

Algorithm 1 :

 1 Enumerates the models of k-DNF with delay 2 O(k) . Data: A k-DNF-formula D begin if D = ∅ then return ∅ ; if D = {C} then Enumerates the models of C ; else d ← Ak 2 2 kα ; Pick C ∈ D ; Every d steps of computation in the next block, output a new model of D[1 C ] ; begin for y ∈ var(C) do D y ← D[0 y C ] ; for y ∈ var(C) do Recursively enumerates 0 y C × sat(D y );to 0 cannot make two terms equal except when one contains x and the other does not contain the variable x. Hence m 0 ≥ max(n 0 , n ). Similarly,m 1 ≥ max(n 1 , n ). Since max(a, b) ≥ (a + b)/2, we have m 0 + m 1 ≥ (n 0 + n )/2 + (n 1 + n )/2 ≥ (n 0 + n 1 + n )/2. By definition m = n 0 + n 1 + n , then we have proved the lemma.

Table 1 :

 1 Overview of the results. In this table, D is a DNF, n its number of models and m its number of terms. New contributions are annotated with ( ).

	Class	Delay

A Proof of Proposition 1

Proof. The models of C on variables X are exactly the assignments of the form 1 C ∪ τ for any τ : X \ var(C) → {0, 1}. Enumerating the models of C intuitively boils down to enumerating all assignments on variables X \ var(C). We use a Gray code to do this in constant time.

More precisely, we represent a model of C in n registers R 1 , . . . , R n where R i holds the value of x i in this model. We initialize the registers such that

∈ var(C). Now let k = |X \ var(C)| and σ : {1, . . . , k} → {1, . . . , n} be such that X \ var(C) = {x σ(1) , . . . , x σ(k) } with σ(1) < • • • < σ(k). We start by storing the values of σ in an array of size k. We then execute Output(1, n) which outputs the first model of C. After that, we run a Gray code enumeration on the subsets of a set of size k. For each new element generated, the Gray code algorithm flips a bit at some position i. We thus switch the bit of R σ(i) and call Output (1, n). This generates all solutions of C since all possible values of the variables not in C are set, while the other