Christophe Rehmet

Frédéric Fabry

Vandad-Julien Rohani

François Cauneau

Laurent Fulcheri

Frédéric Fabry

Mhd

Frederic Fabry

MHD modeling of a 3-phase plasma torch with parallel electrodes, comparison with coplanar electrodes

Keywords: electric arc motion, 3-phase AC plasma torch, MHD modeling

teaching and research institutions in France or abroad, or from public or private research centers.

MHD modeling of a 3-phase plasma torch with parallel electrodes, comparison with coplanar electrodes.

Introduction

3-phase AC power supplies with graphite electrodes have proven their durability and reliability for many years in steel industry for few MW [START_REF] Bakken | Numerical modeling of electric arcs[END_REF]. Plasma technology based on this system, could overcome some problems of reliability and equipment/operating costs of DC arc plasma torches. Indeed, short lifetime of the electrodes [START_REF] Matveev | Guest Editorial Classification of Plasma Systems for Plasma-Assisted Combustion[END_REF] and sensitive electronics due to the rectifier part of the electrical signal give some limits to DC technology.

A solution developed for about 20 years at MINES-ParisTech is a 3-phase AC plasma technology with consumable graphite electrodes. By contrast with 3-phase arcs furnace, where the arcs transit via molten metal as star neutral point, there is not a neutral point which is ground connected in this developed system. Consequently, this technology is based on the generation of rotating arcs entrained in a gas flow generated at the tip of three electrodes. Thus, the 3-phase system is assimilated to a triangle configuration where generally only one arc can exist at one and two arcs can coexist only during a very short duration. Each electrode plays then alternatively two times by period the role of cathode and anode. One period is decomposed into the generation by rotation of 6 arcs. According to Fabry et al [START_REF] Fabry | Waste Gasification by Thermal Plasma: A Review[END_REF], this plasma system could be a solution to for the future development at industrial scale of the Waste-to-Energy gasification processes based on thermal plasma. However, very few studies have been dedicated to 3-phase AC plasma torches. Indeed, most of the industrial applications are based on DC plasma torches [START_REF] Gleizes | Thermal plasma modeling[END_REF][START_REF] Jp Trelles | Arc plasma torch modeling[END_REF].

The 3-Phase plasma system was initially composed of three converging graphite electrodes located on three generating lines of a 15° revolution cone [START_REF] Ravary | Influence of the Electromagnetic Forces on Momentum and Heat Transfer in a 3-Phase ac Plasma Reactor[END_REF]. Recently, a new configuration with three radial electrodes being located in a plan perpendicular to the reactor axis has been developed. The arc behavior has been simulated by a MagnetoHydroDynamic (MHD) model [START_REF] Rehmet | 3D unsteady state MHD modeling of a 3-Phase AC hot graphite electrodes plasma torch[END_REF]. This model automatically simulates the arcs ignition, extinction and produces a particular arc motion between the three electrode tips. Further, the analysis of the 3-phase plasma torch has been carried out using a 100 000 fps high speed video camera and oscilloscope [START_REF] Rehmet | High speed video camera and electrical signal analyses of arcs behavior in a 3-Phase AC arc plasma torch[END_REF]. Results show a fair correlation between MHD modelling and experiments, both regarding the arc behaviour and electrical waveforms [START_REF] Rehmet | Unsteady state analysis of the arc behaviour in a 3-Phase AC arc plasma torch. A comparison between MHD modelling and experimental results[END_REF]. In this paper, the results of MHD modelling with coplanar electrodes configuration are compared with parallel electrodes configuration and experimental images obtained with a high speed camera.

Description of the 3-Phase AC Plasma Torch

The parallel electrodes configuration grid mesh was created using the software Salome 6©. The computational domain refers to a slightly different geometry to the actual converging graphite electrodes located on three generating lines located on a 15° revolution cone. Due to meshing issue, we consider three parallel electrodes (Fig. 1). Nevertheless, we can assume that this 15° has not a significant influence on the 3-phase arc discharge behavior compared to parallel electrodes. A functional fluid area of the arcs evolution has been considered. A non-uniform hexahedral mesh with approximately 800 000 nodes is used and is refined at the vicinity of the electrodes.

Assumptions and equations

The study was carried out for an arc discharge operating with nitrogen. The data of the properties have been taken from [START_REF] Pateyron | T&T Winner, the chemistry of on-line transport properties in interval of 300 K to 20.000 K[END_REF]. The 3-D model of the three phase system is based on the following main assumptions that are commonly used.

(1) The plasma is considered as a single continuous fluid (nitrogen) and at Local Thermodynamic Equilibrium (LTE).

(2) The gas is uncompressible.

(3) The gas flow is laminar and time dependent. (4) Inductive currents are neglected.

(5) Gravitational effects and radiation are not taken into account. [START_REF] Ravary | Influence of the Electromagnetic Forces on Momentum and Heat Transfer in a 3-Phase ac Plasma Reactor[END_REF] The voltage drop in the sheath is not considered in the model Given the assumptions listed above, fluid equations can be written by normal conservation equations including the Lorentz force and the Joule effect, assuming the approximation of Ohm's law in its simplified form and the overall electrical neutrality [START_REF] Rehmet | 3D unsteady state MHD modeling of a 3-Phase AC hot graphite electrodes plasma torch[END_REF].

Boundary Conditions

Boundary conditions are detailed and exposed on a cross section view on the Fig. 2 and Table 1. The boundary conditions could be separated into three different categories: (1) hydrodynamic conditions, (2) thermal conditions and (3) electromagnetic conditions. Regarding the hydrodynamic conditions, the central plasma gas injection is located in the AB, EF zones. The velocity is 0.1 m s -1 , giving a flow rate of 3.8 Nm 3 h -1 .

Current density and temperature assumptions on the boundary condition of the electrode are not suitable to well simulate the complex arc behavior in such system. The electrodes are then incorporate into the computational domain (zone BCDE). The objective is to characterize the magnetic interactions between the 3 graphite electrodes and arcs.

Regarding the electromagnetic conditions, a sinusoidal electrical potential of 680 Hz frequency is imposed on each electrode. This electrical potential is phase shifted by 2π/3 between the three electrodes and is adjusted by the model using a variable called ϕr, in order to control the measured current. Table 1 Boundary conditions applied to a section of the computational domain, with f, t ref and it the frequency, the times step and the iteration number respectively.

Implementation of the Numerical Model on Code Saturne

The software Code_Saturne v. 2.0, electric arc module, was used for the 3D transient model of the three phase plasma torch. A 2 μs time-step has been used to give a reasonable computation time. The study is partitioned on eight processors (Intel Xeon 2.66 GHz). One iteration takes almost 1 min 20 s. To simulate one period (1.5 ms), we need about 750 iterations, i.e. 16 h of calculation. To ignite the first arcs at the first time-steps, a hot zone at 6,000 K is artificially implemented between the three electrodes and the peak current target is set during half of a period (0.75 ms) from 5 to 400 A.

As previously mentioned, the 3-phase system can be assimilated to a triangle configuration [START_REF] Ravary | Influence of the Electromagnetic Forces on Momentum and Heat Transfer in a 3-Phase ac Plasma Reactor[END_REF][START_REF] Weidong | Characterization of a 3-phase a.c. free burning arc plasma[END_REF] where two arcs can coexist just during a short time. We suppose that each new arc disturbs the previous by a repulsive magnetic effect, or the joule effect becomes insufficient to maintain two arcs in same time. Assuming this hypothesis, the theoretical current flowing in the arcs follows the maximum of the phase-to-phase voltage. We impose a constant value for the sum of the currents measured on the three planes between the 3 electrodes. This method shows results close to the experimental form. In the case of a phase delay or a lacking arc, this method gives also more degree of freedom to the model.

Results

In parallel and coplanar electrode configuration, MHD modeling has been successfully implemented and the model automatically simulates the arcs ignition, extinction and motion between the three electrodes tips.

Results show a different arc behavior in both case controlled by hydrodynamic electrode jet forces in coplanar electrode configuration instead of repulsive magnetic forces in parallel electrode configuration, as shown in Fig. 2. This particular arc behavior has been already suggested and estimated by Ravary et al [START_REF] Ravary | Influence of the Electromagnetic Forces on Momentum and Heat Transfer in a 3-Phase ac Plasma Reactor[END_REF]. He demonstrated that, in such geometry, the interaction of the 3 line currents flowing in the electrodes and within the arcs induced strong rotating magnetic repulsive forces on the arc roots leading to a global centrifugal arc motion. Besides, this global arc motion in the centrifugal direction was also suggested and estimated in other papers [START_REF] Bermúdez | Numerical computation of the electromagnetic field in the electrodes of a three-phase arc furnace[END_REF]. Indeed, the 3-phase AC plasma torch has some similarities with electric arc furnaces. Moreover, from a parametric study, it appears that the imposed current will change the arc discharge behavior. By decreasing the current down to 50 A RMS, the arc stay within the inter electrode gap. Conversely, by increasing the current, arc roots move along the surface of the electrodes and finish on the extremity of the electrode tip (Fig. 3). We can then notice that the Lorentz forces produced within the electrode have a strong influence on the arc behavior. A previous model, without the electrode incorporate into the computational domain, has shown that the arcs stay within the inter electrode gap like in coplanar electrodes configuration.

Consequently, by increasing the current above 150 A RMS, the arc motion increase strongly toward the outflow of the computational domain. The radial distance of this domain has also been increased in order to avoid that the arcs impact the torch walls. With the arc elongation, the mass flow emanating from the electrodes jets do not directly contribute to the new arcs ignition. The heat transmitted by this induced mass flow is directed outward to the inter electrode gap. The path of the arc produce a thermal channel, in which arc discharge operate. The increasing of the current raises the Lorentz forces. These forces give a motion to the thermal channel. By increasing the frequency, the motion of this thermal channel decrease and the 3-Phase arc discharge is more stable than previously. At 50 Hz, the 3-Phase system becomes unbalanced and only one arc remains after around 10 ms. Indeed, the temperature field in the others hot channel between the other others electrodes decreases below 5000 K. Arc discharge within this channel is avoided since the electrical conductivity is not enough. Consequently, the evolution of the hot channel is a key feature on the 3-Phase arc discharge stability. This evolution seems to be related to the imposed frequency, which influence the arc motion and then the 3-Phase arc discharge stability.

Regarding the electrical values, the power significantly increases with parallel electrode configuration compared to coplanar configuration. Indeed, heat transmitted by electrode jets in coplanar electrode configuration within the inter electrode gap stabilize the line voltage around 35 V. In fact, the increase of the current stabilizes the arcs within inter-electrode gap. Conversely, in parallel electrode configuration, the arcs are not stabilized within the inter electrode gap. The increase of the current raises the hot channel motion and the velocity of the arc column is faster than the plasma gas velocity injected. Arc elongation induces line voltage up to 100 V RMS. The power increases then significantly around 70-50 KW whereas it was around 8-20 KW in coplanar configuration.

Conclusion

A 3D MHD modeling of a 3-phase AC plasma torch with parallel graphite electrodes has been developed. Significant information that can hardly be achieved experimentally has been obtained with this modeling. Results show a different arc behavior than the previous developed model with coplanar electrode configuration.

In this new MHD model, the arcs are mainly controlled by magnetic forces instead of hydrodynamic electrode jet forces. These forces give a particular arc behavior which has already been suggested and estimated by different studies. Results obtained with this model show that the arc motion is strongly influenced by the current, plasma gas flow rate and frequency. Frequency is then a key feature in electrical arc motion since it modifies the 3phase discharge stability. This point, revealed by this study is interesting for mainly arc applications, as arc furnace systems and also other 3 phase plasma torches. Globally, arcs generated between three graphite electrodes induce a specific heat and mass transfer which tend to homogenize the plasma flow temperature. In other ISPC paper, results obtained with this model are compared to experimental results obtained with a high speed camera. Results show a fair correlation between MHD modeling and experiments, both regarding the arc motion and the dissipated power value.

Fig. 1

 1 Fig.1 Computational domain used in this study composed of three parallel electrodes.

Fig. 2

 2 Fig.2 Grid mesh exposed on a cross section showing one electrode.

Fig. 3

 3 Fig.3 Arc behavior in coplanar electrode configuration (on the top) with temperature field (left) and velocity field with particule tracking (right) in a cross section located in a plan containing the 3 electrode axes. At the bottom, arc behavior in parallel electrode configuration correlated to theoretical arc


Christophe Rehmet

Frederic Fabry

Vandad Rohani

François Cauneau

Laurent Fulcheri

Experimental analysis of arcs behavior in a 3-Phase AC arc plasma torch with hot graphite electrodes: comparison with MHD modeling results

Keywords: electric arc motion, 3-phase AC plasma torch, MHD modeling

Arcs behavior in 3-Phase plasma systems remains poorly studied. In this paper, arc behavior observed with a high speed camera has been compared to some typical results obtained by a MHD modeling. Numerical results highlights the origin of some typical arc shape obtained experimentally.

Introduction

For about 20 years, a 3-phase AC plasma technology with high temperature consumable graphite electrodes is under development at MINES-ParisTech [START_REF] Fabry | Waste Gasification by Thermal Plasma: A Review[END_REF]. As displayed in Fig. 1, this system was initially composed of three converging graphite electrodes located on three generating lines of a 15° revolution cone [START_REF] Ravary | Influence of the Electromagnetic Forces on Momentum and Heat Transfer in a 3-Phase ac Plasma Reactor[END_REF]. Recently, a new configuration with three radial electrodes being located in a plan perpendicular to the reactor axis has been developed. The arc behaviour in such system has been simulated by a MHD model [START_REF] Rehmet | 3D unsteady state MHD modeling of a 3-Phase AC hot graphite electrodes plasma torch[END_REF]. As illustrated in Fig. 2, the analysis of the 3-phase plasma torch using a 100 000 fps high speed video camera [START_REF] Rehmet | High speed video camera and electrical signal analyses of arcs behavior in a 3-Phase AC arc plasma torch[END_REF] has shown a fair correlation between MHD modelling and experiments, both regarding the arc behaviour and the current waveform [START_REF] Rehmet | Unsteady state analysis of the arc behaviour in a 3-Phase AC arc plasma torch. A comparison between MHD modelling and experimental results[END_REF].

This MHD model was faintly modified and adapted for the simulation of the 15° angular electrode configuration. Due to meshing issue, we have simulated this 3-Phase system with parallel electrodes configuration. We assume that this 15° angle has not a significant influence on the arc behaviour compared to the parallel electrode configuration. Besides, the graphite electrodes are now considering into the computational domain as mentioned in other IPSC paper. In this paper, the numerical results are compared with experimental images obtained with a high speed camera. First, the main characteristics of the experimental setup are described. Then, the comparison between experimental and some numerical results are presented. Fig. 2 Comparison of the temperature fields in a cross section located in a plan containing the 3 electrode axes (left) with the arc video sequence (right).

3-Phase AC Plasma Torch Set Up

The studied 3-phase plasma system is composed of three graphite electrodes (Ø 17 mm) located on three generating lines located on a 15° revolution cone, as represented in Fig. 1. This system has no ground connected neutral point and is assimilated to a triangle configuration. Generally, only one arc can exist at time and arc discharge follows the maximal electric gap potential. Each electrode plays alternatively two times the role of cathode and anode by period and one period is decomposed into the generation by rotation of 6 arcs. The frequency was set to 680 Hz by a 263 KW 3-Phase converter (rectifier-chopper-inverter). Video camera analysis was carried out using a high-speed video camera Olympus FS with a 50,000 fps shooting rate, 20 μs exposure time and 180x132 pixels on the raw images recorded with the camera. This camera is located at around 2 m optical distance from the electrodes tips. Experiments were realized for 200A phase current, 3.82 Nm 3 .h -1 nitrogen plasma gas flow rate. Various optical filters were used in order to adjust the image brightness.

Comparison between experimental results and MHD modelling

Video sequences obtained with the high speed camera show that the 3-phase discharge is mainly erratic. The chronological 6 arcs formation in each period is almost verified and the 3-Phase system is sometimes unbalanced. This arc discharge behaviour could not be perfectly reproduced by the simulation since the arc discharge is mainly chaotic. However, this behaviour is partially simulated with the model and it is then analysed with the MHD model in order to improve the understanding of such discharge. Besides, arc motions are observed by a bottom view of the inter electrode gap. The numerical results show the same arc shape with an entire 3D view.

Experimental arcs can have different geometries that can be classified into numerous shapes. Five main classical shapes already reported in the literature for different plasma systems: I, V, U, W or S [START_REF] Weidong | Characterization of a 3-phase a.c. free burning arc plasma[END_REF][START_REF] Schnick | Three-dimensional modelling of arc behaviour and gas shield quality in tandem gas-metal arc welding using anti-phase pulse synchronization[END_REF][START_REF] Colombo | 3D static and time-dependent modelling of a dc transferred arc twin torch system[END_REF][START_REF] Tang | Fang Three-dimensional modelling of a DC arc plasma in a twin-torch system[END_REF] are observed. It is commonly accepted that these shapes are mainly influenced by the electrode jet velocity and magnetic repulsive forces induced by the electric current. Consequently, the arc root position on the electrode has a significant influence on these particularly arc shapes.

The comparison of the results obtained with the high speed camera and the results obtained with MHD model is exposed in Fig. 3. This figure shows the decomposition, in a cross projection view under the inter electrode gap of the arc shapes. The images have been correlated with the numerical results. Numerical results were obtained with a sinusoidal voltage waveform and a 400 A imposed peak current, which gives almost 200 A RMS current value. In the experimental setup, a 150-200 A RMS current is used and the voltage waveform has a non sinusoidal waveform (3-Phase chopper-inverter output waveform). However, it has been demonstrated in a previous article that the voltage waveform shape has a negligible influence on the current waveform [START_REF] Rehmet | Unsteady state analysis of the arc behaviour in a 3-Phase AC arc plasma torch. A comparison between MHD modelling and experimental results[END_REF].

To initialize the simulation, a hot zone at 6,000 K is artificially implemented between the three electrodes and the current target is set during 0.75 ms from 5 to 400 A. With the increasing of the imposed current, the hot zone moves toward the torch outlet. Indeed, the hot zone moves faster than the velocity of the surrounding gas flow. The temperature in this hot zone remains above 5,000 K, which is a privileged area for the arc discharge. Once the arc root and the hot zone reach the bottom of the electrodes, the arc curvature increases; partly induced by the magnetic arc repulsion and the electrode jets. Thus give a centrifugal arc motion, as shown in Fig. 3 (t= 1.46 ms). This arc shape is almost reached in 1-1.5 ms, which correspond to the end of the current increment or one period. Arc shape is mainly influenced by the magnetic forces rather than the plasma gas aerodynamics forces.

Later on, on image t = 1.98 ms, one observe that the arcs are characterized by a particular arc shape with two electrode vapour jets. The mass flow emanating from each electrode vapour jet is almost entirely directed along the normal of its associated electrode tip and the arc has almost a U shape. With the arc column motion toward the torch outlet, the phase voltage increases. Therefore an arc column ignites between the two electrode jets which decreases the phase voltage, as shown at t= 2.02 ms and in Fig. 4. Fig. 3 Comparison of the arc shapes, displayed by a 10 000 K isosurface temperature and hot channels displayed by a 5 000 K isosurface temperature, with the arc video sequence by a bottom view of the inter electrode gap.

With the rotation of the electrical potential difference on the electrodes, the phase voltage increases to maintain the arc. Thus, an arc column ignites with the other electrode in a previous arc column where the gradient of the electrical potential is high. As represented in Fig. 4, the arc shaped by a front view, one part of the initial arc column switches off and the arc discharge remains mainly in the other part of the previous hot channel. Nevertheless, we can observe that the arc does not follow the maximal phase to phase voltage. Concerning the arc shapes, the position of the new arc root will change this shape. Different shapes have been observed, as shown in Fig. 3 and Fig. 4. These shapes modify the position of the hot channels which contributes to the erratic arcs discharge.

Later on, a part of the heat contained in the arc column cross an old hot channel and the arcs ignite on this channel. Voltage between arc column and electrode a key parameter on the 3-phase arc discharge behaviour. Besides, the arc motion around the inactive electrode enhanced the 3-Phase arc discharge. Results of the MHD modelling highlight this behaviour. This mechanism is well reproduced by the simulation and occurs three times per period. The arc rotates by passing from a pair of electrodes to another one with a frequency which is twice that of the power supply. Besides, unbalanced 3-phase behaviour was sometimes observed with the high speed video camera. In this case, the arc column does not cross the periphery of the inactive electrode zone. With the MHD modelling, this behaviour was also observed. However, the calculations diverge with this discharge. As shown on image at t= 3.28 ms or on Fig. 5, the arc is directed toward the torch walls. The temperature within the inter electrode gap decreases and the electrical conductivity becomes insufficient to ignite a new arc with the other electrode. This point highlights the key feature of the hot channels motion near to the inter-electrode gap for the 3-Phase arc discharge stability. Indeed, it has been highlighted that decreasing the frequency down to 50 Hz, the arc elongation increases which produces an unbalanced 3phase arc discharge. Fig. 5 Arc shapes when the 3 phase system is unbalanced, arc is represented by a 8 000 K isosurface temperature and hot channels are represented by a 5 000 K isosurface temperature,.

Conclusion

In order to understand the 3 phase arc discharge behaviour with angular electrode configuration (15°), a MHD model has been developed. This model is based on a previous one, developed for coplanar electrode configuration. The model has been modified by considering parallel electrode configuration. Experimental results obtained with a high speed camera show a fair correlation between the experimental and numerical arc behaviours and arc shapes. These results highlighted the strong influence of the arc motions and the hot channels positions on the 3-phase plasma system operation. The heat contained on the hot channels near the inter-electrode gap contributes to ignite new arc roots by increasing the gradient of the electrical potential in the zone near the inactive electrode. Arc ignition is therefore mainly related to the evolution of the hot channels which is related to the arc motion. This arc motion is mostly influenced by the electrode jet forces (Maecker effect) and Lorentz forces. These forces modify the heat exchanged in the system and are mainly influenced by the current, plasma gas and electrodes configuration. In other ISPC paper, the influence of such parameter on the 3 phase discharge behaviour is discussed.

Fig. 1 3

 1 Fig.1 3-phase plasma torch systems developed at PERSEE Mines-ParisTech with coplanar electrode configuration (right) and angular electrode configuration (left).

Fig. 4

 4 Fig.4 Arc shapes obtained with the MHD modeling, displayed on a front view by a 8 000 K isosurface temperature and the hot channels are displayed by a 5 000 K isosurface temperature. Electrodes are colored by the electrical potential value.


Pierre Lescanne

An exercise on streams: convergence acceleration

Keywords: convergence acceleration, streams, numerical analysis, co-algebra

published or not. The documents may come

Introduction

Assume that we want to compute numerically the limit of a sequence that converges slowly. If we use the sequence itself, we will get significant figures of the limit after a long time. Methods called convergence acceleration have been designed to shorten the time after which we get reasonable amount of significant digits. In other words, convergence acceleration is a set of methods for numerically computing the limit of a sequence of numbers. Those methods are based on sequence transformations and are a nice domain of application of streams [START_REF] Jan | Behavioural differential equations: a coinductive calculus of streams, automata, and power series[END_REF][START_REF] Bonsangue | Defining context-free power series coalgebraically[END_REF], with beautiful higher order functions. This allows us to present elegantly rather complex methods and to code them in Haskell [START_REF] Hutton | Programming in Haskell[END_REF] replacing long, obscure and special purpose Fortran programs by short, generic and arbitrary precision Haskell programs.

In this paper we show, how given a sequence (s n ) n∈N , we can evaluate efficiently lim n→∞ s n . For that we use Levin transforms. There are three kinds of such transforms, which are the result of three sequence transformations labeled traditionally by t, u and v ( [START_REF] Osada | Acceleration Methods for Slowly Convergent Sequences and their Applications[END_REF], p.58).

Presentation of the method

In what follows we speak indistinctly of "sequences" or of "streams" We use Haskell and we work with arbitrary precision reals based on the implementation of David Lester called CReal. In Haskell the type of streams over a type A is written [A].

For the numerical aspect, we follows Naoki Osada [START_REF] Osada | Acceleration Methods for Slowly Convergent Sequences and their Applications[END_REF] (for a historical account see [START_REF] Brezinski | Convergence acceleration during the 20th century[END_REF]) and we show that the stream notation of Haskell makes the presentation simpler. With no surprise, the size of the Haskell code is the same if not shorter than the mathematical description of the reference [START_REF] Osada | Acceleration Methods for Slowly Convergent Sequences and their Applications[END_REF]. Moreover it provides efficient programs.

Levin transformations are somewhat generic in the sense that they are based on elementary transformations. Specialists of convergence acceleration propose three such elementary transformations. Let s be a sequence on CReal, i.e., s :: [CReal].

We define first a basic sequence transformation on which we will found our elementary transformations:

dELTA :: [CReal ] → [CReal ] dELTA s = zipWith (-) (tail s) s
which means that dELT A(s) n = s n+1 -s n . From this basic sequence transformation we define the three elementary other sequence transformations as follows. A unique function depending on a parameter which is either T or U or V is given. For a computer scientist those names T , U , V are meaningless but this terminology (in lower case i.e., u, t, v is traditionally used by mathematicians and we stick to it. It corresponds to the traditional notations of numerical analysis for those sequence transformations.

data Kind = T | U | V delta :: Kind → [CReal ] → [CReal ] delta T s = dELTA s delta U s = zipWith ( * ) (dELTA s) [1..] delta V s = zipWith (/) (zipWith ( * ) (tail $dELTA s) (dELTA s)) (dELTA (dELTA s))
In numerical analysis, people speak about E-algorithm. This is a family of functions eAlg n,k which are also parametrized by a character either T or U or V . It tells which of the basic sequence transformations is chosen. eAlg n,k uses a family of auxiliary functions which we call gAlg n,k for symmetry and regularity. Here is the Haskell code for these functions:

eAlg :: Kind → Int → [CReal ] → [CReal ] eAlg c 0 s = s eAlg c k s = let a = (eAlg c (k -1) s) b = (gAlg c (k -1) k s) in zipWith (-) a (zipWith ( * ) b (zipWith (/) (dELTA a) (dELTA b))) gAlg :: Kind → Int → Int → [CReal ] → [CReal ] gAlg c 0 j s = let nTojMinus1 j = zipWith ( * * ) [1..] (repeat (fromIntegral (j -1))) in zipWith (/) (nTojMinus1 j ) (delta c s) gAlg c k j s = let a = gAlg c (k -1) j s b = gAlg c (k -1) k s in zipWith (-) a (zipWith ( * ) b (zipWith (/) (dELTA a) (dELTA b)))
Here is the formula as it is given in [START_REF] Osada | Acceleration Methods for Slowly Convergent Sequences and their Applications[END_REF]. R n is the generic value of (delta c s) n . gAlg is written g. E (n) k is the n th element of the sequence eAlg c k s, the same for g

E (n) 0 = s n , g (n) 0,j = n 1-j R n , n = 1, 2, ...; j = 1, 2, ..., E (n) k = E (n) k-1 -g (n) k-1,k ∆E (n) k-1 ∆g (n) k-1,k , n = 1, 2, ...; k = 1, 2, ..., g (n) k,j = g (n) k-1,j -g (n) k-1,k ∆g (n) k-1,j ∆g (n) k-1,k , n = 1, 2, ...; k = 1, 2, ..., j > k

Levin's formulas

There is another formula for the E -algorithm:

T (n) k = k j=0 (-1) j k j n + j n + k k-1 s n+j R n+j k j=0 (-1) j k j n + j n + k k-1 1 R n+j T (n) k
is not easily amenable to a Haskell program. 1 We give only the functions for k = 0, 1, 2, which we call levin c 0, levin c 1 and levin c 2. For k = 0, T (n) 0 = s n and levin k 1 and levin k 2 are the result of small calculations.

levin :: Kind → Int → [CReal ] → [CReal ] levin c 0 s = s levin c 1 s = zipWith (-) s (zipWith (/) (zipWith ( * ) (dELTA s) (delta c s)) (dELTA (delta c s)))
levin U 1 s is called Aitken's delta-squared process. One notices that the numerator and the denominator of the above formula differ slightly. Indeed s n+j in the numerator is just replaced by 1 in the denominator. More than any other branch of numerical analysis, convergence acceleration is an experimental science. The researcher applies the algorithm and looks at the results to assess their worth.

formulaForLevinTwo :: Kind → [CReal ] → [CReal ] → [CReal ] formulaForLevinTwo c s ′ s = zipWith (
Dirk Laurie [START_REF] Dirk | Appendix A: Convergence acceleration[END_REF] Given a sequence we use the convergence acceleration to find the main coefficient of the asymptotic equivalent. More precisely given an integer sequence s n , we want to find a number a > 1 such that

s n ∼ a n f (n) where f (n) is subexponential. In other words, for all b ∈ R, b > 1, lim n→∞ f (n)
b n = 0. Actually we compute the limit of the sequence s n+1 /s n . For that we create the function:

expCoeffAC :: ([CReal ] → [CReal ]) → [Integer ] → Int → CReal expCoeffAC transform sequence n = last (transform (zipWith (/) (tail u) u))
where u = map fromIntegral (take n sequence)

Thus expCoeffAC (levin U 2) s 300 gives the approximation of the coefficient one can get after 300 iterations using the sequence transformation levin U 2.

Catalan numbers

A good example to start is Catalan numbers:

catalan = 1 : [let cati = take i catalan in sum (zipWith ( * ) cati (reverse cati)) | i ← [1..]]
We know that catalan!!n ∼ 4 n √ πn 3

Actually we get expCoeffAC (levin U 2) catalan 800 = 4.0000000237 (with 8 exacts digits) and

expCoef f AC (eAlg T 2) catalan 800 ≈ 3.9849561088 expCoef f AC (eAlg U 2) catalan 800 ≈ 3.9773868157 expCoef f AC (eAlg V 2) catalan 800 ≈ 3.9773869346

Counting plain lambda terms

Now we want to use this technique to address a conjecture2 [START_REF] Grygiel | Counting terms in the binary lambda calculus[END_REF], on the asymptotic evaluation of the exponential coefficient of the numbers of typable terms of size n when n goes to ∞. First let us give the recursive definition of the numbers S ∞ of plain lambda terms of size n. This sequence appears on the On-line Encyclopedia of Integer Sequences with the entry number A114851. We assume that abstractions and applications have size 2 and variables have size 1 + k where k is the depth of the variable w.r.t. its binder.

S ∞,0 = S ∞,1 = 0, S ∞,n+2 = 1 + S ∞,n + n k=0 S ∞,k S ∞,n-k .
It has been proved in [START_REF] Grygiel | Counting terms in the binary lambda calculus[END_REF] that

S ∞,n ∼ A n • C n 3/2 ,
where A . = 1.963447954 and C . = 1.021874073. After 300 iterations and using levin U 2 we found 1.9634489522735283291619147713569993355616.

giving six exact digits.

Conclusion

We have shown how streams can be applied to numerical analysis, namely to convergence acceleration. It makes no doubt that they can also be applied to other fields in numerical analysis or elsewhere. For instance, one may imagine applications of acceleration of convergence to the computation of of non numerical sequences.

We did not use the differential equation approach [START_REF] Jan | Behavioural differential equations: a coinductive calculus of streams, automata, and power series[END_REF][START_REF] Winter | QStream: A suite of streams[END_REF], but presenting acceleration of convergence in this framework should be worthwhile.

  . ∆ is the notation for what we write dELTA.

  +) (zipWith (-) (foldl (zipWith ( * )) [2..] [tail $tail s ′ , tail $delta c s, delta c s]) (foldl (zipWith ( * )) [2, 4..] [tail s ′ , tail $tail $delta c s, delta c s])) (foldl (zipWith ( * )) [0..] [s ′ , tail $tail $delta c s, tail $delta c s]) levin c 2 s = zipWith (/) (formulaForLevinTwo c s s) (forLevinTwo c [1, 1..] s) Brezinski [3] proves that the sequences E (n) k and T (n) k are the same, in other words: eAlg = η levin 4 Classic and non classic examples

In particular due to many divisions by 0 and to the complexity of the formula.

This problem is the origin of the interest of the author for this question.

Counting typable lambda terms

The question is now to find the exponential coefficients for the numbers of typable terms. We have no formula for computing those numbers. The only fact we know is the following table of the numbers T ∞,n till 42 which has been obtained after heavy computations (more than 5 days for the 42 nd ). The method consists in generating all the lambda terms of a given size and sieving those that are typable to count them. Therefore the best method to guess exponential coefficient is by acceleration of convergence. After 43 iterations we found 1.8375065809.... Knowing that with the same number 43 we get 1.8925174623... for S ∞ , this is not enough to conclude. But this allows us to speculate that the exponential coefficient for the asymptotic evaluation of T ∞,n could be 1.963447954 like for S ∞,n 's.

Application to divergent series

In his famous paper [START_REF] Euler | De seriebus divergentibus[END_REF] Euler provides a sum to divergent series. See [START_REF] Ramis | Sries divergentes et thories asymptotiques[END_REF][START_REF] Martinet | Elementary acceleration and multisummability[END_REF] and for a light introduction, the reader who understands French is advised to watch the video [START_REF] Ramis | Leonhard euler, ou l'art de donner un sens ce qui n'en avait pas[END_REF] which completes another video [START_REF]+ 4 + 5 +[END_REF] in English.

Among the methods Euler and his followers propose to give a meaning to sum of divergent series there is convergence acceleration. We applied naturally our implementation to some divergent series.

Lets us define the function that associated to a sequence it series.

in map fromIntegral (series s)

Grandi series

Grandi series is also called by Euler, Leibniz series. This is the series

which is sometime written 1 -1 + 1 -1 + . . .. In Haskell it is: grandi = let gr = 1 : (-1) : gr in seq2series gr

We get 1/2 after 3 iterations using eAlg T 2 and this does not change when we increase the number of iterations.


Alexandre Aksenov

A combinatorial approach to rarefaction in b-multiplicative sequences.

par Alexandre AKSENOV Résumé. Pour une suite b-multiplicative donnée et un nombre premier p fixé, l'étude de la p-raréfaction consiste à caractériser le comportement asymptotique des sommes de premiers termes d'indices multiples de p. Sous une hypothèse (dite de ≪finitude≫) sur la suite, les valeurs entières du polynôme homogène ≪norme≫ à 3 variables, N p,i1,i2 (Y 0 , Y 1 , Y 2 ) := p-1 j=1 Y 0 +ζ i1j p Y 1 +ζ i2j p Y 2 , où i 1 , i 2 ∈{1, 2, . . . , p-1}, ζ p est une racine p-ième primitive de l'unité, déterminent ce comportement asymptotique. On montre qu'une méthode combinatoire s'applique à N p,i1,i2 (Y 0 , Y 1 , Y 2 ) qui permet d'établir de nouvelles relations fonctionnelles entre les coefficients de ce polynôme ≪norme≫, diverses propriétés des coefficients de N p,i1,i2 (Y 0 , Y 1 , Y 2 ), notamment pour i 1 =1, i 2 =2, 3; cette méthode fournit des relations entre coefficients binomiaux, de nouvelles preuves des deux identités p-1 j=1 1+ζ j p -ζ 2j p =L p (le p-ième nombre de Lucas) et p . Une méthode algorithmique de recherche des coefficients de N p,i1,i2 est développée.

Abstract. Given a b-multiplicative sequence and a prime p, studying the p-rarefaction consists in characterizing the asymptotic behaviour of the sums of the first terms indexed by the multiples of p. Under the "finiteness" assumption for the sequence, the integer values of the homogeneous "norm" 3-variate polynomial N p,i1,i2 (Y 0 , Y 1 , Y 2 ) := p-1 j=1 Y 0 +ζ i1j p Y 1 +ζ i2j p Y 2 , where i 1 , i 2 ∈{1, 2, . . . , p-1}, and ζ p is a primitive p-th root of unity, determine this asymptotic behaviour. It will be shown that a combinatorial method can be applied to N p,i1,i2 (Y 0 , Y 1 , Y 2 ). The method enables deducing functional relations between the coefficients as well as various properties of the coefficients of N p,i1,i2 (Y 0 , Y 1 , Y 2 ), in particular for i 1 =1, i 2 =2, 3. This method provides relations between binomial coefficients. It gives new proofs of the two identities p-1 j=1 1-ζ j p =p and p-1 j=1 1+ζ j p -ζ 2j p =L p (the p-th Lucas number). The sign and the residue modulo p of the symmetric 2010 Mathematics Subject Classification. 05 A 10, 05 A 18, 11 B 39, 11 R 18.

Introduction

This article deals with a combinatorial method adapted to the coefficients of homogeneous 3-variate "norm" polynomials which determine the asymptotic behaviour of rarified sums of a sub-class of b-multiplicative sequences. The general definition of a b-multiplicative sequence of complex numbers can be written as:

t n+mb i = t n t mb i for each n, m, i ∈ N such that n < b i .
We will be interested in the case where t mb i = t m (this condition will be called the finiteness condition). If a b-multiplicative sequence satisfies the finiteness condition and its values are either 0 or roots of unity, it is b-automatic. An example of such sequence is the Thue-Morse sequence defined by b = 2, t 1 = -1. A survey on the b-multiplicative sequences with values in an arbitrary compact group can be found in [START_REF] Drmota | Generalized Thue-Morse sequence of Squares[END_REF].

Rarified sums (the term is due to [START_REF] Drmota | Rarified sums of the Thue-Morse sequence[END_REF]) of a sequence (t n ) are the sums of initial terms of the subsequence (t pn ) n (the rarefaction step p is supposed to be a prime number in this paper). It is proved in [START_REF] Goldstein | The Fractal Structure of Rarefied Sums of the Thue-Morse Sequence[END_REF] that if (t n ) is the Thue-Morse sequence and b = 2 is a generator of the multiplicative group

F × p , then (1.1)
n<N,p|n

t n = O N log p (p-1) log 2
and this exponent cannot be decreased. For some prime numbers p the rarified sums are always positive, this phenomenon is discussed in [START_REF] Drmota | Rarified sums of the Thue-Morse sequence[END_REF].

The equation (1.1) can be generalized as:

(1.2) n<N,p|n t n = O N log N Q(ζp)/Q ( b-1 j=0 t j ζ j p ) (p-1) log b
for any b-multiplicative sequence (t n ) satisfying the finiteness condition, with values only in {-1, 0, 1}, such that

N Q(ζp)/Q   b-1 j=0 t j ζ j p   > max   ( b-1 j=0 t n ) p-1 , 1  
and such that b < p is a generator of the multiplicative group F × p . The equation (1.2) uses the notation ζ p for a primitive p-th root of unity and N L/K for the norm. This result is proved in the Ph.D. thesis [START_REF] Aksenov | Raréfaction dans les suites b-multiplicatives[END_REF]. Indeed, the equation (1.2) generalizes (1.1) as the Thue-Morse sequence satisfies the conditions of validity of (1.2) and we get the following:

N Q(ζp)/Q   b-1 j=0 t j ζ j p   = N Q(ζp)/Q (1 -ζ p ) = p.
The norms

(1.3) N Q(ζp)/Q   b-1 j=0 t j ζ j p   can also be calculated in a straightforward way if b = 3, t 0 =t 1 =1, t 2 =-1.
Using the resultant of the two polynomials S(X) = X p-1 + . . . + 1 and R(X)=X 2 -X-1, one obtains

(1.4) N Q(ζp)/Q (1 + ζ p -ζ 2 p ) = L p the p-
th term of the Lucas sequence (referred as A000032 by OEIS, cf [START_REF]On-Line Encyclopedia of Integer Sequences[END_REF]) defined recursively by

L 0 = 2, L 1 = 1, L n+2 = L n +L n+1 .
The objective of this article is to study a "norm" expression similar to (1.3) in the case of some b-multiplicative sequences, which is constructed by introducing a finite number of new formal variables Y 0 , Y 1 , . . . , Y d with d p -1. The integer d can be defined as the number of nonzero terms in t 1 , t 2 , . . . , t b-1 . The combinatorial method developed here concerns the case d = 2; in general, if d 3, it leads to too difficult computations 1 . Then, we are interested in the homogeneous polynomial of degree p -1 in 3 variables with integer coefficients

(1.5) N p,i 1 ,i 2 (Y 0 , Y 1 , Y 2 ) = p-1 j=1 Y 0 + ζ i 1 j p Y 1 + ζ i 2 j p Y 2 ,
where i 1 , i 2 are two distinct elements of

F × p . If all numbers t c for c ∈ {1, . . . , b -1} \ {i 1 , i 2 } are zeroes, then the norm (1.3) is recovered as N p,i 1 ,i 2 (1, t i 1 , t i 2 ). By definition, N p,i 1 ,i 2 (Y 0 , Y 1 , Y 2 ) is the norm of (Y 0 + ζ i 1 p Y 1 + ζ i 2 p Y 2 ) as a polynomial in the 4 variables Y 0 , Y 1 , Y 2 , ζ p relative to the extension of fields Q(ζ p )/Q
in the sense of the extended definition of norm introduced in [START_REF] Trager | Algebraic Factoring and Rational Function Integration[END_REF].

The form (1.5) of "norm" polynomial reveals to be common for a large class of b-multiplicative sequences, either by setting the formal variables Y 0 , Y 1 , Y 2 to special values and/or fixing the two residue classes i 1 , i 2 . Since the form (1.5) inherits the properties of its coefficients, any functional relation between these coefficients can be considered as a key result.

Alexandre AKSENOV

In this context, Sections 2 and 3 enunciate a combinatorial interpretation of the coefficients of N p,i 1 ,i 2 and the proof of a functional relation between them, which looks like the recurrence equation of the Pascal's triangle. This relation can be used to find closed formulas for some classes of coefficients (for all of them in the case i 1 = 1, i 2 = 2) and to find the remaining coefficients in a fast algorithmic way. A closed formula for these coefficients is a final goal.

In Section 4 we describe an algorithm in O(p 2 ) additions that calculates the coefficients of (1.5) using this relation. We also re-prove the result (1.4) about the case i 1 = 1, i 2 = 2 and formulate two corollaries of the new proof. We also state some results about the case i 1 = 1, i 2 = 3.

Throughout the paper, |X| and #X will both refer to the size of a finite set X, the symbol # followed by a system of equations, congruences or inequalities will denote the number of solutions; and X, standing for x∈X x, will refer to the sum of a finite subset X of a commutative group with additive notation.

Combinatorics of partitions of a set.

In this section we are going to give an alternative proof of the formula

(2.1) j=p-1 j=1 X -ζ j = 1 + X + . . . + X p-1 ,
and the methods of this proof will be re-used in the proof of the functional equation in Section 3. The new proof uses the properties of the partially ordered sets Π n of partitions of a set of size n (a good reference about the properties of those is the Chapter 3.10.4 of [START_REF] Stanley | Enumerative combinatorics[END_REF]). We are going to prove the following statement, which is equivalent to (2.1). Lemma 2.1. Let p be a prime number and 0 n < p an integer. Define A 0 (n, p) as the number of subsets of F × p of n elements that sum up to 0 modulo p and A 1 (n, p) the number of those subsets that sum up to 1. Then

A 0 (n, p) -A 1 (n, p) = (-1) n .
Let us begin the proof with an obvious observation: if we define similarly the numbers A 2 (n, p), A 3 (n, p), . . . , A p-1 (n, p) , they will all be equal to A 1 (n, p), since multiplying a set that sums to 1 by a constant residue c ∈ F × p gives a set that sums to c, and this correspondence is one-to-one.

Let us deal with a simpler version of the Lemma that allows repetitions and counts sequences instead of subsets, which is formalized in the following

Definition 1. Denote E k 1 ,k 2 ,...,
k i x i = x.
Then we get the following Lemma 2.2. If n is even,

E k 1 ,k 2 ,...,kn 0 (n, p) = (p-1) n +p-1 p and E k 1 ,k 2 ,...,kn 1 (n, p) = (p-1) n -1 p ;
if n is odd,

E k 1 ,k 2 ,...,kn 0 (n, p) = (p-1) n -p+1 p and E k 1 ,k 2 ,...,kn 1 (n, p) = (p-1) n + 1 p .
In both cases,

E k 1 ,k 2 ,...,kn 0 (n, p) -E k 1 ,k 2 ,...,kn 1 (n, p) = (-1) n .
Proof. By induction on n. For n = 0 or n = 1 the result is trivial. For bigger n we always get:

E k 1 ,k 2 ,...,kn 0 (n, p) = (n -1)E k 1 ,k 2 ,...,k n-1 1 (n -1, p) and E k 1 ,k 2 ,...,kn 1 (n, p) = E k 1 ,k 2 ,...,k n-1 0 (n -1, p) + (p -2)E k 1 ,k 2 ,...,k n-1 1 (n -1, p),
since the sequences of length n of linear combination (with coefficients k i ) equal to x are exactly expansions of sequences of length n -1 of linear combination different from x, and this correspondence is one-to-one. Injecting formulas for n -1 concludes the induction. Now we are going to prove Lemma 2.1 for small n.

If n = 0 or n = 1, Lemma is clear. For n = 2, there is one more sequence (x, y) ∈ F × p 2
that sums up to 0, but that counts the sequences of the form (x, x) which should be removed. Since p is prime, these sequences contribute once for every nonzero residue modulo p, and removing them increases the zero's "advantage" to 2. Now, we have to identify (x, y) and (y, x) to be the same, so we get the difference 1 back, establishing Lemma 1 for n = 2.

For n = 3, counting all the sequences (x, y, z) ∈ F × p gives a difference E 0 -E 1 = -1. The sequences (x, x, z) contribute one time more often to the sum equal to 0, so removing them adds -1 to the total difference. The same thing applies to sequences of the form (x, y, y) and (x, y, x). After removing them, we get an intermediate difference of -4, but the triples of the form (x, x, x) have been removed 3 times, which is equivalent to saying they count -2 times. Therefore, they should be "reinjected" with coefficient 2. As p is prime and bigger than 3, the redundant triples contribute once for each nonzero residue; therefore we accumulate the difference of -4 -2 = -6. We have then to identify permutations, that is to divide the score by 6 which gives the final result -1.

Here is the explicit calculation for the case n = 4:

1 (corresponds to E 0 (4, p) -E 1 (4, p))

+6

(for removing (x, x, y, z), (x, y, x, z), (x, y, z, x), (x, y, y, z), (x, y, z, y), (x, y, z, z) )

+2 × 4 (for re-injecting (x, x, x, y), (x, x, y, x), (x, y, x, x) and (x, y, y, y) )

+1 × 3 (for re-injecting (x, x, y, y), (x, y, x, y) and (x, y, y, x))

+6 × 1 (for removing (x, x, x, x)) = 24
which is 4!, therefore Lemma 2.1 is proved for n = 4.

For a general n we can calculate the difference between the number of sequences that sum up to 0 and the number of those that sum up to 1 by assigning to all sequences in F × p n an intermediate coefficient equal to one, then by reducing it by one for each couple of equal terms, then increasing by 2 for each triple of equal terms, and so on, proceeding by successive adjustments of coefficients, each step corresponding to a "poker combination" of n cards. If after adding the contributions of all the steps and the initial (-1) n , we get (-1) n n!, Lemma 2.1 is valid for n independently from p provided that p > n is prime.

Let us introduce a formalization of these concepts using the notions exposed in [START_REF] Kung | Combinatorics: The Rota Way[END_REF]. Call a partition of the set {1, 2, . . . , n} a choice of pairwise disjoint nonempty subsets B 1 , B 2 , . . . , B c of {1, 2, . . . , n} of non-increasing sizes |B i |, and such that B 1 ∪ B 2 ∪ . . . ∪ B c = {1, 2, . . . , n}. The set Π n of all partitions of {1, 2, . . . , n} is partially ordered by reverse refinement: for each two partitions τ and π, we say that τ π if each block of π is included in a block of τ . We define the Möbius function µ( 0, x) on Π n recursively by: if x = {{1}, {2}, . . . , {n}} = 0, then µ( 0, x) = 1; if x is bigger than 0, then

µ( 0, x) = - y ∈ Π n y < x µ( 0, y).
By the Corollary to the Proposition 3 section 7 of [START_REF] Rota | On the Foundations of Combinatorial Theory I. Theory of Möbius Functions[END_REF] and the first Theorem from the section 5.2.1 of [START_REF] Kung | Combinatorics: The Rota Way[END_REF]

, if x is a subdivision of type (λ 1 , λ 2 , . . . , λ n ), then (2.2) µ( 0, x) = n i=1 (-1) λ i -1 (λ i -1)!
This formula will be useful in Section 3.

We are also going to use the following definition: let x = (x 1 , x 2 , . . . , x n ) be a sequence of n nonzero residues modulo p seen as a function

x : {1, 2, . . . , n} → F × p . Then the coimage of x is the partition of {1, 2, . . . , n}, whose blocks are the nonempty preimages of elements of F × p . Now we can prove the following proposition that puts together all the previous study.

Lemma 2.3. The difference A 0 (n, p) -A 1 (n, p)
does not depend on p provided that p is a prime number bigger than n.

Proof. We are going to describe an algorithm that computes this difference (which is the one applied earlier for small values of the argument). For each subdivision x ∈ Π n , denote by r 0 (x, p) the number of sequences (x 1 , x 2 , . . . , x n ) of elements of F × p of coimage x that sum up to 0, and denote by r 1 (x, p) the number of those sequences of coimage x that sum up to 1 and denote r(x, p) = r 0 (x, p)r 1 (x, p). Then, where c(y) is the number of blocks in the subdivision y. By the Möbius inversion formula (see [START_REF] Kung | Combinatorics: The Rota Way[END_REF]),

n!(A 0 (n, p) -A 1 (n, p)) = r( 0, p).
(2.4) r( 0, p) = y∈Πn µ( 0, y)s(y, p) = y∈Πn (-1) c(y) µ( 0, y).
If we compute this sum, we get the value of A 0 (n, p) -A 1 (n, p) in a way that does not depend on p.

The last move consists in proving that (2.5) y∈Πn (-1) c(y) µ( 0, y) = (-1) n n! in a way that uses the equivalence with Lemma 2.1. This proof may seem to be artificial because it is no longer used in the Section 3, and a purely combinatorial and more general proof exists: see the final formula of Chapter 3.10.4 of [START_REF] Stanley | Enumerative combinatorics[END_REF].

Remark that A 0 (n, p) = A 0 (n, p-1-n) since saying that the sum of some subset of F × p is 0 is equivalent to saying that the sum of its complement is 0. For the same kind of reason,

A 1 (n, p) = A -1 (n, p-1-n) = A 1 (n, p-1-n).
Now we can prove Lemma 2.1 by induction on n. It has already been proved for small values of n. If n > 4, by Bertrand's postulate, there is a prime number p ′ such that n < p ′ < 2n. Replace p by p ′ (by the proposition 2.3 this leads to an equivalent statement), then n by p ′ -1n (using the above remark). As p ′ -1n < n, the step of induction is done.

This proof can be analysed from the following point of view: how fast does the number of steps of induction grow as function of n? Suppose that one step of induction reduces Lemma 2.1 for n to Lemma 2.1 for the number f (n) and denote by R(n) the number of steps of induction needed to reach one of the numbers 0 or 1 (the formal definitions will follow). We can prove then the following upper bound on R(n). This definition makes sense, for f (n) < n for each n > 1 by the Bertrand's postulate.

The function R(n) satisfies the estimation

(2.6) R(n) = O(log log n).
Proof. Denote θ = 0.525. By Theorem 1 of [START_REF] Baker | The difference between consecutive primes -II[END_REF], there is a constant N 0 such that for all n > N 0 , the interval [nn θ , n] contains a prime number. We are going to deduce from this the following result: for each θ ∈]0.525, 1[ there exists a constant

N 1 such that n > N 1 implies f (n) < n θ.
Indeed, suppose n > N 0 and denote p = nextprime(n) -1. Then, by the result cited above, (2.7)

n ppθ .

The function

u : [N 0 , +∞[ → [u(N 0 ), +∞[ x → x -x θ
is strictly increasing, continuous and equivalent to x. Therefore, the same is valid for its inverse u -1 . By (2.7), p u -1 (n), therefore

(2.8) ∀n > N 1 f (n) = p -n pθ (u -1 (n)) θ < n θ
for each θ ∈]θ, 1[ and for a bound N 1 N 0 that may depend on θ. The end of the proof is analogous to that of Theorem 1.1 of [START_REF] Luca | On an arithmetic function considered by Pillai[END_REF]. Denote by l the integer such that f l+1 (n) < N 1 f l (n). Then:

n θl N 0 therefore l log θ + log log n log log N 1 which implies l - log log n log θ . Put b = max 1 m N 0 R(m)
, it is a constant. We get:

R(n) l + 1 + b - log log n log θ + 1 + b
which proves our claim.

3. Pascal's equation.

We are going to prove the functional equation verified by the coefficients of the polynomial P i 1 ,i 2 (Y 0 , Y 1 , Y 2 ) (introduced in (1.5)). To do this, we are going to describe a combinatorial interpretation of these numbers. Definition 2. Let p, i 1 , i 2 be fixed as in Introduction and n 1 , n 2 be nonnegative integers such that n 1 + n 2 p -1. Define

C i 1 ,i 2 i (n 1 , n 2 , p) = (3.1) #      (x 1 , . . . , x n 1 +n 2 ) ∈ F × p n 1 +n 2 x k = x l if k = l, i 1 n 1 k=1 x k + i 2 n 1 +n 2 k=n 1 +1 x k = i      and A i 1 ,i 2 i (n 1 , n 2 , p) = #    (X 1 , X 2 ) ∈ P(F × p ) 2 |X 1 | = n 1 , |X 2 | = n 2 , X 1 ∩ X 2 = ∅, i 1 X 1 + i 2 X 2 = i    . (3.2)
Alexandre AKSENOV Definition 2 matches with the notations from the previous section because of the identity

A i 1 ,i 2 i (n, 0, p) = A i 1 ,i 2 i (0, n, p) = A i (n, p) (indepen- dently from i 1 , i 2 ).
From this definition one can see that

C i 1 ,i 2 1 (n 1 , n 2 , p) = . . . = C i 1 ,i 2 p-1 (n 1 , n 2 , p), p-1 i=0 C i 1 ,i 2 i (n 1 , n 2 , p) = (p -1) . . . (p -n 1 -n 2 ),
and for any i,

A i 1 ,i 2 i (n 1 , n 2 , p) = C i 1 ,i 2 i (n 1 ,n 2 ,p) n 1 !n 2 !
. Only one linear equation should be added to these in order to be able to determine all the numbers defined by (3.1) and (3.2). Proposition 3.1 below suggests to research the value of

△ i 1 ,i 2 (n 1 , n 2 , p) = A i 1 ,i 2 0 (n 1 , n 2 , p) -A i 1 ,i 2 1 (n 1 , n 2 , p) = X 1 , X 2 ⊂ F × p |X 1 | = n 1 , |X 2 | = n 2 , X 1 ∩ X 2 = ∅ ζ i 1 X 1 +i 2 X 2 p .
We can express the symmetric polynomials of the quantities (Y

0 +ζ i 1 j p Y 1 + ζ i 2 j
p Y 2 ) in terms of the previously defined numbers via the following Proposition 3.1. Let i 1 , i 2 be two different elements of F × p and denote by σ v,(j=1,...,p-1) the elementary symmetric polynomial of degree v in quantities that depend on an index j varying from 1 to p -1. Then we have the following formal expansion:

(3.3) σ p-1-δ,(j=1,...,p-1) Y 0 + ζ i 1 j p Y 1 + ζ i 2 j p Y 2 = 0 n 0 , n 1 , n 2 p -1 n 0 + n 1 + n 2 = p -1 n 0 δ n 0 δ △ i 1 ,i 2 (n 1 , n 2 , p)Y n 0 -δ 0 Y n 1 1 Y n 2 2 .
In particular,

(3.4) N p,i 1 ,i 2 (Y 0 , Y 1 , Y 2 ) = 0 n 0 , n 1 , n 2 p -1 n 0 + n 1 + n 2 = p -1 △ i 1 ,i 2 (n 1 , n 2 , p)Y n 0 0 Y n 1 1 Y n 2 2 .
Proof. The symmetric polynomial develops as:

σ p-1-δ,(j=1,...,p-1) Y 0 + ζ i 1 j p Y 1 + ζ i 2 j p Y 2 = X ⊂ F × p , |X| = p-1-δ j∈X Y 0 + ζ i 1 j p Y 1 + ζ i 2 j p Y 2 = X ⊂ F × p , |X| = p-1-δ X 0 , X 1 , X 2 , X 0 ∪ X 1 ∪ X 2 = X, X 0 , X 1 , X 2 disjoint ζ i 1 X 1 +i 2 X 2 p Y |X 0 | 0 Y |X 1 | 1 Y |X 2 | 2 = X ′ 0 , X 1 , X 2 , X ′ 0 ∪ X 1 ∪ X 2 = F × p X ′ 0 , X 1 , X 2 disjoint |X ′ 0 | δ ζ i 1 X 1 +i 2 X 2 p Y |X ′ 0 |-δ 0 Y X 1 1 Y X 2 2 .
When we group the terms of this sum by sizes

n 0 = |X ′ 0 |, n 1 = |X 1 |, n 2 = |X 2 | we obtain (3.3).
The following closed formula for the coefficients △ i 1 ,i 2 (n 1 , n 2 , p) is valid under a condition: the multiset consisting of i 1 with multiplicity n 1 and of i 2 with multiplicity n 2 should have no nonempty subset of sum multiple of p (this holds, for example, if the smallest positive representatives of i 1 and i 2 verify n 1 i 1 + n 2 i 2 < p). The proof of Lemma 2.1 can then be generalized to get:

(3.5) △ i 1 ,i 2 (n 1 , n 2 , p) = (-1) n 1 +n 2 n 1 + n 2 n 1 .
The complete proof of this statement can be established as a corollary of the Theorem 3.3 below. Without this condition (3.5) becomes false: for example, △ 2,3 (1, 1, 5) = -3. For the general case, we are going to replace the closed formula by a recursive equation in which the parameters i 1 , i 2 , p are fixed, and the recursion is on different values of n 1 , n 2 . The equation is similar to the equation of the Pascal's triangle, and can be formulated as follows: Theorem 3.2 ("Colored" Pascal's equation). Let p be an odd prime, and

i 1 , i 2 ∈ F × p , n 1 , n 2 ∈ {1, . . . , p -2} such that i 1 = i 2 and n 1 + n 2 < p.
Then,

(3.6) C i 1 ,i 2 0 (n 1 , n 2 , p) -C i 1 ,i 2 1 (n 1 , n 2 , p) ≡ n 1 C i 1 ,i 2 1 (n 1 -1, n 2 , p) + n 2 C i 1 ,i 2 1 (n 1 , n 2 -1, p)- -n 1 C i 1 ,i 2 0 (n 1 -1, n 2 , p) -n 2 C i 1 ,i 2 0 (n 1 , n 2 -1, p) mod p and if p ∤ n 1 i 1 + n 2 i 2 , the equality (3.7) C i 1 ,i 2 0 (n 1 , n 2 , p) -C i 1 ,i 2 1 (n 1 , n 2 , p) = n 1 C i 1 ,i 2 1 (n 1 -1, n 2 , p) + n 2 C i 1 ,i 2 1 (n 1 , n 2 -1, p)- -n 1 C i 1 ,i 2 0 (n 1 -1, n 2 , p) -n 2 C i 1 ,i 2 0 (n 1 , n 2 -1, p) holds.
Proof. Define (3.8)

f k = i 1 if k ∈ {1, . . . , n 1 } i 2 if k ∈ {n 1 + 1, . . . , n 1 + n 2 }
We are going to call a hindrance a subset X of {1, . . . , n 1 + n 2 } such that m∈X f m ≡ 0 mod p. If there are no hindrances, then, by following the proof of Lemma 2.1, we get

C i 1 ,i 2 0 (n 1 , n 2 , p) -C i 1 ,i 2 1 (n 1 , n 2 , p) = (-1) n 1 +n 2 (n 1 + n 2 )!
and this number is the opposite of n times (-1) n-1 (n -1)! .

In general, the formula (2.4) should be replaced by:

(3.9) s(y, p) = (1 -p) d(y) (-1) c(y)
if the partition y of {1, . . . , n 1 +n 2 } contains d(y) blocks that are hindrances. Indeed, suppose that the blocks of y are B 1 , . . . , B c , and for each bloc B j we denote f B j = m∈B j f m ∈ F p . Then, the solutions (x 1 , . . . , x n 1 +n 2 ) of

n 1 +n 2 k=1 f k x k = i,
such that the coimages x of (x 1 , . . . , x n 1 +n 2 ) satisfy x y (as partitions), are in one-to-one correspondence to solutions (x B 1 , . . . , x Bc ) of

c j=1 f B j x B j = i
where x B j ∈ F × p are no longer required to be distinct. Proposition 2.2 states that if we pay no attention to the indices j that correspond to hindrances (i.e. such that f B j = 0), the difference between numbers of solutions of j f B j x B j = 0 and j f B j x B j = 1 is (-1) c-d(y) . Moreover, the values of x B j where B j are hindrances can be chosen arbitrarily (from p -1 options each). The product of these contributions leads to (3.9).

The formula (3.9) can be rewritten as

s(y, p) = d(y) l=0
X 1 , X 2 , . . . , X l hindrances contained in y (-1) c(y)-l p l where the order of X 1 , X 2 , . . . , X l is irrelevant in the sum. Then we get:

C i 1 ,i 2 0 (n 1 , n 2 , p) -C i 1 ,i 2 1 (n 1 , n 2 , p) = y∈Πn µ( 0, y)s(y, p) (3.10) = X 1 , X 2 , . . . , X l disjoint hindrances y ∈ Π n y contains X 1 , . . . , X l as blocks (-1) c(y)-l µ( 0, y)p l (3.11) = X 1 ,...,X l (-1) |X 1 |+|X 2 |+...+|X l |-l (|X 1 | -1)!(|X 2 | -1)! . . . (|X l | -1)!p l (3.12) × y ∈ Π n y contains X 1 , . . . , X l µ( 0, y-X 1 -X 2 -. . . -X l )(-1) c(y-X 1 -X 2 -...-X l )
by factoring µ( 0, y) according to the formula (2.2). In the last sum, (y -X 1 -X 2 -. . . -X l ) denotes the partition y, where the blocs X 1 , . . . , X l are removed (which is a partition of (n 1 + n 2 -|X 1 | -. . . -|X l |) elements). By applying (2.5) to the last sum of (3.12), we get

(3.13) C i 1 ,i 2 0 (n 1 , n 2 , p) -C i 1 ,i 2 1 (n 1 , n 2 , p) = X 1 ,...,X l (|X 1 |-1)!(|X 2 |-1)! . . . (|X l |-1)!(-1) n 1 +n 2 -l p l (n 1 +n 2 -|X 1 |-. . . -|X l |)!. From (3.13), (3.14) C i 1 ,i 2 0 (n 1 , n 2 , p) -C i 1 ,i 2 1 (n 1 , n 2 , p) ≡ (-1) n 1 +n 2 (n 1 + n 2 )! mod p, which implies (3.

6).

Alexandre AKSENOV Suppose that {1, . . . , n 1 + n 2 } is not a hindrance. In order to prove (3.7), remark that the sum (3.13) can be split as

X 1 ,...,X l (|X 1 |-1)!(|X 2 |-1)! . . . (|X l |-1)!(-1) n 1 +n 2 -l p l (n 1 +n 2 -|X 1 |-. . . -|X l |)! = - n 1 +n 2 m=1 X 1 , X 2 , . . . , X l disjoint hindrances not containing m (|X 1 |-1)!(|X 2 |-1)! . . . (|X l |-1)! (-1) n 1 +n 2 -l-1 p l (n 1 +n 2 -|X 1 |-. . . -|X l |-1)!
then gathered into two parts according to the values of f m :

C i 1 ,i 2 0 (n 1 , n 2 , p) -C i 1 ,i 2 1 (n 1 , n 2 , p) = -n 1 X 1 , X 2 , . . . , X l disjoint hindrances not containing 1 (|X 1 |-1)!(|X 2 |-1)! . . . (|X l |-1)! (-1) n 1 +n 2 -l-1 p l (n 1 +n 2 -|X 1 |-. . . -|X l |-1)! -n 2 X 1 , X 2 , . . . , X l disjoint hindrances not containing n 1 + 1 (|X 1 |-1)!(|X 2 |-1)! . . . (|X l |-1)! (-1) n 1 +n 2 -l-1 p l (n 1 +n 2 -|X 1 |-. . . -|X l |-1)!
By identifying each sum in the last formula to the right-hand side of (3.13) with one of the arguments n 1 or n 2 decreased by 1, we get (3.7).

The numbers △ i 1 ,i 2 (n 1 , n 2 , p) satisfy a similar equation. Theorem 3.3 ("Uncolored" Pascal's equation). Let p be an odd prime, and

i 1 , i 2 ∈ F × p , n 1 , n 2 ∈ {1, . . . , p -2} such that i 1 = i 2 and n 1 + n 2 < p. Then, (3.15) △ i 1 ,i 2 (n 1 , n 2 , p) ≡ -△ i 1 ,i 2 (n 1 -1, n 2 , p) -△ i 1 ,i 2 (n 1 , n 2 -1, p) mod p and if p ∤ n 1 i 1 + n 2 i 2 , the equality (3.16) △ i 1 ,i 2 (n 1 , n 2 , p) = -△ i 1 ,i 2 (n 1 -1, n 2 , p) -△ i 1 ,i 2 (n 1 , n 2 -1, p) holds.
Proof. Division of both sides of (3.6) by n 1 !n 2 ! (which is not multiple of p) gives (3.15) and division by the same number of (3.7) gives (3.16).

4. Some properties of finite Pascal's triangles.

4.1. Algorithm. Let us define formally

△ i 1 ,i 2 (n 1 , n 2 ) = 0 when one of n 1 , n 2 is negative or n 1 + n 2 p. Then (3.16) is valid for any n 1 , n 2 ∈ N 2 such that p ∤ n 1 i 1 + n 2 i 2 . Indeed: if n 1 = 0 or n 2 = 0, the identification A i 1 ,i 2 (n 1 , n 2 , p) = A(max(n 1 , n 2 ), p) implies (3.16) via Lemma 2.1. When n 1 + n 2 = p -1, one can use the hypothesis X 1 ∪ X 2 = F × p and the identity F × p = 0 to prove i 1 X 1 + i 2 X 2 = (i 1 -i 2 ) X 1 , which implies A i 1 ,i 2 i (n 1 , n 2 , p) = A i 1 -i 2 ,i 2 i (n 1 , 0, p), therefore △ i 1 ,i 2 (n 1 , n 2 ) = (-1) n 1 . The equation (3.16) is valid, therefore, when n 1 + n 2 = p.
We can now prove that the functional relation (3.16), together with these border values, characterizes the function △ i 1 ,i 2 (•, •, p) as a function defined on Z 2 -1 , with values in Z.

Theorem 4.1. Let p be an odd prime, let i 1 , i 2 be two distinct elements of {1, . . . , p -1}, and let d : Z 2 -1 → Z be a function such that d(0, 0) = 1, (4.1)

d(n 1 , n 2 ) = 0 if n 1 < 0, n 2 < 0 or n 1 + n 2 p, (4.2) d(n 1 , n 2 ) + d(n 1 -1, n 2 ) + d(n 1 , n 2 -1) = 0 if p ∤ n 1 i 1 + n 2 i 2 . (4.3) Then, d(n 1 , n 2 ) = △ i 1 ,i 2 (n 1 , n 2 , p). Proof. Define δ(n 1 , n 2 ) = d(n 1 , n 2 ) -△ i 1 ,i 2 (n 1 , n 2 , p).
Then the function δ satisfies (4.2), (4.3) and δ(0, 0) = 0. In order to prove the theorem we should prove that δ = 0.

By applying (4.3) successively to n 2 = 0 and n 1 = 1, . . . , p -1 one proves that δ(0, 0) = -δ(1, 0) = δ(2, 0) = . . . = δ(p-1, 0). By applying it to n 1 = 0 and n 2 = 1, . . . , p -1 one proves that δ(0, 0) = -δ(0, 1) = . . . = δ(0, p -1).

Let us prove the identity δ(n 1 , n 2 ) = 0 by induction on ñ := p-n 1 -n 2 ∈ {0, . . . , p -2}. If ñ = 0, then δ(n 1 , n 2 ) = 0 as a part of the hypothesis (4.2).

Suppose that the Theorem is proved for ñ ∈ {0, . . . , p -3}, let us prove it for ñ + 1. Denote (n S 1 , n S 2 ) the solution of

   i 1 n S 1 + i 2 n S 2 ≡ 0 mod p n S 1 + n S 2 = p - ñ (n S 1 , n S 2 ) ∈ {1, . . . , p} 2 .
If one applies the functional relation (4.3) to a point where n 2 = pñn 1 (with the restriction n 1 = n S 1 ), and uses the induction hypothesis, one gets (4.4)

δ(n 1 -1, p -ñ -n 1 ) + δ(n 1 , p -ñ -n 1 -1) = 0.
By applying (4.4) successively to n 1 = 1, . . . , n S 1 -1, we prove δ(n 1 , pñn 1 -1) = 0 for n 1 in the same range 1, . . . , n S 1 -1. If n S 1 pñ -1, this concludes the step of induction. Otherwise, by applying (4.4) successively to n 1 = pñ -1, . . . , n S 1 + 1 (in the decreasing order of values of n 1 ), we prove δ(n 1 , pñn 1 -1) = 0 for n 1 in the range pñ, . . . , n S 1 . This concludes the induction and proves δ(n 1 , n 2 ) = 0 for all (n 1 , n 2 ).

The previous proof corresponds to the Algorithm 1, which computes the values of the function △ a,b (x, y, p) line by line. It executes one addition per number to compute, therefore its execution time is proportional to the size of the answer.

Given an odd prime p and two distinct elements i 1 , i 2 of F × p , we are going to call the array of all values of △ i 1 ,i 2 (n 1 , n 2 , p) for n 1 , n 2 0, n 1 + n 2 < p a finite Pascal's triangle, and we will use geometrical terminology when it seems to make exposition simpler.

We are going to call sources the points (n 1 , n 2 ) such that p|i

1 n 1 + i 2 n 2 . Define (4.5) f i 1 ,i 2 (n 1 , n 2 , p) = △ i 1 ,i 2 (n 1 , n 2 , p)+△ i 1 ,i 2 (n 1 -1, n 2 , p)+△ i 1 ,i 2 (n 1 , n 2 -1, p).
The value of f i 1 ,i 2 (n 1 , n 2 , p) (which we will call force) is nonzero only at sources, where it can be computed using (3.13) combined with the end of the proof of Theorem 3.2:

(4.6) n 1 !n 2 !f i 1 ,i 2 (n 1 , n 2 , p) = X 1 , X 2 , . . . , X l partition of {1, . . . , n 1 + n 2 }, ∀j p| m∈X j f m (|X 1 |-1)!(|X 2 |-1)! . . . (|X l |-1)!(-1) n 1 +n 2 -l p l (n 1 +n 2 -|X 1 |-. . .-|X l |)!.
This formula uses the notation (3.8) in order to describe the fact that summation goes through all partitions of {1, . . . , n 1 + n 2 } into hindrances. The definition (4.5) implies, by linearity of the Pascal's equation: (4.7)

△ i 1 ,i 2 (n 1 , n 2 , p) = 0 k n 1 0 l n 2 p|i 1 k+i 2 l f i 1 ,i 2 (k, l, p)(-1) n 1 +n 2 -k-l n 1 +n 2 -k-l n 1 -k . +1 -1 -1 +1 +2 +1 -1 -3 -3 -1 +1 +4 +6 +4 +1 -1 -5 -10 -10 -5 -1 +1 +6 +15 +20 +15 -5 +1 -1 -7 -21 -35 +20 -10 +4 -1 +1 +8 +28 -21 +15 -10 +6 -3 +1 -1 -9 +8 -7 +6 -5 +4 -3 +2 -1 +1 -1 +1 -1 +1 -1 +1 -1 +1 -1 +1 Figure 1. Coefficients of 10 j=1 X + ζ j 11 Y + ζ 2j 11 Z 4.2. The case i 1 = 1, i 2 = 2.
We can find a closed formula for the numbers △ 1,2 (n 1 , n 2 , p) using the identity

(4.8) △ 1,2 (n 1 , n 2 , p) = △ 1,2 (n 1 , p -1 -n 1 -n 2 , p).
It follows indeed from the fact that for each disjoint couple X 1 , X 2 ⊂ F × p , as in the definition (3.2),

X 1 + 2 X 2 = - X 1 + 2 (F × p \ X 1 \ X 2 ) .
Formula (3.5) applies to at least one side of (4.8) for each (n 1 , n 2 ) (and to both sides of (4.8) if n 1 + 2n 2 = p -1), leading to (4.9)

△ 1,2 (n 1 , n 2 , p) = (-1) n 1 +n 2 n 1 +n 2 n 1 if n 1 + 2n 2 p -1 (-1) n 2 p-1-n 2 n 1 if n 1 + 2n 2 p -1.
Therefore, this Pascal's triangle is symmetric with respect to the axis n 1 + 2n 2 = p -1.

One can deduce (1.4) from (4.9) in the following way: by (3.4), (4.10)

j=p-1 j=1 1 + ζ p -ζ 2 p = n 1 ,n 2 ∈N (-1) n 2 △ 1,2 (n 1 , n 2 , p) = n 1 ,n 2 ∈N (-1) n 2 -1 (△ 1,2 (n 1 -1, n 2 , p) + △ 1,2 (n 1 , n 2 -1, p) -f 1,2 (n 1 , n 2 , p)) = n 1 ,n 2 ∈N ((-1) n 2 -1 △ 1,2 (n 1 , n 2 -1, p) -(-1) n 2 △ 1,2 (n 1 -1, n 2 , p) + (-1) n 2 f 1,2 (n 1 , n 2 , p)) = n 1 ,n 2 ∈N (-1) n 2 f 1,2 (n 1 , n 2 , p)
because massive cancellation occurs in the sum of differences of values of the function (-1) y △ 1,2 (x, y, p).

Suppose n 1 , n 2 > 0 and

n 1 + 2n 2 = p (therefore n 1 is odd). Then (4.11) f 1,2 (n 1 , n 2 , p) = △ 1,2 (n 1 -1, n 2 , p) + △ 1,2 (n 1 , n 2 -1, p) + △ 1,2 (n 1 , n 2 , p) = (-1) n 2 n 1 + n 2 -1 n 1 -1 + 2(-1) n 2 n 1 + n 2 -1 n 1 = (-1) n 2 n 1 + n 2 n 1 + n 1 + n 2 -1 n 1 = (-1) n 2 p -n 2 n 2 + p -n 2 -1 n 2 -1 .
The absolute value of (4.11) can be interpreted as the number of ways to put n 2 identical disjoint dominoes on a discrete circle of length p. Indeed (see also [START_REF] Benjamin | The Fibonacci Numbers -Exposed More Discretely[END_REF]), for any k p-1 2 (4.12) #{k disjoint dominoes on a circle of length p} = #{k disjoint dominoes on a line segment of length p}

+ #{k -1 disjoint dominoes on a line segment of length p -2} = p -k k + p -k -1 k -1 .
The sum (4.10) contains three terms not covered by the hypotheses of (4.11): these correspond to n 1 =n 2 =0, n 1 =p, n 2 =0, n 1 =0, n 2 =p and they equal respectively 1, 1 and -1. The overall contribution of these terms can be identified to the number of ways to put 0 dominoes on a discrete circle of length p. Therefore, the norm (4.10) equals to the number of ways to put any number of identical disjoint dominoes on a discrete circle of length p, which is proved in [START_REF] Benjamin | The Fibonacci Numbers -Exposed More Discretely[END_REF] to be L p . For example, if p = 11, the numbers are those of Figure 1 ( denotes a source).

4.13) 1 = j=p-1 j=1 1 -ζ p + ζ 2 p = n 1 ,n 2 ∈N (-1) n 1 △ 1,2 (n 1 , n 2 , p) = n 1 ,n 2 ∈N (-1) n 1 f 1,2 (n 1 , n 2 , p).
We further get:

(4.14) 1 = 1 + n 1 ,n 2 ∈N * (-1) n 1 f 1,2 (n 1 , n 2 , p) = 1 - n 1 ,n 2 ∈N * f 1,2 (n 1 , n 2 , p).
The formula (4.11) leads to the following combinatorial identity2 :

(4.15) If (n 1 , n 2 ) is a point on the upper line of sources, the value of f 1,3 (n 1 , n 2 , p) has a simple expression given by (4.6): In any point (n 1 , n 2 ) such that p n 1 + 3n 2 < 2p, the formula (4.7) takes the following form: (4.24) △ 1,3 (n 1 , n 2 , p) = (-1) n 1 +n 2 n 1 + n 2 n 1 + 0<k n 1 0<l n 2 p = i 1 k+i 2 l f 1,3 (k, l, p)(-1) n 1 +n 2 -k-l n 1 +n 2 -k-l n 1k .

p-1 2 k=1 (-1) k p -k k + p -k -1 k -1 = 0.
We can also compute a simple expression for the forces of sources on the lower line. Suppose that n 1 , n 2 > 0 and n 1 + 3n 2 = 2p. Then, by (4.5) and (4.18), For example, if p = 11, the numbers are those of Figure 2 ( denotes a source).

p- 1 j=1 1

 11 -ζ j p =p, le signe et le résidu modulo p des polynômes symétriques des 1+ζ-ζ 2

  Denote, for each subdivision y of {1, 2, . . . , n}, s(y, p) = x y r(x, p). Then, by Proposition 2.2, (2.3) s(y, p) = (-1) c(y)

Theorem 2 . 1 .

 21 Let nextprime(n) := min{p > n, p prime} and f (n) := nextprime(n)n -1 for each n ∈ N. Further, denote R(n) := min{k, f k (n) ∈ {0, 1}}.

4. 3 .= 1+ζ 3 p

 33 Application: an identity for binomial coefficients. The formulas (3.4) and (4.10) have another application. As 1ζ p + ζ 2 p 1+ζp , we get in a similar way to (4.10):

(

4. 4 .Theorem 4 . 2 .- 1 +1 - 1 +1Figure 2 . 4 . 5 .(- 1 ) n 1 +n 2 n 1 +n 2 n 1 if n 1 + 1 if n 1 + 4

 44211245111114 Figure 2. Coefficients of 10 j=1 X + ζ j 11 Y + ζ 3j 11 Z

( 4 . 1 !n 2

 412 22) f 1,3 (n 1 , n 2 , p) = (n 1 + n 2 -1)!p n 1 !n 2 !because the sum consists of the single term associated to X={1, . . . , n 1 +n 2 }. Under the same hypotheses, (4.7) implies (4.23) △ 1,3 (n 1 , n 2 , p) = (n 1 + n 2 -1)!p n

( 4 .

 4 25) f 1,3 (n 1 , n 2 , p) = △ 1,3 (n 1 , n 2 , p) + △ 1,3 (n 1 -1, n 2 , p) + △ 1,3 (n 1 , n 2 -1, p) = f 2,3 (n 1 , pn 1n 2 , p).By (4.6) (the sum, once again, consists of a single term because2n 1 + 3(p-n 1 -n 2 )=p), (4.26) f 2,3 (n 1 , pn 1n 2 , p) = (-1) n 2 (pn 2 -1)!p n 1 !(pn 1n 2 )! .

  kn x (n, p) (where x ∈ F p and k 1 , k 2 , . . . , k n ∈ F × p ) the number of sequences (x 1 , x 2 , . . . , x n ) of elements of F × p such that

	n
	i=1

In[START_REF] Aksenov | Raréfaction dans les suites b-multiplicatives[END_REF] it is proved that the functional equation (Theorems 3.2 and 3.3) generalizes to the case d 3

The previous proof implies (4.15) in the case of prime p 5. The Zeilberger's algoritm (implemented in Maple 17, see also Chapter 6 of the book[9]) generalizes it for any p 5 congruent to 1 or 5 modulo 6

resolution at the edges

The identity (4.11) leads to (4.17) σ p-1-δ,(j=1,...,p-1)

and the discussion that follows the formula (4.11) identifies each number (-1) n 2 f 1,2 (n 1 , n 2 , p) as the number of ways to put n 2 disjoint dominoes on a discrete circle of length p.