Vincent Herbert
email: vincent.herbert@inria.frsarkas@ucalgary.ca

Sumanta Sarkar

On the Triple-Error-Correcting Cyclic Codes with Zero Set {1, 2 i + 1, 2 j + 1}

Keywords: Cyclic codes, Weight enumerator, Zero set, Dual codes, Schaub algorithm, Cyclotomic cosets, Spectral immunity

We consider a class of 3-error-correcting cyclic codes of length 2 m -1 over the two-element field F2. The generator polynomial of a code of this class has zeroes α, α 2 i +1 and α 2 j +1 , where α is a primitive element of the field F2m . In short, {1, 2 i + 1, 2 j + 1} refers to the zero set of these codes. Kasami in 1971 and Bracken and Helleseth in 2009, showed that cyclic codes with zeroes {1, 2 ℓ + 1, 2 3ℓ + 1} and {1, 2 ℓ + 1, 2 2ℓ + 1} respectively are 3-error correcting, where gcd(ℓ, m) = 1. We present a sufficient condition so that the zero set {1, 2 ℓ + 1, 2 pℓ + 1}, gcd(ℓ, m) = 1 gives a 3-error-correcting cyclic code. The question for p > 3 is open. In addition, we determine all the 3-error-correcting cyclic codes in the class {1, 2 i + 1, 2 j + 1} for m < 20. We investigate their weight distribution via their duals and observe that they have the same weight distribution as 3-error-correcting BCH codes for m < 14. Further our experiment shows that these codes are not equivalent to the 3-error-correcting BCH code in general. We also study the Schaub algorithm which determines a lower bound of the minimum distance of a cyclic code. We introduce a pruning strategy to improve the Schaub algorithm. Finally we study the cryptographic property of a Boolean function, called spectral immunity which is directly related to the minimum distance of cyclic codes over F2m . We apply the improved Schaub algorithm in order to find a lower bound of the spectral immunity of a Boolean function related to the zero set {1, 2 i + 1, 2 j + 1}.

Introduction

The family of cyclic codes is a well known class of error-correcting codes. We deal with 3-error-correcting binary cyclic codes and their duals. Let F 2 m be the extension field of degree m of the two-element field F 2 and n be an odd integer. Consider a binary cyclic code C of length n and let α be a primitive n-th root of unity in F 2 m . One can describe C as a principal ideal, in the ring F 2 [X]/(X n -1), with a so-called generator polynomial g over F 2 , where n | (2 m -1). Therefore, the zeroes of g can be used to define C. The zero set Z of C is the set of exponents i of the primitive element α such that α i is a root of g. Note that if z is a root of g, so is z 2 i for all i. In other words, Z is a union of 2-cyclotomic cosets modulo n. Most of the time, Z is shortly described by the list of distinct representatives of these cyclotomic cosets. If the length of the code is 2 m -1, i.e. α is a primitive element of F 2 m , then the code is said to be primitive.

BCH codes form an important class of cyclic codes. A cyclic code generated by g(X) = lcm({M (i) (X)} i∈I ) where M (i) (X) is the minimal polynomial of α i with respect to F 2 and where I is a set of δ -1 consecutive integers is a BCH code with a designed distance δ. In particular, if I = {1, 2, . . . , δ -1}, then the BCH code is said to be narrow sense.

Kasami [START_REF] Kasami | The Weight Enumerators for Several Classes of Subcodes of the 2nd Order Binary Reed-Muller Codes[END_REF] introduced some classes of primitive binary 3-error-correcting cyclic codes (i.e. their minimum distance is 7) similar to 3-error-correcting BCH codes, in the sense that their zero set consists of the union of 3 cyclotomic cosets. These codes have a dimension k ≥ n -3m. This also means that, these codes asymptotically have a high information rate and moreover, their minimum distance is optimal. Indeed, Hamming bound implies that long cyclic codes defined by τ distinct cyclotomic cosets have an error-correcting capacity t ≤ τ . One of those classes is {1, 2 ℓ + 1, 2 3ℓ + 1} with gcd(ℓ, m) = 1 and m odd. Later Bracken and Helleseth in [START_REF] Bracken | Triple-Error-Correcting BCH-like codes[END_REF] discovered the class {1, 2 ℓ +1, 2 2ℓ +1} with gcd(ℓ, m) = 1.

In this article, we investigate the class {1, 2 i + 1, 2 j + 1} for i, j > 1, i = j. We attempt to generalize the subclass {1, 2 ℓ + 1, 2 3ℓ + 1} and {1, 2 ℓ + 1, 2 2ℓ + 1} with gcd(ℓ, m) = 1 to {1, 2 ℓ + 1, 2 pℓ + 1} with gcd(ℓ, m) = 1 and p > 1.

Another class of 3-error-correcting cyclic codes with zero set {1, 2 ℓ-1 + 1, 2 ℓ + 1} for m = 2ℓ + 1 was introduced in [START_REF] Macwilliams | The Theory of Error-Correcting Codes (North-Holland Mathematical Library)[END_REF]. A question was raised whether these codes have the same weight distribution as the BCH code. Later it was proved to be true in [START_REF] Van Der | Non-BCH Triple-Error-Correcting codes[END_REF]. They showed that these dual of these codes have the same weight distribution as the dual of the BCH codes. This motivates us to study the weight distribution of the general class {1, 2 i + 1, 2 j + 1}. By computation up to m < 14, we check that the dual of all the 3-error-correcting cyclic codes that has zero {1, 2 i + 1, 2 j + 1} and which are not BCH, have the same weight distribution as the dual of the 3-error-correcting BCH code.

Further we study the minimum distance of their duals over F 2 and over F 2 m . In the literature, numerous theoretical lower bounds on the minimum distance of cyclic codes are known (e.g. [START_REF] Bose | On a class of error correcting binary group codes[END_REF][START_REF] Hartmann | Generalizations of the BCH bound[END_REF][START_REF] Roos | A new lower bound for the minimum distance of a cyclic code[END_REF][START_REF] Van Lint | On the minimum distance of cyclic codes[END_REF][START_REF] Wolfmann | New bounds on cyclic codes from algebraic curves[END_REF]). Schaub [START_REF] Schaub | A linear complexity approach to cyclic codes[END_REF] has investigated an algorithmic approach to compute a lower bound on the minimum distance of a given cyclic code. This idea is particularly efficient for the codes which have few cyclic subcodes. We improve time-complexity of the Schaub algorithm using a pruning criteria based on BCH bound in order to be able to manage codes with more cyclic subcodes. We compare the Schaub bound with the Hartmann-Tzeng bound and the true minimum distance of duals of codes with the zero set {1, 2 i + 1, 2 j + 1}. Augot and Levy-dit-Vehel [START_REF] Augot | Bounds on the minimum distance of the duals of BCH codes[END_REF] had also applied the Schaub algorithm to find a lower bound of the minimum distance of the dual of BCH codes and found this algorithm gives better results than Ross bound and Weil bound on the dual of BCH codes. Our numerical results show a similar behavior of the minimum distance of the duals of 3-error-correcting cyclic codes in the class {1, 2 i + 1, 2 j + 1} for i, j > 1, i = j.

In the end we study the spectral immunity of a Boolean function which is a cryptographic property. High value of spectral immunity is a necessary condition in order to resist algebraic cryptanalysis of filter generators. In [START_REF] Helleseth | Simplifying algebraic attacks with univariate analysis[END_REF] the connection between spectral immunity and minimum distance of a cyclic code was shown. The spectral immunity of a Boolean function f over F 2 m (in univariate form) is equal to the minimal weight of the 2 m -ary cyclic code of length n = 2 m -1 generated by gcd(f (z), z n + 1) or gcd(f (z) + 1, z n + 1). In this paper, we find a lower bound on the spectral immunity of the Boolean function Tr(g) using the Schaub algorithm, where g is the generator polynomial of the code with the zero set {1, 2 i + 1, 2 j + 1}.

In summary, we have two major contributions in this paper. First, we present a sufficient condition for which the zero set {1, 2 ℓ +1, 2 pℓ +1}, where gcd(ℓ, m) = 1 will give a 3-error-correcting code of length 2 m -1. Secondly, we improve the Schaub algorithm which determines a lower bound of the minimum distance of a cyclic code.

2 Triple-error-correcting cyclic code with the zero set {1, 2 i + 1, 2 j + 1}

Let Z = {a, b, c} be the zero set of a cyclic code C, where a, b and c are the representatives of distinct 2-cyclotomic cosets. Then, the parity check matrix of C is a (3m × n) matrix over F 2 of form:

H =   1 α a α 2a . . . α (n-1)a 1 α b α 2b . . . α (n-1)b 1 α c α 2c . . . α (n-1)c   ,
where each entry in the matrix is represented as an m-bit column vector with respect to a fixed F 2 -basis of F 2 m . Then, the code corresponding to Z is the binary kernel of H.

In Table 1, we present the list of known zero sets that correspond to 3-errorcorrecting cyclic codes of length 2 m -1.

Triple-error-correcting cyclic code with the zero set

{1, 2 ℓ + 1, 2 pℓ + 1}

The zero sets considered in [START_REF] Kasami | The Weight Enumerators for Several Classes of Subcodes of the 2nd Order Binary Reed-Muller Codes[END_REF] and [START_REF] Bracken | Triple-Error-Correcting BCH-like codes[END_REF] are of the form {1, 2 ℓ + 1, 2 pℓ + 1}, where p = 3 and p = 2 respectively. Therefore, it will be interesting to characterize p for which the zero set {1, 2 ℓ + 1, 2 pℓ + 1} always gives a 3-errorcorrecting cyclic code.

In this section, we present a sufficient condition for the zero set Z = {1, 2 ℓ + 1, 2 pℓ + 1} with gcd(ℓ, m) = 1 corresponds to a 3-error-correcting cyclic code.

Zero Set

Conditions References

{1, 2 ℓ + 1, 2 3ℓ + 1} gcd(ℓ, m) = 1 [Kas71] m odd {2 ℓ + 1, 2 3ℓ + 1, 2 5ℓ + 1} gcd(ℓ, m) = 1 [Kas71] m odd {1, 2 ℓ-1 + 1, 2 ℓ + 1} m = 2ℓ + 1 [MS83] m odd {1, 2 ℓ + 1, 2 ℓ+2 + 3} m = 2ℓ + 1 [CGG + 00] m odd {1, 2 ℓ + 1, 2 2ℓ + 1} gcd(ℓ, m) = 1 [BH09] any
, m) = 1. The minimum distance of C is 7 if for all β ∈ F * 2 m , γ ∈ F 2 m , x 2 pℓ +1 p-1 i=0 (βx -(2 ℓ +1) ) 2 iℓ = γ (1)
have less than or equal to 5 solutions for

x in F * 2 m .
The proof is given in Appendix A. We now apply Theorem 1 for p = 2 and p = 3 to show that {1, 2 ℓ + 1, 2 2ℓ + 1} [BH09] and {1, 2 ℓ + 1, 2 3ℓ + 1} [START_REF] Kasami | The Weight Enumerators for Several Classes of Subcodes of the 2nd Order Binary Reed-Muller Codes[END_REF] are two zero sets that give 3-error-correcting cyclic code in Theorem 2 and Theorem 3 respectively.

Theorem 1 opens up the scope of getting new p for which {1, 2 ℓ + 1, 2 pℓ + 1} gives a 3-error-correcting cyclic code by investigating Equation 1 for p > 3 (see Remark 1).

Below we present a consequence of a result given in [START_REF] Bluher | On x q+1 +ax+b[END_REF].

Lemma 2. The equation of the form x 2 ℓ +1 + rx 2 ℓ + tx + s = 0 does not have more than three solutions when gcd(ℓ, m) = 1 for all r, s, t ∈ F 2 m .

Theorem 2. For p = 2, the code {1, 2 ℓ + 1, 2 pℓ + 1} with gcd(ℓ, m) = 1 is a 3-error-correcting cyclic code.

Proof. For p = 2, the equation (1) becomes

γ = x 2 2ℓ +1 (βx -(2 ℓ +1) + β 2 ℓ x -2 ℓ (2 ℓ +1) ), = βx 2 2ℓ -2 ℓ + β 2 ℓ x 1-2 ℓ , = βx 2 ℓ (2 ℓ -1) + β 2 ℓ x -(2 ℓ -1) .
Let us remind β = 0. Since gcd(ℓ, m) = 1 then gcd(2 ℓ -1, 2 m -1) = 1 and thus x → x 2 ℓ -1 is a bijection. Then transforming x = x 2 ℓ -1 we get:

βx 2 ℓ + β 2 ℓ x -1 = γ, x 2 ℓ +1 + γ β x + β 2 ℓ -1 = 0. (2)
Lemma 2 tells that (2) does not have more than three solutions. Therefore, the zero set {1, 2 ℓ + 1, 2 2ℓ + 1} gives a 3-error-correcting cyclic code.

⊓ ⊔ Theorem 3. {1, 2 ℓ + 1, 2 pℓ + 1} is a 3-error-correcting cyclic code if p = 3, m is odd and gcd(ℓ, m) = 1.

Proof. For p = 3, the equation (1) becomes:

γ = x 2 3ℓ +1 (βx -(2 ℓ +1) + β 2 ℓ x -2 ℓ (2 ℓ +1) +β 2 2ℓ x -2 2ℓ (2 ℓ +1) ), = βx 2 ℓ (2 2ℓ -1) + β 2 ℓ x (2 2ℓ -1)(2 ℓ -1) +β 2 2ℓ x 1-2 2ℓ .
Since m is odd and gcd(ℓ, m) = 1, gcd(2 2ℓ -1, 2 m -1) = 1 and so x → x 2 2ℓ -1 is a bijection. Let us recall β = 0. Now replacing x by x 2 2ℓ -1 , we obtain:

βx 2 ℓ + β 2 ℓ x (2 ℓ -1) + β 2 2ℓ x -1 = γ, x 2 ℓ +1 + β 2 ℓ -1 x 2 ℓ + γ β x + β 2 2ℓ -1 = 0. (3)
Then, from Lemma 2, we get that (3) can not have more than three solutions. Therefore, {1, 2 ℓ + 1, 2 3ℓ + 1} is a 3-error-correcting cyclic code.

⊓ ⊔ Remark 1. Consider m odd. In that case, gcd(2 ℓ + 1, 2 m -1) = 1 and so x → x 2 ℓ +1 is a bijection. If we assume p is odd, then applying the transformation x → x -(2 ℓ +1) , (1) becomes :

p-1 i=0 (βx) 2 iℓ = γx 2 (p-1)ℓ +2 (p-2)ℓ +...+1 (4)
Then, it is interesting to find out for which values of p, this equation does not have more than 5 nonzero solutions in F 2 m .

Finding 3-error-correcting cyclic codes with the zero set

{1, 2 i + 1, 2 j + 1} by computation

We use an implementation of Chose-Joux-Mitton algorithm [START_REF] Chose | Fast correlation attacks: An algorithmic point of view[END_REF] to look for words of weight w = 6 in codes with a zero set {1, 2 i + 1, 2 j + 1} for m < 20 and for all i, j. In Table 4, we provide the exhaustive list of triple-error-correcting cyclic codes up to m = 13. This algorithm has time complexity O(n

w 2 ) = O(n 3 ) and space complexity O(n ⌈ w 4 ⌉ ) = O(n 2 ).
From the foundations, this algorithm is employed to find low-weight polynomial multiples in stream cipher cryptanalysis. In our context, it is an efficient algorithm to search codewords of weight smaller than 8. We notice that each of those zero sets {1, 2 i + 1, 2 j + 1} can be written in the form {1, 2 ℓ + 1, 2 pℓ + 1}, where p = 2 or p = 3 up to m < 20. Therefore, the class {1, 2 i + 1, 2 j + 1} of 3-error-correcting cyclic code is completely described by two known classes [START_REF] Kasami | The Weight Enumerators for Several Classes of Subcodes of the 2nd Order Binary Reed-Muller Codes[END_REF][START_REF] Bracken | Triple-Error-Correcting BCH-like codes[END_REF] for m < 20.

The weight distributions of the 3-error-correcting cyclic codes

with the zero set {1,

2 i + 1, 2 j + 1} Weight # Codewords 0 1 N ± √ 8N 2 (N 2 -3N + 2)(N ∓ √ 8N ) 96 N ± √ 2N 2 (5N 2 + 3N -8)(N ∓ √ 2N ) 24 N 2 9N 3 -3N 2 + 10N -16 16 Table 2. Weight distribution of dual of 3-error-correcting BCH code of length 2 m -1, odd m, N = 2 m .
Weight distribution of a linear code C and its dual code C ⊥ are related by the MacWilliams identity [START_REF] Macwilliams | The Theory of Error-Correcting Codes (North-Holland Mathematical Library)[END_REF]. Therefore, knowing the weight distribution of C ⊥ one can obtain the weight distribution of C. The dual code C ⊥ of the code C is the annihilator of C. If C is cyclic, then C ⊥ is also a cyclic code. We denote the zero set of C ⊥ as Z ⊥ . The generator polynomial of C ⊥ is the reciprocal polynomial of h(X) = (X n -1)/g(X). Its roots are the inverses of the roots of h. In other words, it is established that z ∈ Z ⊥ if and only if n -z / ∈ Z. The weight distribution of the dual of 3-error-correcting BCH code for odd m was determined in [START_REF] Kasami | Weight Distributions of BCH Codes[END_REF] and we present it in Table 2. However, for even m, explicit formula for the weight distribution of the dual of the BCH code is not known. In [START_REF] Van Der | Non-BCH Triple-Error-Correcting codes[END_REF], it was shown that the dual of the code {1, 2 ℓ-1 + 1, 2 ℓ + 1} has the same weight distribution as the dual of BCH code. Then from the MacWilliams identity, the code {1, 2 ℓ-1 + 1, 2 ℓ + 1} has the same weight distribution as the BCH code. This motivated us to find if the code {1, 2 i +1, 2 j + 1} has the same weight distribution as the BCH code. As shown in [START_REF] Van Der | Non-BCH Triple-Error-Correcting codes[END_REF] we also study the weight distribution of these codes via their duals. Since these duals have fewer codewords, it is possible to compute the weight distribution for higher extension degrees. For this we implement a concurrent algorithm. We first compute the codewords using Gray coding. After that, we determine their Hamming weights with an hardware-accelerated instruction from the SSE4 instruction set.

For m ≤ 13, we check that duals of the 3-error-correcting cyclic code with zero set {1, 2 i + 1, 2 j + 1} have the same weight distribution as the duals of BCH code as given in Table 2. Therefore, we raise the following question.

Problem 1. Prove or disprove that all the 3-error-correcting cyclic codes with the zero set {1, 2 i + 1, 2 j + 1} of length 2 m -1 have the same weight distribution as the 3-error-correcting BCH code of length 2 m -1.

2.4

Non-equivalence of the 3-error-correcting cyclic codes with the zero set {1, 2 i + 1, 2 j + 1} with the 3-error-correcting BCH code

In [START_REF] Macwilliams | The Theory of Error-Correcting Codes (North-Holland Mathematical Library)[END_REF], it was also asked whether the cyclic code with the zero set {1, 2 ℓ-1 + 1, 2 ℓ + 1} is equivalent to the 3-error-correcting BCH code. They conjectured that they are not. Since our computational result shows that every 3-errorcorrecting cyclic codes having the zero set {1, 2 i +1, 2 j +1} with m ≤ 13 have the same weight distribution as the 3-error-correcting BCH codes, we are interested in the question whether these codes are equivalent to the BCH code. We use the MAGMA implementation of Leon's algorithm [START_REF] Leon | Computing automorphism groups of error-correcting codes[END_REF] to prove the nonequivalence for m = 7 and m = 8. In particular, for m = 7, the 3-error-correcting cyclic code with the zero set {1, 2 ℓ-1 +1, 2 ℓ +1}, for ℓ = 3 is not equivalent to the BCH code {1, 3, 5}. This supports the conjecture proposed in [START_REF] Macwilliams | The Theory of Error-Correcting Codes (North-Holland Mathematical Library)[END_REF]. We employ the support splitting algorithm [START_REF] Sendrier | Finding the permutation between equivalent codes: the support splitting algorithm[END_REF] to prove the non-equivalence for m = 10. The weight enumerator of the hull of a code is an invariant by permutation. The hull of a linear code is the intersection of the code with its dual. We notice that the 3-error-correcting cyclic codes with the zero set {1, 2 i + 1, 2 j + 1} are selforthogonal for m < 20, i.e. the hull of the code is the code itself. If we puncture two equivalent codes in each position, the multiset of weight enumerators of each punctured code is the same for the two codes. This object is the signature of the code that we compute to determine the equivalence of two codes. Cyclic codes have a transitive automorphism group. It implies that if we puncture a cyclic code in any position, we obtain the same weight enumerator for each punctured code. Thus, we puncture the dual codes in one fixed position first. Then we puncture them a second time in each position. We compute the signature of these dual codes. We obtain signatures which are different from the signature of the dual of BCH code {1, 3, 5} for m = 10. So that we can conclude on the non-equivalence of non-BCH 3-error-correcting cyclic codes with the zero set {1, 2 i + 1, 2 j + 1} with BCH. For m = 9, we get the same signature as that of BCH code {1, 3, 5}. This signature is not enough discriminant to collect information on a potential permutation.

We conclude this section by stating that every non-BCH 3-error-correcting code with the zero set {1, 2 i +1, 2 j +1} are not equivalent to the 3-error-correcting BCH code for m = 7, m = 8 and m = 10. The question remains open for m = 9.

An algorithmic approach to compute a lower bound on the minimum distance of cyclic codes

While the weight distribution of the code C and its dual C ⊥ are directly related by the Pless power moment identity, there is as such no theoretical result known which combines the minimum distance of C and C ⊥ .

In this section we discuss on the minimum distance of the dual of tripleerror-correcting cyclic codes with zero set {1, 2 i + 1, 2 j + 1} for i, j > 1 and i = j. Note that these duals are also cyclic. If 0 is in the zero set of C, then C is an even weight code. This implies that dual of the cyclic code with zero set {1, 2 i + 1, 2 j + 1}, has an even minimum distance.

Theoretical lower bounds on the minimum distance of cyclic codes are known (e.g. [START_REF] Bose | On a class of error correcting binary group codes[END_REF][START_REF] Hartmann | Generalizations of the BCH bound[END_REF][START_REF] Roos | A new lower bound for the minimum distance of a cyclic code[END_REF][START_REF] Van Lint | On the minimum distance of cyclic codes[END_REF][START_REF] Wolfmann | New bounds on cyclic codes from algebraic curves[END_REF]). They rely either on properties of regular distribution of certain patterns contained in the zero set, or on the number of rational points of algebraic curves over finite fields. Schaub [START_REF] Massey | Linear complexity in coding theory[END_REF] has investigated an algorithmic approach to compute a lower bound on the minimum distance of cyclic codes. This idea is particularly efficient for the codes which have few cyclic subcodes.

To find a lower bound of the minimum distance of a given dual code, we apply the Schaub [START_REF] Schaub | A linear complexity approach to cyclic codes[END_REF] algorithm. We propose an improvement with a pruning criteria based on BCH bound. Then, we compare the tightness of the Schaub bound and Hartmann-Tzeng bound.

Schaub Algorithm Description

In [START_REF] Schaub | A linear complexity approach to cyclic codes[END_REF], Schaub introduced an algorithm which computes a lower bound of the minimum distance of a cyclic code. This algorithm iteratively applies a method called Rank bounding on symbolic matrices. Basically, this method computes the linear complexity of the infinite periodic sequence derived from the Discrete Fourier Transform of an n-length word c over F q m , where q is a prime power. The Rank bounding method is described in Appendix B. In our instance, q = 2 and n = 2 m -1. Its time complexity is O(n 3 ) = O(2 3m ). Blahut's theorem ensures that this quantity is equal to the Hamming weight of c (e.g. [START_REF] Massey | The discrete Fourier transform in coding and cryptography[END_REF]). In matrix terms, it means that the weight of c is equal to the rank of the circulant matrix B c of order n,

B c =      A 0 A 1 . . . A n-2 A n-1 A 1 A 2 . . . A n-1 A 0 . . . . . . . . . . . . A n-1 A 0 . . . A n-3 A n-2     
, where (A i ) 0≤i≤n-1 is the family of coefficients of Mattson-Solomon polynomial of c.

Consider an n-length cyclic code C over F 2 with zero set Z. On the one hand, the minimum distance d of C is equal to the minimum rank of B c , for all c ∈ C. However, it is impractical to employ Berlekamp-Massey algorithm, whose time complexity is O(n 2 ), to compute the minimum distance of a cyclic code. On the other hand, since C is cyclic, the coefficients of B c satisfy, for all c ∈ C, the property: A z = 0 for all z ∈ Z. In addition, for all c ∈ C, the set of integers i such that A i = 0 forms an union of 2-cyclotomic cosets modulo n. The Schaub algorithm computes a lower bound on the rank of symbolic matrices that we describe below.

Schaub [START_REF] Schaub | A linear complexity approach to cyclic codes[END_REF] defined an arithmetic with three symbols 0, 1 and X. In this notation, 0 stands for null element of F 2 m , 1 stands for any nonzero element of F 2 m and X stands for any element of F 2 m whose nullity or non-nullity is not known.

The commutative semiring ({0, 1, X}, +, * ) is defined with tables:

+ 0 1 X 0 0 1 X 1 1 X X X X X X * 0 1 X 0 0 0 0 1 0 1 X X 0 X X
If κ cyclotomic cosets do not belong to Z, then the Schaub algorithm computes in effect a lower bound on the rank of 2 κ circulant matrices in M ({0, 1}). These matrices have zero coefficient only in the positions determined by the 2 κ corresponding unions of cyclotomic cosets. Thus, the Schaub algorithm has time complexity O(2 3m+κ ).

An improved Schaub algorithm

Each matrix can be identified with a non-linear subcode of C. This code is defined by the codewords of C having zeroes only in the form α i , where i belongs to the set of positions of zeroes coefficients in the first row of the matrix. Each of these codes is of the form D \ E, where D is a cyclic subcode of C and E is the union of all strict cyclic subcodes of D. Let us denominate these non-linear subcodes as constant-zero codes of C, since their codewords (as polynomials) all have the same zeroes. Constant-zero codes form a partition of C. In addition, we can associate to each constant-zero subcode of C, the cyclic subcode of C with corresponding zero set. Rank bounding method consists in constructing a set of necessarily independent rows of the matrix and returns its cardinality. This cardinality is a lower bound on the minimum distance of a constant-zero code of C. Thus, the Schaub bound is the minimum cardinality computed among all the considered subcodes.

Note that each circulant matrix of order n over {0, 1} can be identified by the integer between 0 and 2 n -1 whose binary representation is the first row of the matrix. We assume that the integers are distinct from 2 n -1. If two integers are in the same cyclotomic coset modulo (2 n -1), then the corresponding matrices are row equivalent and thus have the same rank. Thus, it is only necessary to consider one representative of each cyclotomic coset modulo (2 n -1). Further, if the cyclotomic coset modulo (2 n -1) of the integer contains p elements, then we have only to consider the submatrix containing the first p rows of the matrix instead of the whole circulant matrix. Indeed, the circulant matrix is a block matrix and can be split by horizontal lines into n/p blocks where each block is the submatrix described above. Note that p necessarily divides n.

A natural data structure to represent the set of the considered subcodes is the tree. A node A of the tree corresponds to a constant-zero code of C or equivalently to a cyclic subcode of C with zero set denoted Z A . For our purpose, the root node corresponds to the subcode of C which contains the codewords with zeroes exactly in the positions given by Z. A node C is the child of a parent node P if and only if Z C ⊃ Z P and |Z C | = |Z P | + 1. The number of studied subcodes grows exponentially in the number of cyclotomic cosets which are not in Z. In order to reduce time complexity, our strategy is to prune the tree using the BCH bound which is easily computable. One can think to introduce Hartmann-Tzeng bound in the pruning. However, our empirical analysis shows that Hartmann-Tzeng bound slows down the process and does not give better results than BCH. If, in a node, the Schaub bound is found to be smaller than BCH bound of the associated cyclic subcode of C, it becomes pointless to apply Rank bounding method to the subtree whose root is the considered node. Indeed, in each node of the pruned tree, BCH bound is always greater than Schaub bound and thus the Schaub bound is not be updated in this subtree.

Schaub algorithm and algebraic cryptanalysis

Boolean functions are important building blocks in the design of stream ciphers. High algebraic immunity [START_REF] Courtois | Algebraic attacks on stream ciphers with linear feedback[END_REF] is a necessary condition to protect the stream cipher from algebraic attack. Spectral immunity is a related concept to algebraic immunity [START_REF] Helleseth | Simplifying algebraic attacks with univariate analysis[END_REF]. If the spectral immunity of a Boolean function is small, then one can find the initial state of a filter generator in which that Boolean function is used. In [START_REF] Helleseth | Simplifying algebraic attacks with univariate analysis[END_REF] the connection between spectral immunity and minimum distance of a cyclic code was shown which is as follows. The spectral immunity of a Boolean function f over F 2 m (in univariate form) is equal to the minimal weight of the 2 m -ary cyclic codes of length n = 2 m -1 generated by the polynomials gcd(f (z), z n + 1) and gcd(f (z) + 1, z n + 1). Therefore, we need an algorithm to efficiently estimate the minimum distance of a cyclic code over F 2 m .

We apply the improved Schaub algorithm on cyclic codes of length n over F 2 m with a zero set derived from triple-error-correcting cyclic codes with zero set {1, 2 i + 1, 2 j + 1}.

We define the Trace function as Tr(•) :

F 2 m → F 2 , Tr(z) := z + z 2 + z 2 2 + . . . + z 2 m-1 .
In Table 3, we provide a lower bound on the spectral immunity of the boolean function Tr(g(.)) in univariate form over F 2 m , where g is the generator polynomial of a triple-error-correcting cyclic code with zero set {1, 2 i + 1, 2 j + 1} for m < 9. We denote G and H, the generator polynomials of the cyclic codes over F 2 m upon which we apply the Schaub algorithm, G(z) := gcd(Tr(g(z)), z n + 1),

H(z) := z n + 1 G(z) .
Note that the zero set of G contains the zero set of g since g(z) divides both Tr(g(z)) and z n + 1. In addition, H(z) = gcd(Tr(g(z)) + 1, z n + 1), since Tr(g(.)) is a boolean function.

The polynomial g has coefficients over F 2 since it is the product of minimal polynomials with respect to F 2 . Therefore, the boolean function Tr(g(.)) and the generator polynomials G and H have binary coefficients. From Theorem 9 in [START_REF] Van Lint | On the minimum distance of cyclic codes[END_REF], the minimum distance of these codes over F 2 m is the same as that of their subfield subcodes over F 2 . As a consequence, we apply the Schaub algorithm on their binary subfield subcodes, since these one have much less cyclic subcodes.

Computational results

In Table 4, we give the Hartmann-Tzeng bound and the Schaub bound of every dual codes of triple-error-correcting cyclic codes of length 2 m -1 with 5 ≤ m ≤ 13 and Z = {1, 2 i + 1, 2 j + 1}. Some subclasses of codes of this form are well known. We observe the Schaub bound is sharper than Hartmann-Tzeng bound on this class of codes.

Remark 2. As well, we consider the dual of triple-error-correcting cyclic codes defined with the same zero set {1, 2 i +1, 2 j +1} over the alphabet F 2 m for m < 9. It is interesting to note that we obtain the same bound for the codes over F 2 m and their subfield subcodes over F 2 with the Schaub algorithm.

Conclusions

In this work we have discussed on the 3-error-correcting cyclic code that has zero set of the form {1, 2 i + 1, 2 j + 1}. We have presented a sufficient condition so that {1, 2 ℓ + 1, 2 pℓ + 1} corresponds to a 3-error-correcting cyclic code. Although p = 2 and p = 3 are known, our result opens the window to obtain a zero set of the form {1, 2 ℓ + 1, 2 pℓ + 1} for p > 3. For this one needs to find p > 3 such that Equation (1) does not have more than 5 solutions. Remark 1 highlights some reduced form of this equation when p is odd, which may be easier to handle with.

Our experimental result shows that {1, 2 i + 1, 2 j + 1} has the same weight distribution as that of 3-error-correcting BCH code. However, these codes are not equivalent to the BCH code in general, this supports the conjecture proposed by Sloane and MacWilliams [START_REF] Macwilliams | The Theory of Error-Correcting Codes (North-Holland Mathematical Library)[END_REF].

We have improved the Schaub algorithm that finds a lower bound of the minimum distance of a cyclic code. We have used this algorithm to find a lower bound of the spectral immunity of Boolean function Tr(g), where g is the generator polynomial of the code with the zero set {1, 2 i + 1, 2 j + 1}. We would like to see the deployment of this algorithm to find a lower bound of the minimum distance of some other class of cyclic codes in future. A Proof of Theorem 1

The parity check matrix of C is of the form:

H =   1 α α 2 . . . α (n-1) 1 α 2 ℓ +1 α 2(2 ℓ +1) . . . α (n-1)(2 ℓ +1) 1 α 2 pℓ +1 α 2(2 pℓ +1) . . . α (n-1)(2 pℓ +1)   .
Suppose, for a contradiction, that the minimum distance is five. From Lemma 1, we can consider a codeword of weight six. Thus, there exists a set of six dependent columns over F 2 in H. In other words, there exist six distinct elements x, y, z, u, v, w in F * 2 m such that:

   x + y + z + u + v + w = 0, x 2 ℓ +1 + y 2 ℓ +1 + • • • + v 2 ℓ +1 + w 2 ℓ +1
= 0, x 2 pℓ +1 + y 2 pℓ +1 + • • • + v 2 pℓ +1 + w 2 pℓ +1 = 0.

(5) These equations are symmetric in variables x, y, z, u, v, w. Thus, the system (5) can be written as:

   x + y + z = u + v + w = a, x 2 ℓ +1 + y 2 ℓ +1 + z 2 ℓ +1 = u 2 ℓ +1 + v 2 ℓ +1 + w 2 ℓ +1 = b, x 2 pℓ +1 + y 2 pℓ +1 + z 2 pℓ +1 = u 2 pℓ +1 + v 2 pℓ +1 + w 2 pℓ +1 = c,
where a, b, c ∈ F 2 m . Note that b = a 2 ℓ +1 , otherwise we would have:

x + y + z + a = 0, x 2 ℓ +1 + y 2 ℓ +1 + z 2 ℓ +1 + a 2 ℓ +1 = 0, which means there is a codeword of weight 4 in the code with zero set {1, 2 ℓ + 1} where gcd(ℓ, m) = 1. This is impossible since this code has minimum distance 5 from the proof of Lemma 1. Now we consider the equations:

x + y + z = a, x 2 ℓ +1 + y 2 ℓ +1 + z 2 ℓ +1 = b, x 2 pℓ +1 + y 2 pℓ +1 + z 2 pℓ +1 = c. We substitute x with x + a, y with y + a and z with z + a.

x + y + z = 0, (x + a) 2 ℓ +1 + (y + a) 2 ℓ +1 + (z + a) 2 ℓ +1 = b, (x + a) 2 pℓ +1 + (y + a) 2 pℓ +1 + (z + a) 2 pℓ +1 = c.

Then, expanding the second and third equations and using the relation x+y+z = 0 we obtain the following:

x + y + z = 0, x 2 ℓ +1 + y 2 ℓ +1 + z 2 ℓ +1 = b + a 2 ℓ +1 , x 2 pℓ +1 + y 2 pℓ +1 + z 2 pℓ +1 = c + a 2 pℓ +1 .

Next, we replace z = x + y and we get:

x 2 ℓ y + y 2 ℓ x = β, x 2 pℓ y + y 2 pℓ x = γ, (6)
where β = b + a 2 ℓ +1 and γ = c + a 2 pℓ +1 . Note x = 0, since β = 0. We replace y by xy. (6) become:   

x 2 ℓ +1 (y + y 2 ℓ ) = β,

x 2 pℓ +1 (y + y 2 pℓ ) = γ.

(7) (8) From (7), we get: y + y 2 ℓ = βx -(2 ℓ +1) , and raising to the power 2 ℓ repeatedly, we obtain:

y 2 ℓ + y 2 2ℓ = (βx -(2 ℓ +1) ) 2 ℓ , y 2 2ℓ + y 2 3ℓ = (βx -(2 ℓ +1) ) 2 2ℓ , • • • = • • • , y 2 (p-1)ℓ + y 2 pℓ = (βx -(2 ℓ +1) ) 2 (p-1)ℓ .
Adding them all, we get:

y + y 2 pℓ = p-1 i=0 (βx -(2 ℓ +1) ) 2 iℓ .
Then, we get from (8):

x 2 pℓ +1 p-1 i=0 (βx -(2 ℓ +1) ) 2 iℓ = γ.
If this equation does not have more than 5 solutions over F * 2 m , then for all distinct elements x, y, z, u, v, w in F * 2 m (5) is not satisfied. Hence there is no word of weight 6. Then, Lemma 1 implies that the minimum distance is 7.

⊓ ⊔ B Rank Bounding method [START_REF] Schaub | A linear complexity approach to cyclic codes[END_REF] Cases Sum Terms 1 Conclusion 1 0 0 . . .

Fig. 1 .

 1 Fig. 1. Description of fives cases which enable to determine either unknown coefficients, or the independence of considered row, in Rank Bounding method.

Table 1 .

 1 m Known classes of 3-error-correcting cyclic codes of length 2 m -1 Lemma 1. Let d be the minimum distance of the cyclic code C given by the zero set Z = {1, 2 ℓ + 1, 2 pℓ + 1} with gcd(ℓ, m) = 1 of length n = 2 m -1. Then d = 5 or d = 7 and there exists a codeword of weight d + 1.

	Proof. It is known from [Cha98, Theorem 4.2] that cyclic codes with Z = {1, a}
	have minimum distance 5 if and only if the mapping x → x a is almost perfect
	nonlinear (APN). Consequently, if Z = {1, a, b} and if x → x a is APN, then the
	code has minimum distance at least 5 since it is a subcode of code with zero set {1, 2 ℓ + 1}. The mapping x → x 2 ℓ +1 is APN (e.g. [Dob99]) if and only if
	gcd(ℓ, m) = 1. Therefore, the minimum distance of C is at least 5. Further, the
	code C contains the cyclic Reed-Muller code R * (m -3, m) (cf. [Bla02]). Indeed,
	R * (m -3, m) has zero set ∪ i∈N {2 i + 1}. R * (m -3, m) has minimum distance
	2 m-(m-3) = 7. Hence, we are ensured that 5 ≤ d ≤ 7.
	Moreover, by Corollary 17 of [MS83, page 237], the minimum distance d is
	necessarily odd and there is a codeword of weight d + 1.	⊓ ⊔
	Theorem 1. Consider C, the cyclic code of length 2 m -1 with zero set Z =
	{1, 2 ℓ + 1, 2 pℓ + 1} where gcd(ℓ

Table 3 .

 3 Lower bound for spectral immunity of Boolean functions Tr(g(.)) where g is the generator polynomial of binary 3-error-correcting cyclic codes with zero set {1, 2 i + 1, 2 j + 1}

	Code	Zero	Lower Bound for
	Length Set Spectral Immunity
	31 {1, 3, 5} 63 {1, 3, 5} 127 {1, 3, 5} {1, 3, 9} {1, 5, 9} 255 {1, 3, 5} {1, 5, 9}	2 8 11 13 12 14 14

Table 4 .

 4 Bounds on the minimum distance of the dual of the 3-error-correcting cyclic code over F2 of length 2 m -1 with zero set {1, 2 i + 1, 2 j + 1}

	Code	Dual	Hartmann-Tzeng Schaub Minimum
	Length Zero Set	Bound	Bound Distance
	31 63 127 255 511 1023 {1, 3, 5} {1, 3, 5} {1, 3, 5} {1, 3, 5} {1, 3, 9} {1, 5, 9} {1, 3, 5} {1, 5, 9} {1, 3, 5} {1, 3, 9} {1, 3, 17} {1, 5, 9} {1, 5, 17} {1, 9, 17} {1, 9, 17} 2047 {1, 3, 5} {1, 3, 9} {1, 3, 17} {1, 3, 33} {1, 5, 9} {1, 5, 17} {1, 5, 33} {1, 9, 17} {1, 9, 33} {1, 17, 33} 4095 {1, 3, 5} {1, 5, 33} 8191 {1, 3, 5} {1, 3, 9} {1, 3, 17} {1, 3, 65} {1, 5, 17} {1, 5, 33} {1, 5, 65} {1, 9, 17} {1, 9, 33} {1, 9, 65} {1, 17, 33} {1, 33, 65}	8 16 32 32 48 64 96 128 128 128 192 192 224 256 448 512 512 512 512 768 768 768 896 896 960 1024 1536 2048 2048 2048 2048 3072 3072 3072 3584 3584 3584 3840 3968	8 16 48 48 48 96 96 216 212 210 218 212 224 446 448 930 906 906 876 936 872 916 902 902 960 1886 1814 3110 3588 3643 3668 3594 3678 3802 3912 3718 3722 3844 3968	8 16 48 48 48 96 96 224 224 224 224 224 224 448 448 960 960 960 960 960 960 960 960 960 960 1920 1920 3968 3968 3968 3968 3968 3968 3968 3968 3968 3968 3968 3968

Acknowledgments: The authors are highly grateful to Daniel Augot and Pascale Charpin for their valuable suggestions and comments on this work. The authors also thank the anonymous reviewers for their comments which has improved the editorial quality of this paper. The second author also would like to acknowledge the financial support by the French Agence National de la Recherche under contract ANR-06-SERI-013-RAPIDE that he received during this work.

Algorithm 1 Rank Bound

Input: a nonzero matrix M of size n × n, with coefficients over {0, 1}. Output: a lower bound on the rank of M .


Frédéric Morlot
email: frederic.morlot@orange-ftgroup.com

Benoît Fourestié
email: benoit.fourestie@orange-ftgroup.com

Salah- Eddine Elayoubi
email: salaheddine.elayoubi@orange-ftgroup.com

Optimizing the date of an upgrading investment in a data network

Keywords: data flows, brownian motion, quality of service, customer satisfaction

Due to the introduction of new services, the volume of data transferred in mobile networks is rapidly growing and operators periodically face the necessity to upgrade their network. Such upgrades allow them to increase the capacity, and provide adequate Quality of Service (QoS). In this paper we propose a general framework for deriving the optimal date for a network upgrade. We show that this date is the result of a compromise betweenthe decrease of upgrade investment cost and the loss of profit generated by insufficient capacity. The upgrade should hence be performed when the loss of profit, derived using analytical capacity expressions, exceeds the expected discount. The model presented herein accounts for the randomness of the demand and upgrading cost functions, and results are given for a HSDPA network.

I. INTRODUCTION

As new high throughput services such as multimedia applications are introduced in the market, the traffic is significantly growing in mobile telecommunications networks. In particular, data flows are growing exponentially, as they have been doing since the middle 90' in the Internet. To face these soaring volumes of data to be transferred, mobile operators must periodically upgrade their equipments to offer higher throughputs lest they should undergo blocking problems. However as the demand does not increase steadily and must be considered as partly random, the upgrade dates are difficult to be forecast.

In this article we consider a HSDPA cellular network and assume for the sake of simplicity that there exists no priority between data flows. As the demand increases, the operator can choose between the following two strategies:

• with admission control: it sets a lower throughput threshold and thus has to prevent some users from entering the network to guarantee this minimum throughput. As the demand rises, more and more users are rejected: the blocking rate increases. • without admission control: there is no throughput lower threshold. All users are admitted in the network, but the communication duration becomes longer and longer for each user until the network is saturated. The global throughput experienced in the network may result very small. In the remainder of the article, we will assume that the operator has chosen the second option. Note that if it had chosen the admission control strategy, main results presented herein would still hold and generalization elements are given in section V. We further assume that the "fair power" sharing algorithm is used [START_REF] Bonald | France Telecom R&D), Wireless Downlink Data Channels: User Performance and Cell Dimensioning[END_REF] since it can be shown that it has good stationary properties. As the demand rises, the operator increases its profit, since each user transfers more data and thus spends more money. When the network starts experiencing saturation problems, the global troughput and profit may fall. The operator must then upgrade its network and faces the following trade-off:

• The later the investment, the lower individual throughputs and customer satisfaction. Permanent insatisfaction will result into churn and additional loss of profit. • The sooner the investment, the more expansive as upgraded element costs decrease with time. This article aims at proposing a global framework for evaluating the optimal date for network upgrading, and its content is the following:

• First, operator's profit is calculated. Analytical values for capacity, individual throughput and satisfaction will be derived as a function of the demand. The randomness of the rising demand will be accounted for; • Second, upgrading costs will be estimated assuming exponential decrease according to time; • Third, elements are given in case that an admission control strategy is used; • Fourth, results are given for a HSDPA network and show that the more random (volatile) the demand, the earlier the investment should be made. To the best of our knowledge, this paper presents the first endeavour to analytically derive optimal dates of investments in mobile radio networks taking into account randomly growing data transfers, capacity and throughput evaluations and exponentially decreasing investment costs. Such an analysis will be of interest for future operator upgrading decisions.

II. OPERATOR'S PROFIT

The operator profit depends on the amount of data flowed by the network, i.e. on the sum of individual throughputs. However individual throughputs depend on a large number of parameters and are difficult to calculate. The profit thus depends on the following parameters:

• the mean daily traffic demand per cell (in Mbits/sec/cell).

We denote it by X t , where the time t = 0, 1...T varies discretely, for example day by day. [0, T ] stands for the considered time period. • the customer satisfaction -before the upgrading: S 1 t = f 1 (X t ) -after the upgrading: S 2 t = f 2 (X t ) • the daily profit per surface unit generated by X t , possibly balanced by S i t

π i t = g(X t , S i t ) = g(X t , f i (X t ))
If we denote by t 0 the investment date (0 < t 0 < T ), the total profit Π T (t 0 ) actualized at the date t = T is calculated as follows:

Π T (t 0 ) = t0-1 t=0 e ζ(T -t) π 1 t + T t=t0 e ζ(T -t) π 2 t , (1)
where ζ is the actualization rate. For simplicity, we assume that ζ is constant during the period [0, T ].

A. Traffic Demand

We assume that the network is formed by circular cells of radius R and that the demand is uniformly distributed in the cell. We have:

X t = λ t × E[ξ],
where λ t is the mean arrival rate per cell at the date t and E[ξ] is the mean size of a typical data flow. In the following we assume that E[ξ] remains stable during [0, T ], so that X t is proportional to λ t , whereas λ t is growing. This model is strictly equivalent to the case where the number of active users is constant, but they initiate connections more often. Furthermore a model where the rate of calls λ t would be constant with E[ξ] increasing would lead to similar results.

To model the evolution of (X t ) t∈N , let us consider it as the daily sampling of a continuous stochastic process ( Xt ) t∈R + . Usually, one monitors the 24 mean demands over one hour each, keeps the second or third highest, and multiply it by a given factor. As many random phenomena related to a social behaviour (e.g. [START_REF] Halluin | Managing Capacity for Telecommunications Networks Under Uncertainty[END_REF][5]), we assume that Xt is growing according to the following stochastic differential equation:

     d Xt Xt = αdt + σdW t X0 = x 0 almost surely, (2)
where W t is a standard brownian motion. In other words, Xt is a geometric brownian motion, whose solution is given by (see [START_REF] Karoui | Introduction au Calcul stochastique, première partie[END_REF]):

Xt = x 0 e (α-σ 2 /2)t+σWt , t ∈ R +
Note that we have here E [X t ] = x 0 e αt . Hence, α is the trend of the demand, whereas σ is its volatility. Figure 1 shows the daily evolution of the traffic according to time obtained using the described model.

B. Customer Satisfaction

To evaluate the customer satisfaction, let us first compute the flow throughput γ t (r) a user can expect at distance r from the center of the cell as carried out in [START_REF] Bonald | France Telecom R&D), Wireless Downlink Data Channels: User Performance and Cell Dimensioning[END_REF]. The resource of a single downlink channel is time-shared between active users. Denote by φ u the fraction of time the BS transmits to user u, with u φ u = 1. The data rate of user u is then

C u (r) = C(r) × φ u ,
where C(r) is the peak data rate, obtained in the absence of any other user in the cell, i.e., for φ u = 1. When there are x users, the "fair power" sharing is defined by φ u = 1

x . We have:

C(r) = min C 0 , Z b × Γ(r) η + I(r)
,

where:

• C 0 is the maximum peak rate (which depends on channel bandwidth and coding efficiency); • Z is the cell chip rate; • b is a lower bound for energy-per-bit to noise density ratio (E b /N 0 ), such that the error probability will not be higher than a target threshold; • Γ(r) is the path loss between the BS and user u;

• η is the thermal noise to received power ratio; • I(r) is the interference to received power ratio . Figure 2 shows the peak data rate in a hexagonal cell. In the interference term, we only take immediate adjacent cells in account. Furthermore, we assume that the path loss Γ(r) is decreasing according to 1/r 4 .

Note that we considered the ideal case where a continuous set of peak rates is achievable. But in practice, because of coding constraints, the set of achievable peak rates is discrete, and locally constant.

We define the cell load as (see [START_REF] Bonald | France Telecom R&D), Wireless Downlink Data Channels: User Performance and Cell Dimensioning[END_REF]):

ρ t = X t πR 2 R 0 2πrdr C(r) ,
where R is the radius of the cell.

• if ρ t > 1, the cell is overloaded and one can show ( [START_REF] Bonald | France Telecom R&D), Wireless Downlink Data Channels: User Performance and Cell Dimensioning[END_REF]) that the number of active users grows indefinitely; any individual data rate tends to zero, and the cell is saturated. • if ρ t < 1, the cell is underloaded and the number of active users tends to a finite stationary regime.

We are naturally led to introduce:

X max = πR 2 R 0 2πrdr C(r) -1 , (3)
such that ρ t = X t /X max . Hence ρ t < 1 ⇐⇒ X t < X max .

Let us denote the flow throughput of users at distance r by γ t (r) (it is the ratio of the mean flow size to the mean flow duration). Then, Bonald and Proutière [START_REF] Bonald | France Telecom R&D), Wireless Downlink Data Channels: User Performance and Cell Dimensioning[END_REF] show that, if ρ t < 1:

γ t (r) = C(r)(1 -ρ t ). (4)
By assumption, the mean flow size (the numerator in γ t (r)) does not vary significantly, whereas the mean flow duration (the denominator) does as the load increases. Thus, if we want to compute the mean flow throughput, the significant number to calculate is the harmonic mean of γ t (r) over the cell (balanced by the proportion of active users between r and r + dr; remember we suppose users are uniformly distributed in the cell. In other words, we have to calculate the arithmetic mean of 1/γ t (r) over the cell. From (4) we deduce:

γ t = R 2 R 0 2rdr γ t (r) -1 = X max -X t ,
and of course, in overload, γ t = 0.

Finally, we can summarize the whole calculation by:

γ t = (X max -X t ) + (5)
where, for a given real x, x + = max(x, 0). Note that γ t = 0 if and only if X t ≥ X max , i.e., if and only if the cell is saturated. Now we can compute the customer satisfaction, which can reasonably be supposed to depend on γ t . Since subjective satisfactions have been shown to be more sensitive to small variations at low throughputs than at high throughputs, Enderlé and Lagrange propose in [START_REF] Enderlé | User Satisfaction Models and Scheduling Algorithms for Packet-Switched Services in UMTS[END_REF] to model the customer satisfaction as a negative exponential function of the throughput:

S t = e -β/(Xmax-Xt) + .
For example, β can be chosen as: β = log(2) • γ 1/2 , where γ 1/2 is the throughput value ensuring a satisfaction of 50%. Once again, note that S t = 0 if and only if X t ≥ X max , i.e., if and only if the cell is saturated.

C. Daily Profit

Let us recall that X t is sampled day by day, for example during the second or the third highest peak hour. Within 24 hours, the operator transmits µmin(X t , X max ) to active users, where µ is a multiplicative factor between the peak hour and the whole day. Typically, µ ≈ 4 • 3600 ≈ 14000 sec. Since taxation is applied to the volume of transfers and not to the duration, the gross daily profit per cell is given by:

π gross = δµmin(X t , X max ) ( 6
)
where δ is the transfer price (say in $/Mbit). However the gross profit should be weighed by the customer satisfaction to account for the quality of the communications. The net profit is thus calculated as the product of π gross by S t π net = δµmin(X t , X max )e -β/(Xmax-Xt) + .

If S t = 0, i.e. if the cell is saturated, the net profit is null. If the satisfaction is maximal, i.e., S t = 1, the net profit is equal to the gross profit [START_REF] Enderlé | User Satisfaction Models and Scheduling Algorithms for Packet-Switched Services in UMTS[END_REF]. To sum up we have:

π t = δµX t e -β/(Xmax-Xt) if X t < X max π t = 0 otherwise
Intuitively, as the demand rises, X t will increase as will the profit (Fig. 3). Then, the profit will decrease because the unsatisfaction effect becomes dominant. In case no upgrading action is taken, the profit will progressively tend to zero. Once the operator decides to upgrade the network, two solutions are possible:

• make the network denser, i.e. install new base stations in the network • install additional transmitters operating on different frequency bands with the view to increase the throughputs.

In such case we obtain a higher value of X max , so that:

t = δµmin(X t , X ′ max )e -β/(X ′ max -Xt) +
In this article we restrict our analysis to the latter case, i.e., X max changes but other variables do not.

D. Total Profit

As we have now derived π t analytically, we can compute numerically the expectation of Π T (t 0 ) given t 0 using the Monte Carlo method. Results for the calculation of the daily profit calculated using Scilab are shown in Fig. 4. We observe that the higher sigmais , the more important the risk is and the sooner the maximal value of the expected profit is reached. Additionally, the higher sigma is, the slower the curve tends to zero. When σ becomes paramount, the curve flattens out.

From the previous simulation, we can deduce the total profit shown in Fig. 5. We observe that the expected profit diminishes as the volatility increases, because of the risk. To summarize the results:

• for a given value of σ, the earlier the investment, the happier the consumers, and the higher the profit; • at a given date, the higher σ, the lower the profit, since the system can be saturated earlier. However, we must now take the investment cost into account. This is the subject of next section.

III. UPGRADING INVESTMENT

The upgrading cost is a decreasing function of time. The decrease of the cost is due to many factors, for example the R&D progress, and also the serialization in the manufacturing chain. Moore's law states that electronic devices' capacity doubles every 18 months. So in this paper we assume it decreases exponentially:

K(t) = K 0 e -ǫt ,
where ǫ is the depreciation rate. The number to maximize is what we call the total Profitminus-Cost (PMC), actualized at date T :

P M C(t 0 ) = Π T (t 0 ) -K(t 0 )e ζ(T -t0) = Π T (t 0 ) -K 0 e ζT -(ǫ+ζ)t0 .
Note that generally ǫ ≫ ζ, so that ǫ + ζ ≈ ǫ. Results when taking the K(t) component in our simulations and results are shown in Fig. 6. Again we show that the higher the volatility is, the lower the expected profit is. This result was expectable because of the increased risk. Furthermore, we see that the optimal investment date decreases according to the volatility, as shown in Fig. 6. However, works on investment under uncertainty obtain the opposite result [START_REF] Salahaldin | Investing in Sustainable Transport to Relieve Air Pollution under Population-Growth Uncertainty[END_REF]. In fact, as will be shown in next section, the significant value to observe is the optimal actualized investment date, and we will observe the opposite trend.

IV. STRATEGY

As time goes by, we have realizations of the daily demand. These values can thus be used to derive a more precise optimal date than the one we found at t = 0 assuming demand randomness. Instead of calculating expectations, we can now take in account what really happens during the period. We can implement the following algorithm: each day t ∈ [0, T ], we perform deterministic calculations over [0, t] and expectations over [t, T ], with an initial value x t of traffic demand which is now a deterministic parameter. Each computation thus provides us with a new investment date t 0 ∈ [t, T ]. In other words, t 0 is daily actualized. When t 0 becomes close to t, it is time to invest. Note that such a situation (t 0 → t) will inevitably happen by continuity, since at the end of the process we have t = t 0 = T . Interestingly we can take the upgrading delay into account: if for example we need 1 month of delay for the deployment, we will plan it as soon as t 0 -t ≤ 30 days. In Fig. 7, we show the results in the deterministic case, i.e., σ = 0. For t = 0, we find t 0 = 41 days as expected by Fig. ??. Since there is no randomness, t 0 is fixed until t = t 0 . Then, t 0 remains "stuck" to t. In the random case, we can observe some fluctuations before t 0 becomes stuck to t (see Fig. 8). In fact, the optimal date we find (60 days) seems to be higher than the one found at the beginning (30 days). Is it a hazard, or is the mean final best result always greater ? To examine that assumption, we simulate several possible scenarios, in other words, several possible trajectories of the brownian demand. For each trajectory, we apply the actualization algorithm and we compute a final result, averaged over all the simulations (see Fig. 9): now, we are in accordance with classic investment studies.

Fig. 8. Best investment date during actualization algorithm in the random case (σ = 0.05), for one given trajectory. Here, investment will take place at around t = 60 days. When Xt does not grow as fast as expected, the actualized investment date is postponed. But as soon as Xt has sufficiently grown, the actualized date is stabilized.

V. GENERALIZATION So far, the actualization algorithm has been illustrated in a simple situation. However it could be applied to a wide range of different situations. For example, we could introduce:

• admission control: the network then guarantees to each user a minimal data throughput, but some users are blocked. In that case, the satisfaction S t depends on the blocking rate, not on the throughput. [START_REF] Bonald | France Telecom R&D), Wireless Downlink Data Channels: User Performance and Cell Dimensioning[END_REF] gives an analytic expression of the blocking rate:

Fig. 1 .

 1 Fig. 1. Daily traffic evolution. α = -log(0.85)/30 day -1 (which corresponds to an average increasing of 15% per month), σ = 0.04, x 0 = 1.2 Mbit/sec/cell.

Fig. 2 .

 2 Fig. 2. Peak data rate against distance to the center of the cell. The distance and the rate are normalized with respect to the cell's radius and to C 0 . Z/b = 3 Mbits/sec, η = 10%.

Fig. 3 .

 3 Fig. 3. Daily Profit generated by demand. Xmax = 3 Mbits/sec/cell, δ = 0.1 $/Mbit, µ = 14000 sec, β = 0.7 Mbit/sec (which corresponds to γ 1/2 ≈ 1 Mbit/sec).

  π

Fig. 4 .

 4 Fig. 4. Expected Daily Profit during the period. Parameters are identical to those used in former figures.

Fig. 5 .

 5 Fig. 5. Total Expected Profit. Xmax = 8 Mbits/sec/cell after the upgrade. ζ = -log(0.95)/365 day -1 , which corresponds to a 5% / year actualization.

Fig. 6 .

 6 Fig.6. Expected PMC during the period. K 0 = 300 k$ and ǫ = -log(0.5)/365, which corresponds to a 50% / year cost decrease.

Fig. 7 .

 7 Fig. 7. Best investment date during actualization algorithm in the deterministic case. Investment will take place at t = 41 days.

Acknowledgment. We are especially grateful to Thomas Bonald for the useful pieces of advice he gave to us.

, where ρ t is the cell load (refer to section II-B) and m is the maximum allowed number of users in one cell. Then we assume that the customer satisfaction decreases exponentially according to blocking rate, and the actualization algorithm still holds. • different power allocations. In that case, we have to compute a mean satisfaction over the cell, balanced by the probability density function of the user position. • coding constraints, which entail discontinuities in the peak data rate curve [START_REF] Bonald | France Telecom R&D), Wireless Downlink Data Channels: User Performance and Cell Dimensioning[END_REF].

• discontinuities and randomness in the investment cost curve to modelize R&D's breakthroughs. This random component is added to the expected PMC's computation.

VI. CONCLUSION

In this paper, we presented an analytical approach for determining the optimal date of an upgrade investment in mobile radio networks. The strategy we adopted is:

• first, analytically model the operator profit and the customer satisfaction; • second, compute the optimal upgrading date with an actualization algorithm. The first step is of utmost importance because it has to reflect as accurately as possible the real situation of the network: Radio Ressource Management (RRM) algorithms (such as admission control, power allocation) must be taken into account to correctly model the capacity of the network. To illustrate this modelling, we consider a HSDPA data network and show how to calculate the cell capacity, the user's throughput and the satisfaction, and the operator's profit. Note that the demand parameters, such as the volatility, the trend and the mean size of data flows, have to be evaluated carefully because results are rather sensitive to input variations, especially to volatility (Fig. ??). Models should therefore be calibrated using measurements taken in the past to avoid too big discrepancies between forecasts and reality. In the second step, we introduced the cost of the upgrading investment and developed an actualization algorithm that calculates the optimal investment time using the above described analytical modeling. The method presented herein should be of interest for operators, allowing them to make optimal decisions as to when upgrade their radio networks.