Tania Richmond
email: tania.richmond@univ-st-etienne.fr

Pierre-Louis Cayrel
email: pierre.louis.cayrel@univ-st-etienne.fr

Viktor Fischer
email: fischer@univ-st-etienne.fr

Pascal V Éron

Towards a Secure Goppa Decoder in Hardware Implementation

The irreducible binary Goppa codes are widely used in code-based cryptography, like in the McEliece cryptosystem. The aim of this work is to design an efficient and secure hardware implementation of a Goppa decoder. We will show how to adapt a common step of all decoding algorithms to obtain a "leakage resistant" variant.

Introduction

There exists several domains in cryptography. Code-based cryptography has interesting properties: Fast computations, a well theorical complexity and no quantum algorithm in polynomial time breaks the underlying NP-hard problem [START_REF] Berlekamp | On the inherent intractability of certain coding problems[END_REF]. The first public-key cryptosystem (PKC) based on error-correcting codes was proposed by McEliece in 1978 in [START_REF] Mceliece | A public-key cryptosystem based on algebraic coding theory[END_REF]. In his paper, he proposed to use the family of Goppa codes. Permuted Goppa codes present some advantages: they look like random codes, have at least one efficient decoding algorithm, and it is a dense family of codes. Goppa decoder can be used as a trapdoor in code-based cryptography, both for encryption scheme (like McEliece public-key encryption scheme [START_REF] Mceliece | A public-key cryptosystem based on algebraic coding theory[END_REF]) and signature scheme (like CFS signature scheme [START_REF] Courtois | How to Achieve a McEliece-Based Digital Signature Scheme[END_REF]). It is in this context that takes place our work. We are interested in the way that Goppa decoding algorithms are translated into circuits in hardware. Furthermore, we propose a partial hardware implementation of them. The remainder of this paper is organized as follows. In Section 2, we will recall some preliminaries about the McEliece PKC and the Goppa codes. In Section 3, we will briefly present the principle of side-channel attacks. In Section 4, we will explain our hardware implementation of a partial Goppa decoder. Finally, we will conclude this paper in Section 5.

McEliece PKC with Goppa codes

The reader who wants more explanations about the McEliece PKC or the Goppa codes can find them in Sections 3 and 4 in [START_REF] Heyse | Code-based cryptography on reconfigurable hardware: tweaking niederreiter encryption for performance[END_REF] or directly in [START_REF] Mceliece | A public-key cryptosystem based on algebraic coding theory[END_REF][START_REF] Goppa | A new class of linear correcting codes[END_REF][START_REF] Berlekamp | Goppa codes. Information Theory[END_REF].

1. Compute C = C • P -1 . 2. Decode C with a Goppa decoder to obtain M = M • S. 3. Compute M = M • S -1 . 4. Return M .

Goppa codes

Definition: In this paper, we are only interested in irreducible binary Goppa codes. A Goppa code C is a linear code generated by the Goppa polynomial g and the support L. g is a monic irreducible polynomial of degree t. L = {α1, . . . , αn} is a subset of n elements of the Galois field F2m with g(αi) = 0 ∀i ∈ 1, n . C is the set of words c of length n which verify the following equivalence:

- n i=1 c i g(α i ) × g(x)-g(α i )
x-α i ≡ 0 mod g(x).

Decoding: The aim of a Goppa decoder is to solve the key equation as follows: Find the so-called error locator polynomial σ such that SC (x) × σ(x) = σ (x) mod g(x). All σ roots correspond to all error positions. That is why, to decode a codeword C into a message M , there are the three following steps:

1. Compute the syndrome polynomial:

SC (x) = - n i=1 C i g(α i ) × g(x)-g(α i )
x-α i . 2. Solve the key equation to obtain the error locator polynomial: σ(x). 3. Find the roots βi ∈ L of the error locator polynomial:

σ(x) = t i=1 (x -βi).
The second step can be executed by extended Euclidean algorithm, Berlekamp-Massey algorithm [START_REF] Massey | Shift-register synthesis and BCH decoding[END_REF] or Patterson algorithm [START_REF] Patterson | The Algebraic Decoding of Goppa Codes[END_REF].

Side-channel attacks

Examples. The main side-channel attacks are: timing attack and power consumption attack. The first one exploits variations of processing time of a program, which depends on the size of the data; the second one the power consumption of the whole system in an immediate way, which depends on treatments made on the data. Study of the McEliece cryptosystem resistance against this type of attacks began only five years ago [START_REF] Strenzke | Side Channels in the McEliece PKC[END_REF]. There exists several side-channel attacks and associated countermeasures [2,[START_REF] Molter | A simple power analysis attack on a McEliece cryptoprocessor[END_REF][START_REF] Shoufan | A Timing Attack against Patterson Algorithm in the McEliece PKC[END_REF][START_REF] Strenzke | A Timing Attack against the secret Permutation in the McEliece PKC[END_REF][START_REF] Strenzke | Timing Attacks against the Syndrome Inversion in Code-based Cryptosystems[END_REF][START_REF] Strenzke | Side Channels in the McEliece PKC[END_REF], we want to provide an implementation secure faces all those attacks. In order to reach this goal, we implemented the most critical step in a Goppa decoder (the root-finding step). Our aim is to find new vulnerabilities, namely those related to information leakages via power consumption.

Hardware implementation of a vulnerable part of the McEliece decryption algorithm

As explained above, our final objective is to implement cryptographic protocols based on errorcorrecting codes on various supports and then to analyze the behavior of these implementations against side-channel attacks.

First of all, we needed to localize the part of the McEliece algorithm, which could be interesting for the attacker. We believe that one of the most important in the McEliece PKC is the decryption, in which the hardware manipulates secret key. We therefore concentrated our effort to this part of the McEliece PKC, which can be implemented in three ways: using the Patterson [START_REF] Patterson | The Algebraic Decoding of Goppa Codes[END_REF], Berlekamp-Massey [START_REF] Massey | Shift-register synthesis and BCH decoding[END_REF] or extended Euclidean algorithm. We selected the so-called error locator polynomial, which is used during the Goppa code decoding, as the most security critical operation. Our next task was to implement it in hardware, in order to be able to observe later dynamically its power consumption, which is related to confidential data being processed. Once the information leakage point will be found, our next objective will be to propose an efficient countermeasure on either algorithmic or implementation level or both.

Implementation of Galois field multiplier in hardware

The main component of all decryption algorithms is the Galois field multiplier. Therefore, we started our work by designing an efficient architecture of the multiplier block.

Algorithm of multiplication in Galois fields:

We decided to use the next Horner's scheme [START_REF] Horner | A new method of solving numerical equations of all orders, by continuous approximation[END_REF] of multiplication algorithm.

Description: Multiplication of two elements α and β of a finite field F2m = F2[X]/Qm(X). The product denoted r by Galois Multiplier (GM) is in F2m . Parameters: m the degree of the polynomial, Qm(X) = m i=0 qiX i a polynomial of degree m on F2 (Version 2). Inputs: Two elements of F2m represented by two polynomials of degree m -1: α(X) = m-1 i=0 αiX i and β(X) = m-1 i=0 βiX i , where αi, βi ∈ {0, 1}, Qm(X) = m i=0 qiX i a polynomial of degree m on F2 (Version 1) Output: Product of the inputs seen as a polynomial of degree m -1: r(X) = m-1 i=0 riX i .

We can represent all polynomials as vectors, as follows:

α = m-1 i=0 αiX i ⇒ α = (αm-1, . . . , α0). GM: r(X) ← αm-1β(X) For i from m -1 downto 1 do r(X) ← r(X)X + αi-1β(X) + rm-1Qm(X) End for Return r Example: m = 5 Qm(X) = x 5 + x 3 + x 2 + x + 1 α(x) = x 4 + x 2 + x and β(x) = x 4 + x x 5 x 4 x 3 x 2 x 1 Qm 1 0 1 1 1 1 α 1 0 1 1 0 β 1 0 0 1 0 r 1 0 0 1 0 (i = 4) 1 0 0 1 0 0 0 0 0 0 0 1 0 1 1 1 1 r 0 0 1 0 1 1 (i = 3) . . .

Implementation of a parameterized Galois field multiplier:

We implemented two versions of the Galois field multiplier algorithm from the preceding section in VHDL. The size of the multiplier is very flexible, since m is given as a parameter. Version 1 of the multiplier block has the polynomial given as an input (a binary vector specifying coefficients of the polynomial) and the second version (which is less flexible, but also less expensive) has the polynomial given as a second parameter. We implemented both versions on the Evariste II FPGA [START_REF]Evariste II -A Modular Hardware System for Development and Evaluation of Cryptographic Functions and Random Number Generators[END_REF] module featuring Altera Cyclone III device EP3C25F256-C8. Implementation results obtained using Altera Quartus II version 9.1 software are presented in Table 1 in term of logic Altera cells (LCELLs) and frequency (freq). As it can be seen (and expected), the flexible solution using polynomial as an input (Version 1) is always slower and larger. It is therefore important to evaluate if the selected decryption algorithm needs such flexibility or not.

Implementation of error locator polynomial in hardware

We implemented two versions of architectures for error vector computation. The first one computes the polynomial σ from right to left [R2L] and the second one from left to right [L2R] (according to the Horner algorithm [START_REF] Horner | A new method of solving numerical equations of all orders, by continuous approximation[END_REF]). The polynomial has the following forms:

σ(X) = σtX t + σt-1X t-1 + . . . + σ1X + σ0 [R2L] σ(X) = (((σtX + σt-1)X + . . .)X + σ1)X + σ0 [L2R]
Implementation of the module computing error locator polynomial from right to left: The first architecture of the module computing error locator polynomial from right to left [R2L] is presented in Figure 1 and the implementation results are presented in Table 2. As it can be seen, this version needs two GMs and two registers for saving intermediate results. Coefficients σi and error vector are saved in separate embedded memory. The maximum frequency of this block is much lower than that of GM alone, because of memory accesses (coefficients σi are read from the memory on the fly). Implementation of the module computing error locator polynomial from left to right: The second architecture of the module computing error locator polynomial from left to right [L2R] according to Horner's scheme is depicted in Figure 2 and the implementation results are presented in Table 3. As it can be seen, this architecture needs only one GM and one register. Memories for coefficients σi and error vector (the result) are the same as in the previous case. Of course, both architecture versions give the same error vector as a result. However, intermediate results, which are saved in registers, are different. Therefore, both circuits will leak information differently. We expect that this difference can be used for constructing some efficient countermeasure in the future.

Figure 1 :

 1 Figure 1: Architecture for [R2L] evaluation

Figure 2 :

 2 Figure 2: Architecture for [R2L] evaluation

Table 1 :

 1 Implementation results for Version 1 (with variable Q m ) vs Version 2 (with fixed Q m )

	Version 1 (with variable Q m )	Version 2 (with fixed Q m )
	m LCELLs	freq (MHz)	m LCELLs freq (MHz)
	5	35	348	5	21	453
	6	54	275	6	29	405
	7	68	207	7	37	455
	8	90	175	8	56	355
	9	117	172	9	59	414
	10	140	135	10	83	321
	11	178	140	11	87	411
	12	209	118	12	123	276

Table 2 :

 2 Implementation results for [R2L] evaluation

	m	t	LCELLs freq (MHz) time (µs)
	5	3	287	120	1.07
	6	5	311	118	3.25
	7 10	341	95	14.82
	11 50	511	78	1339.08

Table 3 :

 3 Implementation results for [L2R] evaluation

	m	t	LCELLs freq (MHz) time (µs)
	5	3	274	110	1.16
	6	5	293	110	3.49
	7 10	311	107	13.16
	11 50	414	85	1228.80

Principle. The principle of side-channel attacks appeared only in 1996[START_REF] Kocher | Timing attacks on implementations of Diffie-Hellman, RSA, DSS, and other systems[END_REF]. A side-channel attack is an attack which exploits the laws of the physics to obtain some information contained in channels associated to an implementation, a circuit. The purpose is to extract secrets manipulated by smart cards or cryptographic components. The side-channel does not aim at transmitting voluntarily some information. It leaks information that an observer can interpret. On the other hand, algorithmic or physic countermeasures can be proposed to stop such a leakage. Implementation of a cryptosystem is a tradeoff between security and efficiency. That is why cryptanalysis and implementation must be simultaneously considered.


Sarah Foley

Angela Staicu

Alexandru Pascu

Mironel Enescu

Photolysis of Bovine Serum Albumin by near UV irradiation

Bovine Serum Albumin (BSA) is the most abundant protein in blood plasma and serves a depot as well as a transport protein for numerous compounds. The structure of BSA is for the most part helical (70%). BSA is also a protein characterised by its high pecentage of lysine residues (10 % of its amino acid content). It is well known that side chains containing Tyrosine and Tryptophan play a dominant role in the near UV electronic absorption and subsequent fluorescence of proteins. However, proteins rich in lysine residues possess an unusual spectral absorption tail between 300 -400 nm [1]. This band can play an important role in understanding protein damage induced by near UV irradiation. We report here the photolysis of BSA due to irradiation at 355 nm and studied using Raman spectroscopy which allows the spectral changes induced by photolysis to be directly identified.

Introduction

[1] L. Homchaudhuri and R. Swaminathan, Bull. Chem. Soc. Jpn., 77 (2004), 765-769.

[2] S.U. Sane, S.M. Cramer and T.M. Przybycien, Anal. Biochem. 269 (1999), 255-272.

[3] C.David, S. Foley, C. Mavon and M. Enescu, Biopolymers, 89 (2008), 623-634.

The photolysis of BSA at 355 nm is followed by a partial reduction fo the disulphide bridges. Disulphide bridge reduction gives rise to a modification of the secondary structure of the protein due to unfolding. The photolysis of BSA also brings about changes in the tertiary structure of the protein.

Conclusions

I. Excitation and Emission Spectra

Fig. 3 shows the emission spectra obtained when a solution of BSA is excited at 400 nm. Emission shows a maxima at 463 nm which is slightly red shifted and increases in intensity after irradiation. Fig. 4 shows the excitation spectra of BSA and highlights two species with maxima situated at 360 nm and 306 nm. Following irradiation, the peak at 360 nm is red shifted and decreases in intensity whilst the peak at 306 disappears, thus indicating the involvement of both species in the photolysis of BSA.

Results & Discussion

II. Raman Spectra

The Raman spectrum for native BSA between 380 and 1750 cm -1 is presented in figure 5a. Figures 5b and5c show the Raman spectra of BSA after irradiation with 355 nm light under aerobic and anaerobic conditions.

The absence of any changes between spectra 5b and 5c suggests that the changes observed following the photolysis of BSA are not mediated by singlet molecular oxygen.

Changes are observed in three regions of the spectra :

1 -Between 480 -600 cm -1 (S-S band) 2 -Between 1200 -1400 cm -1 3 -Between 1520 -1700 cm -1 (Amide I band)

Raman Measurements:

Raman spectra for all solutions were measured in the region of 380 -1750 cm -1 using the second harmonic of an Nd/YAG laser 532 nm with a pulse duration of 5 ns. Solutions were irradiated using 355 nm with a fluence of 10 mW / mm 2 for 30 or 60 minutes. All solutions were measured at pH = 7 using a concentration of 1 mM BSA and the reference solution used throughout the measurements was H 2 O. The wavenumbers of Raman bands were reproducible within ±1 cm -1 .

Emission Spectra:

Fluorescence measurements were made using a Horiba Jobin Yvon Fluorolog ® -3-22 flurometer.

The instrument settings were as follows: Emission slits = 5.0 nm Excitation slits = 5.0 nm

III. Structural Modifications of the Protein

The irradiation of BSA brought about some structural changes which manifested themselves in a number of ways.

Changes in the S-S band

Photolysis of BSA results in significant changes in the Raman S-S band. During the reaction the intensity of the band falls with components 2 and 3 showing a marked decrease (Figure 6). The evolution of the S-S bands under irradiation strongly suggests that the S-S bridges are partially reduced.

Changes in the Amide I band

Analysis of the Raman band provides information about the protein secondary structure and can be analysed using the following attribution [2]. Peak 1 at 1640 cm -1 -short segment chains Peak 2 at 1657 cm -1 -the ordered α-helix Peak 3 at 1671 cm -1 -the unordered α-helix Peak 4 at 1684 cm -1 -vibrations of various turns Following the photolysis of BSA after irradiation at 355 nm the cleavage of the disulphide bridges is associated with a significant evolution of the Amide I band (Table 1 andFigure 7). Changes in the secondary structure suggesting protein unfolding can be summarised as: decrease in the fraction of the ordred α-helix increase in the fraction of the unordered α-helix, protein unfolding.

Fluorescence Measurements

The photolysis of BSA was also followed via the fluorescence spectra of the tryptophan residues.

-red shift in the intrinsic tryptophan fluorescence of BSA -in the intensity of emission Indicative of a change in the tertiary structure of the protein.

Mechanism

The mechanism for the photolysis of BSA at 355 nm is at present unknown, however one can consider the following :

Acknowledgements:

  This work was supported by ANCS through the Brancusi French-Romanian bilateral program (Contract Number 207/2009 ) and by CNMP (Contract Number 61-023/2007).

Fig. 3 .

 3 Fig. 3. Emission spectra of native BSA and photolysed BSA, λ ex = 400 nm.

Fig. 4 .

 4 Fig. 4. Excitation spectra of native BSA and irradiated BSA: a) λ em = 460 nm, b) λ em = 400 nm

Fig. 5 .Fig. 2 .

 52 Fig. 5. Raman spectra of BSA and irraditaed BSA under aerobic and anerobic conditions.

Fig. 6 .

 6 Fig. 6. Analysis of the S-S Raman band for 1 mM BSA as a function of irradiation time.

Fig. 7 .

 7 Fig. 7. Analysis of the Amide I Raman band for 1 mM BSA as a function of irradiation time.

Fig. 8 .

 8 Fig. 8. Emission spectra of native BSA and irradiated BSA, λ ex = 400 nm.

  the chromophore absorbing at 355 nm and by which the reduction of the disulphide bridges results in the unfolding of the protein as previously described [3].

Irradiation time / min Peak 1 Peak 2 Peak 3 Peak 4

	0	29.84	36.77	24.96	8.43
	30	27.19	34.95	27.17	10.69
	60	24.84	26.50	35.69	12.97

Table 1 .

 1 Analysis of the Amide I band showing the % of each peak as a function of irradiation time

	ϕ
	hω
	Fig. 1. Ribbon diagram
	of BSA


Catherine Berho

Nicolas Mazzella

Dominique Munaron

Cecile Miege

AQUAREF intercomparison passive sampling exercise: monitoring of pesticides in surface water

  AQUAREF intercomparison passive sampling exercise: monitoring of pesticides in surface water C. Berho, N. Mazzella, D.Munaron, C Miège Other co authors : B. Lepot, H. Budzinski, C. Gonzalez , A. Gravell, R. Gunold, E. Peinerud, F. Smedes, B. Vrana The Water Framework Directive (WFD) is probably the most significant legislation in the water field that has been introduced for many years. Several publications pointed out that passive samplers (PSs) can be valuable tools as complementary method in monitoring water quality in the context of the WFD However, there is still a lack of quality assurance and control (QA/QC) procedures that demonstrate the reliability and the comparability of results obtained by passive sampling. Up to date, very few in situ intercomparison exercises on PSs have been performed until now. To our knowledge only three pan-European or international intercomparison studies have been organized in situ since 2007 (ICES Passive sampling survey and intercalibration, AQUAREF, IPSIC and NORMAN). Among them, the AQUAREF intercomparison exercise was conducted on several passive samplers for priority pollutants (PAH, pesticides and metals) in surface and coastal waters. Presentation of this exercise and first outputs have been already published (Miège et al., 2012) and showed that in spite of variety of expert laboratories, strategies and tools, we observed a low and satisfactory uncertainty on the estimation of mean TWA concentrations with passive sampling. The objective of this poster is to present specific results obtained during this interlaboratory exercise for pesticide sampling. Two trials were conducted during 14 days in the Charente River (at Beillant, France) and in a marine site, the Thau Lagoon (Hérault, France) with the participation of 11 laboratories at Beillant and 8 laboratories at Thau. Nine pesticides (acetochlor, S-metolachlor, 5 WFD priority substances, i.e., alachlor, atrazine, diuron, isoproturon, simazine and two atrazine metabolites, i.e. desethylatrazine, deisopropylatrazine) were considered. Among the tested PSs, we found Polar Organic Chemical Integrative Sampler (POCIS, pharmaceutical & pesticide configuration), polar Chemcatcher (with SDB-XC, SDB-RPS or C18 phases), silicon rod and silicon sheet, speedisk hydrophylic DVB06, ... We will present here and discuss results on i) Comparison of pesticide quantities accumulated per surface of passive sampler exposed ii) Mean and variability of Time Weighted Average (TWA) water concentrations of pesticides obtained with passive samplers and comparison with water concentration from grab samples iii) Pesticide fingerprints in the various passive samplers compared with those in grab samples


Patrice Loisel

Jean François Dhôte

A tree-growth model to optimize silviculture

Keywords: growth model, optimization, control

In this paper, we present the description of a simplified model of the dynamic of a mono-specific even-aged forest. The model studied is a tree-growth model based on a system of two ordinary differential equations concerning the tree basal area and the number of trees. The analytical study of this model permits us to predict the behavior of the system solutions. We are trying to highlight the influence of economic parameters and growth parameters on the system solutions, in the framework of the optimization of silviculture.

Introduction

The forest management, because of its impact on our environment, is a topic that now involves researchers of many disciplines : forestry, economy, and ecology. These various communities have models adapted to the questions they wish to tackle. The economists usually study the best age at which to cut down a tree or stand of trees, the simultaneous management of several forest stands. As to foresters, they are moreover interested in silviculture at the stand level. We will focus on the models developed by foresters.

The models of growth for silviculture, expanded rapidly these last few years, and represent a significant part of the developed models. Here, we are focusing on a particular type of forest (thus a particular type of model) : a mono-specific even-aged forest where all trees belong to the same species and are the same age. Models built by forest modellers are based on statistical adjustments of dendrometric data [START_REF] Porte | Modelling mixed forest growth : a review of models for forest management[END_REF] : these models accurately describe the evolution of a forest, but the analytical study of those models is made difficult due to their complexity. The analytical study of models allows us to predict the influence of different parameters : for instance, if by modifying parameters to take into account the climate change and analyzing so the potential consequences is provide. To allow analytical studies while having realistic model, we decide in this paper to consider a simplified model.

The models, on which we are focusing here, are tree centered distance independent models where the trees are not spatialized. In this type of model each tree is characterized by its basal area at the height of 1.3 meters : s and eventually by its height h. The model here described is based on the concepts developed in the growth model "Fagacées" [START_REF] Dhôte | Forme de la tige, tarifs de cubage et ventilation de la production en volume chez le Chêne sessile[END_REF] [START_REF] Moguedec | Fagacées a tree-centered growth and yield model for Sessile Oak (Quercus petraea L.) and common Beech (Fagus sylvatica L.). Accepted with revisions in Annals of Forest Science[END_REF] for Oak or Beech forest. In this model, the link between the stand level and the individual tree level is explicit. This modeling allows us to describe the evolution of a forest of high density. "Fagacées" was broadcast through the project Capsis [START_REF] De Coligny | Capsis : Computer-Aides Projection for Strategies In Silviculture, a software platform for forestry modellers[END_REF].

In order to allow analytical studies, we are starting with the simplified hypothesis in which we consider that all trees have the same basal area. We then consider a forest of n trees with basal area s which needs management such as thinning e throughout time.

The tree growth (due to the observed densities) is not independent of its neighbor's growth. There is a competition for the available resources : photosynthesis and access to the light on one hand, and mineral nutrients on the other hand. Thus, cutting a tree implies the increase of its neighbors' growth, and cutting no tree limits individual growth. That shows how a forest is not just a juxtaposition of trees. This phenomenon is considered through an assessment equation that allows us to distribute the energy resources of the forest stand level between the various trees. This equation coupled to an equation describing the evolution of the number of trees leads us to a dynamic system of the forest. The studied model takes into account the characteristics that the foresters consider as the most important and as required : the basal area at 1.3 meters, the height, and the number of trees.

In Section 2 we will present the designing of the model. In Section 3 we will present general results on the behavior of the system solutions, then we will look for strategies which permit to leave the viability domain in minimal or maximal time. Finally in Section 4 we will highlight the influence of economic parameters and growth parameters for silviculture within the resolution of an optimization problem.

Designing the model

The trees density and the RDI of forest stand.

Let's consider a forest with a given area. It is intuitively clear that the tree number which this area can bear is limited. The environmental conditions (type of soil, local climatic conditions) are also factors to be taken into account for the maximum tree capacity. Foresters have established a law called "self-thinning", described hereafter, to evaluate this maximum capacity. Let's note s * the average tree basal area (at the height of 1.3 meters) of the forest.

Reineke [START_REF] Reineke | Perfecting a stand-density index for even-aged forest[END_REF] observed monospecific forests with various densities and various species. Out of these observations he claimed the maximum tree number n max (s * ) that a stand can bear is given by the following self thinning relation :

log n max (s * ) = C 0 - q 2 log s * ,
where C 0 > 0 and 1 < q < 2 are characteristic constant values of the forest species and of its environment, and in particular the ground fertility.

As for a given s * , beyond the number n max (s * ) the trees die, the forest stand (in terms of tree number) has to remain under this limit. To simplify we will here make the assumption that all the trees have the same basal area s. For a forest of where the effective tree number is n, taking into account the relation of self-thinning, the density r is the ratio of the tree number and the maximum tree number of basal area s that the forest can support, r is then defined by :

r(n, s) = n n max (s)
is written this way :

r(n, s) = ns q 2 e C 0 = Ans q 2
where A := e -C 0 . This ratio r is called RDI (Relative Density Index or Reineke Density Index). By definition this ratio is always less than 1.

Competition between trees : from forest stand level to individual tree level

We will now describe the temporal evolution of state variables s, n and now r, of the considered forest.

The growth of a tree depends on its neighbors. There is competition for the resources and the death of a tree, natural or due to cuttings, implies an increased growth for its neighbors . We make the assumption, that in the course of time, silviculture makes it possible to maintain the trees uniformly distributed on the area. The model is characterized by the existence of two levels in the modeling.

At the forest stand level the available energy for the considered forest, is considered globally, due to photosynthesis or due to nutrients in the soil. This supplied energy makes it possible to ensure at the same time the maintenance and the growth of the trees. The share reserved for maintenance increases with the tree height, therefore with time, which limits all the more so the available part for growth. The energy left for growth is therefore a decreasing time function and allows the increase of basal area of the forest stand. The increase of basal area of the forest stand at its peak of density r(n(t), s(t)) = 1 is given by the function V (.). We assume that V (t) verify the following properties :

(H 1 ) : V (.) is a positive, decreasing, convex function of t.
For a lower density (r(n(t), s(t)) < 1), the effective increase of the basal area for the forest is reduced by a factor dependent on this same density : g(r(n(t), s(t))) at any time. Thus the energy actually used at time t is given by :

g(r(n(t), s(t)))V (t).
The function g(.) is supposed to satisfy the following properties :

(H 2 ) : g(.) is an increasing, concave function of r such that g(r) > r for r ∈ (0, 1), g(0) = 0, g(1) = 1.

The concavity of g is related to crown development in relation to basal area.

On the individual tree level, tree growth is characterized by the evolution of tree basal area and therefore by the evolution of the function s(t) : the instantaneous increase is thus ds(t) dt . As mentioned in the hypothesis all trees have the same basal area, the total sum increase of basal areas of all trees is n(t) ds(t) dt . This total increase is obtained from the available energy resources. We thus obtain the equation which describes the link between the forest stand level and the individual tree level :

g(r(n(t), s(t)))V (t) = n(t) ds(t) dt for n(t) > 0, for all t
For any t, this enables us to establish the first dynamic equation of our model :

ds(t) dt = g(r(n(t), s(t))) n(t) V (t)
In addition, the evolution of the tree number depends on several factors. To permit analytical study of the model, we've decided to simplify and we suppose the only cause of tree mortality is due to fallings that foresters could operate. We noted e(t) the instantaneous rate of trees cutting at time t. Thus the evolution of the tree number is given by :

dn(t) dt = -e(t).
We wish to preserve a minimum tree number in the forest stand, which implies n(t) ≥ n > 0, for any t. Technologically and to ensure a provisioning not too irregular, the thinning rate is limited : 0 ≤ e(t) ≤ e, for any t.

The forest is therefore described using the two state variables s, n and its evolution follows the following dynamic :

(S 0 )      ds(t) dt = g(r(n(t), s(t))) n(t) V (t) dn(t) dt = -e(t)
with the constraints 0 ≤ e(t) ≤ e, n(t) ≥ n, r(n(t), s(t)) = An(t)s(t) q 2 ≤ 1 for any t.

Foresters built this type of model from observed forest data. The available data only allows us to validate the model on a limited period of time. The system (S 0 ) has therefore a time limit domain :

t ∈ [0, T * ].
It is a dynamic system in the state variables n and s, controlled by the control variable e. For a cutting policy, i.e. the data of a particular function e(.), and for each initial condition (s(0), n(0)), this system has a single solution : we will suppose the functions g(.), V (.) are regular enough for it to happen. We will specify these trajectories in the following paragraph.

To finish with the model description, the tree height h is supposed to depend only on the tree basal area s and on the dominant height h 0 (average height of the 100 largest trees), h 0 is a concave function of time t and shouldn't depend on silviculture (cuttings in the course of time) and thus depends only on time t. The height h has therefore no influence on s(t) and its evolution, h is consequently an output of the model.

If, as supposed earlier, the basal area s at time t is the same for all trees, the height h is also the same. We therefore deduce h(t) = h 0 (t), for any t.

3 Studiing the solutions

Model properties

The solutions of the dynamic system (S 0 ) must satisfy in particular the constraint r(n(t), s(t)) ≤ 1 for any t. If there is no cutting, i.e. if e(t) = 0 for any t, we deduce that n(t) = n(0), s(.) and r(.) are increasing with respect to time t. Let's suppose there is one time τ < T * such as r(n(τ ), s(τ )) = 1, we deduced t > τ if we apply a control identically null then r(n(t), s(t)) > 1 and the constraint is no longer satisfied. In order to let us know which control we should apply we are led to study the evolution of the density function r(.) :

dr(n(t), s(t)) dt = r ′ s (n(t), s(t)) g(r(n(t), s(t))) n(t) V (t) -r ′ n (n(t), s(t))e(t) = r(n(t), s(t)) n(t) [ q 2 g(r(n(t), s(t))) s(t) V (t) -e(t)] (1)
Out of this last equation we can deduce that in order to respect the constraint r(n(t), s(t)) ≤ 1 for t > τ , we should apply a non-identically null control on the system. Thus the cutting e(t) = q 2

V (t) s(t)
for t > τ respects the constraints by binding, i.e. r(n(t), s(t)) = 1. If we define the function e r (., .) by : e r (s, t) := q 2

V (t) s , for any s > 0, t > 0, the solutions, independently of the cutting function e(.) applied to the system, are only valid if the constraint : e r (s(t), t) ≤ e is satisfied. We are therefore led to formulate the following assumption (H 3 ) :

(H 3 ) : e r (s m (t), t) = q 2 V (t) s m (t) < e for all t ∈ (0, T * )
where s m (t) is the minimal value s(t) can reach at the time t.

Remark 3.1 s m (t) is not specified at the moment but will be specified later on, however we can take an approximate lower bound for now : s m (t) > s(0).

We noted previously that the system (S 0 ) is considered only for t ∈ [0, T * ]. It is advisable to specify now, the behavior of the solutions in this interval. The following Lemma shows us that the system validity field depends on this energy value : Lemma 3.1 Assuming (H 2 ), (H 3 ), then :

(i) if V(0; T * ) is large enough then there exists a time τ < T * such that r(n(τ ), s(τ )) = 1 and n(τ ) = n. The dynamical system is only valid on the interval [0, τ ]. This time τ depends on the evolution of the cutting e(.).

(ii) conversely if V(0; T * ) is small enough then the dynamical system is valid throughout the entire interval [0, T * ].

Proof : (i) From g(r) ≥ r we deduce :

ds(t) dt ≥ r(n(t), s(t)) n(t) V (t) = As(t) q 2 V (t) hence : s(T * ) 1-q 2 ≥ s(0) 1-q 2 + A(1 - q 2 )V(0; T * ) From r(n(T * ), s(T * )) ≤ 1 and n(T * ) ≥ n we deduce : s(T * ) ≤ 1 (An) 2 q
. If V(0; T * ) is large enough, we obtain a contradiction.

(ii) Let's set τ the first period where r(n(t), s(t)) reaches 1, then for 0 < t ≤ τ we deduce : Specific trajectories easily expressed in terms of control, will play an important role, we are introducing them here : let's consider the system of equations (S 0 ), for trajectory E 0 from a fixed initial condition (s(0), n(0)) we apply the maximum cutting e(t) = e until we reach the value n for the tree number, t 0,n is the time needed to go from the tree number n(0) to n. By definition, we therefore have t 0,n = n(0)n e .

ds(t) dt ≤ V (t) n(t) ≤ V (
For trajectory E 0 , starting from the same initial condition (with r(n(0), s(0)) < 1) we apply the minimum cutting e(t) = 0 until reaching the value 1 for the RDI r, then we apply the control e r (s(t), t) until n = n. t 0 is the time needed to go for the RDI from r(n(0), s(0)) to 1 and T 0 the final time. By definition, we therefore have t 0 and T 0 respectively solutions of :

q 2 n(0) 2 q -1 A 2 q V(0; t 0 ) = 1 r(n(0),s(0)) u 2 q -1 g(u) du (at constant n) n 1-2 q = n(0) 1-2 q + A 2 q (1 - q 2 )V(t 0 ; T 0 ) (at constant r)
Notations. To summarize we note the following definitions of the specific trajectories :

E 0 : e(t) = e if t < t 0,n 0 if t > t 0,n E 0 : e(t) = 0 if r(n(t), s(t)) < 1, i.e. t < t 0 e r (s(t), t) if t 0 < t < T 0
For t 0,n < T < t 0 we can also define an intermediate trajectory

E T E T : e(t) = 0 if t < T -t 0,n e if t > T -t 0,n if t 0,n < T < t 0,n + t 0 E T : e(t) =      0 if t < t 0 e r (s(t), t) if t 0 < t < t * e if t * < t < T if t 0,n + t 0 < T < T 0
where t * is defined by : n(t * ) 1-2 q = n(0) 1-q 2 + A(1 -2 q )V(t 0 ; t * ) and (Tt * )e = n(t * )n.

We note that E T 0 = E 0 and by extension if T < t 0,n then E T = E 0 .

Figure 1 : Phase plane in the coordinates s and n.

The functions obtained by just following the trajectories E 0 , E T and E 0 will be noted by the indices 0 , T and 0 .

From the increasing of the basal area s and the non increasing of the number n we deduce that the system solutions have no choice but to move to the right bottom in the phase plane. The Lemma 3.1 (i) has shown that, if V(0; T * ) is large enough, the solutions are not valid throughout the entire interval [0, T * ]. That implies that as from a time τ , the solution doesn't belong to the validity domain defined by the constraints r(n(t), s(t)) ≤ 1 and n(t) ≥ n. The only point which makes it possible to leave this validity domain is the point such as r(n(τ ), s(τ )) = 1 and n(τ ) = n, this point is represented by a square on Figure 1. As the solution remains valid basal area s(t) verifies for all t < T : s(t) ≤ s = 1 (An)

2 q
(is deduced from r(n(t), s(t)) ≤ 1).

We define T (resp. T ) as the minimum (resp. maximum) time necessary to reach the point defined by r(n(T ), s(T )) = 1 and n(T ) = n. Then :

-if T ≤ T the solution remains valid whatever the trajectory (i.e. whatever the evolution of the cutting e(.)).

-if T < T ≤ T the system has a solution on [0, T ] for certain controls e(.).

-if T > T the system has no solution on [0, T ] whatever the controls e(.).

We noted that, from Lemma 3.1, if V(0; T * ) is small enough then T and especially T can no exist.

A particular case

We could consider the particular function g

(r) = g θ (r) = r 1-θ , 0 < θ < 1. In that case ds(t) dt = A 1-θ s(t) q 2 (1-θ) n(t) θ V (t)
. The basal area s is explicitly deduced from the tree number n :

s(t) 1-q 2 (1-θ) = s(0) 1-q 2 (1-θ) + A 1-θ (1 - q 2 (1 -θ)) t 0 V (u) n(u) θ du.
In this class of functions g θ (.) we consider the extreme case (θ = 0) for which some of the properties of the hypothesis (H 2 ) are not satisfied : g 0 (r) = r. In this last case G(r) ≡ 0 and the evolution of the basal area s is independent from the evolution of the tree number n :

s(t) 1-q 2 = s(0) 1-q 2 + A(1 - q 2 )V(0; t).
In this particular case, provided that V(0; T * ) is large enough, T and T are equal, don't depend on the cutting e(.) and are the unique solution of the following equation in T :

s(0) 1-q 2 + A(1 - q 2 )V(0; T ) = s 1-q 2 .
3.2 Minimum and maximum time necessary to reach the point (r, n) = (1, n)

In order to succeed in the conclusion of the study of the minimum and maximum time needed to reach this point, we will need the following properties and definitions related to the function g(.). The increase in basal area of each tree is given by g(r(n(t), s(t))) n(t) V (t). We will to know thereafter the evolution of g(r(n(t), s(t))) n(t) , for the same aim, we will need to define the functions G(.), γ(.) :

Definition 3.2 The function G(.) is defined by : G(r) = d dr [ r g(r) ] = g(r) -rg ′ (r) g 2 (

r)

. The function γ(.) is defined by : γ(r) = rg ′ (r) g(r) .

In the "Fagacées" model,

g(r) = (1 + p)r r + p with p > 0, G is constant G(r) ≡ G = 1 1 + p .
From G(.) and γ(.) definitions, we can establish the following properties for the model : Lemma 3.2 Assuming the hypothesis (H 2 ), then :

(i) The function r g(r)
is an increasing function of r and G(.) satisfies 0 < G(r)g(r) ≤ 1 for any r > 0.

(ii) The function g(r(n, s)) n is a decreasing function of n.

(iii) The function g(r(n, s)) is an increasing function of s.

(iv) The function γ(.) satisfies γ(r) ≤ 1 for any r > 0.

(v) if n(t) and s(t) are solutions of systems (S 0 ) the function g(r(n(t), s(t))) n(t) is an increasing function of t.

Proof : (i) From the concavity of g(.), d dr [g(r)rg ′ (r)] = -rg ′′ (r) > 0 for any r > 0 and from g(0) = 0, we deduce that g(r)rg ′ (r) > 0 and G(r

) > 0. From g ′ (r) > 0, G(r) ≤ 1 g(r) . (ii) ∂ ∂n [ g(r(n, s)) n ] = r(n, s)g ′ (r(n, s)) -g(r(n, s)) n 2 = -G(r(n, s)) g 2 (r(n, s)) n 2 < 0. (iii) ∂g(r(n, s)) ∂s = g ′ (r(n, s))r ′ s (n, s) > 0.
(iv) From g(r)rg ′ (r) > 0 we deduce the result.

(v) d dt [ g(r(n(t), s(t))) n(t) ] = (g ′ (r)r ′ s )(n(t), s(t))) n(t) ds(t) dt + G(r(n(t), s(t))) g 2 (r(n(t), s(t))) n 2 (t) e(t)
and from (i) we deduce the result.

We are thus focusing on the trajectories and also on the strategies which allow us to reach respectively in a minimum and maximum time the point (1, n) in the (r, n) coordinates.

The minimum time T (resp. the maximum time T ) is reached by solving the problem of optimal control : min and s(.) are the state variables governed by the system (S 0 ) of initial condition (n(0), s(0)) and satisfying constraints, for all t ∈ [0, T ), r(n(t), s(t)) ≤ 1, n(t) ≥ n and the right end time constraints r(n(T ), s(T )) = 1, n(T ) = n. Proposition 3.1 Assume (H 2 ), (H 3 ). If n(.) and s(.) are the solutions of the system (S ′ ) for a control e(.) then, ∀t ∈ [0, T ] :

(i) n 0 (t) ≤ n(t) ≤ n 0 (t) (ii) s(t) ≥ s 0 (t) (iii) if g(r) = g θ (r) = r 1-θ , s(t) ≤ s 0 (t)
If we assume that the final tree-number n(T ) is equal to n then :

(iv) n(t) ≤ n T (t) (v) if g(r) = g θ (r) = r 1-θ , s(t) ≤ s T (t)
Proof (i) and (iv) Follows from the definition.

(ii) From Lemma 3.2 (ii) we deduce :

ds(t) dt = g(r(n(t), s(t)) n(t) V (t) ≥ g(r(n 0 (t), s(t)) n 0 (t) V (t).
For t ≤ t 0 we deduce : ds g(r(n(0), s

q 2 )) ≥ V (t) n(0)
dt then by integration of the inequality :

4 Optimization of silviculture

In order to optimize the silviculture, we are interested in problems which consist in seeking the minimal and maximum values of a variable function depending on the state variables n and s.

Preliminary results

We consider the hypothesis (H 4 ) :

(H 4 ) : there exists a constant γ > 0 such as γ ≤ γ(r) for any r ∈ (0, 1).

We obtain the following Lemma (with proof in Annex A) : Lemma 4.1 Assume (H 2 ), (H 3 ). If n(.) and s(.) are the solutions of the system (S ′ ) for a control e(.) then, ∀t ∈ [0, T ] : (with the convention that the inequalities including s 0 (t) are valid only if g(r) = g θ (r) = r 1-θ ) (i) the function g(r(n, s)) n satisfies :

g(r(n 0 (t), s 0 (t))) n 0 (t) ≤ g(r(n(t), s(t))) n(t) ≤ g(r(n 0 (t), s 0 (t))) n 0 (t) .
Morever, assume (H 4 ) :

(iia) if 0 < b < 1 -q 2 γ 1 -γ then n 0 (t)s 0 (t) b ≤ n(t)s(t) b ≤ n 0 (t)s 0 (t) b
in particular, for b = q 2 the RDI r(n, s) satisfies : r(n 0 (t), s 0 (t)) ≤ r(n(t), s(t)) ≤ r(n 0 (t), s 0 (t))

(iib) if b > b * = 1 + q 2 ( 1 g(r(n,s(0))) -γ) 1 -γ then n 0 (t)s 0 (t) b ≤ n(t)s(t) b ≤ n 0 (t)s 0 (t) b (iii) the relative increase ξ of the basal area s satisfies ξ m (t) = s ′0 (t) s 0 (t) q 2 s 1-q 2 ≤ ξ(t). Moreover, if g(r) = g θ (r) = r 1-θ , ξ m (t) = s ′0 (t) s 0 (t) q 2 s 0 (t) 1-q 2 ≤ ξ(t).
Remark : If we assume that the final tree-number n(T ) is equal to n then the result obtained in the Lemma 4.1 remains valid if we replace respectively n 0 , s 0 by n T , s T .

The optimization problem

We are focusing here, on the setting in the wood market of a forest whose evolution is set by the model studied in the previous paragraphs. We are introducing the price (minus the cost of thinning) which depends only on the basal area s and the height h : we noted P 0 (s, h, t). Owing to the fact that the height h does not depend on the basal area s and is equal to a fixed function h 0 of time t, the price can be written in a new function P of s and t : P (s(t), t) := P 0 (s(t), h(t), t).

In other words we will set the price function in the following form : P 0 (s, h, t) = p(s)he -δt where p(.) is a price function for the basal area s and δ is the actualisation parameter. We deduce P (s, t) = p(s)h 0 (t)e -δt and if we define the function δ h (.) by : δ h (t) = δ -h ′ 0 (t) h 0 (t) for any t > 0 then P ′ t (s, t) = -δ h (t)P (s, t). Mostly to simplify we'll assume p(s) = ks α , α > 0.

We are assuming that at each time t a quantity e(t) is taken and that at the end of the period of exploitation T the remaining trees would have been cut. The instantaneous value of the trees that would have been cut is P (s(t), t)e(t) and the final value is P (s(T ), T )n(T ).

The criterion which we suggest to maximize consists of an integral term corresponding to the cuttings that would have occurred during the interval [0, T ] and the final term corresponding to the final cuttings at time T .

The optimization problem, relating to the cuttings e(.), on the interval [0, T ], is therefore written :

(P) : max e(.) T 0 P (s(t), t)e(t)dt + P (s(T ), T )n(T )

with 0 ≤ e(t) ≤ e and n and s solutions of (S 0 ) with initial conditions (n(0), s(0)) and fulfilling the constraints : n(t) ≥ n et r(n(t), s(t)) ≤ 1.

Intuitively, from the fact that the function g(r(n, s)) n is a decreasing function of n (Lemma 3.2 (ii)), we are tempted to suggest the following assertion :

In order for the trees to get the best benefits from the nutrients, one should, from the beginning cut a significative number of trees, so that in the end of the exploitation timescale, one should get Corollary 4.1 If we assume that the final tree-number n(T ) is equal to n then : if α < 1 -q 2 γ 1γ and δ h (t) ≤ αγξ m (t), then the optimal trajectory is E T .

The condition on α in Proposition 4.1 (i) implies p must be sufficiently convex. Thus, under the conditions mentioned in the Proposition 4.1 (i) (if p is sufficiently convex and the parameter of actualization not too high), one may find it beneficial to cut the maximum tree number at the beginning to ensure a high rate for the remaining tree basal area at the end of T as foretold in the stated assertion. Similar results were obtained with the full model "Fagacées" [START_REF] Moguedec | Formalisation of silvicultural schedule for optimisation purpose : an application to Sessile Oak (Quercus petraea L.) using the growth model Fagacées[END_REF]. Moreover, in the studied cases in Proposition 4.1, silviculture, i.e. cuttings policy e(.), does not depend on final time T . The conditions depend on economic parameters : a sufficient convexity of the price function relative to the basal area s and a small enough parameter of actualization. The stated assertion however is no longer satisfied under the conditions of the Proposition 4.1 (ii).

Conclusion

In that article, starting with a tree-growth model governed by the tree basal area and the number of trees, we study the viability properties of the system solutions. We highlighted the importance of the economic parameters and growth parameters on silviculture. Hence for a price (minus the thinning costs) is sufficently convex and a parameter of actualization not too high, it is optimal to cut the trees at the beginning of the period of exploitation.

6 Proof of Lemma 4.1 (i) the result is a consequence of Lemma 3.2 (ii) (iii) and Proposition 3.1 (ii) (iii).

(iia) d(n(t)s b (t)) a dt = abn(t) a-1 s(t) ab-1 g(r(n(t), s(t)))V (t)ae(t)n(t) a-1 s(t) ab .

If we denote y(n, s) the expression of d(ns b ) a dt , then, if we assume 0 < a < 1 : y ′ n =an(t) a-2 s(t) ab-1 (b[rg ′ (r) + (a -1)g(r)](n(t), s(t))V (t) + (1a)e(t)s(t)) ≥abn(t) a-2 s(t) ab-1 (γ + a -1)g(r(n(t), s(t)))V (t) y ′ s =abn(t) a-1 s(t) ab-2 ([ q 2 rg ′ (r) -(1ab)g(r)](n(t), s(t))V (t)ae(t)s(t))

≤abn(t) a-1 s(t) ab-2 ( q 2 γ + ab -1)g(r(n(t), s(t))))V (t)

hence, if we choose a such that 1γ < a < min( 1 -q 2 γ b , 1), we deduce q 2 γ + ab -1 < 0 < γ + a -1 then y ′ n > 0 and y ′ s < 0 and : (iib) From the expression of y ′ n and y ′ s and using e(t) ≤ e r (s(t), t), we deduce that, if a < 1-γ : y ′ n ≤an(t) a-2 s(t) ab-1 (b(γ + a -1)g(r(n(t), s(t)))

+ q 2 (1 -a))V (t)
y ′ s ≥abn(t) a-1 s(t) ab-2 (( q 2 γ + ab -1)g(r(n(t), s(t))) -q 2 a)V (t)

and, if b > b 1 (a) = q 2

1a 1aγ 1 g(r(n, s(0)) then y ′ n < 0.

if b > b 2 (a) = q 2 1 g(r(n, s(0)) + 1 -q 2 γ a then y ′ s > 0.

To obtain the minimal limit value for b, as b 1 is increasing in a, and b 2 is decreasing in a, we choose the value a such that b 1 (a) = b 2 (a), this value is a * = (1γ)(1 -q 2 γ) 1 + q 2 ( (iii) From Lemma 3.2 (i), g(r) r is a decreasing function of r then :

s ′ (t) s(t) = A g(r(n(t), s(t))) r(n(t), s(t))

V (t) s(t) 1-q 2 ≥ A g(r(n 0 (t), s 0 (t))) r(n 0 (t), s 0 (t))

V (t) s M (t) 1-q 2
. From s M (t) ≤ s we deduce the result, if g(r) = g θ (r) = r 1-θ , s M (t) = s 0 (t).

Definition 3 . 1

 31 The function V(.; .) is defined by V(t; T ) = T t V (u)du and represents the energy that has been available for growth in the period [t, T ].

  t) n and therefore s(τ ) ≤ s(0) + V(0; τ ) n . If V(0; T * ) is small enough, we deduce r(n(τ ), s(τ )) < 1 in contradiction with the assumption.

T

  ) with the set of admissible values for the control variable [0, e]. n(.)

d(n 0

 0 (t)s 0 (t) b ) a dt = y(n 0 (t), s 0 (t)) ≤ d(n(t)s(t) b ) a dt ≤ y(n 0 (t), s 0 (t)) = d(n 0 (t)(s 0 (t)) b ) a dtby integration, we obtain the result.

  γ g(r(n,s(0))γ) hence the result if b > b * = b i (a * ).

s(t) s(0) dx g(r(n(0), x q 2 )) ≥ V(0; t) n(0) = s 0 (t) s(0) dx g(r(n(0), x q 2 )) and we deduce s(t) ≥ s 0 (t).

For t > t 0 , ds(t) dt ≥ V (t) n 0 (t) = ds 0 (t) dt and from s(t 0 ) ≥ s 0 (t 0 ) we deduce by integration s(t) ≥ s 0 (t).

(iii) From the previously stated expression of the basal area s and n(t) ≥ n 0 (t) we deduce :

and hence the result.

(v) From the previously stated expression of the basal area s and n(t) ≤ n T (t) we deduce :

If we remark that to reach the point (r, n) = (1, n) in minimal time (resp. in maximal time) is equivalent to reach s = s in minimal time (resp. in maximal time), we deduce the the trajectory that allows to reach the point (r, n) = (1, n) in minimal or maximal time : Corollary 3.1 Assume (H 2 ), (H 3 ). Let T the minimal time (resp. T the maximal time) necessary to reach the point (r, n) = (1, n) using the control e(.). Then : (i) if T is finite and g(r) = r 1-θ , the trajectory that allows to reach the point (r, n) = (1, n) in minimal time T is the trajectory E 0 .

(ii) if T is finite, the trajectory that allows to reach the point (r, n) = (1, n) in maximal time T is the trajectory E 0 . Maximal time T is then the solution of :

a limited tree number of good quality.

We will try to validate or invalidate according to the cases this assertion and we will also try to answer the complementary yet important questions for management : 1) Does optimal silviculture depend on the term T ?

2) Which role the various parameters of the model play : economic parameters p(.), δ and growth parameters g(.), q ?

The optimization problem (P) can be rewritten just by replacing e(t) by -dn(t) dt :

where C is the whole set of curves :

By an integration by part we deduce :

under the same constraints as in the initial problem.

We are here defining the function ξ(.), the relative increase of the basal area s, by ξ(t) = s ′ (t) s(t) .

By applying the results of Lemma 4.1 (iii) we deduce the following Proposition (with proof in Annex B) :

, then the optimal trajectory is E 0 .

From the remark following the Lemma 4.1 we deduce :

7 Proof of Proposition 4.1

(i) Let's consider the auxiliary problem which consists in maximizing, at each time t, the integrand :

dt n(t) =P ′ s (s(t), t)g(r(n(t), s(t)))V (t) + P ′ t (s(t), t)n(t) =kh 0 (t)e -δt (αs(t) α-1 g(r(n(t), s(t)))V (t)δ h (t)n(t)s(t) α )

We denote y = ns b and z(y, s, t) by : z(y, s, t) = αs α-1 g(Ays

. Then z ′ s > 0, from Lemma 4.1 (iib) we deduce the result.

(ii) We then denote y = ns α and z(y, s, t) by : z(y, s, t) = αs α-1 g(Ays q 2 -α )V (t)δ h (t)y.


Double Hashing Thresholds via Local Weak Convergence

M. Leconte Technicolor -INRIA mathieu.leconte@inria.fr Abstract-A lot of interest has recently arisen in the analysis of multiple-choice "cuckoo hashing" schemes. In this context, a main performance criterion is the load threshold under which the hashing scheme is able to build a valid hashtable with high probability in the limit of large systems; various techniques have successfully been used to answer this question (differential equations, combinatorics, cavity method) for increasing levels of generality of the model. However, the hashing scheme analysed so far is quite utopic in that it requires to generate a lot of independent, fully random choices. Schemes with reduced randomness exists, such as "double hashing", which is expected to provide similar asymptotic results as the ideal scheme, yet they have been more resistant to analysis so far. In this paper, we point out that the approach via the cavity method extends quite naturally to the analysis of double hashing and allows to compute the corresponding threshold. The path followed is to show that the graph induced by the double hashing scheme has the same local weak limit as the one obtained with full randomness.

I. INTRODUCTION

The hashing paradigm is as follows: we are given n items and m keys, and we want to assign a key to each item so as to be able to retrieve the items efficiently. In the most basic setting, we want to have exactly one key per item and at most one item per key. Critical performance metrics for such systems are the size of the hashtable (the number of keys) needed to hold a given number of items, and the time it takes to either retrieve or insert items. The multiple-choice hashing strategy is one that guaranties a constant look-up time. It consists in pre-defining a set of d ≤ m keys for each item, which is then only allowed to pick a key within that set. Of course, depending on the choices of the pre-defined sets of keys, it may be impossible to handle simultaneously some sets of items and inserting a new item may not be easy. As for the second issue, cuckoo hashing [START_REF] Pagh | Cuckoo hashing[END_REF] is a simple, randomized way to search for a new valid assignement upon arrival of a new item: one first checks whether one of the keys pre-defined for the new item is available, in which case it suffices to pick one of these; otherwise, one of the predefined keys is chosen at random and re-allocated to the new item. The same procedure is then used for the item that has just been evicted (which is thus treated as a new item). In a fully random context, multiple-choice hashing has been shown to have good performance: if the set of d keys pre-defined for each item are chosen independently and uniformly at random among all sets of d keys, then there exists a threshold τ * such that, if n = τ m with τ < τ * , in the limit of m, n → ∞ cukoo hashing will yield a valid hashtable with probability tending to 1 (which we refer to as w.h.p. for with high probability). Furthermore, the cuckoo hashing update procedure is also quite efficient [START_REF] Frieze | An analysis of randomwalk cuckoo hashing[END_REF], [START_REF] Fountoulakis | On the insertion time of cuckoo hashing[END_REF], in that the insertion time is polylogarithmic in the system size w.h.p. A recent refinement [START_REF] Khosla | Balls into bins made faster[END_REF] of the update algorithm also achieves constant average insertion time w.h.p. when items are only inserted. The focus here being on the space utilization rather than on the insertion time, i.e. on the efficiency of multiplechoice hashing rather than on that of cuckoo hashing, we do not expand on this.

A major drawback of the fully random (multiple-choice) hashing scheme described above is the amount of randomness involved: the standard approach requires n independent, uniform choices of sets of d keys among m. Phrased differently, this means roughly d log m independent random bits per item. Generating unbiased, perfectly independent random bits does not come for free [START_REF] Zuckerman | Simulating bpp using a general weak random source[END_REF], [START_REF] Vembu | Generating random bits from an arbitrary source: Fundamental limits[END_REF], therefore a lot of effort has been put into reducing this amount of randomness needed (see [START_REF] Mitzenmacher | Balanced allocations and double hashing[END_REF] and references therein for an account of some directions investigated so far). A candidate alternative is double hashing [START_REF] Guibas | The analysis of double hashing[END_REF], [START_REF] Lueker | More analysis of double hashing[END_REF], which seems to have similar performances as the fully random one while requiring only roughly 2 log m independent random bits per item: assume m is a prime number and label the m keys with the integers from 1 to m; for each item, independently draw two random numbers f ∈ {1, . . . , m} and g ∈ {1, . . . , m-1 2 }; the predefined set of keys associated with a couple (f, g) are the keys labeled f +ig (mod m), for i ∈ {0, . . . d-1}. Although the reduced randomness of the double hashing scheme is the very reason this method could be prefered over fully random hashing, it also makes the theoretical analyzis of its performance more difficult. In this paper, we show that the load threshold under which double hashing succeeds w.h.p. is essentially the threshold τ * of standard fully random multiple-choice hashing, thereby confirming the hypothesis that double hashing has very similar performance to fully random hashing. More specifically, we show that, if n = τ m with τ < τ * , in the limit of m, n → ∞ double hashing is able to handle all but o(n) items w.h.p. In the case of fully random hashing, one would be able to state a similar result with o(1) instead of o(n). It is very likely that we could obtain the same result using an additional argument (as in [START_REF] Lelarge | A new approach to the orientation of random hypergraphs[END_REF], [START_REF] Leconte | Convergence of multivariate belief propagation, with applications to cuckoo hashing and load balancing[END_REF]), however we believe the current result is sufficient to demonstrate the efficiency of the approach and we do not cross this last gap. It is also of interest to note that there are many extensions of the basic setting considered in this paper (one key per item and at most one item per key). The approach used in this paper is also suited for the analysis of those extensions. For the sake of clarity we do not address them here, although the results go through with only minor modifications.

The remainder of this paper is structured as follows: in the next section, we point out some relevant previous work on the subject. In section III, we define the graphs and random graphs used to model the hashing techniques considered in this paper; in section IV, we recall the concept of local weak limit of random graphs, which is the basis of the cavity method approach for computing the load threshold τ * . In section V, we show that the local weak limit of the random graphs induced by double hashing is identical to that of the random graphs induced by fully random hashing, which implies that the two methods have essentially the same load threshold.

II. PREVIOUS WORK ON PERFORMANCE OF HASHING SCHEMES FOR LOAD BALANCING

The load threshold τ * under which fully random hashing works w.h.p. can be computed by various techniques [START_REF] Fountoulakis | Orientability of random hypergraphs and the power of multiple choices[END_REF], [START_REF] Dietzfelbinger | Tight thresholds for cuckoo hashing via xorsat[END_REF], [START_REF] Frieze | Maximum matchings in random bipartite graphs and the space utilization of cuckoo hash tables[END_REF]. For d ≥ 3, it is given by

τ * = ξ * d(1 -e -ξ * ) d-1 ,
where ξ * is the unique solution of the equation

d = ξ(1 -e -ξ ) 1 -e -ξ -ξe -ξ .
In the basic setting (one key per item and at most one item per key) as well as in extensions (more items allowed per key [START_REF] Dietzfelbinger | Balanced allocation and dictionaries with tightly packed constant size bins[END_REF], [START_REF] Cain | The random graph threshold for k-orientiability and a fast algorithm for optimal multiple-choice allocation[END_REF], [START_REF] Fernholz | The k-orientability thresholds for gn, p[END_REF], more keys required per item [START_REF] Gao | Load balancing and orientability thresholds for random hypergraphs[END_REF], [START_REF] Lelarge | A new approach to the orientation of random hypergraphs[END_REF], restrictions on the number of times a couple (item,key) is used [START_REF] Leconte | Convergence of multivariate belief propagation, with applications to cuckoo hashing and load balancing[END_REF]), the load threshold τ * is related to the emergence of a certain core of the bipartite graph with a large density. It can also be viewed as a problem of orientability of hypergraphs. In all the cases, the expression of the load threshold τ * involves a solution of a fixed-point equation which can be traced back ot the cavity method, as shown in [START_REF] Bordenave | Matchings on infinite graphs[END_REF], [START_REF] Lelarge | A new approach to the orientation of random hypergraphs[END_REF], [START_REF] Leconte | Convergence of multivariate belief propagation, with applications to cuckoo hashing and load balancing[END_REF] as they are related to the fixed-point messages of the belief propagation message passing algorithm.

More generally than hashing systems, this work is related to load balancing. A common formulation for the load balancing problem is that we have m bins, in which we are going to throw n balls. If we throw each balls uniformly at random independently of the previous ones, we end up with the most loaded bin holding log n log log n (1 + o(1)) balls. However, if before throwing each ball we check d bins at random and throw the ball in the least loaded bin, then the most loaded bin will only hold log log n log d + O(1) balls [START_REF] Azar | Balanced allocations[END_REF], which represents an exponential improvement. In the same setup, where bins have no limit capacity, double hashing achieves the same performances [START_REF] Mitzenmacher | Peeling arguments and double hashing[END_REF], [START_REF] Mitzenmacher | Balanced allocations and double hashing[END_REF] as fully random choices. However, in order to guaranty a worst-case constant look-up time, one cannot allow a key to refer to an unbounded number of items, which motivates the hashing setup presented in Section I.

III. HASHING GRAPH AND RANDOM GRAPH MODELS

A multiple-choice hashing system can be represented as a bipartite graph G = (L ∪ R, E) that we call the hashing graph, where the left part L represents the items and the right part R the keys; the edges of the graph indicate which keys can be assigned to which items. In the case of fully random hashing, the d edges adjacent to a left-vertex can be anything, while for double hashing they are restricted to having a certain structure. Given such a hashing graph G, it is possible to build a valid hashtable if we can associate each left-vertex (items) to a right-vertex (keys) without collision, i.e. if there exists a left-perfect matching of G. The items with which the hashing system is presented are random and thus the hashing graph is random, following a different distribution whether we consider fully random or double hashing. Thus, the goal is to determine whether the random graphs obtained have a left-perfect matching w.h.p. Note that for fully random hashing, the graph G is simply a uniform random graph with fixed degree d on the left.

For our analysis of double hashing, it is convenient to see its hashing graph as obtained in the following manner, where we merge all the items which are given the same choices of keys and keep track of the multiplicity of the vertices: we start from the bipartite graph G = ( L ∪ R, E) which contains all the information about the choices of keys offered to every item that can potentially be inserted in the hashtable. We have |R| = m and | L| = m 2 , with a leftvertex in G for each different neighborhood (different choice of keys) allowed by the local structure constraints of the hashing technique. If we label the vertices in R with the integers from 1 to m and the vertices in L with the couples (f, g) with f ∈ {1, . . . , m} and g ∈ {1, . . . , m-1 2 }, there is an edge in G between vertex r ∈ R and l = (f, g) ∈ L if there exists i ∈ {0, . . . , d-1} such that r = f +ig (mod m). To build the hashing graph, f and g are drawn independently at random for each of the n items, which corresponds exactly to a sampling with replacement of n left-vertices in L. We denote by L ′ ⊂ L the set of left-vertices obtained in this way, and by

Z ′ l ∼ Bin(n, 1/ m 2 ) the multiplicity of the left- vertices l. The induced graph G ′ = G[L ′ ∪ R]
is exactly the graph obtained via double-hashing. This main result is stated in Corollary 1.

As we are interested in the regime n = τ m → ∞, it is equivalent for our problem to suppose n ∼ Poi(τ m) instead. Indeed, the Poi(τ m) random variable is concentrated around τ m with only logarithmic fluctuations for m → ∞, therefore the existence and the value of a load threshold τ * are identical in the two models. Thus, we consider instead the set L obtained by taking Z l ∼ Poi( 2τ m-1 ) copies of each vertex in L, independently for each vertex in L, and keeping in L only those vertices with at least one copy, as we did before for L ′ . We will focus mainly on the random graph G = G[L ∪ R]; we denote by G m a sample of this random graph for a particular value of m (τ being fixed and n ∼ Poi(τ m)).

IV. LOCAL WEAK CONVERGENCE AND THE CAVITY METHOD

As we mentioned earlier, we intend to show that the approach via the cavity method for computing the value of the load threshold τ * can be extended to double hashing. We will not explain here how or why the cavity method works; the interested reader can refer to the vast literature on that particular subject. However, we need to explain the notion of local weak convergence of a sequence of graphs, which is used in this paper and is at the basis of the cavity approach. The framework used here is that of [START_REF] Aldous | Processes on unimodular random networks[END_REF].

First, we recall that a sequence of probability measures (ρ m ) m∈N on a certain metric space S converges weakly to a probability measure ρ, which we denote by ρ m ρ, if for every bounded continuous function f : S → R, f dρ m converges to f dρ. We focus on locally finite rooted graphs (or rooted networks, if we want to include marks on the edges). A rooted graph (G, v) is a graph G together with a distinguished vertex v of G which is called the root. A rooted isomorphism of rooted graphs is a graph isomorphism which sends the root of one to the root of the other, i.e. a one-toone mapping of the vertices of a graph to those of another graph which preserves the root, the edges (and their marks). We denote by [G, v] the isomorphism class of (G, v) and by G * the set of rooted isomorphism classes of rooted locally finite graphs. We endow G * with the following metric: let the distance between [G, v] and [G ′ , v ′ ] be 1/(1 + δ), where δ is the supremum of those k ∈ N such that there is some rooted isomorphism of the balls of graph-distance radius k around the roots of G and G ′ . With this metric on G * , the definition of weak convergence applies to probability measures on G * . Given a finite graph G, there is a natural way to form a probability measure on G * : we let U (G) be the distribution over G * obtained by choosing a uniform random vertex of G as a root. Then, we say a sequence (G m ) m∈N converges locally weakly towards a measure ρ on G * when the sequence of distributions (U (G m )) m∈N converges weakly towards ρ.

In the case of fully random hashing and with n = τ m, the hashing graph is a uniform random graph with fixed degree d on the left. As m → ∞, such a random graph, when rooted at a random right-vertex, is known to converge locally weakly to a two-step Galton-Watson tree with distributions Poi(τ d) and Deterministic(d). A sample of such a (potentially infinite) tree is constructed as follows: the root (a key) has degree Poi(τ ); every vertex at odd distance from the root (items) has exactly d -1 children; every vertex at even distance from the root (keys) has a Poi(τ d) number of children.

We should now explain how this notion of local weak convergence is used to obtain the load threshold τ * . We noted previously that the existence of a valid hashtable is equivalent to that of a left-perfect matching of the hashing graph. Clearly, such a matching is of maximum size among all matchings and its size is exactly equal to n. In [START_REF] Bordenave | Matchings on infinite graphs[END_REF], [START_REF] Lelarge | A new approach to the orientation of random hypergraphs[END_REF], [START_REF] Leconte | Convergence of multivariate belief propagation, with applications to cuckoo hashing and load balancing[END_REF], the cavity method was used to compute the asymptotic density of a maximum-size matching for any sequence of random graphs converging in the local weak sense to a Galton-Watson tree. It is then possible to determine the threshold τ * under which a maximum-sized matching has size equivalent to n = τ m as m → ∞. Having τ < τ * thus guaranties the existence of a matching of size no(n) as m → ∞; therefore all but at most o(n) items can be handled. In a second shorter step, it was shown that the two thresholds τ * and τ * are actually equal for fully random hashing.

The advantage of the approach presented in this section is that the threshold τ * is already computed for all graphs converging in the local weak sense towards two-step Galton-Watson trees. Therefore, to compute the load threshold for another hashing scheme, one only needs to perform two smaller steps: [START_REF] Pagh | Cuckoo hashing[END_REF] show the appropriate local weak convergence of the hashing graphs, (2) check the equality of the two threshold τ * and τ * for the particular hashing scheme used. In the next section, we stop before this second step and show that the thresholds τ * are identical for the fully random hashing and double hashing, as both random graph models have the same two-step Galton-Watson tree as a local weak limit.

V. LOCAL WEAK CONVERGENCE OF DOUBLE HASHING GRAPH

A. Main Results and Overview of the Proof

The proof consists in similar steps as in [START_REF] Bordenave | Lecture notes on random graphs and probabilistic combinatorial optimization!! draft in construction!![END_REF], with the difference that in our case the graphs are bipartite and have a more constrained local structure.

Let us call G m = (L m ∪ R m , E m ) the random graph obtained at finite m and let ρ Gm be the distribution of [G m , r 0 ], where r 0 is a random root in R m ; let ρ m be the average of ρ Gm over the random graph G m . Let also ρ be the law of a two-step Galton-Watson tree with laws Poi(τ d) and Deterministic(d), as defined in the previous section. In a first step, we will show that ρ m converges weakly towards ρ as m → ∞. Then, we will show that ρ Gm is concentrated around ρ m so that ρ Gm almost surely converges weakly towards ρ. Explained differently, there are two sources of randomness involved in the sampling of a rooted graph [G m , r 0 ]: the first one is artificial and is due to r 0 being a random right-vertex. It has been introduced to turn any fixed graph into a distribution over rooted graphs, which then allows to consider weak convergence of a sequence of such distributions. The second one is inherent to the hashing scheme used: it comes from the fact the couple (f, g) determining the choice of keys offered to each item is chosen at random upon arrival of the item, and it turns G m itself into a random graph, in a way described in Section III. We first show that the measure ρ m obtained by averaging over both sources of randomness at the same time converges weakly towards the law ρ of a two-step Galton-Watson, before showing that the averaging over the random choices offered to the items is not required for the convergence to hold, which means we do not need to average over different realization of a large hashtable for the results to hold.

We will thus show the following results: the first one is the intermediate result where we show local weak convergence with both types of averaging together.

Proposition 1: The average distribution ρ m of the double hashing random rooted graph [G m , r 0 ] rooted at a random vertex r 0 ∈ R m converges weakly towards the distribution of a two-step Galton-Watson tree. In other words, we have ρ m ρ. The result above, together with some results on concentration of measure, leads to our main theorem, which is almost sure local weak convergence of the hashing graph:

Theorem 1: The distribution ρ Gm of the double hashing random graph G m almost surely converges weakly towards the distribution of a two-step Galton-Watson tree. In other words, ρ Gm ρ a.s. Finally, the announced result on the load threshold of double hashing is contained in the following corollary, which is a direct consequence of Theorem 1 and Theorem 2 from [START_REF] Lelarge | A new approach to the orientation of random hypergraphs[END_REF] or equivalently Theorem 2.1 from [START_REF] Leconte | Convergence of multivariate belief propagation, with applications to cuckoo hashing and load balancing[END_REF].

Corollary 1: The load threshold τ * under which all items can be inserted in the hashtable for fully random hashing w.h.p. is also the load threshold under which all n items but o(n) can be inserted for double hashing w.h.p.

In the next subsection, we introduce a main tool for our proofs: a two-step breadth-first search (BFS) exploration of the hashing graph G m ; and illustrate its analysis through the simple example of upper-bounding the number of short cycles in G m .

B. Breadth-First Search Exploration

When it is clear from the context, we simply write G instead of G m . For any k ∈ N * and any graph G * , we let N k G * (v) denote the k-hop neighborhood of vertex v in the graph G * ; to ease the notation, when G * = G we simply write N k (v), and when k = 1 we write N G * (v).

We consider the breadth-first search (BFS) exploration of the graph G, starting from a random right-vertex r 0 ∈ R. At each step, we select an active right-vertex and explore all its left-neighbors as well as their own right-neighbors. We define the sets C R t and C L t of connected right-and leftvertices at time t, the list of active right-vertices A R t , and the sets of unexplored right-and left-vertices U R t and U L t . We have

C R 0 = C L 0 = ∅, A R 0 = (r 0 ), U R 0 = R \ {r 0 }, U L 0 = L.
At each step t, we proceed as follows:

• we let r t be the first element in the list A R t (so that it is the closest vertex to r 0 in A R t );

• if A R t = ∅, we do the following updates: C R t+1 = C R t ∪ {r t }, C L t+1 = C L t ∪ N (r t ), A R t+1 = (A R t \ {r t }, N 2 (r t ) ∩ U R t ), U R t+1 = U R t \ N 2 (r t ), U L t+1 = U L t \ N G (r t ).
We let X R t be the number of vertices added to A R t \ {r t } at step t, i.e. X R t = |N 2 (r t ) ∩ U R t |, and X L t (resp. X L t ) be the number of left-vertices of G (resp. G) explored at step t, i.e.

X L t = |N (r t ) ∩ U L t | and X L t = |N G (r t ) ∩ U L t |
, where we identify the vertices in G and G with the same label. We also let θ = inf{t ≥ 1 : A t = ∅} be the time at which we complete the exploration of the connected component of r t in G; θ is a stopping time for the filtration

F t = σ C R i , C L i , A R i , U R i , U L i 0≤i≤t
. With these definitions and for all t ≤ θ, it is clear that

|C R t | = t, |C L t | = t-1 i=0 X L i , |A R t | = 1 + t-1 i=0 (X R i -1) ≤ 1 -t + (d -1) t-1 i=0 X L i , |U R t | = m -1 -t-1 i=0 X R i , |U L t | = m 2 - t-1 i=0 X L i ≥ m-1 2 (m -dt).
To get used to reasoning on this model of random bipartite graphs with local structure constraints as well as to understand why double hashing graphs converge to trees, we start off by bounding the number of cycles of a given length in G:

Lemma 1: Let C k be the number of cycles of length k in G. We have

E [C 2k ] ≤ d 2k τ k + O(1/m) Proof: Let r 1 , . . . , r k ∈ R be distinct vertices of R.
For r 1 , . . . , r k to be on a cycle of size 2k of G in this order, there must exist distinct vertices l 1 , . . . , l k ∈ L such that, for all i, l i connects r i-1 and r i , where r 0 = r k . Fixing r and r ′ distinct in R, the number of vertices connecting them in G is exactly d 2 (as setting the set of indices {i, i ′ } such that f + ig = r and f + i ′ g = r ′ or f + i ′ g = r and f + ig = r ′ leaves exactly one possibility for the couple (f, g) ∈ L). Therefore, the probability that r and r ′ are connected by a vertex in L is P r and r ′ are connected by a vertex in

L = 1 -e -2τ m-1 ( d 2 ) = d 2 2τ m + O(1/m 2 ).
Furthermore, it is clear that the probability of existence of distinct vertices l 1 , . . . , l k in L such that, for each i, l i connects r i-1 and r i , is no larger than the product of the probabilities of existence of each l i independently of the others: P ∀i, ∃l i ∈ L connecting r i-1 and r i ; l i = l j , ∀j = i

≤ i P ∃l i ∈ L connecting r i-1 and r i ≤ τ d 2 m k + O(1/m k+1 ).
We conclude by summing over all the m . . . (m-k+1) ≤ m k possible choices of r 1 , . . . , r k .

C. Proofs

We now return to the proof of the main theorem. Lemma 2 shows that the joint distribution of the numbers of children of the left-vertices explored during the first t steps of BFS is asymptotically very close to that we would obtain in a two-step Galton-Watson tree.

Lemma 2: On an enlarged probability space, there exists a sequence (Y L t ) t≥0 of i.i.d. Poi(τ d) variables such that

P (X L 0 , . . . , X L (t-1)∧θ ) = (Y L 0 , . . . , Y L (t-1)∧θ ) ≤ 1 m-1 (τ d 2 t 2 + τ 2 dt + 2τ t).
We denote by d TV (p, q) the distance in total variation between the two distributions p and q. To ease the notation, we will also use d TV for random variables, which will refer to the distance between their distributions. The maximal coupling inequality says that, over an enlarged probability space, there exists a coupling of X and Y such that

P(X = Y ) = d TV (X, Y ).
Proof:

For any t ≤ θ, X L t is a binomial random variable of parameters |U L t ∩ N G (r t )| and 1 -e -2τ m-1 . We have that |N G (r t )| = m-1
2 d and, as any two vertices in R (and in particular r i and r t for any i < t) are connected by exactly d 2 vertices in L, we obtain the bound m-1

2 d -d 2 t ≤ |U L t ∩ N G (r t )| ≤ m-1 2 d. Furthermore, 2τ m-1 -2τ 2 (m-1) 2 ≤ 1 -e -2τ m-1 ≤ 2τ m-1 . It follows P (X L 0 , . . . , X L (t-1)∧θ ) = (Y L 0 , . . . , Y L (t-1)∧θ ) ≤ E (t-1) i=0 1(i ≤ θ)P(X L i = Y L i F i ) and E d TV (X L i , Y L i )|F i ≤ P Bin( d 2 i, 1 -e -2τ m-1 ) = 0 +d TV Bin( m-1 2 d, 1 -e -2τ m-1 ), Bin( m-1 2 d, 2τ m-1 ) +d TV Bin( m-1 2 d, 2τ m-1 ), Poi(τ d) ≤ 2τ i m-1 d 2 + τ 2 m-1 + 2τ m-1 .
We conclude by using the maximal coupling inequality and then summing over i.

Next, Lemma 3 asserts that w.h.p. the graph explored in the first t steps of BFS is a tree.

Lemma 3: For an integer t ≤ θ, the portion of the graph G explored until step t becomes a tree w.h.p. as m, n → ∞. More precisely,

P G[(C R t ∪ A R t ) ∪ C L t ] not a tree ≤ τ 2 d 4 t + τ 2 d 4 t 2 m -1 .
Proof: Loops get formed in the portion of the graph explored at step t when one of the following events occur:

• there exists a vertex r ∈ U R t connected to r t by at least two distinct vertices in L;

• there exists a vertex l ∈ U L t ∩ L which connects r t and some vertex r

∈ A R t . Indeed, assuming G[(C R t ∪ A R t ) ∪ C L t ]
is a tree and none of the two events above occur, it is easy to see that G

[(C R t+1 ∪ A R t+1 ) ∪ C L t+1 ]
P(E t ) ≤ E r∈U R t P(|N (r) ∩ N (r t )| ≥ 2 F t ) ≤ τ 2 d 4 m-1 , P(E ′ t ) ≤ E r∈A R t \{rt} P(N (r) ∩ N (r t ) = ∅ F t ) ≤ τ d 2 m-1 (E[|A R t |] -1). We can bound E[|A R t |] as follows E[|A R t |] -1 ≤ -t + (d -1) t-1 i=0 E[X L i ] ≤ τ d 2 t.
Summing over t yields the desired result.

The hashing graph G obtained may be a tree, but remember that we merged all the items with the same choices of keys into a single left-vertex of G. Lemma 4 shows that each left-vertex explored in the first t steps of BFS actually correspond to a single item, i.e. the items explored were given different choices of keys, w.h.p.

Lemma 4: For an integer t ≤ θ, no left-vertex explored until step t has two copies or more (i.e. Z l ≤ 1 for all l ∈ C L t ) w.h.p. as m, n → ∞. We have

P ∃l ∈ C L t : Z l ≥ 2 ≤ 4τ 3 dt (m -1) 2 .
Proof: The result follows straightforwardly from the union bound:

P ∃l ∈ C L t : Z l ≥ 2 ≤ E[|C L t |](1 -P Poi( 2τ m-1 ) ≤ 1)) = t-1 i=0 E[X L i ] 1 -e -2τ m-1 -2τ m-1 e -2τ m-1 ≤ τ dt 4τ 2 (m-1) 2 .
We are now ready to show that the average distribution ρ m of the k-hop neighborhood of a random root-vertex r 0 in a random graph G m tends to the distribution of a twostep Galton-Watson tree cut at depth k, for any k ∈ N and as m → ∞ (Proposition 1).

Proof: [Proof of Proposition 1] Let k ∈ N. For any finite tree T of depth at most k (or more generally for rooted graphs of radius at most k) and any collection of such trees T , let A T = {[G] ∈ G * , (G) k ≃ T }, where ≃ refers to rooted isomorphism, and A T = ∪ T ∈T A T . For every ǫ > 0 and any finite k, there exists a finite set T of trees of depths at most k in which every odd-depth vertex has degree d and such that ρ(A T ) = T ∈T ρ(A T ) ≥ 1ǫ. Let t be the maximum size of the trees in T ; it means that, with probability at least 1ǫ, a BFS run during t steps on a two-step Galton-Watson tree sampled from ρ will explore at least the k-hop neighborhood of the root.

Let r 0 be a random right-vertex in the random graph G m . According to Lemmas 3 and 4 and with probability at least

1 -τ 2 d 4 t+τ 2 d 4 t 2 m-1 -4τ 3 dt (m-1) 2 → m→∞
1, the graph induced by [G m , r 0 ] on the vertices explored during the first t steps of BFS is a tree, with all left-vertices having exactly one copy. Hence, in particular the offspring of each left-vertex in this tree is of size d-1. Furthermore, according to Lemma 2 and with probability at least 1 -τ d 2 t 2 +τ 2 dt+2τ t m-1 → m→∞ 1, there is a coupling between the numbers of newly explored leftvertices in the graph [G m , r 0 ] during the first t steps of BFS and a sequence of i.i.d. Poi(τ d) random variables. Therefore, on an event of probability tending to 1 as m → ∞, we can couple the exploration during the first t steps of BFS on [G m , r 0 ] and on a two-step Galton-Watson tree sampled from ρ. On the event that these two exploration can be coupled and for any T ∈ T such that the BFS exploration on T will have explored at least the k-hop neighborhood of the root within t steps, it follows the BFS exploration on [G m , r 0 ] will also have explored at least the k-hop neighborhood of r 0 . Then, for any T ∈ T , we have

|P((G m , r 0 ) k ≃ T ) -ρ(A T )| ≤ τ 2 d 4 t+τ 2 d 4 t 2 +τ d 2 t 2 +τ 2 dt+2τ t m-1 + 4τ 3 dt (m-1) 2 . Then, it follows lim m→∞ ρ m (A T ) = ρ(A T ) for any T ∈ T .
For any bounded uniformly continuous function f , there exists k ∈ N such that |f ((G) k )f (G)| ≤ ǫ for all rooted graphs G ∈ G * . For this k we can define a finite collection of trees T as before, such that ρ(T ) ≥ 1ǫ. For m large enough, we have ρ m (T ) ≥ 1 -2ǫ, and

f dρ m -f dρ ≤ ǫ(1 + 3||f || ∞ ) + T ∈T f (T ) ρ m (A T ) -ρ(A T ) .
Letting m → ∞ and then ǫ → 0 completes the proof.

We will now prove almost sure local weak convergence of the random graph G m towards the two-step Galton-Watson tree (Theorem 1). To that end, we use the Azuma-Hoeffding measure-concentration inequality:

Proposition 2: Let M = (M t ) 0≤t≤m be a martingale with respect to a filtration F = (F t ) 0≤t≤m . Suppose there exists constants c 1 , . . . , c m such that, for all 1 ≤ t ≤ m, the following holds:

M t -M t-1 ≤ c t .
Then, for all ǫ > 0, Theorem 1] Let k ∈ N * and H be a rooted graph of radius at most k. We define A H as in the proof of Proposition 1. We let h be the number of rightvertices in H; we focus on m large enough such that ln m ≥ h. Recall that ρ Gm (A H ) = .

P M m -M 0 ≥ ǫ ≤ 2e -2ǫ 2 / 0≤t≤m c 2 m . Proof: [Proof of
The term on the right-hand side is the general term of a convergent series, hence we can conclude using the Borel-Cantelli lemma and then the same argument as for Proposition 1.

VI. CONCLUSION

We have shown that the approach via the cavity method for computing load thresholds of multiple-choice hashing extends naturally from fully random hashing to double hashing. The interest of this approach is that the results on the maximum size of matchings in random graphs can be fully re-used once one proves that the hashing graphs have the same local weak limit for the schemes considered, which leaves only smaller, scheme-dependent steps to check.

The hashing graph of double hashing differs from that of fully random hashing through its local structure, which is constrained in a simple way. However, more complex constraints on the local structure may perhaps be handled via the same approach.

1

 1

m r0∈Rm 1 2 e -τ d m ln ln m τ d → m→∞ 0 . 2 2h 2 d 2 ln 2 m + 2m 2 e -τ d m ln ln m τ d

 1202 ((G m , r 0 ) k ≃ H). For any r ∈ R m = {1, . . . , m}, we define ζ r = (Z l : l ∈ N G (r) and l / ∈ N G (r ′ ), ∀r ′ < r)and we letF t = σ (ζ r ) 1≤r≤t. It is easy to obtain from Chernoff's bound that P(|ζ r | ≥ ln m) ≤ P( l∈N G (r) Z l ≥ ln m) ≤ me -τ d m ln ln m τ d , so that P(∃r ∈ R m such that |ζ r | ≥ ln m) ≤ m We say that G m is valid if |ζ r | < ln m for all r, and that ζ r is valid if |ζ r | < ln m. We will sometimes write G m (ζ) and ρ Gm(ζ) (A H ) to avoid confusion. We define ρ m (A H |G m valid) as the expectation of ρ Gm (A H ) over the random graph G m conditionally on G m valid. Now, we letM t = E[ρ Gm (A H )|F t , G m valid]. We have E[M t+1 |F t ] = E[ρ Gm (A H )|F t , G m valid] = M t , thus M is indeed a martingale with respect to F . Note that M 0 = ρ m (A H |G m valid) and, for G m valid, M m = ρ Gm (A H ). Consider two sequences ζ = (ζ r ) 1≤r≤m and ζ ′ = (ζ ′ r ) 1≤r≤m differing only in one value, say ζ r = ζ ′ r . There are at most (1+|ζ r |(d-1))h right-vertices r 0 in G m (ζ) such that (G m (ζ), r 0 ) k ≃ H and r has a 2-hop neighbor in R ∩ (G m (ζ), r 0 ) k . Assuming ζ r and ζ ′ r are valid, we obtain |ρ Gm(ζ) (A H )ρ Gm(ζ ′ ) (A H )| ≤ 2dh ln m m . Therefore, assuming ζ t is valid for t < r, we have |M r -M r-1 | ≤ E max ζr ,ζ ′ r ρ Gm(ζ) (A H )ρ Gm(ζ ′ ) (A H ) F r-1 , G m valid ≤ 2dh ln m m .Then, the Azuma-Hoeffding inequality yieldsP ρ Gm (A H )ρ m (A H |G m valid) ≥ ǫ G m valid ≤ 2e -mǫ 2 /(2h 2 d 2 ln 2 m) .Furthermore, it is easy to check that ρ m (A H )ρ m (A H |G m valid) ≤ 2P(G m not valid) ≤ 2m 2 e -τ d m ln ln m τ d , and it follows thatP ρ Gm (A H )ρ m (A H ) ≥ ǫ ≤ 2e -mǫ

  is still a tree. At step t, if a loop is created, it must contain a newly explored left-vertex (i.e. a vertex in C L t+1 \C L t ), otherwise it would not involve any new vertex or edge and hence would also be a loop in G[(C R t ∪A R t )∪C L t ]; furthermore, we can always assume that loop contains r t . If the second event does not occur, the newly added vertices are connected to those in (C R t ∪ A R t ) ∪ C L t only through r t , and thus any new loop must be contained in G[{r t }∪(A R t+1 \ A R t )], which is prevented since the first event does not occur either. Let us call E t and E ′ t the two events considered.