Dominique Pastor
email: dominique.pastor@enst-bretagne.fr

Asmaa Amehraye
email: asmaa.amehraye@enst-bretagne.fr

FROM NON-PARAMETRIC STATISTICS TO SPEECH DENOISING

Given some signal additively corrupted by independent white Gaussian noise with unknown standard deviation σ, we present a new estimator of σ. This estimator derives from a theoretical result presented and commented in the paper. Without any preliminary signal detection, the estimate is performed on the basis of the time-frequency components returned by a standard spectrogram where the Discrete Fourier Transform is simply weighted by the square window. No assumption about the signal statistics is made. The signal time-frequency components are assumed to have probabilities of presence less than or equal to one half.

This estimator is suited to speech denoising. It avoids the use of any Voice Activity Detector and is an alternative solution to subspace approaches. Objective performance measurements show that the standard Wiener filtering of speech signals can be tuned with the outcome of this estimator without a significant loss in comparison with the measurements obtained when the noise standard deviation is known.

MOTIVATION

Let s[t], t = 1, . . . , T be the samples of some speech signal and suppose that these T samples are corrupted by additive and independent stationary noise x[t], t = 1, 2, . . . , T so that the samples of the observed signal are

y[t] = s[t] + x[t], t = 1, . . . , T. (1)
We assume that noise is white and Gaussian with null mean and standard deviation σ: for every t ∈ {1, 2, . . . , T }, x[t] ∼ N (0, σ 2 ).

The Wiener filtering of the noisy speech signal y requires prior knowledge of the noise standard deviation σ. A basic and popular solution consists in using a Voice Activity Detector (VAD): the estimate of σ is the square root of the Maximum Likelihood Estimate (MLE) computed on the basis of the samples of the time frames that the VAD has detected as noise alone. Subspace approaches can also be used to estimate σ by computing the smallest eigenvalues of the noisy speech autocorrelation matrix; the model order is difficult to choose and the computation of the eigenvalues may prove unstable. This paper proposes a new estimator of the noise standard deviation. The theoretical foundation of this estimator is proposition 3.1, stated in section 3 after exposing some preliminary material in section 2. This theoretical result is non-parametric in the sense that it makes no assumption about the probability distributions of the signals and assumes neither that these signals are identically distributed nor that they have equal probabilities of presence.

Section 4 then presents the estimator deriving from proposition 3.1 and a few preliminary experimental results. This estimator is then employed in section 5 to adjust the standard Wiener filtering of the noisy speech signal y introduced above. According to objective performance measurements, the denoised speech signals are not significantly more distorted than those achieved when the filtering is tuned with the exact value of the noise standard deviation.

Conclusions and perspectives are given in section 6.

PRELIMINARY MATERIAL

The random vectors and variables are supposed to be defined on the same probability space denoted by (Ω, M, P ) and for every element ω ∈ Ω. As usual, if property P holds true almost surely, we write P (a-s). Given a positive real value σ, a sequence X = (X k ) k∈N of random complex variables is said to be a complex white Gaussian noise (CWGN) with standard deviation σ if the random variables X k , k = 1, 2, . . ., are complex, mutually independent and identically Gaussian distributed with null mean and variance σ 2 . The real and imaginary parts ℜX k and ℑX k of each X k form a two-dimensional random vector such that (ℜX k ℑX k ) ∼ N 0, (σ 2 /2)I where I stands for the 2 × 2 identity matrix.

The minimum amplitude a(S) of a sequence S = (S k ) k∈N of random complex variables is defined by

a(S) = sup {α ∈ [0, ∞] : ∀k ∈ N, |S k | ≥ α (a-s)} . (2)
If f is some map of the set of all the sequences of complex random variables into R, we say that the limit of f is ℓ ∈ R when a(S) tends to ∞ and write that lim a(S)→∞ f (S) = ℓ if, for any positive real value η, there exists some A 0 ∈ (0, ∞) such that, for every A ≥ A 0 and every S such that a(S) ≥ A, |f (S) -ℓ| ≤ η.

The set L 2 (Ω, C) stands for the set of those complex random variables Y such that E[|Y | 2 ] < ∞. We then define ℓ ∞ (N, L 2 (Ω, C)) as the set of those sequences S = (S k ) k∈N of complex random variables such that

S k ∈ L 2 (Ω, C) for every k ∈ N and sup k∈N E[|S k | 2 ] is finite.
Given a random variable Y and a real number τ , I(Y ≤ τ ) stands for the indicator function of the event {Y ≤ τ }. As usual, I 0 is the zeroth-order modified Bessel function of the first kind. Throughout the rest of the text, we say 'independent' instead of 'mutually independent' for brevity.

A THEORETICAL RESULT

Proposition 3.1 stated below is a corollary of a more general theorem established in [START_REF] Pastor | On the detection of signals with unknown distributions and priors in white gaussian noise[END_REF]. As an introduction to proposition 3.1, we begin with an intuitive approach. It makes the reader understand the main ideas behind proposition 3.1 and the results given in [START_REF] Pastor | On the detection of signals with unknown distributions and priors in white gaussian noise[END_REF].

Consider a sequence Y = (Y k ) k∈N of complex random variables where each Y k is either the sum of some signal S k and noise X k or noise X k alone. We assume that X = (X k ) k∈N is a CWGN with standard deviation σ. For every given k ∈ N, the presence of noise alone is the null hypothesis whereas the presence of some signal in noise is the alternative one. We assume that, for every k ∈ N, the index of the true hypothesis is a random variable ε, valued in {0, 1} and independent with S k and X k . We thus can write that Y k = ε k S k + X k . The a priori probabilities of presence and absence of the signal Λ k are then P ({ε k = 1}) and P ({ε k = 0}), respectively. Proposition 3.1 significantly reduces the importance of the choice of these probabilities since they will be assumed to be upper-bounded.

At this stage, assume that the random variables ε k , k ∈ N, are independent and identically distributed (iid) as well as the random signals S k , k ∈ N. It follows that the random vectors Y k , k ∈ N, are iid as well. Let p stand for the common value of the probabilities of presence P ({ε k = 1}). We assume that p ≤ 1/2.

Given some real number T , set

A m (T ) = 1 m m k=1 |Y k |I(|Y k | ≤ T )
According to Kolmogorov's classical strong limit theorem,

lim m→∞ A m (T ) = E [|Y k |I(|Y k | ≤ T )] (a-s) ( 3
)
where k is any element of {1, . . . , m}. An easy computation shows that we can write that

E [|Y k |I(|Y k | ≤ T )] = (1-p)E [|X k |I(|X k | ≤ T )] × 1+ p 1 -p E [|S k +X k |I(|S k +X k | ≤ T )] E [|X k |I(|X k | ≤ T )] . (4)
Let A be a lower bound for the amplitudes of the signals S k , k ∈ N. If A is large enough in comparison with σ, we can reasonably expect the existence of some threshold T that makes it possible to distinguish noisy signals from noise alone with a rather small probability of error. As a matter of fact, regarding the choice for T , we can be very specific as follows.

For any given non negative real number h, let T τ stand for the map defined for every complex value z by

T τ (z) = 1 if |z| ≥ τ 0 otherwise. (5)
Clearly, for every k ∈ N, T τ is a statistical test for making a decision on the value of ε k since the composite map

T τ (U k ) is measurable and T τ (U k ) ∈ {0, 1}.
In what follows, T τ is called the thresholding test with threshold height τ . The error probability of this test is then the probability

P ({T τ (U k ) = ε k }) of the event {T τ (U k ) = ε k }.
This error probability does not depend on k since the observations U k , k ∈ N, are assumed to be iid. This is the reason why we simply denote it by P e {T τ }, without mentioning the observation under consideration.

Since the probability of presence of any signal S k , k ∈ nN, is assumed to be less than or equal to 1/2 and by taking into account that (see [1, Eq. 9.6.47, p. 377]) I 0 (x) = 0 F 1 (1; x 2 /4), [4, Theorem VII.1] tells us the following. For every x ∈ R, set

κ(x) = I -1 0 (e x 2 )/2x (6)
with κ(0) = 1; for making a decision on the value of ε k where k is any natural number, the error probability P e {T σκ(A/σ) } of the thresholding test T σκ(A/σ) with threshold height σκ(A/σ) is less than or equal to Q(A/σ) where, for any given non-negative real number x,

Q(x) = e -x 2 κ(x) 0 e -t 2 tI 0 (2xt)dt + 1 2 e -κ(x) 2 . ( 7
)
We thus can write that P e {T σκ(A/σ) } ≤ Q(A/σ). In equation (4), set now T = σh with h = κ(A/σ). The function Q(x) decreases very rapidly when x increases. Hence, for large values of A, the probabilities P ({|X k | > σh}) and P ({|S k + X k | ≤ σh}) are small and, thus, the expectation

E [|S k + X k |I(|S k + X k | ≤ σh)]
can reasonably be expected to be significantly smaller than

E [|X k |I(|X k | ≤ σh)].
Since p is assumed to be less than or equal to one half, p/(1p) is less than or equal to 1. Consequently, in a certain sense to specify, we should be able to prove that

E [|Y k |I(|Y k | ≤ σh)] ≈ (1 -p)E [|X k |I(|X k | ≤ σh)] .
Without caring about mathematical exactness, we combine this approximation to the almost surely convergence of equation ( 3) to obtain that, in a certain sense,

A m (σh) ≈ (1 -p)E [|X k |I(|X k | ≤ σh)] (8
B m (σh) ≈ (1 -p)E [I(|X k | ≤ σh)] . (9)
Consider now the ratio A m (σh)/B m (σh). This ratio makes it possible to get rid of the unknown prior p.

Moreover, since X k ∼ N c (0, σ 2 ), the distribution of |X k | is known and its density f (x) is that of the square of a Rayleigh distributed variable. Taking into account that the variance of the real and imaginary parts of X k both equal σ 2 /2, this density is given by :

f (x) = (2x/σ 2 )e -x 2 /σ 2 if x ≥ 0, 0 otherwise. ( 10
)
Therefore, we easily obtain that

E [|X k |I(|X k | ≤ σh)] = 2σ h 0 t 2 e -t 2 dt and that E [I(|X k | ≤ σh)] = P ({X k ≤ σh}) = 1-e -h 2 .
According to these equalities, equations ( 3) and ( 8), we conclude that

A m (σh)/B m (σh) ≈ 2σ h 0 t 2 e -t 2 dt/(1 -e -h 2 ). (11)
Once again, this approximation must be understood with respect to a certain convergence criterion. This one is introduced in proposition 3.1. Its more general form is given in [START_REF] Pastor | On the detection of signals with unknown distributions and priors in white gaussian noise[END_REF]. As a matter of fact, the same type of intuitive approach as that presented above can be used to guess part of the results established in [START_REF] Pastor | On the detection of signals with unknown distributions and priors in white gaussian noise[END_REF]. Proposition 3.1 and its extension not only specify the exact meaning of equation ( 11) but also significantly extend the conditions of validity of (11) because they state that the convergence holds true even for non iid signals and non iid priors. Given any natural number m and any pair (x, T ) of positive real numbers, define the random variable D m (x, T ) by

Proposition 3.1 Let Y = (Y k ) k∈N be a sequence of com- plex random variables such that, for every k ∈ N, Y k = ε k S k + X k where S = (S k ) ∈ ℓ ∞ (N, L 2 (Ω, C)), X = (X k ) k∈N
D m (x, T ) = m k=1 |Y k |I(|Y k | ≤ xT ) m k=1 I(|Y k | ≤ xT ) -2x T 0 u 2 e -u 2 du 1 -e -T 2 .
Then, the standard deviation σ is the unique positive real number x such that, for every

β 0 ∈ (0, 1], lim a(S)→∞ lim m D m (x, βκ(a(S)/x)) ∞ = 0 (12) uniformly in β ∈ [β 0 , 1] where, for every x ∈ R, κ(x) = I -1 0 (e x 2 )/2x (13)
with κ(0) = 1.

A NEW ALGORITHM FOR ESTIMATING THE NOISE STANDARD DEVIATION

On the basis of proposition 3.1, we start by introducing a discrete cost. A minimum of this discrete cost can be computed and considered as a first estimate of the noise standard deviation. This estimate will be called the Essential Supremum Estimate of type I (ESE-I) because of the crucial role played by the essential supremum norm in its computation. The term ESE-I will also stand for the estimator itself.

Experimental results aimed at assessing the ESE-I suggest another estimate of the noise standard deviation. This new estimate is hereafter called the Essential Supremum Estimate of type II (ESE-II). The term ESE-II will also designate the estimator itself. According to Monte-Carlo experiments of the same type as those mentioned above, the ESE-II performs better than the ESE-I.

The ESE-I

Let L be some natural number and set β ℓ = ℓ/L, ℓ = 1, 2, . . . , L. Suppose that A is some known lower bound for the amplitudes of the signal. We thus have a(S) ≥ A. These new notations are kept hereafter with the same meaning.

Consider when x ranges over a suitable search interval. However, in practice, no lower bound for the amplitudes of the signals is known. Surprisingly enough since 3.1 states that the larger A the better the estimate, the experimental results presented in [START_REF] Pastor | Un théorème limite et un test pour la détection non paramétrique de signaux dans un bruit blanc gaussien de variance inconnue[END_REF] and [START_REF] Pastor | On the detection of signals with unknown distributions and priors in white gaussian noise[END_REF] suggest that the asymptotic condition on the minimum amplitude of the signals can be relaxed significantly. Therefore, we consider the trivial lower bound A = 0 and the discrete cost we minimize is then

m observations Y 1 , Y 2 , . . . , Y m . If
sup ℓ∈{1,...,L}            m k=1 |Y k |I(|Y k | ≤ xβ ℓ ) m k=1 I(|Y k | ≤ xβ ℓ ) -2x β ℓ 0 u 2 e -u 2 du 1 -e -β 2 ℓ            , (15
) straightforwardly derives from (14) with A = 0 and seeing that κ(0) = 1. Any possibly local minimum σ of (15) can be considered as an estimate of the noise standard deviation. Because of the crucial role played by the essential supremum norm in proposition 3.1, σ will be called the Essential Supremum Estimate of type I (ESE-I).

To compute the ESE-I, we choose L = m as a reasonable trade-off between the expected accuracy of the estimate and the computational load incurred by the minimization routine. However, a better choice can certainly be thought up. This will be made elsewhere. where Q is some value in (0, 1), close to 1 but less than or equal to 1 -m 4(m/2-1) 2 . The reasons of this choice for the search interval are given in [START_REF] Pastor | Un théorème limite et un test pour la détection non paramétrique de signaux dans un bruit blanc gaussien de variance inconnue[END_REF] and [START_REF] Pastor | Estimating the standard deviation of some additive white Gaussian noise on the basis of non signal-free observations[END_REF].

Accuracy of the ESE-I

Let k be some natural number and L k stand for the Minimum-Probability-of-Error (MPE) test ([8, section II.B]) for making a decision on the value of ε k . The null hypothesis is thus ε k = 0 and the alternative one is ε k = 1. For the decision problem under consideration, the likelihood ratio test L k guarantees the smallest possible probability of error amongst all the possible binary hypothesis tests.

Given Y 1 , Y 2 , . . . , Y m , the test I(| • | ≤ σκ(A/σ)) simply consists in substituting the estimate σ to the exact value σ in the expression of T σκ(A/σ) . It assigns the value 1 to any complex value z whose modulus less than or equal to σκ(A/σ) and 0 otherwise. This test is not, strictly speaking, a thresholding test in the sense given above for its "thresholding height" is the random variable σκ(A/σ). However, with a slight abuse of language, we denote it by T σκ(A/σ) .

If σ is a reasonably good estimate of σ, the performance of T σκ(A/σ) can be expected to approach that of the thresholding test T σκ(A/σ) . In other words, the use of the estimate σ instead of the true value σ should not induce a significant performance loss even when the minimum amplitude A is known, provided, of course, that m is large enough. In particular, when the signals S k , k ∈ N, are independent, have their probabilities of presence all equal to 1/2 and are such that S k = Ae iΦ k where Φ k is uniformly distributed in [0, 2π], the error probability of the test T σκ(A/σ) equals Q(A/σ) ( [START_REF] Pastor | A Sharp Upper-Bound for the Probability of Error of the Likelihood Ratio Test for Detecting Signals in White Gaussian Noise[END_REF]); therefore, the error probability of the test T σκ(A/σ) should be close to Q(A/σ) when A and m are both large. Even though the computation of the error probability of the test T σκ(A/σ) is an open issue, this intuitive claim can easily be verified via Monte-Carlo simulations aimed at comparing the Binary Error Rate (BER) of this test to Q(A/σ). To achieve this simulations, we follow the standard experimental protocol adopted by practitioners in telecommunication systems.

Fix σ = 1. We carry out independent trials of m observations each by considering a number J of successive independent random copies of the observations Y 1 , . . . , Y m . These copies are henceforth denoted by Y j,1 , Y j,2 , . . . , Y j,m , j = 1, 2, . . . , J. Of course, they are constructed by using independent random copies ε j,k , S j,k and X j,k of ε k , S k and X k respectively. For every copy j and every given k ∈ {1, . . . , m}, we thus have Y j,k = ε j,k S j,k + X j,k .

For each j = 1, 2, . . . , J, let σj be the ESE-I of σ obtained during the jth trial, that is on the basis of Y j,k , k = 1, . . . , m; denote by n j the number of errors made by the test T σκ(A/σj ) = I(| • | ≤ σj κ(A/σ j ) applied to the m observations Y j,k , k = 1, . . . , m.

Since the decision is made on the same observations as those used to estimate the noise standard deviation, the accuracy of the estimate affects m decisions at one go. To reduce this effect, we proceed as practitioners in telecommunication systems usually do by fixing a minimum number J min of trials to achieve and a minimum number N min of errors to obtain during the experiments.

Trials are thus carried out until the total number J of trials is larger than or equal to J min and the total number of errors J j=1 n j obtained after these J trials is larger than or equal to N min . The BER of the test T σκ(A/σ) is then defined as the ratio J j=1 n j /(J × m). the results presented below were achieved with a minimum number of trials equal to J = 150 and a minimum number of errors equal to N min = 400.

Figure 1 displays the BER of the test T σκ(A/σ) for different values of A and m in comparison with the theoretical value Q(A/σ) of the probability of error. Table 1 gives the empirical mean and empirical Mean Square Error (MSE) of the ESE-I obtained during these experiments.

The ESE-II

With the same notations as those used so far, let Ψ m be the random variable defined by

Ψ m = 1 σ m k=1 |Y k |I(|Y k | ≤ σ) / m k=1 I(|Y k | ≤ σ) .
The empirical mean and standard deviation of Ψ m were computed during the experiments described in the previous section.

σ = 1 K m k=1 |Y k |I(|Y k | ≤ σ) / m k=1 I(|Y k | ≤ σ)
(16) where K = 0.7096 is the average value of the empirical means of Ψ m for m = 100, 200, . . . , 1000. This new estimate σ is called the ESE-II.

We conducted the same type of experiments as those presented in section 4.2. The BERs obtained when the noise standard deviation is estimated by the ESE-II are then those of figure 2. The empirical mean and empirical MSE of this estimate are given in table 3. According to these results, the ESE-II is more accurate than the

APPLICATION TO SPEECH ENHANCEMENT

With the same notations and under the same assumptions as those of section 1, we use the ESE-II to estimate the noise standard deviation and adjust the Wiener filtering of the noisy speech signal y.

Standard deviation estimation via the ESE-II

We split the T available samples y[t], t = 1, 2, . . . , T , into non-overlapping frames of N = 2 q successive samples each. As usual, q is chosen so that N F s ≈ 20ms where F s is the sampling frequency. Let K stand for the number of frames such constructed. The kth frame is then the finite sequence of samples y[(k-1) N 2 +n], n = 0, 1, . . . , N -1. The N -Discrete Fourier Transform (DFT) of this frame is then the sequence Y k,ℓ , n = 0, 1, . . . , N -1, with

Y k,ℓ = C N -1 t=0 y[(k -1)N + t]e -i2πℓt/N , (17)
C being some constant, usually chosen in {1, 1/N, 1/ √ N }. We thus obtain the matrix [Y k,ℓ ] k∈{1,...,K},ℓ∈{0,...,N -1} . Because of the Hermitian symmetry of the DFT, we can restrict attention to half of this matrix, namely the complex values Y k,ℓ , k ∈ {1, . . . , K}, ℓ ∈ {0, . . . , N/2 -1}.

Given a frame k and a bin ℓ, we should write that Y k,ℓ = S k,ℓ + X k,ℓ where, obviously, S k,ℓ and X k,ℓ stand respectively for the speech and noise time-frequency components for the kth frame and the ℓth bin. Since the frames do not overlap, the complex random variables X k,ℓ , k = {1, . . . , K}, ℓ ∈ {0, 1, . . . , N -1}, are iid with X k,ℓ ∼ N c (0, γ 2 ) and γ = σC √ N . Depending on the type of speech signal present during frame k, some speech time-frequency components S k,ℓ can be neglected in comparison with noise and other speech time-frequency components. For instance, high frequency components of voiced speech signals are often negligible in comparison with noise and low-frequency components of the same speech signals; many unvoiced fricative speech signals have low-frequency components significantly smaller than those in high frequency and those due to noise. We model the presence and the absence of the speech time-frequency component S k,ℓ by a discrete random variable ε k,ℓ valued in {0, 1} and write that the observation is Y k,ℓ = ε k,ℓ S k,ℓ + X k,ℓ . With respect to this model, P ({ε k,ℓ = 1}) is the probability that some speech component be present in bin ℓ during the frame k. This probability of presence may be larger than one half for low frequency components; however, for high frequency components, this probability of presence becomes less than or equal to 1/2 and even relatively small.

The In order to deal with m observations that can reasonably be considered as mutually independent, these observations can be chosen randomly amongst the M complex values we have. However, this randomization does not affect significantly the results obtained below.

The Wiener filtering

The T available samples y(t), t = 1, 2, . . . , T , are still split into frames of N = 2 q samples each but, in contrast with the preceding subsection, the frames overlap now by one half and the samples of each frame are weighted. Despite these differences with the foregoing, the notations used above are kept.

The Wiener filtering of the kth frame consists in seeking the complex values W k,ℓ , such that, for every bin ℓ

∈ {0, 1, . . . , N -1}, E[|S k,ℓ -W k,ℓ Y k,ℓ | 2 ] is the least value among all the possible quadratic means E[|S k,ℓ -λY k,ℓ | 2 ]
when λ ranges over the set of complex values. The wellknown solution to this problem is

W k,ℓ = E[|S k,ℓ | 2 ]/E[|Y k,ℓ | 2 ] = E[|S k,ℓ | 2 ] γ 2 + E[|S k,ℓ | 2 ] ( 18
)
since X k,ℓ ∼ N c (0, γ 2 ) and γ = σC √ N . Defining the a priori Signal to Noise Ratio (SNR) by

ρ k,ℓ = k,ℓ | 2 ]/γ 2 , (19)
equation ( 18) can be re-written in the form

W k,ℓ = ρ k,ℓ /(1 + ρ k,ℓ ). ( 20
)
The denoised speech signal in the kth frame is then the inverse DFT of the sequence W k,ℓ , ℓ = 0, 1, . . . , N -1.

The main difficulty in performing an estimate of the a priori SNR is that speech signals are not stationary. According to the standard recursive filtering procedure originally introduced in [3], we estimate W k,ℓ by

W k,ℓ = ρk,ℓ /(1 + ρk,ℓ ) (21)
where

ρk,ℓ = (1 -α)h (ζ k,ℓ -1) + α| W k-1,ℓ Y k-1,ℓ | 2 /γ 2 (22)
can be regarded as an estimate of the a priori SNR ρ k,ℓ . In (22), h(x) = x if x ≥ 0 and h(x) = 0 otherwise, α is some weighting factor such that 0 ≤ α < 1 (we chose α = 0.98 in our experiments commented below), and

ζ k,ℓ = |Y k,ℓ | 2 /γ 2
is the so-called a posteriori SNR. When σ is unknown, the value γ can be estimated by proceeding as described in subsection 5.1. Denoting by γ the estimate returned for γ by the ESE-II, we modify the recursive filtering approach defined by equations ( 21) and (22) as follows. The coefficients W k,ℓ are now estimated by

W k,ℓ = ρ k,ℓ /(1 + ρ k,ℓ ), (23)
where the estimate ρ k,ℓ of the a priori SNR is given by

ρ k,ℓ = (1 -α)h ζ k,ℓ -1 + α| W k-1,ℓ Y k-1,ℓ | 2 / γ 2 , ( 24
) and ζ k,ℓ = |Y k,ℓ | 2 / γ 2 is an estimate of ζ k,ℓ
. The denoised speech signal obtained in frame k is then the inverse DFT of the sequence W k,ℓ , ℓ = 0, 1, . . . , N -1. Since it follows from section 4.3 that the estimate γ should approach significantly well the exact unknown value γ, the performance of the recursive procedure defined by equations (23) and (24) can be expected to be close to that obtained by the filtering approach defined by (21) and (22).

Performance evaluation

We consider twenty sentences of the TIMIT database, downsampled to 8 kHz before adding white Gaussian noise. We estimate the noise standard deviation as described in section 5.1 with frames of N = 256 samples each. A frame corresponds to 32ms of noisy speech signals. For estimating the noise standard deviation, these frames do not overlap and are not weighted. As far as the Wiener filtering is concerned, there is a 50% overlap between two adjacent frames and each frame is weighted by a Hanning window before computing the DFT.

evaluate the quality of the filtered speech signals by means of the standard Segmental Signal to Noise Ratio (SSNR) (see [START_REF] Quackenbush | Objective Measures of Speech Quality[END_REF]) and the Modified Bark Spectral Distortion (MBSD) (see [10]). The SSNR is the average of the SNR values on short segments. The SSNR is not relevant enough to measure the distortion of the denoised speech signals. This is the reason we use the MBSD. The MBSD proves to be highly correlated with subjective speech quality assessment [10].

The average SNNR and MBSD obtained over the twenty sentences randomly chosen within the TIMIT database are presented in figures 3 and 4. The solid curves are the performance measurements achieved with the filtering defined by equations ( 23) and ( 24) where the ESE-II is used to estimate the noise level. The dashed curves are the results obtained when the filtering is achieved along equations (21) and ( 22), that is when the noise standard deviation is known. Clearly, the Wiener filtering adjusted with the noise standard deviation estimate yields results that are significantly close to those obtained when the noise standard deviation is known.

CONCLUSION AND PERSPECTIVE

When signals with unknown probability distributions and priors are additively corrupted by independent WGN, the noise standard deviation can be estimated by the ESE-II. This estimator requires no prior knowledge on the signal probability distributions, does not assume that the signals are iid or that the probabilities of presence of these signals are equal. The observations should be independent; the sample size and the signal amplitudes should be large; however, these conditions are seemingly not so constraining in practice.

A direct application of the ESE-II is the estimation of the noise standard deviation when observations are speech signals in additive and independent WGN. The estimate thus performed serves to adjust a standard Wiener filtering without resorting to any VAD or subspace approaches. The SSNR and the MBSD of the denoised speech signals returned by the resulting filtering are very close to those achieved when the noise standard deviation is known.

A rather natural extension of this work is perceptual filtering. For instance, in [START_REF] Amehraye | On the Application of Recent Results in Statistical Decision and Estimation Theory to Perceptual Filtering of Noisy Speech Signals[END_REF], the same type of estimator is used to adjust some perceptual filtering. Our current work involves the use of the estimator proposed in the present paper to carry out perceptually motivated speech denoising in presence of white and coloured noise.

  ) when m and the amplitudes of the signals are both large. If we now set B m (T ) = 1 m m k=1 I(|Y k | ≤ σh), the same type of intuitive approach suggests that

  is a CWGN with standard deviation σ and ε = (ε k ) k∈N is a sequence of random variables valued in {0, 1} respectively. Assume that (A1) for every k ∈ N, S k , X k and ε k are independent; (A2) the random variables Y k , k ∈ N, are independent; (A3) the random variables ε k , k ∈ N, are independent; (A4) the priors P ({ε k = 1}), k ∈ N, are less than or equal to one half.

  A and m are large enough, proposition 3.1 suggests estimating the noise standard deviation by a possibly local minimum of sup ℓ∈{1,...,L} {D m (x, β ℓ κ(A/x))} (14)

  The search interval used to compute the estimate is|Y [kmin] |, |Y [m] | where Y [k] , k = 1,2, . . . , m stands for the sequence Y k , k = 1, 2, . . . , m sorted by increasing modulus, k min = m/2hm and h = 1/ 4m(1 -Q)

Figure 1 .

 1 Figure 1. BER of the test T σκ(A/σ) the error probability Q(A/σ) for different values of m and A. The signals S k , k ∈ N, are independent, have their probabilities of presence equal to 1/2 and are such that S k = Ae iΦ k where Φ k is uniformly distributed in [0, 2π]. These results suggest the construction of a new estimator, namely the ESE-II, which basically derives from the ESE-I.

Figure 2 .

 2 Figure 2. BER of the test T σκ(A/ σ) versus the error probability Q(A/σ) for different values of m and A. These results were obtained with the same signals and the same experimental protocol as those employed to obtain the results of figure 1.

Figure 3 .

 3 Figure 3. SSNR improvement for speech signals in independent AWGN with various SNRs.

Figure 4 .

 4 Figure 4. MBSD improvement for speech signals in independent AWGN with various SNRs.

Table 1 .

 1 The results are those of table 2. The empirical mean of Ψ m is rather steady when m varies and the empirical standard deviation of this same Empirical mean and empirical MSE of the ESE-I for different values of m

	Sample Size	m = 100 m = 200 m = 500 m = 1000
	Empirical mean 1.2187	1.2289	1.2094	1.2262
	Empirical MSE	0.1275	0.0995	0.0737	0.0756
	Sample Size	m = 100 m = 200 m = 500 m = 1000
	Empirical mean	0.7102	0.7120	0.7069	0.7093
	Empirical standard deviation 0.0255	0.0152	0.0082	0.0064

Table 2 .

 2 Empirical mean and standard deviation of Ψ m for different values of m random variable decreases with the sample size. Even though we only give the results obtained for m = 100, 200, 500 and 1000, the values obtained for other sample sizes less than 1000 are quite the same. The foregoing then suggests defining another estimate σ by setting

Table 3 .

 3 ESE-II is used as follows to estimate γ. We split the observation set {Y k,ℓ } where k ∈ {1, . . . , K}, ℓ ∈ Empirical bias and empirical MSE of the ESE-II(m) for different values of m

	Sample Size	m = 100 m = 200 m = 500 m = 1000
	Empirical bias	1.0029	1.0103	1.0001	1.0041
	Empirical MSE 0.0520	0.0302	0.0159	0.0115
	{0, . . . , N/2 -1}, into subsets of m observations each; each subset is used to perform an estimate of γ via the
	ESE-II; we then compute the average value of the KN/2m
	estimates thus obtained to derive an estimate of γ. Divid-ing this average by C √ N yields an estimate of σ.


Hugues Favreliere
email: hugues.favreliere@univ-savoie.fr

Gaëtan Le Goïc
email: gaetan.le-goic@univ-savoie.fr

Serge Samper
email: serge.samper@univ-savoie.fr

Fabien Formosa
email: fabien.formosa]@univ-savoie.fr

Multi

Multi-scale analysis of surface topographies by modal filtering

come

This work introduces an innova�ve method for the mul�-scale analysis of surface topographies, which consists of applying a method based on a new parameteriza�on. This kind of surface parameteriza�on refers to natural modes of vibra�on, and is therefore named modal parameteriza�on. It allows us to characterize the form, waviness and roughness defects of a surface. This parameteriza�on opens up new fields of analysis, such as the appearance quality of surfaces. It is thereby possible to decompose a measured surface in a vector basis, of which vectors are represented by plane natural eigenmodes sorted by frequency and complexity. Different filtering opera�ons can then be produced, such as extrac�ng the primary form of the surface.


Luca De Feo
email: luca.de-feo@uvsq.fr

Hugues Randriam
email: randriam@enst.fr

Édouard Rousseau
email: erousseau@enst.fr

Standard Lattices of Compatibly Embedded Finite Fields

Keywords: • Mathematics of computing → Mathematical software, • Computing methodologies → Algebraic algorithms, Finite fields, field extensions, Conway polynomials

Lattices of compatibly embedded finite fields are useful in computer algebra systems for managing many extensions of a finite field F p at once. They can also be used to represent the algebraic closure Fp , and to represent all finite fields in a standard manner.

The most well known constructions are Conway polynomials, and the Bosma-Cannon-Steel framework used in Magma. In this work, leveraging the theory of the Lenstra-Allombert isomorphism algorithm, we generalize both at the same time.

Compared to Conway polynomials, our construction defines a much larger set of field extensions from a small pre-computed table; however it is provably as inefficient as Conway polynomials if one wants to represent all field extensions, and thus yields no asymptotic improvement for representing Fp .

Compared to Bosma-Cannon-Steel lattices, it is considerably more efficient both in computation time and storage: all algorithms have at worst quadratic complexity, and storage is linear in the number of represented field extensions and their degrees.

Our implementation written in C/Flint/Julia/Nemo shows that our construction in indeed practical.

INTRODUCTION

Computer algebra systems (CAS) are often faced with the problem of constructing several extensions of a finite field F p in a compatible way, i.e., such that the (subfield) inclusion lattice of the given extensions can be computed and evaluated efficiently.

Concretely, what is sought is a data structure Λ to represent arbitrary collections of extensions of F p , in such a way that elements of F p m are represented in optimal space (i.e., O(m) coefficients), and that arithmetic operations are performed efficiently (i.e., O lcm(l, m) d arithmetic operations to combine an element of F p l and an element of F p m , where d ≤ 3 and, possibly, d = 1 + ε).

To this end, it is useful to set several sub-goals:

Effective embeddings: For any pair of extensions k ⊆ K in Λ, there exists an efficiently computable embedding ϕ : k → K, and algorithms to evaluate ϕ on k, and the section ϕ -1 on K. Compatibility: The embeddings are compatible, i.e., for any triple k ⊆ K ⊆ L in Λ, and embeddings ϕ : k → K, ψ : K → L, χ : k → L, one has χ = ψ • ϕ. Incrementality: The data associated with an extension (e.g., its irreducible polynomial, change-of-basis matrices, . . . ) must be computable efficiently and incrementally, i.e., adding a new field extension to Λ does not require recomputing data for all extensions already in Λ. Uniqueness: Any extension of F p is determined by an irreducible polynomial whose definition only depends on the characteristic p and the degree of the extension. Generality: Extensions of F p can be represented by arbitrary irreducible polynomials.

Some goals, such as incrementality, uniqueness and generality are optional, and it is obvious that uniqueness and generality are even in conflict with each other. An incremental data structure can be used to effectively represent an algebraic closure Fp , with new finite extensions built on the fly as they are needed. Uniqueness is useful for defining field elements in a standard way, portable between different CAS, while generality is useful in a context where the user is left with the freedom of choosing the defining polynomials. Note that any solution can be made unique by replacing all random choices with pseudo-random ones, however one is usually interested in unique solutions that have a simple mathematical description. Also, any solution can be made general by means of an isomorphism algorithm [START_REF] Allombert | Explicit Computation of Isomorphisms between Finite Fields[END_REF][START_REF] Brieulle | Computing isomorphisms and embeddings of finite fields[END_REF][START_REF] Lenstra | Finding isomorphisms between finite fields[END_REF][START_REF] Kumar Narayanan | Fast Computation of Isomorphisms Between Finite Fields Using Elliptic Curves[END_REF]. Other optional goals, such as computing normal bases or evaluating Frobenius morphisms, may be added to the list, however they are out of the scope of this work.

Previous work. The first and most well known solution is the family of Conway polynomials [START_REF] Heath | New algorithms for generating Conway polynomials over finite fields[END_REF][START_REF] Nickel | Endliche Körper in dem gruppentheoretischen Programmsystem GAP[END_REF], first adopted in GAP [START_REF]GAP -Groups, Algorithms, and Programming[END_REF], and then also by Magma [START_REF] Bosma | The MAGMA algebra system I: the user language[END_REF] and Sage [START_REF]The Sage Developers[END_REF]. Conway polynomials yield uniqueness, however computing them requires exponential time using the best known algorithm, thus incrementality is only available at a prohibitive cost; for this reason, they are usually pre-computed and tabulated up to some bound.

Lenstra [START_REF] Lenstra | Finding isomorphisms between finite fields[END_REF] was the first to show the existence of a (incremental, general) data structure computable in deterministic polynomial time. He proved that, besides the problem of finding irreducible polynomials, any other question is amenable to linear algebra. Subsequent work of Lenstra and de Smit [START_REF] Lenstra | Standard models for finite fields[END_REF] tackled the uniqueness problem, albeit only from a theoretical point of view.

In practice, randomized algorithms are good enough for a CAS, then polynomial factorization and basic linear algebra provide an easy (incremental, general) solution, that was first analyzed by Bosma, Cannon and Steel [START_REF] Bosma | Lattices of compatibly embedded finite fields[END_REF], and is currently used in Magma.

All solutions presented so far have superquadratic complexity, i.e., d > 2. Recent work on embedding algorithms [START_REF] De Feo | Fast algorithms for ℓ-adic towers over finite fields[END_REF][START_REF] De Feo | Fast Arithmetic for the Algebraic Closure of Finite Fields[END_REF][START_REF] Doliskani | Computing in degree 2 k -extensions of finite fields of odd characteristic[END_REF] yields subquadratic (more precisely, d ≤ 1.5) solutions for specially constructed (non-unique, non-general) families of irreducible polynomials, and even quasi-optimal ones (i.e., d = 1+ε) if a quasi-linear modular composition algorithm is available. However these constructions involve counting points of random elliptic curves over finite fields, and have thus a rather high polynomial dependency in log p; for this reason, they are usually considered practical only for relatively small characteristic.

Our contribution. In this work we present an incremental, general and/or unique solution for lattices of compatibly embedded finite fields, where all embeddings can be computed and evaluated in quasi-quadratic time. Our starting point is Allombert's [START_REF] Allombert | Explicit Computation of Isomorphisms between Finite Fields[END_REF] and subsequent [START_REF] Brieulle | Computing isomorphisms and embeddings of finite fields[END_REF] improvements to Lenstra's isomorphism algorithm [START_REF] Lenstra | Finding isomorphisms between finite fields[END_REF]. Plugging them in the Bosma-Cannon-Steel framework immediately produces an incremental general solution with quasiquadratic complexity; however we go much further. Indeed, we show that the compatibility requirement can be taken a step further by constructing a lattice of F p -algebras with a distinguished element, which is a byproduct of the Lenstra-Allombert algorithm.

The advantages of our construction over a naive combination of the Lenstra-Allombert algorithm and the Bosma-Cannon-Steel framework are multiple. Storage drops from quadratic to linear in the number of extensions stored in Λ and in their degrees, and the cost of adding a new extension to Λ drops similarly.

Our F p -algebras are constructed by tensoring an arbitrary lattice of extensions of F p with what we call a cyclotomic lattice (see next section). In this work we mostly abstract away from the concrete instantiation of the cyclotomic lattice, only fixing a choice in Section 6, where we use Conway polynomials to analyze the complexity and implement our algorithms. This choice allows us to uniquely represent finite fields of degree exponentially larger than with Conway polynomials alone; however it also has the serious drawback of being generically as hard to compute as Conway polynomials, and thus relatively unpractical. We leave the exploration of other, more practical, instantiations of cyclotomic lattices for future work.

Organization. The presentation is structured as follows. In Section 2 we review some basic algorithms and facts on roots of unity and Conway polynomials. In Section 3 we review the Lenstra-Allombert algorithm and we define and study Kummer algebras, the main ingredient to our construction. In Section 4 we introduce a notion of compatibility for solutions of Hilbert 90 in Kummer algebras, that provides standard defining polynomials for finite fields. Then in Section 5 we again use these compatible solutions to construct standard compatible embeddings between finite fields, from which a lattice can be incrementally constructed. Finally, in Section 6 we give the complexities of our algorithms, and present our implementation.

PRELIMINARIES

Fundamental algorithms and complexities. Throughout this paper we let F p be a finite field. For simplicity, we shall assume that p is prime, which is arguably the most useful case, however our results could easily be extend to non-prime fields. We measure time complexities as a number of arithmetic operations (+, ×, /) over F p , and storage as a number of elements of F p . We let M(m) denote the number of operations required to multiply two polynomials with coefficients in F p of degree at most m, and adopt the usual super-linearity assumptions on the function M (see [START_REF] Von | Modern Computer Algebra[END_REF]Ch. 8.3]).

Any finite extension F p m can be represented as the quotient of F p [X ] by an irreducible polynomial of degree m. The algorithms we present in the next sections need not assume any particular representation for finite fields, however when analyzing their complexities we will assume this representation. Then, multiplications in F p m can be carried out using O(M(m)) operations, and inversions using O(M(m) log(m)). In this work we will also need to perform computations in algebras F p m ⊗ F p n : representing them as quotients of a bivariate polynomial ring, we can multiply elements using O(M(mn)) operations.

We will extensively use a few standard routines, that we recall briefly. Brent and Kung's algorithm [START_REF] Brent | Fast Algorithms for Manipulating Formal Power Series[END_REF] computes the modular composition f (д) mod h of three polynomials f , д, h ∈ F p [X ] of degree at most m using O m (ω+1)/2) operations, where ω is the exponent of linear algebra over F p . The Kedlaya-Umans algorithm [START_REF] Kedlaya | Fast Polynomial Factorization and Modular Composition[END_REF] solves the same problem, and has better complexity in the binary RAM model, however it is widely considered impractical, we shall thus not consider it in our complexity estimates.

By applying transposition techniques [START_REF] Bostan | Tellegen's principle into practice[END_REF][START_REF] Bürgisser | Algebraic Complexity Theory[END_REF] to Brent and Kung's algorithm, Shoup [START_REF] Shoup | Fast construction of irreducible polynomials over finite fields[END_REF][START_REF] Shoup | Efficient computation of minimal polynomials in algebraic extensions of finite fields[END_REF] derived an algorithm to compute minimal polynomials of arbitrary elements of F p m , having the same complexity O m (ω+1)/2) . Kedlaya and Umans' improvements also apply to Shoup's minimal polynomial algorithm, with the same practical limitations.

Shoup's techniques can also be applied to evaluate embeddings of finite fields. Let F p l , F p m be a pair of finite fields related by an embedding ϕ : F p l → F p m , and let a generator α l of F p l and its image ϕ(α l ) in F p m be given. Given these data, for any element x ∈ F p l it is possible to compute its image ϕ(x) using O m (ω+1)/2) operations; similarly, given an element y = ϕ(x) in F p m , it is possible to recover x in the same asymptotic number of operations. The relevant algorithms are summarized in [START_REF] Brieulle | Computing isomorphisms and embeddings of finite fields (extended version)[END_REF]Sec. 6]; note that, for specially constructed generators α l , more efficient algorithms may exist [START_REF] De Feo | Fast algorithms for ℓ-adic towers over finite fields[END_REF][START_REF] De Feo | Fast Arithmetic for the Algebraic Closure of Finite Fields[END_REF][START_REF] De | Fast arithmetics in Artin-Schreier towers over finite fields[END_REF][START_REF] Doliskani | Computing in degree 2 k -extensions of finite fields of odd characteristic[END_REF].

The present work focuses on algorithms to compute embeddings of finite fields, i.e., algorithms that, given finite fields F p l and F p m with l | m, find an embedding ϕ : F p l → F p m , a generator α l of F p l , and its image ϕ(α l ). An extensive review of known algorithms is given in [START_REF] Brieulle | Computing isomorphisms and embeddings of finite fields[END_REF]; here we shall only be interested in the Lenstra-Allombert isomorphism algorithm [START_REF] Allombert | Explicit Computation of Isomorphisms between Finite Fields[END_REF][START_REF] Lenstra | Finding isomorphisms between finite fields[END_REF], and its adaptation to compatible lattices of finite fields.

Conway polynomials and cyclotomic lattices. The algorithms of the next sections will be dependent on the availability of a cyclotomic lattice. By this we mean, formally, a collection

S I = {(K l , ζ l )} l ∈I
over some support set I ⊆ N \ pN. Where K l is an explicitly represented finite extension of F p , and ζ l ∈ K l a generating element that is also a primitive l-th root of unity, so

K l = F p (ζ l ), (ζ l ) l = 1, together with explicit embeddings ι l,m : K l → K m ζ l → (ζ m ) m l whenever l | m.
If I is a finite set of indices, there is an easy randomized algorithm to construct a cyclotomic lattice: compute n = lcm l ∈I (l), construct the smallest field F p a such that n divides p a -1, take a random x (p a -1)/n ∈ F p a and test that it has multiplicative order n; then all roots ζ l in the lattice are constructed as powers of this element, and we can set K l = F p (ζ l ) ⊆ F p a and let ι l,m be natural inclusion.

Nevertheless, the most useful cyclotomic lattices are those where I is the whole set N \ pN. It may seem odd to ask for such data, which in the end provides a representation of Fp , as a prerequisite for a construction whose goal is precisely to represent Fp . However we shall see that a relatively small cyclotomic sub-lattice is enough to construct a much larger lattice of compatibly embedded finite fields; it thus makes sense to assume that a cyclotomic lattice is available, if it can be computed incrementally.

Conway polynomials [START_REF] Nickel | Endliche Körper in dem gruppentheoretischen Programmsystem GAP[END_REF] offer a classic example of cyclotomic lattice. The a-th Conway polynomial C a ∈ F p [X ] is defined as the lexicographically smallest monic irreducible polynomial of degree a that is also primitive (i.e., its roots generate F × p a ) and norm compatible (i.e,

C a X p b -1 p a -1 = 0 mod C b whenever a divides b). The cyclotomic lattice is defined first by letting ζ p a -1 be the image of X in K p a -1 = F p a = F p [X ]/C a for any a; this is then extended to all l ∈ I by setting K l = K p a -1 and ζ l = (ζ p a -1 ) p a -1 l
where F p a is the smallest extension F p containing l-th roots of unity.

The best known algorithm to compute Conway polynomials has exponential complexity [START_REF] Heath | New algorithms for generating Conway polynomials over finite fields[END_REF], hence they are usually precomputed and tabulated up to a certain bound. Most computer algebra systems switch to other ways of representing finite fields when the tables of Conway polynomials are not enough. A notable exception is Sage-Math [START_REF]The Sage Developers[END_REF] (since version 5.13 [START_REF] Roe | Implement pseudo-Conway polynomials[END_REF]), that defines pseudo-Conway polynomials by relaxing the "lexicographically first" requirement; although easier to compute in practice, their computation still requires an exponential amount of work.

Other ways to construct cyclotomic lattices are possible. One may, for example, factor cyclotomic polynomials over F p , being careful to maintain compatibility. In the next sections we shall not suppose any cyclotomic lattice construction in particular, and simply assume that we are given a collection S I that satisfies the properties given at the beginning of this paragraph.

THE LENSTRA-ALLOMBERT ALGORITHM

We now review the theory behind Allombert's adaptation [START_REF] Allombert | Explicit Computation of Isomorphisms between Finite Fields[END_REF] of Lenstra's isomorphism algorithm [START_REF] Lenstra | Finding isomorphisms between finite fields[END_REF]. This will be our stepping stone towards the definition of some standard elements in field extensions of F p with an effective compatibility condition.

The main ingredient of the algorithm is an extension of Kummer theory. Because of this, the algorithm is limited to field extensions F p l of degree l prime to p. The easier case of extensions of degree p e is covered in a similar way using Artin-Schreier theory, and the generic case is solved by separately computing isomorphisms for the power-of-p and the prime-to-p parts, and then tensoring the results together. Due to space constraints, we will not give details for the general case here; see [START_REF] Allombert | Explicit Computation of Isomorphisms between Finite Fields[END_REF][START_REF] Brieulle | Computing isomorphisms and embeddings of finite fields[END_REF][START_REF] Lenstra | Finding isomorphisms between finite fields[END_REF].

For any finite extension of F p , we denote by σ : x → x p the Frobenius automorphism. Let l be an integer not divisible by p. Then σ is an F p -linear endomorphism of F p l with minimal polynomial T l -1, separable but not necessarily split, i.e., F p l is not necessarily a Kummer extension of F p . We extend scalars and work in the Kummer algebra of degree l:

A l = F p l ⊗ F p (ζ l ),
where ⊗ is the tensor product over F p , and ζ l is a primitive l-th root of unity, taken from the given cyclotomic lattice. We call F p (ζ l ) the field of scalars of A l , and we define the level of A l as

ν(l) = ord (Z/l Z) × (p) = [F p (ζ l ) : F p ],
that is, the degree of its field of scalars. Now σ ⊗ 1 is a 1 ⊗ F p (ζ l )-linear endomorphism of A l with l distinct eigenvalues, namely the powers of 1 ⊗ ζ l . Thus, if η = ζ i l is any l-th root of unity in F p (ζ l ), the corresponding eigenspace is defined by the Hilbert 90 equation for η:

(σ ⊗ 1)(x) = (1 ⊗ η)x, (H90)
which plays the role of σ (x) = ηx in classical Kummer theory. The solutions of (H90) in A l form a 1 ⊗ F p (ζ l )-vector space of dimension 1, and if x is such a solution for η, then x j is a solution for η j . In particular, let α l be a nonzero solution of (H90) for ζ l . Then 1, α l , • • • , (α l ) l -1 are eigenvectors for distinct eigenvalues and thus form a basis of

A l over 1 ⊗ F p (ζ l ). Likewise (α l ) l = 1 ⊗ c l
for some scalar c l ∈ F p (ζ l ), that we shall call the Kummer constant of α l . This proves: Proposition 1. Any nonzero solution α l of (H90) for ζ l is a generating element for A l as an algebra over 1 ⊗ F p (ζ l ), inducing an isomorphism

A l ≃ F p (ζ l )[T ]/(T l -c l ).
Since A l is known to be an étale algebra, α l being nonzero implies that c l is nonzero, which in turn implies that α l is invertible in A l , indeed (α l ) -1 = (1 ⊗ c -1 l )(α l ) l -1 . We will make frequent use of the following: Lemma 2. Let K, L be two finite extensions of F p . Then, for any

β ∈ K ⊗ L, we have (σ ⊗ σ )(β) = β p . Proof. If β = u ⊗ v is an elementary tensor we have (σ ⊗ σ )(β) = u p ⊗ v p = β p
. This then extends by linearity since we're in characteristic p. □

In this generality we also introduce the following notation: if η ∈ L has degree d over F p , then any β ∈ K ⊗ F p (η) ⊆ K ⊗ L decomposes uniquely as β = d -1 i=0 y i ⊗ η i , and we set ⌊β⌋ η = y 0 .

In particular, coming back to A l , if we write

α l = a-1 i=0 x i ⊗ ζ i l where a = ν (l), it is shown in [1] that x 0 = ⌊α l ⌋ ζ l
is a generating element for the extension F p l /F p . Moreover, [START_REF] Allombert | Explicit Computation of Isomorphisms between Finite Fields[END_REF] (see also [START_REF] Brieulle | Computing isomorphisms and embeddings of finite fields[END_REF]) provides the following equations that allow, in the opposite direction, to recover α l from x 0 :

x a-1 = σ (x 0 )/b 0 x i = σ (x i+1 ) -b i+1 x a-1 for i = a-2 down to 1 (1)
where

ζ a l = a-1 i=0 b i ζ i l is the minimal equation for ζ l .
Proposition 3. With the notations above, there are precisely l elements x ∈ A l that are solutions of (H90) for ζ l and satisfy

x l = 1 ⊗ c l , namely, these are the (1 ⊗ ζ l ) u α l = (σ u ⊗ 1)(α l ) for 0 ≤ u < l. The corresponding generating elements for F p l /F p are the ⌊(σ u ⊗ 1)(α l )⌋ ζ l = σ u (x 0 )
; they all have the same minimal polynomial, which is a generating polynomial for F p l /F p depending only on c l .

Proof. The solutions of (H90) for ζ l form a 1 ⊗ F p (ζ l )-vector space of dimension 1, thus they all are of the form x = (1 ⊗ ξ )α l . Adding the condition x l = 1 ⊗ c l then forces ξ l = 1, from which all assertions follow. □

Now we consider Kummer algebras of various degrees. Since we assumed that the fields of scalars are defined from a cyclotomic lattice S I , they are compatibly embedded: for l | m prime to p, we have the embedding ι l,m :

F p (ζ l ) → F p (ζ m ) ζ l → (ζ m ) m l .
It is easily shown that, as an F p -algebra, A l is isomorphic to a product of copies of F p l (ζ l ), and A m to a product of copies of F p m (ζ m ). This allows us to describe all F p -algebra morphisms from A l to A m . However here we will focus only on a certain subclass of them: Definition 4. A Kummer embedding of A l into A m is an injective F p -algebra morphism Φ : A l → A m such that:

• Φ extends the scalar embedding 1 ⊗ ι l,m • Φ commutes with σ ⊗ 1.

Proposition 5. Let α l ∈ A l be a nonzero solution of (H90) for ζ l , with Kummer constant c l . Then, there is a 1-to-1 correspondence between Kummer embeddings Φ : A l → A m and solutions α ∈ A m of (H90) for (ζ m ) m l that satisfy ( α) l = 1 ⊗ ι l,m (c l ), given by

Φ ←→ α = Φ(α l ).
Proof. Direct consequence of Proposition 1 and Definition 4. □

Actually, Kummer embeddings are easily characterized: Proposition 6. There is a natural 1-to-1 correspondence between Kummer embeddings Φ : A l → A m and embeddings of finite fields ϕ : F p l → F p m , given by:

• If Φ is a Kummer embedding, then Φ maps F p l ⊗1 into F p m ⊗1.
Thus the restriction of Φ to F p l ⊗ 1 is of the form ϕ ⊗ 1 for some ϕ :

F p l → F p m . • Conversely, if ϕ : F p l → F p m is
an embedding of finite fields, then Φ = ϕ ⊗ ι l,m is a Kummer embedding. Moreover, this correspondence commutes with composition of embeddings.

Proof. Let Φ be a Kummer embedding. Being a F p -algebra morphism, it satisfies Φ(β p ) = Φ(β) p for all β ∈ A l . By Lemma 2, this means that Φ commutes with σ ⊗ σ , and thus, also with

(σ ⊗ 1) -1 • (σ ⊗ σ ) = 1 ⊗ σ . This implies that Φ maps F p l ⊗ 1 into F p m ⊗ 1.
The other assertions are clear. □ Corollary 7. Let α l ∈ A l be a nonzero solution of (H90) for ζ l , with Kummer constant c l , and let α ∈ A m be a solution of (H90) for

(ζ m ) m l that satisfies ( α) l = 1 ⊗ ι l,m (c l ). Then: • α ∈ F p m ⊗ F p ((ζ m ) m l ) ⊆ A m ; • the assignation ⌊α l ⌋ ζ l → ⌊ α⌋ (ζ m ) m l defines an embedding ϕ : F p l → F p m ; • Φ = ϕ ⊗ ι l,m is the unique Kummer embedding such that Φ(α l ) = α.
Proof. By Proposition 5 there is a unique Kummer embedding Φ such that Φ(α l ) = α. By Proposition 6 we have that Φ = ϕ ⊗ ι l,m for some ϕ :

F p l → F p m . Writing α l = a-1 i=0 x i ⊗ ζ i l , it follows that α = Φ(α l ) = a-1 i=0 ϕ(x i ) ⊗ (ζ m ) mi l . Thus ⌊ α⌋ (ζ m ) m l = ϕ(x 0 ) = ϕ(⌊α l ⌋ ζ l )
, and, since ⌊α l ⌋ ζ l generates F p l , this uniquely characterizes ϕ. □

We can now state Allombert's algorithm and prove its correctness; we give below a minor variation on the original algorithm, better adapted to our more general setting.

STANDARD SOLUTIONS OF (H90)

Plugging Algorithm 1 into the Bosma-Cannon-Steel framework immediately gives a way to compatibly embed arbitrary finite fields. However, there are two points in Allombert's algorithm on which we would like to improve:

Uniqueness: As mentioned, the element ⌊α l ⌋ ζ l is a generating element for F p l , or equivalently, it provides a defining irreducible polynomial of degree l. However this polynomial depends on the choice of α l (even though only through c l , cf. Proposition 3). Compatibility: The embedding ϕ depends on a constant κ, which itself depends on the choice of α l , α m (and of a l-th root extraction). Thus, given a certain number of finite fields, in order to ensure compatibility of the various embeddings between them, one has to keep track of these constants κ for all pairs (l, m), which grow quadratically with the number of fields. It would be useful if one could force κ = 1, that is, if α l , α m and Φ : A l → A m could be chosen so that Φ(α l ) = (α m ) m l . From the description of Algorithm 1, this requires c m = ι l,m (c l ). Thus, necessarily c m lies in the subfield F p ((ζ m ) m l ) of F p (ζ m ). Possibly this condition could fail if A l and A m do not have the same field of scalars. This motivates: Definition 9. A Kummer algebra is complete if it is of the largest degree for a given level.

Thus, the complete Kummer algebra of level a is

A p a -1 = F p p a -1 ⊗ F p a with field of scalars F p a = F p (ζ p a -1 )
given by the corresponding ζ p a -1 in our cyclotomic lattice S I , e.g., defined by a (pseudo)-Conway polynomial of degree a.

Lemma 10. All nonzero solutions α p a -1 ∈ A p a -1 of (H90) for ζ p a -1 have the same Kummer constant c p a -1 = (ζ p a -1 ) a .

Proof. From Lemma 2 and the fact that σ a is trivial on F p a ≃ F p (ζ p a -1 ) we get that

(α p a -1 ) p a = (σ a ⊗ σ a )(α p a -1 ) = (σ a ⊗ 1)(α p a -1 ) = (1 ⊗ ζ p a -1 ) a α p a -1 .
We conclude since α p a -1 is invertible. □ Definition 11. Let l be an integer prime to p. We define the standard Kummer constant of order l as

c std l = (ι l,p a -1 ) -1 ((ζ p a -1 ) a ) ∈ F p (ζ l ) where a = ν (l) is the level of A l .
We say a solution α l ∈ A l of (H90) for ζ l is standard if its Kummer constant is standard:

(α l ) l = 1 ⊗ c std l . Then, by a decorated Kummer algebra we mean a pair

(A l , α l )
with such α l standard.

Observe that ι l,p a -1 is an isomorphism when a = ν (l), so c std l is well defined.

For complete algebras, Lemma 10 asserts that all nonzero α p a -1 are standard. Proposition 12. Let l be an integer not divisible by p. Then A l can be decorated, i.e., it admits a standard α l . Moreover, this α l is unique up to a l-th root of unity.

Proof. Let α ′ l be any nonzero solution of (H90) for ζ l . Set a = ν (l), pick any α p a -1 ∈ A p a -1 standard (Lemma 10), and pick any Kummer embedding Φ : A l → A p a -1 (Proposition 6). Then Φ(α ′ l ) and (α p a -1 )

p a -1 l are two nonzero solutions of (H90) for

ι l,p a -1 (ζ l ) = (ζ p a -1 ) p a -1 l in A p a -1 , thus there is a scalar λ ∈ F p (ζ p a -1 ) such that (α p a -1 ) p a -1 l = (1 ⊗ λ)Φ(α ′ l ) = Φ((1 ⊗ λ)α ′ l ),
where λ = (ι l,p a -1 ) -1 (λ) ∈ F p (ζ l ).

Setting α l = (1 ⊗ η λ)α ′ l ∈ A l for η ∈ F p (ζ l ), we get c l = η l (ι l,p a -1 ) -1 (c std p a -1 ) = η l c std l , and thus the standard

α l ∈ A l are the (1 ⊗ ζ u l λ)α ′ l , for 0 ≤ u < l. □ Definition 13. A generating element s ∈ F p l is called standard if it is of the form s = ⌊α l ⌋ ζ l for α l ∈ A l a standard solution of (H90).
The standard defining polynomial P l for F p l is then the minimal polynomial over F p of such a standard s.

By Proposition 3, we note that P l is entirely determined by c std l , and thus, by the given cyclotomic lattice S I , possibly up to order p ν (l ) -1. As an example, we give in Table 1 the first ten standard polynomials induced by the system of Conway polynomials for p = 2 (thus, in this example, P l only depends on the Conway polynomial of degree ν(l)).

We remark that it is easy to extend the definitions of decorated algebras and standard elements to any extension degree, similarly to the way this is done for the basic Lenstra-Allombert algorithm. Use any (standard) construction for Artin-Schreier towers over finite fields (e.g., [START_REF] De | Fast arithmetics in Artin-Schreier towers over finite fields[END_REF]), define decorated algebras by tensoring together Kummer algebras and Artin-Schreier extensions of F p , and define standard elements as, e.g., the product of a solution of multiplicative H90 and one of additive H90. While this solution is simple and effective, it is rather orthogonal to our work, hence we omit the details here.

The decoration of an algebra A l , and the associated standard generating element and polynomial for F p l , can be computed by the simple adaptation of Allombert's algorithm presented below.

Algorithm 2 (Decoration -Standardization)

Input: F p l , for l prime to p, and S I a cyclotomic lattice. Output: (A l , α l ) decorated, P l standard irreducible polynomial of degree l, and s ∈ F p l standard generating element inducing

F p l ≃ F p [T ]/(P l ). 1: Prepare the Kummer algebra A l . 2: Prepare c std l = (ι l,p a -1 ) -1 ((ζ p a -1 ) a ) ∈ F p (ζ l ). 3: Find α ′ l ∈ A l nonzero solution of (H90) for ζ l . 4: Compute its Kummer constant: (α ′ l ) l = 1 ⊗ c ′ l . 5: Compute κ a l-th root of c std l /c ′ l . 6: Set α l = (1 ⊗ κ)α ′ l . 7: Compute P l the minimal polynomial of ⌊α l ⌋ ζ l over F p . 8: Return (A l , α l ), P l , and ⌊α l ⌋ ζ l .
Algorithm 2 is correct, indeed Proposition 12 ensures that a standard α l exists, and thus α ′ l = (1 ⊗ κ -1 )α l for some κ ∈ F p (ζ l ), so c std l /c ′ l = κ l is a l-th power. By design, decorated Kummer algebras of the same level admit standard Kummer embeddings, under which the corresponding standard solutions of (H90) are power-compatible: Proposition 14. Let l | m be integers prime to p and such that ν (l) = ν (m) = a. Let (A l , α l ), (A m , α m ) be decorated Kummer algebras of degree l, m respectively (and of the same level a). Then, there is a unique Kummer embedding

Φ std l,m : A l → A m such that Φ std l,m (α l ) = (α m ) m l .
Proof. Proposition 5 with α = (α m ) m l . □

Most often we will apply Proposition 14 with m = p a -1.

On the other hand, since power-compatibility implies c m = ι l,m (c l ), it cannot be satisfied for a Kummer embedding between decorated Kummer algebras of different levels. However, at least between complete decorated algebras, we can request some normcompatibility instead.

x + 1

x 9 + x 7 + x 4 + x 2 + 1 x 3 + x + 1

x 11 + x 8 + x 7 + x 6 + x 2 + x + 1 x 5 + x 3 + 1

x 13 + x 10 + x 5 + x 3 + 1 x 7 + x + 1

x 15 + x + 1

x 17 + x 11 + x 10 + x 8 + x 7 + x 6 + x 4 + x 3 + x 2 + x + 1 x 19 + x 17 + x 16 + x 15 + x 14 + x 13 + x 12 + x 8 + x 7 + x 6 + x 5 + x 3 + 1 Let A m be a Kummer algebra of level b = ν (m), so

A m = F p m ⊗ F p (ζ m ) ∼ -→ F p m ⊗ F p (ζ p b -1 ),
where the isomorphism is given by 1 ⊗ι m,p b -1 . Then, for an integer a | b, the subalgebra of A m invariant under 1 ⊗ σ a is identified by this isomorphism with

(A m ) 1⊗σ a ≃ F p m ⊗ F p ((ζ p b -1 ) p b -1 p a -1 ),
where

(ζ p b -1 ) p b -1 p a -1 = N F p b /F p a (ζ p b -1 ) = ι p a -1,p b -1 (ζ p a -
N b/a,A n : (A n ) 1⊗σ b → (A n ) 1⊗σ a γ → 0≤j < b a (1 ⊗ σ ja )(γ ).
This is well-defined, i.e., the image of N b/a,A n is invariant under 1 ⊗ σ a as specified. Often the ambient algebra A n will be implicit, and we will write N b/a instead of N b/a,A n .

By construction, N b/a acts on 1 ⊗ F × p b as 1 ⊗ N F p b /F p a . Scalar norms are multiplicative:

N b/a (γγ ′ ) = N b/a (γ )N b/a (γ ′ ), transitive: N c/a = N b/a • N c/b ,
and they commute with σ ⊗ 1.

Proposition 16. Let a | b be integers, and let (A p a -1 , α p a -1 ), (A p b -1 , α p b -1 ) be decorated complete Kummer algebras of level a, b respectively. Then there is a unique Kummer embedding

Φ std p a -1,p b -1 : A p a -1 → A p b -1 such that Φ std p a -1,p b -1 (α p a -1 ) = N b/a (α p b -1 ).
Proof. By Lemma 2 and the properties of the norm,

(N b/a (α p b -1 )) p a = (σ a ⊗ σ a )(N b/a (α p b -1 )) = (σ a ⊗ 1)(N b/a (α p b -1 )) = N b/a ((σ a ⊗ 1)(α p b -1 )) = N b/a ((1 ⊗ (ζ p b -1 ) a )α p b -1 ) = (1 ⊗ ι p a -1,p b -1 (ζ p a -1 ) a )N b/a (α p b -1 ). So α = N b/a (α p b -1 ) satisfies ( α) p a -1 = 1 ⊗ ι p a -1,p b -1 (c p a -1
) and we conclude with Proposition 5. □

STANDARD EMBEDDINGS

In Proposition 14, we saw how to construct a standard powercompatible embedding of a decorated Kummer algebra into its decorated complete algebra, and in Proposition 16, a standard normcompatible embedding between decorated complete algebras of dividing levels.

(A p a -1 , α p a -1 ) Φ std p a -1,p b -1 ---------→ (A p b -1 , α p b -1 ) Φ std l,p a -1     Φ std m,p b -1 (A l , α l ) (A m , α m )
of standard embeddings of decorated algebras.

Lemma 17. In this setting, there exists a unique Kummer embedding

Φ std l,m : A l → A m that makes the diagram commute. Proof. Consider α = Φ std p a -1,p b -1 (Φ std l,p a -1 (α l )) ∈ A p b -1 .
Then α is invariant under σ l ⊗ 1 and under 1 ⊗ σ a (because α l is), thus, a fortiori, invariant under σ m ⊗ 1 and under 1 ⊗ σ b , which means it lies in the image of Φ std m,p b -1 . We can then set α = (Φ std m,p b -1

) -1 ( α).

If Φ std l,m exists, then it necessarily maps α l to α. However, chasing in the diagram, it is easily seen that α is a solution of (H90) for (ζ m ) m l that satisfies ( α) l = 1 ⊗ ι l,m (c l ), and we conclude with Proposition 5. □

This existence result is "constructive", but impractical, since it requires computations in the possibly very large algebra A p b -1 . However, as in Algorithm 1, one should be able to write α = (1 ⊗ κ)(α m ) m l for some κ ∈ F p (ζ m ). Moreover α is uniquely determined by our data, thus, so should κ. Now our aim is to give an explicit expression for this κ = κ l,m . We start with the case of complete algebras.

Proposition 18. In the complete algebra A p b -1 we have

(α p b -1 ) p b -1 p a -1 = (1 ⊗ ζ p b -1 ) (b-a)p b+a -bp b +ap a (p a -1) 2 N b/a (α p b -1 ).
Proof. Using first Lemma 2, and then (H90), we get:

(α p b -1 ) p b -1 p a -1 N b/a (α p b -1 ) = 0≤j < b a (σ ja ⊗ σ ja )(α p b -1 ) (1 ⊗ σ ja )(α p b -1 ) = 0≤j < b a (1 ⊗ σ ja ) (σ ja ⊗ 1)(α p b -1 ) α p b -1 = 0≤j < b a (1 ⊗ σ ja )(1 ⊗ ζ p b -1 ) ja = (1 ⊗ ζ p b -1 ) 0≤j < b a jap j a .
We conclude thanks to the identity 0≤j <n ).

jT j = T d dT T n -1 T -1 = (n -1)T n+1 -nT n + T (T -1
Proof. It suffices to check that Φ std Proof. We have to show that Φ std l,n (α l ) = Φ std m,n (Φ std l,m (α l )); nothing but a pleasant computation with the explicit constants given by Corollary 19.

Alternatively, set a = ν (l), b = ν (m), c = ν (n), and decorate

A p a -1 , A p b -1 , A p c -1 . It suffices to show that the elements Φ std l,n (α l ) and Φ std m,n (Φ std l,m (α l )) have the same image under Φ std n,p c -1 in A p c -1 . Chasing in the diagram A p a -1 ------→ A p b -1 ------→ A p c -1       A l ------→ A m ------→ A n we see that this common image is N c/a (α p c -1 ) p a -1 l .
□ By a decorated finite field (of degree l, an integer prime to p, and relative to a given cyclotomic lattice S I ), we mean a pair (F p l , s l ), where F p l is a finite field, and s l ∈ F p l a standard generating element in the sense of Definition 13.

We can finally state:

Algorithm 3 (Standard compatible embeddings)

Input: S I a cyclotomic lattice, and (F p l , s l ), (F p m , s m ), decorated finite fields, for l | m integers prime to p. Output: t ∈ F p m , such that the assignation s l → t defines a standard embedding ϕ std l,m : F p l → F p m , compatible with composition. 1: Prepare the Kummer algebras A l and A m . 2: Recover α l from s l and α m from s m using equations (1). Proof. Corollary 7 and Proposition 20. □

IMPLEMENTATION

In the previous sections we kept the description of Kummer embeddings abstract, leaving many computational details unspecified. There are various ways in which our algorithms can possibly be implemented, depending on how one chooses to represent finite fields and the cyclotomic lattice S I . A reasonable option is to use (pseudo)-Conway polynomials to represent the fields F p (ζ p a -1 ), and deduce from them the smallest possible representation for any other field F p (ζ l ). Assuming this technique, we can prove a bound on the complexity of our algorithms.

Proposition 22. Given a collection of (pseudo)-Conway polynomials for F p , of degree up to d, standard solutions α l of (H90) can be computed for any l | (p i -1) for any i ≤ d using O M(l 2 ) log(l) + M(l) log(l) log(p) operations. After that, Kummer embeddings F p l ⊆ F p m can be computed using O(M(m 2 ) log(m)) operations.

Proof. Let a = ν (l) be the level of A l . We take the a-th polynomial from the collection of (pseudo)-Conway polynomials, and use it to define ζ p a -1 . Because a ∈ O(l), the cost of multiplications in F p (ζ p a -1 ) will be bounded by O(M(l)).

From ζ p a -1 , we compute ζ l using O(lM(l)) operations, and its minimal polynomial in O(l (ω+1)/2 ). Then, the Kummer constant c std l = (ζ p a -1 ) a is computed in negligible time, and its expression in the power basis of ζ l is computed in O(l (ω+1)/2 ) using the algorithms for evaluating embeddings mentioned in Section 2.

To construct the Kummer algebra A l = F p l ⊗ F p (ζ l ) we need an irreducible polynomial of degree l, not necessarily related to the (pseudo)-Conway polynomials used to represent the fields of scalars. Very efficient, quasi-optimal algorithms for finding such a polynomial are given in [START_REF] Bostan | Fast computation of special resultants[END_REF][START_REF] Couveignes | Fast construction of irreducible polynomials over finite fields[END_REF][START_REF] De Feo | Fast algorithms for ℓ-adic towers over finite fields[END_REF], we can thus neglect this cost.

The cost of computing a solution α ′ l to (H90) was extensively studied in [START_REF] Brieulle | Computing isomorphisms and embeddings of finite fields[END_REF], where it was found to be bounded by O M(l 2 ) log(l)+ M(l) log(p) . Then, the constant c ′ l = (α ′ l ) l is computed using O M(l 2 ) log(l) operations, and the l-th root κ is computed in O(M(l) log(l) log(p)) according to [START_REF] Brieulle | Computing isomorphisms and embeddings of finite fields[END_REF].

α l = (1 ⊗ κ)α ′
l is then computed in a negligible number of operations.

Finally, the projection ⌊α l ⌋ ζ l comes for free, and its minimal polynomial P l is computed again in O(l (ω+1)/2 ) operations. Now, in order to compute a Kummer embedding of (F p l , α l ) into (F p m , α m ), we compute the scalar κ l,m in O(mM(m)) operations, and (α m ) We thus need to evaluate x → ⌊Tr(xη)⌋ ζ l for many values in F p (ζ m ), but this is a F p -linear form, hence we can precompute its vector on the power basis of ζ m . Let h m , h l be the minimal polynomials of ζ m , ζ l , and let b, a be their degrees. Let h 0 be the constant coefficient of h l , and let

τ = - h 0 (ζ m ) m l h ′ m (ζ m ) h ′ l ((ζ m ) m l ) ∈ F p (ζ m ), direct calculation shows that b-1 i=0 Tr(ζ i m ) ζ l Z i = τ (Z -1 ) Zh m (Z -1 ) mod Z b ,
where by τ (Z ) we mean τ ∈ F p (ζ m ) seen as a polynomial in ζ m . Hence, we can compute the vector of the linear form x → ⌊Tr(x)⌋ ζ l using only basic polynomial arithmetic and modular composition, i.e., in O(m (ω+1)/2 ) operations. Finally, we compute (1 ⊗ κη)(α m ) m l , we see it as a polynomial with coefficients in F p (ζ m ), and we apply the map ⌊Tr(x)⌋ ζ l to each coefficient. This costs O(mM(m)) operations. □

We remark that storing the decorated fields (F p l , α l ) requires O(l 2 ) field elements, however, using the formulas in [START_REF] Allombert | Explicit Computation of Isomorphisms between Finite Fields[END_REF][START_REF] Brieulle | Computing isomorphisms and embeddings of finite fields[END_REF], it is possible to only store ⌊α l ⌋ ζ l , and recover all other coefficients of α l in O(lM(l) log(p)) operations.

To demonstrate the feasibility of this approach, we implemented it in the Julia-based CAS Nemo [START_REF] Fieker | Nemo/Hecke: Computer Algebra and Number Theory Packages for the Julia Programming Language[END_REF], with performance critical routines written in C/Flint [START_REF] Hart | FLINT: Fast Library for Number Theory[END_REF]. Our code is available as a Julia package at https://github.com/erou/LatticeGFH90.jl.

We tested Algorithms 2 and 3 for various small primes, using precomputed Conway polynomials available in Nemo. We do not see major differences between different primes. In Figure 1 we report timings obtained for the case p = 3, on an Intel Core i7-7500U CPU clocked at 2.70GHz, using Nemo 0.11.1 running on Julia 1.1.0, and Nemo's current version of Flint. The plot on the left shows timings for Algorithm 2, for degrees l growing from 1 to 200, for every l coprime to p and such that the ν (l)-th Conway polynomial is available in Nemo; the color scale shows the level of the associated algebra A l . The bottleneck of this algorithm appears to be the l-th root extraction routine.

The plot on the right shows timings for Algorithm 3, measured by computing the standard embedding of F p 2 in F p l . As expected, computing the embeddings takes negligible time in comparison to the decoration of the finite fields. We also tested embedding fields larger than F p 2 , and noticed that the running time mostly depends on the size of the larger field.

CONCLUSION AND FUTURE WORK

We presented a new family of standard compatible polynomials for defining finite fields. Its construction being dependent on the availability of Conway polynomials, it has, at the present moment, very little practical impact; its existence is nevertheless remarkable in itself. It is even evident that computing our standard polynomials is essentially equivalent to computing Conway polynomials; indeed from α l one can immediately deduce (ζ p a -1 ) a , and by taking an a-th root (doable in polynomial time in l), deduce ζ p a -1 and the associated Conway polynomial. Hence, an efficient algorithm for computing our polynomials (for arbitrary degrees) would imply an efficient algorithm to compute Conway polynomials, which would be unexpected.

However, our proposed implementation is not the only possible way to exploit our definitions. It would be interesting, indeed, to find some middle ground between the flexibility of the Bosma-Steel-Cannon framework and the rigidity of Conway polynomials, for example by lazily enforcing the conditions required to have a standard solution of (H90), while incrementally constructing the lattice of roots of unity.

Another line of work would be to give a complete implementation of a lattice of finite fields, not limited to extensions of degree coprime to p. We leave these questions for future work.

Proposition 8 . 5 :

 85 Algorithm 1 is correct: it returns elements that define an embedding ϕ : F p l → F p m . Proof. By Propositions 5 and 6, there exists α ∈ A m solution of (H90) for(ζ m ) m l that satisfies ( α) l = 1 ⊗ ι l,m (c l ). On the other hand, (α m ) m l ∈ A m is also a solution of (H90) for (ζ m ) m l , thus α = (1 ⊗ λ)(α m ) m l for some λ ∈ F p (ζ m ). It follows that ι l,m (c l )/c m = λ lis a l-th power, and κ = (ζ m ) um l λ for some integer u. Now we can replace α with (1 ⊗ (ζ m ) um l )( α) = (1 ⊗ κ)(α m ) m l and conclude with Corollary 7. □ From this proof and Proposition 3, it follows that another choice of the l-th root κ only changes ϕ by a power of σ . Algorithm 1 (Allombert's algorithm) Input: F p l , F p m , for l | m integers prime to p, and a cyclotomic lattice S {l,m } . Output: s ∈ F p l , t ∈ F p m , such that the assignation s → t defines an embedding ϕ : F p l → F p m . 1: Prepare the Kummer algebras A l and A m . 2: Find α l ∈ A l and α m ∈ A m , nonzero solutions of (H90) for ζ l and ζ m respectively. 3: Compute their Kummer constants: (α l ) l = 1 ⊗ c l and (α m ) m = 1 ⊗ c m . 4: Compute κ, a l-th root of ι l,m (c l )/c m . Return ⌊α l ⌋ ζ l and (1 ⊗ κ)(α m )

Definition 15 .

 15 1 ), and with N F p b /F p a the norm of the field extension F p b /F p a . Given a Kummer algebra A n , and some integers a | b | ν(n), we define the scalar norm operator

  Now, consider general l | m not divisible by p, set a = ν(l), b = ν (m), and consider the diagram

) 2 .□Corollary 19 .

 219 Let (A l , α l ) and (A m , α m ) be decorated Kummer algebras, of respective degrees l | m prime to p. Then the standard Kummer embedding Φ std l,m : A l → A m is defined by the assignation α l → (1 ⊗ κ l,m )(α m ) m l , where κ l,m = (ι m,p b -1 ) -1 ((ζ p b -1 ) -(b-a)p b+a -bp b +ap a (p a -1)l

p a - 1 ,p b - 1 (Proposition 20 .

 1120 Φ std l,p a -1 (α l )) and the image of the right-hand-side under Φ std m,p b -1 coincide in A p b -1 . However, we have Φ std p a -1,p b -1 (Φ std l,p a -1 (α l )) = N b/a (α p b -1 ) Standard Kummer embeddings are compatible with composition: if (A l , α l ), (A m , α m ), and (A n , α n ) are decorated Kummer algebras with l | m | n, the corresponding standard embeddings satisfy Φ std l,n = Φ std m,n • Φ std l,m .

3 :

 3 Compute κ l,m = (ι m,p b -1 ) -1 ((ζ p b -1 ) -(b-a)p b+a -bp b +ap a (p a -1)l ) where a = ν (l), b = ν(m).

4 :.Proposition 21 .

 421 Return (1 ⊗ κ)(α m ) m l (ζ m ) m l Standard finite field embeddings are compatible with composition: if (F p l , α l ), (F p m , α m ), and (F p n , α n ) are decorated finite fields with l | m | n, the corresponding standard embeddings satisfy ϕ std l,n = ϕ std m,n • ϕ std l,m .

ml

  in O(M(m 2 ) log(m)). We then need to convert (α m ) m l in the power basis of ζ l . Applying a generic change of basis algorithm as before would be too expensive: indeed we have to convert m coefficients from the field of scalars F p (ζ m ) to F p (ζ l ), which would cost O(m (ω+3)/2 ). Instead we notice that we are only interested in the value (1 ⊗ κ)(α m ) m l ζ l , therefore we proceed as follows. Let Tr denote the trace map from F p (ζ m ) to F p (ζ l ), and let η ∈ F p (ζ m ) be such that Tr(η) = 1. Then the map x → Tr(xη) sends (ζ m ) m l to ζ l , and is F p (ζ l )-linear, it thus agrees with the inverse map of ζ l → (ζ m ) m l on the image of F p (ζ l ).

Figure 1 :

 1 Figure 1: Timings for computing decorated fields (F p l , α l ) (left, log scale), and for computing the standard Kummer embedding from F p 2 to F p l (right) for p = 3.

Table 1 :

 1 The first ten standard polynomials derived from Conway polynomials for p = 2.

ACKNOWLEDGMENTS

We would like to thank the anonymous reviewers for their useful comments. We thank Éric Schost for fruitful discussions and for helping bootstrap this work during a visit by two of the authors to the University of Waterloo. We acknowledge financial support from the French ANR-15-CE39-0013 project Manta, the OpenDreamKit Horizon 2020 European Research Infrastructures project (#676541), and from the French Domaine d'Intérêt Majeur Math'Innov.