Maxime Dénès
email: maxime.denes@inria.fr

Anders Mörtberg
email: mortberg@chalmers.se

Vincent Siles
email: siles@chalmers.se

A refinement-based approach to computational algebra in COQ *

We describe a step-by-step approach to the implementation and formal verification of efficient algebraic algorithms. Formal specifications are expressed on rich data types which are suitable for deriving essential theoretical properties. These specifications are then refined to concrete implementations on more efficient data structures and linked to their abstract counterparts. We illustrate this methodology on key applications: matrix rank computation, Winograd's fast matrix product, Karatsuba's polynomial multiplication, and the gcd of multivariate polynomials.

Introduction

In the past decade, the range of application of proof assistants has extended its traditional ground in theoretical computer science to mainstream mathematics. Formalised proofs of important theorems like the fundamental theorem of algebra [START_REF] Barendregt | The "fundamental theorem of algebra" project[END_REF], the four colour theorem [START_REF] Gonthier | Formal proof-the four-color theorem[END_REF] and the Jordan curve theorem [START_REF] Hales | The jordan curve theorem, formally and informally[END_REF] have advertised the use of proof assistants in mathematical activity, even in cases when the pen and paper approach was no longer tractable.

But since these results established proofs of concept, more effort has been put into designing an actually scalable library of formalised mathematics. The Mathematical Components project (developing the SSREFLECT library [START_REF] Gonthier | A Small Scale Reflection Extension for the Coq system[END_REF] for the COQ proof assistant) advocates the use of small scale reflection to achieve a nearly comparable level of detail to usual mathematics on paper, even for advanced theories like the proof of the Feit-Thompson theorem. In this approach, the user expresses significant deductive steps while low-level details are taken care of by small computational steps, at least when properties are decidable. Such an approach makes the proof style closer to usual mathematics.

One of the main features of these libraries is that they heavily rely on rich dependent types, which gives the opportunity to encode a lot of information directly into the type of objects: for instance, the type of matrices embeds their size, which makes operations like multiplication easy to implement. Also, algorithms on these objects are simple enough so that their correctness can easily be derived from the definition. However in practice, most efficient algorithms in modern computer algebra systems do not rely on dependent types and do not provide any proof of correctness. We show in this paper how to use this rich mathematical framework to develop efficient computer algebra programs with proofs of correctness. This is a step towards closing the gap between proof assistants and computer algebra systems.

The methodology we suggest for achieving this is the following: we are able to prove the correctness of some mathematical algorithms having all the high-level theory at our disposal and we then refine them to an implementation on simpler data structures that will be actually running on machines. In short, we aim at formally linking convenient high-level properties to efficient low-level implementations, ensuring safety of the whole approach while enjoying better performance thanks to the separation of proofs and computational content.

In the next section, we describe the methodology of refinements. Then, we give two examples of such refinements for matrices in Section 3, and polynomials in Section 4. In Section 5, we give a solution to unify both examples by describing COQEAL 1 , a library built using this methodology on top of the SSREFLECT libraries.

Refinements

Refinements are commonly used to describe successive steps when verifying a program. Typically, a specification is expressed in Hoare logic, then the program is described in a high-level language and finally implemented in C. Each step is proved correct with respect to the previous one. By using several formalisms, one has to trust every translation step or prove them correct in yet another formalism.

Our approach is similar: we refine the definition of a concept to an efficient algorithm described on high-level data structures. Then, we implement it on data structures that are closer to machine representations, once we no longer need rich theory to prove the correctness. Thus the implementation is an immediate translation of the algorithm, see Fig. 1.

However, in our approach, the three layers can be expressed in the same formalism (the Calculus of Inductive Constructions), though they do not use exactly the same features. On one hand, the high-level layers use rich dependent types that are very useful when describing theories because they allow abuse of notations

Morphism lemma

Figure 1: The three steps of refinement and concise statements which quickly become necessary when working with advanced mathematics. On the other hand, the efficient implementations use simple types, which are closer to standard implementations in traditional programming languages. The main advantage of this approach is that the correctness of translations can easily be expressed in the formalism itself, and we do not rely on any additional external proofs.

In the next sections, we are going to use the following methodology to build efficient algorithms from high-level descriptions:

1. Implement an abstract version of the algorithm using SSREFLECT's structures and use the libraries to prove properties about them. Here we can use the full power of dependent types when proving correctness.

2. Refine this algorithm into an efficient one using SSREFLECT's structures and prove that it behaves like the abstract version.

3. Translate the SSREFLECT structures and the efficient algorithm to the lowlevel data types, ensuring that they will perform the same operations as their high-level counterparts.

Matrices

Linear algebra is a natural first test-case to validate our approach, as a pervasive and inherently computational area of mathematics, which is well covered by the SSREFLECT library [START_REF] Gonthier | Point-Free, Set-Free concrete linear algebra[END_REF]. In this section, we will detail the (quite simple) data structure we use to represent matrices and then review two fundamental examples: rank computation and efficient matrix product. This encoding makes many properties easy to derive, but it is inefficient for evaluation. Indeed, a finite function over 'I_m * 'I_n is internally represented as a flat list of m × n values which has to be traversed whenever the function is evaluated. Moreover, having the size of matrices encoded in their type allows to state concise lemmas without explicit side conditions, but it is not always flexible enough when getting closer to machine-level implementation details.

Representation

To be able to implement efficient matrix operations we introduce a low-level data type seqmatrix representing matrices as lists of lists. A concrete matrix is built from an abstract one by mapping canonical enumerations (enum) of ordinals to the corresponding coefficients in the abstract matrix:

Definition seqmx_of_mx (M : 'M[R]_(m,n)) : seqmatrix := [seq [seq M i j | j <-enum 'I_n] | i <-enum 'I_m].
To ensure the correct behaviour of concrete matrices it is sufficient to prove that seqmx_of_mx is injective (== denotes boolean equality):

Lemma seqmx_eqP (M N : 'M[R]_(m,n)) : reflect (M = N) (seqmx_of_mx M == seqmx_of_mx N).
Operations like addition are straightforward to implement, and their correctness is expressed through a morphism lemma, stating that the concrete representation of the sum of two matrices is the concrete sum of their concrete representations:

Definition addseqmx (M N : seqmatrix) : seqmatrix := zipwith (zipwith (fun x y => add x y)) M N.

Lemma addseqmxE : {morph (@seqmx_of_mx m n) : M N / M + N >-> addseqmx M N}.

Here morph is notation meaning that seqmx_of_mx is an additive morphism from abstract to concrete matrices. It is worth noting that we could have stated all our morphism lemmas with the converse operator (from concrete matrices to abstract ones). But these lemmas would then have been quantified over lists of lists, with poorer types, which would have required a well-formedness predicate as well as premises expressing size constraints. The way we have chosen takes full advantage of the information carried by richer types. Like the addseqmx operation, we have developed concrete implementations of most of the matrix operations provided by the SSREFLECT library and proved the corresponding morphism lemmas. Among these operations we can cite: subtraction, scaling, transpose and block operations.

Computing the rank

Now that the basic data structure and operations have been defined, it is possible to apply our approach to an algorithm based on Gaussian elimination which computes the rank of a matrix A = (a i, j ) over a field K. We first specify the algorithm using abstract matrices and then refine it to the low-level structures.

An elimination step consists of finding a nonzero pivot in the first column of A. If there is none, it is possible to drop the first column without changing the rank. Otherwise, there is an index i such that a i,1 = 0. By linear combinations of rows (preserving the rank) A can be transformed into the following matrix B:

B =           0 a 1,2 - a 1,1 ×a i,2 a i,1 • • • a 1,n - a 1,1 ×a i,n a i,1 0 . . . . . . a i,1 a i,2 • • • a i,n 0 . . . . . . 0 a n,2 - a n,1 ×a i,2 a i,1 • • • a n,n - a n,1 ×a i,n a i,1           =             0 R 1 . . . 0 a i,1 • • • a i,n 0 R 2 . . . 0             Now pose R = R 1 R 2 , since a i,1 = 0, this means that r ank A = r ank B = 1 + r ank R.
Hence the current rank can be incremented and the algorithm can be recursively applied on R.

In our development we defined a function elim_step returning the matrix R above and a boolean b indicating if a pivot has been found. A wrapper function rank_elim is in charge of maintaining the current rank and performing the recursive call on R:

Fixpoint rank_elim (m n : nat) {struct n} : 'M[K]_(m,n) -> nat := match n return 'M[K]_(m,n) -> nat with | q.+1 => fun M => let (R,b) := elim_step M in (rank_elim R + b)%N | _ => fun _ => %N end.
Note that booleans are coerced to natural numbers: b is interpreted as 1 if true and 0 if false. The correctness of rank_elim is expressed by relating it to the \rank function of the SSREFLECT library:

Lemma rank_elimP n m (M : 'M[K]_(m,n)) : rank_elim M = \rank M.
The proof of this specification relies on a key invariant of elim_step, relating the ranks of the input and output matrices:

Lemma elim_step_rank m n (M : 'M[K]_(m, 1 + n)) : let (R,b) := elim_step M in \rank M = (\rank R + b)%N.
Now the proof of rank_elimP follows by induction on n. The concrete version of this algorithm is a direct translation of the algorithm using only concrete matrices and executable operations on them. This executable version (called rank_elim_seqmx ) is then linked to the abstract implementation by the lemma:

Lemma rank_elim_seqmxE : forall m n (M : 'M[K]_(m, n)), rank_elim_seqmx m n (seqmx_of_mx M) = rank_elim M.
The proof of this is straightforward as all of the operations on concrete matrices have morphism lemmas which means that the proof can be done simply by expanding the definitions and applying the translation morphisms.

Fast matrix product

In the context we presented, the naïve matrix product (i.e. with cubic complexity) of two matrices M and N can be implemented by transposing the list of lists representing N and then for each i and j compute k M i,k N T j,k :

Definition mulseqmx (M N : seqmatrix) : seqmatrix := let N' := trseqmx N in map (fun r => map (foldl2 (fun z x y => x * y + z) r) N') M. Lemma mulseqmxE (M : 'M[R]_(m,p)) (N : 'M[R]_(p,n)) : mulseqmx (seqmx_of_mx M) (seqmx_of_mx N) = seqmx_of_mx (M *m N).
*m is SSREFLECT's notation for the matrix product. Once again, the rich type information in the quantification of the morphism lemma ensures that it can be applied only if the two matrices have compatible sizes.

In 1969, Strassen [START_REF] Strassen | Gaussian elimination is not optimal[END_REF] showed that 2 × 2 matrices can be multiplied using only 7 multiplications without requiring commutativity. This yields an immediate recursive scheme for the product of two n × n matrices with (n log2 7 ) complexity. 2 This is an important theoretical result, since matrix multiplication was commonly thought to be intrinsically of cubic complexity, it opened the way to many further improvements and gave birth to a fertile branch of algebraic complexity theory.

However, Strassen's result is also still of practical interest since the asymptotically best algorithms known today [START_REF] Coppersmith | Matrix multiplication via arithmetic progressions[END_REF] are slower in practice because of huge hidden constants. Thus, we implemented a variant of this algorithm suggested by Winograd in 1971 [START_REF] Winograd | On multiplication of 2x2 matrices[END_REF], decreasing the required number of additions and subtractions to 15 (instead of 18 in Strassen's original proposal). This choice reflects the implementation of matrix product in most of modern computer algebra systems. A previous formal description of this algorithm has been developed in ACL2 [START_REF] Palomo-Lozano | Certification of matrix multiplication algorithms. strassen's algorithm in ACL2[END_REF], but it is restricted to matrices whose sizes are powers of 2. The extension to arbitrary matrices represents a significant part of our development, which is to the best of our knowledge the first complete formally verified description of Winograd's algorithm.

We define a function expressing a recursion step in Winograd's algorithm. Given two matrices A and B and an operator f representing matrix product, it reformulates the algebraic identities involved in the description of the algorithm:

= B22 -B12 in let C21 := f X Y in let X := A21 + A22 in let Y := B12 -B11 in let C22 := f X Y in let X := X -A11 in let Y := B22 -Y in let C12 := f X Y in let X := A12 -X in let C11 := f X B22 in let X := f A11 B11 in let C12 := X + C12 in let C21 := C12 + C21 in let C12 := C12 + C22 in let C22 := C21 + C22 in let C12 := C12 + C11 in let Y := Y -B21 in let C11 := f A22 Y in let C21 := C21 -C11 in let C11 := f A12 B21 in let C11 := X + C11 in block_mx C11 C12 C21 C22.
This is an implementation of matrix multiplication that is clearly not suited for proving algebraic properties, like associativity. The correctness of this function is expressed by the fact that if f is instantiated by the multiplication of matrices, winograd_step A B should be the product of A and B (=2 denotes extensional equality):

Lemma winograd_stepP (p : positive) (A B : 'M[R]_(p + p)) f : f =2 mulmx -> winograd_step A B f = A *m B.
This proof is made easy by the use of the ring tactic (the script is two lines long). Since version 8.4 of COQ, ring is applicable to non-commutative rings, which has allowed its use in our context.

Note that the above implementation only works for even-sized matrices. This means that the general procedure has to implement a strategy for handling oddsized matrices. Several standard techniques have been proposed, which fall into two categories. Some are static, in the sense that they preprocess the matrices to obtain sizes that are powers of 2. Others are dynamic, meaning that parity is tested at each recursive step. Two standard treatments can be implemented either statically or dynamically: padding and peeling. The first consists of adding rows and/or columns of zeros as required to get even dimensions (or a power of 2), these lines are then simply removed from the result. Peeling on the other hand removes rows or columns when needed, and corrects the result accordingly.

We chose to implement dynamic peeling because it seemed to be the most challenging technique from the formalisation point of view, since the size of matrices involved depend on dynamic information and the post processing of the result is more sophisticated than using padding. Another motivation is that dynamic peeling has shown to give good results in practice.

The function that implements Winograd multiplication with dynamic peeling is called winograd and it is proved correct with respect to the usual matrix product:

Lemma winogradP : forall (n : positive) (M N : 'M[R]_n), winograd M N = M *m N.
The concrete version is called winograd_seqmx and it is also just a direct translation of winograd using only concrete operations on seq based matrices. In the next section, Fig. 2 shows some benchmarks of how well this implementation performs compared to the naïve matrix product, but we will first discuss how to implement concrete algorithms based on dependently typed polynomials.

Polynomials

Polynomials in the SSREFLECT library are represented as records with a list representing the coefficients and a proof that the last of these is nonzero. The library also contains basic operations on this representation like addition and multiplication and proofs that the polynomials form a commutative ring using these operations. The implementation of these operations use big operators [START_REF] Bertot | Canonical big operators[END_REF] which means that it is not possible to compute with them.

To remedy this we have implemented polynomials as lists without any proofs together with executable implementations of the basic operations. It is very easy to build a concrete polynomial from an abstract polynomial, simply apply the record projection (called polyseq) to extract the list from the record. The soundness of concrete polynomials is proved by showing that the pointwise boolean equality on the projected lists reflects the equality on abstract polynomials: Lemma polyseqP p q : reflect (p = q) (polyseq p == polyseq q). Basic operations like addition and multiplication are slightly more complicated to implement for concrete polynomials than for concrete matrices as it is necessary to ensure that these operations preserve the invariant that the last element is nonzero. For instance multiplication is implemented as: Fixpoint mul_seq p q := match p,q with | [::], _ => [::] | _, [::] => [::] | x :: xs,_ => add_seq (scale_seq x q) (mul_seq xs ( %R :: q)) end.

Lemma mul_seqE : {morph polyseq : p q / p * q >-> mul_seq p q}.

Here add_seq is addition of concrete polynomials and scale_seq x q means that every coefficient of q is multiplied by x (both of these are implemented in such a way that the invariant that the last element is nonzero is satisfied). Using this approach we have implemented a substantial part of the SSREFLECT polynomial library, including pseudo-division, using executable polynomials.

Fast polynomial multiplication

The naïve polynomial multiplication algorithm presented in the previous section requires (n 2 ) operations. A more efficient algorithm is Karatsuba's algorithm [START_REF] Abdeljaoued | Méthodes matricielles -Introduction à la complexité algébrique[END_REF][START_REF] Karatsuba | Multiplication of many-digital numbers by automatic computers[END_REF] which is a divide and conquer algorithm based on reducing the number of recursive calls in the multiplication. More precisely, in order to multiply two polynomials written as aX k + b and cX k + d the ordinary method

(aX k + b)(cX k + d) = acX 2k + (ad + bc)X k + cd
requires four multiplications (as the multiplications by X n can be implemented efficiently by padding the list of coefficients by n zeroes). The key observation is that this can be rewritten as

(aX k + b)(cX k + d) = acX 2k + ((a + b)(c + d) -ac -bd)X k + bd
which only requires three multiplication: ac, (a + b)(c + d) and bd. Now if the two polynomials have 2 n coefficients and the splitting is performed in the middle at every point then the algorithm will only require (n log 23 ) which is better than the naïve algorithm. 3 If the polynomials do not have 2 n coefficients it is possible to split the polynomials at for example ⌊n/2⌋ as the formula above holds for any k ∈ and still obtain a faster algorithm. This algorithm has been implemented in COQ previously for binary natural numbers [START_REF] O'connor | Karatsuba's multiplication[END_REF] and for numbers represented by a tree-like structure [START_REF] Grégoire | A purely functional library for modular arithmetic and its application to certifying large prime numbers[END_REF]. But as far as we know, it has never been implemented for polynomials before. When implementing this algorithm we first implemented it using dependently typed polynomials as:

Fixpoint karatsuba_rec (n : nat) p q := match n with | %N => p * q | n'.+1 => if (size p <= 2) || (size q <= 2) then p * q else let m := minn (size p)./2 (size q)./2 in let (p1,p2) := splitp m p in let (q1,q2) := splitp m q in let p1q1 := karatsuba_rec n' p1 q1 in let p2q2 := karatsuba_rec n' p2 q2 in let p12 := p1 + p2 in let q12 := q1 + q2 in let p12q12 := karatsuba_rec n' p12 q12 in p1q1 * 'X^(2 * m) + (p12q12 -p1q1 -p2q2) * 'X^m + p2q2 end.

Here splitp is a function that splits the polynomial at the correct point using take and drop. There is also a wrapper function named karatsuba that calls karatsuba_seq with the greatest degree of p and q. The correctness of this algorithm is expressed by: Lemma karatsubaE : forall p q, karatsuba p q = p * q.

As p and q are SSREFLECT polynomials this lemma can be proved using all of the theory in the library. The next step is to implement the executable version (karatsuba_seq) of this algorithm which is done by changing all the operations in the above version to executable operations on concrete polynomials. The correctness of the concrete algorithm is then proved by: Lemma karatsuba_seqE : {morph polyseq : p q / karatsuba p q >-> karatsuba_seq p q}.

The proof of this is straightforward as all of the operations have morphism lemmas for translating back and forth between the concrete representation and the high-level ones.

In Fig. 2 The benchmarks have been done by computing the square of integer matrices and polynomials using the COQ virtual machine (i.e. by running vm_compute). It is clear that both the implementation of Winograd matrix multiplication and Karatsuba polynomial multiplication is faster than their naïve counterparts, as expected.

gcd of multivariate polynomials

An important feature of modern computer algebra systems is to compute the greatest common divisor (gcd) of multivariate polynomials. The main idea of our implementation is based on the observation that in order to compute the gcd of elements in R[X 1 , . . . , X n ] it suffices to show how to compute the gcd in R[X ] given that it is possible to compute the gcd of elements in R. So, for example, to compute the gcd of elements in [X , Y ] we model it as ( [X ])[Y ], i.e. as univariate polynomials in Y with coefficients in [X ], and then use that there is a gcd algorithm in .

The algorithm that we implemented is based on the presentation of Knuth in [START_REF] Knuth | The art of computer programming, volume 2: seminumerical algorithms[END_REF] which uses that in order to compute the gcd of two multivariate polynomials it is possible to instead consider the task of computing the gcd of primitive polynomials, i.e. polynomials where all coefficients are coprime. Using that any polynomial can be split in a primitive part and a non-primitive part by dividing by the gcd of its coefficients (this is called the content of the polynomial) we get an algorithm for computing the gcd of any two polynomials. Below is our implementation of this algorithm together with explanations of the operations:

Fixpoint gcdp_rec (n : nat) (p q : {poly R}) := let r := modp p q in if r == then q else if n is m.+1 then gcdp_rec m q (pp r) else pp r.

Definition gcdp p q := let (p1,q1) := if size p < size q then (q,p) else (p,q) in let d := (gcdr (gcdsr p1) (gcdsr q1))%:P in d * gcdp_rec (size (pp p1)) (pp p1) (pp q1).

• modp p q computes the remainder after pseudo-dividing p by q.

• pp p computes the primitive part of p by dividing it by its content.

• gcdsr p computes the content of p.

• gcdr (gcdsr p1)(gcdsr q1) computes the gcd (using the operation in the underlying ring) of the content of p1 and the content of q1.

The correctness of this algorithm is now expressed by: Lemma gcdpP : forall p q g, g %| gcdp p q = (g %| p) && (g %| q).

Here p %| q computes whether p divides q or not. As divisibility is reflexive this equality is a compact way of expressing that the function actually computes the gcd of p and q.

Our result is stated in constructive algebra [START_REF] Mines | A Course in Constructive Algebra[END_REF] as: If R is a gcd domain then so is R[X ]. Our algorithmic proof is different (and arguably simpler) than the one in [START_REF] Mines | A Course in Constructive Algebra[END_REF]; for instance, we do not go via the field of fractions of the ring.

As noted in [START_REF] Knuth | The art of computer programming, volume 2: seminumerical algorithms[END_REF], this algorithm may be inefficient when applied on the polynomials over integers. The reference [START_REF] Knuth | The art of computer programming, volume 2: seminumerical algorithms[END_REF] provides a solution in this case, based on subresultants. This would be a further refinement of the algorithm, which would be interesting to explore since subresultants have been already analysed in COQ [START_REF] Mahboubi | Proving formally the implementation of an efficient gcd algorithm for polynomials[END_REF].

The executable version (gcdp_seq) of the algorithm has also been implemented and is linked to the abstract version above by: Lemma gcdp_seqE : {morph polyseq : p q / gcdp p q >-> gcdp_seq p q}.

But when running the concrete implementation there is a quite subtle problem: the polyseq projection links the abstract polynomials with the concrete polynomials of type seq R where R is a ring with a gcd operation. Let us consider multivariate polynomials, for example R[x, y]. In this case the concrete type will be seq (seq R), but seq R is not a ring so our algorithm is not applicable! The next section explains how to resolve this issue so that it is possible to implement computable algorithms of the above kind that rely on the computability of the underlying ring.

Algebraic hierarchy of computable structures

As noted in the previous section there is a problem when implementing multivariate polynomials by iterating the polynomial construction, i.e. by representing

R[X , Y ] as (R[X ])[Y ].
The same problem occurs when considering other structures where the computation relies on the computability of the underlying ring as is the case when computing the characteristic polynomial of a square matrix for instance. For this, one needs to compute with matrices of polynomials which will require a concrete implementation of matrices with coefficients being a concrete implementation of polynomials.

However, both the list based matrices and polynomials have something in common: we can guarantee the correctness of the operations on a subset of the lowlevel structure. This can be used to implement another hierarchy of computable structures corresponding to the SSREFLECT algebraic hierarchy.

Design of the library

We have implemented computable counterparts to the basic structures in this hierarchy, e.g. -modules, rings and fields. These are implemented in the same manner as presented in [START_REF] Garillot | Packaging mathematical structures[END_REF] The type czmodType is the computable -module type parametrized by amodule. The trans function is the translation function from SSREFLECT structures to the computable structures and the only property that is required of it is that it is injective, so we are sure that different high-level objects are mapped to different computable objects.

This way we can implement all the basic operations of the algebraic structures the way we want (for example using fast matrix multiplication as an implementation of *m instead of a naïve one), and the only thing we have to prove is that the implementations behave the same as SSREFLECT's operations on the subset of "wellformed terms" (e.g. for polynomials, lists that do not end with 0). This is done by providing the corresponding morphism lemmas.

The operations presented in the previous sections can then be implemented by having computable structures as the underlying structure instead of dependently typed ones. This way one can prove that polynomials represented as lists is a computable ring by assuming that the coefficients are computable and hence get ring operations that can be applied on multivariate polynomials built by iterating the construction.

It is interesting to note that the equational behavior of an abstract structure is carried as a parameter, but does not appear in its computable counterpart, which depends only on the operations to be implemented. For instance, the same computable ring structure can implement a commutative ring or an arbitrary one, only its parameter varies.

Example: computable ring of polynomials

Let us explain how the list based polynomials can be made a computable ring. First, we define:

Variable R : comRingType. Variable CR : cringType R.

This says that CR is a computable ring parametrized by a commutative ring which makes sense as any commutative ring is a ring. Next we need to implement the translation function from {poly R} to seq CR and prove that this translation is injective: The two polynomials are 2 3 + 2 3 x + 1 2 x y + 1 2 x 2 y and 2 3 + 2 3 y + 1 2 x y + 1 2 x y 2 . The resulting gcd should be interpreted as 1 3 + 1 4 x y.

Conclusions and Further Work

In this paper, we showed how to use high-level libraries to prove properties of algorithms, while retaining good execution capabilities by providing efficient lowlevel implementations. The need of modularity of the executable structure appears naturally and the methodology explained in [START_REF] Garillot | Packaging mathematical structures[END_REF] works quite well. The only thing a user has to provide is a proof of an injectivity lemma stating that the translation behaves correctly. The methodology we suggest has already been used in other contexts, like the CoRN library, where properties of real numbers described in [START_REF] O'connor | Certified exact transcendental real number computation in coq[END_REF] are obtained by proving that these real numbers are isomorphic to an abstract, pre-existing but less efficient version. We tried to show that this approach can be applied in a systematic and modular way.

The library we designed also helps to solve a restriction of SSREFLECT: due to a lot of computations during deduction steps, some of the structures are locked to allow type-checking to be performed in a reasonable amount of time. This locking prevents full-scale reflection on some of the most complex types like big operators, polynomials or matrices. Our implementation restores the ability to perform full-scale reflection on abstract structures, and more generally to compute. For instance, addition of two fully instantiated polynomials cannot be evaluated to its actual numerical result but we can refine it to a computable object that will reduce. This is a first step towards having in the same system definitions of objects on which properties can be proved and some of the usual features of a computer algebra system.

However, in its current state, the inner structure of our library is slightly more rigid than necessary: we create a type for computable -modules, but in practice, all the operations it contains could be packaged independently. Indeed, on each of these operations we prove only a morphism lemma linking it to its abstract counterpart, whereas in usual algebraic structures, expressing properties like distributivity require access to several operations at once. This specificity would make it possible to reorganise the library and create independent structures for each operation, instead of creating one of them for each type. Also, we could use other packaging methods, like type classes [START_REF] Sozeau | First-Class type classes[END_REF], to simplify the layout of the library. However, modifying the library to use type classes on top of SSREFLECT's canonical structures is still on-going work, since we faced some incompatibilities between the different instance resolution mechanisms.

  Matrices are represented by finite functions over pairs of ordinals (the indices): (* 'I_n *) Inductive ordinal (n : nat) : predArgType := Ordinal m of m < n. (* 'M[R]_(m,n) *) Inductive matrix R m n := Matrix of {ffun 'I_m * 'I_n -> R}.

  Definition winograd_step {p : positive} (A B : 'M[R]_(p + p)) f := let A11 := ulsubmx A in let A12 := ursubmx A in let A21 := dlsubmx A in let A22 := drsubmx A in let B11 := ulsubmx B in let B12 := ursubmx B in let B21 := dlsubmx B in let B22 := drsubmx B in let X := A11 -A21 in let Y :

Figure 2 :

 2 Figure 2: Benchmarks of Winograd and Karatsuba multiplication

  using canonical structures. Here are a few examples of the mixins we use: Record trans_struct (A B: Type) : Type := Trans { trans : A -> B; _ : injective trans }. (* Mixin for "Computable" Z-modules *) Record mixin_of (V : zmodType) (T: Type) : Type := Mixin { zero : T; opp : T -> T; add : T -> T -> T; tstruct : trans_struct V T; _ : (trans tstruct) = zero; _ : {morph (trans tstruct) : x / -x >-> opp x}; _ : {morph (trans tstruct) : x y / x + y >-> add x y} }. (* Mixin for "Computable" Rings *) Record mixin_of (R : ringType) (V : czmodType R) : Type := Mixin { one : V; mul : V -> V -> V; _ : (trans V) 1 = one; _ : {morph (trans V) : x y / x * y >-> mul x y} }.

  Definition trans_poly (p : {poly R}) : seq CR := > Eval compute in gcdp_seq p q. = [:: [:: 1 # 3]; [:: ; 1 # 4]]

1

  Documentation available at http://www-sop.inria.fr/members/Maxime.Denes/coqeal/

	Abstract definitions
	Correctness proof
	Algorithmic refinement
	Implementation
	2

log 2 7 is approximately 2.807

log 2 3 is approximately 1.585.

The research leading to these results has received funding from the European Union's 7th Framework Programme under grant agreement nr. 243847 (ForMath).

map (@trans R CR) (polyseq p).

Lemma inj_trans_poly : injective trans_poly.

Assuming that computable polynomials already are an instance of the computable -module structure it is possible to prove that they are computable rings by implementing multiplication (exactly like above) and then prove the corresponding morphism lemmas:

Lemma mul_seqE : {morph trans_poly : p q / p * q >-> mul_seq p q}.

At this point, we could also have used the karatsuba_seq implementation of polynomial multiplication instead of mul_seq since we can prove its correctness using the karatsubaE and karatsuba_seqE lemmas. Finally this can be used to build the CRing mixin and make it a canonical structure.

Examples of computations

This computable ring structure has also been instantiated by the COQ implementation of and which means that they can be used as basis when building multivariate polynomials. To multiply 2 + x y and 1 The result should be interpreted as (2 + 2x) + (3x + x 2 ) y + 3x 2 y 2 + x 3 y 3 . The gcd of 1 + x + (x + x 2 ) y and 1 + (1 + x) y + x y 2 in [x, y] can be computed by:


R Coulon

S Normand

G Ban

V Dumarcher

H P Brau

L Barbot

T Domenech

V Kondrasovs

G Corre

A M Frelin

T Montagu #

T Dautremer #

E Barat

Sodium fast reactor power monitoring using 20 F tagging agent

This work deals with the use of gamma spectrometry to monitor the fourth generation of sodium fast reactor (SFR) power. Simulation part has shown that power monitoring in short response time and with good accuracy is possible measuring liquid sodium delayed gamma emitters produced in-core. An experimental test is under preparation at French SFR Phénix experimental reactor to validate simulation studies. Physical calculations have been done to correlate gamma activity to the released thermal power. Gamma emitter production rate in the reactor core was calculated with technical and nuclear data as sodium velocity, atomic densities, neutron spectra and incident neutron cross-sections of reactions producing gamma emitters. Then, a thermal hydraulic transfer function was used for taking into account primary sodium flow in our calculations and gamma spectra were determined by Monte-Carlo simulations. For power monitoring problematic, use of a short decay period gamma emitter will allowed to obtain a very fast response system without cumulative and flow distortion effects. The experiment will be set during the reactor "end of life testing". The Delayed Neutron Detection (DND) system cell has been chosen as the best available primary sodium sample for gamma power monitoring on Phénix reactor due to its short transit time from reactor core to measurement sample and homogenized sampling in the reactor hot pool. The main gamma spectrometer is composed of a coaxial high purity germanium diode (HPGe) coupled with a transistor reset preamplifier. The signal will be then processed by a digital signal processing system (called Adonis) which is optimum for high count rate and various time activity measurements. To limit statistical problems of the signal, an analytical pileup correction method using duration variable given by our spectrometry system (Adonis) and a nonparametric Bayesian inference for photopeack deconvolution will be used.

I. INTRODUCTION

The main requirement for power monitoring systems is the ability to measure a large range of power with good accuracy and fast response. On SFR, power is usually measured by heat balance on water to vapour heat exchangers which is correlated with the ex-core neutron flux measurements. Heat balance power measurement gives the thermal power with a response time about ten minutes. The ex-core neutron power measurement is done by fission chambers (five on Phénix reactor) with several activity ranges settled at the bottom of the primary vessel to cover the entire neutron flux range of more than 11 decades. This measurement instantly estimates the power but it has some shift problems due to change in chamber fissile nuclei concentration.

Gamma spectrometry systems are not used for this kind of measurement. Although power monitoring using gamma measurement is already used in Pressurised Water Reactors (PWR) by gamma counting measurement using the 16 N tagging agent.

Gamma spectrometry is studying as a power monitoring system for the fourth generation SFR but conventional gamma spectrometry systems are not designed for time varying activity measurement. The Adonis signal processing for gamma spectrometry developed by the CEA-LIST is able to go further through this instrumental limitation [START_REF] Morel | Metrological characterization of the Adonis system used in gammaray spectrometry[END_REF] .

II. METHODS

Gamma activity of each radionuclide is linked to reactor parameters as released thermal power, sodium flow rate in the fast neutron flux and other sodium flow data.

In the first part we will introduce the Adonis signal processing system; the second part the primary sodium gamma source term calculations; the third part deals with Phénix experiment and the fourth part gamma spectra simulation on the Phénix experiment.

II.A. Adonis signal processing for gamma spectrometry

Conventional digital signal processing systems use the same matched filters method as in analogue signal processing. Parameters of this linear filter are adjusted to pre-requirement pulses shape to optimize spectral resolutions and dynamic range. This trapezoidal shaping method gives short range activity abilities. Some pile-up rejection and life time correction methods are settled to increase the dynamic range but they produced shifts on signal measurement above 6x10 5 s -1 of input count rate.

For time varying activity measurement, loss free counting methods are implemented. Count rates are corrected in real time with instant evaluation of the system dead time but metrological performances of this kind of method are intrinsically limited by the dead time estimation accuracy 2 .

Adonis system has a totally different approach for signal processing using a bimodal Kalman smoother [START_REF] Barat | A Bimodal Kalman Smoother for Nuclear Spectrometry[END_REF] . This nonlinear digital signal processing leads us to take into account the stochastic nature of the pulses arrival. An event stream composed, for each pulse of precise value of energy, duration, and separation with previous event is given. This triplet allows an accurate analysis for time varying activities and the dynamic range is the largest without setting on a prejudged input signal. The Adonis electronic is shown on Fig. 1. The Fig. 2 shows 241 Am piled-up spectra done with Adonis and commercial Lynx system which is the most advance DSP spectrometry system of Canberra Company. At 1.3x10 6 s -1 , Adonis gives a less noise spectrum than the Lynx without pile-up rejection using.

Energy (keV)

Count rate (s-1)

Adonis system at 6e6 s-1

Lynx system at 6e6 s-1

Adonis system at 1.3e6 cps Lynx system at 1.3e6 s-1 Fig. 2. Adonis and Lynx 241 Am spectra at high count rate

This system is able to provide on-line gamma spectra for input pulse rates from 0 to 5x10 6 s -1 with still high metrology requirement and optimum spectral resolution (Fig. 3). Furthermore the analytical pileup correction method is only available on Adonis system and is able to remove all spectral the pileup distortion [START_REF] Trigano | Traitement Statistique du Signal Spectrométrique : Etude du Désempilement de Spectre en Energie pour la Spectrométrie Gamma[END_REF] . Input count rate (s-1) Photopeack count rate / Input count rate Adonis system Lynx system Fig. 3. Adonis and Lynx metrological performance at very high input count rate for 59 keV of gamma energy.

So Adonis system will allow a new application field for online gamma spectrometry, as for example reactor monitoring application.

II.B. SFR primary coolant gamma source term calculation

The gamma source term is calculated as a function of power and sodium flow rate in the fast neutron flux. Equations are an adaptation to SFR from PWR equations developed by B. Papin and P. Bernard 4 . Gamma source terms are calculated with neutron reactions introduced in Table I for each activation radionuclide. One should notice that the calculation also took into account dissolved argon (concentration of 0.019 ppm) contained into primary sodium [START_REF] Devaux | Activité de l'Argon 41 dans le Ciel de Pile de Phénix[END_REF] .

First, reaction rate is calculated by discrete nuclear data values as shown in Eq. ( 1):

( ) ( )

j x j E R N E E σ φ = ∑ (1)
Where:

• R j is the reaction rate of the reaction number "j".

• N x is the atomic density of the target nucleus "x" in the homogenized core with non-irradiated fuel (data extracted from technical report [START_REF] Sicard | Données Nécessaires au Calcul Neutronique Simplifié des 3 Premiers Cycles du Coeur Phénix[END_REF] ). • σ j (E) is the partial cross section of the reaction numerber j" on the nucleus "x". Nuclear data evaluations are used and extracted in logarithmic scale with 100 energy groups by decade (see Fig. 4). • Ф(E) is the mean neutron Phenix reactor spectrum computed in 33 energy groups using a MCNP5 Neutron transport simulations with a surface flux tally (F4) and homogeneous configuration (Fig. 5 and Fig. 6). The fast neutron velocity induced a high production of threshold reactions as 23 Na(n,2n) 22 Na, 23 Na(n,p) 23 Ne and 23 Na(n, α) 20 F. Neutron shielding around the core limits ex-core neutron activation (bleu spectrum on Fig. 6). The production rate is then calculated as a function of the thermal power by this formula: . 0.9 . .

act y f fiss c P R X E R V = (2)
Where: • X y is the production rate producing "y" radionuclide.

• P is the instant thermal power released by the core (the residual power is not include in this value). • R act is the reaction rate of the reaction creating "y" as activating isotope. We consider only the fissile zone of the core where 90 % of the thermal power is released and where the homogeneous atomic density was calculated. Therefore, a factor is added to compensate this approximation.

• R fiss is the fission rate on 235 U and 239 Pu nucleus calculated using fission cross-sections (Eq. 1 and Fig. 5). • E f is the mean energy released per fission. • V c is sodium active volume submitted to the fast neutron flux.

Where:

• N y is the activating isotope "y" concentration at reactor core outlet. • X y is the production rate of the activating isotope "y".

• λ y is the decay constant of the activating isotope "y".

• α is an hydraulic parameter which depends of the sodium velocity in active volume. It is calculated as the ratio between sodium flow rate and the sodium active volume of a fuel assembly.

The fuel burn-up is considered as a secondary parameter which changes the fuel enrichment (changing R fiss parameter in Eq. 2) and control rod position (changing V c parameter in Eq.2 and Eq. 3). The sodium temperature is also a second order parameter which induced a nonlinear production rate behaviour (changing R act parameter in Eq. 2).

During the sodium coolant cycle, radionuclides are diluted in primary coolant and could perform a complete coolant cycle if their decay period is long enough. This effect is neglected for short decay radionuclides as 20 F and 23 Ne but it is predominant for 24 Na and 22 Na. To take into account this cumulative effect, activating isotope concentrations are calculated by adding a cumulative correction as shown in Eq. 4.

( ) _ ' 1 0 1 y n t y y c n N N t e dt λ δ ∞ ∞ =   = +     ∑ ∫ (4)
Where:

• N' y is the cumulative concentration of the radionuclide "y" at the outlet of fuel assembly. • n is the coolant cycle number.

• δ c (t) is concentration distribution during a coolant cycle.

This cumulative effect increases time to estimate a power value for long period tagging agent. Indeed a short decay period gamma emitter is better to delete this effect and to immediately obtain the power measurement. We can notice that sodium flow in the hot pool is very heterogeneous due to the control plug. Sodium has a complex behavior due to the mixture of mixed jet effect, diffusion effect and dilution effect. One of the six responses shows a high jet effect on Fig. 11 which is important to obtain short decay period radionuclides in DND sodium sample. Flow heterogeneities will be corrected by the 6 sodium sample mixture measurement which is essential to obtain global power estimation.

A tridimensional thermal hydraulic simulation project using Trio code is under consideration in collaboration with the thermal hydraulic department of the CEA/DEN to improve the sodium flow distribution knowledge inside Phenix primary hot pool [START_REF] Bieder | Qualification of the CFD Code Trio U for Full Scale Reactor Applications[END_REF] . By convolving a more accurate dilution distribution with radionuclides decay probability, a more precise thermal hydraulic transfer function will be determined.

Thanks to the tagging agent concentration in the DND sample, the incident gamma spectrum is simulated by MCNP5 using a flux point tally (F5) and reduction variance technique due to shield thickness (17 cm of lead around the sodium sample). Fig. 12 is the MCNP model of experiment. • S e (E) is the incident spectrum on sensor surface due to source photons of energy "e"(listed in Table. 1). • λ y is the decay constant of the radionuclide "y"

producing "e" energy photons.

• η e is the branching ratio of "e" energy photons produced by the radionuclide "y". • N y " is the "y" radionuclide concentration in the sodium sample calculated in Eq. 5.

• s is the sensor input surface.

• Ф e ∆ (E) is the MCNP5 simulation response for the source photons of "e" energy emitted in the sampling part "∆" (Fig. 13).

• V ∆ is the volume of the sampling part "∆".

To simulate the sensor response, MCNP5 is used with a pulse high tally (F8) in photon-electron mode. Deposited energy during each photon interaction in the sensor is then calculated and finally gives energy deposited spectra by convolving by incident spectra (Eq. 9).

( ) (

)

( ) e e E S E S E E e dE = Η - ∑ ∫ (9)
Where:

• S(E) is the response spectrum without noise.

• S e (E) is the incident spectrum on sensor due to source photons of energy "e". • H(E) is efficiency spectra of sensor.

The electronic noise is taken into account by a Gaussian noise convolution on the previews spectrum (Eq. 10).

( )

( ) ( ) ( ) 2 2 2 , '( ) 2
,

CR E T CR e S E S d T ε σ ε ε ε ε πσ ε - - = ∫ (10)
Where:

• S'(E) is the noisy response spectrum.

• S(E) is the deposited energy spectrum.

• σ(E,T CR ) is the resolution parameter of the Gaussian which depends of the pulse energy "E" and the total input count rate "T CR " (The resolution evolution as a function of the total input count rate is an Adonis particularity).

In Adonis, piled-up pulses are never rejected. Pile-up contribution is then finally simulated by a recursive calculation using pile-up probability. A summation is done

Proceedings of Global 2009 Paris, France, September 6-11, 2009 Paper 9416 when a pulse arrives during treatment of a preview pulse (about 350 ns).

III. RESULTS

This study links reactor power parameters to primary sodium tagging agent concentration which principally depends on the instant released power and the sodium flow rate in the neutron flux. First part will introduced the Phénix experiment simulation results. Second part will introduced simulation results of an optimum SFR power measurement system.

III.A. Simulation results for the Phénix experiment

Experimental Phénix Spectra simulations show the possibility to measure 24 Na, and 20 F tagging agents during our experiment. Our study is then focused on those radionuclides to measure reactor power because 23 Ne exhibits too low photon energy of 440 keV regarding to the DND shielding thickness. The LaBr3:Ce sensor could be settled next to a sodium sample without any lead shielding (see in Fig. 12). The direct to diffuse photon ratio will be better than in the HPGe diode case.

Sodium thermal hydraulic behaviour is the most influent phenomena on the calculation. Simulations are made for 2 extremes cases where the jet effect in the outlet of the assembly is maximal (case 1 in Table 2) and where this effect is minimal (case 2 in Table 2). This simulation shows the possibility to measure the 20 F gamma activity in the Phénix experiment.

During this experiment, reactor power will be increased by scale of 100 MW to compare at each scale the gamma signal to heat balance power measurement. In a first time, a relative relation between 20 F activity and reactor power will be obtained. In a second time, this simulation results will be used as numerical calibration to estimate 20 F nuclei concentration after acquisitions and link it to the power by the Eq. 3.

III.B. Simulation results of an optimum SFR power measurement system

The tagging agent core production is directly link to the instant released power but the sodium flow to the measurement cell induced interferences. Sampling has to be done close to the core outlet to limit flow heterogeneity phenomena. Gamma emitters should also have a short decay period (as 20 F and 23 Ne) to limit flow accumulation phenomena.

To design a power monitoring system measuring short decay period activation products, the transit time from core outlet to measurement sample has to be set as short as possible (Fig. 16). This calculation shows that a short time response power monitoring system using 20 F or 23 Ne tagging agent is easily achievable if the sodium sampling volume and the transit time to measurement sample can be designed. Fig. 17 is a spectrum simulated with 2 s of transit time and 1.5x10 6 s -1 of total incident input count rate. An Adonis system could cover 3 decades of SFR power range. Five systems could cover the entire power range measuring 23 Ne, 20 F and 24 Na gamma activities.

IV. DISCUSSION

This simulation study shows the possibility of 20 F measurement but statistical problem due to pileup effects or too low 20 F detection limits could be met.

Post-processing softwares are developed at CEA-LIST to improve statistics for gamma spectrometry. A new nonparametric Bayesian method is currently designed and used to extract the 20 F signal [START_REF] Barat | Nonparametric Bayesian Inference in Nuclear Spectrometry[END_REF] . In case of piled-up spectra, we will use an analytical pileup correction method using the duration variable given by Adonis 10 .

On-line gamma spectrometry measurement will allow several activation products analysis. This system will permit to design a more autonomous power measurement system with temperature and flow rate correction done only by the gamma spectrometric signal.

IV. CONCLUSIONS

This study deals with a new method for power measurement based on sodium activation calculations; thermal hydraulic and gamma photons transport simulations to link gamma signal of primary sodium tagging agents to the released in-core reactor power.

Simulation work has shown that Adonis spectrometry system could be settled as generation IV SFR power monitoring systems matching a large range of power measurement in a short response time.

The Phénix experiment is currently under preparation to validate the method and determine the experimental accuracy of this power measurement system.

Fig. 1 .

 1 Fig. 1. Adonis electronic

Fig. 4 .

 4 Fig. 4. JEFF 3.1.1 Activation cross-sections of 23 Na

Fig. 5 .Fig. 6 .

 56 Fig. 5. MCNP5 schematic Phénix primary vessel model

Fig. 7 .

 7 Fig. 7. JEFF 3.1.1 Fission cross-sections on 239 Pu and 235 U nuclei

Fig. 11 .

 11 Fig. 11. Typical concentration distribution functions for one fuel assembly outlet

Fig. 12 .

 12 Fig. 12. MCNP model -Phénix experimental test

Fig. 19 .

 19 Fig. 19. Simulated mixed gamma spectra of the HPGe signal.

Fig. 16 .

 16 Fig. 16. Tagging agent concentration as a function of transit time to core outlet

Fig. 17 .

 17 Fig. 17. Simulated mixed gamma spectrum designed for power monitoring system.

TABLE I Neutron

 I

			reactions in primary sodium
	Nuclear	Induced	Decay	Photon	Branching
	reaction	emitter	period	energy	ratio
			(s)	(keV)
	23 Na(n,2n)	22 Na	8.22x10 7	1274.5	1
				1368.63	1
				2754.03	1
	23 Na(n,γ)	24 Na	5.39x10 4	2870	2.0x10 -6
				3867.2	5.2x10 -4
				4238	7.5x10 -6
	23 Na(n,γ)	24m Na	2.02x10 -2	472.0	1
	23 Na(n,α)	20 F	1.10x10 1	1633.7 3334.4	1 1.4x10 -4
				439.9	3.3x10 -1
				1636.4	1.0x10 -2
	23 Na(n,p)	23 Ne	3.72x10 1	2076.3	1.0x10 -3
				2542.3	2.7x10 -4
				2982.2	3.8x10 -4
				1293.6	1
	40 Ar(n,γ)	41 Ar	6.58x10 3	1677.0	5.2x10 -4
				3967.0	1.5x10 -5

TABLE 2 20

 2 F photopeak signal estimation for the Phénix experiment

					Time to
	Detector	case	Signal 20 F (s -1 )	Noise (keV -1 .s -1 )	exceed detection
					limit (s)
		1	1.75	0.70	37.2
	HPGe
		2	0.73	0.56	71.4
		1	1873	318	221
	LaBr3
		2	643	270	547

TABLE 3

 3 Power dynamic range computed for an Adonis spectrometer

	Tagging agents	23 Ne	22 Na	20 F	24 Na
	Energy (keV)	440	1275	1334	2754
	Decay period (s)	3.72x10 1 8.22x10 7	1.1x10 1	5.39x10 4
	Signal (s-1)	2.97x10 4 6.00x10 2 4.14x10 4 2.39x10 4
	Noise (keV-1.s-1)	3.5x10 2	4.28x10 2 1.97x10 2 1.06x10 2
	Dynamic range	1.2x10 3	1.8x10 1	1.8x10 3	1.3x10 3

ACKNOWLEDGMENTS

The authors would like to thank CEA-DRT and CEA-DAM for financial support.

Proceedings of Global 2009

Paris, France, September 6-11, 2009 Paper 9416

II.C. Phénix experiment

The experimental system will be set during the Phénix reactor "end of life testing". The Delayed Neutron Detection (DND) cell has been chosen as the best available location on this Phénix reactor (Fig. 9). This DND system has been set in order to detect clad failure by delayed neutron emitter measurement. This location is optimum for global power measurement due to short sodium transit time (about 30 s) and homogenized measurement done by 6 points in hot pool near each intermediary heat exchanger. Then, short decay radionuclides measurement and thermal power correlation is potentially achievable.

The main spectrometer is composed of a coaxial high purity germanium (HPGe) diode coupled with a reset transistor preamplifier and a hybrid cryostat (electrical LN2 Cryo-cycle) which is convenient for deported measurement. The HPGe diode signal will be processed by the Adonis signal processing.

The secondary spectrometer is a LaBr3:Ce scintillator coupled with photoelectric electron-multiplier tube. A new numerical electronic coupled with a digital triangular filter will be also implemented. LaBr3:Ce scintillators have the highest resolution of scintillators but it still gives about 2.2 % energy resolution where HPGe diode gives 0.1 % of resolution at the 1.17 MeV 60 Co peak. Despite of its resolution, it will give another measurement of the primary sodium sample with a more gamma dose robust detector.

Detectors are shown in Fig. 10.


L Brasco
email: brasco@cmi.univ-mrs.fr

G Carlier
email: carlier@ceremade.dauphine.fr

Congested traffic equilibria and degenerate anisotropic PDEs

Keywords: traffic congestion, anisotropic and degenerate PDEs, regularity

 to the anisotropic case.

Introduction

Congested traffic equilibrium models on finite networks have received a lot of attention since the early 50's because of applications to road traffic as well as to communication networks. In such problems, the notion of Wardrop equilibrium [START_REF] Wardrop | Some theoretical aspects of road traffic research[END_REF] plays a distinguished role. The standard congested network model consists of:

• a finite oriented connected graph G = {N, E} (N is the set of nodes of G and E the set of its edges) modelling the network;

• distribution of sources f G -:= x∈N f G -(x) δ x and sinks f G + := y∈N f G + (y) δ y which are discrete measures with same total mass (which can be assumed to be 1) on the set of nodes N that capture respectively the amount of traffic emitted and absorbed at the nodes of the network;

• for each edge e ∈ E, there is a function h e : R + → R + that is typically increasing and which captures the congestion effect in the sense that h e (m) represents the travelling time of edge e when the traffic flow on e is m.

One then looks for an equilibrium flow configuration, that is a flow configuration that is compatible with the distribution of sources and sinks and such that paths on which there is a positive flow minimize the travelling time between their endpoints, given the congestion effects. More precisely, a flow configuration is a probability measure Q on the set C G of (loop-free) paths on the network that satisfies the mass conservation constraints

Q {γ ∈ C G : γ starts at x} = f G -(x), for every x ∈ N, Q {γ ∈ C G : γ ends at y} = f G + (y),
for every y ∈ N.

A flow configuration Q induces a collection of arc-flows through m Q (e) := Q {γ ∈ C G : e ∈ γ} , for every e ∈ E, and then a collection of arc-travelling times h e (m Q (e)) as well as the total travelling time of a path γ ∈ C G :

T Q (γ) = e∈γ h e (m Q (e)).
A Wardrop equilibrium is then by definition a flow configuration Q such that for every x and y and every γ in the set of paths C G x,y starting at x and ending at y, if Q({γ}) > 0 then one shoud have

T Q (γ) = min τ ∈C G x,y T Q (τ ).
In other words, a Wardrop equilibrium requires that users behave rationally by choosing the shortest available paths, taking congestion into account, i.e. the fact that travel times increase with the flow. Finding Wardrop equilibria is a fixedpoint problem in nature that presents some analogies with the mean-field games theory of Lasry and Lions [START_REF] Lasry | Mean-Field Games, Japan[END_REF], eventhough it is purely stationary. Soon after the work of Wardrop, it was observed by Beckmann, McGuire and Winsten [START_REF] Beckmann | Studies in Economics of Transportation[END_REF] that Wardrop equilibria coincide with the minimizers of the convex functional

e∈E H e (m Q (e))
where

H e (m) = m 0 h e (s) ds,
among flow configurations Q. This variational characterization is nice both from a theoretical point and for numerical computations. Unfortunately, the minimization problem above has one flow variable per admissible path on the network, it may therefore quickly become untractable for realistic road or communication networks.

A natural question at this point is whether the situation somehow simplifies on very dense networks and if one can find some continuous counterpart to the discrete theory. This issue was recently addressed in [START_REF] Baillon | From discrete to continuous Wardrop equilibria[END_REF] in the case of a two-dimensional cartesian grid, where

G = G ε = {N ε , E ε } with N ε = εZ 2 ∩ Ω and E ε = {[x, x + εv k ] : x ∈ N ε , k ∈ {1, . . . , 4}} ,
where (v 1 , v 2 , v 3 , v 4 ) := ((1, 0), (0, 1), (-1, 0), (0, -1)) denote the directions of the grid and Ω is some bounded two-dimensional domain. Under suitable assumptions1 , it is shown in [START_REF] Baillon | From discrete to continuous Wardrop equilibria[END_REF] that Wardrop equilibria at scale ε converge as ε → 0 + to solutions of inf

Q∈Q(f + ,f -) Ω 4 k=1 H k (x, m Q k (x)) dx, (1.1)
where

Q(f + , f -) is the set of Borel probability measures Q on C([0, 1], R 2 ) satis- fying the mass conservation conditions e 0# Q = f -, e 1# Q = f + , where e t (γ) = γ(t), t ∈ [0, 1], (1.2)
and such that

Q(L) = 1. Here L = W 1,∞ ([0, 1], Ω) is the set of Lipschitz curves. For k = 1, . . . , 4, the measure m Q k is defined by: Ω ϕ(x) dm Q k (x) := L 1 0 ϕ(γ(t)) ( γ(t) • v k ) + dt dQ(γ), (1.3)
for every ϕ ∈ C(Ω, R), while the functions H k are related to h k by H k (x, m) = m 0 h k (x, s) ds. They are therefore convex and nondecreasing in their second argument, since h k are typically nonnegative and nondecreasing in their second argument. We shall suppose in the sequel that there is some p > 1 such that for k = 1, . . . , 4, the functions H k satisfies the growth condition

1 λ (m p -1) ≤ H k (x, m) ≤ λ (m p + 1). (1.4) So that in (1.1), it is intended that m Q k ∈ L p (Ω), k = 1, . . . , 4
and thus, we shall also assume that

Q p (f + , f -) := Q ∈ Q(f + , f -) : m Q k ∈ L p (Ω), k = 1, . . . , 4 = ∅. (1.5)
This assumption is satisfied for instance as soon as f + and f -belong to L p (Ω) and Ω is convex, as follows from the results of De Pascale, Evans and Pratelli [START_REF] Pascale | Integral estimates for transport densities[END_REF][START_REF] Pascale | Regularity properties for Monge transport density and for solutions of some shape optimization problem[END_REF][START_REF] Pascale | Sharp summability for Monge Transport density via Interpolation[END_REF] and Santambrogio [START_REF] Santambrogio | Absolute continuity and summability of transport densities: simpler proofs and new estimates[END_REF].

Let us note that (1.1) is very similar to the optimal transport with congestion studied in [START_REF] Carlier | Optimal transportation with traffic congestion and Wardrop equilibria[END_REF], except that it keeps track of the anisotropy of the network. However, at first glance, it is absolutely not clear that problem (1.1) which is posed over probability measures over curves (two layers of infinite-dimensions!) offers any simplication with respect to the discrete problem. In the present paper, we want to go one step further with respect to the convergence results of [START_REF] Baillon | From discrete to continuous Wardrop equilibria[END_REF] and we shall see that solving (1.1) roughly amounts to solve a single PDE. Following a similar approach as the one developed in [START_REF] Brasco | Congested traffic dynamics, weak flows and very degenerate elliptic equations[END_REF] for the isotropic case, we shall investigate the relationship between (1.1) and the simpler problem

inf σ∈L p (Ω,R 2 ) Ω 4 k=1 H k (x, (σ(x) • v k ) + ) dx : -div σ = f , (1.6)
where f := f + -f -and the equation -div (σ) = f should be understood as

Ω ∇u • σ = Ω u df, for every u ∈ C 1 Ω ,
so that it incorporates in the weak sense the homogeneous Neumann boundary condition σ • ν Ω = 0 on ∂Ω. Problem (1.6) belongs to a class of problems introduced by Beckmann [START_REF] Beckmann | Studies in Economics of Transportation[END_REF] for the design of an efficient commodity transport program2 . Solving (1.6) can be done by first solving the Euler-Lagrange equation of its dual formulation and then by using the primal-dual optimality conditions. However, in typical congestion models, the functions H k (x, •) have a positive derivative at zero. Indeed, recall that this derivative is h k (x, 0) and it should be positive since one cannot go at infinite speed even when there is no congestion. This creates a singularity in the integrand in (1.6) which, in turn, makes the Euler-Lagrange equation of the dual extremely degenerate. Because of the anisotropic feature of the problem, the kind of degeneracies we are facing is even worse than the one in [START_REF] Brasco | Congested traffic dynamics, weak flows and very degenerate elliptic equations[END_REF]. Indeed, the prototypical equation of [START_REF] Brasco | Congested traffic dynamics, weak flows and very degenerate elliptic equations[END_REF] was

-div (|∇u| -1) q-1 + ∇u |∇u| = f, with q = p = p/(p -1)
, whereas here we shall rather deal with anisotropic equations of the form

- 2 i=1 ∂ i (|∂ i u| -δ i ) q-1 + ∂ i u |∂ i u| = f,
which, in contrast with the isotropic case, degenerates in an unbounded set of values of the gradient. Even in the less degenerate case where all the δ i 's are zero, the previous equation (refered to as the pseudo q-Laplacian equation in [START_REF] Belloni | The pseudo p-Laplace eigenvalue problem and viscosity solutions as p → ∞[END_REF]) is more delicate than its isotropic counterpart, which has been much more studied and for which more regularity results are available.

The paper is organized as follows. In Section 2, we investigate some relationship between (1.1) and (1.6) in a similar way as in [START_REF] Brasco | Congested traffic dynamics, weak flows and very degenerate elliptic equations[END_REF]. In Section 3, we formulate the optimality conditions for (1.6) in terms of the solutions of the dual and emphasize the kind of PDEs realistic anisotropic models of congestion lead to. Section 4 is devoted to some regularity results for such degenerate anisotropic elliptic PDEs. In this final section, we shall work in any dimension d.

Relationship with Beckmann problem

We consider again problems (1.1) and (1.6), our aim is to emphasize some precise connections between those two problems. We assume that Ω is an open bounded connected subset of R 2 with a Lipschitz boundary, that (1.5) holds and that there is a neighbourhood U of Ω such that each function H k is Carathéodory on U × R + , convex nondecreasing in its second argument with H k (x, 0) = 0 a.e. x ∈ U and satisfies for some λ > 1 the growth condition (1.4), for a.e. x ∈ U and every m ∈ R + .

For Q ∈ Q p (f + , f -), let us define the vector-measure σ Q by:

Ω F dσ Q := L 1 0 F (γ(t)) • γ(t) dt dQ(γ), for every F ∈ C(Ω, R 2 ), (2.1)
where we recall that L = W 1,∞ ([0, 1], Ω) has full mass for Q. For u ∈ C 1 (Ω), we then have by the very definition of

Q(f + , f -) Ω ∇u • dσ Q = L u(γ(1)) -u(γ(0)) dQ(γ) = Ω u df, i.e. -div σ Q = f := f + -f -. Now for ϕ ∈ C(Ω) with ϕ ≥ 0, let us observe that Ω ϕ d(σ Q • v k ) = L 1 0 ϕ(γ(t)) γ(t) • v k dt dQ(γ) ≤ L 1 0 ϕ(γ(t)) ( γ(t) • v k ) + dt dQ(γ) = Ω ϕ dm Q k , i.e. for k = 1, . . . , 4, one has σ Q • v k ≤ m Q k in the sense of measures. Since Q ∈ Q p (f + , f -), this implies that σ Q ∈ L p (Ω, R 2 ) and thus (σ Q • v k ) + ≤ m Q k a.e. in Ω.
(2.2)

Note then that the assumption (1.5) asserting that Q p (f + , f -) is not empty implies that -div σ = f can be solved in L p (Ω, R 2 ), so that the zero-mass signed measure f belongs to the dual of W 1,q (Ω), where q = p/(p -1) is the conjugate exponent of p. The existence of solutions to (1.6) directly follows; as for the existence of solutions to (1.1), it follows from the same arguments as in [START_REF] Baillon | From discrete to continuous Wardrop equilibria[END_REF] or [START_REF] Carlier | Optimal transportation with traffic congestion and Wardrop equilibria[END_REF].

Theorem 2.1. Under the assumptions above, we have

inf (1.1) = inf (1.6). Proof. Let Q ∈ Q p (f + , f -), since H k is nondecreasing in its second argument, we deduce from (2.2) that Ω 4 k=1 H k (x, m Q k ) dx ≥ Ω 4 k=1 H k (x, (σ • v k ) + ) dx,
and since -div

(σ Q ) = f := f + -f -, we immediately get inf (1.1) ≥ inf (1.6).
For the converse inequality, we shall use Moser's flow argument (see [START_REF] Dacorogna | On a partial differential equation involving the Jacobian determinant[END_REF][START_REF] Moser | On the volume elements on a manifold[END_REF]), together with some standard regularization as follows. Let σ be a solution of (1.6), extend it by 0 outside Ω, let then ρ ∈ C ∞ c (R 2 ) be a positive function, supported in the unit ball B 1 and such that

R 2 ρ = 1. For ε 1 so that Ω ε := Ω + ε B 1 U , we set ρ ε (x) := ε -2 ρ(ε -1 x) and σ ε := ρ ε * σ. By construction, we thus have that σ ε ∈ C ∞ (Ω ε ) and -div σ ε = f ε + -f ε - in Ω ε , σ ε = 0, on ∂Ω ε , (2.3)
where

f ε ± = ρ ε * (f ± 1 Ω ) + ε. Now comes Moser's flow construction: for t ∈ [0, 1] and x ∈ Ω ε , we set g ε (t, x) := (1 -t) f ε -(x) + t f ε + (x)
, since σ ε and g ε are smooth and the latter is bounded from below by ε > 0, one may define the flow X ε of the vector field

v ε := σ ε /g ε , i.e. Ẋε t (x) = v ε (t, X ε t (x)) = σ ε (X ε t (x)) (1 -t) f ε -(X ε t (x)) + t f ε + (X ε t (x)) , X ε 0 (x) = x, (t, x) ∈ [0, 1] × Ω ε .
Since by construction g ε satisfies the continuity equation

∂ t g ε + div (g ε v ε ) = 0, with initial datum g ε (0, •) = f ε -and v ε is smooth, we have X ε t # f ε -= g ε (t, •) for every t ∈ [0, 1]. Let then L ε := W 1,∞ ([0, 1], Ω ε ) and consider the measure Q ε on L ε defined by Q ε := Ωε δ X ε • (x) df ε -(x).
We then have e

t# Q ε = X ε t # f ε -= g ε (t, •) for every t ∈ [0, 1] and in particular e 0# Q ε = f ε - and e 1# Q ε = f ε + . (2.4)
Now let us define σ Q ε and m Q ε k as in (2.1) and (1.3) respectively, by using testfunctions that are now defined on Ω ε . For every ϕ ∈ C(Ω ε ), by the definition of Q ε , the fact that X ε t # f ε -= g ε (t, •) and that v ε g ε = σ ε and by using Fubini's formula, we then have

Ωε ϕ dm Q ε k = Ωε 1 0 ϕ(X ε t (x)) (v ε (t, X ε t (x)) • v k ) + dt df ε -(x) = 1 0 Ωε ϕ(y) (v ε (t, y) • v k ) + d(X ε t # f ε -)(y) dt = 1 0 Ωε ϕ(y) (v ε (t, y) • v k ) + g ε (t, y) dy dt = Ωε ϕ(y) (σ ε (y) • v k ) + dy, so that m Q ε k = (σ ε •v k ) + .
A similar computation (with a vector valued test-function) similarly shows that σ Q ε = σ ε . Now, let us remark that the definitions of m Q k and σ Q are invariant by arbitrary injective reparameterization of curves. In particular, if for any Lipischitz curve γ, we denote by γ its constant speed reparametrization and if we denote by Q the push forward of Q through the map γ → γ, we have

m Q k = m Q k and σ Q = σ Q . Arguing as in [8, Lemma 2.8], the uniform L p bound on σ Q ε = σ Q ε = σ ε implies that the family of Borel measures Q ε on C([0, 1], R 2
) is tight and thus admits a (not relabeled) subsequence that * -weakly converges to some measure Q. Observe that this limit Q is a probability measure, since the total mass of Q ε equals that of f ε + i.e. 1 + ε and, arguing again as in [8, Lemma 2.8], one can show that Q(L) = 1 (more precisely for fixed ε > 0, one proves exactly as in [START_REF] Carlier | Optimal transportation with traffic congestion and Wardrop equilibria[END_REF]Lemma 2.8

] that Q(L ε ) = 1 so that Q(L) = lim ε↓0 + Q(L ε ) = 1). Moreover, thanks to the * -weak convergence of Q ε to Q, passing to the limit in (2.4) gives Q ∈ Q(f + , f -). We finally remark that m Q ε k = m Q ε k = (σ ε • v k ) + converges strongly in L p to (σ • v k ) + and
then, thanks to the same semicontinuity argument as in [8, Lemma 2.9], we have

m Q k ≤ (σ • v k ) + in the sense of measures. This implies that m Q k ∈ L p i.e. Q ∈ Q p (f + , f -). Using the monotonicity of H k , we obtain Ω 4 k=1 H k (x, m Q k (x))dx ≤ Ω 4 k=1 H k (x, (σ(x) • v k ) + )dx = inf(1.6)
and thus we can finally infer inf (1.1) ≤ inf (1.6).

Remark 2.2. The convexity assumption on the H k 's is essential for the existence of minimizers, but it played no role in the proof that inf (1.1) = inf (1.6) where only monotonicity matters.

Remark 2.3. Note that with the previous proof, we have in fact also shown slightly more precise results. Firstly

Q solves (1.1) =⇒ σ Q solves (1.6).
Secondly, if in addition each function H k is increasing in its second argument, it easily follows from inequality (2.2) that we have the equivalence

Q solves (1.1) ⇐⇒ σ Q solves (1.6) and m Q k = (σ Q • v k ) + , k = 1, . . . , 4.
Lastly, we have also shown how to build a minimizing sequence for (1.1) from a regularization of a solution σ of (1.6) thanks to Moser's flow method. Also notice that if σ and f + and f -are regular enough, so that the flow can be defined (in the Ambrosio-Di Perna-Lions sense [START_REF] Ambrosio | Transport equation and Cauchy problem for BV vector fields[END_REF][START_REF] Diperna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF], for instance), there is no need to do such a regularization. By the way, the regularity of σ will be addressed in section 4.

Link with anisotropic elliptic PDEs

To shorten notations, let us set

H(x, σ) := 4 k=1 H k (x, (σ • v k ) + ), (3.1)
which can also be written in the separable form H(x, σ) = H 1 (x, σ 1 ) + H 2 (x, σ 2 ), where

H 1 (x, σ 1 ) = H 1 (x, σ 1+ ) + H 3 (x, σ 1-), H 2 (x, σ 2 ) = H 2 (x, σ 2+ ) + H 4 (x, σ 2-).
The dual problem of (1.6) then reads sup

u∈W 1,q (Ω) Ω u df - Ω H * (x, ∇u) dx , (3.2)
and since f = f + -f -has zero mean, we may as well restrict the optimization in (3.2) to zero-mean W 1,q (Ω) functions. Here q = p = p/(p -1). As usual, H * (x, •) denotes the Legendre transform of H(x, •) and it is separable as well

H * (x, z) = H * 1 (x, z 1 ) + H * 2 (x, z 2
). Let us recall that in congestion models the functions H k (x, •) are the primitives of functions h k (x, •), where h k (x, m) represents the time per unit of length to move at x in the direction v k when the intensity of traffic in this direction is m. The functions h k therefore represent a sort of anisotropic metric that captures the congestion effect through some relationship between time per unit of length and mass per unit of length. Having this interpretation in mind, it is reasonable to assume that h k (x, •) are continuous increasing functions on R + which satisfy some p -1 growth condition (so that (1.4) holds) and that h k (x, 0) > 0 (the metric is everywhere positive even when there is no traffic). This implies that H is strictly convex in its second variable (hence, by duality H * is C 1 ) but not differentiable on the axes of coordinates, more precisely, the subdifferentials of H 1 (x, .) and H 2 (x, .) are intervals having 0 in their interior. The Legendre transforms H * 1 (x, .) and H * 2 (x, .) therefore vanish on such intervals which makes the Euler-Lagrange equation of (3.2) degenerate close to the axes of coordinates.

Recall that by standard convex duality (see [START_REF] Ekeland | Convex analysis and Variational problems[END_REF] for example), the values min (1.6) and max (3.2) coincide and the primal-dual optimality conditions characterize the minimizer (unique by strict convexity) σ of (1.6) by

σ(x) = ∇H * (x, ∇u(x)), a.e. x ∈ Ω
where u is a solution of (3.2). This is equivalent to the requirement that u is a weak solution of the Euler-Lagrange equation:

-div ∇H * (x, ∇u) = f, in Ω, ∇H * (x, ∇u) • ν Ω = 0, on ∂Ω, (3.3)
in the sense that

Ω ∇H * (x, ∇u(x)) • ∇ϕ(x) dx = Ω ϕ(x) df (x), for every ϕ ∈ W 1,q (Ω).
Let us note that eventhough u is non unique, σ is.

A prototypical example is h k (x, m) = a k (x) m p-1 + δ k with δ k > 0 and the weights a k are smooth and bounded away from zero. A direct computation then gives

H * 1 (x, z 1 ) = b 1 (x) q (z 1 -δ 1 ) q + + b 3 (x) q (-z 1 -δ 3 ) q + , and
H * 2 (x, z 2 ) = b 2 (x) q (z 2 -δ 2 ) q + + b 4 (x) q (-z 2 -δ 4 ) q + ,
where b k = a

-1 p-1 k
. In this case (3.3) takes the form

-∂ 1 b 1 (∂ 1 u -δ 1 ) q-1 + -b 3 (-∂ 1 u -δ 3 ) q-1 + -∂ 2 b 2 (∂ 2 u -δ 2 ) q-1 + -b 4 (-∂ 2 u -δ 4 ) q-1 + = f. (3.4)
Observe that in the symmetric homogeneous case, i.e. when δ 1 = δ 3 , δ 2 = δ 4 and a k ≡ 1, the previous simplifies to

-∂ 1 (|∂ 1 u| -δ 1 ) q-1 + ∂ 1 u |∂ 1 u| -∂ 2 (|∂ 2 u| -δ 2 ) q-1 + ∂ 2 u |∂ 2 u| = f.
Since H * 1 (x, z 1 ) and H * 2 (x, z 2 ) vanish whenever z 1 ∈ [-δ 3 , δ 1 ] and z 2 ∈ [-δ 4 , δ 2 ] respectively, any u whose gradient belongs to the rectangle [-δ 3 , δ 1 ] × [-δ 4 , δ 2 ] solves the previous equation with f = 0. Hence, there is no hope to recover estimates on the second derivatives of u or even oscillation estimates on ∇u from (3.4). The best one can hope for is that u is Lipschitz. However, we shall see in the next section how to obtain some regularity results directly for the vector field σ = (σ 1 , σ 2 ) that solves (1.6), i.e.

σ 1 = b 1 (∂ 1 u -δ 1 ) q-1 + -b 3 (-∂ 1 u -δ 3 ) q-1 + , σ 2 = b 2 (∂ 2 u -δ 2 ) q-1 + -b 4 (-∂ 2 u -δ 4 ) q-1 + .
Typical traffic congestion problems are two-dimensional, however degenerate equations of the form (3.4) may arise in other contexts, as in the relaxation of some nonconvex variational problems. For this reason, in the next section we will actually work in arbitrary dimension d.

Regularity 4.1 A general result

Let q ≥ 2, Ω be an open bounded subset of R d , F : Ω×R d → R d and G : R d → R d be some continuous vector fields such that there exist ν > 0 and µ ≥ 0 such that for every (x, z) and (y, w) in Ω × R d , one has

|F (x, z)| ≤ µ|z| q-1 , (4.1)
(F (x, z) -F (x, w)) • (z -w) ≥ ν |G(z) -G(w)| 2 , (4.2)
|F (x, z) -F (x, w)| ≤ µ |G(z)| q-2 q + |G(w)| q-2 q |G(z) -G(w)|, (4.3)
F (•, z) is differentiable for every z and satisfies

|∇ x F (x, z)| ≤ µ |G(z)| 2q-2 q (4.4) |∇ x F (x, z) -∇ x F (x, w)| ≤ µ |G(z)| q-2 q + |G(q)| q-2 q |G(z) -G(w)| (4.5) |∇ x F (x, z) -∇ x F (y, z)| ≤ µ |x -y| |G(z)| 2q-2 q (4.6)
Let finally f ∈ W 1,p loc (Ω) with p = q and let us consider the equation

-div F (x, ∇u) = f. (4.7)
Using Nirenberg's method of incremental ratios, we then get the following result.

Theorem 4.1. Let q ≥ 2, p = q , f ∈ W 1,p loc (Ω), F and G be vector fields that satisfy conditions (4.1)-(4.6) above and let u ∈ W 1,q loc (Ω) be a local weak solution of (4.7). Then G := G(∇u) ∈ W 1,2 loc (Ω).

Proof. To shorten notations, set

F := F (• , ∇u(•)) (note that F ∈ L p loc (Ω) and G ∈ L 2 loc (Ω) thanks to (4.1)
) and denote by τ h ϕ := ϕ(• + h) the translate of the function ϕ by the vector h. Let ϕ ∈ W 1,p (Ω) be compactly supported in Ω and

h ∈ R d \ {0} be such that |h| < dist(supp(ϕ), R d \ Ω), we then have Ω τ h F -F |h| • ∇ϕ dx = Ω τ h f -f |h| ϕ dx. (4.8)
Let

ω ω 0 Ω, let ξ ∈ C ∞ c (Ω) be such that supp(ξ) ⊂ ω 0 , 0 ≤ ξ ≤ 1 and ξ = 1 on ω and h ∈ R d \ {0} such that |h| ≤ r 0 < 1 2 dist(ω 0 , R d \ Ω).
In what follows, c will denote a nonnegative constant that does not depend on h, but may vary from one line to another. We then insert the test function

ϕ = ξ 2 |h| -1 (τ h u -u),
into (4.8). Using the fact that u ∈ W 1,q loc (Ω), f ∈ W 1,p loc (Ω) and Hölder's inequality and defining ω := ω 0 + B(0, r 0 ), we get

|h| -2 Ω (τ h F -F) • ξ 2 (τ h ∇u -∇u) + 2 ξ ∇ξ (τ h u -u) ≤ ∇f L p (ω ) ∇u L q (ω ) .
Let us now write the left-hand side of the previous inequality as the sum of four terms I 1 + I 2 + I 3 + I 4 , where

I 1 := |h| -2 Ω ξ 2 (F (x + h, ∇u(x + h)) -F (x + h, ∇u(x)) • (τ h ∇u -∇u), I 2 := |h| -2 Ω (F (x + h, ∇u(x + h)) -F (x + h, ∇u(x)) • ∇ξ ξ (τ h u -u), I 3 := |h| -2 Ω (F (x + h, ∇u(x)) -F (x, ∇u(x)) • ∇ξ ξ (τ h u -u),
and

I 4 := |h| -2 Ω (F (x + h, ∇u(x)) -F (x, ∇u(x)) • (τ h ∇u -∇u) ξ 2 .
Thanks to (4.2), the first term satisfies satisfies:

I 1 ≥ ν ξ |h| -1 (τ h G -G) 2 L 2 . (4.9)
As for the second term, if q > 2 using (4.3) and Hölder's inequality with exponents 2, q and 2q q-2 yields

|I 2 | ≤ |h| -2 Ω |ξ∇ξ ||τ h u -u| |τ h G -G| |τ h G| q-2 q + |G| q-2 q ≤ c ∇ξ L ∞ |h| -1 (τ h u -u) L q (ω 0 ) ξ|h| -1 (τ h G -G) L 2 ω 0 |G| 2 + |τ h G| 2 q-2 2q ≤ c ξ |h| -1 (τ h G -G) L 2 ,
and if q = 2, simply using Cauchy-Schwarz inequality, we similarly obtain:

|I 2 | ≤ c ξ |h| -1 (τ h G -G) L 2 .
We now come to the term I 3 . Thanks to (4.4), we have

|h| -1 |F (x + h, ∇u(x)) -F (x, ∇u(x))| ≤ µ |G| 2q-2 q ∈ L p loc ,
using again Hölder's inequality and the fact that |h| -1 (τ h u -u) L q (ω 0 ) is bounded independently of h, we get |I 3 | ≤ C. To handle the last term we first write

F (x + h, ∇u(x)) -F (x, ∇u(x)) |h| = 1 0 ∇ x F (x + sh, ∇u(x)) h |h| ds.
By defining

Φ h (x) := 1 0 ∇ x F (x + sh, ∇u(x)) -∇ x F (x + (1 + s)h, ∇u(x + h)) |h| h |h| ds,
and rearranging terms, we can rewrite

I 4 = Ω ξ 2 τ h ∇u • Φ h dx - Ω τ h ξ 2 -ξ 2 |h| τ h ∇u • 1 0 ∇ x F (x + (1 + s)h, ∇u(x + h)) h |h| ds dx.
Thanks to (4.4), we get a bound on the absolute value of the second term exactly as we did for the term I 3 . To treat the remaining term, using (4.5) and (4.6), we first have

|Φ h | ≤ 1 0 |∇ x F (x + sh, ∇u) -∇ x F (x + (1 + s)h, ∇u)| |h| ds + 1 0 |∇ x F (x + (1 + s)h, ∇u) -∇ x F (x + (1 + s)h, τ h ∇u)| |h| ds ≤ µ |G| 2q-2 q + |h| -1 |τ h G -G| |G| q-2 q + |τ h G| q-2 q
.

Proceeding as for the term I 2 , we obtain

|I 4 | ≤ c (1 + ξ|h| -1 (τ h G -G) L 2 ).
Putting all the estimates together, we then get

ξ τ h G -G h 2 L 2 ≤ c 1 + ξ τ h G -G h L 2
, from which we can finally infer

τ h G -G h 2 L 2 (ω) ≤ c,
for some constant c depending on q, f W 1,p , u W 1,q and the distance between ω and ∂Ω, but not on h. This finally proves the desired result, namely that G ∈ W 1,2 loc (Ω).

Applications

We shall now apply Theorem 4.1 to our model equation

- d i=1 ∂ i (|∂ i u| -δ i ) q-1 + ∂ i u |∂ i u| = f, (4.10)
which corresponds to

F (z) = (F 1 (z 1 ), . . . , F d (z d )), with F i (z i ) := (|z i | -δ i ) q-1 + z i |z i | , in (4.7
). We then define

G(z) = (G 1 (z 1 ), . . . , G d (z d )), with G i (z i ) := (|z i | -δ i ) q 2 + z i |z i | ,
and again we assume that q ≥ 2. It is immediate to check that condition (4.1) holds and since F does not depend on x, (4.4), (4.5), (4.6) are straightforward. As for (4.2) and (4.3), we have the following.

Lemma 4.2. Let F and G be defined as above with q ≥ 2, then for every (z, w) ∈ R d × R d , the following inequalities hold,

|F (z) -F (w)| ≤ (q -1) |G(z)| q-2 q + |G(w)| q-2 q |G(z) -G(w)|, (4.11)
and

(F (z) -F (w)) • (z -w) ≥ 4 q 2 |G(z) -G(w)| 2 . (4.12)
Proof. First recall (see for instance [START_REF] Lindqvist | Notes on the p-Laplace equation[END_REF]), that for any (a, b) ∈ R 2 , one has

|a| q-2 a -|b| q-2 b ≤ (q -1) |a| q-2 2 + |b| q-2 2 |a| q-2 2 a -|b| q-2 2 b . (4.13) Taking a = (|z i |-δ i ) + z i /|z i | and b = (|w i |-δ i ) + w i /|w i |, in the previous inequality, we thus get |F i (z i ) -F i (w i )| ≤ (q -1) |G(z)| q-2 q + |G(w)| q-2 q |G i (z i ) -G i (w i )|,
and then (4.11) directly follows.

We shall now prove that for every i and z i , w i in R, there holds

(F i (z i ) -F i (w i ))(z i -w i ) ≥ 4 q 2 (G i (z i ) -G i (w i )) 2 , (4.14)
from which (4.12) will directly follow. First note that (4.

(F i (z i ) -F i (w i ))(z i -w i ) = (|z i | -δ i ) q-1 (|z i | - z i |z i | w i ) ≥ (|z i | -δ i ) q = G i (z i ) 2 .
As for the case |z i | > δ i and |w i | > δ i , again taking a = (|z i | -δ i ) + z i /|z i | and b = (|w i | -δ i ) + w i /|w i | in the following inequality (for which we again refer to [START_REF] Lindqvist | Notes on the p-Laplace equation[END_REF])

(|a| q-2 a -|b| q-2 b)(a -b) ≥ 4 q 2 |a| q-2 2 a -|b| q-2 2 b 2 , we get 4 q 2 (G i (z i )-G i (w i )) 2 ≤ (F i (z i )-F i (w i ))(z i -w i )-δ i (F i (z i )-F i (w i )) z i |z i | - w i |w i | .
We then observe that we then obtain the following Sobolev regularity for the (unique) minimizer. Proof. By duality, we know that σ is related to any solution of the dual problem u through the componentwise relation is in W 1,2 loc (Ω) and we obtain the desired result exactly as in [START_REF] Brasco | Congested traffic dynamics, weak flows and very degenerate elliptic equations[END_REF], by observing that

(F i (z i )-F i (w i )) z i |z i | - w i |w i | = (|z i | -δ i ) q-1 + (|w i | -δ i ) q-1 1 - z i w i |z i w i | ≥ 0,
σ i = (|∂ i u| -δ i ) q-1 + ∂ i u |∂ i u| , i =
σ i = |G i | q-2 q G i .
Let us observe that when d = 2 (the case which is relevant for applications to network congestion), the previous result implies that σ ∈ L s loc for every s > 1. As shown recently in [START_REF] Brasco | On certain anisotropic elliptic equations arising in congested optimal transport: local gradient bounds[END_REF], this higher integrability result is still valid for d ≥ 2 and can be proved without appealing to the Sobolev result of Corollary 4.3, by proving directly that any solution of the dual problem has a gradient in L s loc , for every s > 1. We point out that the proof in [START_REF] Brasco | On certain anisotropic elliptic equations arising in congested optimal transport: local gradient bounds[END_REF] requires f to be in L ∞ loc , an hypothesis which is not directly comparable with that of Corollary 4.3.

Remark 4.4. The previous Sobolev result can be easily generalized to equations with weights such as

- d i=1 ∂ i b i (x)(|∂ i u| -δ i ) q-1 + ∂ i u |∂ i u| = f,

which finally proves ( 4 . 1 p

 41 [START_REF] Ekeland | Convex analysis and Variational problems[END_REF].If we go back to the variational problem of Beckmann type:inf σ∈L p (Ω) Ω d i=1 |σ i (x)| p + δ i |σ i (x)| dx : -div σ = f (4.15)

Corollary 4 . 3 .if q = 2 ,

 432 If q ≥ 2 and f ∈ W 1,p loc (Ω) then the solution σ of (4.15) belongs to the Sobolev space W 1,r loc (Ω), any value < 2, if q > 2 and d = 2, dq dq-(d+q)+2 , if q > 2 and d > 2.

1 , 2 +

 12 . . . , d.Since u ∈ W 1,p (Ω) is a weak solution of the Euler-Lagrange equation (4.10), it follows from Theorem 4.1 and Lemma 4.2 that the vector field whose components areG i := (|∂ i u| -δ i ) q ∂ i u |∂ i u| , i = 1, . . . , d,

In particular, the fact that the travelling time functions on arcs [x, x + εv k ] scale like ε h k (x, m/ε) and that the discrete measures f Gε + and f Gε -weakly converge to some f + and f -.

Interestingly, the connection with the Monge-Kantorovich theory was realized much later by Robert McCann.

Acknowledgements. We thank Filippo Santambrogio for some useful discussions. This work has been supported by the ANR through the projects ANR-09-JCJC-0096-01 EVAMEF and ANR-07-BLAN-0235 OTARIE, as well as by the ERC Advanced Grant n. 226234.

as well as to non symmetric equations of the form

and finally to equations of the form (3.4) which both have weights and are nonsymmetric.