Pierrick Gaudry

Alexander Kruppa

Paul Zimmermann

A GMP-based Implementation of Schönhage-Strassen's Large Integer Multiplication Algorithm

Keywords: I.1.2 [Computing methodologies]: Algorithms-Symbolic and algebraic manipulation Algorithms, Performance Integer multiplication, multiprecision arithmetic

Schönhage-Strassen's algorithm is one of the best known algorithms for multiplying large integers. Implementing it efficiently is of utmost importance, since many other algorithms rely on it as a subroutine. We present here an improved implementation, based on the one distributed within the GMP library. The following ideas and techniques were used or tried: faster arithmetic modulo 2 n + 1, improved cache locality, Mersenne transforms, Chinese Remainder Reconstruction, the √ 2 trick, Harley's and Granlund's tricks, improved tuning.

INTRODUCTION

Since Schönhage and Strassen presented in 1971 a method to multiply two N -bit integers in O(N log N log log N ) time [START_REF] Schönhage | Multiplikation großer Zahlen[END_REF], several authors have shown how to reduce other operations -inverse, division, square root, gcd, base conversion, elementary functions -to multiplication, possibly with log N multiplicative factors [START_REF] Bernstein | Fast multiplication its applications[END_REF][START_REF] Brent | Multiple-precision zero-finding methods and the complexity of elementary function evaluation[END_REF][START_REF] Knuth | The analysis of algorithms[END_REF][START_REF] Moenck | Fast modular transforms via division[END_REF][START_REF] Schönhage | Berechnung von Kettenbruchentwicklungen[END_REF][START_REF] Steel | Reduce everything to multiplication[END_REF]. It has now become common practice to express complexities in terms of the cost M (N ) to multiply two N -bit numbers, and many researchers tried hard to get the best possible constants in front of M (N ) for the above-mentioned operations (see for example [START_REF] Bernstein | Removing redundancy in high-precision Newton iteration[END_REF][START_REF] Karp | High-precision division and square root[END_REF]).

Strangely, much less effort was made for decreasing the implicit constant in M (N ) itself, although any gain on that constant will give a similar gain on all multiplication-based operations. Some authors reported on implementations of large integer arithmetic for specific hardware or as part of a number-theoretic project [START_REF] Bailey | The computation of π to 29,360,000 decimal digits using Borwein's quartically convergent algorithm[END_REF][START_REF] Crandall | Discrete weighted transforms and large-integer arithmetic[END_REF]. In this article we concentrate on the question of an optimized implementation of Schönhage-Strassen's algorithm on a classical workstation.

In the last few years, the multiplication of large integers has found several new applications in "real life", and not only in computing billions of digits of π. One such application is the segmentation method (called Kronecker substitution in [START_REF] Von Zur Gathen | Modern Computer Algebra[END_REF]) to reduce the multiplication of polynomials with integer coefficients to one huge integer multiplication; this is used for example in the GMP-ECM software [START_REF] Zimmermann | 20 years of ECM[END_REF]. Another example is the multiplication or factorization of multivariate polynomials [START_REF] Steel | Reduce everything to multiplication[END_REF][START_REF] Van Der Hoeven | The truncated Fourier transform and applications[END_REF].

In this article we detail several ideas or techniques that may be used to implement Schönhage-Strassen's algorithm (SSA) efficiently. As a consequence, we obtain what we believe is the best existing implementation of SSA on current processors; this implementation might be used as a reference to compare with other algorithms based on the Fast Fourier Transform, in particular those using complex floating-point numbers.

The paper is organized as follows: §1 revisits the original SSA and defines the notation used in the rest of the paper;

§2 describes the different ideas and techniques we tried; finally §3 provides timing figures and graphs obtained with our new GMP implementation, and compares it to other implementations.

THE ALGORITHM OF SCHÖNHAGE AND STRASSEN

Throughout the paper we use w for the computer word size in bits -usually 32 or 64 -and denote by N the number of bits of the numbers we want to multiply.

Several descriptions of SSA can be found in the literature, see [START_REF] Crandall | Prime Numbers: A Computational Perspective[END_REF][START_REF] Schönhage | Multiplikation großer Zahlen[END_REF] for example. We recall it here to establish the notations.

Let R + N -or simply RN -be the ring of integers modulo 2 N + 1. SSA reduces integer multiplication to multiplication in RN , which reduces to polynomial multiplication in Z[x] mod (x K + 1), which in turn reduces to polynomial multiplication in Rn[x] mod (x K + 1), which finally reduces to multiplication in Rn. The reason for choosing RN as the ring to map the input integers to is that the multiplications of elements of Rn can use SSA recursively, skipping the first step of mapping from integers to RN again. The first reduction -from Z to RN -is simple: to multiply two non-negative integers of u and v bits, it suffices to compute their product mod 2 N + 1 for N ≥ u + v.

Z ⇒ R N ⇒ Z[x] mod (x K + 1) ⇒ Rn[x] mod (x K + 1) ⇒ Rn
The second step -a map from RN to Z[x] mod (x K + 1) -works as follows. Assume N = 2 k M for integers k and M , and define K := 2 k . An integer a ∈ [0, 2 N ] can be uniquely written a = P K-1 i=0 ai2 iM , with 0 ≤ ai < 2 M for i < K -1, and 0 ≤ aK-1 ≤ 2 M , that is, we cut a into K pieces of M bits each, except the last piece can be equal to 2 M . Now the integer a is the value at x = 2 M of the polynomial A(x) = P K-1 i=0 aix i . Assume we decompose an integer b ∈ RN in the same manner, and let C(x) be the product

A(x)B(x) over Z[x]: C(x) = P 2K-2 i=0 cix i . One now has ab = A(2 M )B(2 M ) = C(2 M ), thus ab = P 2K-2 i=0 ci2 iM . Now what we really want is ab mod (2 N + 1), i.e., (c0 -cK)+• • •+(cK-2 -c2K-2)2 (K-2)M +cK-12 (K-1)M (1) which comes from C + (x) := P K-1 i=0 cix i = A(x)B(x) mod (x K + 1), since x = 2 M and N = KM . To determine C + (x)
, one uses a negacyclic convolution over the ring Rn, i.e., modulo 2 n + 1, where n is taken large enough so that the ci can be recovered exactly. For 0

≤ i ≤ K -1, one has 0 ≤ ci = P i j=0 ajbi-j < (i + 1)2 2M . Similarly for K ≤ i ≤ 2K -3, one has 0 ≤ ci < (2K -1 -i)2 2M
and finally 0 ≤ c2K-2 ≤ 2 2M . With the convention that c2K-1 = 0, according to (1), we have

((i + 1) -K)2 2M ≤ ci = ci -ci+K < (i + 1)2 2M (2)
for 0 ≤ i < K. Hence each coefficient of C(x) mod (x K + 1) is confined to an interval of length K2 2M , and so it suffices to have 2 n + 1 ≥ K2 2M , i.e., n ≥ 2M + k 1 . The negacyclic convolution A(x)B(x) mod (x K + 1) can be performed efficiently using the Fast Fourier Transform (FFT). More precisely, SSA uses here a simple case of the Discrete Weighted Transform (DWT) [START_REF] Crandall | Discrete weighted transforms and large-integer arithmetic[END_REF]. Assume ω = θ 2 is a primitive Kth root of unity in Rn. (All operations in this paragraph are in Rn.) Given (ai) 0≤i<K , the weight signal is (a ′ i := θ i ai) 0≤i<K . The forward transform computes (âi := P K-1 j=0 ω ij a ′ j ) 0≤i<K , and similarly for ( bi). One then multiplies âi and bi together in Rn (pointwise products): let ĉi = âi bi. The backward transform computes (c ′ i := P K-1 j=0 ω -ij ĉj)0≤i<K :

c ′ i = K-1 X j=0 ω -ij âj bj = K-1 X j=0 ω -ij K-1 X ℓ=0 ω jℓ a ′ ℓ ! K-1 X m=0 ω jm b ′ m ! = K-1 X ℓ,m=0 a ℓ bmθ ℓ+m K-1 X j=0 ω j(ℓ+m-i) .
Since ω is a primitive Kth root of unity, P K-1 j=0 ω j(ℓ+m-i) is zero unless ℓ + m -i ≡ 0 mod K, which holds for ℓ + m = i 1 One might use n ≥ 2M + k + 1 to get a lifting algorithm from Rn to Z which is independent of i.

or ℓ + m = i + K. Since θ K ≡ -1 mod (2 n + 1), it follows: c ′ i = Kθ i K-1 X ℓ=0 (a ℓ b i-ℓ -a ℓ b i+K-ℓ ) = Kθ i (ci -ci+K ),
where bm is assumed zero for m outside the range [0, K -1]. SSA thus consists of five consecutive steps, where all computations in steps (2) to (4) are done modulo 2 n + 1:

(1) the "decompose" step extracts from a the M -bit parts ai, and multiplies them by the weight signal θ i , obtaining a ′ i (similarly for bi);

(2) the "forward transform" computes (â0, . . . , âK-1) from (a ′ 0 , . . . , a ′ K-1 ) (similarly for bi);

(3) the "pointwise product" step computes ĉi = âi bi, for 0 ≤ i < K;

(4) the "backward transform" computes (c ′ 0 , . . . , c ′ K-1 ) from (ĉ0, . . . , ĉK-1);

(5) the "recompose" step divides c ′ i by 2 k θ i , and constructs the final result as c0 +c12 M +• • •+cK-12 (K-1)M . Some ci, defined in Eq. ( 2), may be negative, but the sum is necessarily non-negative.

For a given input bit-size N , several choices of the FFT length K may be possible. SSA is thus a whole family of algorithms: we call FFT-K -or FFT-2 k -the algorithm splitting the inputs into K = 2 k parts. For a given input size N , one of the main practical problems is how to choose the best value of the FFT length K, and thus of the bit-size n of the smaller multiplies (see §2.6).

Choice of n and Efficiency

SSA takes for n a multiple of K, so that ω = 2 2n/K is a primitive Kth root of unity, and θ = 2 n/K is used for the weight signal. This ensures that all FFT butterflies only involve additions/subtractions and shifts on a radix 2 computer (see §2.1).

In practice one may additionally require n to be a multiple of the word size w, to make the arithmetic in 2 n + 1 simpler. Indeed, a number from Rn is then represented by n/w machine words, plus one additional bit of weight 2 n . We call this a semi-normalized representation, since values up to 2 n+1 -1 can be represented.

For a given bit size N divisible by K = 2 k , we define the efficiency of the FFT-K scheme:

2N/K + k n ,
where n is the smallest multiple of K larger than or equal to 2N/K + k. For example for N = 1, 000, 448 and K = 2 10 , we have 2N/K + k = 1964, and the next multiple of K is n = 2048, therefore the efficiency is 1964 2048 ≈ 96%. For N = 1, 044, 480 with the same value of K, we have 2N/K + k = 2050, and the next multiple of K is n = 3072, with an efficiency of about 67%. The FFT scheme is close to optimal when its efficiency is near 100%.

Note that a scheme with efficiency below 50% does not need to be considered. Indeed, this means that 2N/K + k ≤ 1 2 n, which necessarily implies that n = K (remember n has to be divisible by K). Then the FFT scheme of length K/2 can be performed with the same value of n, since 2(N/(K/2)) + (k -1) < 4N/K + 2k ≤ n, and n is a multiple of K/2.

From this last remark, we can assume 2N/K ≥ 1 2 nneglecting the small k term -, which together with n ≥ K gives:

K ≤ 2 √ N . (3)

OUR IMPROVEMENTS

We describe in this section the ideas and techniques we have tried to improve the GMP implementation of SSA. We started from the GMP 4.2.1 implementation, and used the graph of the multiplication time up to 1, 000, 000 words on an Opteron as benchmark. After encoding each idea, if the new graph was better than the old one, the new idea was validated, otherwise it was discarded. Each technique saved only 5% up to 20%, but all techniques together saved a factor of about 2 with respect to GMP 4.2.1.

Arithmetic Modulo 2 n + 1

Arithmetic operations modulo 2 n +1 have to be performed during the forward and backward transforms, when applying the weight signal, and when unapplying it. Thanks to the fact that the primitive roots of unity are powers of two, the only needed operations are additions, subtractions, and multiplications by a power of two. Divisions by 2 k can be reduced to multiplications by 2 2n-k .

We recall that we desire n to be a multiple of the number w of bits per word. Since n must also be a multiple of K = 2 k , this is not a real constraint, unless k < 5 on a 32bit computer, or k < 6 on a 64-bit computer. Let m = n/w be the number of computer words corresponding to an nbit number. A residue mod 2 n + 1 has a semi-normalized representation with m full words and one carry of weight 2 n : a = (am, am-1, . . . , a0), with 0 ≤ ai < 2 w for 0 ≤ i < m, and 0 ≤ am ≤ 1.

The addition of two such representations is done as follows (we give here the GMP code):

c = a[m] + b[m] + mpn_add_n (r, a, b, m); r[m] = (r[0] < c); MPN_DECR_U (r, m + 1, c -r[m]);
The first line adds (am-1, . . . , a0) with (bm-1, . . . , b0), puts the low m words of the result in (rm-1, . . . , r0), and adds the out carry to am + bm; we thus have 0 ≤ c ≤ 3. The second line yields rm = 0 if r0 ≥ c, in which case we simply subtract c from r0 at the third line. Otherwise rm = 1, and we subtract c -1 from r0: a borrow may propagate, but at most to rm. In all cases r = a + b mod (2 n + 1), and r is semi-normalized. The subtraction is done in a similar manner.

The multiplication by 2 e is more tricky to implement. However this operation mainly appears in the butterflies [a, t] ← [a + b, (a -b)2 e ] of the forward and backward transforms, which may be performed as follows: Step 3 means that the most significant words from t are formed with al -bl, and the least significant words with bh -ah, where we assume that borrows are propagated, so that t is semi-normalized. Thus the only real multiplication by a power of two is that of step 5, which may be efficiently performed with GMP's mpn_lshift routine.

If one has a combined addsub routine which computes simultaneously x + y and x -y faster than two separate calls, then step 4 can be written a <-(bh + ah, al + bl), which shows that t and a may be computed with two addsub calls.

Cache Locality During the Transforms

When multiplying large integers with SSA, the time spent in accessing data for performing the Fourier transforms is non-negligible. The literature is rich with papers dealing with the organization of the computations in order to improve the locality. However most of these papers are concerned with contexts which are different from ours: usually the coefficients are small and most often they are complex numbers represented as a pair of double's. Also there is a variety of target platforms, from embedded hardware implementations to super-scalar computers.

We have tried to apply several of these approaches in our context where the coefficients are modular integers that fit in at least a few cache lines and the target platform is a standard PC workstation.

In this work, we concentrate on multiplying large, but not huge integers. By this we mean that we consider only 3 levels of memory for our data: L1 cache, L2 cache, and standard RAM. In the future we might consider also the case where we have to use the hard disk as a 4th level.

Here are the orders of magnitude for these memories, to fix ideas: on a typical Opteron, a cache line is 64 bytes; the L1 data cache is 64 kB; the L2 cache is 1 MB; the RAM is 8 GB. The smallest coefficient size (i.e., n-bit residues) we consider is about 50 machine words, that is 400 bytes. For very large integers, a single coefficient hardly fits in the L1 cache.

The very first FFT algorithm is the iterative one. In our context this is a really bad idea. The main advantage of it is that the data is accessed in a sequential way. In the case where the coefficients are small enough so that several of them fit in a cache line, this saves many cache misses. But in our case, contiguity is irrelevant due to the size of the coefficients compared to cache lines.

The next very classical FFT algorithm is the recursive one. In this algorithm, at a certain level of recursion, we work on a small set of coefficients, so that they must fit in the cache. This version (or a variant of it) was implemented in GMP up to version 4.2.1. This behaves well for moderate sizes, but when multiplying large numbers, everything fits in the cache only at the tail of the recursion, so that most of the transform is already done when we are at last in the cache. The problem is that before getting to the appropriate recursion level, the accesses are very cache unfriendly.

In order to improve the locality for large transforms, we have tried three strategies found in the literature: the Belgian approach, the radix-2 k transform, and Bailey's 4-step algorithm. Step 2

Step 3 a7 a3 a5 a1 a6 a2 a4 a0

Figure 1: FFT circuit of length 8 and butterfly tree of depth 3.

The Belgian Transform

In [START_REF] Brockmeyer | Parametrizable behavioral IP module for a data-localized low-power FFT[END_REF], Brockmeyer et al. propose a way of organizing the transform that reduces cache misses. In order to explain it, let us first define a tree of butterflies as follows (we don't mention the root of unity for simplicity):

TreeBfy(A, index, depth, stride) Bfy(A[index], A[index+stride]) if depth > 1 TreeBfy(A, index-stride/2, depth-1, stride/2) TreeBfy(A, index+stride/2, depth-1, stride/2)
An example of a tree of depth 3 is given on the right of Figure 1. Now, the depth of a butterfly tree is bounded by a value that is not the same for every tree. For instance, on Figure 1, the butterfly tree that starts with the butterfly between a0 and a4 has depth 1: one can not continue the tree on step 2. Similarly, the tree starting with the butterfly between a1 and a5 has depth 1, the tree starting between a2 and a6 has depth 2 and the tree starting between a3 and a7 has depth 3. More generally, the depth can be computed by a simple formula.

One can check that by considering all the trees of butterflies starting with an operation at step 1, we cover the complete FFT circuit. It remains to find the right ordering for computing those trees of butterflies. For instance, in the example of Figure 1, it is important to do the tree that starts between a3 and a7 in the end, since it requires data from all the other trees.

One solution is to perform the trees of butterflies following the BitReverse order. For an integer i whose binary representation fits in at most k bits, the value BitReverse(i,k) is the integer one obtains by reading the k bits (maybe padded with zeros) in the opposite order. One obtains the following algorithm, where ord_2 stands for the number of trailing zeros in the binary representation of an integer (together with the 4-line TreeBfy routine, this is a recursive description of the 36-line routine from [9, Code 6.1]):

BelgianFFT(A, k) K = 2^{k-1} for i := 0 to K-1 TreeBfy(A, BitReverse(i, k-1), 1+ord_2(i+1), K)
Inside a tree of butterflies, we see that most of the time, the butterfly operation will involve a coefficient that has been used just before, so that it should still be in the cache. Therefore an approximate 50% cache-hit is provided by construction, and we can hope for more if the data is not too large compared to the cache size.

We have implemented this in GMP, and this saved a few percent for large sizes, thus confirming the fact that this approach is better than the classical recursive transform.

Higher Radix Transforms

The principle of higher radix transforms is to use an atomic operation which groups several butterflies. In the book [START_REF] Arndt | Algorithms for programmers[END_REF] the reader will find a description of several variants in this spirit. The classical FFT can be viewed as a radix-2 transform. The next step is a radix-4 transform, where the atomic operation has 4 inputs and 4 outputs (without counting roots of unity) and groups 4 butterflies of 2 consecutive steps of the FFT.

We can then build a recursive algorithm upon this atomic operation. Of course, since we perform 2 steps at a time, the number of steps in the recursion is reduced by a factor of 2, and we have to handle separately the last step when the FFT level k is odd.

In the literature, the main interest for higher radix transforms comes from the fact that the number of operations is reduced for a transform of complex numbers (this is done by exhibiting a free multiplication by i). In our case, the number of operations remains the same. However, in the atomic block each input is used in two butterflies, so that the number of cache misses is less than 50%, just as for the Belgian approach. Furthermore, with the recursive structure, just as for the classical recursive FFT, at some point we deal with a number of inputs which is small enough so that everything fits in the cache.

We have tested this approach, and this was faster than the Belgian transform by a few percent.

The next step after radix 4 is radix 8 which works in the same spirit, but grouping 3 levels at a time. We have also implemented it, but this saved nothing, and was even sometimes slower than the radix 4 approach. Our explanation is that for small numbers, radix 4 is close to optimal with respect to cache locality, and for large numbers, the number of coefficients that fit in the cache is rather small and we have misses inside the atomic block of 12 butterflies. Further investigation is needed to validate this explanation.

More generally, radix 2 t groups t levels together, with a total of t2 t-1 butterflies, over 2 t residues. If all those residues fit in the cache, the cache miss rate is less than 1/t. Thus the optimal strategy seems to choose for t the largest integer such that 2 t n bits fit in the cache (either L1 or L2, in fact the smallest cache where a single radix 2 butterfly fits).

Bailey's 4-step Algorithm

The algorithm we describe in this section can be found in a paper by Bailey [START_REF] Bailey | FFTs in external or hierarchical memory[END_REF]. In there, the reader will find earlier references tracing back the original idea. For simplicity we stick to the "Bailey's algorithm" denomination.

A way of seeing Bailey's 4-step transform algorithm is as a radix-√ K transform, where K = 2 k is the length of the input sequence. In other words, instead of grouping 2 steps as in radix-4, we group k/2 steps. To be more general, let us write k = k1 + k2, where k1 and k2 are to be thought as close to k/2, but this is not really necessary. Then Bailey's 4-step algorithm consists in the following phases:

1. Perform 2 k 2 transforms of length 2 k 1 ; 2. Multiply the data by weights; 3. Perform 2 k 1 transforms of length 2 k 2 .

There are only three phases in this description. The fourth phase is usually some matrix transposition, but this is irrelevant in our case: the coefficients are large so that we keep a table of pointers to them, and this transposition is just pointer exchanges which are basically for free, and fit very well in the cache.

The second step involving weights is due to the fact that in the usual description of Bailey's 4-step algorithm, the transforms of length 2 k 1 are exactly Fourier transforms, whereas the needed operation is a twisted Fourier transform where the roots of unity involved in the butterflies are different (since they involve a 2 k -th root of unity, whereas the classical transform of length 2 k 1 involves a 2 k 1 -th root of unity). In the classical FFT setting this is very interesting, since we can then reuse some small-dimension implementation that has been very well optimized. In our case, we have found it better to write separate code for this twisted FFT, so that we merge the first and second phases.

The interest of this way of organizing the computation is again not due to a reduction of the number of operations, since they are exactly the same as with the other FFT approaches mentioned above. The goal is to help locality. Indeed, assume that √ K coefficients fit in the cache, then the number of cache misses is at most 2K, since each call to the internal FFT or twisted FFT operates on √ K coefficients. Of course we are interested in numbers for which √ K coefficients do not fit in the L1 cache, but for all numbers we might want to multiply, they do fit in the L2 cache. Therefore the structure of the code follows the memory hierarchy: at the top level of Bailey's algorithm, we deal with the RAM vs L2 cache locality question, then in each internal FFT or twisted FFT, we can take care of the L2 vs L1 cache locality question. This is done by using the radix-4 variant inside our Bailey-algorithm implementation.

We have implemented this approach (with a threshold for activating Bailey's algorithm only for large sizes), and combined with radix-4, this gave us our best timings. We have also tried a higher dimensional transform, in particular 3 steps of size 3 √ K. This did not help for the sizes we considered.

Mixing Several Phases

Another way to improve locality is to mix different phases of the algorithm in order to do as much work as possible on the data while they are in the cache. An easy improvement in this spirit is to mix the pointwise multiplication and the backward transform, in particular when Bailey's algorithm is used. Indeed, after the two forward transforms have been computed, one can load the data corresponding to the first column, do the pointwise multiplication of its elements, and readily perform the small transform of this column. Then the data corresponding to the second column is loaded, multiplied and transformed, and so on. In this way, one saves one full pass on the data. Taking the idea one step further, assuming that the forward transform for the first input number has been done already (or that we are squaring one number), after performing the column-wise forward transform on the second number we can immediately do the point-wise multiply and the backward transform on the column, so saving another pass over memory.

Following this idea, we can also merge the "decompose" and "recompose" steps with the transforms, again to save a pass on the data. In the case of the "decompose" step, there is more to it since one can also save unnecessary copies by merging it with the first step of the forward transform.

The "decompose" step consists of cutting parts of M bits from the input numbers, then multiplying each part ai by θ i modulo 2 n + 1, giving a ′ i . If one closely looks at the first FFT level, it will perform a butterfly between a ′ i and a ′ i+K/2 with θ 2i as multiplier. This will compute a ′ i + a ′ i+K/2 and a ′ i -a ′ i+K/2 , and multiply the latter by θ 2i . It can be seen that the M non-zero bits from a ′ i and a ′ i+K/2 do not overlap, thus no real addition or subtraction is required: the results a ′ i + a ′ i+K/2 and a ′ i -a ′ i+K/2 can be obtained with just copies and ones' complements. As a consequence, it should be possible to completely avoid the "decompose" step and the first FFT level, by directly starting from the second FFT level, which for instance will add a ′ i + a ′ i+K/2 to (a ′ j -a ′ j+K/2 )θ 2j ; here the four operands a ′ i , a ′ i+K/2 , a ′ j , a ′ j+K/2 will be directly taken from the input integer a, and the implicit multiplier θ 2j will be used to know where to add or subtract a ′ j and a ′ j+K/2 . This example illustrates the kind of savings obtained by avoiding trivial operations like copies and ones' complements, and furthermore improving the locality. This idea was not used in the results from §3.

Fermat and Mersenne Transforms

The reason why SSA uses negacyclic convolutions is because the algorithm can be used recursively: the "pointwise products" modulo 2 n + 1 can in turn be performed using the same algorithm, each one giving rise to K ′ smaller pointwise products modulo 2 n ′ + 1. (In that case, n must satisfy an additional divisibility condition related to K ′ .) A drawback of this approach is that it requires a weighted transform, i.e., additional operations before the forward transforms and after the backward transform. However, if one looks carefully, power-of-two roots of unity are needed only at the "lower level", i.e., in R + n . Therefore one can replace RN by R - Ni.e., the ring of integers modulo 2 N -1 -in the original algorithm, and replace the weighted transform by a classical cyclic convolution, to compute a product mod 2 N -1. This works only at the top level of the algorithm, and not recursively. We call this a "Mersenne transform", whereas the original SSA performs a "Fermat transform"2 . This idea of using a Mersenne transform is already present in [START_REF] Bernstein | Multidigit multiplication for mathematicians[END_REF] where it is called "cyclic Schönhage-Strassen trick".

Despite the fact that it can be used at the top level only, the Mersenne transform is nevertheless very interesting for the following reasons: • a Mersenne transform modulo 2 N -1, combined with a Fermat transform modulo 2 N + 1 and CRT reconstruction, can be used to compute a product of two N -bit integers; • a Mersenne transform can use a larger FFT length K = 2 k than the corresponding Fermat transform. Indeed, while K must divide N for the Fermat transform, so that θ = 2 N/K is a power of two, it only needs to divide 2N for the Mersenne transform, so that ω = 2 2N/K is a power of two. This improves the efficiency for K near √ N , and enables one to use a value of K close to optimal. (The constraint on the FFT length can still be decreased by using the " √ 2 trick", see §2.4.)

The above idea can be generalized to a Fermat transform mod 2 aN + 1 and a Mersenne transform mod 2 bN -1 for small integers a, b.

Lemma 1. Let a, b be two positive integers. Then at least one of gcd(2 a + 1, 2 b -1) and gcd(2 a -1, 2 b + 1) is 1.

Proof. Let g = gcd(a, b), r = 2 g , a ′ = a/g, b ′ = b/g.
Denote by ordp(r) the multiplicative order of r (mod p). In the case of b

′ odd, p | r b ′ -1 ⇒ ordp(r) | b ′ ⇒ 2 ∤ ordp(r), and p | r a ′ + 1 ⇒ ordp(r) | 2a ′ and ordp(r) ∤ a ′ ⇒ 2 | ordp(r), hence no prime can divide both r b ′ -1 and r a ′ + 1.
In the other case of b ′ even, a ′ must be odd, and the same argument holds with the roles of a ′ an b ′ exchanged, so no prime can divide both r a ′ -1 and r b ′ + 1.

It follows from Lemma 1 that we can use one Fermat transform of size aN (respectively bN ) and one Mersenne transform of size bN (respectively aN ). However this does not imply that the reconstruction is easy: in practice we used b = 1 and made only a vary (see §2.6.2).

The √ 2 Trick

Since all prime factors of 2 n + 1 are p ≡ 1 (mod 8) if 4 | n, 2 is a quadratic residue (mod n), and it turns out that √ 2 is of a simple enough form to make it useful as a root of unity with power-of-two order. Specifically,

" 2 3n/4 -2 n/4 " 2 ≡
2 (mod 2 n + 1), which is easily checked by expanding the square. Hence we can use √ 2 = 2 3n/4 -2 n/4 as a root of unity of order 2 k+2 in the transform to double the possible transform length for a given n. In the case of the negacyclic transform, this allows a length 2 k+1 transform, and √ 2 is used only in the weight signal. For a cyclic transform, √ 2 is used normally as a root of unity during the transform, allowing a transform length of 2 k+2 . This idea is mentioned in [4, §9] where it is credited without reference to Schönhage, but we have been unable to track down the original source. In our implementation, this √ 2 trick saved roughly 10% on the total time of integer multiplication.

Unfortunately using higher roots of unity for the transform is not feasible as prime divisors of 2 n + 1 are not necessarily congruent to 1 (mod 2 k+3 ), deciding whether they are or not requires factoring 2 n + 1, and even if they are as in the case of the eighth Fermat number F8 = 2 256 + 1 [START_REF] Brent | Factorization of the eighth Fermat number[END_REF], there does not seem to be a simple form for 4 √ 2 which would make it useful as a root on unity in the transform.

Harley's and Granlund's Tricks

Rob Harley [START_REF] Harley | Personal communication[END_REF] suggested the following trick3 to improve the efficiency of a given FFT scheme. Assume 2M + k is just above an integer multiple of K, say λK. Then we have to use n = (λ + 1)K, which gives an efficiency of only about λ λ+1 . Harley's idea is to use n = λK instead, and recover the missing information from a CRT-reconstruction with an additional computation modulo the machine word 2 w .

A drawback of Harley's trick is that when only a few bits are missing, the K 2 word products may become relatively expensive. When only a few bits are missing, we can multiply A(x) and B(x) over Z[x] modulo a small power of 2 using the segmentation method. That way, if h ≤ w bits are missing, one trades K 2 word products for the product of two large integers of (2h + k)K/w words, which can in turn use fast multiplication 4 .

Torbjörn Granlund [12] found that this idea -combining computations mod 2 n + 1 with computations mod 2 h -can also be used at the top-level for the plain integer multiplication, and not only at the lower-level as in Harley's trick. Assume one wants to multiply two integers u and v whose product has m bits, where m is just above an "optimal" Fermat scheme (2 N +1, K), say m = N +h. Then first compute uv mod (2 N +1), and second compute uv mod 2 h , by simply computing the plain integer product (u mod 2 h )(v mod 2 h ), again possibly in turn with fast multiplication. The exact value of uv can be efficiently reconstructed by CRT from both values. We call this idea "Granlund's trick".

Let us denote M (N ), M + (N ) and M -(N ) the cost of the multiplication of two N -bit integers, multiplication modulo 2 N + 1 and multiplication modulo 2 N -1 respectively. Granlund's trick can be written M (N + h) = M + (N ) + M (h), or M (N + h) = M + (N ) + M + (2h) if one reduces the plain product modulo 2 h to a modular product modulo 2 2h + 1. Marco Bodrato (personal communication) discovered that Granlund's trick can be applied simultaneously to the low and high ends of the product, giving M (N + h) = M + (N ) + 2M (h/2).

We use neither Harley's nor Granlund's trick in our current implementation. We believe Granlund's trick is less efficient than the generalized Fermat-Mersenne scheme we propose ( §2.3), which yields M (a + b) = M + (a) + M -(b), with M -(N ) the cost of a multiplication modulo 2 N -1, if a good efficiency is possible for M + (a) and M -(b). As for Harley's trick, we tried it only at the word level, i.e., for λK < 2M + k ≤ λK + w, which happens in rare cases only, and it made little difference. However, when multiplying two numbers modulo a Fermat number 2 2 n + 1, Harley's trick becomes very attractive, since 2M is a power of two in that case.

Improved Tuning

We found that significant speedups could be obtained with better tuning schemes, which we describe here. All examples given in this section are related to an Opteron.

Tuning the Fermat and Mersenne Transforms

Until version 4.2.1, GMP used a naive tuning scheme for the FFT multiplication. For the Fermat transforms modulo 2 N + 1, an FFT of length 2 k was used for t k ≤ N < t k+1 , where t k is the smallest bit-size for which FFT-2 k is faster than FFT-2 k-1 . For example on an Opteron, the default gmp-mparam.h file uses k = 4 for a size less than 528 machine words, then k = 5 for less than 1184 words, and so on: A special rule is used for the last entry: here k = 14 is used for less than m = 12582912 words, k = 15 is used for less than 4m = 50331648 words, and then k = 16 is used. An additional single threshold determines from which size upward -still in words -a Fermat transform mod 2 n + 1 is faster than a full product of two n-bit integers:

#define MUL_FFT_TABLE {

#define MUL_FFT_MODF_THRESHOLD 544

For a product mod 2 n + 1 of at least 544 words, GMP 4.2.1 therefore uses a Fermat transform, with k = 5 until 1183 words according to the above MUL_FFT_TABLE. Below the 544 words threshold, the algorithm used is the 3-way Toom-Cook algorithm, followed by a reduction mod 2 n + 1. This scheme is clearly not optimal since the FFT-2 k curves intersect several times, as shown by Figure 2

Tuning the Plain Integer Multiplication

Up to GMP 4.2.1, a single threshold controls the plain integer multiplication:

#define MUL_FFT_THRESHOLD 7680

This means that SSA is used for a product of two integers of at least 7680 words, which corresponds to about 148, 000 decimal digits, and the Toom-Cook 3-way algorithm is used below that threshold.

We now use the generalized Fermat-Mersenne scheme described in §2.3 with b = 1 (in our implementation we found 1 ≤ a ≤ 7 was enough). Again, for each size, the best value of a is determined by our tuning program: #define MUL_FFT_FULL_TABLE2 {{16, 1}, {4224, 2}, {4416, 6}, {4480, 2}, {4608, 4}, {4640, 2}, ...

For example, the entry {4608, 4} means that to multiply two numbers of 4608 words -or whose product has 2×4608 words -the new algorithm uses one Mersenne transform modulo 2 N -1 and one Fermat transform modulo 2 4N + 1.

Reconstruction is easy since 2 aN + 1 ≡ 2 mod (2 N -1).

RESULTS AND CONCLUSION

On July 1st, 2005, Allan Steel wrote a web page [START_REF] Steel | Magma V2.12-1 is up to 2.3 times faster than GMP 4.1.4 for large integer multiplication[END_REF] entitled "Magma V2.12-1 is up to 2.3 times faster than GMP 4.1.4 for large integer multiplication". This was actually our first motivation for improving GMP's implementation. Magma V2.13-6 takes 2.22s to multiply two numbers of 784141 words, whereas our GMP development code takes only 0.96s. Thus our GMP-based code is clearly faster than Magma by a factor of 2.3. Note that this does not mean that we have gained a factor 2. cases, 2.3 is the maximal ratio between Magma V2.12-1 and GMP 4.1.4, and between our code and Magma V2.13-6 respectively, following the well known "benchmarketing" strategy5 (both versions of Magma give very similar timings).

We have tested other freely available packages providing an implementation for large integer arithmetic. Among them, some (OpenSSL/BN, LiDiA/libI) do not go beyond Karatsuba algorithm, some do have some kind of FFT, but are not really made for really large integers: arprec, Miracl. Two useful implementations we have tested are apfloat and CLN. They take about 4 to 5 seconds on our test machine to multiply one million-word integers, whereas we need about 1 second. Bernstein mentions some partial implementation Zmult of Schönhage-Strassen's algorithm, with good timings, but right now, only very few sizes are handled, so that the comparison with our software is not really possible.

A program that implements a complex floating-point FFT for integer multiplication is George Woltman's Prime95. It is written mainly for testing large Mersenne numbers 2 p -1 for primality in the in the Great Internet Mersenne Prime Search [START_REF] Woltman | The Great Internet Mersenne Prime Search[END_REF]. It uses a DWT for multiplication mod a2 n ± c, with a and c not too large, see [START_REF] Percival | Rapid multiplication modulo the sum and difference of highly composite numbers[END_REF]. We compared multiplication modulo 2 2wn -1 in Prime95 version 24.14.2 with multiplication of n-word integers using our SSA implementation on a Pentium 4 at 3.2 GHz, and on an Opteron 250 at 2.4 GHz, see Figure 4. It is plain that Prime95 beats our implementation by a wide margin, in fact usually by more than a factor of 10 on a Pentium 4, and by a factor between 2.5 and 3 on the Opteron. Some differences between Prime95 and our implementation need to be pointed out: due to the floating point nature of Prime95's FFT, rounding errors can build up for particular input data to the point where the result are incorrectly rounded to integers. The floating point FFT can be made provably correct, see again [START_REF] Percival | Rapid multiplication modulo the sum and difference of highly composite numbers[END_REF], but at the cost of using larger FFT lengths. For example, for a length 2 25 FFT, [START_REF] Percival | Rapid multiplication modulo the sum and difference of highly composite numbers[END_REF] allows 9 bits per double, whereas Prime95 uses up to 17.76. To eliminate any chance of fatal round-offerror, the transform length and hence run-time would need to be about doubled. Also, the implementation of the FFT in Prime95 is done in hand-optimized assembly for the x86 family of processors, and will not run on other architectures.

Another implementation of complex floating point FFT is Guillermo Ballester Valor's Glucas. The algorithm it uses is similar to that in Prime95, but it is written portably in C. This makes it slower than Prime95, but still faster than our code on both the Pentium 4 and the Opteron, as shown in Figure 4.

  Bfy(a, b, t, e) 1. Write e = d*w + s with 0 <= s < w, where w is the number of bits per word 2. Decompose a = (ah, al), where ah contains the upper d words Idem for b 3. t <-(al -bl, bh -ah) 4. a <-a + b 5. t <-t * 2^s

  3 2 = 5.29 over GMP 4.1.4. In both

Figure 3 :

 3 Figure 3: Comparison of GMP 4.1.4, GMP 4.2.1, Magma V2.13-6 and our new code for the plain integer multiplication on a 2.4Ghz Opteron (horizontal axis in 64-bit words, vertical axis in seconds).

Figure 4 :

 4 Figure 4: Time in seconds for multiplication of different word lengths with our implementation, Prime95 and Glucas on a 3.2 GHz Pentium 4 and a 2.4 GHz Opteron.

  . Figure 2: Time in milliseconds needed to multiply numbers modulo 2 n +1 with an FFT of length 2 k for k = 5, 6, 7. On the right, the zoom (with only k = 5, 6) illustrates that two curves can intersect several times.To take into account those multiple crossings, the new tuning scheme determines word-intervals [m1, m2] where the FFT of length 2 k is preferred for Fermat transforms:

	0.9	mpn_mul_fft.5				0.2	mpn_mul_fft.5
	0.8	mpn_mul_fft.6 mpn_mul_fft.7				0.19	mpn_mul_fft.6
	0.7
	0.6					0.18
	0.5					0.17
	0.4					0.16
	0.3					0.15
	0.2
	0.1					0.14
	0					0.13
	500	1000	1500	2000	2500	700	750	800	850	900
	#define MUL_FFT_TABLE2 {{1, 4 /*66*/}, {401, 5 /*96*/},
	{417, 4 /*98*/}, {433, 5 /*96*/}, {865, 6 /*96*/}, ...
	The entry {433, 5 /*96*/} means that from 433 words -and up to the next size of 865 words -FFT-2 5 is pre-
	ferred, with an efficiency of 96%. A similar table is used for
	Mersenne transforms.

In the whole paper, a Fermat transform, product, or scheme is meant modulo 2 N +1, without N being necessarily a power of two as in Fermat numbers.

Bernstein attributes a similar idea to Karp in [4, §9].

Harley's trick extends to h > w: together with the segmentation method, the exact same reasoning holds.

The word "benchmarketing" has been suggested to us by Torbjörn Granlund.

Acknowledgments.

This work was done in collaboration with Torbjörn Granlund, during his visits as invited professor at INRIA Lorraine in March-April and November-December 2006; we also thank him for proof-reading this paper. This work would probably not have been achieved without the initial stimulation from Allan Steel; the authors are very grateful to him. Thanks to Markus Hegland for pointing to the Belgian paper, and to William Hart for finding a typo in a preliminary version. We also thank the anonymous referees for many valuable remarks on the manuscript.


Wim Van Ackooij
email: wim.van-ackooij@edf.fr

Jérôme Malick
email: jerome.malick@inria.fr

Second-order differentiability of probability functions

Keywords: Stochastic optimization, probabilistic constraints, joint-chance-constraints, differentiability of probabilistic function, nonlinear optimization Mathematics Subject Classification (2000) 49M37, 65K05, 90C15

In this paper, we study second-order differentiability properties of probability functions. We present conditions under which probability functions involving nonlinear systems and Gaussian (or Student) multi-variate random vectors are twice continuously differentiable. We provide an expression for their Hessian that can be useful in numerical methods for solving probabilistic constrained optimization problems.

Introduction

Given a mapping g : R n × R m → R k and a random vector ξ ∈ R m , a probability function ϕ : R n → [0, 1] is defined by ϕ(x) := P[g(x, ξ) ≤ 0], for all x ∈ R n .

(

) 1
Such functions are used to model uncertainty as probability constraints ϕ(x) ≥ p, with a user-defined safety level p ∈ (0, 1). Probability constrained optimization problems appear in fields such as energy, telecommunications, network expansion, mineral blending, chemical engineering; see, e.g., [START_REF] Arnold | A mixed-integer stochastic nonlinear optimization problem with joint probabilistic constraints[END_REF][START_REF] Henrion | Optimization of a continuous distillation process under random inflow rate[END_REF][START_REF] Morgan | Aquifer remediation design under uncertainty using a new chance constraint programming technique[END_REF][START_REF] Van Ackooij | Decomposition approaches for block-structured chance-constrained programs with application to hydro-thermal unit commitment[END_REF][START_REF] Van Ackooij | Joint chance constrained programming for hydro reservoir management[END_REF]. For a comprehensive overview on the theory and the applications of probabilistic constraints, we refer to [START_REF] Prékopa | Stochastic Programming[END_REF]18,[START_REF] Shapiro | Lectures on Stochastic Programming. Modeling and Theory, volume 9 of MPS-SIAM series on optimization[END_REF].

The numerical treatment of the probabilistic constraints by optimization algorithms opens several questions, as convexity of the level set {x ∈ R n : ϕ(x) ≥ p} and smoothness of ϕ. Both (generalized) convexity (see e.g. [START_REF] Henrion | Convexity of chance constraints with independent random variables[END_REF][START_REF] Henrion | Convexity of Chance Constraints with Dependent Random Variables: the use of Copulae[END_REF][START_REF] Prékopa | Stochastic Programming[END_REF]18,[START_REF] Van Ackooij | Eventual convexity of chance constrained feasible sets[END_REF]) and smoothness of ϕ have been the topic of several investigations. Regarding differentiability, the important results summarized in [21] prove that ϕ is continuously differentiable under the assumptions that {z ∈ R m : g(x, z) ≤ 0} is compact in a neighbourhood of x. Weaker conditions are shown to be sufficient in presence of special structure, for example separability of g(x, ξ) = ξh(x), see [START_REF] Henrion | A gradient formula for linear chance constraints under Gaussian distribution[END_REF]. For the general case of a nonlinear function g, the compactness assumption can also be replaced by growth conditions on the gradient ∇ x g, as recently developed by [START_REF] Van Ackooij | Gradient formulae for nonlinear probabilistic constraints with Gaussian and Gaussianlike distributions[END_REF] in presence of Gaussian multivariate random vectors. To our knowledge, second-order differentiability has not been investigated yetexcept for a specific case, given in the appendix of [START_REF] Van Ackooij | Joint chance constrained programming for hydro reservoir management[END_REF].

In this work, we study second-order differentiability of the probabilistic functions in the general framework of [START_REF] Van Ackooij | Gradient formulae for nonlinear probabilistic constraints with Gaussian and Gaussianlike distributions[END_REF]. We prove that ϕ is twice continuously differentiable around a trial point x, under growth conditions on ∇ 2 g(x) the Hessian of g at x. We derive a general expression of ∇ 2 ϕ(x) the Hessian of ϕ at x, as an integral with respect to the uniform distribution on the sphere. By Deàk's sampling scheme of the sphere [START_REF] Deák | Computing probabilities of rectangles in case of multinormal distribution[END_REF][START_REF] Deák | Subroutines for computing normal probabilities of sets -computer experiences[END_REF], we can then approximate ∇ 2 ϕ(x) as well as ϕ(x) and ∇ϕ(x) simultaneously. This expression of Hessians of probability functions could be useful to develop and analyze second-order methods for solving probabilistically constrained optimization problems. As shown for example in [START_REF] Van Ackooij | Constrained bundle methods for upper inexact oracles with application to joint chance constrained energy problems[END_REF], first order optimization methods dealing with probabilistic constraints spend most time evaluating the gradients of the probability functions. Since we show here that Hessians are available at the same computational cost, second order methods (see e.g. [START_REF] Izmailov | Newton-Type Methods for Optimization and Variational Problems[END_REF]) could be considered, and this could speed up the overall resolution scheme. This may be particularly so for SQP methods that have only been investigated recently in this context, in [START_REF] Bremer | Probabilistic constraints via SQP solver: Application to a renewable energy management problem[END_REF].

Here is the outline of this short paper. Section 2 presents notation, background, and our second-order differentiability result. Then section 3 develops the technical lemmas leading to the proof of our result; this part relies heavily on the analysis of [START_REF] Van Ackooij | Gradient formulae for nonlinear probabilistic constraints with Gaussian and Gaussianlike distributions[END_REF] regarding first-order differentiability. Finally section 4 presents examples and extensions.

(Second-order) differentiability under growth conditions

In this paper, we consider the general framework of [START_REF] Van Ackooij | Gradient formulae for nonlinear probabilistic constraints with Gaussian and Gaussianlike distributions[END_REF] where g is a continuously differentiable function, convex with respect to the second argument, and ξ is a multivariate non-degenerate Gaussian random vector in R m . Without loss of generality 1 , we can assume that ξ ∼ N (0, R) is centered and with an arbitrary positive definite correlation matrix R. Note that the set M (x) = {z ∈ R m : g(x, z) ≤ 0} is convex so (Lebesgue) measurable and that ξ admits a density (with respect to the Lebesgue measure). Consequently the probability function ϕ is well-defined by [START_REF] Arnold | A mixed-integer stochastic nonlinear optimization problem with joint probabilistic constraints[END_REF]. Without additional assumptions though, ϕ is not differentiable in general (see a counterexample in [24, section 2]). The approach of [START_REF] Van Ackooij | Gradient formulae for nonlinear probabilistic constraints with Gaussian and Gaussianlike distributions[END_REF] to get differentiability roughly requires to control the derivatives of g when z escapes to infinity. Specifically, the key assumption taken from [START_REF] Van Ackooij | Gradient formulae for nonlinear probabilistic constraints with Gaussian and Gaussianlike distributions[END_REF] is the following first-order exponential growth condition, that we extend here to a second-order growth condition. This framework allows for many applications; see the examples in [START_REF] Van Ackooij | Gradient formulae for nonlinear probabilistic constraints with Gaussian and Gaussianlike distributions[END_REF] and in the forthcoming section 4. Let us formalize here the two growth condition that will appear repeatedly in this paper.

Assumptions 1 (exponential growth conditions) (i) The mapping g satisfies the first order exponential growth condition at x:

∇ x g (x, z) ≤ δ exp ( z ) for all x ∈ U and all z such that z ≥ C, (2)
for a neighbourhood U of x and constants C, δ > 0. The norm • is the usual euclidian norm in R n . (ii) The mapping g satisfies the second order exponential growth condition at x if it is twice continuously differentiable at x and

∇ 2 g (x, z) ≤ δ exp ( z ) for all x ∈ U and all z such that z ≥ C, (3)
for a neighbourhood U of x and constants C, δ > 0. Here the norm of the Hessian refers to the operator norm induced by the usual euclidian norm of R n .

Throughout this paper, we consider the Cholesky decomposition of R, i.e., R = LL T . We consider now the spherical-radial decomposition of ξ (e.g., [START_REF] Genz | Computation of multivariate normal and t probabilities[END_REF][START_REF] Naor | Projecting the surface measure of the sphere of ℓ n p[END_REF])

ξ = ηLζ
with two independent random vectors: ζ having a uniform distribution over the Euclidian unit sphere S m-1 in R m and η having a chi-distribution with m degrees of freedom. We denote by µ ζ the uniform measure over the unit sphere S m-1 and by µ η the measure related to the chi-distribution with m degrees of freedom. With the factor κ := 2 1-m 2 /Γ ( m 2 ) > 0, the density associated to µ η and its first derivative are

χ (y) = κy m-1 e -y 2 /2 and χ ′ (y) = κ (m -1) -y 2 y m-2 e -y 2 /2 , (4)
and the associated cumulative function is denoted F χ . The mapping ϕ of (1) admits the following description as an integral

ϕ(x) = {r≥0,v∈S m-1 : g(x,rLv)≤0} χ(r) dr dµ ζ (v).
Given this description, the importance of studying the rays {r ≥ 0, v ∈ S m-1 : g(x, rLv) ≤ 0} becomes apparent. One can visualize a given v ∈ S m-1 as a direction along which we will investigate the length of a beam cutting the set M (x). Since the set M (x) is convex, such a beam can intersect the boundary of M (x) at most twice. Under the additional assumption that 0 ∈ int M (x), such a beam intersects the boundary at most once. We will set aside these directions, called finite or relevant directions by introducing the following set-valued mapping F : R n ⇒ S m-1 ,

F (x) = v ∈ S m-1 : ∃r > 0, g(x, rLv) = 0 . (5)
For an appropriate couple x ∈ R n and v ∈ F (x), the equation g(x, rLv) = 0 with respect to r can be solved locally around (x, v). This is exploited in the following nice alternative representation of ϕ, obtained as a consequence of Lemmas 3.1(4) and 3.3(1) of [START_REF] Van Ackooij | Gradient formulae for nonlinear probabilistic constraints with Gaussian and Gaussianlike distributions[END_REF].

Theorem 1 (Integral representation) In the above-described context, consider a point x such that g (x, 0) < 0. Then ϕ can be expressed in a neighbourhood of x as

ϕ (x) = v∈F (x) F χ (ρ x,v (x, v)) dµ ζ (v) + v ∈F (x) dµ ζ (v). ( 6
)
where F is the set-valued function defined in (5), and ρ x,v is the real-valued function defined implicitly around (x, v) by the equation g(x, rLv) = 0 with respect to r.

We notice that the condition g(x, 0) < 0 is not very restrictive as it holds in particular if ϕ(x) ≥ 1 2 under a Slater assumption on g; see [START_REF] Van Ackooij | Gradient formulae for nonlinear probabilistic constraints with Gaussian and Gaussianlike distributions[END_REF]Proposition 3.11]. From the integral expression of this theorem, we see that differentiating ϕ involves differentiating under the integral a function defined implictly. The main result of [START_REF] Van Ackooij | Gradient formulae for nonlinear probabilistic constraints with Gaussian and Gaussianlike distributions[END_REF] gives general conditions that allow this differentiation.

Theorem 2 (First-order differentiability) Let the assumption of Theorem 1 hold. Assume furthermore that g satisfies the first order exponential growth condition at x (Assumption 1(i)). Then, ϕ is continuously differentiable on U , and for all x ∈ U ∇ϕ

(x) = - v∈F (x) w=ρ x,v (x,v)Lv χ (ρ x,v (x, v)) ∇ x g (x, w) α(x, v, w) dµ ζ (v).
where α is the real-valued function defined by α(x, v, w) := 1/ ∇ z g (x, w) , Lv .

We mention that, by [START_REF] Van Ackooij | Gradient formulae for nonlinear probabilistic constraints with Gaussian and Gaussianlike distributions[END_REF]Theorem 3.14], the growth condition of this theorem can be replaced by the condition that the set {z : g(x, z) ≤ 0} is bounded (and in this case, ρ x,v is in fact independent from (x, v) and F (x) = S m-1 by [24, Lemma 3.12]). In this short paper, we establish that a similar differentiability result holds at the second-order; more precisely here is our main second-order differentiability result.

Theorem 3 (Second-order differentiability) Let the assumption of Theorem 1 hold. Assume furthermore that g satisfies the exponential growth conditions at x (Assumption 1 (i) and (ii)). Then ϕ is twice continuously differentiable on a neighbourhood U of x and for all x ∈ U

∇ 2 ϕ (x) = v∈F (x) w=ρ x,v (x,v)Lv χ ′ (ρ x,v (x, v)) ∇ x g(x, w)∇ x g(x, w) T + χ (ρ x,v (x, v)) ∇ x g(x, w)(∇ xz g(x, w)Lv) T -χ (ρ x,v (x, v)) α(x, v, w) v T L T ∇ zz g(x, w)Lv∇ x g(x, w)∇ x g(x, w) T (7) -χ (ρ x,v (x, v)) α(x, v, w) -1 ∇ xx g(x, w) + χ (ρ x,v (x, v)) v T L T ∇ zx g(x, w)∇ x g(x, w) T α(x, v, w) 2 dµ ζ (v).
Evaluating the Hessian ∇ 2 ϕ (x) with the above integral expression can be done while evaluating ϕ(x) and ∇ϕ(x) by Deàk's sampling method [START_REF] Deák | Subroutines for computing normal probabilities of sets -computer experiences[END_REF]. For each sampled point v ∈ S m-1 , we try to solve the equation g(x, rLv) = 0 explicitly or approximately (by a Newton-Raphson algorithm for instance). If there is a solution r > 0 (i.e., v ∈ F (x)) then we evaluate the integrand; else (i.e., v ∈ F (x)) the term does not contribute to the approximated integral. This procedure does not add significant computational cost on top of evaluating ϕ(x) and ∇ϕ(x). Note that this is quite different from [25, Lemma 2], where the given expression of the Hessian (of a simple probability function) has a computing cost O(m) times higher than the computing cost of the gradient.

Proof of second-order differentiability

This section provides several auxiliary developments that lead us to the proof of Theorem 3. We consider a given fixed, but arbitrary, point x ∈ R n such that g(x, 0) < 0. The bulk of the work consists in studying the second-order differentiability of the function e : R n × S m-1 → [0, 1] defined by for all x ∈ R n and v ∈ S m-1 e(x, v)

:= µ η ({r ≥ 0 : g(x, rLv) ≤ 0}) = {r≥0: g(x,rLv)≤0} χ(r) dr (8)
We start by giving more details on the implicit ray function ρ x,v and the expressions of its first and second order derivatives (when v ∈ F (x)).

Lemma 1 (Implicit ray function) If v ∈ F (x) (where F is defined as in (5)), then there exist neighbourhoods U of x and V of v as well as a twice continuously differentiable function

ρ x,v : U × V → R + satisfying the equivalence g(x, rLv) = 0 ⇐⇒ r = ρ x,v (x, v) for all (x, v, r) ∈ U × V × R + .
Moreover for all (x, v) ∈ U × V and for w = ρ x,v (x, v)Lv

∇ x ρ x,v (x, v) = -α(x, v, w)∇ x g(x, w) (9)
∇ xx ρ x,v (x, v) = α(x, v, w) 2 ∇ x g(x, w)(∇ xz g(x, w)Lv) T -α(x, v, w) 3 v T L T ∇ zz g(x, w)Lv∇ x g(x, w)∇ x g(x, w) T (10) -α(x, v, w)∇ xx g(x, w) + α(x, v, w) 2 v T L T ∇ zx g(x, w)∇ x g(x, w) T .
Proof The existence and smoothness of ρ x,v comes from applying the implicit function theorem. The details are in the proof of [START_REF] Van Ackooij | Gradient formulae for nonlinear probabilistic constraints with Gaussian and Gaussianlike distributions[END_REF]Lemma 3.2], which also provides the expression of the gradient. We focus here on the calculation of the second order derivative. To this end let us consider h(x, v) = α(x, v, w) -1 = ∇ z g(x, ρ x,v (x, v)Lv), Lv . Then for any k = 1, ..., n we have the following derivative formula:

∂h ∂x k (x, v) = ∂ ∂x k ( m j=1 ∂g ∂z j (x, ρ x,v (x, v)Lv)(Lv) j ) = m j=1 ∂ 2 g ∂x k ∂z j (x, ρ x,v (x, v)Lv)(Lv) j + m j=1 (Lv) j m ℓ=1 ∂ 2 g ∂z ℓ ∂z j (x, ρ x,v (x, v)Lv)(Lv) ℓ ∂ρ x,v ∂x k (x, v) = (∇ xz g(x, w)Lv) k -α(x, v, w)v T L T ∇ zz g(x, w)Lv ∂g ∂x k (x, w).
Consequently for any k, i = 1, ..., n we get:

∂ 2 ρ x,v ∂x i ∂x k (x, v) = h(x, v) -2 ∂h ∂x k (x, v) ∂g ∂x i (x, ρ x,v (x, v)Lv) -h(x, v) -1 ∂ 2 g ∂x i ∂x k (x, ρ x,v (x, v)Lv) + h(x, v) -1 m j=1 ∂ 2 g ∂z j ∂x i (x, ρ x,v (x, v)Lv)(Lv) j ∂ρ x,v ∂x k (x, v) ,
which by substitution gives

∂ 2 ρ x,v ∂x i ∂x k (x, v) = h(x, v) -2 ∂g ∂x i (x, ρ x,v (x, v)Lv)(∇ xz g(x, ρ x,v (x, v)Lv)Lv) k -h(x, v) -3 v T L T ∇ zz g(x, ρ x,v (x, v)Lv)Lv ∂g ∂x i (x, ρ x,v (x, v)Lv) ∂g ∂x k (x, ρ x,v (x, v)Lv) -h(x, v) -1 ∂ 2 g ∂x i ∂x k (x, ρ x,v (x, v)Lv) + h(x, v) -2 (v T L T ∇ zx g(x, ρ x,v (x, v)Lv)) i ∂g ∂x k (x, ρ x,v (x, v)Lv).
This can be written as [START_REF] Henrion | A gradient formula for linear chance constraints under Gaussian distribution[END_REF].

⊓ ⊔

The second-order differentiability of e when v ∈ F (x) follows immediately from the previous lemma together with Corollary 3.5 of [START_REF] Van Ackooij | Gradient formulae for nonlinear probabilistic constraints with Gaussian and Gaussianlike distributions[END_REF].

Lemma 2 (In relevant directions) If v ∈ F (x), then there exist neighbourhoods U of x and V of v such that, for (x, v) ∈ U × V , the function e defined in (8) satisfies e(x, v) = F χ (ρ x,v (x, v)) as well as

∂e ∂x k (x, v) = -α(x, v, w)χ(ρ x,v (x, v)) ∂g ∂x k (x, w)
and is moreover twice continuously differentiable with partial derivative:

∂ 2 e ∂x k ∂x i (x, v) = χ ′ (ρ x,v (x, v))α(x, v, w) 2 ∂g ∂x k (x, w) ∂g ∂x i (x, w) + χ ρ x,v (x, v) ∂ 2 ρ x,v ∂x k ∂x i (x, v),
where k, i ∈ {1, ..., n} are arbitrary and w = ρ x,v (x, v)Lv.

We now show that Assumptions 1 yield that the exponentials of the χ-distribution make the Hessian of e vanish when getting close to irrelevant directions (v ∈ F (x)).

Lemma 3 (Toward irrelevant directions) Let v ∈ F (x) and consider a sequence

(x k , v k ) → (x, v) with v k ∈ F (x k ). If g satisfies Assumption 1 (i) and (ii) at x, then lim k→∞ ∇ xx e(x k , v k ) = 0.
Proof By [24, Lemma 3.3], we have first that

ρ k := ρ x k ,v k (x k , v k ) → ∞ when k → ∞, since x k → x and v k ∈ F (x k ).
Recall also that, by [24, Lemma 3.1(3)], there exists δ 1 > 0 such that, on an appropriately small neighbourhood of x,

∇ z g (x k , ρ k Lv k ) , Lv k ≥ δ 1 [ρ x k ,v k (x k , v k )] -1 = δ 1 ρ -1 k > 0.
Now note that, by Lemma 2, we have the upper bound:

∇ xx e(x k , v k ) ≤ χ ′ (ρ k ) ∇ z g(x k , ρ k Lv k ), Lv k 2 ∇ x g(x k , ρ k Lv k ) 2 + χ(ρ k ) ∇ xx ρ x k ,v k (x k , v k ) . ( 11
)
By using the growth condition for k sufficiently large and the expressions of (4), we can bound the first term (which we call A k for short) by

A k ≤ δ 2 δ -2 1 ρ 2 k κ (m -1) -ρ 2 k [ρ k ] m-2 e -[ρ k ] 2 /2 e 2 L ρ k = δκδ -2 1 κ (m -1) -ρ 2 k [ρ k ] m e 2 L ρ k e -[ρ k ] 2 /2 ,
which allows us to show that lim k→∞ A k = 0 since ρ k := ρ x k ,v k (x k , v k ) → ∞ and y α e y e -y 2 /2 → 0 for y → ∞, where α > 0 is arbitrary. We bound the second term in (11) by using the second order growth condition:

χ(ρ x k ,v k (x k , v k )) ∇ xx ρ x k ,v k (x k , v k ) ≤ κδδ -1 1 δ -1 1 ρ 2 k L e 2 L ρ k [ρ k ] m-1 e -[ρ k ] 2 /2 + δ -2 1 ρ 3 k L 2 e 3 L ρ k [ρ k ] m-1 e -[ρ k ] 2 /2 +ρ k e L ρ k [ρ k ] m-1 e -[ρ k ] 2 /2 + δ -1 1 ρ 2 k L e 2 L ρ k [ρ k ] m-1 e -[ρ k ] 2 /2 .
The right-hand side can be seen to converge to zero following the above evoked arguments.

⊓ ⊔

We are now in position to establish the twice differentiability of e, which is the main technical lemma in our way to Theorem 3.

Lemma 4 (Second derivative of the function e) If g satisfies the first and second order exponential growth conditions at x, then there is a neighbourhood U of x such that the mapping e is twice continuously differentiable with respect to x on U . For x ∈ U and v ∈ F (x), the gradient and Hessian of e are null, and for v ∈ F (x), the gradient is given in Lemma 2 and the Hessian by

∇ xx e(x, v) = α(x, v, w) 2 χ ′ (ρ x,v (x, v)) ∇ x g(x, w)∇ x g(x, w) T + χ (ρ x,v (x, v)) ∇ x g(x, w)(∇ xz g(x, w)Lv) T -χ (ρ x,v (x, v)) α(x, v, w) v T L T ∇ zz g(x, w)Lv∇ x g(x, w)∇ x g(x, w) T -χ (ρ x,v (x, v)) /α(x, v, w) ∇ xx g(x, w) + χ (ρ x,v (x, v)) v T L T ∇ zx g(x, w)∇ x g(x, w) T ,
where w = ρ x,v (x, v)Lv. We have moreover that ∇ xx e is continuous with respect to the couple (x, v).

Proof Fix i, ℓ ∈ {1, ..., n}. For any v ∈ F (x), the provided formula and differentiability statement follows from Lemma 2. We focus here on the case v ∈ F (x). We will show, by contradiction, that,

lim t↑0 ∂e ∂x ℓ (x + tu i , v) -∂e ∂x ℓ (x, v) t = 0, (12)
where u i is the i-th canonical unit vector in R n . In exactly the same way, one can also show that the limit for t ↓ 0 equals zero too. Altogether, this will prove that e is twice differentiable at x and that

∂ 2 e ∂x i ∂x ℓ e(x, v) = 0.
Notice first that, by [24, Corollary 3.8], we have ∂e ∂x ℓ (x, v) = 0. Negating [START_REF] Henrion | Convexity of Chance Constraints with Dependent Random Variables: the use of Copulae[END_REF] implies the existence of some ε > 0 and a subsequence t k ↑ 0 such that

∂e ∂x ℓ (x + t k u i , v) t k ≥ ε. (13)
This implies in particular that ∂e ∂x ℓ (x + t k u i , v) = 0 and then v ∈ F (x + t k u i ), again by [24, Corollary 3.8]. We also have that g(x + t k u i , 0) < 0 for all large k. Now, for an arbitrary k, we define (recall that t k < 0)

T := inf τ ∈ [t k , 0] : ∂e ∂x ℓ (x + τ u i , v) = 0 (≤ 0).
Due to ∂e ∂x ℓ (x, v) = 0 we have that T ≤ 0. On the other hand, ∂e ∂x ℓ (x + t k u i , v) = 0 and the continuity of ∇ x e established in [24, Corollary 3.9] provides T > t k . We infer that ∂e ∂x ℓ (x+τ u i , v) = 0 for all τ ∈ [t k , T ) and, hence,

v ∈ F (x + τ u i ) ∀τ ∈ [t k , T ). ( 14
)
But then, the function

β(τ ) := ∂e ∂x ℓ (x + τ u i , v)
is differentiable for all τ ∈ (t k , T ) by virtue of Lemma 2 and its derivative is given by

β ′ (τ ) = n j=1 ∂ 2 e ∂x j ∂x ℓ (x + τ u i , v)u i (j) = ∂ 2 e ∂x i ∂x ℓ (x + τ u i , v).
Therefore, the mean value theorem guarantees the existence of some τ * k ∈ (t k , T ) such that

β ′ (τ * k ) = β(T ) -β(t k ) T -t k .
By continuity of ∇ x e of [24, Corollary 3.9], we have that β(T ) = ∂e ∂x ℓ (x + T u i , v) = 0, whence, by

t k < T ≤ 0, |β ′ (τ * k )| = -∂e ∂x ℓ (x + t k u i , v) T -t k ≥ ∂e ∂x ℓ (x + t k u i , v) t k ≥ ε,
where the last inequality is [START_REF] Hunter | Applied Analysis[END_REF]. Now, since k was arbitrarily fixed, we have constructed a sequence

τ * k such that t k < τ * k ≤ 0 such that ∂ 2 e ∂x i ∂x ℓ (x + τ u i , v) ≥ ε ∀k. ( 15
)
Since t k ↑ 0, we also have that [START_REF] Izmailov | Newton-Type Methods for Optimization and Variational Problems[END_REF]. Due to our assumption that g satisfies the growth conditions at x, Lemma 3 yields that lim k→∞ ∇ xx e(x k , v) = 0 which contradicts [START_REF] Morgan | Aquifer remediation design under uncertainty using a new chance constraint programming technique[END_REF] and thus concludes about twice differentiability. Finally, the fact that the second order partial derivative ∇ xx e is also continuous at (x, v) for any v ∈ S m-1 can be shown as follows. For a given x and v ∈ F (x), the provided formula of the Hessian holds true locally and continuity is evident. In contrast, if v ∈ F (x), assuming that ∇ xx e(x k , v k ) does not tend to zero along a sequence (x k , v k ) → (x, v), means in particular that ∇ xx e(x k , v k ) ≥ ε for a given ε > 0. Consequently by the arguments above, v k ∈ F (x k ). But then this leads to a contradiction with Lemma 3.

τ * k ↑ 0. Moreover, v ∈ F (x + τ * k u i ) by

⊓ ⊔

The proof of Theorem 3 about second-order differentiability of ϕ now follows from the previous result after justifying the derivation under the integral.

Proof (of Theorem 3) Theorem 1 gives the integral representation (6) of ϕ around x. Theorem 2 allows us to obtain an integral representation for ∇ϕ(x) involving e. Indeed for a given k = 1, ..., n, we have

∂ϕ ∂x k (x) = v∈S m-1 ∂e ∂x k (x, v)dµ ζ (v), (16)
as a consequence of Theorem 2 and expressions of derivatives of e of Lemmas 2 and 4. Invoking again Lemma 4, the second order partial derivative ∇ xx e of the function e exists and is continuous on U ×S m-1 . The continuity of ∇ xx e on U × S m-1 and compactness of S m-1 guarantee that the function γ :

U → R γ(x) := max v∈S m-1 ∇ xx e(x, v)
is well-defined and continuous. Restricting U if necessary, we may assume that γ(x) ≤ 2γ(x) for all x ∈ U . This together with µ ζ (S m-1 ) = 1 for the law µ ζ of the uniform distribution on S m-1 we infer that the constant 2γ(x) is an integrable function on S m-1 dominating ∇ xx e(x, v) on S m-1 for all x ∈ U . Now, Lebesgue's dominated convergence theorem [13, Theorem 12.35, Corollary 12.36] allows to differentiate ( 16) under the integral sign and provides continuity with respect to x. We have established

∇ 2 ϕ (x) = v∈S m-1 ∇ xx e(x, v)dµ ζ
which gives the desired expression by Lemma 4.

⊓ ⊔ 4 Examples and extensions

We illustrate now our second-order differentiability result on some special examples: for chi-squared random vectors in section 4.1 and for lognormal random vectors in section 4.2. Finally we show in section 4.3 that results similar to Theorem 3 also hold with Student random vectors.

Chi-squared example

We consider here a linear probabilistic function

ϕ(x) = P ( η, x ≤ b) (17)
with b > 0 and a random vector η ∈ R n whose components η i are independent and have a χ 2 -distribution with m i degrees of freedom. For each i = 1, . . . , n, we can write

η i = m i k=1 ξ 2 i,k
, where m i is the degree of freedom, and ξ i,k ∼ N (0, 1) (for k = 1, . . . , m i ) are independent random variables. We consider the Gaussian random vector of

m = m 1 + • • • + m n dimensions ξ := (ξ 1,1 , . . . , ξ 1,m 1 , . . . , ξ n,1 , . . . , ξ n,m n ) ∼ N (0, I) .
We also introduce the function f : R m → R n defined by its components

f i (z) = m i k=1 z 2
i,k (for z partitioned as ξ above), and the function g : R n × R m → R defined by g(x, z) = x, f (z)b. Thus, we are in the situation of this paper as the probability function [START_REF] Prékopa | Stochastic Programming[END_REF] writes

ϕ(x) = P ( η, x ≤ b) = P ( f (ξ), x ≤ b) = P (g(x, ξ) ≤ 0) .
We observe that g is C 2 , satisfies the first and second order growth conditions and is convex with respect to z as soon as x i ≥ 0 for all i = 1, ..., n. We apply Theorem 3 at x (such that xi > 0 for all i = 1, ..., n) to establish that ϕ of ( 17) is C 2 around x, and we can specify the integral expression of the Hessian for this case as follows. We observe that the set {z ∈ R m : g(x, z) ≤ 0} is bounded, so that F (x) = S m-1 by [START_REF] Van Ackooij | Gradient formulae for nonlinear probabilistic constraints with Gaussian and Gaussianlike distributions[END_REF]Lemma 3.12]. We also dispose of an explicit ray function (for all (x, v)) as the equation f (rv), x = b in r admits the solution

ρ(x, v) = b/ f (v), x .
We readily compute the terms that appeared in the expression of Theorem 3, including

α(x, v, ρ(x, v)v) = 1/ 2 b f (v), x and v T ∇ zz g(x, ρ (x, v) v)v = 2 f (v), x .
Substituting these expressions in [START_REF]Encyclopedia of Optimization[END_REF] and simplifying yields:

∇ 2 ϕ(x) = v∈S m-1 b χ ′ b/ f (v), x + 3 b f (v), x χ b/ f (v), x f (v) T f (v) 4 f (v), x 3 dµ ζ (v).

Multivariate log-normal distribution

We consider here the probabilistic function

ϕ(x) = P ( η, x ≤ h(x)) (18)
with a twice continuously differentiable h : R m → R and a random vector η with a multivariate lognormal distribution, i.e., the component-wise logarithm log η has a multivariate Gaussian distribution). Without loss of generality we assume that ξ = log η ∼ N (0, R) for some correlation matrix R = LL T . Thus we remark that the probabilistic function of (18) can be written as ϕ(x) = P[g(x, ξ) ≤ 0] with the function g(x, z) = x, e zh(x), which is C 2 , satisfies the first and second order growth conditions, and is convex with respect to z as soon as x i ≥ 0 for all i = 1, ..., m.

Let x be such that xi > 0 for i = 1, . . . , m and h(x) > m i=1 xi . This implies that g(x, 0) < 0 and we can therefore apply Theorem 3, to get that ϕ of (18) is C 2 around x (more precisely, in a neighbourhood U of x such that g (x, 0) < 0 and x i > 0 for i = 1, . . . , m). The expression of α and F instantiates as

α(x, v, w) = m i=1 x i e ρ x,v (x,v)(Lv) i (Lv) i -1 and F (x) = {v ∈ S m-1 : ∃ i such that (Lv) i > 0}. (19)
The integral expression of the Hessian then simplifies to

∇ 2 ϕ (x) = v∈F (x)
χ ′ (ρ x,v (x, v)) e ρ x,v (x,v)Lv -∇h(x) e ρ x,v (x,v)Lv -∇h(x)

T + χ (ρ x,v (x, v)) e ρ x,v (x,v)Lv -∇h(x) diag(e ρ x,v (x,v)Lv )Lv T - χ (ρ x,v (x, v)) α x, v, ρ x,v (x, v) v T L T diag(xe ρ x,v (x,v)Lv )Lv e ρ x,v (x,v)Lv -∇h(x) e ρ x,v (x,v)Lv -∇h(x) T + χ (ρ x,v (x, v)) α x, v, ρ x,v (x, v) ∇ 2 h(x) + χ (ρ x,v (x, v)) v T L T diag(e ρ x,v (x,v)Lv ) e ρ x,v (x,v)Lv -∇h(x) T α(x, v, w) 2 dµ ζ (v).
Here, e z for a vector z has to be understood component-wise, diag(y) denotes the diagonal matrix with vector y on the diagonal, and xe z has to be understood as the components wise product. The expression of F (x) in ( 19) is the main technical point to applying Theorem 3 to the example of this section. For sake of completeness, let us provide a quick proof of this expression: we show by double inclusion that, for all

x ∈ U , S m-1 \ F (x) = v ∈ S m-1 : Lv ≤ 0 . (20)
First, let x ∈ U and v ∈ S m-1 with Lv ≤ 0; then for all r > 0 g(x, rLv) = e rLv , xh(x) ≤ e 0 , xh(x) = g (x, 0) < 0, whence v / ∈ F (x). Conversely, let x ∈ U and v / ∈ F (x). Then, e rLv , x < h(x) for all r > 0. Define J := {i : (Lv) i > 0} and observe that i∈J x i e r(Lv) i is bounded from above independently of r. If J = ∅, then this sum would tend to ∞ for r → ∞ which is a contradiction. Consequently, J = ∅, proving Lv ≤ 0 and, thus, the reverse inclusion of (20).

Second-order differentiability with Student's law

The rationale of the proof of section 3 can be adapted to a similar context under slight changes in conditions and analysis. In this section, we show for example that the second-order differentiability still holds when the random vector ξ in (1) follows a multi-variate Student distribution, under a slightly adapted growth condition on the Hessian of g. Consider the probability function ϕ defined in [START_REF] Arnold | A mixed-integer stochastic nonlinear optimization problem with joint probabilistic constraints[END_REF] where ξ ∼ T (0, R, ν) has a standard multi-variate Student distribution (with correlation matrix R = LL T and ν degrees of freedom). By [START_REF] Van Ackooij | Gradient formulae for nonlinear probabilistic constraints with Gaussian and Gaussianlike distributions[END_REF]Theorem 4.4], ϕ can be expressed as

ϕ(x) = v∈S m-1 ẽ(x, v)dµ ζ (v) with ẽ(x, v) = F m,ν (m -1 [ρ x,v (x, v)] 2 ) if v ∈ F (x) 1 else , (21)
where F m,ν (t) is the Fisher-Snedecor distribution function with density:

f m,ν (t) = Γ (m/2+ν/2) Γ (m/2)Γ (ν/2) m m/2 ν ν/2 t m/2-1 (mt + ν) -(m+ν)/2 t ≥ 0 0 t < 0 ( 22
)
Theorem 4 Let the mapping g : R n × R m → R be twice continuously differentiable and convex in the second argument. Assume moreover that g satisfies a first and second order polynomial growth condition at x (i.e., the exponential term in (3) is replaced with z ℓ , and the one of (2) by z κ and these coefficients are related as follows ℓ + 2κ < ν -2). Then ϕ of (21) is twice continuously differentiable on a neighbourhood U of x and for all x ∈ U ∇ 2 ϕ (x) = v∈F (x) w=ρ x,v (x,v)Lv

4 ρ x,v (x, v) 2 f ′ m,ν (m -1 ρ x,v (x, v) 2 ) mf m,ν (m -1 ρ x,v (x, v)
2 ) ∇ x g(x, w)∇ x g(x, w) T + 2∇ x g(x, w)∇ x g(x, w) T + 2ρ x,v (x, v)∇ x g(x, w)(∇ xz g(x, w)Lv) T -2ρ x,v (x, v)α(x, v, w) v T L T ∇ zz g(x, w)Lv∇ x g(x, w)∇ x g(x, w) T -2ρ x,v (x, v)/α(x, v, w)∇ xx g(x, w) + 2ρ x,v (x, v)v T L T ∇ zx g(x, w)∇ x g(x, w) T f m,ν ρ x,v (x, v) 2 m α(x, v, w)

2 m dµ ζ (v)
Proof The continuous differentiability is established in [START_REF] Van Ackooij | Gradient formulae for nonlinear probabilistic constraints with Gaussian and Gaussianlike distributions[END_REF]Theorem 4.7]. We focus on the proof of the second-order differentiability and, as remarked in [START_REF] Van Ackooij | Gradient formulae for nonlinear probabilistic constraints with Gaussian and Gaussianlike distributions[END_REF], we only need to derive the equivalent of Lemma 3 for the Student setting. Following the arguments of Lemma 3, we have

∇ xx ẽ(x k , v k ) ≤ 4ρ x,v (x, v) 2 f ′ m,ν (m -1 [ρ x,v (x, v)] 2 ) m 2 ∇ z g(x k , ρ x k ,v k (x k , v k )Lv k ), Lv k 2 ∇ x g(x k , ρ x k ,v k (x k , v k )Lv k ) 2 + 2f m,ν (m -1 [ρ x,v (x, v)] 2 ) m ∇ z g(x k , ρ x k ,v k (x k , v k )Lv k ), Lv k 2 ∇ x g(x k , ρ x k ,v k (x k , v k )Lv k ) 2 + 2 m ρ x,v (x, v) f m,ν (m -1 [ρ x,v (x, v)] 2 ) ∇ xx ρ x k ,v k (x k , v k ) .
Using the expression of f ′ m,ν (t) (for t ≥ 0)

f ′ m,ν (t) =
Γ (m/2 + ν/2) Γ (m/2)Γ (ν/2) m m/2 ν ν/2 t m/2-1 (mt + ν)

-(m+ν)/2 (( m 2 -1)t -1 -m m + ν 2 (mt + ν) -1 )
and the polynomial growth condition, we can bound the first term (A k ) as:

A k ≤ C1 ρ m+2κ k (ρ 2 k + ν) -m+ν 2 + C2 ρ m+2κ+2 k (ρ 2 k + ν) -m+ν+2 2 ,
for appropriate constants C1 , C2 the last term tends to zero whenever κ < ν 2 and where ρ k is short for ρ k := ρ x k ,v k (x k , v k ). The second term (B k ) can be bounded by:

B k ≤ C3 ρ x k ,v k (x k , v k ) m+2κ (ρ x k ,v k (x k , v k ) 2 + ν) -m+ν 2
which again tends to zero when k → ∞ as long as κ < ν 2 . Finally the last term (C k ) can be bounded by using the second order polynomial growth condition as follows:

C k ≤ C4 ρ m k (ρ 2 k + ν) -m+ν 2 (2ρ 1+κ+ℓ k + ρ 2+2κ+ℓ k + ρ ℓ k
). When ℓ + 2κ < ν -2, we also have ℓ < ν, 1 + κ + ℓ < ν and κ < ν 2 , so that this condition entails that also C k → 0, when k → ∞. We have thus established the equivalent of Lemma 3: the mapping ẽ has the same properties as e. Consequently the remaining proofs carry forth.

⊓ ⊔

[2] [19] [START_REF]Encyclopedia of Optimization[END_REF]