F Orieux

O Féron

J.-F Giovannelli

Sampling high-dimensional Gaussian distributions for general linear inverse problems

Keywords: Stochastic sampling, high-dimensional sampling, inverse problem, Bayesian strategy, unsupervised, myopic

This paper is devoted to the problem of sampling Gaussian distributions in high dimension. Solutions exist for two specific structures of inverse covariance: sparse and circulant. The proposed algorithm is valid in a more general case especially as it emerges in linear inverse problems as well as in some hierarchical or latent Gaussian models. It relies on a perturbationoptimization principle: adequate stochastic perturbation of a criterion and optimization of the perturbed criterion. It is proved that the criterion optimizer is a sample of the target distribution. The main motivation is in inverse problems related to general (non-convolutive) linear observation models and their solution in a Bayesian framework implemented through sampling algorithms when existing samplers are infeasible. It finds a direct application in myopic / unsupervised inversion methods as well as in some non-Gaussian inversion methods. An illustration focused on hyperparameter estimation for super-resolution method shows the interest and the feasibility of the proposed algorithm.

I. INTRODUCTION

This work deals with simulation of high-dimensional Gaussian and conditional Gaussian distributions. The difficulty of the problem is directly related to handling high-dimensional covariances R and precision matrices Q = R -1 . The problem has already been investigated and solutions exist in two cases.

• When Q is sparse, two strategies are available. The first one [1, chap. 8], relies on a parallel Gibbs sampler based on a chessboard-like decomposition. It takes advantage of the sparsity of Q to update simultaneously large blocks of variables. The second strategy [START_REF] Rue | Fast sampling of Gaussian Markov random fields[END_REF][START_REF] Lalanne | Stochastic artificial retinas: algorithm, optoelectronic circuits, and implementation[END_REF] relies on a Cholesky decomposition Q = L t L: a sample x is obtained by solving the linear system Lx = ε, where ε is a zeromean white Gaussian vector. The sparsity of Q ensures feasible numerical factorization and the sparsity of L ensures feasible numerical resolution of the linear system.

• [START_REF] Chellappa | Classification of textures using Gaussian Markov random fields[END_REF][START_REF] Chellappa | Markov Random Fields: Theory and Application[END_REF] propose a solution for circulant matrix Q, even non-sparse. In this case, the covariance is diagonal in the Fourier domain: the sampling is based on independent sampling of the Fourier coefficients. Finally, the sample is computed by FFT and it has been used in [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF][START_REF] Charbonnier | Deterministic edge-preserving regularization in computed imaging[END_REF][START_REF] Giovannelli | Unsupervised Bayesian convex deconvolution based on a field with an explicit partition function[END_REF][START_REF] Orieux | Bayesian estimation of regularization and point spread function parameters for Wiener-Hunt deconvolution[END_REF][START_REF] Helgason | Fast and exact synthesis of stationary multivariate Gaussian time series using circulant embedding[END_REF]. To our knowledge there is no solution for more general structure in high dimension because factorization (Cholesky, QR, square root,. . . ), diagonalization and inversion of Q and R are numerically infeasible. The obstacle is due to both computational cost and memory footprint. The proposed F. Orieux is with Pasteur Institute, 75015 Paris, France, orieux@pasteur.fr. O. Féron is with EDF Research & Developments, OSIRIS, 92140 Clamart, France, olivier-2.feron@edf.fr. J.-F. Giovannelli is with Univ. Bordeaux, IMS, 33400 Talence, France, Giova@IMS-Bordeaux.fr.

algorithm overcomes this obstacle when Q is of the form

Q = K k=1 M t k R -1 k M k (1)
as it appears in inverse problems [START_REF] Idier | Bayesian Approach to Inverse Problems[END_REF]. Indeed, let us consider the general linear forward model y = Ax+n, where y, n and x are the observation, the noise and the unknown image and A is a linear operator. Consider, again, two prior distributions for n and x that are Gaussian conditionally on a parameter θ. This framework is very general: it includes linear inverse problems [START_REF] Idier | Bayesian Approach to Inverse Problems[END_REF] as well as some hierarchical or latent Gaussian models [START_REF] Rue | Approximate bayesian inference for latent gaussian models by using integrated nested laplace approximations[END_REF] and it can be used in many applications. In image reconstruction, it covers a majority of current problems, e.g. unsupervised [START_REF] Giovannelli | Unsupervised Bayesian convex deconvolution based on a field with an explicit partition function[END_REF] or myopic (semi-blind) [START_REF] Orieux | Bayesian estimation of regularization and point spread function parameters for Wiener-Hunt deconvolution[END_REF] inverse problems, by including acquisition parameters and hyperparameters in θ.

Moreover, the framework also includes non-linear models, based on conditional linear models such as bilinear or multilinear ones (see Section III-B). The framework also covers some non-stationary or inhomogeneous Gaussian priors and non-Gaussian priors involving auxiliary / latent variables [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF][START_REF] Giovannelli | Unsupervised Bayesian convex deconvolution based on a field with an explicit partition function[END_REF][START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF][START_REF] Féron | Microwave imaging of piecewise constant objects in a 2D-TE configuration[END_REF][START_REF] Ayasso | Joint NDT image restoration and segmentation using Gauss-Markov-Potts prior models and variational Bayesian computation[END_REF]] (e.g., location or scale mixtures of Gaussian), by including these variables in θ.

Let us focus on the joint estimation of x and θ from the posterior p(x, θ|y). It commonly requires the handling of the conditional posterior p(x|θ, y) that is Gaussian with precision matrix Q of the form (1), as will be shown in section II-B. In the general case, Q is neither sparse nor circulant so existing sampling algorithms fail when the dimension of x is very large while the proposed one handles this case. It relies on a perturbation-optimization principle: adequate stochastic perturbation of a quadratic criterion and optimization of the perturbed criterion. A recent paper [START_REF] Tan | Efficient sparse Bayesian learning via Gibbs sampling[END_REF] briefly describes a similar algorithm for compressed sensing in signal processing. Our paper deepens and generalizes this contribution.

Subsequently, Section II presents the proposed algorithm and its direct application to linear inverse problems. Section III gives an illustration through an academic problem in superresolution. Section IV presents conclusions and perspectives.

II. PERTURBATION-OPTIMIZATION ALGORITHM

A. Description

We focus on the problem of sampling from a target Gaussian distribution whose precision matrix Q is in the form [START_REF] Winkler | Image Analysis, Random Fields and Markov Chain Monte Carlo Methods[END_REF]. When Q is neither sparse nor circulant, existing algorithms fail in high dimension because of an excessive memory footprint as illustrated in section III. We propose a solution based on the Perturbation-Optimization (PO) algorithm described hereafter, whose memory footprint is far smaller.

Proposition 1: The optimizer x of criterion (5) resulting from Algorithm 1 is Gaussian

x ∼ N Q -1 K k=1 M t k R -1 k m k , Q -1 . (2)
Proof: The optimizer x of criterion ( 5) is explicit:

x = K k=1 M t k R -1 k M k -1 K k=1 M t k R -1 k η k = Q -1 K k=1 M t k R -1 k η k . (3)
It is clearly a Gaussian vector as a linear combination of K Gaussian vectors. Its expectation and covariance are calculated below using elementary algebra: from ( 4) and ( 3), we have

E [ x] = Q -1 k M t k R -1 k E [η k ] = Q -1 k M t k R -1 k m k E x x t = Q -1   k,k ′ M t k R -1 k E η k η t k ′ R -1 k ′ M k ′   Q -1 = Q -1 k M t k R -1 k R k + m k m t k R -1 k M k Q -1 = Q -1 + E [ x] E [ x] t
that completes the proof.

The feasibility of Step P clearly depends on the capability to sample from Gaussian distributions N (m k , R k ). It is usually the case in inverse problems and it will be actually the case in super-resolution applications shown in section III-A and in other contributions shortly described in section III-B.

Regarding

Step O, J being quadratic, a large literature [17] is available about its numerical optimization, e.g. gradient procedure (standard, corrected, conjugate, optimal step size. . . ). Such algorithms require the computation of criterion ( 5) and its gradient. The feasibility of Step O clearly depends on the capability to compute that without the storage of large matrices. It is usually the case in inverse problems and it will be actually the case in applications shown in section III-A and described in section III-B.

However, the desired sample is the exact optimizer, so, Step O could require N iterations of a conjugate gradient algorithm for a problem of dimension N . Therefore the complexity could be O(N 3 ) that is equivalent to the one of a Cholesky decomposition. However, the optimization procedure can be stopped earlier without practical loss of precision and the complexity falls down to O(P N 2 ) for P iterations. In addition, for a band matrix, the complexity of the proposed algorithm becomes O(M P N ) and the one of the Cholesky decomposition becomes O(M N 2 ). Anyway, the main advantage of the proposed algorithm is its reduced memory footprint: it avoids the storage of neither Q nor its (Cholesky, QR, square root,. . . ) factors.

Algorithm 1 : Perturbation-Optimization algorithm.

1:

Step P (Perturbation): Generate independent vectors

η k ∼ N (m k , R k ), for k = 1, . . . K (4) 2:
Step O (Optimization): Compute x as the minimizer of

J(x) = K k=1 (η k -M k x) t R -1 k (η k -M k x) (5)
Remark 1: Still regarding Step O, it would be awkward if Q was badly scaled, but it is not the case here for the following reason. In usual ill-conditioned inverse problems, A is badlyscaled but the aim of regularization is precisely to overcome this difficulty and to produce a well-scaled matrix Q.

B. Application to inverse problems

The purpose is to solve an inverse problem, stated by the forward model y = Ax + n, in a Bayesian framework based on the following models:

• A describes any observation system that can depend on unknown acquisition parameters,

• priors for the noise n and the object x are Gaussian N (m n , R n ) and N (m x , R x ), conditionally on a set of hyperparameters and auxiliary variables.

In a general statement, acquisition parameters, hyperparameters and auxiliary variables are collected in θ. The general inverse problem then consists in estimating x and θ through the posterior p(x, θ|y). Its exploration can be achieved by means of a Gibbs sampler which iteratively samples from p(θ|x, y) and p(x|θ, y). The conditional posterior p(x|y, θ) is a correlated Gaussian distribution:

N (m post x , R post x ) with R post x = A t R -1 n A + R -1 x -1 m post x = R post x A t R -1 n [y -m n ] + R -1 x m x
where θ is embedded in A, R n and R x for simpler notations.

If A has no particular properties, Q = (R post x ) -1 is neither sparse nor circulant, and existing sampling algorithms are not applicable. The PO algorithm makes it possible to sample from

N (m post x , R post x ) by applying Algorithm 1 with K = 2, M 1 = A, M 2 = I, R 1 = R n , R 2 = R x , m 1 = m n and m 2 = m x .
In this context, it can be said that the optimization procedure converts prior samples into a posterior one.

III. ILLUSTRATION

The proposed PO algorithm makes it possible to resort to stochastic sampling algorithms in inverse problems providing two main advances:

• capability to jointly estimate extra unknowns included in θ (acquisition parameters, hyperparameters, . . . ),

• access to the entire unknown distribution providing uncertainties (standard deviation, credibility interval,. . . ).

These advances are illustrated in the present section.

A. Unsupervised super-resolution

We detail an application of the proposed PO algorithm to the super-resolution (SR) academic problem: several blurred and down-sampled (low resolution) images of a scene are available in order to retrieve the original (high resolution) scene [START_REF] Park | Super-resolution image reconstruction: a technical overview[END_REF][START_REF] Rochefort | An improved observation model for SR under affine motion[END_REF]. It is shown that the crucial novelty, enabled by the proposed PO algorithm, is to allow the use of sampling algorithms in SR methods and thus to provide hyperparameter estimation.

We resort to the standard forward model in SR: y = Ax + n = P Hx + n. In this equation, y ∈ R M collects the pixels of the low resolution images, here 5 images of 128×128 pixels (M = 81920) and x ∈ R N collects the pixels of the original image, here 256 × 256 pixels (N = 65536). The noise n ∈ R M accounts for measurement and modeling errors. H is a N × N circulant convolution matrix that accounts for the convolution part of the observation system. Practically, the impulse response is a Laplace shape with FWHM of 4 pixels. P is a M × N decimation matrix: it is a binary matrix indicating which pixel is observed. Finally, A is a M × N matrix (that is to say 81920 × 65536). The prior distribution for n is N (0, γ -1 n I) and the one for x is N (0, γ -1 x D t D) where D is the N × N circulant convolution matrix of the Laplacian filter. The hyperparameters γ n and γ x are unknown and their prior law are Jeffreys'. The posterior [START_REF] Orieux | Bayesian estimation of regularization and point spread function parameters for Wiener-Hunt deconvolution[END_REF] is

p(x, γ n , γ x |y) ∝ γ M/2-1 n γ (N -1)/2-1 x exp -γ n y -P Hx 2 /2 -γ x Dx 2 /2 . (6)
It is explored by a Gibbs sampler: iteratively sampling γ n , γ x and x under their respective posterior conditional distribution

p(γ (k) n |x, γ x , y) = G 1 + M/2, 2/ y -P Hx (k-1) 2 p(γ (k) x |x, γ n , y) = G 1 + (N -1)/2, 2/ Dx (k-1) 2 p(x (k) |γ x , γ n , y) = N (m post x , R post x ) with R post x = γ (k) n H t P t P H + γ (k) x D t D -1 and m post x = γ (k) n R post x P t H t y.
The conditional posteriors for the hyperparameters are Gamma distributions so they are easy to sample. The conditional posterior for x is Gaussian, but the use of existing algorithms is impossible due to the structure and the size of R post x . Regarding the structure, according to Section II-B, with A = P H: A is non-circulant due to the decimation and A is non-sparse due to large support of the impulse response. Regarding the size, R post x (and its Cholesky factor) is a huge N × N matrix, that is to say 65536×65536 and its footprint in memory would be 32 GB. As a consequence, neither the precision matrix nor its Cholesky factor can be stored on standard computers.

On the contrary, the proposed PO algorithm only requires the storage of four 256 × 256 matrices and its footprint in memory is only 2 MB that is easy to manage on standard computers. Regarding the computational cost:

• Step P requires a sample under each prior distribution: x is computed by FFT (see item 2 of Section I) and n is trivially computed since it is a white noise.

•

Step O is achived by a conjugate gradient procedure with optimal step size. It only requires computations of convolutions (by FFT), decimation and zero-padding. So, the proposed PO algorithm is feasible and it easily provides a desired sample. Practically, it takes 1 about one second (i.e. around P = 50 gradient iterations) to obtain one sample. Fig. 1 shows the iterates and illustrates the operation and convergence. After a burn-in period of about 25 iterations, the algorithm is in its converged state and the total number of iterations is 59 to ensure a good exploration of the distribution. Histograms approximating marginal posteriors are also given and the posterior means are γ n ≈ 7.7 and γ x ≈ 2.2 × 10 -3 .

Concerning the images themselves, results are shown in Fig. 2: the estimated image in 2(c) clearly shows a better resolution than the data in 2(b) and it is visually close to the original image in 2(a). Nevertheless, it is important to keep in mind that, w.r.t. other SR methods, the proposed PO algorithm does not improve image quality itself but the crucial novelty is to allow for hyperparameter estimation. In this sense, it is clear that the approach produces correct hyperparameters i.e. correct balance between data and prior. Moreover, uncertainties are derived from the samples through the posterior standard deviation. It is illustrated in Fig. 2(d): the true image is inside the 99% credibility interval around the estimate. As a conclusion, the proposed PO algorithm makes it possible to resort to sampling algorithms in SR method whereas it was not possible before. It then enables hyperparameter estimation while other SR methods require hand-made hyperparameter tuning. In addition, it enables to compute uncertainties based on posterior standard deviation.

B. Three other examples

The PO algorithm has been used in three other contexts: electromagnetic inverse scattering [START_REF] Féron | Microwave imaging of piecewise constant objects in a 2D-TE configuration[END_REF], fluorescent microscopy trough structured illumination [START_REF] Orieux | Bayesian estimation for optimized structured illumination microscopy[END_REF] and superresolution from data provided the Herschel observatory in astronomy [?]. The problems are tackled in a Bayesian framework and implemented by means of stochastic sampling. In these contexts, the distribution for the object given the other variables is Gaussian with large size precision matrix. Its structure is neither sparse nor circulant making the use of existing algorithms impossible. This is due to non-linearity and label variables in [START_REF] Féron | Microwave imaging of piecewise constant objects in a 2D-TE configuration[END_REF] and non-invariance of the observation model in [20, ?]. Nevertheless, the precision matrix is in the form (1), so, the proposed PO is applicable.

IV. CONCLUSION

The paper presents an algorithm for sampling high-dimensional Gaussian distributions when existing algorithms are infeasible. It relies on a perturbation-optimization principle: adequate stochastic perturbation of a criterion and optimization of the perturbed criterion. It is shown that the criterion optimizer is a sample of the target distribution. The algorithm is applicable for a particular decomposition of the precision matrix that emerges in general linear inverse problems.

There is a wide class of applications, in particular any processing problem based on a conditional linear forward model and conditional Gaussian priors for noise and object. The interest and the feasibility of the proposed algorithm have been illustrated in [14, 20, ?] and in this paper on a more academic super-resolution problem allowing automatic tuning of hyperparameters.

An interesting perspective deals with the case of stopped optimization procedure. It is a question under consderation to prove that, embedded in a Gibbs loop, a finite number (maybe one) of iteration of the optimization step is enough to guarantee convergence towards the target law.

Fig. 1 .

 1 Fig. 1. Chains and histograms of hyperparameters γn and γx.

Fig. 2 .

 2 Fig. 2. Image reconstruction: true image 2(a), one of the low resolution images 2(b) and the proposed estimate 2(c). The plot 2(d) is a true image slice inside the 99% credibility interval around the estimate.

The algorithm is implemented within the computing environment Matlab on a PC with a 3 GHz CPU and 3 GB of RAM.

V. ACKNOWLEDGMENT

The authors would like to thank J. IDIER (IRRCyN) for inspiration of this work [21], T. RODET and A. DJAFARI (L2S), for fruitful discussions, C. VACAR (IMS) and P. SZACHERSKI (IMS and CEA-LETI) for carefully reading the paper.


Houcein El Abdalaoui
email: elhoucein.elabdalaoui@univ-rouen.fr

Mohamed Dahmoune
email: mohamed.dahmoune@u-pec.fr

Djelloul Ziadi
email: djelloul.ziadi@univ-rouen.fr

On the transition reduction problem for finite automata

We are interested in the problem of transition reduction of nondeterministic automata. We present some results on the reduction of the automata recognizing the language L(En) denoted by the regular expres-

These results can be used in the general case of the transition reduction problem.

Introduction

Minimizing the number of states of an automaton is a subject that has been studied extensively since the 1950s, both in the deterministic case and the nondeterministic case [START_REF] Moore | Gedanken-experiments on sequentiel machines[END_REF][START_REF] Hopcroft | An n log n algorithm for minimizing states in a finite automaton[END_REF]. However, works on the minimization of the number of transitions have appeared recently.

In 1997, J. Hromkovic et al. [START_REF] Hromkovic | Translating regular expressions into small epsilon-free nondeterministic finite automata[END_REF] have proposed an algorithm based on the concept of Common Follow Set of a regular expression, that converts a regular expression of size n into a finite state automaton with O(n) states, O(n log n) transitions as lower bound and O(n log 2 n) transitions as upper bound. Muscholl et al. [START_REF] Hagenah | Computing epsilon-free nfa from regular expressions in O(n log 2 (n)) time[END_REF], showed that this algorithm can be implemented in time O(n log 2 n). In [START_REF] Khorsi | Fast equation automaton computation[END_REF] Ouardi and Ziadi, based on the ZPC structure [START_REF] Champarnaud | An efficient computation of the equation k-automaton of a regular k-expression[END_REF], gave an O(n log 2 n) algorithm to convert a weighted regular expression of size n into a weighted automaton having O(n) states and O(n log 2 n) transitions. In [START_REF] Geffert | Translation of binary regular expressions into nondeterministic ǫfree automata with transitions[END_REF], Viliam Geffert showed that every regular expression of size n over a fixed alphabet of s symbols can be converted into a nondeterministic ε-free finite state automaton with O(sn log n) transitions.

Lower bound was improved by Yuri Lifshits [START_REF] Yuri | A lower bound on the size of e-free nfa corresponding to a regular expression[END_REF] to Ω( n log 2 n log log n ), after, Schnitger [START_REF] Schnitger | Regular expressions and nfas without ǫ-transitions[END_REF] improved it to Ω(n log 2 n) transitions.

In [START_REF] Cox | Minimal number of edges in e-free non-deterministic finite automata (nfa) for regular expression (1 + ε) • (2 + ε) • (3 + ε) • • • (n + ε). The On-Line Encyclopedia of Integer Sequences[END_REF], R. Cox has done an exhaustive search to find the transition minimal automata of L(E n ) for n = 1 to 7. He has also used an heuristic approach that construct transition reduced automata for n = 8 to 10.

Here, we are able to produce an algorithm for which the number of transitions is minimal for L(E n ) languages class, in the sense that, asymptotically, this number of transitions is equivalent to n log 2 n (see Section 6).

We mention that most of complexity results mentioned above are obtained from the study of L(E n ) languages class. This class of languages corresponds to a simple class of automata, in which, the minimization of the number of transitions is difficult and not obvious. The study of this class of languages, can also find its application in bioinformatics, since that L(E n ) is exactly the set of all sub-sequences of the word 1.2.3 . . . n.

Our approach to reduce the number of transitions of a nondeterministic homogeneous finite state automaton is based on the decomposition of the transition table of the automaton into blocks. This decomposition is based on the concept of Common Follow Sets. From a block decomposition we construct an automaton with less transitions than the initial automaton. See Figure 1.

The main problem in our approach is to find a good block decomposition of the transition table (even the best one). In the case where this matrix is lower triangular or upper triangular, finding a minimal decomposition block leading to a minimal transition automaton, is not evident. Our study is focused on the upper triangular matrix, which corresponds to the transition table of the deterministic minimal automaton recognizing the L(E n ) language. The case of lower triangular matrix can be obtained in a similar manner. In this paper we present the following results: At first, in Section 3, we extend the concept of Common Follow Sets to homogeneous automata. Then, in Section 4 we introduce particular decompositions called Z-partitions associated with expressions E n . Then, in Section 5 we introduce the notion of Z-tree to represent any Z-partition by a binary tree. Then, we propose an algorithm of O(n log n) time complexity to generate the Z-minimal trees. We finish our study by experimental results and a last section in which we show that our algorithms construct automata with a number of transitions equivalents to n log 2 (n) which is the minimal lower bound according to Schnitger.

Notation and terminology

We recall the basics of regular expressions, languages and finite state machines and introduce the notation that we use. Let Σ be a non-empty finite set of symbols, called alphabet. The set of all the words over Σ is denoted by Σ * . The empty word is denoted by ε. A language over Σ is a subset of Σ * . A finite automaton over Σ is a 5-tuple A = (Q, Σ, I, δ, F ) where Q is a set of states, I is a subset of Q whose elements are the initial states, F is a subset of Q whose elements are the final states, δ is a subset of the cartesian product Q × Σ × Q whose elements are the transitions. A transition (q, a, p) ∈ δ goes from the head q to the tail p. A path in A is a sequence of transitions (q i , a i , q i+1 ), i = 1 to n, of consecutive transitions. Its label is the word w = a 1 a 2 • • • a n . A word w ∈ Σ * is recognized by the automaton A if there is a path with label w such that q 1 ∈ I and q n+1 ∈ F . The language recognized by the automaton A is the set of words that are recognized by A. The automaton A is homogeneous if for all (q, a, p), (q ′ , a ′ , p ′ ) ∈ δ, p = p ′ implies that a = a ′ , in this case we write h(p) = a. The function h assigns to each non-initial state q of an homogeneous automaton the symbol that is the unique label of all the transitions having q as tail.

In Appendix A we recall the basics of asymptotic notations.

CFS for homogeneous automata

J. Hromkovic et al. [START_REF] Hromkovic | Translating regular expressions into small epsilon-free nondeterministic finite automata[END_REF] have given an elegant algorithm based on the notion of Common Follow Sets, to convert a regular expression of size n into a nondeterministic finite automaton having O(n) states and O(n log 2 n) transitions. This notion can be easily extended to homogeneous automata. Let A = (Q, Σ, {q 0 }, δ, F ) be an homogeneous automaton. In order to capture the final states in the A, we introduce a dummy state denoted by # which is not in Q. We define over Q the function f ollow as follows:

f ollow(q) = {p | (q, a, p) ∈ δ} ∪ {#} if q ∈ F, {p | (q, a, p) ∈ δ} otherwise.
Let q ∈ Q be a state in A, we denote by dec(q

) = {Q 1 , Q 2 , • • • , Q k } (where Q i ⊆ f ollow(q)
) any decomposition of the set f ollow(q), i.e. f ollow(q) = Qi∈dec(q) Q i .

In the case where dec(q) is a partition of the set f ollow(q), the decomposition dec(q) will be called a partition decomposition. Figure 2 provides examples of decompositions.

Definition 1 (Common Follow Sets System).

Let A be a homogeneous automaton. A CFS system for A is given as S(A) = (dec(q)) q∈Q , where each dec(q) ⊆ 2 Q is a decomposition of f ollow(q). Definition 2 (CFS automaton). Let A = (Q, Σ, {q 0 }, δ, F ) be a homogeneous automaton and S(A) an associated Common Follow Sets system. The Common Follow Sets automaton associated with S(A) is defined by

C S(A) = (Q ′ , Σ, I ′ , δ ′ , F ′ ) where -Q ′ = q∈Q dec(q) -I ′ = dec(q 0 ) -For Q 1 ∈ Q ′ , Q 1 ∈ F ′ if and only if # ∈ Q 1 -δ ′ = {(Q 1 , a, Q 2 ) | ∃q ∈ Q 1 s.t. h(q) = a and Q 2 ∈ dec(q)}.
Theorem 1. Let A be a homogeneous automaton, S(A) be a Common Follow Sets System associated with A and C S(A) its Common Follow Sets automaton. Then C S(A) and A recognize the same language. This theorem can be proved in the same way as Theorem 5 of the paper of J. Hromkovic et al. [START_REF] Hromkovic | Translating regular expressions into small epsilon-free nondeterministic finite automata[END_REF].

To evaluate the number of transitions in the automaton C S(A) we define over the states of A two functions, -a(q) = |dec(q)| the size of the decomposition of the set f ollow(q) -b(q

) = |{Q 1 ∈ Q ′ | q ∈ Q 1 }| the number of states in Q ′ that contain the state q. Lemma 1. The number of transitions T C in C S(A) is such that T C ≤ q∈Q a(q)b(q).
It is easy to see that if for all p, q ∈ Q\{q 0 } such that p = q, we have h(p) = h(q) then the equality holds. From Lemma 1, we can see that the number of transitions in a CFS automaton depends on the decomposition system. A decomposition which is not a partition will induce more transitions than a partition decomposition. Therefore in the following we are interested only in partition decompositions. As it was mentioned in the introduction our study will focus on the CFS automata associated with the family of automata (A n ) n≥1 . The automaton A n = (Q, Σ, I, δ, F ) is defined by:

-Σ = {1, 2, . . . , n} -Q = Σ ∪ {0} -F = Q -I = {0}, δ = {(p, q, q) ∈ Q × Σ × Q | q > p}.
Figure 3 shows two CFS automata associated with the automaton A 3 .

In the next sections we present two algorithms that construct particular CFS systems which correspond to CFS automata with a reduced number of transitions. In the last section we give comparative and experimental results.

{1, 2, 3, #} {2, 3, #} {3, #} {#} 2 1 2 3 2 3 3 q a(q) b(q) a(q)b(q) dec(q) 0 1 0 0 {{1, 2, 3, #}} 1 1 1 1 {{2, 3, #}} 2 1 2 2 {{3, #}} 3 1 3 3 {{#}} TC = 6 {3, #} {1, 2} {2} {3} {#} 1 1 2 2 2 2 3 3 q a(q) b(q) a(q)b(q) dec(q) 0 2 0 0 {{1, 2}, {3, #}} 1 2 1 2 {{2}, {3, #}} 2 2 2 4 {{3}, {#}} 3 1 2 2 {{#}} TC = 8

A n Reduction Algorithm

The notion of Z-partition is nowhere introduced formally. In the following we are interested in reducing the number of transitions in the automaton A n . The following algorithm computes particular CFS systems S(A n ) that provide C S(An) automata with small number of transitions and having n + 1 states. Let

E n = (1 + ε) • (2 + ε) • (3 + ε) • • • (n + ε
) be a regular expression, it is easy to see that the language denoted by the expression E n is exactly the language recognized by the automaton A n . We have:

Proposition 1. Each transition minimal automaton that recognizes L(E n ) has ex- actly n + 1 states.
This proposition can be proved using properties of the universal automaton [START_REF] Lombardy | The universal automaton[END_REF] of L(E n ).

For fixed n, a CFS system produced by the following algorithm will be denoted by Z(A n ). The set of all Z(A n ) will be denoted CF SZ(n). Our aim in this section is to compute all minimal decompositions in CF SZ(n).

Algorithm 1 CFSPartitions(n)

Require: n ∈ N Ensure: Z(An) 1: Q ← {0, 1, 2, 3, 4, . . . , n} 2: for i = 0 to n do 3: f ollow(i) ← {j ∈ Q|j > i} ∪ {#} 4: dec(i) ← φ 5:
Qi ← φ 6: end for 7: for i = 0 to n do 8:

Choose

j in Q 9: Q ← Q\{j} 10: Qj ← f ollow(j) 11: for all k ∈ Q do 12: if (Qj ⊆ f ollow(k)) then 13: dec(k) ← dec(k) ∪ {Qj} 14: f ollow(k) ← f ollow(k)\Qj 15:
end if 16:

end for 17: end for Proposition 2. The number of all Z(A n ) CFS partition systems is the n th Catalan number:

|CF SZ(n)| = 1 n+1 2n n .
The successive choice of values of j (line 8) leads to a permutation of size n. So, each CFS partition system Z(A n ) can be associated with at least one permutation of size n.

The following algorithm is a recursive version of Algorithm 1. Its first call is done by RecursiveDecomposition(0, n). Without loss of generality we associate in this last algorithm the dummy state # to the number n + 1. At each call, Algorithm 2 constructs one block from the transition matrix M , for the call RecursiveDecomposition(n 1 , n 2 ) and the choice of j (line 2), it produces the block B j which is the submatrix M

[j..n 1 ; j..n 2 ]. Algorithm 2 RecursiveDecomposition (n 1 , n 2 ) Require: n1, n2 ∈ N Ensure: Z(An) when n1 = 0 and n2 = n 1: if n1 ≤ n2 then 2:
Choose an integer j between n1 and n2, j ∈ {n1, . . . , n2} 3:

Qj ← {j + 1, . . . , n2 + 1} 4:
for k = n1 to j do 5:

dec(k) ← dec(k) ∪ {Qj } 6:
end for 7:

Bj = M [j..n1; j..n2] 8:
RecursiveDecomposition(n1, j -1) 9:

RecursiveDecomposition(j + 1, n2) 10: end if Example 1. In this example we shows the CFS partition systems associated with permutations (0, 2, 1, 3) and permutation (1, 0, 2, 3).

dec(0) = { { 1 , 2 , 3 , # } } dec(1) = { { 2 } , { 3 , # } } dec(2) = { { 3 , # } } dec(3) = { { # } } 0 : 1 2 3 # 1 : 2 3 # 2 : 3 # 3 : # dec(0) = { { 1 } , { 2 , 3 , # } } dec(1) = { { 2 , 3 , # } } dec(2) = { { 3 , # } } dec(3) = { { # } } 0 : 1 2 3 # 1 : 2 3 # 2 : 3 # 3 : # For permutation (0, 2, 1, 3) the first block B 0 is 1 2 3 # , the second block B 2 is 3 # 3 #
, the third block B 1 is 2 and B 3 = # is the last one.

Proposition 3. The computation of all minimal partition system

Z(A n ) in CF SZ(n) can be done in time O(n!).
Proof. This can be done by calling the nondeterministic Algorithm 1 or 2 for each possible execution.

Remark 1. By the use of the dynamic programming, we can improve the exponential brute force method to a polynomial algorithm as shown in Algorithm ??.

In the following sections we will introduce our second algorithm which is based on trees. It computes efficiently the reduced Z(A n ) systems.

Tree based reduction

A binary tree is a structure defined on a finite set of nodes that either contains no nodes, or is made of three disjoint sets of nodes:

a root node a binary tree called its left subtree a binary tree called its right subtree.

The binary tree that contains no nodes is called the empty tree. If the left subtree is non-empty, its root is called the left child of the root of the entire tree. Likewise, the root of a non-empty right subtree is the right child of the root of the entire tree. Therefore, in a full binary tree each node is either a leaf or has degree exactly 2, there is no degree-1 nodes. In the following we call a n-tree a full binary tree with n leaves. There is a unique n-tree for n = 0 to 2.

Let t be a n-tree and let π be a path in t. The left weight (resp. right weight) a π (resp. b π ) is defined as the number of left (resp. right) edges in the path π. The length of π denoted by l π = a π + b π is the length of the path π. Denote by w π = a π b π the weight of π. The cost c π of π is the sum of its weight and its length. So we have c π = w π + l π .

Let ν be a node in t. Denote by π ν the path from the node ν to the root of t. Denote by ν l (resp. ν r ) the left child of ν (resp. the right child of ν). Denote by f ν the father4 of ν. If π is a path from the node ν to the root of t then we denote by f π the path from the node f ν to the root of t. We also associate a πν , b πν , w πν , l πν and c πν to the node ν and we denote them respectively by a ν , b ν , w ν , l ν and c ν . The set of leaves of a tree t will be denoted by L t . The weight w(t) of the tree t is defined as the sum of the weight of its leaves, that is w(t) = ν∈Lt w ν . Proposition 4. Each Z(A n ) partition system corresponds to a unique n-tree.

Proof. The idea is that if we follow the execution trace of the recursive Algorithm 2 we can see that it corresponds to a binary tree whose weight is the number of transitions of the reduced automaton. And by induction we can prove that for each state q we have a νq = a(q) and b νq = b(q). See Figure 4. So, finding a minimal Z(A n ) partition system is reduced to finding a n-tree having minimal weight. Let us denote it by Z-tree of rank n.

Let Split(t) be the function that returns the tree obtained from t by replacing a leaf having minimal cost in t by the unique 2-tree. See Figure 5.

Proposition 5. The set of Z-trees can be computed inductively as follows:

-1-tree is the Z-tree of rank one if t is a Z-tree (of rank i) then Split(t) is a Z-tree (of rank i + 1).

Proof. Let t n be a Z-tree of rank n. To get a tree t n+1 of rank n + 1 from t n we have to split a leaf µ. The weight of t n+1 is:

w(t n+1 ) = ν∈Lt n+1 w ν = ( ν∈Lt n w ν ) -w µ + w lef t-child(µ) + w right-child(µ) = w(t n ) -a µ b µ + (a µ + 1)b µ + a µ (b µ + 1) = w(t n ) + c µ
If µ is the leaf of t n which have the minimal coast, then, the tree t n+1 will have minimal weight.

So, this inductive construction allows us to have the minimal weight tree. The difference of weights between two consecutive minimal trees is exactly the cost of the split leaf. All Z-trees of rank less than n, can be generated by the following Algorithm 3.

Algorithm 3 MinZtree (n)

Require: n ∈ N Ensure: Z-tree of rank less than n 1: t ← 1-tree 2: for i = 1 to n do 3: t ← Split(t) 4: end for Theorem 2. Algorithm 3 computes one Z-tree of rank i for all i = 1 to n in O(n log n) time.

Proof. At each step of this algorithm we look for a minimal cost leaf and then we split it. We can maintain the costs of the leaves in a dynamic structure which allow us a logarithmic time search for the minimal cost leaf and also a logarithmic time insertion of the two leaves obtained from the split function.

It is clear that for a given n, there may exist several Z-trees of rank n.

In the following we introduce a subclass of full binary trees (called P -trees), for which the Z-trees are unique. We do that in order to study the size-complexity of the reduced automata (the number of transitions).

P -Trees

We denote by M t the set of leaves having minimal cost in t, that is: M t = arg min ν∈Lt c ν . The function SplitAll(t) returns the tree obtained from t by replacing every leaf in M t by the unique 2-tree. See Figure 5. Fig. 5. A 8-tree t with a Split tree and its SplitAll tree. Values in the nodes are costs.

Definition 3. The class (t n ) n>0 of P -trees is defined inductively as follows:

t 1 = 1-tree and t (n+1) = SplitAll(t n ).

See Figure 6.

Remark 2. Notice that if ν ∈ M tn then c ν = (n -1). This can be established by induction on n. Therefore, to get t (n+1) from t n we split leaves of cost (n -1).

t1 t2 t3 t4 t5
Fig. 6. The first four P -trees. Leaves with big circle have minimal cost.

Eratosthenes-Pascal's Triangle

The Eratosthenes-Pascal's Triangle is constructed from Pascal's Triangle as follows: we interleave each column k of Pascal's Triangle with (k -1) zeros. The element of the Eratosthenes-Pascal's Triangle at the n th row and the k th column is denoted by T k n with n ≥ 1 and k ≥ 1.

1 2 3 4 5 6 7 8 9 10 • • • 1 1 2 1 1 3 1 0 1 4 1 2 0 1 5 1 0 0 0 1 6 1 3 3 0 0 1 7 1 0 0 0 0 0 1 8 1 4 0 4 0 0 0 1 9 1 0 6 0 0 0 0 0 1 10 1 5 0 0 5 0 0 0 0 1 11 1 0 0 0 0 0 0 0 0 0 1 Let n be a natural number. We denote by D n the set of divisors of n. Proposition 6.

T k n =      n k + k -2 k -1 if k ∈ D n 0 otherwise.
Proof. The element of the Pascal's Triangle at the n th row and the j th column is i-1 j-1 . This element is moved in the Eratosthenes-Pascal's Triangle to the row r = (i -j + 1)j in the same column j. We have then

T j (i-j+1)j = i -1 j -1 Thus T j r = r j + j -2 j -1
Let S n be the sum of the n th row in the Eratosthenes-Pascal's Triangle.

S n = n k=1 T k n = k∈Dn n k + k -2 k -1

P -trees and Eratosthenes-Pascal's Triangle

In this section we will describe the link between the Eratosthenes-Pascal's Triangle and the set of P -trees.

Theorem 3. The sum of the elements of the n th row of Eratosthenes-Pascal's Triangle's, S n , is exactly |M tn | the number of leaves of minimal cost in the P -tree t n .

To prove this theorem we introduce a family F n which is in bijection with both the Eratosthenes-Pascal's Triangle's rows and with P -trees. We focus on the following question: Given a natural number n what are all the possible paths that have (n -1) as cost? To answer this question, we define F (n-1) as the set of all paths of cost (n -1):

F (n-1) = {π | c π = (n -1)} Lemma 2. For all n ≥ 1, |F (n-1) | = S n .
Proof. Let F i (n-1) for 0 ≤ i ≤ (n -1) be the set of paths of F (n-1) having i left edges. We have

F (n-1) = (n-1) i=0 F i (n-1)
where

F i (n-1) = {π ∈ F (n-1) | a π = i} = {π | (a π b π + a π + b π = (n -1)) ∧ (a π = i)} = {π | (b π = n i + 1 -1) ∧ (a π = i)} So, we get |F i (n-1) | =      n i + 1 + i -1 i if (i + 1) ∈ D n 0 otherwise.
This corresponds to the different ways to arrange i left edges in a path of length

n i + 1 + i -1. Let k = i + 1 then |F k-1 (n-1) | =      n k + k -2 k -1 if k ∈ D n 0 otherwise. Thus |F k-1 (n-1) | = T k n . That is for 0 ≤ i < n we get |F i (n-1) | = T i+1 n .
Finally

|F (n-1) | = (n-1) i=0 |F i (n-1) | = (n-1) i=0 T i+1 n = n i=1 T i n = S n
Therefore, we can associate the n th row of the Eratosthenes-Pascal's Triangle to Proof. We will show that each set F n is associated with the P -tree t n . We proceed by induction. Assume that for all k ≤ n, the P -tree t n contains all paths of cost less or equal to n, that is, for all k ≤ n, F k is in the P -tree t n . From this hypothesis we should show that t (n+1) = SplitAll(t n ) contains all the paths of cost less or equal to (n + 1). Of course, t (n+1) contains all the paths of cost less or equal to n because the tree t (n+1) is obtained from t n . However, does t (n+1) contain all paths of cost equal to (n + 1)? To answer this question, we proceed by absurd. Let π be a path with cost c π = (n + 1) which is not in t (n+1) . It is clear that the cost of the father of π is:

F (n-1) . (a) (b) (c) (d)
c fπ = (a π -1)b π + (a π -1) + b π if π is a left child, a π (b π -1) + a π + (b π -1) if π is a right child.
As c fπ < c π then c fπ ≤ n. From this we deduce that when c fπ < n, the father f π of π was split by the SplitAll() function within an earlier or it will be split in the current tree t n in the case where c fπ = n (see Remark 2). In both cases, the path π is necessarily in t (n+1) .

From the two last lemmas, we construct a bijection between the P -trees and the Eratosthenes-Pascal's Triangle's rows, and we claim following corollary:

Corollary 1. Let t n be a P -tree. Then |M tn | = |F (n-1) | = S n .
From Corollary 1, we have

S n = k|n n k + k -2 k -1 .
Proposition 7. Let t n be a P -tree. Its size s(t n ) = |L tn | (the number of leaves), and its weight w(t n ) are:

s(t n ) = 1 + n-1 i=1 S i and w(t n ) = n i=2 (i -2)S i-1 .
Proof. From Corollary 1 and Remark 2, we have, the size of the tree t n is the sum of the size of the tree t n-1 and the number of all split minimal leaves, that is, s(t n ) = s(t n-1 ) + |M tn-1 |. We have also, the weight of the tree t n is the sum of the weight of the tree t n-1 and the costs of all split minimal leaves, that is,

w(t n ) = w(t n-1 ) + (n -2)|M tn-1 |.
From Proposition 4, the P -tree t n (which is a Z-tree) corresponds to a Z(A s(tn) ) partition system. The CFS automaton associated with Z(A s(tn) ) has w(t n ) transitions.

With our construction the CFS automaton associated with a Z(A n ) CFS partition system may contain several initial states. However, in order to compare the number of minimal automata and their number of transitions, with those obtained by R. Cox [START_REF] Cox | Minimal number of edges in e-free non-deterministic finite automata (nfa) for regular expression (1 + ε) • (2 + ε) • (3 + ε) • • • (n + ε). The On-Line Encyclopedia of Integer Sequences[END_REF] (seen Table 1), we must restrict our study to CFS partition systems leading to a unique initial state automata. It is easy to verify that in this case a P -tree t n has s(t [START_REF] Cox | Minimal number of edges in e-free non-deterministic finite automata (nfa) for regular expression (1 + ε) • (2 + ε) • (3 + ε) • • • (n + ε). The On-Line Encyclopedia of Integer Sequences[END_REF] (ii) Number of minimal automata estimated by R. Cox [START_REF] Cox | Minimal number of edges in e-free non-deterministic finite automata (nfa) for regular expression (1 + ε) • (2 + ε) • (3 + ε) • • • (n + ε). The On-Line Encyclopedia of Integer Sequences[END_REF] (iii) Number of transitions in our reduced automaton (iv) Number of reduced automata estimated by our approach.

n ) = 1 + n-1 i=1 S i+2 2 and w(t n ) = n i=1 iS i+1 2

Asymptote behavior of the number of transitions

Appendix A contains the basics of asymptotic notations. In this section we shall establish one of our main result which concern the asymptotically result on the behavior of the number of transitions w(t n ), where s(t n ) is the number of states. Namely, we will show that the weight of our automata is asymptotically equivalent to s(t n )log 2 s(t n ) up to constant which means that the number of transitions is minimal in the sense that we reach the lower bounded of Shnitger [START_REF] Schnitger | Regular expressions and nfas without ǫ-transitions[END_REF]. Indeed, we have

Theorem 4. log 2 (4) w(t n ) ∼ s(t n ) log 2 s(t n ).
As a consequence we obtain the following Corollary 2. For a large n, we have

ω(t n ) < s(t n ) log 2 s(t n ).
It is also easy to deduce the following

Corollary 3. s(t n ) log 2 s(t n ) log log s(t n ) = o w(t n ) .
Before starting the proof of Theorem 4. Observe that according to the Corollary 1 combined with Proposition 7, one may consider that ω(t n ) and s(t n ) are given by

s(t n ) = 1 + n-1 1 S i and ω(t n ) = n-2 i=1 iS i+1 .
As usual in the number theory, any arithmetical function f : Z -→ R can be extended to the real line by putting, for any x ∈ R, f (x) = f (⌊x⌋). Therefore, for any x ≥ 2, we have

s(t x ) = 1 + ⌊x⌋-1 1 S i and ω(t x ) = ⌊x⌋-2 i=1 iS i+1 .
But the function x ∈ [3, +∞[ → x log(x) is increasing function. It follows that we have, For all u ≥ 3,

s(t u ) ≥ 1 3 u -2 log(u -2) ,
Which achieve the proof of the proposition.

For any x ≥ 2 and any positive integer n, let

Li n (x) = ˆx 2 dt log n (t)
.

Let us summarize in the following proposition a classical well-known results on Li n (x) that we shall used.

Proposition 9. For every x ≥ 2 and integer n ≥ 1, we have

Li 1 (x) = x log(x) + Li 2 (x) - 2 log(2) , (4)
Li n (x) = O x log n (x) . (5)
Now, we are able to formulate our key estimation of ω(t n ) with respect to s(t n ) in the following proposition.

Proposition 10. lim sup ω(t n ) (n -2) s(t n ) ≤ 1.
Proof. Applying Lemma 4 we have, for any n ≥ 3,

ω(t n ) (n -2) s(t n ) = 1 - 3 (n -2) s(t n ) - 1 (n -2) s(t n ) ˆn-2 2 s(t u )du.
Therefore, for any n ≥ 5, we have

ω(t n ) (n -2) s(t n ) ≤ 1 - 3 (n -2) s(t n )
From Proposition 8 we deduce, that for any n ≥ 5 we have

3 (n -2) s(t n ) ≤ 1 log(n -2)
Hence, by letting n goes to ∞ we obtain

3 (n -2) s(t n ) -----→ n→∞ 0.
Which implies that

lim sup ω(t n ) (n -2) s(t n ) ≤ 1,
and this finish the proof of the proposition.

We are reduced to compare the sequences (n -2)s(t n ) and s(t n ) log 2 s(t n ). For that we shall estimate s(t x ). Precisely, we argue that we have the following Theorem 7. For a large x > 0 we have

x 3 4 4 √ x-1 √ π ≤ s(t x ) ≤ x 3 4 log(x) 4 √ x-1 √ π .
The proof of Theorem 7 will be given later. For instance, using Theorem 7 holds we shall extended Proposition 10 as follows.

Proposition 11. The sequences ω(t n ) and (n -2)s(t n ) two sequences be equivalent. That is,

ω(t n ) (n -2)s(t n ) -----→ n→∞ 1.
Proof. By Lemma 4, write

ω(t n ) (n -2)s(t n ) = 1 - 3 (n -2)s(t n ) - 1 (n -2)s(t n ) ˆn-2 2 s(t x )dx.
Let ε > 0 and x sufficiently large. Then, by Theorem 7, for a large x, we have

x 3 4 4 √ x-1 √ π ≤ s(t x ) ≤ x 3 4 log(x) 4 √ x-1 √ π . (6)
But ˆn-2

2 4 √ x dx u= √ x = ˆ√n-2 √ 2 2 u 4 u du = 2u log(4) 4 u √ n-2 √ 2 - 2 log 2 (4) 4 u √ n-2 √ 2 = 2 √ n -2 .4 √ n-2 log(4) - 2 √ 2 .4 √ 2 log(4) - 2 .4 √ n-2 log 2 (4) + 2 .4 √ 2 log 2 (4) . (7) Since 1
(n-2)s(tn) vanishes at the infinity and s(t n ) is equivalent to s(t n ) we may assume that (6) holds starting from 2 and Theorem 7 is valid for s(t n ). Therefore

1 (n -2) s(t n ) ˆn-2 2 s(t x )dx ≤ 1 (n -2)(n -2) 3 4 4 √ n ˆn-2 2 x 3 4 log(x) 4 √ x dx ≤ (n -2) 3 4 log(n -2) (n -2)(n -2) 3 4 4 √ n ˆn-2 2 4 √ x dx ≤ log(n -2) (n -2)4 √ n ˆn-2 2 4 √ x dx. (8)
From ( 7) combined with (8) it follows that

1 (n -2) s(t n ) ˆn-2 2 s(t x )dx ≤ 2 √ n -2 .4 √ n-2 log(4) × log(n -2) (n -2)4 √ n + 2 .4 √ 2 log 2 (4) × log(n -2) (n -2)4 √ n -----→ n→∞ 0. We conclude that ω(t n ) (n -2) s(t n ) -----→ n→∞ 1.
which proves the proposition.

It remains to prove Theorem 7. For that we shall need the following classical lemma. The proof of it can be found in [START_REF] Flajolet | Analytic Combinatorics[END_REF]. Nevertheless we include the proof for the sake of completeness.

Lemma 5. 2n n ∼ 4 n √ π √ n .
Proof. By Stirling formula we have

n! ∼ n n e -n √ 2πn
Hence 2n n = 2n! (n!) 2 ∼ 4 n √ π √ n .
This finishes the proof of the lemma.

Proof (Proof of Theorem 7). For x ≥ 2, Write

s(t x ) = n≤x d|n n d + d -2 d -1 = dq≤x q + d -2 d -1 = ⌊x⌋ d=1 ⌊ x d ⌋ q=1 q + d -2 d - 1
From this we see that

s(t x ) ≤ ⌊x⌋ d=1 ⌊ x d ⌋ 2(⌊x⌋ -1) ⌊x⌋ -1 ≤ x log(x) 2(⌊x⌋ -1) ⌊x⌋ -1 , (9)
and

s(t x ) ≥ ⌊x⌋ 2(⌊x⌋ -1) ⌊x⌋ -1 , (10)
Using the relation ⌊x⌋ = x + O(1) combined with ( 9) and ( 10), we obtain

x 2(⌊x⌋ -1) ⌊x⌋ -1 ≤ s(t x ) ≤ x log(x) 2(⌊x⌋ -1) ⌊x⌋ -1 .
By Lemma 5, this gives

x 3 4 4 √ x-1 √ π ≤ s(t x ) ≤ x 3 4 log(x) 4 √ x-1 √ π ,
which proves the theorem.

Now we are able to give the proof of Theorem 4.

Proof (of Theorem 4). By Proposition 11, it is sufficient to show that

s(t n ) log 2 (s(t n )) ∼ (n -2)s(t n ).
For that, observe that we have

s(t n ) log 2 (s(t n )) (n -2)s(t n ) = log 2 (s(t n )) n -2
Applying Theorem 3, we deduce that

log(s(t n )) ∼ log(4) √ n. Whence log 2 (s(t n )) ∼ log 2 (4) n. Hence log 2 (s(t n )) (n -2) -----→ n→∞ log 2 (4). We deduce that ω(t n ) s(t n ) log 2 s(t n ) -----→ n→∞ 1 log 2 (4) < 1.
This finishes the proof of the theorem.

Conclusion

In this paper we show how binary trees can be used to design a fast algorithm for computing an automaton with a reduced5 number of transitions recognizing the language L(E n ). We have verify that our algorithm gives the minimal number of transitions for n = 1 to 7 (see Table 1) and we have shown that our reduction is asymptotically a minimization. Hence, we conjecture that Algorithm 3 computes the minimal transition automaton.

Z >0 {+∞} with s 0 = +∞, S = R with s 0 any point in R; S = C or a subset of C with s 0 = 0, and so on. Two functions f and g from S \ {s 0 } to R or C are given.

-O-notation: write f (s) One also says that "f is of order at most g", or "f is big-Oh of g"(as s tends to s 0 ). -o-notation: write f (s) s→s0

= o(g(s)),

if the ratio f (s) g(s)

tends to 0 as s → s 0 in S. In other words, for any (arbitrarily small) ε > 0, there exists a neighborhood V ε of s 0 (depending on ε), such that |f (s)| < ε|g(s)|, s ∈ V ε , s = s 0 .

One also says that "f is of order smaller than g, or f is little-oh of g" (as s tends to s 0 ). -∼-notation: write f (s) s→s0 ∼ g(s), if the ratio f (s) g(s)

tends to 1 as s → s 0 in S. One also says that "f and g are asymptotically equivalent" (as s tends to s 0 ). -Ω-notation: write f (s) s→s0

= Ω(g(s)), if the ratio f (s) g(s) stays bounded from below in modulus by a non-zero quantity, as s → s 0 in S. Which means that there exists k > 0 and a neighborhood V of s 0 , such that f (s) ≥ k.g(s), s ∈ V.

One then says that f is of order at least g. -θ-notation: if f (s) = O(g(s)) and f (s) = Ω(g(s)), write f (s) s→s0 = θ(g(s)). This implies that there exits k, C > 0 and a neighborhood V of s 0 , such that k.g(s) ≤ f (s) ≤ C.g(s), s ∈ V.

One then says that f is of order exactly g.

At this point we are able to make a parallel between the history of our contribution and the history of the famous Prime Number Theorem (PNT) which we shall use later in its weaker form. The PNT Theorem concerns the asymptotic behavior of the prime-counting function π(x) = |{p ≤ x, p prime}|. Using asymptotic notation the PNT can be restated as π(x) ∼ x ln x .

The behavior of π(x) has been the object of intense study by many celebrated mathematicians ever since the eighteenth centry. Inspection of tables of primes led Gauss (1792) and Legendre (1798) to conjecture the PNT. In 1808 Legendre published the formula π(x) = x/(logx + A(x)), where A(x) tends to a constant B = -1.08366 as x -→ +∞, which means that π is Ω(x/log(x)).

According to Bateman and Diamond [START_REF] Harold | A hundred years of prime numbers[END_REF], The first person to establish the true order of π(x) was P. L. Chebyshev. Indeed, in two papers from 1848 and 1850, Chebychev prove that π(x) is θ(x/log(x). This result is known in nowadays as Chebychev Theorem. Finally, in 1896 the PNT was first proved by Hadamard and de la Vallée Poussin. Their proofs were long and intricate. A simplified modern presentation is given on pages 41-47 of Titchmarsh's book on the Riemann Zeta function [START_REF] Titchmarsh | The theory of the Riemann zeta-function[END_REF].

7 Fig. 1 .

 71 Fig. 1. The reduced automaton (at right) is obtained from a decomposition transition table (at right bottom) of the homogeneous automaton (at left). The automaton A(L5) is the part of the reduced automaton which represents the triangle (in the transition table decomposition).

Fig. 2 .

 2 Fig. 2. We have f ollow(0) = {1, 2, 3, #}. Here are three possible decompositions of f ollow(0). The two first ones are partition decompositions. (i) dec(0) = {{1, 2}, {3, #}} (ii) dec(0) = {{1}, {2}, {3, #}} (iii) dec(0) = {{1, 2}, {2, 3, #}} .

Fig. 3 .

 3 Fig. 3. Two CFS Automata constructed from the automaton A3 shown in Figure 2.

Fig. 4 .

 4 Fig. 4. A full binary 5-tree and an associated Z(A5) partition (left edges are represented by dotted lines and right edges with solid lines).

Fig. 7 . 3 1 = 3 .Lemma 3 .

 7333 Fig. 7. All paths of cost 5. (a) The unique path of cost 5 with no left edges which corresponds to 5 0 = 1. (b) The three paths of cost 5 with one left edges which corresponds to 3 1 = 3. (c) The three paths of cost 5 with two left edges which corresponds to 3 2 = 3. (d) The unique path of cost 5 with five left edges which corresponds to 5 5 = 1. There is no paths of cost 5 with three or four left edges that correspond to zeros in the 6 th row of the Eratosthenes-Pascal's Triangle.

s→s0=

  O(g(s)), if the ratio f (s) g(s) stays bounded as s → s 0 in S. In other words, there exists a neighborhood V of s 0 and a constant C > 0 such that |f (s)| < C|g(s)|, s ∈ V, s = s 0 .

Table 1 .

 1 . Comparison table. (i) Minimal transition number estimated by R. Cox

	n 1	2	3	4	5	6	7	8	9	10	11	12	13	14
	(i) 1	3	6	9	13	18	23	≤ 28 ≤ 34 ≤ 41 ?	?	?	?
	(ii) 1	1	2	1	1	4	6	≥ 1 ≥ 1 ≥ 1 ?	?	?	?
	(iii) 1	3	6	9	13	18	23	28	33	39	46	53	60	67
	(iv) 1	1	2	1	1	4	6	4	1	1	5	10	10	5

The first ancestor.

Asymptotically minimal.

Acknowledgments. Special thanks to Saïd Abdeddaïm, Alexis Bès, Patrick Cégielski, Jean-Marc Champarnaud and Yuri Matiyasevich.

We recall that the classical arithmetical function π(x) denote the number of primes not exceeding x. We shall need also the following classical identity due to Abel Theorem 5 (Abel's identity, [START_REF] Tom | Introduction to analytic number theory[END_REF]). For any arithmetical function a(n) let

where 0 < y < x. Then we have

We deduce easily from the Abel's identity the following lemma Lemma 4. For any integer n ≥ 1, we have

where s(t u ) = s(t u ) -1.

We need to estimate ω(t n ) with respect to s(t n ). For that, we shall need the following weaker form of the Prime Number Theorem (WPNT for short) due to Cheyshev

Theorem 6 ([1]

). For every integer n ≥ 2 we have

We deduce from the WPNT the following crucial proposition Proposition 8. For all u ≥ 4 we have

Proof. We can show that for any prime number p ≥ 2, we have S p = 2. Hence

Consequently, by WPNT we get

.