Vincent Lefèvre

Paul Zimmermann

Optimized Binary64 and Binary128 Arithmetic with GNU MPFR

Keywords: floating point, correct rounding, faithful rounding, binary64, binary128, GNU MPFR

We describe algorithms used to optimize the GNU MPFR library when the operands fit into one or two words. On modern processors, a correctly rounded addition of two quadruple precision numbers is now performed in 22 cycles, a subtraction in 24 cycles, a multiplication in 32 cycles, a division in 64 cycles, and a square root in 69 cycles. We also introduce a new faithful rounding mode, which enables even faster computations. Those optimizations will be available in version 4 of MPFR.

I. INTRODUCTION

The IEEE 754-2008 standard [START_REF]Revision of ANSI-IEEE Standard[END_REF] defines -among others -the binary floating-point formats binary64 and binary128. They provide respectively a precision of 53 and 113 bits. Those standard formats are used in many applications, therefore it is very important to provide fast arithmetic for them, either in hardware or in software.

GNU MPFR [START_REF] Fousse | MPFR: A multiple-precision binary floating-point library with correct rounding[END_REF] (MPFR for short) is a reference software implementation of the IEEE 754-2008 standard in arbitrary precision. MPFR guarantees correct rounding for all its operations, including elementary and special functions [6, Table 9.1]. It provides mixed-precision operations, i.e., the precision of the input operand(s) and the result may differ.

Since 2000, several authors have cited MPFR, either to use it for a particular application, or to compare their own library to MPFR [START_REF] De Dinechin | Towards the postultimate libm[END_REF], [START_REF] Graillat | Accurate summation, dot product and polynomial evaluation in complex floating point arithmetic[END_REF], [START_REF] Johansson | Efficient implementation of elementary functions in the medium-precision range[END_REF], [START_REF] Johansson | Efficient arbitrary-precision midpoint-radius interval arithmetic[END_REF]. Most of those comparisons are in the case when all operands have the same precision, which is usually one of the standard binary64 or binary128 formats.

Since its 4.3.0 release, GCC uses MPFR in its middleend to generate correctly rounded compile-time results regardless of the math library implementation or floating-point precision of the host platform. Several interval arithmetic libraries depend on MPFR: MPFI [START_REF] Revol | Motivations for an arbitrary precision interval arithmetic and the MPFI library[END_REF], GNU Octave interval package 1 , libieeep1788 2 , Moore [START_REF] Mascarenhas | Interval arithmetic in modern C++[END_REF]. The Sage computer algebra system 3 uses MPFR with a default precision of 53 bits for its numerical computations.

The contributions of this article are: (i) a description of new algorithms for basic floating-point arithmetic in one-or two-word precision, (ii) an implementation in MPFR of those algorithms with corresponding timings for basic arithmetic and mathematical functions, (iii) a description of a new faithful rounding mode, with corresponding implementation and timings in MPFR.

Notations: We use p to denote the precision in bits. A limb, following GMP terminology [START_REF] Granlund | GNU MP: The GNU Multiple Precision Arithmetic Library[END_REF], is an unsigned integer that fits in a machine word, and is assumed to have 64 bits here. Some algorithms are also valid when the radix is not 2 64 , we then denote β the radix, assumed to be a power of two. We call left shift (resp. right shift) a shift towards the most (resp. least) significant bits. We use a few shorthands: HIGHBIT denotes the limb 2 63 ; ONE denotes the limb 1; GMP_NUMB_BITS, which is the number of bits in a limb, is replaced by 64. In the source code, a 0 and e min are written a0 and emin respectively.

II. BASIC ARITHMETIC A. The MPFR Internal Format

Internally, a MPFR number is represented by a precision p, a sign s, an exponent e, denoted by EXP(•), and a significand m. The significand is a multiple-precision natural integer represented by an array of n = p/64 limbs: we denote by m[0] the least significant limb, and by m[n -1] the most significant one. For a regular number -neither NaN, nor ±∞, nor ±0 -the most significant bit of m[n -1] must be set, i.e., the number is always normalized 4 . The corresponding number is:

(-1) s • (m • 2 -64n ) • 2 e ,
i.e., the rational significand m/2 64n lies in [1/2, 1[. When the precision p is not an exact multiple of 64, the least sh = 64n -p bits of m[0] are not used. By convention they must always be zero, like the 3 bits below (drawing the most significant limbs and bits on the left):

1xxxxxxxxxxx m[n-1] xxxxxxxxxxxx m[n-2] • • • xxxxxxxxx000 m[0]
In the description of the algorithms below, those sh trailing bits will be represented by 3 zeros.

In this section, we describe basic algorithms used to perform arithmetic on 1 or 2 limbs (addition and subtraction, multiplication, division and square root) for regular inputs, after non-regular inputs have been dealt with. Those algorithms assume the number sh of unused bits is non-zero. As a consequence, on a 64-bit processor, if efficiency is required, it is recommended to use precision p = 63 or p = 127 instead of p = 64 or p = 128, whenever possible. In the following, we assume that all input and output arguments have the same precision p. Note that MPFR functions have to take into account that input and output variables can be the same, for example mpfr_sub (a, b, a, rnd_mode) will round b -a and store the result in a.

In addition to correct rounding, all MPFR functions return a ternary value, giving the sign of the rounding error: 0 indicates that the computed correctly rounded result y exactly equals the infinite precision value f (x) (no rounding error did occur); a positive value indicates that y > f (x), and a negative value indicates that y < f (x). Determining that ternary value is sometimes more expensive than determining the correct rounding: for example if two high-precision numbers x and 1 -x are added with a low target precision, the rounding error will usually be less than 1 2 ulp, thus we can easily decide that the correctly rounded result to nearest is 1, but more work is needed to determine the ternary value.

B. Addition

We describe here the internal mpfr_add1sp1 function for 0 < p < 64. Let b and c be the input operands, and a the result. By "addition", we mean that b and c have same sign. Their significands consist of one limb only: bp[0] and cp[0], and their exponents are bx and cx. Since b and c are regular numbers, we have 2 63 ≤ bp[0], cp[0] < 2 64 . If bx = cx, the addition of the significands always produces a carry:

2 64 ≤ bp[0] + cp[0] < 2 65 : bp[0] 1xxxxxxxx000 1yyyyyyyy000 cp[0]
We thus simply add the two shifted significands -since p < 64, we lose no bits in doing this -and increase bx by 1 (this will be the result exponent):

a0 = (bp[0] >> 1) + (cp[0] >> 1); bx ++;
Since sh = 64 -p is the number of unused bits of a 0 , the round bit is bit sh -1 from a 0 , which we set to 0 before storing in memory, and the sticky bit is always 0 in this case: ). An overflow occurs when bx > emax. Otherwise we set the exponent of a to bx. If rb = sb = 0, we return 0 as ternary value, which means the result is exact, whatever the rounding mode. If rounding to nearest, we let ap[0] unchanged when either rb = 0, or rb = 0 and sb = 0 and bit sh of ap[0] is 0 (even rule); otherwise we add one ulp, i.e., 2 sh to ap[0]. If rounding toward zero, ap[0] is unchanged, and we return the opposite of the sign of a as ternary value, which means that the computed value is less than b+c when b, c > 0 (remember the exact case rb = sb = 0 was treated before). If rounding away from zero, we add one ulp to ap[0]; while doing this, we might have a carry, which might in turn give an overflow.

C. Subtraction

We detail here the mpfr_sub1sp2 function, in the case where the exponent difference d = bx -cx satisfies 0 < d < 64. In that case, the significand of c overlaps with the upper word of the significand of b. We align c on b, giving 3 limbs cp [START_REF] Brent | Modern Computer Arithmetic[END_REF] >> d, t and sb:

bp[1] 1xxxxxxxxxxx bp[0] xxxxxxxxx000 ← d bits → 1yyyyyyyyyyy cp[1] yyyyyyyyy000 cp[0] yyyyyyyyyyyy cp[1] >> d yyyyyyyyyyyy t yyyyyy000000 sb
The goal is to subtract these 3 limbs from b, word by word. While this 3-word subtraction could be done with only 3 instructions on a typical processor, we need to write a portable C code, which is more complex to express the borrow propagation in C (and unfortunately, compilers are not yet able to detect patterns to generate only these 3 instructions):

t = (cp[1] << (64 -d)) | (cp[0] >> d); a0 = bp[0] -t; a1 = bp[1] -(cp[1] >> d) -(bp[0] < t); sb = cp[0] << (64 -d); if (sb) {
a1 -= (a0 == 0); a0 --; sb = -sb; / * 2^64 -sb * / } At this point the exact subtraction corresponds to a 1 + 2 -64 a 0 + 2 -128 sb, where a 1 and a 0 cannot be both zero: if d ≥ 2 then a 1 ≥ 2 62 ; if d = 1, since bit 0 of cp[0] is 0 because p < 128, necessarily sb = 0. However we can have a 1 = 0, in which case we shift a by one word to the left, and decrease the exponent by 64 (this can only occur for d = 1, in which case sb = 0): We now compute the round and sticky bit, and set to zero the last sh bits of a 0 before storing it in memory, where mask = 2 sh -1: The rounding is done exactly like for the addition, except in the subtraction we cannot have any overflow since here all arguments have the same precision, but we can have an underflow.

if (a1 == 0) { a1 = a0; a0 = 0; bx -= 64; } Now a 1 = 0,

D. Multiplication

The multiplication is easy to implement with the internal format chosen for MPFR. We detail here the mpfr_mul_1 function, for 0 < p < 64. We first add the two exponents and multiply the two input significands bp[0] and cp[0] using GMP's umul_ppmm macro, which multiplies two 64-bit words, and stores the upper and lower words in a 0 and sb here:

ax = EXP(b) + EXP(c); umul_ppmm (a0, sb, bp[0], cp[0]); Since 2 63 ≤ bp[0], cp[0] < 2 64 , we have 2 126 ≤ bp[0]cp[0] = a 0 • 2 64 + sb < 2 128
. The upper word a 0 therefore satisfies 2 62 ≤ a 0 < 2 64 , and in case 2 62 ≤ a 0 < 2 63 , we have to shift a 0 and sb by one bit to the left and decrease the output exponent:

if (a0 < HIGHBIT) { ax --; a0 = (a0 << 1) | (sb >> 63); sb = sb << 1; }
The round and sticky bits are computed exactly like in the 2-limb subtraction, and the sign of the result is set:

rb = a0 & (ONE << (sh -1)); sb |= (a0 & mask) ^rb; ap[0] = a0 & ~mask; SIGN(a) = MULT_SIGN (SIGN (b), SIGN (c));
For the multiplication, both overflow and underflow can happen. Overflow is easy to handle. For underflow, MPFR signals underflow after rounding: this means that underflow occurs when the rounded result -with an unbounded exponent range -has an exponent smaller that e min . Thus for example, when ax = e min -1, rb or sb is non-zero, ap[0] = 111...111000 and we round away from zero, there is no underflow. Apart from this special case, the rounding is the same as for the subtraction. (This special case where ax = e min -1 and nevertheless no underflow occurs cannot happen in the subtraction, since b and c are multiples of 2 emin-p , likewise for b -c; thus if the exponent ax of the difference is smaller than e min , the difference b -c is exact, and rb = sb = 0.)

In the mpfr_mul_2 function, which multiplies b 1 •2 64 +b 0 by c 1 • 2 64 + c 0 , the product b 0 c 0 of the low limbs contributes to less than 1 to the second most significant limb of the full 4-limb product (before a potential left shift by 1 bit to get a normalized result). Thus we only perform 3 calls to umul_ppmm, to compute b 1 c 1 , b 1 c 0 and b 0 c 1 . Moreover, we ignore the lower limbs of b 1 c 0 and b 0 c 1 , which yields a 127or 128-bit approximate product with error less than 3.

E. Division

The division code for 1 or 2 limbs first computes an approximation of the quotient, and if this approximation is not sufficient to determine the correct rounding, an exact quotient and remainder are computed, starting from that approximation. The algorithms use GMP's invert_limb function (cur-

rently re-implemented in MPFR since it is not yet in GMP's public interface). Given a limb v such that β/2 ≤ v < β, invert_limb returns the limb (β 2 -1)/v -β, which is called the reciprocal of v.
All division functions first check whether the dividend significand u is larger or equal to the divisor significand v, where β n /2 ≤ u, v < β n . If that occurs, v is subtracted from u, which generates an extra leading quotient bit. After that subtraction, we get u -v < v since u < β n ≤ 2v. Therefore in the following we assume u < v.

We detail in Algorithm DivApprox1 how we compute an approximate quotient q for the mpfr_div_1 function, which divides two one-limb significands. For β = 2 64 , the approx-

Algorithm 1 DivApprox1 Input: integers u, v with 0 ≤ u < v and β/2 ≤ v < β Output: integer q approximating uβ/v 1: compute an approximate reciprocal i of v, satisfying i ≤ (β 2 -1)/v -β ≤ i + 1 (1)
2: q = iu/β + u imate reciprocal i in step 1 is obtained using the variable v 3 from Algorithm 2 (RECIPROCAL_WORD) in [START_REF] Möller | Improved division by invariant integers[END_REF]; indeed Theorem 1 from [START_REF] Möller | Improved division by invariant integers[END_REF] proves that -with our notation -

0 < β 2 -(β + v 3 )v < 2v, from which the inequalities (1)
follow.

Theorem 1: The approximate quotient returned by Algorithm DivApprox1 satisfies

q ≤ uβ v ≤ q + 2.
Proof.

Step 2 of DivApprox1 is simply step 1 of Algorithm 1 (DIV2BY1) from [START_REF] Möller | Improved division by invariant integers[END_REF]. If i is the exact reciprocal i 0 := (β 2 -1)/v -β, and q 0 := i 0 u/β +u is the corresponding quotient, it is proven in [START_REF] Möller | Improved division by invariant integers[END_REF] that the corresponding remainder r 0 = βu -q 0 v satisfies 0 ≤ r 0 < 4v. However the upper bound 4v includes 2v coming from a lower dividend term which is zero here, thus we have r 0 < 2v. In the case i = i 0 -1, then q ≤ q 0 ≤ q+1, thus r = βu-qv satisfies r < 3v.

As a consequence of Theorem 1, we can decide the correct rounding except when the last sh-1 bits from q are 000...000, 111...111, or 111...110, which occurs with probability less than 0.15% for the binary64 format. If the rare cases where we cannot decide the correct rounding, we compute r = βu -qv, subtract v and increment q at most two times until r < v, and deduce the round and sticky bits from the updated quotient q and remainder r.

The upper bound q + 2 is tight: it is attained for (u, v) = (10, 12) in radix β = 18, and for (u, v) = (21, 24) in powerof-two radix β = 32, taking i = i 0 -1 in step 1.

For the 2-limb division, we use a similar algorithm (Div-Approx2), which approximates the integer quotient of the significands. Since this algorithm is not specific to radix 2 64 , we describe it for a general radix β.

Algorithm DivApprox2 first computes a lower approximation q 1 of the upper quotient word, from the upper word u 1 of the dividend and the approximate reciprocal of the upper word

Algorithm 2 DivApprox2 Input: integers u, v with 0 ≤ u < v and β 2 /2 ≤ v < β 2 Output: approximate quotient q of uβ 2 /v 1: write u = u 1 β +u 0 and v = v 1 β +v 0 with 0 ≤ u i , v i < β 2: compute x satisfying x ≤ β 2 /(v 1 + 1) ≤ x + 1 3: q 1 = u 1 x/β 4: r = u -q 1 v/β
(equals (uβ -q 1 v)/β ) 5: q 0 = rx/β 6: return q = q 1 β + q 0 v 1 of the divisor. The end of Algorithm DivApprox2 can be seen as an integer version of Karp-Markstein's algorithm for the division [START_REF] Karp | High-precision division and square root[END_REF], or more simply, once the remainder has been approximated by r, we use the reciprocal again to approach the lower quotient word.

In step 2, x can be computed as follows: if

v 1 + 1 = β, take x = β, otherwise (assuming β is a power of two) v 1 + 1 cannot divide β 2 exactly, thus β 2 /(v 1 +1) = (β 2 -1)/(v 1 + 1)
, and we can take the same approximate reciprocal as in Algorithm DivApprox1 (with v 1 + 1 instead of v). 5 Write r = r 1 β +r 0 with 0 ≤ r 0 < β. As we will see later, the upper word satisfies r 1 ≤ 4. Step 6 becomes q = q 1 β + r 1 x + r 0 x/β , where the computation of r 1 x is done using a few additions instead of a multiplication.

Theorem 2: For β ≥ 5, the approximate quotient returned by Algorithm DivApprox2 satisfies

q ≤ uβ 2 v ≤ q + 21.
Proof. The input condition u < v implies u 1 ≤ v 1 , thus since x ≤ β 2 /(v 1 +1), it follows q 1 ≤ u 1 β/(v 1 +1) < β, thus q 1 fits in one limb. We have

β 2 = (v 1 +1)x+κ with 0 ≤ κ ≤ 2v 1 +1. Step 3 gives u 1 x = q 1 β + t, with 0 ≤ t < β. Then uβ -q 1 v = uβ -q 1 (v 1 β + v 0 ) = uβ -q 1 (v 1 + 1)β + q 1 (β -v 0 ) = uβ -(v 1 + 1)(u 1 x -t) + q 1 (β -v 0 ) = uβ -(β 2 -κ)u 1 + (v 1 + 1)t + q 1 (β -v 0 ) = u 0 β + κu 1 + (v 1 + 1)t + q 1 (β -v 0 ). Write v 1 = αβ with 1/2 ≤ α < 1. Thus uβ -q 1 v ≤ (β -1)β + (2αβ + 1)(αβ) +(αβ + 1)(β -1) + (β -1)β < β 2 (2 + α + 2α 2 ).
This proves that r < 5β at step 4. From r < β(2 + α + 2α 2 ) it follows q 0 ≤ rx/β < (2 + α + 2α 2 )x, thus since x ≤ β 2 /(v 1 + 1) < β/α, q 0 < (2/α + 1 + 5 In practice, like in Algorithm DivApprox1, x = β + x is split between one implicit upper bit β and a lower part x , that for simplicity we do not distinguish here. In the real code, step 3 therefore becomes q 1 = u 1 + u 1 x /β , and step 5 becomes q 0 = r + rx /β . 2α)β (note that q 0 might exceed β). Let rx = q 0 β + δ q , and uβ -q 1 v = βr + δ r , with 0 ≤ δ q , δ r < β. Then

uβ 2 -qv = β(βr + δ r ) -q 0 v = β 2 r + βδ r -q 0 (v 1 + 1)β + q 0 (β -v 0 ) = β 2 r + βδ r -(rx -δ q )(v 1 + 1) + q 0 (β -v 0 ) = β 2 r + βδ r -r(β 2 -κ) + δ q (v 1 + 1) + q 0 (β -v 0 ) = βδ r + rκ + δ q (v 1 + 1) + q 0 (β -v 0 ).
Bounding βδ r by β 2 -β, rκ by β(2 + α + 2α 2 )(2αβ + 1), δ q (v 1 + 1) by (β -1)(αβ + 1), and q 0 (β -v 0 ) by (2/α + 1 + 2α)β 2 , it follows:

0 ≤ uβ 2 -qv < β 2 ( 2 α + 2 + 7α + 2α 2 + 4α 3 ) + β(2 + 2α 2 ).
Dividing by v ≥ αβ 2 we get:

0 ≤ uβ 2 v -q < 2 α 2 + 2 α + 7 + 2α + 4α 2 + 1 β ( 2 α + 2α). (2)
The right-hand side is bounded by 21 + 5/β for 1/2 ≤ α < 1, thus for β ≥ 5 we have q ≤ uβ 2 /v < q + 22.

The largest error we could find is 20, for example for β = 4096, u = 8298491, v = 8474666. 6 With β = 2 64 , on 10 7 random inputs, the proportion of inputs for which

r 1 = 0 is 29%, r 1 = 1 is 59%, r 1 = 2 is 12%, r 1 = 3 is 0.2%,
and r 1 = 4 is 0%; and the average difference between q and uβ 2 /v is 2.8.

F. Square Root

Like for the division, we first compute an approximate result for the square root, and if we are not able to get the correctly rounded result, we compute an exact square root with remainder from that approximation. We first outline the exact algorithm for 1 limb. In case the input exponent is odd, we shift the input significand u 0 by one bit to the right (this is possible without any bit lost because p < 64). Then given the (possibly shifted) limb u 0 , one computes an integer square root:

u 0 • 2 64 = r 2 0 + s, 0 ≤ s ≤ 2r 0 .
Since 2 62 ≤ u 0 < 2 64 , we have 2 126 ≤ u 0 • 2 64 < 2 128 , thus 2 63 ≤ r 0 < 2 64 : r 0 has exactly 64 significant bits. The round and sticky bits are then computed from the 64 -p bits of r 0 and from the remainder s.

The expensive part is the computation of the (approximate or exact) integer square root. Since MPFR has to be independent of the machine floating-point hardware, we should use integer operations only. GMP's mpn_sqrtrem function implements the algorithm described in [START_REF] Zimmermann | Karatsuba square root[END_REF], which relies on integer division. We can do better by first computing an approximate reciprocal square root.

1) Approximate Reciprocal Square Root: Algorithm 3 uses an integer variant of Newton's iteration for the reciprocal square root: v 0 is a 11-bit value such that x 0 := v 0 /2 10 approximates the root of a 0 := d 10 /2

d 37 = 2 -27 d + 1 4: e 0 = 2 57 -v 2 0 d 37 5: v 1 = 2 11 v 0 + 2 -47 v 0 e 0 6: e 1 = 2 79 -v 2 1 d 37 7: v 2 = 2 10 v 1 + 2 -70 v 1 e 1 8: e 2 = 2 126 -v 2 2 d 9: v 3 = 2 33 v 2 + 2 -94 v 2 e 2
value such that x 1 := v 1 /2 21 approximates the root of a 1 := d 37 /2 37 , v 2 is a 32-bit value such that x 2 := v 2 /2 31 approximates the root of a 2 := d 37 /2 37 , and v 3 is a 65bit value such that x 3 := v 3 /2 64 approximates the root of a 3 := d/2 64 .

Theorem 3: The value v 3 returned by Algorithm 3 differs by at most 8 from the reciprocal square root:

v 3 ≤ s := 2 96 / √ d ≤ v 3 + 8.
Lemma 1: Assume positive real numbers x 0 , x 1 , ..., x n are computed using the following recurrence, with x 0 ≤ a -1/2 0 :

x k+1 = x k + x k 2 (1 -a k+1 x 2 k ), (3)
where a 0 ≥ a 1 ≥ • • • ≥ a n ≥ a > 0, then:

x 1 ≤ x 2 ≤ • • • ≤ x n ≤ a -1/2 .
Proof. It follows from [1, Lemma 3.14] that x k+1 ≤ a -1/2 k+1 . Together with x 0 ≤ a -1/2 0 , this gives 1 -a k x 2 k ≥ 0 for all k. Since a k+1 ≤ a k , it follows that 1 -a k+1 x 2 k ≥ 1 -a k x 2 k ≥ 0. Put into (3) this proves x k+1 ≥ x k , and thus the lemma since x k+1 ≤ a -1/2 k+1 and a k+1 ≥ a imply x k+1 ≤ a -1/2 . This lemma remains true when the correction term x k 2 (1a k+1 x 2 k ) in Eq. ( 3) -which is thus non-negative -is rounded down towards zero. By applying this result to Proof of Theorem 3. The construction of the lookup table ensures that v 0 < 2 11 and e 0 ≥ 0 at step 4 of the algorithm. By exhaustive search on all possible d 10 values, we get:

x 0 = v 0 /2 10 , x 1 = v 1 /2 21 , x 2 = v 2 /2 31 , x 3 = v 3 /2 64 , with a 0 = d 10 /2 10 ≥ a 1 = d 37 /2 37 = a 2 ≥ a 3 = d/2 64 , it follows 1 ≤ v 0 /2 10 ≤ v 1 /2 21 ≤ v 2 /2 31 ≤ v 3 /2 64 ≤ 2 32 / √ d ≤ 2, thus 2 10 ≤ v 0 < 2 11 , 2 21 ≤ v 1 < 2 22 , 2 31 ≤ v 2 < 2 32 , 2 64 ≤ v 3 < 2 65
0 ≤ e 0 < 2 49.263 .

Thus at step 4, we can compute 2 57 -v 2 0 d 37 in 64-bit arithmetic, which will result in a number of at most 50 bits, and step 5 can be done in 64-bit arithmetic, since the product v 0 e 0 will have at most 61 bits.

Let δ 1 be the truncation error in step 5, with 0 ≤ δ 1 < 1:

2 -47 v 0 e 0 = 2 -47 v 0 e 0 + δ 1 .

Then: 37 , and:

e 1 = 2 79 -v 2 1 d 37 = 2 79 -d 37 ((v 1 + δ 1 ) -δ 1 ) 2 = 2 79 -d 37 (v 1 + δ 1 ) 2 + d 37 δ 1 (2v 1 + δ 1 ). Since v 2 1 d 37 ≤ 2 79 -because e 1 ≥ 0 -, then v 1 d 37 ≤ 2 79 /v 1 ≤ 2 58 because 2 21 ≤ v 1 . It follows that γ 1 := d 37 δ 1 (2v 1 + δ 1 ) satisfies 0 ≤ γ 1 < 2 59 + 2
e 1 -γ 1 = 2 79 -d 37 (v 1 + δ 1 ) 2 = 2 79 -d 37 (2 11 v 0 + 2 -47 v 0 e 0 ) 2 = 2 79 -d 37 v 2 0 (2
2 -70 v 1 e 1 = 2 -70 v 1 e 1 + δ 2 .
Then: By doing again this error analysis for a given value of d 37 , we get a tighter bound involving the δ i truncation errors. For example we can bound d 37 δ 1 (2v 1 + δ 1 ) by d 37 (2v 1,max + 1), where v 1,max = 2 79 /d 37 . This yields a finer bound. By exhaustive search on all possible d 37 values, we found the maximal error is at most 8.

e 2 = 2 126 -v 2 2 d = 2 126 -d((v 2 + δ 2 ) -δ 2 ) 2 = 2 126 -d(v 2 + δ 2 ) 2 + dδ 2 (2v 2 + δ 2 ). Since v 2 ≥ 2 31 and v 2 2 d ≤ 2 126 , v 2 d ≤ 2 95 , thus dδ 2 (2v 2 + δ 2 ) < 2 96 + 2 64 . e 2 -(2 96 + 2 64 ) ≤ 2 126 -d(v 2 + δ 2 ) 2 = 2 126 -d(2 10 v 1 + 2 -70 v 1 e 1 ) 2 = 2 126 -dv 2 1 (2
e 3 = 2 192 -v 2 3 d = 2 192 -d((v 3 + δ 3 ) -δ 3 ) 2 = 2 192 -d(v 3 + δ 3 ) 2 + dδ 3 (2v 3 + δ 3 ). Since v 2 3 d ≤ 2 192 , and v 3 ≥ 2 64 , it follows v 3 d ≤ 2 128 , thus γ 3 := dδ 3 (2v 3 + δ 3 ) < 2 129 + 2 64 . e 3 -γ 3 = 2 192 -d(v 3 + δ 3 ) 2 = 2 192 -d(2 33 v 2 + 2 -94 v 2 e 2 ) 2 = 2 192 -
The bound of 8 is optimal, since it is attained for d = 4755801239923458105.

Remark 1: one can replace 2 -94 v 2 e 2 in step 9 by 2 -29 v 2 e 2 , where e 2 = 2 -65 e 2 has at most 32 bits, like v 2 , and thus the product v 2 e 2 can be computed with a low 64-bit product, which is faster than a full product giving both the low and high words. If we write e 2 = 2 65 e 2 + r with 0 ≤ r < 2 65 , then:

2 -94 v 2 e 2 -2 -29 v 2 e 2 = 2 -94 v 2 r < 2 3 .
This increases the maximal error from 8 to 15, since for the only value d 37 = 35433480442 that can yield c ≥ 8 with the original step 9, we have c < 8.006 and v 2 = 4229391409, thus c + 2 -94 v 2 r < 15.9.

2) The mpfr_sqrt1 function: The mpfr_sqrt1 function computes the square root of a 1-limb number n using Algorithm 4 (SqrtApprox1). In step 1, it computes a 32-bit approximation of 2 63 / √ n using the value v 2 of Algorithm 3 (called with d = n), then uses this approximation to deduce the exact integer square root y of n (step 2) and the corresponding remainder z (step 3), and finally uses again x to approximate the correction term t (step 4) to form the approximation s in step 5. Steps 2 and 3 can be implemented as follows. First compute an initial y = 2 -32 x 2 -31 n and the corresponding remainder z = n -y 2 . Then as long as z ≥ 2y + 1, subtract 2y + 1 to z and increment y by one. Note that since x 2 n < 2 126 , xn/2 31 < 2 95 /n ≤ 2 32 , thus x 2 -31 n fits on 64 bits.

Algorithm 4 SqrtApprox1

Input: integer n with 2 62 ≤ n < 2 64 Output: integer s approximating √ 2 64 n 1: compute an integer x approximating 2 63 / √ n with

x ≤ 2 63 / √ n 2: y = √ n (using the approximation x) 3: z = n -y 2 4: t = 2 -32 xz 5: s = y • 2 32 + t Theorem 4: If the approximation x in step 1 is the value v 2 of Algorithm 3, then Algorithm SqrtApprox1 returns s satisfying s ≤ √ 2 64 n ≤ s + 7. Proof Write n = α2 64 with 1/4 ≤ α < 1, x = 2 63 / √ n -δ x , y =
√ n -δ y and t = 2 -32 xz -δ t with 0 ≤ δ y , δ t < 1. Since x 2 n > 2 126 -2 96.338 from Eq. ( 4), x √ n > 2 63 -2 32.339 , thus

2 63 -x √ n = δ x √ n < 2 33.339 : 2 64 n -s 2 = 2 64 n -(y • 2 32 + t) 2 = 2 64 (n -y 2 ) -2 33 yt -t 2 = 2 64 z -2 33 yt -t 2 . (5)
We first bound 2 64 n-s 2 by above, assuming it is non-negative:

2 33 yt = 2xyz -2 33 δ t y = (2 64 / √ n -2δ x )( √ n -δ y )z -2 33 δ t y (6)
thus since δ x √ n < 2 33.339 and δ y , δ t < 1:

2 64 z -2 33 yt ≤ (2δ x √ n + 2 64 / √ n)z + 2 33 y = 2 32 z(2 2.339 √ α + 1/ √ α) + 2 65 √ α,
where we used y ≤ √ n = 2 32 √ α. Also:

z = n -( √ n -δ y ) 2 ≤ 2 √ nδ y ≤ 2 33 α 1/2 .
Substituting this in the bound for 2 64 n -s 2 gives:

2 64 n -s 2 ≤ 2 64 z -2 33 yt ≤ 2 65 f (α), with f (α) = 2 2.339 α + 1 + α 1/2 . Let c ≥ 0 such that 2 64 n = (s+c) 2 . Then 2 64 n-s 2 = 2sc+c 2 , which implies 2sc < 2 65 f (α), thus c < 2 64 f (α)/s. (7)
Since s ≥ 2 63 , and the maximum of f (α) for 1/4 ≤ α ≤ 1 is less than 7.06 (attained at α = 1), we get c < 14.12. Now this gives s >

√ 2 64 n -15 > √ α2 64 /1.01, therefore c < 1.01 • f (α)/ √ α. Now the function f (α)/ √
α is bounded by 7.06 for 1/4 ≤ α ≤ 1, the maximum still attained at α = 1. Therefore c < 1.01 • 7.06 < 7.14, which proves the upper bound 7.

Now assume 2 64 n -s 2 < 0. Since 2 64 z -2 33 yt ≥ 0 by Eq. ( 6), we have 2 64 n -s 2 ≥ -t 2 by Eq. ( 5). Since x ≤ 2 63 / √ n and z ≤ 2 33 √ α:

t ≤ 2 -32 xz < 2 32 .
(The last inequality is strict because we can have t = 2 32 only when x = 2 63 / √ n, which can only happen when n = 2 62 , but in that case z = 0.) This proves that the product xz at step 4 can be computed modulo 2 64 . Now write 2 64 n = (s -c) 2 with c > 0:

2 64 n -s 2 ≥ -t 2 > -2 64 (s -c) 2 -s 2 > -2 64 2sc < 2 64 + c 2 .
This inequality has no solutions for 2 63 ≤ s -c < 2 64 . Indeed, since we assumed 2 64 n -s 2 < 0, this implies s > 2 63 , because for s = 2 63 we have s 2 ≤ 2 64 n. But then if s = 2 63 + u with u ≥ 1 we would need 2 64 c + 2uc < 2 64 + c 2 , which we can rewrite as 2u < 2 64 /c+c-2 64 . For 1 ≤ c ≤ 2 64 , the expression 2 64 /c + c is bounded by 2 64 + 1, which yields 2u < 1, having no integer solutions u ≥ 1.

3) The mpfr_sqrt2 function: The mpfr_sqrt2 function first computes a 128-bit approximation of the square root using Algorithm 5, where in step 1 we can use Algorithm 3, with the variant of Remark 3, achieving the bound δ = 16. Al- III. FAITHFUL ROUNDING In addition to the usual rounding modes, faithful rounding (MPFR_RNDF) will be supported in MPFR 4 by the most basic operations. We say that an operation is faithfully rounded if its result is correctly rounded toward -∞ or toward +∞. The actual rounding direction is indeterminate: it may depend on the argument, on the platform, etc. (except when the true result is exactly representable, because in this case, this representable value is always returned, whatever the rounding mode). From this definition, the error on the operation is bounded by 1 ulp, like in the directed rounding modes. Moreover, the ternary value is unspecified for MPFR_RNDF.

Algorithm 5 SqrtApprox2 Input: integer n = n 1 • 2 64 + n 0 with 2 62 ≤ n 1 <
The advantage of faithful rounding is that it can be computed more quickly than the usual roundings: computing a good approximation (say, with an error less than 1/2 ulp) and rounding it to the nearest is sufficient. So, the main goal of MPFR_RNDF is to speed up the rounding process in internal computations. Indeed, we often do not need correct rounding at this level, just a good approximation and an error bound.

In particular, MPFR_RNDF allows us to completely avoid the Table Maker's Dilemma (TMD), either for the rounding or for the ternary value.

MPFR provides a function mpfr_can_round taking in entry: an approximation to some real, unknown value; a corresponding error bound (possibly with some given direction); a target rounding mode; a target precision. The goal of this function is to tell the caller whether rounding this approximation to the given target precision with the target rounding mode necessarily provides the correct rounding of the unknown value. For correctness, it is important that if the function returns true, then correct rounding is guaranteed (no false positives). Rare false negatives could be acceptable, but for reproducibility, it was chosen to never return a false negative with the usual rounding modes. When faithful rounding was introduced, a new semantic had to be chosen for the MPFR_RNDF rounding mode argument: true means that rounding the approximation in an adequate rounding mode can guarantee a faithful rounding of the unknown value. Reproducibility was no longer important as in general with faithful rounding. One reason to possibly accept false negatives here was to avoid the equivalent of the TMD for this function. However, a false negative would often mean a useless recomputation in a higher precision. Since it is better to spend a bit more time (at most linear) to avoid a false negative than getting a very costly reiteration in a higher precision, it was chosen to exclude false negatives entirely, like with the other rounding modes.

IV. EXPERIMENTAL RESULTS

We used a 3.2Ghz Intel Core i5-6500 for our experiments (Skylake microarchitecture), with only one core enabled to decrease turbo-boost effects, and revision 11249 of MPFR (aka MPFR 4.0-dev), and GCC 6.2.1 under Debian GNU/Linux testing (stretch). MPFR was configured with GMP 6.1.2 [START_REF] Granlund | GNU MP: The GNU Multiple Precision Arithmetic Library[END_REF] and --disable-shared. The number of cycles were measured with the tools/mbench utility distributed with MPFR. To get stable results, we called the mbench binary with numactl --physcpubind=0, that binds it to cpu 0, and we ensured no other processes were running on the machine.

Table I compares the number of cycles for basic arithmetic operations -with rounding to nearest -between MPFR 3 and the upcoming version 4, with 53 and 113 bits of precision, corresponding to the binary64 and binary128 formats respectively. The speedup goes from a factor 1.7 (for the 113bit division) to a factor 3.6 (for the 53-bit squaring). This improvement of basic arithmetic automatically speeds up the computation of mathematical functions, as demonstrated in Table IV (still for rounding to nearest). Indeed, the computation of a mathematical function reduces to basic arithmetic for the argument reduction or reconstruction, for the computation of Taylor series or asymptotic series, ... Here the speedup goes from 15% (for the 113-bit arc-tangent) to a factor 2.2 (for the 53-bit power).

Timings for the mpfr_sum function in the faithful branch7 are done with a special test (misc/sum-timings) because the input type, an array of MPFR numbers, is specific to this function, and we need to test various kinds of inputs. The inputs are chosen randomly with fixed parameters, corresponding to the first 5 columns 8 some cancellation, 2 even more cancellation), the precision of the input numbers (here, all of them have the same precision), the precision of the result, and an exponent range e (1 means that the input numbers have the same order of magnitude, and a large value means that the magnitudes are very different: from 1 to 2 e-1 ). A test consists in several timings on the same data, and we consider the average, the minimum and the maximum, and if the difference between the maximum and the minimum is too large, the result is ignored. Due to differences in timings between invocations, 4 tests have been run. On some parameters, MPFR_RNDF is faster than the other rounding modes (its maximum timing is less than the minimum timing for the usual rounding modes): for a small array size, small output precision and no cancellation (first line of Table III), there is a small gain; when the input precision is small and the range is large (so that the input values do not overlap in general, meaning that the TMD occurs in the usual rounding modes), we can typically notice a 25% to 40% gain when there is no cancellation (lines 2, 5 and 6), and a bit less when a cancellation occurs (lines 3 and 4). MPFR_RNDF is never slower. Note that the kind of inputs have a great influence on a possible gain, which could be much larger than the one observed here (these tests have not been specifically designed to trigger the TMD).

V. CONCLUSION

This article presents algorithms for fast floating-point arithmetic with correct rounding on small precision numbers. The implementation of those algorithms in version 4 of MPFR gives a speedup of up to 3.6 with respect to MPFR 3. As a consequence, the computation of mathematical functions is greatly improved.

  rb = a0 & (ONE << (sh -1)); ap[0] = a0 ^rb; sb = 0; The case bx > cx is treated similarly: we have to shift cp[0] by d = bx -cx bits to the right. We distinguish three cases: d < sh where the shifted cp[0] fully overlaps with ap[0], sh ≤ d < 64, where the shifted cp[0] partly overlaps with ap[0], and d ≥ 64, where the shifted cp[0] does not overlap with ap[0]. In the last case bx < cx, we simply swap the inputs and reduce to the bx > cx case. Now consider the rounding. All rounding modes are first converted to nearest, toward zero or away from zero. At this point, ap[0] contains the current significand, bx the exponent, rb the round bit and sb the sticky bit. No underflow is possible in the addition, since the inputs have same sign: |b + c| ≥ min(|b|, |c|

  rb = a0 & (ONE << (sh -1)); sb |= (a0 & mask) ^rb; ap[0] = a0 & ~mask;

  (the right bounds are strict because in the only case where 2 32 / √ d = 2, i.e., d = 2 62 , we have v 3 = 2 65 -3).

  and we shift a to the left by the number of leading zeros of a 1 :

	count_leading_zeros (cnt, a1);
	if (cnt) {
	ap[1] = (a1 << cnt) | (a0 >> (64-cnt));
	a0	= (a0 << cnt) | (sb >> (64-cnt));
	sb <<= cnt;
	bx -= cnt; }
	else
	ap[1] = a1;

  10 , v 1 is a 22-bit Algorithm 3 RecSqrtApprox1 Input: integer d with 2 62 ≤ d < 2 64 Output: integer v 3 approximating s = 2 96 / √ d 1: d 10 = 2 -54 d + 1 2: v 0 = 2 30 /d 10 (table lookup) 3:

  11 + 2 -47 e 0 ) 2 . Now, using v 2 0 d 37 = 2 57 -e 0 :e 1 -γ 1 = 2 79 -(2 57 -e 0 )(2 11 + 2 -47 e 0 ) 2Since e 0 < 2 49.263 we deduce 3/4 + 2 -59 e 0 < 2 -0.412 and 0 ≤ e 1 < 2 -35 2 98.526 2 -0.412 + 2 59 + 2 37 < 263.196 .

	= 2 79 -(2 57 -e 0 )(2 22 + 2 -35 e 0 + 2 -94 e 2 0 )
	= 2 79 -(2 79 -2 -35 e 2 0 + 2 -37 e 2 0 -2 -94 e 3 0 )
	= 2 -35 e 2 0 (3/4 + 2 -59 e 0 ).

Therefore e 1 < 2 64 and e 1 can be computed using integer arithmetic modulo 2 64 . Since d 10 = 2 -27 (d 37 -1) + 1, e 1 only depends on d 37 , so that we can perform an exhaustive search on the 2 37 -2 35 possible values of d 37 . By doing this, we find that the largest value of e 1 is obtained for d 37 = 33 • 2 30 + 1, which corresponds to d 10 = 265; this gives e 1 = 10263103231123743388 < 2 63.155 . Let δ 2 be the truncation error in step 7, 0 ≤ δ 2 < 1:

  10 + 2 -70 e 1 ) 2 . Now writing d = 2 27 d 37 -ρ with 0 < ρ ≤ 2 27 , using d 37 v 2 1 = 2 79 -e 1 , and writing ε = 2 96 + 2 64 + ρv 2 1 (2 10 + 2 -70 e 1 ) 2 : e 2 -ε ≤ 2 126 -2 27 (2 79 -e 1 )(2 10 + 2 -70 e 1 ) 2 Here again, an exhaustive search is possible, since for a given value of d 37 , e 2 is maximal when d = 2 27 (d 37 -1): the largest value of e 2 is obtained for d 37 = 132607222902, corresponding to d 10 = 989, with e 2 = 81611919949651931475229016064 < 2 96.043 . The final error is estimated using the truncation error δ 3 in step 9: 2 -94 v 2 e 2 = 2 -94 v 2 e 2 + δ 3

	= 2 126 -(2 79 -e 1 )(2 47 + 2 -32 e 1 + 2 -113 e 2 1 )
	= 2 126 -(2 126 -2 -32 e 2 1 + 2 -34 e 2 1 -2 -113 e 3 1 )
	= 2 -32 e 2 1 (3/4 + 2

-81 e 1 ).

Since e 1 < 2 63.155 we deduce 3/4 + 2 -81 e 1 < 2 -0.415 , and:

ρv 2 1 (2

10 + 2 -70 e 1 ) 2 ≤ 2 71 (2 20 + 2 5 + 2 -12 ). Thus ε ≤ 2 96 + 2 64 + 2 91 + 2 76 + 2 59 and: 0 ≤ e 2 < 2 -32 2 126.31 2 -0.415 + ε < 2 96.338 . (4)

  dv 2 2 (2 33 + 2 -94 e 2 ) 2 . Now since dv 2 2 = 2 126 -e 2 : e 3 -γ 3 ≤ 2 192 -(2 126 -e 2 )(2 33 + 2 -94 e 2 ) 2 = 2 192 -(2 126 -e 2 )(2 66 + 2 -60 e 2 + 2 -188 e 2 2 ) = 2 192 -(2 192 -2 -60 e 2 2 + 2 -62 e 2 2 -2 -188 e 3 2 ) = 2 -60 e 2 2 (3/4 + 2 -128 e 2 ). Since e 2 < 2 96.043 we deduce 3/4 + 2 -128 e 2 < 2 -0.415 : 0 ≤ e 3 < 2 -60 2 192.086 2 -0.415 + 2 129 + 2 64 < 2 131.882 . Again by exhaustive search, restricted on the values of d 37 that give a sufficient large value of e 2 , we found the maximal value of e 3 satisfies e 3 < 2 131.878 . Now, given v 3 returned by Algorithm 3, let c ≥ 0 be the real number such that (v 3 + c) 2 d = 2 192 . Since 2 192 -v 2 3 d < 2 131.878 , it follows 2v 3 cd < 2 131.878 , thus c < 2 131.878 /(2v 3 d). Since v 2 3 d = 2 192 -e 3 , and v 3 < 2 65 , we have v 3 d ≥ (2 192 -2 131.878 )/2 65 > 2 126.999 . Therefore c < 2 3.879 < 15.

  ≤ n 1 x/2 64 ≤ nx/2 128 < √ n. As a consequence, we have z > 0 at step 3.The proof of Theorem 5 is very similar to that of Theorem 4, and can be found in the appendix.The largest difference we found is 24, for example with n = 2 64 • 18355027010654859995 -1, taking x = 2 96 / √ n 1 -15, whereas the actual code might use a larger value. The smallest difference we found is -1, for example with n = 2 64 • 5462010773357419421 -1, taking x = 2 96 / √ n 1 .

				By construction, we have x ≤ 2 96 / √	n 1 + 1 < 2 128 / √	n,
				therefore y
				2 64 and
	0 ≤ n 0 < 2 64 Output: integer s approximating 1: compute an integer x approximating 2 96 / √ 2 128 n √ n 1 with
	0 ≤ 2 96 /	√	n 1 -x < δ := 16
	2: y = 2 32 √ n 1			(using the approximation x)
	3: z = n -y 2
	4: t = xz/2 65
	5: s = y • 2 64 + t
	gorithm SqrtApprox2 is an integer version of Karp-Markstein's
	algorithm for the square root [9], which incorporates n/2 128
	in Newton's iteration for the reciprocal square root.
	Theorem 5: Algorithm SqrtApprox2 returns s satisfying
	s -4 ≤	√	2 128 n ≤ s + 26.

  of Table III: the size of the array (number of input MPFR numbers), the number of cancellation terms in the test (0 means no probable cancellation, 1 means

		Parameters		RND *	RNDF
	10 1	0 10 7	10 1	1	411	399
	10 3	0 10 1	10 5	10 8	27216	20366
	10 3	1 10 1	10 5	10 8	39639	32898
	10 3	2 10 1	10 5	10 8	44025	35276
	10 5	0 10 1	10 1	10 8	1656988 1034802
	10 5	0 10 1	10 3	10 8	1393447	833711
			TABLE III
	NUMBER OF CYCLES FOR M P F R_S U M WITH THE USUAL ROUNDING
		MODES AND MPFR_RNDF.
			MPFR 3.1.5	MPFR 4.0-dev
		precision	53	113	53	113
		mpfr_exp	3996	6787	2352	4253
		mpfr_sin	3746	4907	2684	3838
		mpfr_cos	2853	3553	1872	2819
		mpfr_log	4473	7086	2325	5993
		mpfr_atan	13937 20182	7175	11385
		mpfr_pow 12915 18393	5704	11222
			TABLE IV
	NUMBER OF CYCLES FOR MATHEMATICAL FUNCTIONS.

https://octave.sourceforge.io/interval/

https://github.com/nehmeier/libieeep1788

http://www.sagemath.org

MPFR does not have subnormals, but provides a function to emulate them for IEEE 754 support.

If β 2 + 1 has no divisors in [β/2 + 1, β], as for β = 2 32 and β = 2 64 , then we cannot have κ = 2v 1 + 1, the 1/β term disappears in Eq. (2), and the bound becomes uβ 2 /v < q + 21.

branches/faithful revision 11121.

These are the arguments of position 2 to 6 of sum-timings.

ACKNOWLEDGEMENTS The authors thank Patrick Pélissier who designed the mbench utility, which was very useful in comparing different algorithms and measuring their relative efficiency.

APPENDIX

Proof of Theorem 5. Write n 1 = α2 64 with 1/4 ≤ α < 1. Write x = 2 96 / √ n 1 -δ x , y = 2 32 √ n 1 -δ y and t = xz/2 65δ t with 0 ≤ δ x < δ, and 0 ≤ δ y , δ t < 1.

We first bound 2 128 n-s 2 by above, assuming 2 128 n-s 2 ≥ 0:

thus

where we used y ≤

Substituting Eq. ( 11) in Eq. ( 10) yields:

Substituting in Eq. ( 8) gives

Let c be the real number such that 2 128 n = (s + c) 2 . Then since 2 128 n -

Since s ≥ 2 127 , and the maximum of f (α) for 1/4 ≤ α ≤ 1 is 26. 9) that 2 128 z -2 65 yt ≥ 0. Eq. ( 8) therefore yields 2 128 n -s 2 ≥ -t 2 , and we have to bound t 2 . Using x ≤ 2 96 / √ n 1 , Eq. ( 11) and

.

Writing 2 128 n = (s -c) 2 with c > 0 yields:

This implies for all 1/4 ≤ α < 1: c < c 2 2 128 + 4, which proves c < 4.01, thus the lower bound is s -4.


Yong Han
email: hanyong@amss.ac.cn

Yuefei Wang
email: wangyf@math.ac.cn

Michel Zinsmeister

On The Brownian Loop Measure

Keywords: Brownian loop, SLE bubble, Brownian bubble, Disconnect from boundary

In 2003 Lawler and Werner introduced the Brownian loop measure and studied some of its properties. Cardy and Gamsa has predicted a formula for the total mass of the Brownian loop measure on the set of simple loops in the upper half plane and disconnect two given points from the boundary. In this paper we give a rigorous proof of the formula.

Introduction

The conformally invariant scaling limits of a series of planar lattice models can be described by the one-parameter family of random fractal curves SLE(κ), which was introduced by Schramm. These models include site percolation on the triangular graph, loop erased random walk, Ising model, harmonic random walk, discrete Gaussian free field, FK-Ising model and uniform spanning tree. Using SLE as a tool, many problems related to the Theorem 1. Denote by µ loop H the Brownian loop measure on the upper half plane and by γ a sample of the Brownian loop. Given two points z = x + iy, w = u + iv ∈ H, let E(z, w) denote the event that γ disconnects both z and w from the boundary of H. Then we have

µ loop H [E(z, w)] = - π 5 √ 3 - 1 10 η 3 F 2 (1, 4 3 , 1; 5 3
, 2; η) -1 10 log(η(η -1)) + Γ( 2 3 ) 2 5Γ( 43 ) (η(η -1))

1 3 2 F 1 (1, 2 3 ; 4 3
, η).

(1.1)

where η = η(z, w) = - (x -u) 2 + (y -v) 2 4yv , (1.2)
and 3 F 2 , 2 F 1 are the hypergeometric functions.

The formula (1.1) was first given by Cardy and Gamsa using conformal field theory which assumes that O(n) model has the scaling limit. In fact (1.1) has a nicer form:

µ loop [E(z, w)] = - 1 10 [log σ + (1 -σ) 3 F 2 (1, 4 3 , 1; 5 3 , 2; 1 -σ)], (1.3)
where

σ = σ(z, w) = |z -w| 2 |z -w| 2 = (x -u) 2 + (y -v) 2 (x -u) 2 + (y + v) 2 , (1.4)
and 3 F 2 is the hypergeometric function. Remark 2. By the conformal invariance of Brownian loop measure (see [START_REF] Gregory | The brownian loop soup[END_REF]), for any simply connected domain D ⊂ C with z, w ∈ D, we can get the total mass of the Brownian loop in D that disconnect both z and w from ∂D by the conformal map from D to H. In particular, if D = D, we choose the conformal map ϕ(z) = i 1+z 1-z from D onto H. Then the total mass of the Brownian loop measure in

D that disconnects z, w ∈ D from ∂D is - 1 10 [log σ + (1 -σ) 3 F 2 (1, 4 3 , 1; 5 3 , 2; 1 -σ)],
where σ = σ(z, w)

= |z-w| 2 |1-z w| 2 .

Background

In this section very brief introductions on the chordal SLE, Brownian loop measure and SLE bubble measure are given.

Chordal SLE process

The chordal SLE process from 0 to ∞ in H is a random family of conformal maps (g t : t ≥ 0) that satisfies

∂ t g t (z) = 2 g t (z) - √ κB t , 0 ≤ t < τ (z), g 0 (z) = z. (2.1)
where B is a standard Brownian motion and τ (z) is the blow-up time for the differential equation (2.1). The SLE(κ) process is generated by a continuous curve γ, (see [START_REF] Rohde | Basic properties of SLE[END_REF] and [START_REF] Gregory F Lawler | Conformal invariance of planar loop-erased random walks and uniform spanning trees[END_REF]). That is, the following limits exists:

γ(t) = lim y↓0 g -1 t (iy + √ κB t ).
The curve γ is a simple curve if and only if κ ∈ (0, 4], (see [START_REF] Rohde | Basic properties of SLE[END_REF]). The SLE(κ) curve satisfies the scaling invariance property (see [START_REF] Rohde | Basic properties of SLE[END_REF]), i.e for any r > 0, rγ has the same distribution as γ (with a time rescaling). So given any triple set (D, a, b), where D is a simply connected domain with two given boundary points a and b, we can define the SLE process on D from a to b by the conformal map from H to D that sends 0 to a and ∞ to b.

The following lemma will be useful in defining the SLE-bubble measure.

Lemma 3 (see [START_REF] Schramm | A percolation formula[END_REF]). Given z = x + iy ∈ H, suppose κ ∈ (0, 4] and γ is the SLE(κ) curve from 0 to ∞ in H. Then the probability that γ passes the left of z is

p(z) = C ∫ x y -∞ (1 + t 2 ) -4 κ dt, (2.2)
where C = C(κ) is the constant that make the total integral above equal to 1.

Brownian loop measure and bubble measure

In this section, we will introduce several measures on the space of continuous curves in the plane. To keep the present article short, we will not provide the detailed discussions but instead refer the reader to the fifth chapter of Lawler's book [START_REF] Gregory | Conformally invariant processes in the plane[END_REF] and [START_REF] Gregory | The brownian loop soup[END_REF]. Let µ(z; t) be the law of a complex Brownian motion (B s : 0 ≤ s ≤ t) starting from z. And µ(z; t) can be written as

µ(z; t) = ∫ C µ(z, w; t)dw
where the above integral can be regarded as the integral of measure-valued functions. Here µ(z; t) and µ(z, w; t) are regarded as measures on the space of curves in the plane. Using the density function of the complex Brownian motion, we can see that the total mass of µ(z, w; t) is 1 2πt exp{-1 2t |z -w| 2 }. Let µ(z, w) be the measure defined by µ(z, w) = ∫ ∞ 0 µ(z, w; t). This is a σ-finite infinite measure. If D ⊂ C is a simply connected domain with nice boundary and z, w ∈ D, we can define µ D (z, w) be the restriction of µ(z, w) on the space of curves that lie inside D.

If z ̸ = w, the total measure of µ D (z, w) is πG D (z, w), where G D (z, w) is the Green function on D.
If D is a simply connected domain with nice boundary, let B be a complex Brownian motion starting from z ∈ D and τ D the exit time. Denote µ D (z, ∂D) the law of (B t : 0 ≤ t ≤ τ D ), we can write

µ D (z, ∂D) = ∫ ∂D µ D (z, w)dw.
Here µ D (z, w) can be regarded as a measure on the space of curves in D from z to w ∈ ∂D, and the total mass of µ D (z, w) is the the Poisson kernel H D (z, w). For z ∈ D, w ∈ ∂D, µ D (z, w) can also be equivalently defined by the limits

µ D (z, w) = lim ϵ→0 1 2ϵ µ D (z, w + ϵn w ),
where n w is the inner normal at w. And similarly for z, w ∈ ∂D, we can also define

µ D (z, w) = lim ϵ→0 1 2ϵ 2 µ D (z + ϵn z , w + ϵn w )
It can be showed that the above limits exist in the sense of Prohorov convergence(see Chapter 5 of [START_REF] Gregory | Conformally invariant processes in the plane[END_REF]).

Given z ∈ ∂D, the Brownian bubble measure µ bub D (z) is defined as the limit

µ bub D (z) := lim w∈∂D,w→z πµ D (z, w).
The Brownian loop measure is defined as following:

µ loop C := ∫ C 1 t γ µ(z, z)dz = ∫ C ∫ ∞ 0 1 t γ µ(z, z; t)dtdz.
Since µ(z, z) is a measure defined on loops with z as a marked point(called a root), the Brownian loop measure should be understood as the above integral of measures by forgetting the root. For any domain D, let µ loop D be the restriction of the Brownian loop measure on the space of loops inside D.

For any a ∈ R, define H a := {x + iy ∈ C : y > a}, according to the lowest point of the Brownian loop, the Brownian loop can be decomposed into the following integral of Brownian bubbles(see [START_REF] Gregory | The brownian loop soup[END_REF]):

µ loop C = 1 π ∫ C µ bub Hy (x + iy)dxdy. (2.3) (2.
3) will be very important in our computation.

SLE bubble measure

In this section we will define the SLE bubble measure and give the relation between SLE( 83 ) bubble and Brownian bubble measure. Suppose κ ∈ (0, 4], ϵ > 0 and γ ϵ is the SLE(κ) curve from 0 to ϵ in the upper half plane. Let µ ϵ denote the law of γ ϵ . Lemma 4. The limit of the following limit exists:

µ bub SLE(κ) (0) = lim ϵ→0 ϵ 1-8 κ µ ϵ . (2.4)
We call µ bub SLE(κ) (0) the SLE(κ)-bubble measure.

Proof. We only need to show that the limit restricted to some generated algebras that consist of finite mass exists. Here we choose the measurable sets {γ : γ disconnects z from ∞} for fixed z ∈ H. By the definition of the SLE(κ) from 0 to ϵ, we choose the auto-conformal map F ϵ (z) = ϵz z+1 that sends ∞ to ϵ and fixes 0. We have

P[γ disconnects z from ∞] = p(F -1 ϵ (z)),
where p(z) is defined in (2.2). Therefore

lim ϵ→0 ϵ 1-8 κ p(F -1 ϵ (z)) = Γ( 4 κ ) √ πΓ( 8-κ 2κ )( 8 κ -1) ( x 2 + y 2 y ) 1-8 κ (2.5)
So for fix z ∈ H, if we denote µ ϵ (z) the restriction of µ ϵ restricted to the curves that disconnect z from ∞, then by the above equation, we know that the limit

µ bub SLE(κ) (0, z) := lim ϵ→0 ϵ 1-8 κ µ ϵ (z)
exists and therefore we can define µ bub SLE(κ) (0) as the limit of µ bub SLE(κ) (0, z) as z tends to zero.

If κ = 8
3 , from (2.5), we get that the total mass of the SLE( 83 )-bubble that disconnects a given point z 2 which corresponds to the part (a) of proposition 3.1 in [START_REF] Beliaev | Some remarks on SLE bubbles and schramm's two-point observable[END_REF]. In fact, [START_REF] Beliaev | Some remarks on SLE bubbles and schramm's two-point observable[END_REF] also gives the measure of the SLE( 83 )-bubble that disconnects two points z, w ∈ H from ∞ which we will state as the following lemma.

= x + iy ∈ H is 1 4 ( y x 2 +y 2 ) 2 = 1 4 (Im 1 z )
Lemma 5 (see [START_REF] Beliaev | Some remarks on SLE bubbles and schramm's two-point observable[END_REF]). Let E(z, w) be the event that two points z, w ∈ H are disconnected from ∞ by a SLE( 83 ) curve from 0 to ϵ, then

µ ϵ [E(z, w)] = 1 4 Im ( 1 z )Im ( 1 w )G(σ(z, w))ϵ 2 + O(ϵ 3 ). (2.6)
where σ is defined as (1.4) and

G(t) = 1 -t 2 F 1 (1, 4 3 ; 5 3
; 1 -t).

(2.7)

Here 2 F 1 is the hypergeometric function.

Notice that when κ = 8 3 , it holds that 1 -8 κ = -2. so we have

µ bub SLE (0)[E(z, w)] = 1 4 Im ( 1 z )Im ( 1 w ) = 1 4 yv (x 2 + y 2 )(u 2 + v 2 )
G(σ(z, w)). (2.8) SLE( 83 )-bubble measure is closely related to Brownian bubble, in fact they only differ by a constant. Given a loop γ such that γ(0) = γ(t γ ) = 0 and γ(0, t γ ) ⊂ H, call the complement of the unbounded connected component of H γ the hull enclosed by γ. In the following we will let γ denote both the loop and hull enclosed by it. Lemma 6 (see [START_REF] Lawler | Conformal restriction: the chordal case[END_REF]). As measures on the hulls enclosed by loops, the following holds:

µ bub H (0) = 8 5 µ bub SLE( 83 ) (0). Proof. By the construction of the Brownian bubble measure at 0 (see the Chapter 5 of Lawler's book [START_REF] Gregory | Conformally invariant processes in the plane[END_REF]), it is the unique measure on loops(or hulls enclosed by loops) in H rooted at 0 such that the total mass that the sample intersect |z| = r is 1 r 2 for any r > 0. Since

{ {γ : γ ∩ {|z| = r} = ∅} : r > 0 }
is a algebra that generate the σ-algebra of the space of loops, we only need to show that the total mass of the SLE( 83 )-bubble sample intersecting |z| = r is 5 8r 2 . Define F ϵ (z) = z ϵ-z , the imagine of the circle |z| = r under F ϵ is a circle with center c 0 = -r 2 r 2 -ϵ 2 and radius ρ = ϵr r 2 -ϵ 2 . Define the conformal map ϕ ϵ (z) = z -c 0 + ρ 2 z-c 0 which maps H B(c 0 , ρ) onto H with the derivative at ∞ equaling to 1. By the conformal restriction property of SLE( 83 ), we have

µ ϵ [γ ∩ |z| = r = ∅] = µ ∞ [γ ∩ B(c 0 , ρ) = ∅] = ϕ ′ ϵ (0) 5 8 .
Therefore we can check that

µ bub SLE(κ) (0)[γ ∩ {|z| = r} ̸ = ∅] = lim ϵ→0 1 ϵ 2 (1 -ϕ ′ ϵ (0)) 5 8 = 5 8r 2 .

Proof of the main theorem

In this section we will give a detailed proof of our main theorem.

Proof of the nicer form

Given two points z 0 = x 0 + iy 0 and w 0 = u 0 + iv 0 ∈ H. By the symmetry property of the the Brownian loop measure, we may assume without loss of generality y 0 ≤ v 0 , u 0 ≥ x 0 . By (2.3),

µ loop H [E(z 0 , w 0 )] = 1 π ∫ H µ bub Hy (x + iy)[E(z 0 , w 0 )]dxdy = 1 π ∫ y 0 0 ∫ R µ bub Hy (x + iy)[E(z 0 , w 0 )]dxdy.
Here E(z 0 , w 0 ) denotes the event that the Brownian loop sample in H disconnects both z 0 and w 0 from the boundary of H.

By the translation invariance of the Brownian bubble measure, we have

µ bub Hy (x + iy)[E(z 0 , w 0 )] = µ bub H (0)[E(z 0 -z, w 0 -z)].
By Lemma 6 we have µ bub H (0) = 8 5 µ bub SLE(κ) (0). Therefore by (2.8),

µ loop H [E(z 0 , w 0 )] = 8 5π ∫ y 0 0 ∫ R µ bub SLE(κ) (0)[E(z 0 -x -iy, w 0 -x -iy)]dxdy = 8 5π ∫ y 0 0 ∫ R 1 4 Im ( 1 z 0 -x -iy ) Im ( 1 w 0 -x -iy ) G(σ(z 0 -x -iy, w 0 -x -iy))dxdy. (3.1)
So in order to prove the theorem, we only need to compute above integral. Define two functions as follows:

f (x, y) := (y 0 -y)(v 0 -y) [(x 0 -x) 2 + (y 0 -y) 2 ] × [(u 0 -x) 2 + (v 0 -y) 2 ] . (3.2) g(y) := (x 0 -u 0 ) 2 + (y 0 -v 0 ) 2 (x 0 -u 0 ) 2 + (y 0 + v 0 -2y) 2 2 F 1 (1, 4 3 ; 5 3
; 4(y 0 -y)(v 0 -y) (x 0 -u 0 ) 2 + (y 0 + v 0 -2y) 2 ). (3.3)
Lemma 7. Take the notations as above, for fixed y > 0,

∫ R f (x, y)dx = 2(y 0 -y) + v 0 -y 0 (x 0 -u 0 ) 2 + ( 2(y 0 -y) + v 0 -y 0 ) 2 π. (3.4)
Proof. For fixed y > 0, denote

a = y 0 -y, b = v 0 -y, c = u 0 -x 0 , d = v 0 -y 0 .
Then we have

f (x, y) = ab [(x 0 -x) 2 + a 2 ][(u 0 -x) 2 + b 2 ]
.

By standard calculus, ∫ R f (x, y)dx = ∫ R ab [(x 0 -x) 2 + a 2 ][(u 0 -x) 2 + b 2 ] dx =ab ∫ R 1 [x 2 + a 2 ][(x + c) 2 + b 2 ] dx = abπ ab ( a 4 -2a 2 (b 2 -c 2 ) + (b 2 + c 2 ) 2 ) ( b(b 2 + c 2 -a 2 ) arctan[ x a ] + a [ (a 2 + c 2 -b 2 ) arctan[ c + x b ] + bc log b 2 + (c + x) 2 a 2 + x 2 ] ) | ∞ -∞ =π b(b 2 + c 2 -a 2 ) + a(a 2 + c 2 -b 2 ) a 4 -2a 2 (b 2 -c 2 ) + (b 2 + c 2 ) 2 . Replace b by a + d, we get ∫ R f (x, y)dx = (2a + d)π c 2 + (2a + d) 2 ,
which is what we want.

By (3.1), we have

µ loop H [E(z 0 , w 0 )] = 8 5π ∫ y 0 0 ∫ R 1 4 f (x, y)(1 -g(y))dxdy = 8 5π ∫ y 0 0 π 4 2a + d c 2 + (2a + d) 2 [1 -g(y)]dy = 2 5 (A -B), (3.5)
where

A = A(z 0 , w 0 ) = ∫ y 0 0 2(y 0 -y) + v 0 -y 0 (x 0 -u 0 ) 2 + ( 2(y 0 -y) + v 0 -y 0 ) 2 dy, (3.6)
and

B = B(z 0 , w 0 ) = ∫ y 0 0 2(y 0 -y) + v 0 -y 0 (x 0 -u 0 ) 2 + ( 2(y 0 -y) + v 0 -y 0 ) 2 g(y)dy. (3.7) Lemma 8. A = 1 4 log 1 σ ,
where σ is defined as (1.4).

Proof. By (3.6) we have

A = ∫ y 0 0 2(y 0 -y) + d c 2 + (2(y 0 -y) + d) 2 dy = ∫ y 0 0 2y + d c 2 + (2y + d) 2 dy = 1 2
∫ 2y 0 +d c d/c y 1 + y 2 dy = 1 4 log c 2 + (2y 0 + d) 2 c 2 + d 2 .
In the second equation we used the change of variable y → y 0 -y and in the last equation we used the change of variable y → 2y+d c . Notice that

c 2 + (2y 0 + d) 2 c 2 + d 2 = (u 0 -x 0 ) 2 + (y 0 + v 0 ) 2 (u 0 -x 0 ) 2 + (v 0 -y 0 ) 2 = 1 σ . Lemma 9. B = 1 4 (1 -σ) 3 F 2 (1, 4 3 , 1; 5 3 , 2; 1 -σ),
where σ is defined as (1.4).

Proof. By (3.7) and the definition of g(y), we have

B = ∫ y 0 0 2(y 0 -y) + d c 2 + (2(y 0 -y) + d) 2 • c 2 + d 2 c 2 + (2(y 0 -y) + d) 2 • 2 F 1 (1, 4 3 ; 5 3 ; 4(y 0 -y)((y 0 -y) + d) c 2 + (2(y 0 -y) + d) 2 )dy = ∫ y 0 0 2y + d c 2 + (2y + d) 2 • c 2 + d 2 c 2 + (2y + d) 2 • 2 F 1 (1, 4 3 ; 5 3 ; 4y(y + d) c 2 + (2y + d) 2 )dy = ∫ 2y 0 +d c d c cy c 2 + c 2 y 2 • c 2 + d c 2 + c 2 y 2 • 2 F 1 (1, 4 3 ; 5 3
;

c 2 y 2 -d 2 c 2 + c 2 y 2 ) • c 2 dy = 1 2 c 2 + d 2 c 2 ∫ 2y 0 +d c d c y (1 + y 2 ) 2 • 2 F 1 (1,
(y 0 + d) c 2 + (2y 0 + d) 2 • 3 F 2 (1, 4 3 , 1; 5 3 , 2; 4y 0 (y 0 + d) c 2 + (2y 0 + d) 2 ) = 1 4 (1 -σ) 3 F 2 (1, 4 3 , 1; 5 3 , 2; 1 -σ).
Here the second equation used the change of variable y → y 0 -y, the third equation used the change of variable y → 2y+d c , the fifth equation used the change of variable c 2 y 2 -d 2 c 2 +c 2 y 2 → y and the sixth equation used the equation about hypergeometric functions below:

∫ x 0 2 F 1 (a, b; c, y)dy = x 3 F 2 (a, b, 1; c, 2, x).
Now by (3.5) and Lemma 8 and Lemma 9 we get (1.3).

Proof of Cardy's formula

In this section we will prove the equivalence of formula (1.1) and (1.3). First we will recall some identities for the hypergeometric functions which will be used in our proof. We will assume that our hypergeometric functions are all well defined. And they satisfies the following identities (see Chapter 8 of [START_REF] Beals | Special functions: a graduate text[END_REF]):

2 F 1 (a, b; c; x) = (1 -x) -b 2 F 1 (c -a, b; c; x x -1
).

(3.8)

2 F 1 (a, b; c; x) = Γ(c)Γ(c -a -b) Γ(c -a)Γ(c -b) 2 F 1 (a, b; a + b + 1 -c; 1 -x) + Γ(c)Γ(a + b -c) Γ(a)Γ(b) (1 -x) c-a-b 2 F 1 (c -a, c -b; c + 1 -a -b; 1 -x). (3.9) 2 F 1 (a, b; c; x) = (1 -x) c-a-b 2 F 1 (c -a, c -b; c; x). ( 3.10)
Notice that η = σ σ-1 and σ ∈ (0, 1). We define a function ϕ on [0, 1] as following. ). (3.11) To prove that (1.1) and (1.3) are equivalent, we only need to show that ϕ(t) ≡ 0. Notice that ϕ(0) = 2π √ 3 -3 F 2 (1, 4 3 , 1; 5 3 , 2; 1) = 0, it is left to show that ϕ ′ (t) ≡ 0. Take the following notations. ).

ϕ(t) = 2π √ 3 + t t -1 3 F 2 (1, 4
I(t) := t t -1 3 F 2 (1, 4
above lemma and the definition of the Brownian bubble measure and lemma 6, we can get

µ bub H (0)(E 1 (z, w)) = 1 10 [( y x 2 + y 2 ) 2 - y x 2 + y 2 v u 2 + v 2 G(σ(z, w))]. µ bub H (0)(E 2 (z, w)) = 1 10 [( v u 2 + v 2 ) 2 - y x 2 + y 2 v u 2 + v 2 G(σ(z, w))].
µ bub H (0)(E 3 (z, w)) = ∞. By relation (2.3) and calculating the integral on the upper half plane, we can see that the total mass of the Brownian loop measure on these three sets are infinite. In fact, we can see intuitively that these three cases all contain the loops with arbitrary small diameter, while the event E(z, w) in the main theorem exclude these small loops.

c 2 y 2 -d 2 c 2 + c 2 y 2

 22

Define f (x) = x 3 F 2 (1, 4 3 , 1; 5 3 , 2; x). It is easy to check that

; x).

So

; t).

(3.12)

The last equation follows from (3.8) by assigning

; 1 -t).

By letting a = 1 3 , b = 2 3 , c = 1 3 in (3.10), the following holds

Lastly we deal with the derivative of K(t) with respect to t. By letting

; t).

Consequently,

And

Combining (3.12),(3.13) and (3.14), we have

; t).

Lemma 10.

; t) = 0.

Proof. By definition we have

3 )Γ( 53 ) Γ( 43 )Γ(n + 5 3 )

Therefore

3 )Γ( 53 ) Γ( 43 )Γ(n + 2 3 )

3 )Γ( 53 ) Γ( 13 )Γ(n + 5 3 )

By using the relation Γ(x + 1) = xΓ(x), we can see that the coefficient of t n in the sum is Γ(n + 1 3 )Γ( 53 ) Γ( 43 )Γ(n + 2 3 )

3 )Γ( 53 ) Γ( 43 )Γ(n + 5 3 )

From Lemma 10, we have ϕ ′ (t) ≡ 0, and therefore ϕ ≡ ϕ(0) = 0. This completes the proof of the equivalence between (1.1) and (1.3).

The other cases

Given z, w ∈ H, and γ the sample of the Brownian loop in the upper half plane. According to the property of Brownian path, almost surely, z, w ̸ ∈ γ. So except the case that γ disconnects both z and w from the boundary, there are three other cases:

(1) γ disconnects z from the boundary but does not disconnect w from the boundary;

(2) γ disconnects w from the boundary but does not disconnect z from the boundary;

(3) γ does neither disconnects z from the boundary nor disconnects w from the boundary.

We will show that the total measure of above three cases are infinite. In fact, using the same method as [START_REF] Beliaev | Some remarks on SLE bubbles and schramm's two-point observable[END_REF], we can show the following lemma.

Lemma 11. Suppose that γ is the sample of the SLE( 83 ) from 0 to ϵ and denote above three cases by E 1 (z, w), E 2 (z, w) and E 3 (z, w) respectively. Then

3)

The proof of this lemma is the same as in [START_REF] Beliaev | Some remarks on SLE bubbles and schramm's two-point observable[END_REF]. We only need to prove that for SLE( 83 ) γ from 0 to ∞, the following holds. P[γ passes the left of z and the right of w] where G(σ) is the same as (2.7). Then using the conformal map F ϵ (z) = ϵz 1+z to convert the SLE( 83 ) from 0 to ∞ into the SLE( 83 ) from 0 to ϵ. Combing