Lingli Deng
email: denglingli@chinamobile.com

Ziyao Xu
email: ziyao.xu@alibaba-inc.com

Wei Chen
email: chenweiyj@chinamobile.com

Lu Lu
email: lulu@chinamobile.com

Xiaodong Duan
email: duanxiaodong@chinamobile.com

Efficient Integrity Protection for P2P Streaming

Keywords: hash, chunk, integrity, P2P streaming

This paper proposes an efficient checksum consolidation method, where the content source publishes only a few consolidated checksums (as the direct verification proof) for the media content. A tree-based and a chain-based scheme are proposed to consolidate chunk checksums respectively for VoD and live streaming application scenarios.

Introduction

In a typical P2P streaming system [START_REF] Vu | Measurement of a large-scale overlay for multimedia streaming[END_REF], there exist several well-known portals for the content sources to publish their local media content to potential audience accessing the network through the other peers. For a published media program, its responsible tracker server as well as other human readable description information, are returned to a querying peer by the portal, as the response to the latter's program selection request. The responsible tracker holds the IP/port address list for correspondent swarm (i.e. peer group uploading/downloading the same program) members. By submitting a swarm joining request to the responsible tracker, a downloading peer receives the peer list (containing sharing peers' IP/Port addresses, etc.) and hands in its own IP address and Port number to the tracker to update the swarm's peer list. From this point, a downloading peer can initiate media downloading request to other intermediate peers in the list.

To exploit the many-to-many sharing pattern of the P2P protocol, a large media file is divided into a group of chunks and each chunk is distributed independently afterwards [START_REF] Zhang | Problem Statement of P2P Streaming Protocol (PPSP)[END_REF]. In other words, a downloading peer needs only a single chunk, rather than the complete chunk set for the original media file, to become a server for uploading the local chunk to other peers. It is desirable from the respects of both the service provider and the downloading peers, to ensure that the media content is exactly the same as published and not manipulated by any intermediate party, especially when the provider holds certain reputation/authority/responsibility for the information's authenticity/validity it delivers to the public. This work is supported by Chinese National Grant 2012ZX03002008.

However, previous P2P streaming systems put little concern on integrity protection. Moreover, particular requirements for streaming applications to achieve chunk-level integrity protection make previous solutions for P2P file-sharing networks inapplicable to streaming settings. On the contrary, the consolidated method proposed in this paper is highly applicable to P2P streaming systems, for its cost-efficiency in achieving chunk-level integrity protection for chunk sets, achieved by employing hash consolidation methods (i.e. hash tree based scheme for VoD applications and hash chain based scheme for live streaming applications). It is expected to bring down the storage and verification cost for downloading peers, the publishing cost for content sources, and the storage cost for system portal/trackers. The rest of this paper is organized as follows: Section 2 gives a brief review on related work. Section 4 analysis integrity protection requirements for P2P Streaming applications. Section 3 introduces the hash tree and hash chain mechanisms to be utilized by our proposals in Section 5. In section 6, a preliminary analysis is conducted. Section 7 concludes.

Related Work

Previous P2P file-downloading systems take measures to ensure file integrity from the publisher all the way to the ultimate downloading peer. For example, BitTorrent users can publish a file's hash checksum to the system portal/tracker, to be used later by the downloading peer to verify a file's integrity after downloading the complete chunk set and rebuild the original file. However, in BitTorrent's scheme, publisher only provides the single checksum for the whole file, a downloading peer cannot locate the specific chunk(s) if only part(s) of the chunk set get manipulated, which means a single manipulated chunk would cause the searching and downloading the complete file all over again. Therefore, eMule enhances BitTorrent's scheme allowing content sources to incorporate the hash checksums for each individual chunk into the published content handle, so that chunk-level integrity verification and re-transmission are viable for the downloading peer. Another drawback of BitTorrent's scheme is that the total verification work is performed after the whole set of chunks are downloaded completely, indicating an insufferable delay to the completion of file downloading transaction, indicating devastated user experience. To conquer the problem, [START_REF] He | Latency Hiding Algorithm for P2P File Integrity Verification[END_REF] proposed to exploite the iterative feature of the hash algorithm, in order to hide the check delay into the file downloading process, where the downloading peer uses the hash algorithm to compute a chunk's checksum and updates the file's hash vector whenever it received a new chunk. After downloading the last chunk, a peer only has to compute a single chunk's hash checksum before building the file's checksum and is able to complete integrity verification rapidly.

Background

A single hash tree [START_REF] Merkle | Protocols for Public Key Cryptosystems[END_REF] is able to protect the integrity of a large number of unsafely stored data objects, once the tree's root hash is safely stored. A leaf node is created for each data object, and each internal node is the hash value of the concatenation of its direct children. It is computationally infeasible for an attacker to construct another hash tree out of some modified leaf nodes to produce the same root hash, as the internal nodes (including the root hash) are constructed with collision-free hash algorithms [START_REF] Goldreich | Foundations of cryptography: Basic applications[END_REF]. See Algorithm 1, where tree(S, H, d) calculates the root hash r for the d-degreed tree constructed from an array S using a hash functions from H.

Algorithm 1 tree(S, H, d) k ← 1; s ← S.length while s > 1 do O ← ∅; n ← s/d ; temp ← S[k] for i = 1 to n do for j = 1 to d -1 do temp ← temp||S[k + 1]; k ← k + 1 end for O[i] ← H(temp) end for s ← n; S ← O end while return O[1]
To check for the integrity of a leaf node, one needs to repeatedly compute the hash checksum for the parent node by concatenating the values of the node in question and its siblings and move upwards along the tree until the root is arrived. If the computed value is the same as the securely stored root hash, the leaf's integrity is intact. To update the value of a leaf node, one needs first to check the leaf's integrity (as stated above). If the value is intact, he/she can go further to modify the leaf's value, and update the leaf's parent, its grandparent, until the root hash is updated accordingly.

Intuitively, a hash chain is a hash tree, where "the height of each right-side sub-tree of any internal node is 1". Hash chains can be used to protect the integrity of a sequence.

Problem Statement

Definition 1 (P2P Chunk Integrity). A chunk set has integrity, if each of its element maintains hash-integrity during the transmission chain from the content provider all the way to the downloading peer (including all the intermediate peers who downloaded the media chunk and uploaded to the next peer along the way).

To ensure chunk-level integrity, previous P2P file-downloading systems allow a content source to provide the hash checksums for each chunk (in the form of a hash vector) as part of the descriptive information of the intended file. A downloading peer is able to verify a chunk's integrity by locally calculating its hash checksum and check the value against correspondent element in the hash vector published by the portal.

However, as summarized by Table 1, streaming applications differ from the file downloading applications in dictating the following specific requirements for chunk integrity protection. First of all, in streaming applications, where the downloading peer is always holding a small number of chunks (not necessarily the whole chunk set for the media program) and is required to perform chunk-level integrity verification before playing the media program locally. In other words, they dictate real-time verification requirement that the local verification of individual/subset of media chunks cannot rely on the program's global information (e.g. the complete chunk set or the whole vector of hash checksums).

Secondly, the non-stop feature of some 7 × 24 live programs makes their content sources themselves is incapable of computing and publishing the global checksum for the media beforehand statically. Indicating partial publication requirement that the publication of checksum information cannot rely on the program's global information.

Thirdly, it is intuitively that until the portal delete the relevant downloading link(s) from its resource sharing list out of management reasons, the system should support the efficient downloading and correct verification of media chunks of a VoD media file, once there is a live peer group who together own a complete chunk set of the file, even when the content source leaved the system temporarily or permanently. In short, VoD applications also demand off-line verification requirement that the local verification for individual/subset of media chunks of a VoD media program cannot rely on the program's source's presence.

Consolidated Integrity Verification

The heart of our proposal, is the introduction of a peer-to-peer paradigm for the hash checksum generation and distribution, for the purpose of reducing the publication and storage signaling expenditures for potentially centralized portal/tracker in ensuring chunk-level integrity protection.

The proposed method mainly involves three procedures in a P2P streaming system, including3 : (1) a source-portal interaction for media and direct checksum publication, when the content source to compute the root hash and its signature locally and publishes them to the system portal; (2) a peer-peer interaction for indirect checksum; (3) generation and distribution, when the intermediate sharing peer to compute the verification paths for locally verified chunks and share them with other downloading peers; and (4) a local verification procedure, when the downloading peer verifies a newly downloaded chunk's integrity employing the indirect checksums (verification paths) from other swarm members as well as the direct checksums (root hashes) from the system portal.

Chain-based scheme for live streaming

The root hash of a published live channel s(t), is published to the portal and updated at time t on a periodic basis (e.g. for every u chunks), by the content source using the hash chain algorithm, while the previous value s(t -1) is also stored by the portal to be used by the downloading peer to verify chunk(s)'s integrity within the latest updating period (t-1,t]. The intermediate sharing peer, after downloading a media chunk x within the latest updating period, computes its hash checksum h .u] from any intermediate sharing peer, as well as the root hash pair s(t) and s(t -1) from the portal, to verify chunk x's integrity. Root hash's formulation and publication Assume that a content source P eer -S wishes to publish a live channel X to the portal. The correspondent process of the channel's root hash's formulation and publication is as follows.

-Step1: P eer -S locally produces the descriptive information m(X) for the live channel X, and submits it, as well as the first set of u chunks x 1u 's root hash s(0) (computed through the hash chain algorithm), to the portal for publication.

s(0) = chain(X[1..u], H, d) (1)
-Step2: For each updating period (i.e. the time interval for playing u media chunks), P eer -S locally computes the current root hash s(t), and submits it to the portal for update. To respond, the portal stores s(t) and s(t -1) for the channel.

s(t) = chain(S(t -1)||X[1..u], H, d) (2)
Media chunks' integrity verification Assume that P eer-D wishes to downloads the xth media chunk of the latest updating period from an intermediate sharing peer P eer -U and checks the chunk's integrity locally. Two types of local verification processes are provided: an eager downloading peer may choose Process 2, while a prudent peer may choose Process 1.

1. Local verification process 1: P eer-D checks a chunk's verification path's validity before downloading it and perform the correspondent integrity verification immediately afterwards.

-Step1: P eer -D checks if the local verification vector D contains x's verification path s x : if so, executes Step3; or Step2, otherwise. -Step2: P eer -D asks P eer -U for verification path

s x = D [1..u],
and verifies the validity of the latter's response s x through the following equation, according to the root hash pair s(t) and s(t -1) from the portal:

s(t) = chain(S(t -1)||D [1..u], H, d) (3)
If the equation holds, the path is valid, update the local verification vector D = s x and executes Step3; otherwise, give up P eer -U and choose another intermediate sharing peer for chunk x. -Step3: After downloading chunk x from P eer -U , P eer -D computes its hash checksum h x and verifies its integrity through the following equation, according to the locally stored verification path s x :

h x = s x [x] (4)
If the equation holds, chunk x is downloaded successfully; otherwise, give up P eer -U and choose another intermediate sharing peer for downloading x.

Local verification process 2: P eer -D performs the integrity verification

after successfully downloading all the chunks of the latest updating period.

-Step1: After downloading chunk x from P eer -U , P eer -D computes its hash checksum h x and store it to the correspondent element D[x] of the verification vector D. -Step2: If all the chunks within the latest updating period are downloaded locally, execute Step3; otherwise, terminate the current chunk downloading transaction. -Step3: P eer -D performs the integrity verification for the chunk set locally downloaded for the latest updating period, through using the local verification vector D as well as the root hash pair s(t) and s(t -1) from the portal, according to the following equation:

s(t) = chain(S(t -1)||D[1..u], H, d) (5)
If the equation holds, terminate the chunk downloading transaction for the latest updating period. Otherwise, starting from the first media chunk of the current period, request the correspondent intermediate sharing peer for the chunk's verification path, until a valid path, satisfying the following equation, is received or no other peer to request:

s(t) = chain(S(t -1)||D [1..u], H, d) (6)
If the received D is valid, P eer -D locates a set of integrity-corrupted chunks S, by comparing h xi and D [i]. Otherwise, S is set of the complete set of chunks within the latest updating period. -Step4: P eer-D through the "local verification process 1" for relocating, downloading and verifying each chunk in set S.

Verification path's computation and distribution Assume that P eer -U wishes to share its locally downloaded chunk x, and the correspondent verification path for other downloaders. Intuitively, there are two possibilities:(a) P eer -U downloaded the valid verification path s x as well as x itself from some other intermediate sharing peer P eer -X through the "Local verification process 1"; or (b)P eer -U downloaded the complete set of media chunks through "Local verification process 2" and acquire correspondent verification paths by local hashing and comparison to the root hash pair. As a result, P eer -U may directly share its locally stored verification vector D[1.

.u] to another downloading peer as the verification path for any chunk x of the latest updating period.

No extra processing is required.

Tree-based scheme for VoD

The root hash of a published VoD file s, is published to the portal by the content source using the hash tree algorithm based on the complete chunk set of the file. The intermediate sharing peer, after downloading a media chunk x, computes its hash checksum h x and stores the value into the local verification tree D as the xth leaf on the tree. The hash values of the sibling nodes along the tree path from the correspondent leaf node to the root constitute the verification path for chunk x, i.e. s x = tree -path(X, d, n, x). (See Algorithm 2). The downloading peer needs only to acquire the verification path s x from any intermediate sharing peer, as well as the root hash s from the portal, to verify chunk x's integrity through Algorithm 3.

Root hash's formulation and publication Assume that a content source P eer -S wishes to publish a VoD file X to the portal. P eer -S locally produces the descriptive information m(X) for the file X, and submits it, as well as the root hash s(0), to the portal for publication. s(0) is computed through the hash tree algorithm where the complete chunk set of the file constitute the leaves of the tree and its degree d is determined locally based on the file's scale and average storage and computation expenditure for verifying the integrity of each chunk. (See Section 6 for reference.)

s(0) = tree(X[1..n], H, d) (7)
Intuitively, there are two possibilities: (a) P eer -U downloaded the valid verification path s x as well as x itself from another intermediate sharing peer P eer -X, and checked x's integrity; or (b) before downloading x, P eer -U finds that a valid path is available from the current hash tree D stored locally for file X (i.e all the hash values that constitute x's verification path have been computed and stored for some other chunks downloaded earlier). As a result, P eer -U may directly share its locally stored verification path s x to another downloading peer.

Analysis

In Table 2 and3, variable n stands for the number of chunks contained by an average VoD media file, which is proportional to the length of the file. m represents the average number of chunks downloaded by a participating peer, which is correlated to the coding rate of the media file and the average on-line duration of a typical user. System parameters u and d are constants determined by the system administrator, with the former depicts the number of chunks of an updating period, and the latter the degree of a VoD hash tree. In comparison, System A in the tables refers to a VoD application using a previous method with each chunk's checksum included into a vector as part of the publication data stored at portal/tracker, while system B (live streaming scenario) and C (VoD scenario) employ the proposed chain-based and tree-based schemes, respectively. In system A, a downloading peer may verify the integrity of a newly downloaded chunk by using the global hash vector from the portal/tracker, therefore fulfilling the requirements of real-time verification and off-line verification. However, it is inapplicable to the live streaming scenario at all since it cannot meet the partial publication requirement as dictated in live streaming scenarios. The content source needs to compute the hash checksums for each media chunk in the form of a hash vector with length n (where n depicts the average number of chunks for a VoD file), resulting in O(n) computation and communication costs for the provider in the content publication procedure, O(n) storage and communication costs for the portal/tracker in checksum distribution procedure, and O(n) computation and/or communication costs for the downloading peers in integrity verification procedure. On the other hand, the proposed two schemes respectively satisfy the requirements in terms of live and VoD scenarios. In both settings, the portal/tracker only needs to store the root hashes for a published program, fulfilling the partial publication requirement and reducing publication expenditures for both the provider and portal/tracker to O(1), with the exception that for live streaming scenario where the root hash of a program is updated periodically (i.e. once every u chunks) by the provider the publication cost is reduced to O(n/u).

In terms of the worse case computational cost for the downloading peer to perform integrity verification of a given chunk, system B and C seem to degrade the O(1) performance respectively to O(u) and O(log d (n) in order to verify the first chunk of a media content. However, the average verification cost for each chunk will decrease because the internal nodes' values are filled into the chain/tree as the downloaded chunk group grows, and finally approaches O(1).

Conclusion

This paper proposes an efficient checksum consolidation, distribution and verification scheme to ensure chunk-level integrity protection in P2P streaming systems. The content source publishes only consolidated checksums to the network. A peer securely verifies a downloaded chunk's integrity by combining the checksums from other intermediate peers as well as the one from the content source. Specifically, a tree-based and a chain-based scheme are proposed respectively for VoD and live streaming application scenarios. Preliminary analysis shows that the proposal is highly applicable to the streaming settings and more efficient compared to previous solutions.

  x and stores the value into a u-length verification vector D as D[x]. The verification vector D[1..u] constitutes the verification path s x for chunk x, i.e. s x . = D[1..u]. The downloading peer needs only to acquire the verification path D[1.

Table 1 .

 1 Special Requirements for Media Chunk Integrity Verification

	Application Real-Time Verification Partial Publication Off-line Verification
	Live	•	•
	VoD	•	•

Table 2 .

 2 Requirement-Satisfaction Analysis

	Application Real-Time Verification Partial Publication Off-line Verification A (VoD-vector) √ √ × B(Live-chain) √ √ -C(VoD-tree) √ √ √

Table 3 .

 3 Verification expenditures for each chunk

	Application Publication Verification Storage
	A(VoD-vetor)	O(n)	O(1)	O(n)
	B(Live-chain)	O(n/u)	O(1)	O(1)
	C(VoD-tree)	O(1)	O(1)	O(1)

For systems employing centralized tracker architecture, the tracker may take the responsibility of checksum publication and our proposal also applies by replacing the portal's role with tracker.

Algorithm 2 tree -path(T, x) P ← ∅; temp ← getleaf (T, x) while temp = T.root do AddSibling(T, temp)toP ; temp ← P arent(T, temp) end while return P Algorithm 3 tree -verif y(P, r, x, d)

Media chunks' integrity verification Assume that P eer-D wishes to downloads a media chunk x of a VoD file X from an intermediate sharing peer P eer-U and checks x's integrity locally.

-Step1: P eer -D checks if the local hash tree D contains x's verification path s x : if so, executes Step3; or Step2, otherwise. -Step2: P eer -D asks P eer -U for verification path s x , and verifies the validity of the latter's response s x through the following equation, according to the root hash pair s(0) from the portal:

If the equation holds, the path is valid, update the local hash tree D with the hash values contained by s x and executes Step3; otherwise, give up P eer -U and choose another intermediate sharing peer for chunk x. -Step3: After downloading chunk x from P eer -U , P eer -D computes its hash checksum h x and verifies its integrity through the following equation, according to the locally stored verification path s x :

If the equation holds, chunk x is downloaded successfully; otherwise, give up P eer -U and choose another intermediate sharing peer for downloading x.

Verification path's computation and distribution Assume that P eer -U wishes to share its locally downloaded chunk x, and the correspondent verification path for other downloaders.


Gwenaël Moysan
email: gwenael.moysan@globms.com

Mehdi Senouci
email: mehdi.senouci@centralesupelec.fr

Antoine D'autume

Daniel Cohen

Rodolphe Dos

Santos Ferreira

Patrick Pintus

Gilles Saint-Paul

Robert M Solow

Bertrand Wigniolle

A note on 2-input neoclassical production functions *

Keywords: Production function, elasticity of substitution, capital share, labor share, Solow model

In this short note, we show how the space of elasticity of substitution functions maps into the space of 2-input neoclassical production functions. In doing so we derive a general analytical formula for every 2-input neoclassical production function of class 2 . We present a simple set of sufficient conditions for the Inada conditions to hold; and prove that the Solow model under capital-augmenting (or investment-specific) technical change is asymptotically balanced if and only if the capital share converges to a non-degenerated limit as the capital-labor ratio tends to infinity.

Introduction

Production functions constitute the cornerstone of supply side economics; still, perhaps surprisingly, the economist's standard toolbox contains very few of these. The economist's toolbox contains the often used Cobb-Douglass specification as well as the more recently popular CES functions. Some variable elasticity of substitution (VES) production functions have come up -see for example [START_REF] Revankar | A Class of Variable Elasticity of Substitution Production Functions[END_REF] -but their use has remained marginal in theoretical as well as in empirical literature. Recent evidence on the existence of medium-run growth regimes makes the design of new production functions an issue of particular interest to growth theory and empirics. 1 The prime goal of this paper is to derive an analytical formula for every 2-input neoclassical production function (f ) of class 2 . The method that we present rests on the elasticity of substitution functions (σ). We prove that for any continuous function of * + into itself σ there exists a 2-input neoclassical production function (f ) of class 2 -which is unique up to the choice of two constants -such that σ is the elasticity of substitution function associated to f . In short, we prove that the space of elasticity of substitution functions is exactly the space of continuous functions of * + into itself, and we map the space of elasticity of substitution functions into the space of neoclassical production functions through an integral formula.

The rest of the paper is organized as follows. In section 1 we derive the general formula for 2-input production functions of class 2 . In section 2 we focus on two types of production functions that are important for pure neoclassical growth theory: the production functions that satisfy the Inada conditions, and the production functions for which the Solow model under capitalaugmenting (or investment-specific) technical change is asymptotically balanced. Section 3 concludes.

A general formula for neoclassical production functions

A neoclassical production function with n ≥ 2 inputs is a constant-returns to scale function of class 2 F : * + n → * + , (X 1 , X 2 , ..., X n ) → F (X 1 , X 2 , ..., X n ) satisfying the assumption of strictly positive and strictly decreasing marginal returns to each input: ∀i ∈ {1, ..., n },

∂ F ∂ X i > 0, ∂ 2 F ∂ X 2 i < 0.
In the case of two inputs (n = 2), and by labelling X 1 = K and X 2 = L , the property of constant returns to scale yields the intensive form of

F : ∀K , L > 0, F (K , L ) = L f K L , with f (•) = F (•, 1
). These assumptions of strictly positive and strictly diminishing marginal returns on F are equivalent to the following assumptions on the intensive form f of F : f > 0, f < 0. Further these assumptions on F are sufficient to guarantee the concavity of F . This preliminary result is proven in appendix A.1; for the rest of the paper we focus on intensive forms. We denote by the set of 2-input neoclassical production functions in intensive form:

= f : * + -→ * + , 2 f > 0, f < 0 .
For any function f ∈ , at any capital-labor ratio k , we can define the elasticity of substitution of f at k by:

σ(k ) = - f (k ) f (k ) -k f (k ) k f (k )f (k ) (1)
Let S be the following space of functions: S = σ : * + -→ * + , 0 . We prove below that S is precisely the space of admissible elasticity of substitution functions.

The relationship (1) can be inverted. Let f ∈ and α(•) be the corresponding capital share function and Π(•) the relative factor shares function:

∀k > 0, α(k ) = k f f ∈ (0, 1) and Π(k ) = α(k ) 1-α(k ) = k f f -k f > 0.
The derivative function of Π(•) is linked to the elasticity of substitution through the following equation:

Lemma 1. Π (k ) Π(k ) = σ(k ) -1 σ(k ) 1 k (2) Proof. d ln Π(k ) d k = d ln k f f -k f d k = 1 k + f f - -k f f -k f = 1 k 1 - f -k f f 1 σ - k f f 1 σ = σ -1 σ 1 k .
As is well known, when σ(k ) < 1 (resp. σ(k ) > 1), α and Π decrease (resp. increase) in response to capital deepening around k .

Let's take k > 0 to be some reference capital-labor ratio -for example, k = 1 -and let's denote by Π = Π( k ) the corresponding relative factor shares ratio. Then, by integrating equation (2) between k and k , it comes that for all k > 0:

Π(k ) = Π exp k k σ(k )-1 σ(k ) dk k
. Hence, f is a solution to the differential equation:

k f f -k f = Π exp k k σ(k )-1 σ(k ) dk k
which can be re-arranged like:

f f = 1 k 1 1 + exp - k k σ(k )-1 σ(k ) dk k Π (3)
Notice that the term 1/ 1 + exp -

k k σ(k )-1 σ(k ) dk k / Π = Π(k ) 1+Π(k ) is equal to the capital share α(k ).
Integrating (3), we conclude that there exists a constant A > 0 such that, for all k > 0:

f (k ) = A exp     k k dk k 1 + exp - k k σ(k )-1 σ(k ) dk k Π     (4)
We can state and prove our first main result, which maps the space of elasticity of substitution functions into the space of neoclassical production functions :

Theorem 1. Let k > 0 be some reference capital-labor ratio, and let Π > 0 and A > 0 be two constants.

Then, for any function σ ∈ S, there exists one and only one production function f ∈ such that (i) σ is the elasticity substitution function associated to f , and (ii) Π( k ) = Π and A = f k . The formula for this function is given by (4).

Proof. Let σ ∈ S and let f be the function of * + into itself defined by equation [START_REF] Greenwood | Long-Run Implications of Investment-Specific Technological Change[END_REF]. By construction, if f ∈ then σ is the elasticity of substitution function corresponding to function f ; so there remains only to demonstrate that f defined in (4) belongs to the set .

Since σ is continuous, f is 2 . f is also nonnegative from (4) and for all k > 0:

f (k ) = f (k ) k 1+ exp ( -k k σ-1 σ dk k ) Π > 0 f (k ) = -1 σ(k ) f (k ) exp - k k σ-1 σ dk k Π k 1+ exp ( -k k σ-1 σ dk k ) Π < 0 which proves that f ∈ .
Thus, all well-behaved neoclassical production functions have a representation in terms of their elasticity of substitution function, which can be any continuous function of * + into itself. 2 For example, if σ ≡ 1, equation (4) becomes:

f (k ) = A exp k k α k dk = A k k α , where α = Π 1+ Π ∈ (0, 1).
If σ ≡ σ, with σ = 1, then:

k k σ-1 σ dk k = σ-1 σ ln k k .
In this case, formula (4) translates into:

f (k ) = A exp   k k dk k 1+ (k / k ) -σ-1 σ Π   = A exp σ σ-1 ln k k k 1 σ + Πk -1 σ-1 ln k k k = k = A Π 1+ Π k k σ-1 σ + 1 1+ Π σ σ-1
Let us define ᾱ = Π 1+ Π . The above equation yields:

f (k ) = A ᾱ k k σ-1 σ + 1 -ᾱ σ σ-1
, which is the expression of the only CES production function of elasticity parameter σ, and such that product at k is equal to A and capital share at k is ᾱ. 3

Two classes of production functions

We study how certain restrictions on function σ ensure that the resulting production functions have some properties that are often required in economic growth models.

The elasticity of substitution function and the Inada conditions

The integral formula (4) permits to link the Inada conditions to the behavior of the capital share near zero and infinity.

Let σ ∈ S, and k , A and Π be strictly positive constants and let f ∈ be the production function such that σ is the elasticity of substitution function associated to f and such that f ( k

) = A and k f ( k ) f ( k )-k f ( k ) = Π. Theorem 1 proves that function f is unique.
f ∈ , so f is strictly increasing and strictly positive over * + while f is strictly decreasing and strictly positive over * + . Hence, f and f admit some limits at the borders of * + . Let's denote these limits by l 0 , l ∞ , l 0 , l ∞ :

l 0 = lim 0 f ∈ + l ∞ = lim +∞ f ∈ * + ∪ {+∞} l 0 = lim 0 f ∈ * + ∪ {+∞} l ∞ = lim +∞ f ∈ +
The Inada conditions hold when l 0 = l ∞ = 0 and l ∞ = l 0 = +∞. 2 Remark that equation ( 4) is equivalent to:

f (k ) = A exp k k α(k )
k dk , where α(•) is the capital share function. This seems to be a simpler representation than (4); however, some restrictions have to be imposed on α(•) to ensure that function A exp

k k α(k )
k dk is a neoclassical production function.

3 See Arrow, Chenery, Minhas and Solow (1961), who first derived the analytical expression of two-inputs CES production function through to the same method we use here.

Call

Π(k ) = k f f -k f the relative factor shares function and α(k ) = k f
f the capital share function. We have seen in section 1 that for all k > 0, Π(k

) = Π exp k k σ(k )-1 σ(k ) dk k , and that α(k ) = Π(k ) 1+Π(k ) . Lemma 2. l 0 > 0 =⇒ lim 0 α = 0 l ∞ < +∞ =⇒ lim +∞ α = 0 l 0 < +∞ and l 0 = 0 =⇒ lim 0 α = 1 l ∞ > 0 and l ∞ = +∞ =⇒ lim +∞ α = 1 Proof. See appendix A.2. 4
This result can be stated in plain English the following way:

Let f be a neoclassical production function. If factor shares do not tend to degenerated values in the neighborhood of infinity, then f fulfills the Inada conditions near infinity. If factor shares do not tend to degenerated values in the neighborhood of zero, then f fulfills the Inada conditions near zero.

Remarkably, the Inada conditions are compatible with ever-fluctuating factor shares. Indeed, asymptotic fluctuations in factor shares constitute a sufficient condition for f to meet the Inada conditions -which also follows if factor shares converge to non-degenerated limits as the capitallabor ratio tends to infinity. In the next section we lay some examples of production functions pertaining to these different cases.

Let's call the set of neoclassical production functions that satisfy the Inada conditions:5

= f ∈ lim 0 f = lim ∞ f = 0, lim ∞ f = lim 0 f = ∞
We also define a special subset of S:

S * = σ ∈ S k 0 σ -1 σ dk k and ∞ k σ -1 σ dk k do not diverge to + ∞ nor to -∞
The results of lemma 2 show that:

Theorem 2. Let f ∈ , and let σ be the corresponding elasticity of substitution function. Then:

σ ∈ S * =⇒ f ∈
The reciprocal of theorem 2 does not hold: there exists some production functions that fulfill the Inada conditions, with factor shares converging to degenerated values. 6 For example, if

σ 1 (k ) = ln k ln k -1 (for k > exp(1)), then σ 1 -1 σ 1 = 1 ln k . A primitive of 1/k ln k is ln(ln k )
, which tends to infinity as k tends to infinity, and so lim k →+∞ α(k ) = 1. In virtue of theorem 1, the corresponding functions f 1 all take the form:

f 1 (k ) = A exp   k k dk k 1+ exp - k k dk k ln k Π   , with k > 0, A > 0 and Π > 0. It is possible to
analytically derive all functions f 1 :

f 1 (k ) = A exp k k dk k 1+ exp(-ln(ln k )+ln(ln k )) Π = A exp k k dk k ( 1+ λ ln k ) = A exp (ln k -λ ln(λ + ln k ) -c ) with λ = ln k / Π > 0 and c = ln k -λ ln(λ + ln k ). Call A = A exp(-c ) > 0, then f 1 admits the fol- lowing analytical expression: f 1 (k ) = A k (λ+ln k ) λ . So f 1 (k ) ∞ -→ ∞
as a consequence of power functions beating powers of logarithms near infinity. By differentiating the above expression, it comes: [START_REF] Robinson | The Classification of Inventions[END_REF] and [START_REF] Uzawa | Neutral Inventions and the Stability of Growth Equilibrium[END_REF] have shown that, in the canonical neoclassical growth framework, if the growth rate of real variables are constant then technical change has a labor-augmenting (LATC) representation, i.e. there exists a production function F (•, •) such that output can be ex-

f 1 (k ) = A 1 (λ+ln k ) λ - λ (λ+ln k ) λ+1 = A ln k (λ+ln k ) λ+1 , which proves that f 1 ∞ -→ 0.

The asymptotically-quasi-Cobb-Douglas production functions

pressed like Y t = F (K t , B t L t )
, where B t = B 0 exp(g B t ) is a steadily-increasing productivity term (g B > 0). Then, output and capital both grow at the constant rate of g B along the stable equilibrium path. This result is known as the 'steady growth theorem'. 7 The Robinson-Uzawa theorem is probably the reason for the widespread representation of technical change as labor-augmenting. 8 However, a whole strand of empirical literature, initiated by Greenwood et al. (1997), has firmly established that technical change was not correctly represented as labor-augmenting in the USA over the postwar period, and even less so since the early 1980'sinstead, technical change seems to be mostly embedded in more efficient (and/or cheaper) capital goods. In addition, it has been noted that the labor share has decreased steadily in most advanced countries between the late 1970's and the early 2010's, which also contradicts the LATC assumption. 9 These empirical biases question the relevance of LATC-driven growth and raises the issue of whether the neoclassical growth framework can account for medium-run imbalances as well as for long-run regularities.

To do so, we investigate the conditions under which a simple neoclassical growth model yields asymptotically (rather than absolutely) constant growth rates. More precisely, we consider be- 7 Jones and Scrimgeour (2008) interpret this result the following way: if there are several inputs to production, some accumulable and others not (e.g., capital vs. labor), then economic growth makes the accumulable inputs grow more quickly than the non-accumulable ones, except when technical change precisely increases the effective quantities of non-accumulable inputs. In the capital and labor inputs framework, technical change then has to be labor augmenting. 8 Barro and Sala-i-Martín (2004, p. 53) reflect the common perception of the range of the Robinson-Uzawa theorem: "If we want to consider models that possess a steady state, we have to assume that technological progress takes the labor-augmenting form. (...) the long-term experiences of the United States and some other developed countries indicate that per capita growth rates can be positive and trendless over long periods of time (...). This empirical phenomenon suggests that a useful theory would predict that per capita growth rates approach constants in the long run; that is, the model would possess a steady state." We prove below that some neoclassical growth models can yield asymptotically constant ratios and growth rates without assuming that technological change is labor-augmenting nor that the production function is Cobb-Douglas. 9 See de La Grandville (2009, chapter 5 co-written with Robert M. Solow) and Neiman and Karabarbounis (2014).

low capital-augmenting (KATC) and investment-specific technical change (ISTC) -which both increase the effective quantity of capital -in a standard, continuous-time Solow model. We prove that the set of production functions for which growth is asymptotically balanced is exactly the set of production functions whose labor share function tends to some limit as the capital-labor ratio tends to infinity; a class of functions that we might dub the asymptotically-quasi-Cobb-Douglas production functions.

We denote by the set of production functions such that the capital share tends to some constant in [0, 1) as the effective capital-labor ratio tends to infinity:

= f ∈ lim x →+∞ x f (x ) f (x ) exists and ∈ [0, 1)
For f ∈ , f ∈ if and only if the corresponding elasticity of substitution function σ f belongs to the set Theorem 3. Let s ∈ (0, 1), δ > 0 and n ≥ 0 and let f ∈ .

S A = σ ∈ S ∞ k σ(k )-1 σ(k ) dk k ∈ [-∞, +∞) . Definition 1. Let f t t ≥0 ∈ [0,
• Let (k t , y t ) t ≥0 be the growth path corresponding to the Solow model with steady KATC and LATC:

y t = B t f A t k t B t kt = s y t -(δ + n )k t ( 5
)
with Ȧt /A t = g A ≥ 0, Ḃt /B t = g B ≥ 0, k 0 > 0, A 0 > 0 and B 0 > 0 given. Then:

(k t , y t ) is asymptotically balanced ⇐⇒      g A = 0 or g A > 0 and f ∈
The asymptotic growth rates of k and y are then respectively: g k = g y = g B in the case g A = 0;

and g k = g y = g B + α ∞ 1-α ∞ g A in the case g A > 0 and f ∈ , where α ∞ = lim x →+∞ x f (x ) f (x ) .
• Let (k t , y t ) be the growth path corresponding to the Solow model with steady ISTC and LATC:

y t = B t f k t B t kt = s q t y t -(δ + n )k t (6)
with Ḃt /B t = g B ≥ 0, qt /q t = g q ≥ 0, k 0 > 0, A 0 > 0 and q 0 > 0 given. Then:

(k t , y t ) is asymptotically balanced ⇐⇒      g q = 0
or g q > 0 and f ∈ The asymptotic growth rates of k and y are then respectively: g k = g y = g B in the case g q = 0;

and g k = g B + 1 1-α ∞ g q , g y = g B + α ∞ 1-α ∞ g q in the case g q > 0 and f ∈ , where α ∞ = lim x →+∞ x f (x ) f (x ) .
Proof. See appendix A.3.

Theorem 3 highlights the set of production functions, as well as the corresponding elasticity of substitution functions, such that the Solow growth model under KATC or ISTC yields asymptotically balanced growth. It comes out that these functions are exactly the functions in . As such, theorem 3 constitutes an extension of the Robinson-Uzawa theorem, which implies that balanced (i.e. constant-rate) growth under KATC or ISTC is only possible if the production function is Cobb-Douglas. 10 Notice that there are some production functions in such that the associated elasticity of substitution function does not tend to one near infinity. More precisely, it is possible that f ∈ with σ oscillating around 1. Take, for instance,

σ 2 (k ) = 1 1-sin k 2 for all k ≥ 0, so that σ 2 (k )-1 σ 2 (k ) = sin k 2 .
Since function x → sin x /x is integrable near infinity, for any function f 2 corresponding to elasticity of substitution σ 2 , the relative factor shares ratio Π(k ) tends to a strictly positive constant as k tends to infinity. Still, σ 2 does not tend to one, as σ 2 oscillates periodically between 2/3 and 2.

The effective inputs ratio tends to a constant under pure LATC while it diverges under KATC/ISTC. Consequently, the production function influences growth dynamics in the KATC/ISTC-driven models, in which case theorem 3 proves that fluctuations of the labor share around a long-run mean is compatible with asymptotically balanced growth.

Conclusion

In this paper we have investigated the links between the set of 2-inputs neoclassical production function and the set of elasticity of substitution functions. In doing so we have unveiled an integral formula for all these production functions of class 2 . We linked this formula to the Inada conditions and to the asymptotic balance of a standard growth model.

The generalization of this formula to the case of n ≥ 3 inputs -which today is beyond our means -would be an answer to the question: 'how to mathematically represent constant-returnsto-scale aggregation processes?' and would as such represent a powerful tool for neoclassical analysis. terms: f : * + → * + defined by ∀k ∈ * + , f (k ) = F (k , 1). Then:

A Appendix

∀(K , L ) ∈ * + 2 ,          ∂ F ∂ K (K , L ) > 0 ∂ F ∂ L (K , L ) > 0 ∂ 2 F ∂ K 2 (K , L ) < 0 ∂ 2 F ∂ L 2 (K , L ) < 0 ⇐⇒ ∀k ∈ * + , f (k ) > 0, f (k ) < 0.
Proof. =⇒ Assume that the left-hand side conditions are fulfilled.

By the definition of f ,

f (•) = F (•, 1) so ∀k ∈ * + , f (k ) = ∂ F ∂ K (k , 1) > 0 and ∀k ∈ * + , f (k ) = ∂ 2 F ∂ K 2 (k , 1) < 0.
⇐= Assume that the right-hand side conditions are fulfilled.

Let (K , L ) ∈ * + 2 . By the definition of f , F (K , L ) = L f K L .
Differentiating this identity with respect to K twice yields:

∂ F ∂ K (K , L ) = f K L > 0 ∂ 2 F ∂ K 2 (K , L ) = 1 L f K L < 0 Differentiating identity F (K , L ) = L f K L with respect to L yields: ∂ F ∂ L (K , L ) = f K L -K L f K L .
Since f is assumed to be strictly positive and strictly increasing, f admits a limit near 0 + which we call f (0) ∈ + . Then:

f (K /L ) = f (0)+ K /L 0 f (u
) du ; so the latter identity can be re-written:

∂ F ∂ L (K , L ) = f (0) + K L 0 f (u ) du - K L f K L = f (0) + K L 0 f (u) -f K L du Since f is strictly decreasing, then f (u ) -f K L > 0 for all u < K /L , so K L 0 f (u ) -f K L du > 0. Since f (0) ≥ 0, we conclude that ∂ F ∂ L (K , L ) > 0. Differentiating the identity F (K , L ) = L f K L twice with respect to L yields: ∂ 2 F ∂ L 2 (K , L ) = - K L 2 f K L + K L 2 f K L - K L - K L 2 f K L = K 2 L 3 f K L < 0.
It should be noted that the above computations show that for any neoclassical production functions in extensive terms of class 2 F , the Hessian matrix of F have a zero determinant and a strictly negative trace. The non-diagonal elements are:

∂ 2 F ∂ K ∂ L = ∂ ∂ L f K L = - K L 2 f K L This leads to det H = ∂ 2 F ∂ K 2 ∂ 2 F ∂ L 2 -∂ 2 F ∂ K ∂ L 2 = 0 and Tr H = ∂ 2 F ∂ 2 K + ∂ 2 F ∂ 2 L < 0.
The Hessian matrix is negative semi-definite and so function F is globally concave. 11Constant-returns-to-scale functions in extensive terms are never strictly concave. Let n ∈ * and F be a CRS function of n into R . Let x ∈ n \{(0, • • • , 0)} and let µ ∈ \{1} and y = µx . Then, for any λ ∈ (0, 1):

F λx + (1 -λ)y = F (λ + µ(1 -λ))x = (λ + µ(1 -λ))F (x ) = λF (x ) + (1 -λ)F (y )
which proves that F is not strictly concave.

A.2 Proof of lemma 2

• If l 0 > 0, then by the definition of f in (4), function k →

1 k 1+ exp - k k σ-1 σ dk k Π
is a (strictly) positive function that is integrable in the neighborhood of 0. Since function 1 k is not integrable around 0 + , in particular it must hold that:

1 1+ exp - k k σ-1 σ dk k Π k →0 --→ 0. Hence: lim k →0 Π(k ) = 0 and lim k →0 α(k ) = 0. • If l ∞ < +∞, then function: k → 1 k 1+ exp - k k σ-1 σ dk k Π
is positive and integrable around +∞ and so, in particular, must be negligible compared to function 1 k , which is not integrable near 0 + . So

1 1+ exp - k k σ-1 σ dk k Π k →+∞
----→ 0, and then lim k →+∞ α(k ) = 0.

• If l 0 < +∞ and l 0 = 0, then by applying de L'Hôpital's rule we conclude that

f (k ) k k →0 --→ l 0 . But from equation (3), this implies that 1 1+ exp - k k σ-1 σ dk k Π k →0 --→ 1. Hence, lim k →0 k k σ-1 σ dk k = +∞
and so lim k →0 α(k ) = 1.

• If l ∞ > 0 and l ∞ = +∞, then by de l'Hôpital's rule

f (k ) k k ∞
-→ l ∞ and so from equation ( 3):

1 1+ exp - k k σ-1 σ dk k Π k →+∞ ----→ 1. So: lim k →+∞ k k σ-1 σ dk k = +∞ and lim k →+∞ α(k ) = 1.

A.3 Proof of theorem 3

• Let's denote by x t the effective inputs ratio: x t = A t k t /B t . We rewrite equation ( 5) in terms of (x t , y t ):

y t = B t f (x t ) ẋt = s A t f (x t ) -(δ + n -g A + g B )x t (7)
Differentiating the first equation yields:

ẏt

y t = g B + ẋt f (x t ) f (x t ) = g B + α t ẋt x t ( 8
)
where

α t = α(x t ) = x t f (x t )
f (x t ) denotes the relative capital share at date t . Then, by differentiating the second equation in [START_REF] Karabarnounis | The Global Decline of the Labor Share[END_REF] and re-injecting the same equation and the definition of the capital share, it comes that ẋ /x follows the process:

• ẋt x t = g A -(1 -α(x t )) ẋt x t ẋt x t + δ + n -g A + g B ( 9
)
Note that the second term is positive, since kt /k t ≥ -δn at all dates. The differential equation of ẋt /x t (9) is not autonomous in ẋ /x due to the term α(x t ).

We begin by proving a useful lemma:

Lemma 3. If g A > 0, then x t t →+∞ ---→ +∞.
Proof. We first show that x t is ultimately increasing, that is: ∃τ ≥ 0 |∀t ≥ τ, ẋt ≥ 0 -If δ + ng A + g B < 0, then by the second equation in [START_REF] Karabarnounis | The Global Decline of the Labor Share[END_REF], ẋ is positive at each date.

-If δ + ng A + g B > 0, let's define for all t the ratio x * t such that ẋt 0 ⇐⇒ x t x * t * First case: suppose that x 0 < x * 0 , so that ẋ0 > 0. Suppose that there exists some T > 0 such that ẋT < 0, so x * T < x T . Since ẋ is a continuous function of time, we can define τ = sup t 1 < T ẋt 1 ≥ 0 . Then, by continuity of ẋt , ẋτ = 0 [START_REF] Revankar | A Class of Variable Elasticity of Substitution Production Functions[END_REF] ∀t ∈ (τ, T ], ẋt < 0

s A t f (x * t ) = (δ + n -g A + g B )x *
From [START_REF] Revankar | A Class of Variable Elasticity of Substitution Production Functions[END_REF], [START_REF] Robinson | The Classification of Inventions[END_REF] proves that x τ > x T , which brings a contradiction. So when ẋ0 > 0, ẋt ≥ 0 at all dates t ≥ 0. * Second case: suppose now that x 0 ≥ x * 0 , so that ẋ0 ≤ 0. If ẋ stays negative forever, then x t (which is then decreasing and positive) tends to a constant, while ẋt then tends to zero; which is inconsistent with the second equation in [START_REF] Karabarnounis | The Global Decline of the Labor Share[END_REF] taken to the limit when g A > 0. So there exists τ > 0 such that ẋτ > 0. Then, by the same reasoning as in the first case above, ẋt ≥ 0 for all t ≥ τ.

x τ = x * τ . But since τ < T , it comes that x * τ < x * T . So x τ = x * τ < x * T < x T . But since τ < T ,
Then, since x t is ultimately increasing, x t either converges to a constant or diverges to +∞. If x t converges to some constant x l , then it converges monotonically to x l and so ẋt must have limit 0. But these two facts taken together are inconsistent with the second equation in [START_REF] Karabarnounis | The Global Decline of the Labor Share[END_REF] when A t increases steadily, so x t t ∞ -→ +∞.

We now proceed to the proof of the first part of theorem 3 concerning the Solow model under KATC and/or LATC.

=⇒ Suppose that ẏt y t and kt k t tend to constants g y ≥ 0 and g k ≥ 0, then by definition of x , ẋt x t has limit g k + g Ag B as t tends to infinity.

-If g k = g Bg A then ẋt x t t →+∞ ---→ 0. Equation (9) taken to the limit proves that

• ẋt x t t →+∞ ---→ g A (δ + n -g A + g B ) = g A (δ + n + g k ). If g A = 0 then the limit of • ẋt x t
is stricly positive, which is inconsistent with the convergence of ẋt x t to 0. Thus, in this case,

g A = 0. -If g k = g B -g A then ẋt x t t →+∞ ---→ g k + g A -g B .
Taking equation ( 8) to the limit proves that

α t t →+∞ ---→ g y -g B g k +g A -g B .
Consequently, if g A > 0 then from lemma 2 x t t →+∞ ---→ +∞, so the function x → α(x ) admits a limit when x → +∞. Call this limit α ∞ . Equation (9) taken to the limit proves that • ẋt x t t →+∞ ---→ g A -(1 -α ∞ )(g k + g Ag B ) g k + δ + n . This quantity must be equal to zero, otherwise ẋt x t does not converge. The only possibility is that we have g A = (1 -α ∞ )(g k + g Ag B ), which means that α ∞ = 1 and g k = g B + α ∞ 1-α ∞ g A . We get g y from equation ( 8): g y = g B + α ∞ (g k + g Ag B ) = g k . Finally, the limit value of α is g y -g B g k +g A -g B = α ∞ ∈ [0, 1), which proves that f ∈ . ⇐= -If g A = 0, then A t = A 0 for all t ≥ 0, so the second equation in [START_REF] Karabarnounis | The Global Decline of the Labor Share[END_REF] is autonomous in x t . Since f satisfies the Inada conditions, x t t →+∞ ---→ x * characterized by s A f (x * ) = (δ + n + g B )x * and ẋt t →+∞ ---→ 0. From the definition of x t it comes that k t /B t tends to a constant, and does monotonically, which proves that kt /k t t →+∞ ---→ g B . Equation [START_REF] Meade | A Neo-Classical Theory of Economic Growth[END_REF] proves that ẏt /y t t →+∞ ---→ g B .

-If g A > 0 and f ∈ , let's call α ∞ = lim x ∞ x f (x )/ f (x ) ∈ [0, 1). By lemma 3, α t t →+∞ ---→ α ∞ . The differential equation ( 9) then converges to the following one:

• ẋt x t = g A -(1 -α ∞ ) ẋt x t ẋt x t + δ + n -g A + g B = Γ ẋt x t ( 12
)
which intercepts the horizontal axis at g A /(1 -α ∞ ). Let be a small, strictly positive number. Since α t converges to α ∞ ∈ [0, 1), then after some date T 1 α t ∈ (α ∞ -, α ∞ +

).

Let's respectively denote by Γ -and Γ + the following functions:

Γ -(χ) = g A -(1 -α ∞ + )χ χ + δ + n -g A + g B Γ + (χ) = g A -(1 -α ∞ -)χ χ + δ + n -g A + g B
With these definitions, ∀t ≥ T 1 : Γ - ẋt x t ≤ Γ ẋt x t ≤ Γ + ẋt x t . So after T 1 , ẋt /x t increases more than if it was guided by Γ -and less than if it was guided by Γ + . But the processes χ = Γ -(χ) and χ = Γ + (χ) converge respectively to N = g A /(1 -α ∞ + ) and O = g A /(1α ∞ -). So there exists a date T 2 ≥ T 1 such that:

∀t ≥ T 2 , M = N -≤ ẋt x t ≤ O + = P which proves that ẋt /x t t →+∞ ---→ g A /(1 -α ∞ )
. By the definition of x , g k = g B + α ∞ 1-α ∞ g A and by equation [START_REF] Meade | A Neo-Classical Theory of Economic Growth[END_REF]

g y = g B + α ∞ 1-α ∞ g A .
• We denote by u t the effective inputs ratio: u t = k t /B t . Differentiating the first equation in [START_REF] Jones | A New Proof of Uzawa's Steady-State Growth Theorem[END_REF] yields: ẏt

y t = g B + ut u t α(u t ) (13)
where α(u t ) = u t f (u t ) f (u t ) . The evolution of u t is guided by: ut = s q t f (u t ) -(δ + n + g B )u t (14)

A.1 2 -

 2 input production functions in extensive and intensive forms Proposition. Let F be a 2-input, constant returns to scale production function of class 2 in extensive terms: F : * + 2 → * + and f be the corresponding production function of class 2 in intensive 10 Another extension found in the literature is due to Meade (1961, appendix II) who proved that a 2-sector growth model with Cobb-Douglas production functions for consumption and investment, when animated by steady consumption-and investment-specific TFP growth, is (strongly) balanced.

  * t is the limit that x t would take if A stayed constant at A t from date t on. If g A > 0 then x * t increases with time and lim t ∞ x * t = ∞. By the second equation in (7):

t x

We are grateful to an anonymous referee for pointing out a mistake in the proof of this lemma in an earlier version of this paper

Clearly, ⊂ but = , since the CES, non-Cobb-Douglas production functions belong to but not to .

In the same vein, Palivos and Karagiannis (2010) proved that σ → 1 is not a sufficient condition for Inada conditions to hold.

We acknowledge a referee for pointing out this result.

Differentiating equation ( 14) and re-injecting the same equation and the definition of the capital share yields:

The same way as in the KATC case, we can prove the following lemma:

Proof. Same method as in the proof of lemma (2) applied to (14).

We know prove the assertions of the second part of theorem 3 concerning the Solow model under ISTC and/or LATC with the same method as in the KATC/LATC case.

=⇒ Suppose that ẏ /y → g y ≥ 0 and k /k → g k ≥ 0, then by equation (13

-If g k = g B then u /u → 0 and so by equation ( 15)

The only possibility is that this quantity is zero, so g q = 0.

-If

Consequently, if g q > 0 then from lemma 4 u t t →+∞ ---→ +∞, so function u → α(u ) admits a limit as u tends to infinity. Call this limit α ∞ . Equation (15) taken to the limit proves that

But since ut u t converges, the limit of • ut u t must be zero, which proves that g q -(1α ∞ )(g k -g B ) = 0, so α ∞ < 1, which achieves to prove that f ∈ . Also g k = g B + 1 1-α ∞ g q and equation (13) taken to the limit proves that g y = g B + α ∞ 1-α ∞ g q . ⇐= -If g q = 0, u t follows the autonomous, stable process ut = s q 0 f (u t ) -(δ + n + g B )u t , so u t converges to some limit u * > 0 and the growth path is asymptotically balanced, with g k = g y = g B .

-If g q > 0 and f ∈ , the the process for ut /u t converges to the following one:

With the same method than in the Solow model under KATC above, this is straightforward to prove that ut /u t → g q 1-α ∞ . So kt /k t t ∞ -→ g B + g q 1-α ∞ and from equation (13), ẏt /y t t ∞ -→ g B + g q α ∞ 1-α ∞ .


Alfons Salden

Satellite Image Analysis and Processing Tools to be used in Air-Pollution Forecast and Simulation Systems

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


François Févotte
email: francois.fevotte@edf.fr

VERROU: Assessing Floating-Point Accuracy Without Recompiling

Keywords: Monte Carlo Arithmetic, floating-point arithmetic, numerical verification, round-off errors

As an industrial facility relying on numerical simulation to improve the safety and efficiency of its electricity production units, EDF is committed to ensure that all the numerical simulation codes it develops and uses are correctly validated and verified. Within this context, the accuracy of floating-point operations has progressively become one of the important topics to study, especially since computing codes are exploited on ever more powerful hardware to solve ever larger problems. The Verification and Validation (V&V) process should therefore include the monitoring of inaccuracies introduced by floating-point arithmetic, as well as the verification that they are kept within acceptable limits. In this paper, we present the verrou tool, a valgrind-based system which uses Monte Carlo Arithmetic (MCA) to monitor the accuracy of floating-point operations without needing to instrument the source code or even recompile it. This tool is therefore well-suited to be part of an industrial V&V process. It has been successfully tested both on small-scale, well understood numerical applications, and on large-scale, more complex industrial computing codes.

Introduction

As numerical simulation tools take an ever more important place within industrial processes, Verification & Validation (V&V) has become a major concern both for code developers and users. Within this context, and seeing the large amount of work devoted to numerical verification in the academic world, it might come as a surprise that floating-point accuracy assessment only recently started making its way to the industry. Indeed, most simulation codes manipulate floating-point values and perform computations on them using the rules promoted by standard IEEE-754 [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF]. As an ever growing number of such computations are routinely performed, rounding errors have become one of the key sources of inaccuracies in the results, along with discretization and iterative convergence issues. which were until then the only focus of verification studies.

Furthermore, numerous tools exist to diagnose floating-point problems, among which cadna [START_REF] Lamotte | CADNA_C: A version of CADNA for use with C or C++ programs[END_REF] is one of the most advanced. It allows instrumenting the source code, replacing standard floating-point arithmetic by Discrete Stochastic Arithmetic (DSA), in order to assess floating-point inaccuracies, detect their origin in the source code, and follow their propagation throughout the computation. cadna has already been successfully used on large industrial simulation codes [START_REF] Montan | Sur la validation numérique des codes de calcul industriels[END_REF], which makes it a natural point of comparison throughout this paper.

Other tools exist, such as FpDebug [START_REF] Benz | A dynamic program analysis to find floating-point accuracy problems[END_REF], which performs the Dynamic Binary Instrumentation (DBI) of an input executable, in order to compare each floating-point value in the program with one obtained with higher precision. FpDebug outputs a graph indicating, for each floating-point operation, whether inaccuracies come from input data or from the operation itself.

A third tool worth mentioning here is Craft HPC [START_REF] Lam | Dynamic floating-point cancellation detection[END_REF], another tool performing DBI in order to detect cancellation errors, and furthermore optimize the use of single-and double-precision variables throughout the source code in order to balance speed and accuracy.

However, despite the availability of these tools and the growing part taken by round-off errors in results inaccuracies, only a few industrial V&V processes include numerical verification of simulation codes. A few reasons might be proposed to explain this observation: first, some development teams are not aware of the ubiquity of rounding errors. Some, on the contrary, are convinced that floating-point inaccuracies are a fatality against which nothing can be done. Second, and more importantly, analyzing an industrial code is a far more difficult task, when compared to smaller programs:

• instrumenting the source code with cadna can be hard, and there generally are numerous dependencies to third-party software libraries of which the development team only has limited understanding.

• Industrial simulation codes routinely produce gigabytes of results, an in-depth analysis of which is in itself often a challenge. In this context, the size of graphs output by FpDebug or of lists of cancellations detected by Craft HPC can quickly become overwhelming. As a consequence, few developers are willing to initiate such work before being certain that the code quality will improve.

In this paper, we introduce verrou,1 a tool attempting to bridge the gap between simulation code developers and in-depth numerical verification tools. As such, verrou aims at providing the following features:

• Low entry cost: no need for source code instrumentation, nor even recompilation [START_REF] Févotte | VERROU: L'arithmétique stochastique sans recompiler[END_REF]. To this end, verrou performs DBI using the Valgrind framework.

• Reliable detection of floating-point problems, using Monte Carlo Arithmetic (MCA).

• Easy and seamless introduction within industrial V&V processes to analyze floating-point errors.

• As low as possible run-time overhead.

The rest of this paper is organized as follows. In part 2, we briefly describe the principles of MCA, which verrou uses. The implementation of the tool is presented in part 3, where we first review the Valgrind Dynamic Binary Instrumentation framework, before showing how it is used to implement MCA. Finally, part 4 is devoted to experimental tests of verrou. Small-scale, well understood numerical applications allow validating the implementation of MCA, while the study of a large-scale, more complex industrial code showcases how verrou meets its objectives of integrating into an industrial V&V process.

Monte Carlo Arithmetic

Monte Carlo Arithmetic (MCA) is a floating-point arithmetic system in which inexact values and round-off errors are respectively modeled with randomization of the input values and the output results. Systematically replacing all floating-point operations in a program with their respective MCA counterparts therefore turns any set of program executions into a sample of a random variable representing the result and its floating-point accuracy.

In this section, we try to briefly describe the parts of MCA that are relevant to verrou. The interested reader can refer to [START_REF] Parker | Monte Carlo arithmetic: exploiting randomness in floating-point arithmetic[END_REF] for the details.

Inexact values

The MCA approach proposes to model inaccurate values as random variables. If x is an inaccurate value approaching value x to s significant digits in base β, then x can be represented as the random variable

x = x + β e-s ξ,
where e is the order of magnitude of x in base β, and ξ is a random variable on the interval [-1 2 , 1 2 ]. Conversely, knowing such a random variable, it is possible to retrieve the number s of exact significant digits of the inaccurate value it represents as

s = log β E(x) S(x)
,

where E(x) and S(x) respectively denote the expected value and the standard deviation of random variable x. Therefore, assuming that repeated executions of a program under MCA produce a sample of the (inexact) result expressed as such a random variable, it is possible to perform a statistical analysis to compute estimates μ and σ of E(x) and S(x). From there, μ can be considered to represent the computation result with a number of significant digits approximately given by

s log β μ σ . (1)
In the following, the term "inexact value" will refer to this definition of a random variable modelling a value and the inaccuracies attached to it.

Assessment of inaccuracies due to floating-point arithmetic

We now describe the Monte Carlo arithmetic which can be used by a program to output inexact values modelling the inaccuracies resulting from floating-point arithmetic. In what follows, we will suppose that all floating-point values are defined using a given, known precision (which in practice will be either the single or double binary precision of IEEE-754).

The full MCA approach can be summarized as follows: if x and y are floating point values and * is a binary operator, we introduce the MCA floating-point implementation of * , denoted by and defined by

x y = random_round value(x) * value(y) ,
where value is a function which maps values to inexact values, and is defined as

value(v) = v if v can be expressed exactly, inexact(v) otherwise,
and inexact is a function randomizing value v. The function random_round computes a random rounding of a given value:

random_round = round • value.
Using inexact in the expressions above makes the result of an MCA operation a random variable, unless operands are exact values and the result can itself be represented exactly with the considered precision. Moreover, a careful choice of the definition of inexact can lead the output result to be an inexact value in the sense of section 2.1.

Calling value on inputs x and y is known as "input precision bounding", and helps detecting catastrophic cancellations by introducing random lower-order digits which will track the loss of significance resulting from cancellation.

On the other hand, calling random_round on the result of the operation is known as "random rounding", and models forward errors. verrou does not implement the full MCA approach: it only considers random rounding, which is thus described in more details in the following.

For any value v, random_round(v) can have different definitions (corresponding to different definitions of inexact), including:

• round_nearest: rounds v to the nearest exactly representable value in the given precision. In this case, random_round is identical to the ordinary round function.

• round_random: rounds v up or down with probability 1 2 . It is similar to the random rounding mode used by the cestac method upon which cadna is built.

• round_average: rounds v up or down, with probabilities proportional to the distance from v to either of these two alternatives. This type of rounding, initially introduced and studied in [START_REF] Forsythe | Reprint of a note on rounding-off errors[END_REF], is also referred to as "uniform_absolute output randomization". As shown in [START_REF] Parker | Monte Carlo arithmetic: exploiting randomness in floating-point arithmetic[END_REF], this type of random rounding produces an inexact value with zero expected bias, i.e. it eliminates the expected error from rounding.

VERROU: floating-point accuracy assessment

We now describe the implementation of verrou. In particular, we briefly introduce Valgrind, which is used to perform Dynamic Binary Instrumentation. We then describe how MCA has been implemented within this framework. In a third part, we detail the use of verrou for numerical verification, and particularly the ways it can be introduced into an industrial V&V process. Finally, we describe more advanced instrumentation features implemented in verrou, as well as the limitations which still remain to overcome.

Dynamic binary analysis based on Valgrind

In communities of developers, Valgrind2 is a well-known tool used mostly for memory debugging. Less known is the fact that it is actually a whole framework allowing to build Dynamic Binary Instrumentation (DBI) tools [START_REF] Nethercote | Valgrind: A framework for heavyweight dynamic binary instrumentation[END_REF]. As the tool is already heavily used in numerous development and V&V toolchains, it was a logical choice for the base component upon which to build verrou.

Describing the full internal workings of Valgrind is well outside the scope of this paper. In this section, we will try to briefly describe how DBI works from the perspective of a Valgrind tool writer and user. Figure 1 presents a schematic overview of Valgrind's dynamic binary re-compilation process. On the left part of the figure is an example C source file, which has been compiled to produce a binary executable file. This will be the analyzed program. A Valgrind tool, represented by the large central box, is composed of the Valgrind core -which implements all the heavy DBI operations -and a tool plugin -which is the only part left for the tool developer to implement. For example, Valgrind's default and most popular tool is Memcheck,3 a memory error detector which is often confused with Valgrind itself.

Valgrind uses a disassemble-and-resynthetize strategy, which means that the first task performed by the core consists in converting machine code from the instrumented binary to an Intermediate Representation (IR, represented as the left column and labelled "input" in the Valgrind box on figure 1). The tool plugin is then responsible to transform the original IR into an instrumented IR (labeled "output" on the figure). In the end, the IR is converted back to machine code to be executed.

In the example given on figure 1, the first line in the C code is converted to an integer instruction in the IR, which is left untouched by the tool plugin (it is copied to the output IR). The second line in the C code is converted to two floating-point instructions in the IR. Each of these instructions will be erased, and replaced by other instructions in the output IR. In this example, custom functions are called to count instrumented floating-point instructions (pure side effect), and replace the original instruction with another implementation (which means putting the desired result back where it would have been).

Implementation of MCA operations

verrou is built as a Valgrind tool, which allows forgetting details about DBI and focusing instead on the implementation of MCA operations. However, developing within Valgrind has its own constraints, most of which are related to the fact that the Valgrind code will be executed in the same process as the instrumented code. In particular, two main limitations are impacting the development of verrou:

1. No calls to external libraries are allowed (even the C standard library is not available; the Valgrind core proposes its own implementation of a subset of the C library). This makes it difficult to re-use available MCA implementations such as MCALIB [START_REF] Frechtling | MCALIB: Measuring sensitivyty to rounding error with monte carlo programming[END_REF], since they would need to be made compatible with Valgrind's constraints, which in turn would require adapting all their dependencies to Valgrind. Among these dependencies is for example MPFR [START_REF] Fousse | MPFR: a multiple-precision binary floatingpoint library with correct rounding[END_REF], whose adaptation is feasible [START_REF] Benz | A dynamic program analysis to find floating-point accuracy problems[END_REF] but extremely heavy.

2. As for floating-point operations, only rounding to the nearest is supported. This prevents changing the CPU rounding mode to easily implement the round_random output randomization as proposed in [START_REF] Kahan | The improbability of probabilistic error analyses for numerical computations[END_REF].

MCA operations in verrou are therefore reimplemented from scratch, using as few dependencies as possible. The approach retained is based on error-free transformations (EFT). Following the notations and definitions of section 2, given two floating-point operands x and y, the objective is to compute x y = random_round(x * y), using only rounded to nearest floating point operations. Supposing that there exists an error-free transformation for binary operator * , one can write

x * y = round(x * y) + δ, (2)
where round(x * y) is the rounded result of the floating-point operation, and δ is a floating-point value which can be computed exactly using standard floating-point arithmetic. In the following, let us denote by z the exact operation result, and z , z and z its approximate floating-point values respectively rounded downwards, to nearest, and upwards as defined by the IEEE-754 standard. Knowing z, the sign of δ allows determining the downwards and upwards-rounded values of z:

if δ > 0, if δ < 0, z = z, z = z + ulp(z), z = z -ulp(z), z = z,
where ulp(z) represents the difference between two consecutive representable floating-point numbers for values of the order of magnitude of z with the considered precision.

The choice of an error-free transformation in equation ( 2) must be made in accordance with the transformed operator, and should ideally be as cost-effective as possible. In verrou, Knuth's twosum algorithm [START_REF] Knuth | The Art of Computer Programming[END_REF] is used for the sum and subtraction operations, and Dekker's twoprod [START_REF] Dekker | A floating-point technique for extending the available precision[END_REF] for the product. The case of the division is a bit more difficult to handle and is developed in further details in the next section.

The MCA result is then computed as follows: if δ = 0, then z = z is returned. Otherwise, the result is chosen according to a user-determined strategy:

• NEAREST, DOWNWARD, UPWARD: the relevant IEEE-754 rounding mode is emulated, and z, z or z is returned accordingly.

• RANDOM: this mode corresponds to the round_random strategy presented above. z or z is returned randomly, with equiprobability.

• AVERAGE: this mode corresponds to the round_average strategy presented above. z or z is returned randomly, with probabilities given by:

if δ > 0, if δ < 0, p( z ) = 1 - δ ulp(z) , p( z ) = δ ulp(z) , p( z ) = δ ulp(z) , p( z ) = 1 - δ ulp(z) .

Implementation of the division

If a and b are two floating-point numbers, we try to develop in this section a transformation without too much error for the division operation:

a b q + r, (3)
with q = round(a/b).

Our objective here is thus to evaluate, at least approximately, the remainder r, so that it can be used to compute an MCA version of the division operation, in the same way as value δ in the previous section. In order to evaluate an approximation to r, we propose algorithm 1.

Input: a, b

Output: r such that a/b round(a/b) + r. 1 q ← round(a/b) ; 2 (p, s) ← twoprod(b, q) ; 3 t ← round(a -p) ;

Let us try and compute b r. Since floating-point numbers q and p are computed with exact rounding,

q = a b (1 + 1 )
where

| 1 | < m , p = b q (1 + 2 ) where | 2 | < m ,
where m is machine precision. This implies

p = a (1 + 1 ) (1 + 2 ) ,
and thus p 2 a 2 p.

From Sterbenz' lemma, a -p is a floating-point number, so that

t := round(a -p) = a -p.
Injecting this result in equation ( 4) yields b r = t -s, which will then be computed with exact rounding:

u := round(t -s) = round(b r) = b r (1 + 3 ) where | 3 | < m ,
which in turns yields

r := round(u/b) = u b (1 + 4 ) where | 4 | < m , = r (1 + 3 ) (1 + 4 ) .

Use and introduction in an industrial V&V process

Using Valgrind's tooling makes it very easy to run any executable binary within an MCA context using verrou. The command line only needs a slight modification:

valgrind --tool=verrou --rounding-mode=random EXECUTABLE [ARGS]

When run within this MCA context, the instrumented program executes in a normal way, except output results are inexact values samples. verrou also appends its own output, indicating how many floating-point instructions were detected and successfully instrumented.

As underlined in the introduction, it is of crucial importance that verrou be compatible with the strong requirements of industrial V&V processes.

Figure 2 presents a schematic view of non-regression testing within a typical industrial Verification & Validation process. The code is compiled from sources, then run on any number of test cases to produce results. Non-regression testing tools are then used to compare these results with previously stored reference results, and produce discrepancy measurements.

As seen above, in order to perform numerical validation of a code using verrou, two major steps are needed:

1. run the code multiple times under MCA using verrou;

2. perform a statistical analysis of the results sample to produce accuracy estimates.

verrou itself allows easily tackling the problem of running a program under MCA. However, it does not help with the statistical analysis, which is entirely left for the user to implement. While this might seem to be a simple enough task for small-scale codes, it is often not so easy in the case of large-scale, industrial codes, who routinely produce several gigabytes of results data on each run.

Fortunately, large-scale, industrial codes often have non-regression testing suites as described above. And it is precisely the aim of such tools to be able to compare several executions of the computing code and produce discrepancy measurements. Figure 3 presents two ways of introducing verrou into a non-regression testing process to allow for numerical verification. On the left (figure 3a), the regression test suite is simply run in the execution context of verrou, so that results become MCA inexact values. Non-regression testing tools are then used as usual to compare

Non-regression testing tools execution compilation

Discrepancies

Reference results perturbed results to reference results which are obtained through a standard execution of the code. The measured discrepancies give an indication as to whether accumulated errors in each test case of the suite are acceptable. The major advantage of such a scenario is its ease of implementation: no modification is needed. However, its main drawback is the lack of statistics. This can be mitigated by running the test suite several times, or relying on the fact that the test suite will be run every day or so during all the code development phase.

✓ OK x KO

Figure 3b presents a slightly more complex way of introducing verrou in an industrial non-regression testing process. In this case, each test case is run several times within the MCA execution context provided by verrou. Non-regression testing tools need to be modified to compare all MCA results between them and produce statistics, instead of relying on a single reference result. This produces more reliable estimates for the accuracy of the results.

Advanced instrumentation features

The main advantage of Valgrind and verrou lies in the fact that no source modification is required. This, however, comes with a downside: it is difficult to relate an unexpected behavior to its origins in the source code (when they are available). To partially mitigate this problem, we use Valgrind's client request mechanism: the user can annotate parts of the code to generate dummy instructions in the executable. These instructions will then be interpreted by verrou during the run, but have no effect when the binary is executed outside Valgrind.

A first feature taking advantage of this mechanism is the possibility to activate or deactivate MCA to delimit sections in the program which are run in standard floating-point arithmetic. This feature is very useful in practice, since it allows roughly locating instability sources by trial and error. Another benefit offered by this feature is to avoid instabilities which are not caused by round-off errors per se, but rather by the stochastic approach itself. One such problem arises when a pure function is called twice with the same data input. Under standard floating-point arithmetic, whatever round-off errors happen in the code, the same result is always computed, which is a property on which developers might rely. However, determinism disappears under MCA, which may lead to incorrect program execution.

Rewriting programs so that they don't rely on determinism is not always a good option in terms of performance. On the other hand, although locally deactivating MCA restores determinism, it also prevents real instabilities in the function to be detected. Therefore, in order to correctly account for deterministic constraints, another client request in verrou allows annotating functions so that they compute deterministic results, but still propagate round-off errors using MCA. This feature is implemented by seeding the random generator in a deterministic way on entry of such functions.

Delta-debugging and localization of errors

The features mentioned above allow to implement coarse-grained localization of the sources of floating-point errors by performing a bisection on the portions of the binary which are instrumented. Such a search is best conducted using specialized procedures like the Delta-Debugging algorithm [START_REF] Zeller | Yesterday, My Program Worked. Today, It Does Not. Why?[END_REF][START_REF] Zeller | Simplifying and isolating failure-inducing input[END_REF][START_REF] Zeller | Why Programs Fail[END_REF].

verrou comes with a utility tool, called verrou-DD, which simplifies the localization of errors by automating the following process:

1. A first run of the test case is performed, without instrumentation. Results of this run are considered valid. A list of all functions (or rather, symbols) encountered during the run is generated. This list will become the list of code parts to check for floating-point errors (so-called "deltas", in the Delta-Debugging terminology).

2. The DD_max variant of the Delta-Debugging algorithm is performed on the list of symbols. The output is a maximal set of failure-inducing symbols, in the sense that perturbating floating-point rounding-modes in the instructions of each of these symbols produces failing results. Conversely, perturbating floating-point rounding modes in all instructions of the symbols not listed in the DD_max set is guaranteed to yield a passing result.

3. A second Delta-Debugging pass is performed at the grain-size of source code lines: for each symbol listed in the DD_max set, the following steps are performed:

(a) A first reference run is performed to generate the list of source code lines encountered in the symbol.

(b) The DD_max algorithm is performed on the list of source code lines. As for the first pass, this generates a maximal set of failure-inducing source code lines.

Limitations of verrou

Calling a mathematical function such as cos is either translated to a specific assembly operation, or is delegated to a software library. In the former case, verrou sees the special instruction, which it does not handle yet. MCA is not applied in this case. In the latter case, all verrou sees is a series of floating-point computations implemented by the mathematics library. However, these algorithms have been designed only for the nearest rounding mode and are not compatible with MCA. For example, with random rounding mode in double precision, cos may return a value greater than 42. Implementing a way to automatically circumvent this problem is one of the most important perspectives for verrou. But for now, users have to use client requests to deactivate MCA during each call to a trigonometric function.

For each floating point operation (+,-,*,/), three types of instructions can be found in an executable binary file:

1. scalar instructions of the SSE set;

2. SIMD instructions of the SSE and AVX sets;

3. scalar instructions of the x87 set.

SSE scalar instructions are fully instrumented by verrou. However, only double precision vector instructions of the SSE set are currently handled. These limitations should be removed in the short term.

Valgrind -and consequently verrou-is unable to properly instrument x87 scalar instructions. Nowadays, these instructions are no more generated by the compiler,4 except for 80-bit instructions, 5 which are used for long double variables. With verrou, a long double is treated as a double (but instructions are counted separately, so that the user is able to know exactly what has been instrumented).

Experimental evaluation

In this part, we evaluate verrou in various contexts. Well-understood numerical benchmarks are first used to validate our implementation of MCA. In a second step, we present the use of verrou in the real-world numerical verification of an industrial simulation code.

Numerical benchmarks

We first assess the usability of verrou on a few classical test cases, which are well enough understood for a complete analysis to be performed. Following Kahan [START_REF] Kahan | A survey of error analysis[END_REF], we first consider trinomial a x 2 -b x + c with a = 7169, b = -8686 and c = 2631, and define its roots

{r 1 , r 2 } = -b ± √ b 2 -4 a c 2 a .
Exact solutions to this problem are reproduced in table 1, along with floating-point results obtained using a simple C program. The table shows results from a standard execution in IEEE-754 single precision (float) with rounding to nearest, as well as results obtained with the same program (i.e. the exact same executable binary) running under single-precision MCA using verrou. In this latter case, statistics are collected from 5 samples, and bold digits are considered exact using estimation [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF]. For reference, results obtained with cadna are also shown in the table. In this case, the program sources have been instrumented, and only exact significant digits (as evaluated by cadna) are shown.

It is interesting to note that even for such a simple program, performing as few as 14 binary operations and 2 square roots calculations, a few MCA samples are enough to detect the occurrence of a problem. Obtaining a correct and reproducible estimation of inaccuracies would however have required far more sampling.

A second test case is the second-order recurrent sequence (u k ) defined by [START_REF] Bajard | Introduction to the special issue "Real Numbers and Computers[END_REF]

u 0 = 2, u 1 = -4, u k+1 = 111 - 1130 u k + 3000 u k u k-1
, ∀k > 0.

Exact calculations show that lim k u k = 6. However, when computing iterates of this sequence using IEEE-754 floatingpoint arithmetic with default rounding to nearest, we get very different results, converging to the value 100. Table 2 presents results from such a computation, when run under single precision MCA using verrou. The originality of this experience comes from the fact that the instrumented program was written in Python. No recompilation of any sort has had to take place: the verrou command line was simply prepended to the call to the Python interpreter to produce the results shown.

Statistics are collected from 5 samples. Using formula [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF], one can estimate that digits written in bold face are exact. Rows in slanted face in the table indicate intermediate iterates for which the estimated number of exact significant digits under MCA is less than 1. These results show that, when observing only a late iterate (such as u 20 ), MCA does not allow to detect any problem. It is interesting to note here that, although running a similar program under Discrete Stochastic Arithmetic (DSA) using cadna would have yielded essentially the same results, the self-validation feature would at least have warned about the (in)validity of the results.

However, the problem appears very clearly when observing all iterates: the standard deviation starts increasing rapidly, before vanishing as the sequence converges to its attractive fixed point at 100.

Verification & Validation of an industrial code

Ultra-sonic non-destructive testing is widely used to control steel welds. Athena [START_REF] Chassignole | Ultrasonic and structural characterization of anisotropic austenitic stainless steel welds: Towards a higher reliability in ultrasonic non-destructive testing[END_REF] is a computing code designed to simulate the wave propagation inside the transducer with the Stepanishen equation [START_REF] Stepanishen | Transient radiation from pistons in an infinite planar baffle[END_REF] and in the steel with a finite element method [START_REF] Becache | Fictitious domains, mixed finite elements and perfectly matched layers for 2-d elastic wave propagation[END_REF]. After 15 years of development, Athena-2D (resp. Athena-3D) is now a complex software with 27 (resp. 55) thousand lines of code (95% of Fortran), depending on external libraries Blas, Lapack, MPI (3D only), and Umfpack (3D only). As the V&V process is crucial for Athena, a large number of comparisons to experimental results have been performed and non-regression testing tools have been put in place. However, due to legacy code and external libraries, it is far from easy to instrument sources using cadna or to perform comparisons with higher precision results.

Using verrou required only a few modifications to Athena: some small parts had to be uninstrumented using client requests in order to avoid breaking pre-conditions required by the algorithm (such as the |cos(x)| > 1 or non-determinism problems mentioned in sections 3.5 and 3.7).

Following the scheme described in figure 3a, four executions (two with RANDOM and two with AVERAGE mode) of a short selection of the non-regression test suite of Athena-2D are reported in Table 3. Each output file contains a vector of results, of which the relative error is reported in 2-norm. Result files which are binary equal for all executions do not appear in the table. For example, Case6 and Case8 do not reveal any discrepancy. These two test cases output their results in ascii format with 4 decimal places, which eliminates numerical noise less than 10 -4 . While such treatments might be useful in some instances, they should be avoided when implementing accuracy assessment as described in figure 3b.

Case0, Case1, Case3, Case4, Case7 and Case10 reveal differences in the output files but the discrepancies are negligible (inferior to 10 -5 relative difference). Case2 and Case5 reveal significant discrepancies which have to be analyzed. A close look of the execution traces and the output results shows that a threshold effect in the coupling between the transducer and the steel leads to a time shift (which can be easily visualized with existing tools: cf. figure 4). Since the amplitude of signals and the time between echos are not modified, users of Athena are not concerned by these discrepancies. However, developers are more disturbed as they have to manually check and understand such cases. To avoid this fastidious work, Athena's developers can improve the transducer-steel coupling to avoid threshold effects, or develop norms adapted to user requirements. The first approach seems easier.

To be applicable in an industrial process, MCA also has to feature a reasonable run-time overhead. For Athena, compiled as used in production with gfortran -O3, the slow-down for one sample is usually between ×12 and ×40. The run-time overhead of Case8 is surprisingly high in AVERAGE rounding mode, and even higher in RANDOM rounding mode. This case has two specificities. First there is an intensive use of trigonometric functions and, consequently, of client requests. Second, small matrices are assembled and stored only if the difference with neighboring matrices is not too small. As the criterion is conservative, around 40 times more matrices are stored and factorized in the RANDOM or AVERAGE rounding modes than in the default rounding to the nearest.

In Athena-3D, only 89% of floating point operations are instrumented due to the presence of single-precision SSE instructions. We therefore don't report numerical results in this paper. However, the experiment is still useful to demonstrate that external libraries such as MPI and Umfpack are taken into account without any difficulties, and that larger test cases (10 7 cells and 600 time steps) can be run with similar run-time overheads (×25).

The Athena experiment demonstrates that verrou and non-regression tools can really easily be used together. Still, a strong knowledge of the industrial code is required to be able to analyze the results.

Conclusions

In this paper, we have shown that the MCA approach is adapted to diagnose inaccuracies due to floating-point arithmetic in numerical simulations. We proposed verrou, an MCA implementation built upon the Valgrind framework, which does not require re-compiling the code nor its external dependencies. We showed that verrou was suitable for the analysis of both small scale programs and large scale simulation codes. Thanks to a low entry cost, this new tool will be easier to promote in the industry as a preliminary diagnostic tool allowing detection of numerical errors at a coarse grain. This can then lead users to extract smaller parts of the application which reproduce the problem, and analyze them in more depth, either with verrou itself, with tools previously cited in introduction such as cadna, or with other arithmetics (Rational Number or Interval Arithmetic).

Figure 1 :

 1 Figure 1: Schematic description of dynamic binary analysis with Valgrind

4 u 1 :

 41 ← round(t -s) ; 5 r ← round(u/b) ; Algorithm Approximate evaluation of the remainder in a floating-point division. Multiplying equation (3) by b yields a = b q + b r.An algorithm like twoprod allows computing an EFT for the product b q: b q = p + s, where p = round(b q). Thus, we obtain b r = a -p -s.

Figure 2 :

 2 Figure 2: Schematic view of an industrial non-regression testing process.

Figure 3 :

 3 Figure 3: Introduction of verrou into an industrial Verification & Validation process: non-regression testing tools can be reused for numerical verification.

Figure 4 :

 4 Figure 4: Comparaison of an ascan file for each execution of Case2.

Table 1 :

 1 Roots of 7169x 2 -8686x + 2631, computed by different means.

		Computation method		r 1	r 2
		exact (rounded)		0.6062438663 0.6053616575
		IEEE-754 (float, nearest)	0.6061973	0.6054083
		error		0.0000466	0.0000466
		verrou average (5 samples) 0.6062421	0.6053803
		standard deviation	0.0000397	0.0000228
		average error		0.0000018	0.0000186
		cadna (float_st)		0.6062	0.6053
		error		0.0000439	0.0000617
	k u k average	k u k average	k	u k average
	0	2.000000	6	6.805962	12 100.777983
	1 -4.000000	7	6.580517	13 100.047690
	2 18.500000	8	6.265057	14 100.002459
	3 4 5	9.378378 7.801148 7.154356	9 10 -83.174968 3.400501 11 114.316190	. . . 20 100.000000 . . .

Table 2 :

 2 Initial iterates of sequence u k computed by a Python program running under single precision MCA (5 samples) with verrou.

Table 3 :

 3 Output of the non-regression test suite of Athena2D for 2 runs with RANDOM mode and 2 with AVERAGE mode.

		Reference Time	RANDOM 1	RANDOM 2	AVERAGE 1	AVERAGE 2
	Case0	0.55s	13.89s (×25)	13.88s (×25)	15.04s (×27)	15.05s (×27)
	ascan.dat		0.000000	0.000000	0.000000	0.000000
	Case1	4.75s	91.44s (×19)	91.12s (×19)	96.75s (×20)	95.94s (×20)
	ins1.dat		0.000006	0.000006	0.000006	0.000006
	ascan.dat		0.000000	0.000000	0.000000	0.000000
	Case2	0.45s	14.17s (×31)	14.20s (×31)	11.88s (×26)	11.79s (×26)
	ins1.dat		0.280809	0.355187	0.105049	0.000000
	enerloc.dat		0.234995	0.227002	0.086034	0.000000
	sismo.dat		0.160344	0.196354	0.008062	0.000000
	ascan.dat		0.290362	0.194622	0.344207	0.016042
	Case3	11.33s	316.73s (×27)	315.42s (×27)	237.13s (×20)	237.52s (×20)
	ascan.dat		0.000000	0.000000	0.000000	0.000000
	ins2.dat		0.000000	0.000000	0.000000	0.000000
	Case4	29.91s	1048.10s (×35) 1052.27s (×35) 1120.39s (×37) 1114.05s (×37)
	sismo.dat		0.000000	0.000000	0.000000	0.000000
	ascan.dat		0.000000	0.000000	0.000000	0.000000
	Case5	10.48s	184.01s (×17)	183.46s (×17)	194.64s (×18)	193.14s (×18)
	ins1.dat		0.000000	0.232258	0.000000	0.000000
	enerloc.dat		0.000000	0.055204	0.000000	0.000000
	ascan.dat		0.000000	0.244644	0.000000	0.000000
	Case6	2.04s	24.97s (×12)	25.44s (×12)	26.67s (×13)	26.79s (×13)
	Case7	0.80s	17.97s (×22)	18.11s (×22)	18.95s (×23)	18.76s (×23)
	ins1.dat		0.000000	0.000000	0.000000	0.000000
	ascan.dat		0.000000	0.000000	0.000000	0.000000
	Case8	0.39s	147.77s (×379) 144.15s (×370) 60.65s (×155)	59.93s (×154)
	Case9	10.31s	149.74s (×14)	147.56s (×14)	154.88s (×15)	155.58s (×15)
	ascan.dat		0.000000	0.000000	0.000000	0.000000
	Case10	21.28s	352.96s (×16)	347.29s (×16)	378.13s (×17)	374.82s (×17)
	ins1.dat		0.000000	0.000005	0.000000	0.000005
	sismo.dat		0.000000	0.000000	0.000000	0.000000
	ascan.dat		0.000000	0.000000	0.000000	0.000000

The name verrou was chosen as a tribute to cadna, since cadenas and verrou are two French words respectively translating to "padlock" and "bolt".

Valgrind: http://www.valgrind.org/

Memcheck: http://valgrind.org/docs/manual/mc-manual.html

unless explicitly asked using switches such as -mfpmath=387

It should be noted that double rounding problems can occur with 80-bit instructions, due to the different size of the mantissa in the register and in memory.

Acknowledgement

The authors feel greatly indebted to Josef Weidendorfer, who helpfully gave invaluable advice on valgrind tools development.