Sarah De Nigris

Timoteo Carletti

Renaud Lambiotte

Onset of anomalous diffusion from local motion rules

Keywords: numbers: 05.40.Fb, 02.50.-r, 05.60.Cd

Anomalous diffusion processes, in particular superdiffusive ones, are known to be efficient strategies for searching and navigation by animals and also in human mobility. One way to create such regimes are Lévy flights, where the walkers are allowed to perform jumps, the "flights", that can eventually be very long as their length distribution is asymptotically power-law distributed. In our work, we present a model in which walkers are allowed to perform, on a 1D lattice, "cascades" of n unitary steps instead of one jump of a randomly generated length, as in the Lévy case, where n is drawn from a cascade distribution pn. We show that this local mechanism may give rise to superdiffusion or normal diffusion when pn is distributed as a power law. We also introduce waiting times that are power-law distributed as well and therefore the probability distribution scaling is steered by the two PDF's power-law exponents. As a perspective, our approach may engender a possible generalization of anomalous diffusion in context where distances are difficult to define, as in the case of complex networks, and also provide an interesting model for diffusion in temporal networks.

I. INTRODUCTION

Diffusion processes, when seen as the continuous limit of a random walk, are well known to display uncanny properties when the associated probability distribution of length or duration steps for a walker possesses diverging moments. Among these unusual diffusion processes, Lévy flights have been extensively studied on lattices and continuous media [1,2] as they can display superdiffusion, so that the variance of the distance covered during the process grows superlinearly x 2 (t) ∝ t γ with γ > 1 at odds with linear diffusion for the Brownian motion [3,[START_REF] Klages | Anomalous transport: foundations and applications[END_REF]. This enhanced diffusion entails an efficient exploration of the space in which the diffusion process takes place: thus both in natural contexts and in artificial ones Lévy flights have emerged as a strategic choice for such an exploration and for search strategies [START_REF] Yang | Proceedings of the World Congress on Nature & Biologically Inspired Computing[END_REF][START_REF] Sharma | [END_REF][7][8][9][10][START_REF] Viswanathan | The physics of foraging: an introduction to random searches and biological encounters[END_REF][START_REF] Méndez | Stochastic foundations in movement ecology: anomalous diffusion, front propagation and random searches[END_REF][START_REF] Brockmann | [END_REF][14][15][16][17][18][19][20][21]. In the case of Lévy flights, the whole process relies on the divergence of the second moment of the jump probability distribution P ( ), i.e. the probability to perform a jump of length . Therefore the walker is allowed to perform very long jumps, the flights, which give, as macroscopic effect, the aforementioned superlinear growth of the total displacement variance x 2 (t) ∝ t γ [START_REF] Klafter | First steps in random walks: from tools to applications[END_REF][START_REF] Shlesinger | [END_REF].

On the other hand, if we focus on the temporal properties of the diffusion, we can introduce for the walker a waiting time probability distribution ψ(t) determining the probability of jumping after a time t has elapsed since the last move. It is straightforward to see that, assuming its first moment is divergent, a subdiffusive behaviour can emerge due to the occurrence of very long waiting times that slow down the dynamics, i.e. x 2 (t) ∝ t β with β < 1 [START_REF] Klafter | First steps in random walks: from tools to applications[END_REF][START_REF] Shlesinger | [END_REF]. These two ingredients, the jump length and the waiting time distributions can be blended to create a richer phenomenology as it is possible to steer from the subdiffusive regime to the superdiffusive one by tuning the power law distribution exponents of the jump and waiting time probabilities [START_REF] Klafter | First steps in random walks: from tools to applications[END_REF]24].

In the framework we just described, anomalous diffusion arises from such a choice of the probability distributions for jumps and rest times but it could be unleashed by other properties of the walkers' motion. In this work, we adopt precisely this perspective: in our model we rely on setting microscopic rules for the walker's displacement so that each "flight" is seen as the result as a series of n unitary very small hops, as in Fig. 1. Anomalous diffusion will therefore stem without the need of an a priori knowledge of the jump length distribution, as in the canonical Lévy flight frame, but it shall be the macroscopic manifestation of such a fragmented and microscopic walk.

The fundamental pivot for the analysis will thus be to relate these microscopic displacements with a macroscopic jump probability distribution P ( ). For a sake of simplicity, we investigate this relation on a 1-D chain where we derive an analytical form for the P ( ) distribution as well as an explicit formula for the displacement variance x 2 (t) . However, we would like to stress that our results could be extended to a more general setting of higher dimensional regular lattices. Our main result will be that, under suitable conditions on the elementary micro-steps distribution p n , the walker can indeed exhibit nonlinear diffusion.

The paper is structured as follows: in Sec. II we introduce the model and we demonstrate that the probability distribution of the jumps P ( ) can display a divergent second moment. Then, in Sec. III, we calculate the probability distribution for the walker that, having in the asymptotic limit a stretched exponential form, leads to a superdiffusive behaviour. In Sec. V we show some numerical simulations to display the Lévy form of the probability distributions and we conclude in Sec. VI with some final remarks.

II. THE MODEL

In our model, we consider walkers moving on a 1-D lattice able of performing elementary steps of unitary length, say +1 and -1, both with equal probability 1/2, as shown in Fig. 1. At each time step, the walker is able to perform n such elementary steps, where n is extracted by some probability distribution function p n . In the following we will assume the latter to follow a power law distribution of exponent γ > 1:

p 0 ∈ [0, 1] and p n = C γ n γ ∀n ≥ 1 , (1)
being C γ = (1 -p 0 )/ζ(γ) a normalising factor and ζ(γ) the Riemann ζ-function. If the probability of not performing any elementary jumps p 0 > 0 then the walker can remain stuck in its current position without doing any elementary steps; on the other hand if p 0 = 0 the walker, at each time step, always performs some elementary jumps, whose possible outcome may eventually be returning to its starting position. As we sketched in the Introduction, the pivotal passage for the analysis is to determine the probability π(m) to perform a total jump of length m for some m ∈ Z in a time step. Assume the walker performs n elementary steps, then the probability of making k ≥ 0 steps in the positive direction, and thus n -k in the negative one, is given by a binomial process 1 2 n n k , hence the total length will result to be m = k -(n -k) = 2k -n. In conclusion we can found:

π(m) = n≥1 p n 2 n n n+m 2 + p 0 δ m,0 , (2)
being the last term the probability of performing a total jump of length m = 0 because the walker did not move at all. The probability to have m = 0 is composed by this term and an additional one, given by n≥1 pn

2 n n n 2 , 1/2 -1 +1 1/2
Elementary move (a) which accounts for the case the walker makes an even number of elementary steps whose total sum is equal to 0. Let us observe that, given m, not all the values of n do contribute to the sum: to ensure the positivity of the binomial coefficient, we must require n ≥ |m| and their sum should be an even number, n + m = 2k, i.e. they should be both odd or even at the same time. The function π(m) is even, as we demonstrate in Appendix A; we can thus restrict ourselves to m ≥ 0 and rewrite Eq. ( 2) for even integers m = 2l as follows:

π(2l) = h≥l p 2h 2 2h 2h h + l ∀l ≥ 1 , (3)
(note that h = 0 is not allowed in the sum because it is taken into account thanks to the term p 0 ) and the case m = 0 reads

π(0) = h≥1 p 2h 2 2h 2h h + p 0 . (4)
For odd integers m = 2l -1 we obtain

π(2l -1) = h≥l p 2h-1 2 2h-1 2h -1 h + l -1 ∀l ≥ 1 . (5)
Having computed the probability π(m), we now focus on its momenta, in particular the second one as its divergence is known to cause the departure from normal diffusion [START_REF] Klafter | First steps in random walks: from tools to applications[END_REF]. Let (X i ) i≥1 be independent random variables such that P (X i = m) = π(m), that is X i is the displacement of the walker at the i-th jump, then one can define T k = X 1 + • • • + X k to be the walker position after k time steps. Because of the parity property of π(m) one gets X k = 0, an thus T k = 0 for all k ≥ 0 (see Appedix A). Using this last remark one can compute the mean square deviation (MSD) as

E(T 2 k ) = i≤k E(X 2 i ) and thus E(X 2 i ) = m m 2 π(m) = 2 m≥1 m 2 n≥m p n 2 n n n+m 2 = n≥1 p n 2 n-1 n m=1 m 2 n n+m 2 (6)
where we used the definition of π(m) and we rearranged the terms in the sum. This latter expression acquires a far simpler form (see Lemma 1 in Appendix A and the probability distribution Eq. ( 1)):

E(X 2 i ) = n≥1 np n = n≥1 C γ 1 n γ-1 (7)
and thus

E(X 2 i ) = C γ /C γ-1 < +∞ if γ > 2 +∞ if 1 < γ ≤ 2 . (8)
In conclusion, if γ > 2 the walker undergoes a linear diffusion process, E(T 2 k ) = kC γ /C γ-1 . On the other hand, if 1 < γ ≤ 2 we cannot conclude anything using the previous analysis; to overcome this difficulty we will consider separately the case 1 < γ ≤ 2 in the next section. Before proceeding, we would like to stress that the existence of an interval for the γ parameter in which the second moment diverges is a crucial passage: in our model the walker performs only local moves without any a priori knowledge of the length it is meant to cover with a jump. This fact paves the way to a generalization to contexts in which the space underneath the walker is highly inhomogeneous as we not necessarily require a metric to define the π(m). Therefore, the second moment divergence in our case of study emerges from the interplay of functional form for the p n and the topology, in this case a 1-D lattice and, as we show in the next section, this divergence reverberates on the probability distribution itself.

III. DISCRETE TIME L ÉVY FLIGHTS

Although Eq. ( 8) has proven the divergence of the MSD when the p n ∼ 1/n γ with 1 < γ ≤ 2, we do not possess so far any information on how, from a functional perspective, this divergence impacts on the probability distribution. Let us define P k (d) the probability for the walker to be at distance d from the starting position after exactly k time steps, that is P k (d) = P (T k = d). Then using the independence of each jump one can derive the following relation:

P k+1 (d) = m P k (d -m)π(m) , (9)
that is the probability to be at distance d at step k + 1 is given by the probability to be one step before at some position d -m and then make a jump of length m. To disentangle this convolution is customary to pass in the Fourier space:

Pk (θ) = d P k (d)e idθ and λ(θ) = m π(m)e imθ .
(10) Hence using Eqs. ( 9) and (10), we obtain

Pk+1 (θ) = Pk (θ)λ(θ) , (11)
from which, by iteration, the following expression results

Pk+1 (θ) = P0 (θ) (λ(θ)) k+1 , (12)
and, applying the inverse Fourier Transform, one can recover P k (d) from

P k (d) = 1 2π ∞ -∞ (λ(θ)) k e -2πidθ dθ . (13)
Eq. ( 13) illustrates how, from the behaviour of λ(θ) for θ → 0, one can deduce the behaviour of Pk (θ) and thus of P k (d) in the asymptotic limit of large |d|. We lever here this standard result to circumvent the divergence in Eq. ( 8) and, in order to unveil the divergence rate of the MSD, we shall focus on the behaviour of λ(θ) for small θ in the following. As we detail in Appendix B, we are able to explicitly cast it in the form

λ(θ) = n≥1 p n (cos θ) n + p 0 . (14)
Let us define s = cos θ and µ(s) = λ(cos(θ)), then, using the chosen form for p n , we can rewrite Eq. ( 14) as:

µ(s) = C γ n≥1 s n n γ + p 0 . ( 15
)
To determine the dependence on s in the sum we use the following approximation:

n≥1 s n n γ ∼ ∞ 1 s x x -γ dx , (16)
for any s ∈ (0, 1)-let us remember that we are interested in θ → 0 and thus s → 1 --we can define y = -x log s > 0 and thus change the integration variable form x to y:

∞ 1 s x x -γ dx = = (-log s) γ-1 ∞ -log s e -y y -γ dy = (-log s) γ-1 I γ (s) , (17)
where I γ (s) is defined by the last equality. We note that for s → 1 -then I γ (s) → I γ (1) = Γ(1 -γ). So, in conclusion we obtain

µ(s) -1 ∼ C γ Γ(1 -γ)(-log s) γ-1 for s → 1 -, (18)
where we used the fact that µ(1) = 1. Back to λ(θ) we obtain for θ → 0

λ(θ)-1 ∼ C γ Γ(1-γ)(-log cos θ) γ-1 ∼ C γ Γ(1 -γ) 2 γ-1 θ 2(γ-1) (19) being λ(0) = 1, cos θ ∼ 1 -θ 2 /2 and -log(1 -θ 2 /2) ∼ θ 2 /2.
We thus have, from Eq. ( 12), for small θ

Pk (θ) = λ k (θ) ∼ 1 -A γ θ 2(γ-1) k , (20)
where A γ = -C γ Γ(1 -γ)/2 γ-1 > 0. Therefore, in the limit of large k, the above expression tends to the stretched exponential form typical of Lévy flights characteristic function:

Pk (θ) ∼ e -kAγ θ 2(γ-1) . ( 21
)
The inverse Fourier Transform of the characteristic function in Eq. ( 21) does not have a straightforward analytical expression and, being non-analytic, the evaluation of the MSD using the standard rule

d 2 = ∂ 2 ∂ 2 θ Pk (θ)
θ=0 is impeded because the latter expression diverges. It is nevertheless possible to exploit the self-similarity property of the distribution (21) in order to obtain a scaling relation showing the impact of the local exponent γ on the probability distribution. Using P ak (a 1/g d) = a -1/g P k (d), we can thus recast the inverse Fourier Transform

P k (d) = 1 2π ∞ -∞ e -2πidθ-kDγ θ 2(γ-1) dθ , (22)
from which we get

P k (d) = k -1/2(γ-1) Π d k 1/2(γ-1) , ( 23
)
where Π is a function of the sole variable d/k 1/2(γ-1) . We have thus found that the PDF is shaped by the dynamical exponent µ = 1 2(γ-1) ; therefore, as we anticipated in the previous section, the rule governing the size of the walker's microsteps cascade resonates in the overall diffusion process. As a closure to the present section we would like to make a remark on the probability of not making any microscopic move p 0 . Let us observe that the walker will always perform Lévy flights for any p 0 ∈ [0, 1), being the impact of p 0 only on C γ , more precisely C γ → 0 when p 0 → 1, but not on the exponent 2(γ -1). Only in the extremal case p 0 = 1 the walk degenerates into an absence of movement.

IV. CONTINUOUS TIME APPROACH

In the previous section we considered a discrete time process in which the steps occurred at a regular pace. In this section we extend our analysis introducing in our description the waiting time probability distribution, which allows the walker to wait after n micro-steps at the reached position for a time interval t before hopping again. Thus our process is now composed by two moves: a waiting time, whose length is weighted by a distribution ψ(t), and a "dynamic"phase in which n elementary steps are instantaneously performed. In our approach, we consider the probability distributions π(m) and ψ(t) as independent and the dynamic phase can be interpreted as the flights in our model since it does not take time, similarly to the classical Lévy flights. With these hypotheses, the derivation of the final probability distribution in Fourier-Laplace space P (θ, s) is straightforward in the Continuous Time Random Walk (CTRW) frame [START_REF] Klafter | First steps in random walks: from tools to applications[END_REF], but we detail here the passages for the sake of completeness. We thus assume that the walker starts at t = 0 and let ψ k (t) be the probability distribution function of the occurrence of the k-th jump at time t = t 1 + • • • + t k where t i is the waiting time drawn at the i-th jump.One clearly has

ψ k (t) = t 0 ψ k-1 (t )ψ(t -t ) dt . ( 24
)
This equation leads, passing in Laplace space (the complete derivation can be found in Ref. [START_REF] Klafter | First steps in random walks: from tools to applications[END_REF]), to an expression for χ k (t), which is the probability to make exactly k jumps up to the time t. In Laplace space, χk (s) reads:

χk (s) = ψ(s) k 1 -ψ(s) s . ( 25
)
Now that the distribution χ k accounts for the non-linear relation between steps and time, we can proceed to include it within the definition of the P (d, t), i.e. the probability for the walker to be at distance d from the origin (the initial point at time t = 0) at time t. We observe that this probability is a generalisation of the previously defined P k (d): of course, in case all the waiting times are equal to the duration of the rest period, τ , then P (d, t)

reduces to P k (d)
where k = t/τ , as we had in Eq. ( 23). Using our starting hypotheses, i.e. that the jumps are costless in time and the waiting time is uncorrelated with the jumps, we can write

P (d, t) = k≥0 P k (d)χ k (t) , (26)
meaning that the probability P (d, t) is the probability to be at d after exactly k steps times the probability to have performed k steps in the time interval t. Using once again the Laplace transform for time, Fourier for space and the result of the previous section we arrive at the classical result [START_REF] Klafter | First steps in random walks: from tools to applications[END_REF]:

P (θ, s) = k≥0 λ k (θ) χk (t) = k≥0 λ k (θ) ψk (s) 1 -ψ(s) s (27) = 1 -ψ(s) s 1 1 -λ(θ) ψ(s)
.

Macro movements (a) We thus have that the asymptotic behaviour of P (d, t) shall be governed, in the d, t → ∞ limit, by the moments of the λ(θ) and ψ(s) in the corresponding limit s, θ → 0 in the Fourier-Laplace space. Therefore we combine the approximation of λ(θ) in Eq. ( 19) with a waiting time distribution assumed to have a power-law form ψ(t) ∼ τ α /t 1+α with 0 < α < 1 for t → ∞. In order to investigate if the waiting time distribution interferes with the PDF's profile, we focus, for the jump part, on the interesting case where 1 < γ ≤ 2 as we demonstrated in the previous section that it leads the second spatial moment to diverge. Passing to Fourier-Laplace space, the Laplace transform of ψ(t) reads ψ(s) ∼ 1 -τ α s α for small s by virtue of the Tauberian theorem and, substituting the approximation of λ(θ) and ψ(s) in Eq. [START_REF] Klafter | First steps in random walks: from tools to applications[END_REF] we obtain P (θ, s) ∼ τ α s α-1 τ α s α + A γ θ 2(γ-1) .

(

) 28
We then extrapolate the scaling behaviour in the same fashion we derived Eq. ( 23), where now both the exponents α and γ intervene in the temporal scaling [25] P (d, t) = t -α/2(γ-1) Π d t α/2(γ-1) .

(29)

V. NUMERICAL SIMULATIONS

The aim of this section is to present some numerical results to support the theory presented above. We are left now with the numerical evaluation of the probability distribution to confirm the impact brought by the local exponents α and γ on the overall diffusion process. As for the discrete case, in Fig. 3, we show how the γ exponent governs the behaviour of the probability distribution: indeed as soon as γ > 2, the second moment of the π(d) becomes finite and the PDF tends to a Gaussian distribution (Fig. 3a). On the other hand, in the γ 2 regime, the PDF clearly exibits the fat-tailed Lévy functional form (Fig. 3b). In the Lévy case, the variable rescaling that induces the curves collapse in Fig. 3 bears the mark of the γ exponent since we rescale with respect to ξ = d k 1/2(γ-1) , as obtained in Eq. 23. For the continuous time regime, both the exponents α and γ intervene in the shape of the PDF as shown by Eq. 29; therefore in Fig. 4 the superposition of the PDF at different times emerges in the same fashion as before once the rescaling is done with respect to the variable ξ = d t α/2(γ-1) . As a closing note, it is worth mentioning that other methods exists to tame the numerical instability and investigate, albeit indirectly, the theoretical scaling of the P (d, t) moments such as computing the fractional moments x δ with 0 < δ < µ ≤ 2 and µ = 2(γ -1) [1], the mean of the displacements x i (t) logarithm, called the geometric mean rg [START_REF] Lubashevsky | [END_REF] and, finally, computing the probability density averaged within a box with time depending bounds [L 1 t 1/µ , L 2 t 1/µ ] [28].

VI. CONCLUSION

Concluding, in this work we introduced a random walk model igniting a Lévy flight type of behaviour and leading to superdiffusion on a one dimensional lattice. The specificity of this model is to impose a microscopic condition on the walk, with no need for an a priori knowledge of the topology in order to perform the jumps. In our approach, one jump event corresponds to an "avalanche" of n elementary steps, whose size n is distributed according to a probability distribution p n . We then demonstrated that a power law form p n ∼ 1/n γ entails the divergence of the second moment of the jumps length distribution π(m) when 1 < γ 2. Starting from this divergence, we derived, in Sec. III, the probability distribution Pk (θ) in Fourier space which has the characteristic stretched exponential form. We furthermore introduced the possibility for the walker to stay put on a in Sec. IV, showing how the α exponent of the waiting time PDF determines, along with the one of the avalanches γ, the form of the probability distribution P (d, t). Finally, in Sec.V we confirmed through direct numerical simulation the analytical behaviour of the latter, displaying the tails' scaling. On a closing note, we would like to stress that the approach itself is independent of the p n functional form and that it could be generalised to other distributions. The actual meaningful information carried by the p n is the creation of a divergence in the jumps second moment computed using the π(m) distribution. It is worth of note that this divergence stems from the interplay of both the p n shape and the 1D topology; therefore a careful choice of the former might be a way to create anomalous diffusion in more general network topologies. Widening further our perspective, the walk described in this paper could be used when the underlying space does not possess a proper metric and is small-world, as in the case of complex networks [START_REF] Newman | Networks: an introduction[END_REF], such that the probability to perform a walk at a certain distance is not univocally defined. In that case, adopting a local perspective for the walker dynamics might prove useful to test the notion of anomalous diffusion [START_REF] Riascos | [END_REF]. Another possible application would the modelling of diffusion on temporal networks [START_REF] Newman | The structure and dynamics of networks[END_REF], especially in the presence of burstiness [START_REF] Holme | [END_REF] and the number of events within a time window can be broadly distributed, possibly under the form of trains of events [33].

FIG. 1 .

 1 FIG.1. The walker is allowed to perform at each step, n elementary moves each one having equal probability 1/2 to be in the positive or negative direction (Panel a). In panel b) we present a possible avalanche, the walker performs n = 6 hops, two in the positive direction and 4 in the negative one, for a total effective displacement of m = -2.

FIG. 2 .

 2 FIG. 2.In the continuous time frame, at each time ti the walker can, as before, perform n microscopic moves: in panel a) these macro-movements result in displacements of m = 2, m = 5 and m = 2. If we introduce the waiting times, the displacements are interspersed by waiting intervals: in panel b) the walker stays put from t1 to t 1 and from t2 to t 2 .

FIG. 3 .

 3 FIG. 3. Probability distribution P k (d) in the discrete case for (a) γ = 2.1 and (b) γ = 1.9 for N = 10 5 where µ = 1 2(γ-1)

FIG. 4 .

 4 FIG. 4. Probability distribution P (t) in the continuous time case for (a) α = 0.8, γ = 1.5 and (b) α = 0.3, γ = 1.7 for N = 10 5 where µ = α 2(γ-1) and ξ =

ACKNOWLEDGMENTS

The work of T.C. and R.L. presents research results of the Belgian Network DYSCO (Dynamical Systems, Control, and Optimization), funded by the Interuniversity Attraction Poles Programme, initiated by the Belgian State, Science Policy Office.

Appendix A: Walk properties

The symmetry of the walk reflects in the parity of π(m), i.e. π(-m) = π(m):

Therefore, considering the first moment is trivially E(X i ) = 0 for all i ≥ 1:

Hence on average the walker doesn't move from the initial position E(T k ) = i≤k E(X i ) = 0. On the other hand, for the MSD, the last equality in Eq. ( 6) gives

where

In order to elucidate its behaviour, we shall use the following Lemma Lemma 1. Let us define for all n ≥ 1

Then one has

Proof. Let us consider separately the case n = 2l (even) and n = 2l -1 (odd).

From the definition of c n and using the parity assumption on m and n we can rewrite m = 2h, for some 1 ≤ h ≤ l, and thus

The following relations hold to be true

Let us develop the definition of c 2l to be able to use the previous relations:

where Eqs. (A8) have been used to pass from the first line to the second one. Let us rewrite the rightmost term using the change of summing index j -l = -h:

where we used the fact that 2l l-h = 2l l+h . In conclusion we have thus found

The case n = 2l -1 (odd) can be handled in the same fashion, thus concluding that

From this equality and the definition of the probability distribution (Eq. ( 1)), it can be obtained

In this appendix we detail the derivation of Eq. ( 14) for the function λ(θ). Firstly, we observe that using the parity of π(m) one can write its Fourier Transform λ(θ) as

but it is not possible to use the Taylor development cos(mθ) = 1-m 2 θ 2 /2+. . . because in the present case, 1 < γ ≤ 2, we already know that m≥1 π(m)m 2 diverges. We thus turn to definition of π(m) and write

where g n (θ) is defined using the last equality. For this function g n (θ) holds the following lemma Lemma 2. Let g n (θ) be defined by Eq. (B2), then ∀l ≥ 1.

Proof. Let us consider once again separately the case n = 2l (even) and n = 2l -1 (odd) for some l ≥ 1. Then one can rewrite

Let us rewrite the sum for the even case using the variable j = l + h and the sum for the odd case with the variable j = l + h -1. Then one has, for the even case (the odd case can be treated exactly in the same manner):

Let add and remove in both sums the number of terms up to h = 0:

Rewriting cos x = (e ix + e -ix )/2

and using the definition of binomial coefficient (1 + x) n = n k=0 n k x k we get:

and some manipulations give

replacing in both sums j = l -h (and isolating the term j = 0 in the sum for even n) we get

and using the defintion of g 2l and g 2l-1 we obtain: and obtain Eq. ( 14).


Franc ¸ois Boulier

¸ois Franc

Lemaire
email: lemaire@lifl.fr

COMPUTING CANONICAL REPRESENTATIVES OF REGULAR DIFFERENTIAL IDEALS

Keywords: differential algebra, Rosenfeld-Gröbner, canonical representatives, lextriangular, regular differential ideal, characteristic sets, characteristic presentations

In this paper, we give three theoretical and practical contributions for solving polynomial ODE or PDE systems. The first one is practical: an algorithm which improves the purely algebraic part of Rosenfeld-Gröbner (the polynomial ODE or PDE systems simplifier which is the core of the Maple 5.5 diffalg package). It is a variant of lextriangular but does not need any Gröbner basis computation. The second one is theoretical: a characterization of the output of Rosenfeld-Gröbner and a clarification of the existing relationship between algebraic and differential characteristic sets. The third one is theoretical as well as practical: an algorithm to compute canonical representatives of differential polynomials modulo regular differential ideals without any use of Gröbner bases. This algorithm simplifies the theory (somehow a "pedagogic" contribution) but permits us also to perform easily linear algebra over the base field in the factor differential ring defined by a regular differential ideal.

INTRODUCTION

We make precise in next sections some of the terms used in this introduction.

Rosenfeld-Gröbner [START_REF] Boulier | Étude et implantation de quelques algorithmes en algèbre différentielle[END_REF][START_REF] Boulier | Representation for the radical of a finitely generated differential ideal[END_REF][START_REF] Boulier | Computing representations for radicals of finitely generated differential ideals[END_REF] is a simplifier for systems of polynomial differential equations (ordinary or with partial derivatives). It solves a theoretical problem of differential algebra (deciding membership in the radical of finitely generated differential ideals) and furnishes as a byproduct tools to solve systems of polynomial ODE and PDE. The situation is similar to that of the Buchberger algorithm [START_REF] Buchberger | An Algorithm for Finding a Basis for the Residue Class Ring of a Zero-Dimensional Polynomial Ideal[END_REF], which solves a theoretical problem of commutative algebra (deciding membership in polynomial ideals) and furnishes as a byproduct tools to solve systems of polynomial equations.

It is the heart of the diffalg package which is part of Maple 5.5 standard library.

Technically, given a system Σ of differential polynomials, Rosenfeld-Gröbner outputs finitely many characteristic presentations C1, . . . , Ct which are particular cases of triangular systems. Each characteristic presentation Ci represents a differential ideal, denoted [Ci] : H ∞ C i . The radical p [Σ] of the differential ideal generated by Σ is the intersection

p [Σ] = [C1] : H ∞ C 1 ∩ • • • ∩ [Ct] : H ∞ C t .
Roughly speaking, Rosenfeld-Gröbner consists in two steps.

The first step is the differential step which transforms Σ as finitely many systems A = 0, S = 0 of equations and inequations, called regular differential systems.

p [Σ] = [A1] : S ∞ 1 ∩ • • • ∩ [A t ′ ] : S ∞ t ′ .
These systems satisfy the hypotheses of the key Rosenfeld's lemma which "reduces differential problems to purely algebraic ones". Each regular differential system A = 0, S = 0 represents a regular differential ideal denoted [A] : S ∞ . Its set of equations A is triangular but A is not necessary a characteristic presentation of [A] : S ∞ .

The second step is the purely algebraic step which transforms a regular differential system A = 0, S = 0 into finitely many characteristic presentations C1, . . . , C t ′′ satisfying

[A] : S ∞ = [C1] : H ∞ C 1 ∩ • • • ∩ [C t ′′ ] : H ∞ C t ′′ .
It is applied separately over each regular differential system produced by the first step.

The fact that computing characteristic presentations from regular differential systems is a purely zerodimensional algebraic problem is only claimed in [8, page 35] with no proofs. The first complete proofs of that claim are given in [START_REF] Hubert | Factorization free decomposition algorithms in differential algebra[END_REF].

The first result of this paper is an efficient algorithm, called regCharacteristic which performs the second step without computing any Gröbner basis.

The regCharacteristic algorithm first reduces the problem to a dimension zero problem 1 . It puts in the coefficient field all the derivatives occuring in A and S which are not leaders of any element of A, giving sets A and S. It then calls the subalgorithm satTriangular which computes a set of normalized triangular sets {T1, . . . , Tt} which (nearly immediately) give the desired characteristic presentations {C1, . . . , Ct}.

The satTriangular subalgorithm is a very close variant of the lextriangular algorithm [20, page 129, algorithm D5lextriangular] [START_REF] Maza | Polynomial gcd computations over towers of algebraic extensions[END_REF] and [25, page 133] which applies the D5 [START_REF] Dora | About a new method for computing in algebraic number fields[END_REF] process. Let's quote [20, page 129]: lextriangular is given a Gröbner basis B of a zero-dimensional ideal of polynomials in X1, . . . , Xn, sorted by increasing leading monomials for the lexicographical ordering such that X1 < • • • < Xn and produces on output a finite family of normalized triangular sets T1, . . . , Tt such that

V (B) = V (T1) ∪ • • • ∪ V (Tt) (1)
where V (Ti) denotes the set of the common zeroes of the elements of Ti in the algebraic closure of the base field.

Our implementation of satTriangular is directly inspired from [START_REF] Maza | Polynomial gcd computations over towers of algebraic extensions[END_REF][START_REF] Maza | Calculs de Pgcd au-dessus des Tours d'Extensions Simples et Résolution des Systèmes d' Équations Algébriques[END_REF] (we do not claim any algorithmic improvement w.r.t. these papers): it performs the subresultant algorithm in the factor ring defined by the already built triangular sets instead of performing it generically and specializing its result afterwards.

There are however important theoretical differences between satTriangular and lextriangular: we start from a system of polynomial equations A = 0, S = 0 which is not a Gröbner basis2 ; more important, we want much stronger properties for our output than (1) i.e. properties P1 to P5 stated in section 4. These properties just do not hold for the output of lextriangular in general and we were led to write complete proofs.

We have implemented two versions of Rosenfeld-Gröbner based on regCharacteristic. One in Maple 5.5 and one in C++. We compare existing methods (which are written in Maple) with our Maple version.

Let us compare our work with existing methods.

Wang and Li solve the second step in [START_REF] Li | Coherent, regular and simple systems in zero decompositions of partial differential systems[END_REF] by using the SimSys [START_REF] Wang | Decomposing polynomial systems into simple systems[END_REF] algorithm which handles general polynomial systems.

They claim [21, p. 59] that a more specialized algorithm than SimSys could be applied. Our regCharacteristic is such an algorithm.

In [START_REF] Boulier | Computing representations for radicals of finitely generated differential ideals[END_REF], this second step is solved by computing first a Gröbner basis of the localized ideal S -1 (A). This is expensive, does not take into account the fact that A is already triangular and explicits the inverses of the elements of S which are not needed at all. Algorithm 7.1 in [START_REF] Hubert | Factorization free decomposition algorithms in differential algebra[END_REF] applies exactly the same principles as regCharacteristic. The only difference is that it computes a Gröbner basis of S -1 (A) instead of calling satTriangular.

It suffers therefore of the drawbacks mentioned above for [START_REF] Boulier | Computing representations for radicals of finitely generated differential ideals[END_REF]. Note also that [START_REF] Hubert | Factorization free decomposition algorithms in differential algebra[END_REF] is the first to prove completely that computing characteristic presentations from regular differ-ential systems is a purely algebraic (zerodimensional) problem. This was only claimed in [8, page 35].

Observe that testing the invertibility of differential polynomials modulo triangular sets in order to build characteristic sets was already considered in [START_REF] Maârouf | Étude de Quelques Problèmes Effectifs en Algèbre Différentielle[END_REF][START_REF] Bouziane | Unmixed-Dimensional Decomposition of a Finitely Generated Perfect Differential Ideal[END_REF]. The method of these authors is different from ours, at least because it is based on Gröbner bases computations [9, page 7] and [23, page 29].

The second result of this paper consists in two theorems (theorems 2 and 3) which clarify the relationship between characteristic presentations, characteristic sets and regular differential ideals. Our third result is a new method to compute canonical forms of differential polynomials modulo a regular differential ideal.

Computing canonical forms of differential polynomials modulo a given differential ideal a is a real issue for the set of canonical forms modulo a often forms (here it is the case) a vector space over the base field K of the equations. Using canonical forms we thus can look for linear dependencies over K modulo a by easily performing linear algebra in the factor ring. This is one of the main ideas carried out by the important FGLM algorithm [START_REF] Faugère | Efficient computation of Gröbner bases by change of orderings[END_REF] in the context of Gröbner bases and lifted to regular differential systems in [START_REF] Boulier | Efficient computation of regular differential systems by change of rankings using Kähler differentials[END_REF].

This problem of computing canonical forms modulo polynomial differential equations was, as far as we know, only addressed in [START_REF] Boulier | Efficient computation of regular differential systems by change of rankings using Kähler differentials[END_REF] and required the computation of a Gröbner basis. The method we give in section 8 is based on triangular sets and pseudo reduction only. It solves a problem left open in [START_REF] Boulier | Efficient computation of regular differential systems by change of rankings using Kähler differentials[END_REF].

Our method applies for ordinary (non differential) polynomials modulo regular sets [START_REF] Kalkbrener | A Generalized Euclidean Algorithm for Computing Triangular Representations of Algebraic Varieties[END_REF][START_REF] Aubry | On the theories of triangular sets[END_REF] of polynomial equations too and we believe it could be interesting in this context also. Computing canonical forms imposes to compute the inverses of the initials of the elements of the set. These al-gebraic inverses computations can be performed only once by making the set strongly normalized. This is probably very CPU expensive but is interesting at least for pedagogic reasons.

DEFINITIONS AND NOTATIONS 2.1 Commutative algebra

Let X be an ordered alphabet (possibly infinite). A term over X is a power product of elements of X.

Let R = K[X] be a polynomial ring where K is a field. Let p ∈ R \ K be a polynomial. The leader of p, denoted ld p, is the greatest indeterminate x which occurs in p. The polynomial p can be written as

p = a d x d + • • • + a1 x + a0
where d = deg(p, x) and the polynomials ai are free of x.

The polynomial ip = a d is the initial of p. The rank of p is the monomial x d . If x d and y e are two ranks then x d < y e if x < y or x = y and d < e. The separant of p is the polynomial

sp = ∂p ∂x •
The polynomial p is said to be monic if its initial is equal to 1. The set iter(p) of the iterated initials of p is defined as follows: if p ∈ K then iter(p) = ∅ otherwise iter(p) = {p} ∪ iter(ip).

Let A ⊂ R \ K be a set of polynomials. Then IA (resp. SA) denotes the set of the initials (resp. the separants) of its elements. We denote HA = IA ∪ SA. The set A is said to be triangular if its elements have distinct leaders.

Let A be a triangular set. A polynomial p is said to be normalized w.r.t. A if the set of leaders of iter(p) is disjoint from the set of leaders of A. The set A is said to be normalized if every p ∈ A is normalized w.r.t. A \ {p}.

A polynomial p is said to be strongly normalized w.r.t. A if no leader of A occurs in the initial of p. The set A is said to be strongly normalized if every p ∈ A is strongly normalized w.r.t. A \ {p}.

Every strongly normalized triangular set is normalized.

We denote prem(p, A) the pseudo-remainder [17, volume 2, page 407] of p by all the elements of A viewed as univariate polynomials in their leaders.

If R is a unique factorization domain and p ∈ R[X] then p can be written:

p = a0 t0 + • • • + a k t k
where the ti are terms over X and the ai ∈ R. The content of p over R is the gcd of its coefficients:

cont(p) = gcd(a0, . . . , a k )
The primitive part of p over R is the polynomial

pp(p) = p cont(p) •
A polynomial is said to be primitive if it is equal to its primitive part.

If A is a subset of a ring R then (A) denotes the ideal generated by A. Let a be an ideal of R. Then √ a denotes the radical of a. If S = {s1, . . . , st} is a finite family of elements of R then the saturation a : S ∞ of a by S is the ideal:

a : S ∞ = {p ∈ R | ∃a1, . . . , at ∈ N s.t. s a 1 1 • • • s a t t p ∈ a}.

Differential algebra

We only provide a short presentation. The reference books are [START_REF] Ritt | Differential Algebra[END_REF] and [START_REF] Kolchin | Differential Algebra and Algebraic Groups[END_REF]. We also refer to the Maple 5.5 diffalg package and thus to the articles [START_REF] Boulier | Representation for the radical of a finitely generated differential ideal[END_REF][START_REF] Boulier | Computing representations for radicals of finitely generated differential ideals[END_REF] which present it.

A derivation over a ring R is a map δ : R → R which satisfies, for every

a, b ∈ R δ(a + b) = δa + δb, δ(a b) = (δa)b + a(δb).
A differential ring is a ring endowed with finitely many derivations which commute pairwise. The commutative monoid generated by the derivations is denoted Θ. Its elements are the derivation operators

θ = δ a 1 1 • • • δ am m
where the ai are nonnegative integer numbers. The sum of the exponents ai, called the order of the operator θ, is denoted ord θ. The identity operator is the unique operator with order 0. The other ones are called proper.

If φ = δ b 1 1 • • • δ bm m then θφ = δ a 1 +b 1 1 • • • δ am+bm m . If ai > bi for each 1 ≤ i ≤ m then θ/φ = δ a 1 -b 1 1 • • • δ am-bm m . A differential ideal a of R is an ideal of R stable under deriva- tion i.e. such that a ∈ a ⇒ δa ∈ a.
Let A be a nonempty subset of R. We denote [A] the differential ideal generated by A which is the smallest differential ideal which contains A.

Differential polynomials

Let U = {u1, . . . , un} be a set of differential indeterminates. Derivation operators act over differential indeterminates giving derivatives θu. We denote ΘU the set of all the derivatives. Let K be a differential field. The differential ring of the differential polynomials built over the alphabet ΘU with coefficients in K is denoted R = K{U }.

A ranking is a total ordering over the set of the derivatives [18, page 75] satisfying the following axioms 1. δv > v for each derivative v and derivation δ, 2. v > w ⇒ δv > δw for all derivatives v, w and each derivation δ.

Fix a ranking. The infinite alphabet ΘU gets ordered. Consider a polynomial p ∈ R \ K. Then the leader, initial, separant . . . of p are well defined. Axioms of rankings imply that the separant of p is the initial of every proper derivative of p.

Let rank p = v d . A differential polynomial q is said to be partially reduced w.r.t. p if no proper derivative of v occurs in q. It is said to be reduced w.r.t. p if it is partially reduced w.r.t. p and deg(q, v) < d.

A set A of differential polynomials is said to be differentially triangular if it is triangular and if its elements are pairwise partially reduced. It is said to be autoreduced if its elements are pairwise reduced.

Every autoreduced set is differentially triangular.

If A is a set of differential polynomials and v is a derivative then Av = {θp | p ∈ A, ld θp ≤ v}. Thus Rv denotes the set of all the differential polynomials having leader less than or equal to v.

Ritt's reduction algorithms

They are generalizations of the Euclidean division algorithm for differential polynomials. One distingues the partial reduction algorithm, denoted partial rem from the full reduction algorithm, denoted full rem. We only give specifications of these algorithms. See [18, page 77] for a more precise description. Let q be a differential polynomial and A be a set of differential polynomials. Let v = ld q and A = A ∩ Rv.

If q = partial rem(q, A) denotes the partial remainder of q by A then q is partially reduced w.r.t. all the elements of A and there exists a power product h of elements of S A such that h q ≡ q mod (Av).

If q = full rem(q, A) denotes the full remainder of q by A then q is reduced w.r.t. all the elements of A, there exists a power product h of elements of H A such that h q ≡ q mod (Av).

Critical pairs

A pair {p1, p2} of differential polynomials is said to be a critical pair 3 if the leaders of p1 and p2 are derivatives of some same differential indeterminate u (say ld p1 = θ1u and ld p2 = θ2u). Assume A is differentially triangular. Then critical pairs(A) denotes all the critical pairs that can be formed with any two elements of A. Let {p1, p2} ∈ critical pairs(A) be a critical pair. Denote θ12 the least common multiple between θ1 and θ2. The ∆-polynomial

∆(p1, p2) is ∆(p1, p2) = s2 θ12 θ1 p1 -s1 θ12 θ2 p2
where s1, s2 denote the separants of p1 and p2. Let A = 0, S = 0 be a system of differential polynomial equations and inequations. The critical pair {p1, p2} is said to be solved by A = 0, S = 0 if there exists a derivative v < θ12u such that

∆(p1, p2) ∈ (Av) : (S ∩ Rv) ∞ .
If full rem(∆(p1, p2), A) = 0 and HA ⊂ S then the critical pair {p1, p2} is solved by A = 0, S = 0.

Regular differential systems

Definition 1. (regular differential systems)4

A differential system A = 0, S = 0 of a differential polynomial ring R is said to be a regular differential system (for a ranking R) if C1 A is differentially triangular, C2 S contains the separants of the elements of A and is partially reduced w.r.t. A, C3 all the critical pairs {p, p ′ } ∈ critical pairs(A) are solved by the system A = 0, S = 0 (coherence property 5 ).

If A = 0, S = 0 is a regular differential system then the ideal [A] : S ∞ (resp. (A) : S ∞ ) is called the regular differential ideal (resp. regular algebraic ideal) defined by the system.

If A = 0, S = 0 is a regular differential system, we call derivatives under the stairs of A the elements of ΘU which are not derivatives of any leader of element of A.

Regular systems enjoy the following properties.

Theorem 1. Let A = 0, S = 0 be a regular differential system of R = K{U }. Let L denote the set of leaders of A and N the set of the derivatives under the stairs of A. Then Proof. [8, Lazard's lemma, the lifting of Lazard's lemma and Rosenfeld's lemma]. See also [START_REF] Seidenberg | An elimination theory for differential algebra[END_REF][START_REF] Rosenfeld | Specializations in differential algebra[END_REF][START_REF] Boulier | Representation for the radical of a finitely generated differential ideal[END_REF][START_REF] Schicho | A construction of radical ideals in polynomial algebra[END_REF][START_REF] Morrison | The Differential Ideal [P ] : M ∞[END_REF][START_REF] Hubert | Factorization free decomposition algorithms in differential algebra[END_REF][START_REF] Li | Coherent, regular and simple systems in zero decompositions of partial differential systems[END_REF].

Characteristic sets

Definition 2. Let a be a differential ideal of R. A set C ⊂ a is said to be a characteristic set of a if C is autoreduced and a contains no nonzero polynomial reduced w.r.t. C.

In the following theorem, the equivalence between the two first items is well-known. This theorem is very close to [27, Théorème 2, page 94] though Ollivier only considers differential prime ideals. The equivalence with the third item may be essentially proved already in [START_REF] Bouziane | Unmixed-Dimensional Decomposition of a Finitely Generated Perfect Differential Ideal[END_REF]. The correspondence with the recent [2, Theorem 6.1] is interesting in itself anyway.

Characteristic presentations

The following definition is different from that of [8, Definition 6.1] but we do believe both definitions are equivalent (we do not prove this claim). Hubert weakens the definition

6 A triangular set C = f1 < • • • < fn is regular in the sense of [2] if the initial of f k does not divide zero modulo (f1, . . . , f k-1 ) : (i1 • • • i k-1
) ∞ for every 1 ≤ k ≤ n. 7 A regular (in the sense of [START_REF] Aubry | On the theories of triangular sets[END_REF]) triangular set C is said to be squarefree if (C) : I ∞ C is radical (in [1, def 4.5.11], the qualifier separable is used in place of squarefree).

of characteristic presentations in [START_REF] Hubert | Factorization free decomposition algorithms in differential algebra[END_REF] (she only imposes D1 and D2) thus looses canonicity properties.

Definition 3. A set C ⊂ K{U } is said to be a charac- teristic presentation of the differential ideal [C] : H ∞ C if D1 the differential system C = 0, HC = 0 is regular, D2 if p ∈ R then p ∈ [C] : H ∞ C ⇔ full rem(p, C) = 0, D3 C is a strongly normalized autoreduced set of K[L, N ]
such that the elements of C are primitive over K[N ] where L denotes the set of leaders of the elements of C and N denotes the other derivatives occuring in C.

Corollary 1. A set C of differential polynomials is a characteristic presentation of [C] : H ∞ C if and only if C is a characteristic set of [C] : H ∞ C which satisfies D3.

THE PROBLEM

We are given a regular differential system A = 0, S = 0 of R. We assume moreover that HA ⊂ S. We want to compute sets C1, . . . , Ct of differential polynomials such that the following conditions hold

A1 each Ci is a characteristic presentation of the differential ideal [Ci] : H ∞ C i , A2 [A] : S ∞ = [C1] : H ∞ C 1 ∩ • • • ∩ [Ct] : H ∞ C t , A3 the intersection is not redundant: if p is a differential prime component of [A]:S ∞ then p is a minimal differ- ential prime of exactly one differential ideal [Ci] : H ∞ C i .

REGCHARACTERISTIC

Denote X the set of the derivatives occuring in A ∪ S and L ⊂ X the set of the leaders of the elements of A and N = X \ L. Denote G = K(N ) the ring obtained by putting the elements of N in the base field of the differential polynomial ring. The algorithm involves three steps.

1. Transform the system A = 0, S = 0 into a system A = 0, S = 0 of G[L]. That step is purely formal.

2. Apply the algorithm satTriangular (described later) over A = 0, S = 0. This algorithm returns a possibly empty set of squarefree normalized 8 autoreduced sets {T1, . . . , Tt} of G[L] satisfying the following properties:

P1 if the set is empty then (A) : S ∞ = G[L], P2 (A) : S ∞ = (T1) ∩ • • • ∩ (Tt), P3 if i = j then (Ti) + (Tj) = (1), P4 for every 1 ≤ i ≤ t we have (Ti) = (Ti) : H ∞ T i , P5 ld Ti = ld A.

Proof of the algorithm

The part of regCharacteristic that we prove in this section is shared together with [START_REF] Hubert | Factorization free decomposition algorithms in differential algebra[END_REF]Algorithm 7.1]. Proofs can thus be found in [15, Theorem 3.10 and 6.2]. We give them anew to make this article selfcontained.

Let φ be the canonical ring homomorphism

K[X] → G[L]. Lemma 1. (A) : S ∞ = (C1) : H ∞ C 1 ∩ • • • ∩ (Ct) : H ∞ C t .
Proof. The relation above is obtained by applying φ -1 over P2 componentwise. Indeed, φ -1 preserves intersections by [ 8 Observe in this case, every normalized set is strongly normalized for its elements are monic. A1)

Let 1 ≤ i ≤ t be an index. The set Ci is a characteristic presentation of the ideal [Ci] : H ∞ C i .
Proof. The set Ci is autoreduced. All its elements are strongly normalized for they are obtained by multiplying monic polynomials by elements of K[N ]. They are primitive over K[N ] by construction. Thus Ci satisfies D3. Condition D1 holds by lemma 3. Condition D2 thus holds by theorem 1 (3) and lemma 4.

Lemma 5. A differential prime p is minimal over [A]:S ∞ iff it is minimal over some [Ci] : H ∞ C i .
Proof. By lemmas 2, 3 and theorem 1 (5).

Proposition 2. (condition A2)

[A] :

S ∞ = [C1] : H ∞ C 1 ∩ • • • ∩ [Ct] : H ∞ C t .
Proof. The proposition is a corollary of lemma 5.

Proposition 3. (condition A3) If p is a differential prime component of [A] : S ∞ then p is a differential prime component of exactly one differential ideal [Ci] : H ∞ C i . Proof. Let p be a minimal differential prime of [Ci]:H ∞ C i and of [Cj] : H ∞ C j . Then b = p ∩ K[X]
is a minimal prime of (A) : S ∞ by lemma 5. Then (φb) is a minimal prime of (Ti) and of (Tj). Property P3 implies i = j.

THE INVERT SUBALGORITHM

Computing the normalized sets Ti amounts to normalizing the polynomials of A, which consists in inverting the initials of the polynomials of A.

Denote L = {X1, . . . , Xn}.

Specification of invert

We describe precisely the inputs and outputs of the algorithm invert:

Inputs of invert:

• p ∈ G[L], a non zero polynomial • T = {p1, . . . , p k } a normalized triangular set of the ring G[X1, . . . , X k ] with k ≤ n.
Outputs of invert:

• either the inverse q ∈ G[L] such that pq = 1 (mod (T )).

We say invert has found the inverse of p modulo the ideal (T ).

• or a triple (j, g, h) (if invert could not compute an inverse) with j ≤ k such that ld(g) = ld(h) = Xj g and h are monic pj = g h (mod (p1, . . . , pj-1))

Algorithmic scheme of invert

The algorithms invert and ExtEuclid (see below) respectively are simplified versions of the functions QuasiRecipElseSplit and extendedSubResGcdElseSplit taken from [START_REF] Maza | Calculs de Pgcd au-dessus des Tours d'Extensions Simples et Résolution des Systèmes d' Équations Algébriques[END_REF] and based on a splitting process à la D5 [START_REF] Dora | About a new method for computing in algebraic number fields[END_REF]. See [START_REF] Lazard | Solving Zero-dimensional Algebraic Systems[END_REF][START_REF] Aubry | Ensembles triangulaires de polynômes et résolution de systèmes algébriques[END_REF] too. The implementation is far more complicated and optimizations are detailed in section 7.

Denote p and q two nonconstant polynomials of the ring G[X1, . . . , Xn] such that q is monic and has leader X k .

The polynomial quo(p, q, X k ) (resp. rem(p, q, X k ) ) denotes the quotient (resp remainder) of the Euclidean division of p by q. The division does not raise any problem since q is monic.

If T = {p1, . . . , p k } is a normalized triangular set of the ring G[X1, . . . , X k ], then, in the algorithms, p mod T denotes rem(. . . rem(p, p k , X k ) . . . , p1, X1), (i.e. the remainder of the Euclidian division of p by all polynomials of T ). Remark: T can be considered as a Gröbner basis and p mod T is equal to the normal form [3, page 199] of p by T .

invert(p, {p1, . . . , p k }) if p ∈ G return 1/p else let 1 ≤ j ≤ k such that ld p = ld pj (g, u, v) := ExtEuclid(p, pj, {p1, . . . , pj-1}) if g = 1 then return u else
return to top level the triple (j, g, quo(pj, g, X k ) mod {p1, . . . , pj-1}) fi fi

The algorithm invert is based on the algorithm ExtEuclid we describe below.

Specification of ExtEuclid

T = {p1, . . . , p k-1 } is a normalized triangular set of the ring G[X1, . . . , X k-1 ] and p and q are two polynomials of G[X1, . . . , X k ] with ld(p) = ld(q) = X k and q is monic. ExtEuclid tries to computes a triple (g, u, v) of three polynomials of G[X1, . . . , X k ] such that:

• up + vq = g (mod (T )) (Bézout [4] identity)

• lcoeff(g, X k ) = 1, i.e. g is either the constant polynomial 1 or a monic polynomial with leader X k .

• g divides both p and q modulo the ideal (T )

Algorithmic scheme of ExtEuclid ExtEuclid(p, q, T ) if q = 0 then ıp :=invert(ip, T ) return (ıp p mod T, ıp, 0) else a := quo(p, q, X k , T ) b := rem(p, q, X k , T )

ı b :=invert(i b , T ) (g, u, v) := ExtEuclid(q, ı b b mod T , T ) return (g, ı b v mod T, u -ı b v a mod T ) fi
Remark: in case of splitting, the algorithm ExtEuclid does not compute the expected triple because it is interrupted by the algorithm invert at the line "return to top level".

SATTRIANGULAR

Recall H A ⊂ S and L = ld(A). The algorithm satTriangular builds a finite sequence (Fi) 1≤i≤r of r sets of systems of equations and inequations.

Initially, take F0 = {(A = 0, S = 0)}. We suppose we have built the set Fi. Two cases may arise:

• Denote Fi = {(T1 = 0, S1 = 0), . . . , (Tt = 0, St = 0)}.

If each Ti is a normalized triangular set and each Si is empty, the algorithm then stops and outputs the set {T1, . . . , Tt}.

• Fi contains a system A = 0, S = 0 such that A contains a non monic polynomial or S is not empty. Then transform A = 0, S = 0 with one of the two rules R1 or R2 (defined below). Both rules compute a set F of zero, one or two systems of equations and inequations.

Take Fi+1 = Fi \ {A = 0, S = 0} ∪ F.

R1 : try to make a polynomial monic. If there exists some non monic p k ∈ A s.t. p1, . . . , p k-1 are monic, three cases are possible: Same as R1.1.

• Case R2.2

The inequation s admits an inverse we call q. The property s q = 1 (mod (p1, . . . , p k )) clearly implies that s does not vanish on the zeros of (A) : S ∞ . Thus, by lemma 6, we have (A) : S ∞ = (A) : S ′∞ and Fi+1 satisfies the three invariants.

• Case R2.3 Same as R1.3

This ends the proof of proposition 4

Proposition 5. The algorithms stops.

Proof. To each system Σ = (A = 0, S = 0), associate the sum of the degrees of the elements of A in their leader, the number of the non monic elements of A and the number of elements of S. Denote v(Σ) this positive integer.

When the algorithm rewrites a system Σ into a set F , then v(Σ ′ ) < v(Σ) for each Σ ′ ∈ F.

By [START_REF] König | Theorie der endlichen und unendlichen Graphen[END_REF]Satz 6.6] (i.e. every infinite locally finite tree contains a branch of infinite length), the algorithm stops. Proposition 6. The output {T1, . . . , Tt} of the algorithm satisfies P1, . . . , P5 Proof. This is an immediate consequence of the invariants I1, I2 and I3 satisfied by Fr.

IMPLEMENTATION

The codes of the algorithms invert and ExtEuclid, in section 5, just are algorithmic schemes. Most of the optimizations below are implemented in our private version of the Maple package diffalg. They were suggested to us by Marc Moreno Maza, who is very familiar with triangularization algorithms.

• One can compute pairs instead of triples in ExtEuclid following [17, Volume 2, page 325]: using the final pair and the Bézout identity, one can recover the missing coefficient by a mere Euclidean division modulo a triangular set.

• For efficiency reasons, one can implement different versions of ExtEuclid, one which computes a full Bézout identity, one which only computes one of the Bézout coefficients (for computing an inverse), and one which only computes the gcd (for testing invertibility). For example, Marc Moreno Maza observed that the last division of the former item can be very expensive.

• The classical scheme of the extended Euclidean algorithm can be performed with the subresultants algorithm which controls the growth of the coefficients of the intermediate remainders. The recent Ducos [START_REF] Ducos | Optimizations of the subresultant algorithm[END_REF] optimization of the subresultants algorithm is worth being implemented. It replaces some pseudoremainders of the traditional subresultant algorithm by a few more computations which involve smaller datas. The alternative [START_REF] Lombardi | New structure theorem for subresultants[END_REF] algorithm could be used as well.

• The rule R2.3 can be specialized when j = k. In this case, the system Ag = 0, S = 0 is inconsistent. The same optimization applies for R1.3 when j = k -1.

• It is interesting to invert separants of the equations as soon as possible. Once the separants of all considered equations are inverted, the ideal generated by these equations is radical (Lazard's lemma). This gives another optimization of the rule R2.3 when j = k. In this case, s is necessarily invertible modulo (A h ) and can be removed from S.

• A good strategy consists in applying R2 as soon as possible in the case both R1 and R2 apply. Indeed, this may split the systems into smaller systems which are easier to handle.

Experiments

Some other algorithms which have the same specifications as regCharacteristic are implemented in Maple [START_REF] Boulier | Computing representations for radicals of finitely generated differential ideals[END_REF] and [START_REF] Hubert | Factorization free decomposition algorithms in differential algebra[END_REF] but are superseded by ours: the first step of the other algorithms consists in computing a Gröbner basis of (A) : S ∞ . This is certainly not the best way since the Buchberger algorithm does not take into account the fact that A is already triangular and explicits the inverses of the inequations, which are not needed.

The advantage of our method only appears for differential systems which make Rosenfeld-Gröbner spend a lot of time on the purely algebraic treatment of regular differential systems.

Let's try the following system which has no physical significance but makes Rosenfeld-Gröbner spend most of the time in the purely algebraic part (second step):

v uxx + u 2 xx + ux, uyy + uy for the orderly ranking

• • • > uxx > uxy > uyy > vxx > vxy > vyy > ux > uy > vx > vy > u > v.
We perform our comparisons in Maple 5.5 over a Sun Ultra 5 at 333Mhz with 128Mb memory.

Our private version of Rosenfeld-Gröbner produces six regular differential ideals in 75 seconds, including 62 seconds for regCharacteristic. One of the regular differential systems A = 0, S = 0 to deal with is quite large. We do not give its characteristic presentation.

The other implementations of Rosenfeld-Gröbner (the one in Maple 5.5 and its variant by Hubert) cannot carry this example out: Maple 5.5 does not succeed in computing the Gröbner basis of the localized ideal S -1 (A). To any differential polynomial q we may associate a fraction NF(q, C) = a b satisfying 1. b q = a mod a, 2. NF(q, C) is a canonical form of the equivalence class of q in R/a (the fraction only depends on the ranking, the differential ideal a and the equivalence class of q). In particular q = q ′ mod a ⇒ NF(q, C) = NF(q ′ , C). The following proposition is straightforward but very important for it permits us to perform easily linear algebra over K in R/a. Proposition 7. The set {NF(q, C) | q ∈ R} forms a vector space over K.

CANONICAL FORMS

We now show how to compute normal forms. First consider the case q ∈ K[L, N ] for which things are easy since C is strongly normalized. Let q = prem(q, C). For some power product h of initials of C we have h q = q mod (C). We define NF(q, C) as the fraction a/b obtained by making q/h irreducible. The differential polynomial q is reduced w.r.t. C. Since C is strongly normalized, h ∈ K[N ].

Lemma 7. b q = a mod a.

Proof. We have h q = q mod a. Factors of h are not zero divisors in R/a thus we may factor them out from above relation. Lemma 8. NF(q, C) = NF(q ′ , C) for any q ′ ∈ K[L, N ] such that q = q ′ mod a. We now consider the general case. The problem is due to the separants which do not belong to K[N ]. It is overruled by inverting them using the extended Euclidean algorithm.

Let q be a differential polynomial. Let q = partial rem(q, C). For some power product h of separants of elements of C we have h q = q mod a.

Using the invert algorithm, we may compute10 the inverse of h modulo (C) : H ∞ C in K(N ) [L]. Multiplying by some element of K[N ] to clear denominators we find a differential polynomial h such that h h = g mod a and g ∈ K[N ]. Therefore g q = h q mod a. Let NF(h q, C) = a/b. We define NF(q, C) to be the fraction a/b obtained by making the fraction a/g b irreducible.

Proofs of the general case are easy variants of the former ones.

Canonicity properties are consequences of theorem 3 and lemma 8.

CONCLUSION

The algorithm presented in this paper is a first step in merging the efficient solvers of polynomial equations implemented by [START_REF] Maza | Calculs de Pgcd au-dessus des Tours d'Extensions Simples et Résolution des Systèmes d' Équations Algébriques[END_REF][START_REF] Aubry | Ensembles triangulaires de polynômes et résolution de systèmes algébriques[END_REF] and the Rosenfeld-Gröbner algorithm. Merging these solvers, we expect to obtain an efficient solver for differential polynomial equations in which there would not be any distinction between differential and algebraic parts. In particular, such a new solver would be able to handle the purely algebraic subproblems which also arise in the differential part of the current Rosenfeld-Gröbner. It would be very interesting for systems of DAE and PDAE.

1 .

 1 the regular algebraic ideal (A) : S ∞ is radical (Lazard's lemma) ; 2. if b denotes a prime ideal minimal over (A) : S ∞ then the set N furnishes a transcendence basis of the field of fractions of R/b over K (Lazard's lemma) ; 3. we have [A] : S ∞ ∩ K[L, N ] = (A) : S ∞ (Rosenfeld's lemma) ; 4. the regular differential ideal [A] : S ∞ is radical (lifting of Lazard's lemma) ; 5. if (A) : S ∞ has t minimal primes bi then [A] : S ∞ has t minimal differential primes pi defined by (lifting of Lazard's lemma) pi ∩ K[L, N ] = bi. Actually pi = {p ∈ R | partial rem(p, A) ∈ bi}.

Theorem 2 . 1 .

 21 Let C be an autoreduced set of differential polynomials and a = [C] : H ∞ C . The following conditions are equivalent. C is a characteristic set of a, 2. p ∈ a ⇔ full rem(p, C) = 0, 3. C = 0, HC = 0 is a regular differential system such that C is regular 6 in the sense of[START_REF] Aubry | On the theories of triangular sets[END_REF] and squarefree 7 in the sense of[1, def 4.5.11].Proof. 1 ⇒ 2. If p ∈ a then full rem(p, C) ∈ a, is reduced w.r.t. C, thus is zero. If full rem(p, C) = 0 then, for some power product h of elements of HC we have h p ∈ [C] ⊂ a whence p ∈ a.2 ⇒ 3. The differential system C = 0, HC = 0 satisfies conditions C1 and C2. The set C reduces to zero all ∆(p, p ′ ) such that {p, p ′ } ∈ critical pairs(C) (for C reduces a to zero) thus C = 0, HC = 0 satisfies C3 and is a regular differential system. C reduces to zero (C) : I ∞ C ⊂ a whence is regular in the sense of [2] by [2, Theorem 6.1]. C reduces to zero (C) : H ∞ C ⊂ a without differentiating any element of C whence (C):H ∞ C ⊂ (C):I ∞ C . Since the converse inclusion obviously holds too we have (C) : I ∞ C = (C) : H ∞ C thus C is squarefree by theorem 1 (1). 3 ⇒ 1. Let f be a differential polynomial reduced w.r.t. C. We must prove f = 0. We have f ∈ K[L, N ]. Since C = 0, HC = 0 is a regular differential system, theorem 1 (3) applies and f ∈ (C) : H ∞ C . The set C is regular in the sense of [2] thus C is a characteristic set of (C) : I ∞ C by [2, Theorem 6.1]. Since C is squarefree, (C) : I ∞ C is radical and (C) : I ∞ C = (C) : H ∞ C by [15, Proposition 3.3]. Therefore C is a characteristic set of (C) : H ∞ C and f = 0.

Theorem 3 .

 3 A characteristic presentation C is a canonical representative of the regular differential [C]:H ∞ C (it only depends on the ideal and on the ranking). Proof. Let C and C ′ be two characteristic presentations of [C] : H ∞ C . Both sets have the same rank for they are characteristic sets of the same ideal. Consider any f ∈ C and f ′ ∈ C ′ having the same rank. Denote i and i ′ their initials. The polynomial i ′ fi f ′ ∈ (C) : H ∞ C . Since C and C ′ are autoreduced and strongly normalized, this polynomial is reduced w.r.t. both C and C ′ whence is zero. Since f and f ′ are primitive over K[N ] we have f = f ′ thus C = C ′ . In practice, in order to have canonicity properties, we impose also the coefficients in K of the elements of C to be normalized. In the case K = Q(Y ) is a pure transcendental field extension of the field of the rational numbers (Y is an alphabet of indeterminates) then the elements of C are polynomials in Z[Y ∪N ∪L] primitive over the ring Z[Y ∪N ].

3 .

 3 Transform each triangular set Ti of G[L] as a triangular set Ci of K[X] by replacing every polynomial p/s occuring in the Ti systems by the primitive part of p over K[N ]. The obtained systems C1, . . . , Ct are the characteristic presentations we are looking for.

13 ,Lemma 3 .

 133 Proposition 2.2, (a)] and maps (A) : S ∞ (resp. (Ti)) to (A) : S ∞ (resp. (Ci) : H ∞ C i ) by P4, [13, Proposition 2.2] and the fact that the nonzero elements of K[N ] belong to none of the minimal primes of (A) : S ∞ (resp. (Ci) : H ∞ C i ) by theorem 1 (2) (resp. P5 and theorem 1 (2)). Lemma 2. A prime ideal b is minimal over (A) : S ∞ iff b is minimal over some (Ci) : H ∞ C i . Proof. By P5 and theorem 1 (2), the prime ideals which are minimal over (A) : S ∞ and the ones which are minimal over the (Ci) : H ∞ C i all have the same dimension. If b ⊂ b ′ are two prime ideals having the same dimension then b = b ′ . The lemma follows now from lemma 1. Each system Ci = 0, HC i = 0 is a regular differential system. Proof. It suffices to prove if {p, p ′ } ∈ critical pairs(Ci) then full rem(∆(p, p ′ ), Ci) ∈ b where b is any prime ideal minimal over (Ci) : H ∞ C i . By lemma 2 and theorem 1 (5) b is the intersection with K[X] of some differential prime ideal p minimal over [A] : S ∞ . Since p, p ′ , Ci belong to p we have full rem(∆(p, p ′ ), Ci) ∈ b.

Lemma 4 . 0 for 9

 409 Let 1 ≤ i ≤ t be an index. The set Ci is a characteristic set, in the sense of Ritt, of the ideal (Ci):H ∞ C i Proof. Let us prove that p ∈ (Ci) : H ∞ C i iff prem(p, Ci) = 0 for any p ∈ K[X]. The implication from right to left is clear. The converse one comes from the following: first, p ∈ (Ci) : H ∞ C i iff φp ∈ (Ti) by P4 ; second, prem(p, Ci) = 0 iff prem(φp, Ti) = 0 ; third prem(φp, Ti) = 0 iff φp * ----→ T i the elements of Ti are monic ; last, Ti is a Gröbner basis of (Ti) [3, Lemma 5.66 (Buchberger's first criterion)] thus reduces to zero all the elements of (Ti) [3, Proposition 5.38]. Proposition 1. (condition

Let

  C be the characteristic presentation of the differential ideal a = [C] : H ∞ C in the differential ring R = K{U }. De-note L the set of leaders of C and N the set of the other derivatives occuring in C.

3 .

 3 a ∈ K[L, N ] is reduced w.r.t. C and b ∈ K[N ] (in particular b does not divide zero in R/a).

  Proof. Denote NF(q ′ , C) = a ′ /b ′ . The differential polynomial a b ′a ′ b ∈ a, is reduced w.r.t. C for b, b ′ ∈ K[N ]and a, a ′ are reduced w.r.t. C. It is zero for C is a characteristic set of a. Since the fractions are reduced a = a ′ and b = b ′ .

As[START_REF] Boulier | Computing representations for radicals of finitely generated differential ideals[END_REF][START_REF] Hubert | Factorization free decomposition algorithms in differential algebra[END_REF] already do.

This is quite anecdotic for the algorithmic consequences of the Gianni and Kalkbrener theorem do not apply in our case.

This definition was introduced for the first time in[8, sec-

tion 4], under the name "pair".

[START_REF] Boulier | Computing representations for radicals of finitely generated differential ideals[END_REF] Definition 4.4].

If IA ⊂ S and full rem(∆(p, p ′ ), A) = 0 for every critical pair {p, p ′ } ∈ critical pairs(A) then C3 is satisfied.

This denotes the reduction in the sense of the Gröbner basis theory[3, page 199].

Observe that, though h is invertible mod (C) : H ∞ C the invert algorithm may fail to compute the inverse. In that case however a splitting of the ideal is exhibited and computations can be restarted over each of its branches. One avoids this problem by precomputing the inverses of the separants of the elements of C

We would like to thank Marc Moreno Maza for many explanations about his algorithms and triangular sets. We would like to thank both Marc Moreno Maza and François Ollivier for their comments on a preliminary version of this paper.

R1. [START_REF] Aubry | Ensembles triangulaires de polynômes et résolution de systèmes algébriques[END_REF] The initial of p k is zero modulo (p1, . . . , p k-1 ).

Take F = ∅.

R1.2 invert finds the inverse q of ip k modulo the ideal (p1, . . . , p k-1 ). Take F = {A ′ = 0, S = 0} where A ′ = A \ {p k } ∪ {p k } where p k = q p k mod {p1, . . . , p k-1 }.

R1.3 invert does not find the inverse of ip k modulo (p1, . . . , p k-1 ), but a triple (j, g, h) such that 1 ≤ j < k and pj = g h (mod (p1, . . . , pj-1)). Take F = {(Ag = 0, S = 0), (A h = 0, S = 0)} where Ag = A \ {pj} ∪ {g} and A h = A \ {pj} ∪ {h}.

R2 : try to get rid of an inequation. If there is some s ∈ S such that ld s = ld p k and p1, . . . , p k are monic then three cases are possible: R2.1 s is zero modulo (p1, . . . , p k ). Take F = ∅.

R2.2 invert finds the inverse of s modulo (p1, . . . , p k ). Take F = {A = 0, S ′ = 0} where S ′ = S \ {s}.

R2.3 invert does not find the inverse of s modulo the ideal (p1, . . . , p k ) but a triple (j, g, h) such that 1 ≤ j ≤ k and pj = g h (mod (p1, . . . , pj-1)). Take F = {(Ag = 0, S = 0), (A h = 0, S = 0)} where Ag = A\{pj}∪{g} and A h = A\{pj}∪{h}.

Proof

Proof. (sketched) The proof follows from the two following points: first, every zero of I is an irreducible component of the algebraic variety of I (since I is zerodimensional); second, the algebraic variety of I : S ∞ is the union of all the irreducible components of I which do not annihilate any element of S.

Corollary 2. For each 1 ≤ l ≤ t, the ideal (A l ) : S ∞ l has dimension zero and is radical.

Proof. This is a corollary to invariant I2 and Lazard's lemma (theorem 1 (1,2)).

Proof. of proposition 4

The proof is an induction on i. The basis of the induction is clear. Suppose Fi verifies the invariants. We prove that Fi+1 does too. We distinguish several cases, corresponding to the way Fi+1 is built. • Case R1. [START_REF] Aubry | On the theories of triangular sets[END_REF] We have (A) = (A ′ ) since q p k ∈ (A) and p = ip k q p k (mod (p1, . . . , p k-1 )). Therefore, by saturating, (A) : S ∞ = (A ′ ) : S ∞ . Thus, Fi+1 satisfies I1 and I3.

It remains to prove I2. We have (A) : S ∞ = (A) : (S ∪ HA) ∞ (invariant I2 over Fi). Lemma 6 implies that the polynomials of HA do not vanish on the zeros of (A) : S ∞ . Since S A ′ = SA \ {sp k } ∪ {sp k } and sp k = q sp k (mod (p1, . . . , p k-1 )), the polynomials of H A ′ do not vanish either on the zeros of (A) : S ∞ . Thus, by lemma 6, (A) :

We have pj = g h (mod (p1, . . . , pj-1)). Since pj, g and h have the same leader, sp j = sgh + gs h (mod (p1, . . . , pj-1)).

Denote I = (A) : S ∞ and Ig = (Ag) : S ∞ and

We have V (I) ⊃ V (Ig) ∪ V (I h ) for I ⊂ Ig and I ⊂ I h .

The converse inclusion. Let z ∈ V (I).

g(z)h(z) = 0 (4)

Thus, if g(z) = 0 then sg(z) = 0 and h(z) = 0 ; if g(z) = 0, then h(z) = 0 and s h (z) = 0. This implies z is either a zero of Ig or a zero of I h , and it can't be zero of both ideals. Thus 2 and 3 are proved and Fi+1 satisfies I3.

Proof of I2 :

By invariant I2 on Fi and lemma 6, the polynomials of HA do not vanish on V (I). By relation 2, they do not vanish on V (Ig). The polynomial sg does not vanish on V (Ig) (consequence of relation 5). We have HA g = HA \ {sp j } ∪ {sg}. Therefore, by lemma 6, (Ag) : S ∞ = (Ag) : (S ∪ HA g ) ∞ .

The same proof holds for I h . This ends the proof of I2.

By the theorem of zeros, √ I = p Ig ∩ √ I h . By invariant I2 and theorem 1(1), all these ideals are radical. Thus Fi+1 satisfies I1.

• Case R2.1