10.6092/issn.1972-5787/9730
Large-integer arithmetic algorithms are used in contexts where both their performance and their correctness are critical, such as cryptographic software. The fastest algorithms are complex enough that formally verifying them is desirable but challenging. We have formally verified a comprehensive arbitrary-precision integer arithmetic library that implements many state-of-the-art algorithms from the GMP library. The algorithms we have verified include addition, subtraction, Toom-Cook multiplication, division and square root. We use the Why3 platform to perform the proof semi-automatically. We obtain an efficient and formally verified C library of low-level functions on arbitrary-precision natural integers. This verification covers the functional correctness of the algorithms, as well as the memory safety and absence of arithmetic overflows of their C implementation. Using detailed pseudocode, we present the algorithms that we verified and outline their proofs.
application/pdf
1972-5787
9730
1
en
Raphael Rieu-Helft
Journal of Formalized Reasoning
arbitrary-precision integer arithmetic
A Why3 proof of GMP algorithms
Text.Serial.Journal
https://jfr.unibo.it/
12
Articles
2020-02-13
2019-07-24
2019-02-19
2020-02-19
53 - 97
JFR
Raphael Rieu-Helft
2019
10.6092/issn.1972-5787/9730
53
TrustInSoft/Inria
1972-5787
1
arbitrary-precision integer arithmetic
en
97
Jean Karim Zinzindohoué, Karthikeyan Bhargavan, Jonathan Protzenko, and Benjamin Beurdouche. HACL*: A verified modern cryptographic library. Cryptology ePrint Archive, Report 2017/536, 2017. https://eprint.iacr.org/2017/536
A Why3 proof of GMP algorithms
https://jfr.unibo.it/article/download/9730/10489
12
A Why3 proof of GMP algorithms
Raphaël Rieu-Helft
TrustInSoft
Inria


Large-integer arithmetic algorithms are used in contexts where both their performance and their
correctness are critical, such as cryptographic software. The fastest algorithms are complex enough
that formally verifying them is desirable but challenging. We have formally verified a compre-
hensive arbitrary-precision integer arithmetic library that implements many state-of-the-art al-
gorithms from the GMP library. The algorithms we have verified include addition, subtraction,
Toom-Cook multiplication, division and square root. We use the Why3 platform to perform the
proof semi-automatically. We obtain an efficient and formally verified C library of low-level func-
tions on arbitrary-precision natural integers. This verification covers the functional correctness of
the algorithms, as well as the memory safety and absence of arithmetic overflows of their C im-
plementation. Using detailed pseudocode, we present the algorithms that we verified and outline
their proofs.




1.   INTRODUCTION
The GNU Multi-Precision library, or GMP,1 is a very widely used library for arbitrary-
precision arithmetic. It provides state-of-the-art algorithms for basic arithmetic operations
on integer, rational and floating-point numbers as well as number-theoretic primitives. It is
used in safety-critical contexts such as cryptography and security of Internet applications.
For performance reasons, many of these algorithms are very intricate, which makes the
presence of correctness bugs a real possibility. Moreover, even proving the absence of
crashes and memory issues sometimes relies on functional correctness properties, such as
the fact that a carry propagation cannot go past the bounds of an array.
   The library is extensively tested, but some parts of the code are visited only in very un-
likely cases, such as a 1/264 probability (e.g. a case in long division that occurs when two
machine words have the same value, assuming uniformly distributed inputs). This makes
random testing very challenging. Correctness bugs occuring with very low probability
have in fact been found in the past.2 In order to ensure the absence of correctness bugs,
formal verification for all possible inputs is therefore desirable.
   GMP has multiple layers, which handle different types of numbers (natural, relative,
rational, floating-point) at different levels of abstraction. This work focuses only on the
lowest-level layer of GMP, called mpn, which provides arithmetic primitives for arbitrary-
sized natural integers. We did not verify the entire mpn layer. Indeed, for each operation,
mpn implements many different algorithms, each most suitable for different input sizes.
The fragment of mpn that we have verified includes at least one algorithm for each of addi-
tion, subtraction, multiplication, division and square root (fast modular exponentiation is
a work-in-progress). This paper mostly focuses on the functional correctness of the algo-
rithms from a mathematical point of view. Indeed, most implementation-level concerns are

1 https://gmplib.org/
2 Look   for ’division’ at https://gmplib.org/gmp5.0.html.

                                                             Journal of Formalized Reasoning Vol.12, No.1, 2019, Pages 53–97.
54       ·      Raphaël Rieu-Helft

automatically proved with very little user input on top of the program code. This leaves the
proof of higher-level mathematical facts as the only non-trivial part of the proof. Further-
more, the Why3 mechanized proofs of these mathematical facts are very similar to paper
proofs. Therefore, the proofs we present in this work deal with faithful transcriptions of the
algorithms in detailed pseudocode. We give the specifications and invariants, and provide
step-by-step explanations of the algorithms and proofs that the implementations indeed
match the specifications.
   Our approach is to implement the GMP algorithms (originally implemented in C and
handwritten assembly) in WhyML, the high-level programming language supported by the
Why3 verification environment. We give the algorithms a formal specification, and verify
that they fulfill these specifications with the help of Why3 and automated theorem provers.
Our functions are implemented on top of a Why3 model of the C language that accounts
for memory safety and lack of runtime errors. Thus, they can be implemented using the
low-level expressivity of C that we need to mirror GMP’s own implementation. The model
also permits a direct compilation from WhyML to C. Thus, we obtain a formally verified
C library that is both compatible with GMP and almost as efficient for inputs smaller
than about 1,000 bits (100,000 for multiplication). The full development is available at
http://toccata.lri.fr/gallery/multiprecision.en.html. We have found no bugs in
GMP while performing this verification, but the developers modified a section of a Toom-
Cook algorithm (Toom-2) after we pointed out that its correction relied on much more
intricate reasonings than what seemed needed. Section 2 provides some context on the
process of Why3 proofs. As our Why3 model of C has already been presented [24], the
rest of this paper focuses on the higher-level, mathematical aspects of the verification.
   There are several challenges to overcome when carrying out such a verification work.
The first is to understand why the algorithms are in fact correct. This is easy for the algo-
rithms that we discuss in Section 3 (comparison, addition, subtraction, schoolbook multi-
plication), but represents a significant part of the work for most algorithms. Some of the
optimizations layered over the algorithms obfuscate them somewhat, such as the division
algorithm in Section 4. Moreover, seemingly innocuous statements sometimes require a
complex proof in order to show that, for example, an operation or a carry propagation does
not overflow. Section 5, which focuses on the Toom-Cook divide-and-conquer multipli-
cation, provides good examples of this. Textbooks such as Brent and Zimmermann’s [6]
explain the broad strokes of most of the algorithms, but do not mention any of these im-
plementation tricks in their pseudocode, not to mention that some algorithms are simply
more recent than the textbooks. Many of these algorithms do not seem to have been pub-
lished by their original authors. The second challenge is to structure the proof in such a
way that Why3 and the automated provers can replay it. Due to the nature of the weakest-
precondition calculus used by Why3, all loops must be annotated with inductive invariants,
and all functions must carry extensive specifications. Furthermore, the automated provers
we use do not handle logical goals that involve non-linear arithmetic and numbers with
unknown size very well. In practice, the user sometimes has to provide extensive proof
hints in order to prove seemingly simple facts.
   We discuss GMP’s square root algorithm in Section 6. The divide-and-conquer part of it
has already been verified in Coq by Bertot et al. [3], and our proof is largely adapted from
theirs. Related work, covered in Section 7, otherwise mostly comprises formal proofs of
single state-of-the-art algorithms, and verified libraries of simpler, less efficient algorithms.
To the best of our knowledge, this work is the first formal verification of a comprehensive
Journal of Formalized Reasoning Vol.12, No.1, 2019
                                                     A Why3 proof of GMP algorithms                 ·      55

state-of-the-art arbitrary-precision integer library.

2.     WHY3 AND GMP
The goal of this section is to provide some context on the process of verifying GMP algo-
rithms. We will first go over GMP’s number representation, and then explain the process
of Why3 proofs, using a very simple GMP routine as an example.

2.1     GMP number representation
In GMP, natural integers are represented as arrays of unsigned integers called limbs. We
set a radix β = 264 (also calledP   radix in the formal development). Any natural number
                                      n−1
N has a unique decomposition k=0 a[k]β k in base β such that a[n − 1] 6= 0, and is
represented as the buffer a[0]a[1] . . . a[n − 1] (with the least significant limb first).
   For efficiency, there is no memory management in GMP’s low-level functions, so the
caller code has to keep track of number sizes. Instead, operands are specified by a pointer
to their least significant limb and a limb count of type int32.
 type limb = uint64


 type t = ptr limb

    A pointer is considered valid over a size s if it is not null and points to a zone of size at
least s.
    Let us now establish the link between mathematical integers and arrays of machine in-
tegers. If a pointer a is valid over a size n, we denote value(a, n) = a[0] . . . a[n − 1] =
Pn−1          k
    k=0 a[k]β . We also denote value(incr(a, k), p) as a[k]...a[k + p − 1].
    We use a few lemmas to express what happens to the value of a big integer when part of
it is modified. In particular, loops tend to change only one end of a big integer (usually by
increasing its length), and being able to separate what changed from what did not is crucial
to prove that loop invariants are preserved.
     L EMMA 1. Let p a pointer valid over a sufficiently large length.
∀k, p[0] . . . p[n − 1] = p[0] . . . p[k − 1] + β k p[k] . . . p[n − 1]           [value_sub_concat]
p[0] . . . p[n] = p[0] . . . p[n − 1] + β n p[n]                                     [value_sub_tail]
p[0] . . . p[n − 1] = p[0] + βp[1] . . . p[n − 1]                                    [value_sub_head]
∀k, v, p[0] . . . p[k − 1]vp[k + 1] . . . p[n − 1]
= p[0] . . . p[k] . . . p[n − 1] + β k (v − p[k])                                 [value_sub_update]
   We also need some bounds for the value of an integer of size n. These follow easily
from the fact that the values a[i] lie between 0 and β.
     L EMMA 2. Let p a pointer valid over the size n.
       0 ≤ p[0] . . . p[n − 1]                                     [value_sub_lower_bound]
       p[n − 1]β n−1 ≤ p[0] . . . p[n − 1]                [value_sub_lower_bound_tight]
                                 n
       p[0] . . . p[n − 1] < β                                     [value_sub_upper_bound]
       p[0] . . . p[n − 1] < (p[n − 1] + 1)β n−1          [value_sub_upper_bound_tight]

                                                             Journal of Formalized Reasoning Vol.12, No.1, 2019
56          ·    Raphaël Rieu-Helft

2.2     From program to logical goals
Let us now go over the process of Why3 proofs with the example of GMP’s copy routine,
shown in Alg. 1. It takes two numbers r and x valid over a size n and copies the contents
of x into r, starting with the least significant limbs. Note that x and r are allowed to point
to zones longer than n, or even to point to the middle of a number. The only requirement
is that there are at least n valid limbs to the right of each pointer.

Algorithm 1 Copy of a long integer
  function COPYI(r, x, n)
     i←0
     while i < n do
         r[i] ← x[i]
         i←i+1

   The Why3 deductive program verification platform provides a language for specification
and programming called WhyML, the programming part of which is a dialect of ML. The
first step of a Why3 proof is to translate the program that we want to prove in WhyML. We
also use WhyML to write a contract for the program. For the COPYI function, the (slightly
simplified) Why3 implementation is shown in Figure 1.

  let wmpn_copyi (r x: ptr uint64) (n: int32) : unit
  =
      let ref i = 0 in
      while (Int32.(<) i n) do
        r[i]   ← x[i];
        i    ← i+1;
      done



                               Fig. 1.   WhyML transcription of the COPYI function.

   WhyML functions are annotated with contracts based on Hoare logic [17, 15]. Given an
expression e, a precondition P and a postcondition Q, where P and Q are assertions in first-
order logic, one can write the Hoare triple {P }e{Q}. Its intuitive reading is: Whenever
the state before the execution of e is such that P holds, then the computation terminates,
there is no runtime error, and the final state satisfies Q. Hoare triples can be derived from
a set of inference rules based on program syntax. A soundness theorem implies that any
derivable triple is correct. Some examples are shown in Figure 2.

                                                                     {P }s1 {Q}     {Q}s2 {R}
                     {P }skip{P }        {P [x ← e]}x ← e{P }             {P }s1 ; s2 {R}

                                     Fig. 2.   Some Hoare triple inference rules.

   The contract of a function forms a Hoare triple with the function body. The precondition
P is the conjunction of the requires clauses, and the postcondition Q is the conjunction
of the ensures clauses. For example, the contract of COPYI is as follows:


Journal of Formalized Reasoning Vol.12, No.1, 2019
                                                     A Why3 proof of GMP algorithms                     ·   57

 let wmpn_copyi (r x: ptr uint64) (n: int32) : unit
   requires { valid x n∧ valid r n }
   ensures               ≤ i < n → r[i] = x[i] }
             { forall i. 0
   ensures { forall i. i < 0 ∨ n ≤ i → r[i] = old r[i] }

   Using this contract, Why3 produces a logical goal that implies that the program satisfies
this specification. The aim is to find out whether the Hoare triple defined by the contract
is valid. Using the inference rules naively would be far too tedious. Indeed, the rule
for sequences of instructions (see Figure 2) requires an intermediate assertion between
each pair of statements. The code would need to be very heavily annotated for the naive
approach to work.
   Instead, Why3 uses Dijkstra’s weakest-precondition calculus [10]. We define the pred-
icate transformer WP(., .). Where e is an expression and Q a postcondition, WP(e, Q)
computes the weakest precondition P such that {P }e{Q} holds. This is a way to auto-
mate the search for intermediate assertions in the derivations. A subset of the computation
rules for WP(., .) can be found in Figure 3.

                                      WP(skip, Q) = Q
                             WP(e1 ; e2 , Q) =WP(e1 , WP(e2 , Q))
                             WP(assert R, Q) = R ∧ (R → Q)

                                         Fig. 3.   Some WP rules.

   The soundness theorem of WP states that for any e and Q, the triple {WP(e, Q)}e{Q}
is valid. Therefore, to show that a contract {P }e{Q} is valid, it suffices to prove that
P → WP(e, Q). Why3 computes WP(e, Q) and outputs P → WP(e, Q) as a logical
goal for the user to prove.
   One major difficulty remains. The previous figures do not show the rules for deriving
Hoare triples or computing WP for loops. This problem is undecidable in general, as
the loop would need to be unrolled until termination. Therefore, WhyML loops must be
annotated with an inductive loop invariant. The corresponding WP computation makes
sure that this invariant is valid at the start of the loop, maintained through one iteration of
the loop, and that at the end of the loop, it implies the required postcondition. A simplified
version of this rule is shown in Figure 4. It assumes that the only modification of the
memory state in the loop occurs on a variable v, and does not handle termination.

          WP(while t invariant J do e done, Q) ≡
              J∧                                                           [invariant initialisation]
              ∀v.
                (J ∧ t → WP(e, J))∧                                        [invariant preservation]
                (J ∧ ¬t → Q)                                        [loop exit and postcondition]



                               Fig. 4.     Simplified WP rule for while.

   The only problem that remains is proving termination, which would also be undecidable
in general. Why3 requires that proofs and recursive functions be annotated with a variant.


                                                              Journal of Formalized Reasoning Vol.12, No.1, 2019
58          ·    Raphaël Rieu-Helft

This variant is a program expression that is decreasing at each loop iteration (or function
call) on a well-founded order. In the end, the annotated version of COPYI is as follows:
  let wmpn_copyi (r x: ptr uint64) (n: int32) : unit
      requires { valid x n∧ valid r n }
      ensures               ≤ i < n → r[i] = x[i] }
                 { forall i. 0
      ensures { forall i. i < 0 ∨ n ≤ i → r[i] = old r[i] }
  = let ref i = 0 in
      while (Int32.(<) i n) do
        variant { n - i }
        invariant { forall j. 0       ≤ j < i → r[j] = x[j] }
        invariant { forall j. j       < 0 ∨ i ≤ j → r[j] = old r[j] }
        r[i]   ← x[i];
        i    ← i+1;
      done

   Note that it is not necessary to specify that the memory block pointed by x was not
modified, which would happen if x and r were aliased. Indeed, Why3’s region-based type
system imposes strong constraints upon WhyML programs, in such a way that all aliases
are known statically and encoded in the types [16]. We used this fact to design our C
memory model in such a way that pointers used as function parameters are forced to be
separated by typing (unless they are read-only). Moreover, typing enforces that separated
pointers cannot reach each other through pointer arithmetic. This makes proofs much
simpler, both in terms of specification length (many more annotations would otherwise be
necessary) and in terms of size of the proof context.
2.3     Proving the goals
Once Why3 has generated the logical goals that imply the program meets its specification,
they should be proved by the user. Why3 offers several ways for them to do so. The
most important one is to dispatch the proof obligations to a variety of automated theorem
provers. Why3 has drivers for various provers, such as SMT solvers Alt-Ergo, Z3, and
CVC4, as well as superposition-based provers such as E, SPASS and Vampire. The default
way of proving a goal is to send it to one or more SMT solvers. If at least one of them
returns that the goal is correct, Why3 considers it proved. Why3 offers task transformations
that can be used, for example, to split a very large goal such as the weakest-precondition
of a large program into many tractable ones.
   While SMT solvers are very powerful tools in many contexts, sometimes proof obli-
gations fall outside their scope and they do not manage to prove them in any reasonable
timeframe. For example, the theory of nonlinear integer arithmetic is undecidable. While
solvers do manage to solve some nonlinear arithmetic goals using various proof search
techniques, many goals that naturally arise during the verification of GMP algorithms are
not solved by any of the provers that we tried. Furthermore, the fact that GMP integers
are arbitrary-sized and the large size of the proof contexts make these goals even harder to
solve using non-complete heuristics.
   When simply sending the goal to automated solvers is not enough to prove it in a rea-
sonable time, Why3 offers several tools for the user to provide extra guidance. First, the
user can write a proof in an interactive theorem prover, such as Coq, Isabelle or PVS. This
approach is generally quite tedious, and we used it only for a handful of goals in this work.
Second, the Why3 user can declare and prove lemmas. These lemmas are passed to the
solvers alongside the rest of the proof context. The user may also call them like functions
Journal of Formalized Reasoning Vol.12, No.1, 2019
                                                 A Why3 proof of GMP algorithms                   ·      59

in the code, in order to force their application on specific arguments. Figure 5 shows an
example of this. The first lemma is declared and proved trivially, and the second lemma
is proved using the first lemma twice. Note that arbitrary control structures can be used in
the proof of lemmas, as they are simply ghost WhyML functions. For example, the third
lemma is proved by induction over n, using an if statement to split between the base case
and the inductive case. The right column shows the interpretation of these functions as
lemmas in Why3’s logical context.

 let lemma prod_compat_r (a b c:int)
   requires { 0 ≤ a ≤ b }
   requires { 0 ≤ c }                               ∀a, b, c ∈ Z, 0 ≤ a ≤ b ∧ 0 ≤ c =⇒ ca ≤ cb
   ensures { c * a ≤ c * b }
 = ()

 let lemma prod_compat_lr (a b c d:int)
   requires { 0≤ a ≤ b }
   requires { 0≤ c ≤ d }
   ensures { a * c ≤ b * d }                 ∀a, b, c, d ∈ Z, 0 ≤ a ≤ b ∧ 0 ≤ c ≤ d =⇒ ac ≤ bd
 = prod_compat_r c d a;    (* ac   ≤ ad *)
   prod_compat_r a b d     (* da   ≤ db *)

 let rec lemma pow2_gt (n:int)
   requires { 0    ≤ n }
   ensures { n    ≤ power 2 n }                                       ∀n ∈ Z, 0 ≤ n =⇒ n ≤ 2n
   variant { n }
 = if n   > 0 then pow2_gt (n-1)


                            Fig. 5. Declaring, proving and using lemmas.

   A third way to prove difficult goals is to write assertions inside the code. The assertions
are logical formulas that are interpreted as proof cuts, using the WP rules for sequences
and assert in Figure 3. Why3 checks that the assertion is correct using the same methods
as for the other goals, and then inserts it in the logical context. Figure 6 features a Why3
lemma and its proof, which consists in a large assertion. In our case, the assertion is a full-
fledged proof, written using the logical connectives by and so [7]. The proof is essentially
a sequence of implications, which Why3 checks using automated provers. The end result
is essentially what we would expect from a paper proof. Much like in a paper proof, when
developing such a Why3 assertion, we simply add extra steps when the automated provers
do not manage to check an implication. Many of the mathematical proofs that are given in
this paper were proved using this method in our Why3 development.
   Finally, another way for the Why3 user to prove difficult goals is computational reflec-
tion. Why3 allows a user to define and verify dedicated decision procedures as WhyML
programs, and to execute these procedures in order to prove the goals [20]. We used this
approach in our development to prove goals that previously required hundreds of lines of
handwritten assertions.
   Using these various tools, it is usually not too difficult to go from a detailed paper proof
to a Why3 mechanized proof. Aside from a few goals related to the absence of overflow in
some operations (which we do discuss), most implementation-level concerns are automat-
ically discharged by the SMT solvers with next to no work on our part. The proof effort

                                                           Journal of Formalized Reasoning Vol.12, No.1, 2019
60       ·      Raphaël Rieu-Helft


  let lemma fact_div (x y z:int)
      requires { y   > 0 }
      ensures   { div (x + y * z) y = (div x y) + z }
  =
      assert { div (x + y * z) y = (div x y) + z
               by x + y * z = y * (div (x + y * z) y) + mod (x + y * z) y
               so mod (x + y * z) y = mod (y * z + x) y = mod x y
                so x + y * z = y * (div (x + y * z) y) + mod x y
                so x = y * div x y + mod x y
                so x + y * z =      y * div x y + mod x y + y * z
                so y * (div (x      + y * z) y) + mod x y
                   = y * div x      y + mod x y + y * z
                so y * (div (x      + y * z) y)
                     = y * div x y + y * z
                     = y * ((div x y) + z)
                so div (x + y * z) y = div x y + z }



                                   Fig. 6.   A Why3 proof using a large assertion.

was overwhelmingly concentrated on the mathematical correctness of the algorithms. In
the rest of this paper, we will therefore focus on proving the functional correctness of
the various algorithms on paper, while abstracting some language-specific concerns away.
Furthermore, the Why3 proofs of mathematical correctness are extremely similar to paper
proofs, as exemplified by the proof in Figure 6. Therefore, we will focus on proving the
correctness of the algorithms on paper, rather than showing the Why3 translation of these
paper proofs.

3.     BASIC ALGORITHMS
3.1     Comparisons
The mpn layer of GMP exposes a single comparison function, which compares two integers
of same length (Alg. 2). The algorithm is very straightforward: it simply iterates both
operands until it finds a difference, starting at the most significant limb.

Algorithm 2 Comparison of two integers of identical length.
  function CMP(x, y, n)
Require: valid(x, n), valid(y, n)
Ensure: result > 0 ⇔ value(x, n) > value(y, n)
Ensure: result = 0 ⇔ value(x, n) = value(y, n)
Ensure: result < 0 ⇔ value(x, n) < value(y, n)
     for i = n − 1 downto 0 do
         if x[i] 6= y[i] then
             if x[i] > y[i] then return 1
             else return −1
     return 0

   Since the function involves a loop, we must provide a loop invariant. Here, the loop
invariant is that both source operands are identical from offsets i + 1 to n.

Journal of Formalized Reasoning Vol.12, No.1, 2019
                                                A Why3 proof of GMP algorithms                  ·      61

  An additional lemma is required to prove the invariant and complete the proof, it simply
says that two big integers with equal limbs at all offsets are equal:
  L EMMA value_sub_frame.
Let a[0], . . . , a[n − 1], b[0], . . . , b[m − 1] such that for all u ≤ k ≤ v, a[k] = b[k]. Then
a[u] . . . a[v] = b[u] . . . b[v].
  The lemma is proved by a straightforward induction, which translates well into a Why3
lemma function as the recursive call takes care of the inductive case:
 let rec lemma value_sub_frame (a b:map int limb) (u v:int)
      requires { forall i. u   ≤ i < v → a[i] = b[i] }
      variant
           { m - n }
   ensures { value_sub a u v = value_sub b u v }
 = if u < v then value_sub_frame a b (u+1) v else ()

   This lemma shows that the numbers are equal if no difference was found by the end of
the loop.
3.2     Addition, Subtraction
We use the schoolbook algorithms for the addition and subtraction of big integers, repre-
sented by their decomposition in base β.
   We first need to give the specifications of basic operations on limbs. The following three
primitives can be used to add limbs.
   The first primitive + is the defensive addition: it requires that the sum of the two inputs
does not overflow.
 val (+) (a b:limb) : limb
      requires { "expl:integer overflow" Limb.min ≤ to_int a + to_int b     ≤ Limb.max }
      ensures { to_int result = to_int a + to_int b }

   The second primitive, add_mod, has the semantics of the + operator on unsigned integers
in C: if there is an overflow, the result wraps around. When compiling WhyML programs
to C, both (+) and add_mod are translated to the C operator +. The former has a stronger
postcondition, so using it simplifies the proofs when its (also stronger) precondition is met.
The latter captures the full semantics of the addition, so it can be used in the remaining
cases.
 val add_mod (x y:limb) : limb
   ensures { to_int result = mod (to_int x + to_int y) radix }

  Finally, the third primitive accepts a carry to be added to the other two operands, and
outputs both the carry and the (potentially wrapped-around) result of the addition.
 val add_with_carry (x y:limb) (c:limb) : (limb,limb)
   requires { 0 ≤ to_int c ≤ 1 }
      returns { (r,d) →
        to_int r + radix * to_int d =
        to_int x + to_int y + to_int c
        ∧ 0 ≤ to_int d ≤ 1 }

   Similar primitives are used for limb subtraction.
   Our library implements many variants of addition and subtraction (depending on whether
the operation is done in place, the operation may overflow, the operands are known to be
of same length, etc.) As an example, let us examine the general-case addition (Alg. 3).
                                                         Journal of Formalized Reasoning Vol.12, No.1, 2019
62       ·      Raphaël Rieu-Helft

Algorithm 3 Addition of two integers
Require: 0 ≤ n ≤ m
Require: valid(a, m), valid(b, n), valid(r, m)
Ensure: value(r, m) + β m · result = value(a, m) + value(b, n)
Ensure: 0 ≤ result ≤ 1
  function ADD(r, a, m, b, n)
     i←0
     while i < n do                                                      . Add b to a.
          x ← a[i]
          y ← b[i]
          (z, c) ← ADD _ WITH _ CARRY(x, y, c)
          r[i] ← z
          i←i+1
     if c 6= 0 then                . Keep copying a into r while propagating the carry.
          while i < m do
              x ← a[i]
              z ← ADD _ MOD(x, 1)                                             . c = 1.
              r[i] ← z
              i←i+1
              if z 6= 0 then                    . No overflow: there is no more carry.
                   c←0
                   break
     while i < m do                                    . No more carry: copy a into r.
          r[i] ← a[i]
          i←i+1
     return c

   The algorithm is schoolbook addition with a few optimisations. There are three main
steps.
   The first step is to add together the two operands over the length of the shorter one. This
corresponds to the first while loop. Its loop invariants are as follows:
(1) value(r, i) + cβ i = value(a, i) + value(b, i)
(2) 0 ≤ i ≤ n.
   At the end of the first loop, we have value(r, n) + cβ n = value(a, n) + value(b, n) and
i = n. What remains to be done is to copy the last m − n limbs of a into r and propagate
the carry.
   The second loop copies a into r while propagating the carry. It is skipped if c = 0. Its
loop invariants are:
(1) value(r, i) + cβ i = value(a, i) + value(b, n),
(2) n ≤ i ≤ m,
(3) i = m ∨ c = 1.
We break out of the loop whenever the carry c becomes 0 (or if we have finished copying a
into r, in which case i = m). Note that we are only adding the carry to a limb x of a. There
is an overflow if and only if x = β − 1. There is no need to use ADD _ WITH _ CARRY. It
is more efficient to instead use ADD _ MOD and check if the result is 0, in which case there
Journal of Formalized Reasoning Vol.12, No.1, 2019
                                              A Why3 proof of GMP algorithms                  ·      63

was an overflow. This is relatively unlikely (probability 1/β if the operands are randomly
drawn from a uniform distribution), so the loop typically only runs for zero or one iteration.
  Finally, in the third loop (which is skipped if we already have i = m), we only have to
copy the last limbs of a into r. The loop invariants are:
(1) value(r, i) + cβ i = value(a, i) + value(b, n),
(2) n ≤ i ≤ m,
(3) i = m ∨ c = 0.
  At the end of the loop, we can return c and easily see that the postcondition is verified.
3.3   Schoolbook multiplication
We have implemented several algorithms for integer multiplication. For smaller integers
(fewer than 30 limbs), the fastest is the schoolbook one, which has complexity O(n2 ). For
larger operands, Toom-Cook multiplication is used (Section 5), as it has better asymptotic
complexity.

Algorithm 4 Multiply-and-add
Require: valid(r, m), valid(a, m)
Ensure: value(r, m) + β m · result = value(old r, m) + value(a, m) × y
Ensure: ∀j.j < 0 ∨ m ≤ j ⇒ r[j] = (old r)[j]
  function ADDMUL _1(r, a, m, y)
     c, i = 0
     while i < m do
          x ← a[i]
          z ← r[i]
          (l, h) ← MUL _ DOUBLE(x, y)                            .l+β·h=x·y
          v, c0 ← ADD 3(z, l, c)                           . v + β · c0 = z + l + c
          r[i] ← v
          c ← c0 + h
          i←i+1
     return c

   Let us first consider the auxiliary function ADDMUL _1 (Alg. 4). It multiplies a big
integer a by a limb y and adds the result to r, without modifying the contents of r outside
the area of the addition.
   The loop invariants of ADDMUL _1 are:
(1) 0 ≤ i ≤ m,
(2) value(r, i) + cβ i = value(old r, i) + value(a, i) · y,
(3) ∀j.j < 0 ∨ m ≤ j ⇒ r[j] = (old r)[j].
  Let us now consider the main function that implements schoolbook multiplication (Alg. 5).
At each loop iteration, one limb of the second operand is multiplied by the entire first
operand, and the product is added to the result (shifted appropriately).
  The loop invariants are as follows:
(1) 1 ≤ i ≤ n,
(2) value(r, m + i) = value(a, m) · value(b, i),
                                                       Journal of Formalized Reasoning Vol.12, No.1, 2019
64       ·         Raphaël Rieu-Helft

Algorithm 5 Schoolbook multiplication of two integers
Require: 0 < n ≤ m
Require: valid(a, m), valid(b, n), valid(r, m + n)
Ensure: value(r, m + n) = value(a, m) · value(b, n)
  function MUL _ BASECASE(r, a, b, m, n)
     y ← b[0]
     r[m] ← MUL _1(r, a, y, m)                 . value(r, m + 1) = value(a, m) · b[0].
     p←r+1
     i←1
     while i < n do
         y ← b[i]
         p[m] ← ADDMUL _1(p, a, m, y)                                   . See Alg. 4.
         i←i+1
         p←p+1

(3) p = r + i.
   It is easy to see that the postcondition follows from the invariants. The fact that the
invariants are maintained follows from the specification of the auxiliary function ADD -
MUL _1.
   Indeed, if we pose r0 the state of r at the beginning of the loop, we have at the end of the
loop:
value(r, m + i + 1)

= value(r, i) + β i · value(r + i, m + 1)                                        decomposition
               0          i
= value(r , i) + β · value(r + i, m + 1)                                     no writes in r(0, i)
               0          i               0
= value(r , i) + β · (value(r + i, m) + value(a, m) · y) postcondition of ADDMUL _1
= value(r0 , m + i) + β i · value(a, m) · y                                      recomposition
                                              i
= value(a, m) · value(b, i) + β · value(a, m) · y                                 loop invariant
                                                  i
= value(a, m) · (value(b, i) + β · y)
= value(a, m) · value(b, i + 1).                                                 recomposition

4.    DIVISION
Long division consists in computing the quotient and remainder of the division of big
integers of arbitrary sizes. It is a significantly more complex problem than long addition
and multiplication. While GMP’s algorithm is a variation on the schoolbook algorithm, it
is thoroughly optimized to the point of making it hard to understand and prove. The GMP
general-case long division function is about 30-line long,3 and the code extracted from our
implementation has about the same length. However, our proof for it is about 2000-line
long.
   Let us first review a more naïve algorithm: Knuth’s Algorithm D [19, p.257] (see
also [28]), shown in Alg. 6. We did not use this algorithm in our development, but it
is simple enough to explain the core ideas more easily.

3 mpn/generic/sbpi1_div_qr.c       in GMP 6.1.2

Journal of Formalized Reasoning Vol.12, No.1, 2019
                                              A Why3 proof of GMP algorithms                  ·      65

Algorithm 6 Knuth’s Algorithm D
Require: m ≥ n > 0, valid(a, m), valid(d, n), valid(q, m − n), valid(r, n)
Require: d[n − 1] ≥ β/2
Require: value(a + m − n, n) < value(d, n)           . Otherwise an extra quotient limb is
    needed.
Ensure: value(a, m) = value(d, n) × value(q, m − n) + value(r, n)
Ensure: value(r, n) < value(d, n)
 1: function A LGORITHM D(q, r, a, d, m, n)
 2:    for j = m − n − 1 downto 0 do
 3:        q̂ ← DIV _2 BY 1(a[j + n − 1], a[j + n], d[n − 1]) . Candidate quotient limb.
 4:        r̂ ← βa[j + n] + a[j + n − 1] − d[n − 1] × q̂          . Candidate remainder.
         adjust:
 5:         if q̂ ≥ β or q̂ × d[n − 2] > β × r̂ + a[j + n − 2] then
 6:               q̂ ← q̂ − 1                                . Quotient is too large; adjust.
 7:               r̂ ← r̂ + d[n − 1]
 8:               if r̂ < β then goto adjust                       . Happens at most once.
 9:         b ← SUBMUL _ IN _ PLACE(a + j, d, q̂, n)               . Subtract d × q̂ from a.
10:         q[j] ← q̂
11:         if b > 0 then                  . There was a borrow, the quotient was too large.
12:               q[j] ← q[j] − 1
13:               c ← ADD _ N _ IN _ PLACE(a + j, d, n)
14:               a[j + n] ← a[j + n] + c                             . Propagate the carry.
15:     for i = 0 to n − 1 do                   . The remainder is written in a, copy it to r.
16:         r[i] ← a[i]
        return

   We assume a primitive DIV _2 BY 1 that divides a 2-limb integer by a 1-limb integer and
returns the quotient. It has no WhyML code and we assume that the hardware provides
such a function. It is the only division primitive used by the functions in this section.
 val div_2by1 (l h d:limb) : limb
   requires { to_int h < to_int d }
   ensures { to_int result = div (to_int l + (max_uint64+1) * to_int h) (to_int d) }

   The algorithm consists in computing the limbs of the quotient one by one, starting with
the most significant. The numerator is overwritten at each step to contain the partial re-
mainder.
   At each iteration of the loop, we compute a quotient limb and subtract from the current
remainder the product of that quotient limb and the denominator, left-shifted appropriately
to cancel out the most significant limb of the current remainder.
   To compute a quotient limb, a candidate value is first guessed by dividing the two most
significant limbs from the current remainder by the most significant limb of the denomina-
tor.
   This candidate value is then adjusted to match the correct value of the quotient (lines
5-8 and 11-14). This process is called “adjustment step” throughout this section. The
algorithms that are actually implemented in GMP are variants of Algorithm D that try to
minimize the number of adjustments that occur.

                                                       Journal of Formalized Reasoning Vol.12, No.1, 2019
66       ·      Raphaël Rieu-Helft

  This is where the requirement that d[n − 1] ≥ β/2 comes into play. When this is the
case (we call such a denominator normalized) then the initial 2-by-1 division gives a good
approximation of the target quotient limb.
  D EFINITION 1. An integer p[0] . . . p[n − 1] is said to be normalized when p[n − 1] ≥
β/2.

  predicate normalized (x:ptr limb) (sz:int32) =
      valid x sz   ∧ x[x.offset + sz - 1] ≥ div radix 2

   More precisely, as shown by Knuth [19, p.257, Theorem B], the candidate quotient
is at most too large by 2, under the condition that the denominator is normalized. The
denominator being normalized is therefore a precondition of Knuth’s algorithm, and of the
other division algorithms in this paper for similar reasons.
   In the general case, we remark that an integer is normalized if and only if its most
significant bit is set to 1. The denominator can therefore be normalized by counting the
leading zeros in the denominator, shifting the numerator and denominator by that amount
(the denominator is normalized), calling a division procedure, and correcting the output by
shifting the remainder to the right by the same amount.
   This normalization is done by a wrapper around the main division primitive. This wrap-
per is the function that is exposed to the user. We verified a version of the wrapper that is
very simple and only performs this normalization, so we will not discuss it in this paper.
GMP’s version also implements an alternative algorithm that is not needed for correctness,
but that is faster than the default one when the divisor is close in length to the dividend.
We have not implemented and verified this alternate algorithm yet; the more general algo-
rithm is used in all cases. In the rest of this paper, we will continue to assume the divisor
normalized.
4.1     General case algorithm
GMP does not use Knuth’s algorithm, but a similar one that uses a 3-by-2 division to
compute each quotient limb (Alg. 7). Let us now discuss the differences between this
algorithm and Algorithm D.
   First, there is an extra local variable x. It is really a proxy for the most significant limb
in the current remainder, in the sense that whenever we would read from a[n+i], we take x
instead. Thus, instead of being stored in a[i] . . . a[n + i], the current remainder is stored
in a[i] . . . a[n + i − 1] with the most significant limb stored separately in x. This saves a
few memory accesses, and it can be done because the most significant limb of the current
remainder is no longer needed after the current loop iteration.
   We are now better equipped to express the main loop invariants.
   Let X be equal to the initial value of value(a, m). The following invariants are main-
tained in the main loop:
  X = value(d, n) · β i · value(q + i, m − n − i) + value(a, n + i − 1) + β n+i−1 × x
                                   n−1
  value(a + i, n − 1) + β                · x < value(d, n)
  d[n − 2] + βd[n − 1] ≥ a[n + i − 2] + βx                   (implied by the previous two).
  An important difference with Algorithm D is that instead of dividing the 2 most signifi-
cant limbs of the numerator by the most significant limb of the denominator, we divide the
3 most significant limbs of the former by the 2 most significant limbs of the latter.
Journal of Formalized Reasoning Vol.12, No.1, 2019
                                               A Why3 proof of GMP algorithms                  ·      67

Algorithm 7 General case long division
Require: m ≥ n ≥ 3, valid(a, m), valid(d, n), valid(q, m − n)
Require: d[n − 1] ≥ β/2
Require: value(a + m − n, n) < value(d, n)            . Otherwise an extra quotient limb is
    needed
Ensure: value(old a, m) = value(q, m − n) × value(d, n) + value(a, n)
Ensure: value(a, n) < value(d, n)
 1: function DIV _ SB _ QR (q, a, d, m, n)
 2:    v ← RECIPROCAL _ WORD _3 BY 2(d[n − 1], d[n − 2])
 3:    x ← a[m − 1]
 4:    i=m−n
 5:    while i > 0 do
 6:        i←i−1
 7:        if x = d[n − 1] and a[n + i − 1] = d[n − 2] then                      . Unlikely.
 8:            q̂ ← β − 1
 9:             SUBMUL _1(a + i, d, n, q̂)                . We know the result is d[n − 1].
10:            x ← a[n + i − 1]
11:        else
12:            (q̂, x, l) ← DIV 3 BY 2_ INV(x, a[n+i−1], a[n+i−2], d[n−1], d[n−2], v)
13:            b ← SUBMUL _1(a + i, d, n − 2, q̂)
14:            b1 ← (l < b)                  . Last two steps of the subtraction are inlined.
15:            a[i + n − 2] ← l − b mod β
16:            b2 ← (x < b1 )
17:            x ← x − b1 mod β                                        . Finish subtraction.
18:            if b2 6= 0 then                                      . Unlikely, and b2 = 1.
19:                 q̂ ← q̂ − 1                              . We only need to adjust by 1.
20:                 c ← ADD _ IN _ PLACE(a + i, d, n − 1) . Add only over n − 1 limbs.
21:                 x ← x + d[n − 1] + c mod β                  . The carry out is always 1.
22:        q[i] ← q̂
23:    a[n − 1] ← x

  This means that the adjustment step is much more efficient than that of Algorithm D.
Indeed, the candidate quotient obtained this way is very likely to be correct (Lemma 4).
  L EMMA 4. Let q̂ and r[0]r[1] the quotient and remainder of the division of
  a[n + i − 2]a[n + i − 1]x by d[n − 2]d[n − 1].
  If q̂ × d[0] . . . d[n − 1] > a[i] . . . a[n + i − 1]x, then r[1] = 0.
  P ROOF. We have a[n + i − 2]a[n + i − 1]x = d[n − 2]d[n − 1]q̂ + r[0]r[1].
  We also have a[i] . . . a[n + i − 1]x ≥ β n−2 a[n + i − 2]a[n + i − 1]x.
  If q̂ × d[0] . . . d[n − 1] > a[i] . . . a[n + i − 1]x, this implies:
               q̂ × d[0] . . . d[n − 1] > β n−2 (d[n − 2]d[n − 1]q̂ + r[0]r[1])
  However,
         q̂ × d[0] . . . d[n − 1] = q̂ × d[0] . . . d[n − 3] + β n−2 d[n − 2]d[n − 1]q̂
                                < β n−1 + β n−2 d[n − 2]d[n − 1]q̂.
                                                        Journal of Formalized Reasoning Vol.12, No.1, 2019
68       ·      Raphaël Rieu-Helft

   Therefore we must have β n−2 r[0]r[1] < β n−1 , which implies r[1] = 0.
   This means that the borrow at line 18 is only non-zero when the high limb of the re-
mainder returned by the 3-by-2 division at line 12 is zero, which is intuitively rare (if the
outcomes were evenly distributed, the probability would be 1/β). This lemma does not
have an equivalent in the Why3 proof, as it is not needed to prove functional correctness.
   Not only is the initial guess for the quotient very likely correct, but when it is not, it
is only too large by 1 and we can correct it with a single incrementation. Compare this
to Algorithm D, where two separate blocks were dedicated to adjusting the quotient, and
one of which could be executed twice. The following lemma justifies that the candidate
quotient is at worst too large by 1, which justifies that the adjustment step at line 18 is
needed at most once.
   L EMMA 5. Let q̂ and r[0]r[1] the quotient and remainder of the division of
   a[n + i − 2]a[n + i − 1]x by d[n − 2]d[n − 1].
   Then (q̂ − 1) × d[0] . . . d[n − 1] ≤ a[i] . . . a[n + i − 1]x.
     P ROOF. We have q̂ × d[n − 2]d[n − 1] ≤ a[n + i − 2]a[n + i − 1]x, so:
(q̂ − 1)d[0] . . . d[n − 1] ≤ (q̂ − 1)d[0] . . . d[n − 3] + β n−2 (q̂ − 1)d[n − 2]d[n − 1]
                                < β n−1 + β n−2 a[n + i − 2]a[n + i − 1]x − β n−2 d[n − 2]d[n − 1]
                                < β n−2 a[n + i − 2]a[n + i − 1]x (we have d[n − 1] ≥ β/2)
                                < a[i] . . . a[n + i − 1]x.
   The remainder of the 3-by-2 division is also used: instead of a simple long subtraction
over a length n like in Algorithm D, we perform a long subtraction over a length n − 2 only
and inline the last two steps. These last two steps consist simply in propagating the borrow
from the previous subtraction, as the result of the 3 most significant limbs of subtraction are
known to be lx0 in the absence of borrow (the postcondition of the division is exactly that
a[n + i − 2]a[n + i − 1]x = q̂ × d[n − 2]d[n − 1] + lx). We then propagate the borrow
on lx0. Hence, lines 13 to 17 are equivalent to computing the subtraction
                              a[i] . . . a[n + i − 1]x − q̂ × d[0] . . . d[n − 1]
returning b2 as borrow and writing the result in a[i] . . . a[n + i − 2]x (one limb fewer).
   If b2 = 0, the first invariant is maintained. Otherwise, there is an adjustment to make
(line 18): if the subtraction overflows, our candidate quotient q̂ was too large, we subtract
1 from it and add back value(d, n) to the remainder. This addition is done at lines 20-21.
The last limb is added separately because we save a few operations by ignoring the carry
out (we know it is equal to 1).
   It is efficient to do the subtraction first and potentially backtrack on it later because it
happens very rarely (Lemma 4), and evaluating the subtraction makes it easy to check if
the candidate quotient was correct with a simple integer comparison.
   The specification of the 3-by-2 division primitive ensures that lx < d[n − 2]d[n − 1],
hence the second and third invariants are maintained if there is no adjustment. If there is
an adjustment, the addition at line 20 may overflow, but the value of the top two limbs of
the remainder will still be lx + 1 ≤ d[n − 2]d[n − 1]. This shows that the third invariant
still holds in this case. For the second invariant, we use the fact that value(a + i, n − 1) +
β n−1 x − β n b2 is negative if b2 6= 0, so adding back value(d, n) maintains the second
invariant.
Journal of Formalized Reasoning Vol.12, No.1, 2019
                                                A Why3 proof of GMP algorithms                  ·      69

   There is an extra contingency in the main loop: when the two most significant limbs
of the denominator and the current numerator are identical, then we can skip the 3-by-2
division and adjustment step. Indeed, the candidate quotient output by the division would
necessarily be β, but the third invariant ensures that the long subtraction that would follow
would have a non-zero borrow, and that the adjustment step would knock the candidate
quotient down to β − 1. Therefore, we immediately write in β − 1 as quotient limb. With
this case out of the way, we may write the 3-by-2 division in such a way that its quotient
output is always a single limb, which would otherwise not be true in general.
4.2     3-by-2 division
Let us now take a closer look at the 3-by-2 division subroutine used by our division algo-
rithm. It was introduced by Möller and Granlund [22] and is used in GMP 6.1.2.
 predicate reciprocal_3by2 (v dh dl:limb) =
      v = div (radix*radix*radix -1) (dl + radix * dh) - radix


 let div3by2_inv (uh um ul dh dl v: limb) : (limb,limb,limb)
      requires { dh   ≥ div radix 2 }
      requires { reciprocal_3by2 v dh dl }
      requires { um + radix * uh < dl + radix * dh }
      returns { q, rl, rh → q * (dl + radix * dh) + rl + radix * rh
                             = ul + radix * (um + radix * uh) }
      returns { _q, rl, rh → 0 ≤ rl + radix * rh < dl + radix * dh }

  The algorithm takes a precomputed pseudo-inverse v of the denominator dl dh as an extra
parameter. More precisely,
                                     3        
                                       β −1
                               v=                − β.
                                      dl + βdh
The reason it is precomputed and passed as a parameter rather than computed on the fly is
that the caller function (Alg. 7) always uses the same denominator over a long division, so
it is much more efficient to compute the pseudo-inverse only once.
    The algorithm consists essentially in multiplying dl dh by its pseudo-inverse v and then
performing some simple adjustments. Remarkably, this means that no division primitive is
used. As division primitives tend to be much more expensive than additions or multiplica-
tions, this makes the algorithm a very efficient way to perform a 3-by-2 division.
    The “trick” is that the computation of the pseudo-inverse itself does use a division prim-
itive. In fact, computing the pseudo-inverse is about as complex as the 3-by-2 division
proper because of this division. However, over an m-by-n long division, m − n short divi-
sions are performed, all with the same precomputed pseudo-inverse, which amortizes that
cost.
    While the algorithm itself is short, its proof is non-trivial. The hardest part was directly
taken from Möller and Granlund’s on-paper proof [22, Theorem 3] and adapted for Why3.
The algorithm that computes the pseudo-inverse with a single 2-by-1 division primitive
was also taken from their paper.
4.3     Smaller cases: n = 1 and n = 2
The general case algorithm only handles the case where the denominator has length 3 or
more. Different algorithms are used for smaller denominators. When the divisor is exactly
one limb long, the schoolbook algorithm is used (Alg. 8).
                                                         Journal of Formalized Reasoning Vol.12, No.1, 2019
70       ·      Raphaël Rieu-Helft

Algorithm 8 Schoolbook division by a 1-limb number
Require: m ≥ 1, valid(q, m), valid(a, m), d ≥ β/2
Ensure: value(a, m) = result + d × value(q, m)
Ensure: result < d
 1: function DIVREM _1(q, a, m, d)
 2:    v ← INVERT _ LIMB(d)
 3:    r←0
 4:    i←m−1
 5:    while i ≥ 0 do
 6:        (q̂, r̂) ← DIV 2 BY 1_ INV(r, a[i], d, v)                      . Divide a[i]r by d
 7:        q[i] ← q
 8:        r ← r̂
 9:        i←i−1
10:    return r

   The variable r is the partial remainder, in the sense of the following loop invariants:
(1) value(a + i + 1, m − i − 1) = d × value(q + i + 1, m − i − 1) + r
(2) r < β.
Let us prove that the first invariant is maintained. Assume value(a + i + 1, m − i − 1) =
d × value(q + i + 1, m − i − 1) + r and r < β. Let q̂, r̂ such that q̂d + r̂ = a[i] + βr.
Then after q[i] ← q,
      value(a + (i − 1) + 1, m − (i − 1) − 1)
       = value(a + i, m − i)
       = a[i] + β value(a + i + 1, m − i − 1)                       [value_sub_head]
       = a[i] + βd × value(q + i + 1, m − i − 1) + βr
       = q̂d + r̂ + βd × value(q + i + 1, m − i − 1)
       = d × value(q + i, m − i) + r̂
       = d × value(q + (i − 1) + 1, m − (i − 1) − 1) + r̂.
   Finally, similarly to the general algorithm, we precompute a pseudo-inverse v of d for
the 2-by-1 division, and subsequently
                                    j 2useka 2-by-1 division algorithm that uses no division
primitive. This time, we have v = β d−1 − β.

  predicate reciprocal (v d:limb) =
     v = (div (radix*radix - 1) (d)) - radix


  let div2by1_inv (uh ul d v:limb) : (limb,limb)
     requires { d    ≥ div radix 2 }
     requires { uh    < d }
     requires { reciprocal v d }
     returns { q, r     → q * d + r = ul + radix * uh }
     returns { q, r     → 0 ≤ r < d }

  When n = 2, a very similar algorithm is used (Alg. 9). The only difference is that 3-
by-2 division is used at each loop iteration. Similarly to the n = 1 case, lh is the current

Journal of Formalized Reasoning Vol.12, No.1, 2019
                                              A Why3 proof of GMP algorithms                  ·      71

Algorithm 9 Schoolbook division by a 2-limb number
Require: m ≥ 2, valid(q, m − 2), valid(a, m), valid(d, 2), d[1] ≥ β/2
Require: value(a + m − 2, 2) < value(d, 2)              . Otherwise an extra quotient limb is
    needed
Ensure: value(old a, m) = value(d, 2) × value(q, m − 2) + value(a, 2)
Ensure: value(a, 2) < value(d, 2)
 1: function DIVREM _2(q, a, m, d)
 2:    v ← RECIPROCAL _ WORD _3 BY 2(d[1], d[0])
 3:    h ← a[m − 1]
 4:    l ← a[m − 2]
 5:    i←m−2
 6:    while i > 0 do
 7:        (q̂, l, h) ← DIV 3 BY 2_ INV(h, l, a[i − 1], d[1], d[0], v)
 8:        i←i−1
 9:        q[i] ← q̂
10:    a[1] ← h
11:    a[0] ← l

remainder, and the loop invariants are as follows:
     value(a, m)   = value(a, i) + β i (value(q + i, m − i − 2) × value(d, 2) + l + βh)
          l + βh < value(d, 2).
   Let us now discuss the differences between these two special cases on one hand, and the
general case algorithm in the other hand. The adjustment step that is found in Algorithm
D and other general long division algorithm is not needed in any of these two special
cases. Indeed, the source of error in Algorithm D (and in GMP’s algorithm) that makes the
adjustment step necessary is that the candidate quotient is computed using an incomplete
denominator. When the denominator has length 1 or 2, the division that occurs at each loop
iteration uses the full denominator. Even though the numerator still has arbitrary length,
this does not cause errors on the candidate quotient.
   Another difference is that no subtraction is needed to compute the partial remainders in
the two short cases. The corresponding operation in Alg. 8 is simply the r ← r̂ assignment
at line 8. Indeed, we have dq̂ + r̂ = a[i]r from the 2-by-1 division, and the long subtraction
that we would perform with Algorithm D would be exactly a[i]r−dq̂ = r̂, so decrementing
i and assigning r to r̂ does the same thing and saves a subtraction. Again, this is due to the
fact that we use the whole denominator at each loop iteration.

5.    TOOM-COOK MULTIPLICATION
The schoolbook multiplication algorithm from Section 3.3 has quadratic complexity and
is only optimal for numbers shorter than about 30 limbs, or 2000 bits. For larger numbers
(between 2000 and about 100,000 bits), GMP uses a family of recursive multiplication
algorithms initially introduced by Toom [27] and Cook [8]. These algorithms for inte-
ger and polynomial multiplication can be viewed as solving a multipoint evaluation and
polynomial interpolation problem.
   The general principle of Toom-Cook algorithms is to choose a base B, typically a power
of 264 , and to view the digits of the factors in base B as coefficients of polynomials a
                                                       Journal of Formalized Reasoning Vol.12, No.1, 2019
72       ·      Raphaël Rieu-Helft

and b. We then evaluate those polynomials at well-chosen points vi , compute the products
a(vi )b(vi ) by calling the algorithm recursively, and interpolate to obtain the coefficients of
the product polynomial c. The product is then obtained by evaluating c(B).
   We have verified two Toom-Cook algorithms: Toom-2 (Sec. 5.1), which is similar to
Karatsuba multiplication [18], and its unbalanced variant Toom-2.5 (Sec. 5.2), introduced
by Bodrato and Zanoni [5].
   Toom-2 (also called toom22 in the code) splits each of the operands into two parts
roughly equal in length, and Toom-2.5 (or toom32) splits the largest operand into three
parts and the smallest into two. Toom-2 is called on operands of roughly equal length
and Toom-2.5 is called when one of the operands is about 1.5 times as long as the other.
This way, after splitting, we are left with parts that have roughly equal length. A general
case algorithm, which we describe in Sec. 5.3, reduces all cases to applications of the two
former ones.
5.1     Toom-2

Algorithm 10 Toom-2 multiplication
Require: 2 ≤ n ≤ m < 30 × 2k
Require: valid(r, m + n), valid(a, m), valid(b, n), valid(s, 2(m + k))
Require: 4 · m < 5 · n
Ensure: value(r, m + n) = value(a, m) × value(b, n)
 1: function TOOM 22_ MUL (r, a, b, s, m, n, k)
 2:    µ←m1
 3:    λ←m−µ
 4:    ν ←n−λ
 5:    (a0 , a1 ) ← (a, a + λ)
 6:    (b0 , b1 ) ← (b, b + λ)
 7:    Compute |a(−1)| in r, |b(−1)| in r + λ, sign in  (see Alg. 12)
 8:    (c0 , c∞ ) ← (r, r + 2λ)
 9:    s0 ← s + 2λ
10:    TOOM 22_ MUL _ REC (s, r, r + λ, s0 , λ, λ, k − 1) . Compute |c(−1)| recursively.
11:    TOOM 22_ MUL _ REC (c∞ , a1 , b1 , s0 , µ, ν, k − 1) . Compute c(+∞) recursively.
12:    TOOM 22_ MUL _ REC (c0 , a0 , b0 , s0 , λ, λ, k − 1)    . Compute c(0) recursively.
13:    v ← ADD _ N(c∞ , c0 + λ, c∞ , λ)                                      . H0 + L∞
14:    v2 ← v + ADD _ N(c0 + λ, c∞ , c0 , λ)                            . L0 + H0 + L∞
15:    v ← v + ADD(c∞ , c∞ , λ, c∞ + λ, µ + ν − λ)                    . H0 + L∞ + H∞
16:    if  = 1 then
17:         v ← v + ADD(r + λ, r + λ, c−1 , 2λ)
18:    else v ← v − SUB(r + λ, r + λ, c−1 , 2λ) mod β              . v ∈ {0, 1, 2, β − 1}.
19:    INCR (c∞ , v2 )
20:    if v ≤ 2 then                                                     . Implies 0 ≤ v.
21:         INCR (r + 3λ, v)
22:    else DECR(r + 3λ, 1) . v = β − 1 instead of −1, due to integer representation.




Journal of Formalized Reasoning Vol.12, No.1, 2019
A Why3 proof of GMP algorithms                                                             73


Algorithm 11 Recursive call in Toom-2
Require: valid(r, m + n), valid(a, m), valid(b, n), valid(s, 2(m + k))
Require: 4m < 5n
Ensure: value(r, m + n) = value(a, m) × value(b, n)
  function TOOM 22_ MUL _ REC(r, a, b, s, m, n, k)
     if b < 30 then                    . Operands are small, use the schoolbook algorithm.
          MUL _ BASECASE(r, a, b, m, n)
     else
          if 4 · m < 5 · n then
               TOOM 22_ MUL (r, a, b, s, m, n, k)
          else                                   . Operands are unbalanced, use Toom-2.5.
               TOOM 32_ MUL (r, a, b, s, m, n, k)




Algorithm 12 Computation of |a(−1)| in r and |b(−1)| in r + λ
Ensure: value(r, λ) = |a(−1)|
Ensure: value(r + λ, λ) = |b(−1)|
Ensure:  · value(r, λ) · value(r + λ, λ) = a(−1)b(−1)
  ←1                                                 . Will hold the sign of a(−1)b(−1).
  if λ = µ then                                                           . Compute a(−1).
      if COMPARE(a0 , a1 , λ) < 0 then                                            . A1 > A0
           SUB (r, a1 , a0 , λ)
            ← −1
      else SUB(r, a0 , a1 , λ)
  else                                                                           .λ=µ+1
      if a0 [µ] = 0 ∧ COMPARE(a0 , a1 , µ) < 0 then                               . A1 > A0
           SUB (r, a1 , a0 , µ)
           r[µ] ← 0
            ← −1
      else
           t ← SUB(r, a0 , a1 , µ)
           r[µ] ← a0 [µ] − t                              . No carry, as we know a0 ≥ a1 .
  if λ = ν then                                                           . Compute b(−1).
      if COMPARE(b0 , b1 , λ) < 0 then                                            . B1 > B0
           SUB (r + λ, b1 , b0 , λ)
            ← −                                       . Change the sign of a(−1)b(−1).
      else SUB(r + λ, b0 , b1 , λ)
  else
      if IS _ ZERO(b0 + ν, λ − ν) ∧ COMPARE(b0 , b1 , ν) < 0 then                 . B1 > B0
           SUB (r + λ, b1 , b0 , ν)                          . b0 also has length at most ν.
           ZERO (r + λ + ν, λ − ν)       . We still have to initialize the rest of (r + λ, λ).
            ← −
      else SUB(r + λ, b0 , λ, b1 , ν)
74       ·      Raphaël Rieu-Helft

   The parameters of the toom22_mul function (Alg. 10) are as follows: r is the destination
buffer, a and b are the source operands, m and n are their lengths in limbs, and s is an extra
buffer to store temporary results.
   The amount of space needed for s is approximately 2(m + log2 (m)) limbs. Rather
than explicitly talking about logarithms in the specification, we use an extra parameter
k ≥ log2 (m). The variable k is ghost [12], which means that it is never used in the
computations, but only in the proof. The first precondition ensures that k is a suitable
bound. In practice, the caller of Toom-2 can give k = 64 (and allocate 2m + 128 limbs as
scratch space), as integer lengths are machine integers, so smaller than 264 .
   The last precondition makes sure that the operands are sufficiently close in size. It could
be a bit looser without breaking the algorithm, but the unbalanced version Toom-2.5 is
more efficient when the operand sizes are too unbalanced.
   The algorithm is organized in four steps. First we split the operands into two parts
of roughly equal length (Sec. 5.1.1). Then we evaluate the product polynomial c at three
points (Sec. 5.1.2). We then recompose the coefficients of c through interpolation (Sec. 5.1.3).
Finally, we propagate the remaining carries (Sec. 5.1.4).
  5.1.1 Splitting (Alg. 10, lines 2-6). We pose µ = b m
                                                      2 c, λ = m − µ, ν = n − λ. The
preconditions ensure the following:
                                              0<ν≤µ≤λ<m
                                              λ−1≤µ≤λ
                                                     λ<µ+ν
We can split a into two subwords a0 and a1 such that value(a, m) = value(a0 , λ) + β λ ×
value(a1 , µ). Similarly we have b0 and b1 such that value(b, n) = value(b0 , λ) + β λ ×
value(b1 , ν).
  We denote A0 := value(a0 , λ) and so on. Likewise, we define the polynomials a(X) :=
A0 + A1 X and b(X) := B0 + B1 X. The goal is to compute c(β λ ), where c(X) =
a(X)b(X) is a degree-2 polynomial. We pose c(X) = C0 + C1 X + C2 X 2 .

   a              A0                    A1             b     B0        B1
                   λ                     µ                   λ          ν


   5.1.2 Evaluation (Alg. 10, lines 7-12). The first step is to obtain three values of c(X)
for interpolation. GMP chooses to evaluate c at 0, −1 and +∞ (where c(+∞) is defined as
C2 ). We first evaluate a(−1) = A0 − A1 and b(−1) = B0 − B1 (Alg. 12). To avoid carry
propagations, we first check which of A0 and A1 is larger to compute |a(−1)| and store
its sign separately in a variable . If µ = λ − 1, we can optimise slightly by performing a
subtraction of length µ instead of λ.
   Similarly, we compute |b(−1)| and update  to contain the sign of a(−1)b(−1). We
store |a(−1)| in the first λ limbs of r (we denote this subarray r(0, λ)) and |b(−1)| in the
next λ limbs (which we denote r(λ, 2λ)). We then call Toom-2 recursively to compute
|a(−1)b(−1)| and store the result in s(0, 2λ) (Alg. 10, line 10). We use s(2λ, . . .) as
scratch space (there is enough space because k decreased by one).
   The constant 30 in the recursive call function TOOM 22_ MUL _ REC (Alg. 11) is the min-
imum integer length for which our library calls Toom-Cook multiplication algorithms. For

Journal of Formalized Reasoning Vol.12, No.1, 2019
                                              A Why3 proof of GMP algorithms                      ·      75

smaller numbers, the schoolbook multiplication is called instead. The value 30 was picked
experimentally, and the optimum can vary from machine to machine.
   After the recursive call, the memory layout is as follows. Throughout this section, the
memory layout will be illustrated with such diagrams (in addition to the formulas) to help
visualise the current state of the computation.

  r        |a(−1)|            |b(−1)|
              λ                  λ                     µ+ν


  s               |a(−1)b(−1)|
                       2λ                              remaining scratch space


   The next step is to compute c(+∞), that is, A1 B1 . This is done with a simple recursive
call to Toom-2 (Alg. 10, line 11), using s(2λ, . . .) as scratch space again. The result has
size µ + ν, which fits in r(2λ, 2λ + µ + ν).
   Finally, we compute c(0), that is, A0 B0 (line 12). We use the second half of s as scratch
space again, and store the result in r(0, 2λ), writing over |a(−1)| and |b(−1)| (but their
product still is in the first half of s).

  r                  A0 B 0                         A1 B 1
                       2λ                              µ+ν


  s               |a(−1)b(−1)|
                       2λ                              remaining scratch space


  We further decompose c(0), c(−1) and c(+∞) in halves of size λ or less:
                            A0 B0 = c(0) = L0 + β λ × H0
                            A1 B1 = c(+∞) = L∞ + β λ × H∞ .
At the end of the evaluation step, we have the following memory layout:

  r          L0                  H0              L∞                 H∞
              λ                  λ                 λ              µ+ν−λ


  s               |a(−1)b(−1)|
                       2λ                              remaining scratch space


   5.1.3 Recomposition (Alg. 10, lines 13-15). We first add H0 to L∞ in place, storing
the carry in a variable v (line 13). The outgoing up-right arrows in the diagrams that follow
represent the carry out.
                                                       v
                                                  →




             L0                                  L∞                 H∞
  r
              +                                   H0
              λ                  λ                 λ              µ+ν−λ



                                                           Journal of Formalized Reasoning Vol.12, No.1, 2019
76       ·      Raphaël Rieu-Helft

   We then add the result to L0 , writing over the original location of H0 (line 14). We store
the sum of v and the new carry in v2 . At that point, value(r+λ, λ)+β λ v2 = L0 +H0 +L∞ .




                                               2
                                               v




                                                                        v
                                          →




                                                                   →
                  L0                     H0                        L∞       H∞
   r               +                     L0                        H0
                   +                     L∞
                   λ                       λ                        λ       µ+ν−λ


   Finally, we add H∞ to r(2λ, 3λ) in place, incrementing v if needed (line 15).
   At that point, we have:
       value(r    + λ, 2λ) + β λ v2 + β 2λ v = (H0 + L0 + L∞ ) + β λ (H0 + L∞ + H∞ )
                                                          = A0 B0 + A1 B1 + H0 + β λ L∞ .   (1)
                                               2
                                               v




                                                                        v
                                          →




                                                                   →

                  L0                     H0                        L∞       H∞
   r               +                     L0                        H0
                   +                     L∞                        H∞
                   λ                       λ                        λ       µ+ν−λ


  Finally, we subtract c(−1) from r(λ, 3λ) by adding or subtracting |a(−1)||b(−1)|, de-
pending on the stored sign (lines 16-18).
  At that point, we have −1 ≤ v ≤ 3, and:
value(r      + λ, 2λ) + β λ v2 + β 2λ v = A0 B0 + A1 B1 + H0 + β λ L∞ − (A0 − A1 )(B0 − B1 )
                                                   = A0 B1 + A1 B0 + H0 + β λ L∞ .
                                               2
                                               v




                                                                        v
                                          →




                                                                   →




                  L0                     H0                        L∞       H∞
   r
                   +                           A0 B 1 + A1 B 0
                   λ                                      2λ                µ+ν−λ


   Therefore,
                                                     2λ
                       value(r, m + n) + β                v2 + β 3λ v
                       = L0 + β λ (value(r + λ, 2λ) + β λ v2 + β 2λ v) + β 3λ H∞
                       = A0 B0 + β λ (A0 B1 + A1 B0 ) + β 2λ A1 B1
                       = (A0 + β λ A1 )(B0 + β λ B1 ).
   The only thing left to do is propagating the carries.


Journal of Formalized Reasoning Vol.12, No.1, 2019
                                              A Why3 proof of GMP algorithms                  ·      77

   5.1.4 Carry propagation (Alg. 10, lines 19-22). We first propagate v2 at line 19, then
v (lines 20-22). There is an if statement because the case v2 = −1 (represented as the
unsigned integer β − 1) requires special treatment, as INCR(·, β − 1) is not the same thing
as DECR(·, 1). Another thing to look out for is that the functions INCR and DECR that
propagate the carries never check whether they reach the bounds of the array that holds the
number to be incremented. Rather, they perform an addition (or subtraction) and propagate
the carry until there is none. Their preconditions include the fact that this computation
should not overflow. Calling these functions incorrectly could result in buffer overflows.
This caused various bugs in past versions of GMP where this precondition was mistakenly
violated in rare cases. Notably, this is a situation where the memory safety of the program
(absence of buffer overflow) directly depends on its functional correctness (the number that
is being incremented fits between certain bounds).
   The most difficult part of the whole Toom-2 proof is ensuring that the propagation of the
carries v and v2 does not overflow. This is easy in the case where v ≥ 0, as the total product
ab (obtained after propagating both carries) is certain to fit in m + n and the intermediate
value obtained after propagating one of the carries is certain to be between 0 and ab.
   The only nontrivial case is v = −1, v2 6= 0. In this case, it is not obvious that the first
propagation, that of v2 , does not overflow out of r, as we have:
             value(r, m + n)   = (A0 + β λ A1 )(B0 + β λ B1 ) − β 2λ v2 + β 3λ ,
and we might have (A0 + β λ A1 )(B0 + β λ B1 ) + β 3λ ≥ β m+n .
  It is sufficient to show that H∞ < β µ+ν−λ − 1, that is, that the binary representation of
H∞ is not all ones. Indeed, the carry is absorbed somewhere in H∞ if it is the case.
  If v = −1, the first addition (of H0 and L∞ ) cannot overflow, and we must have v2 ≤ 1.
The only case to consider is therefore v = −1, v2 = 1.
  If v = −1, the substraction of c−1 from r + λ necessarily underflows, that is,
              A0 B0 + A1 B1 + H0 + β λ L∞ − β λ < (A0 − A1 )(B0 − B1 )
(Equation (1) with v2 = 1, v = 0).
   Noticing that 0 ≤ H0 and 0 ≤ A0 B1 + A1 B0 = A0 B0 + A1 B1 − (A0 − A1 )(B0 − B1 ),
we are left with β λ L∞ − β λ < 0, which implies L∞ = 0.
   Let us pose x, y the 2-adic valuations of A1 and B1 , and a0 , b0 odd integers such that
A1 = 2x a0 , B1 = 2y b0 . We have 2x+y a0 b0 = A1 B1 = β λ H∞ = 264λ H∞ . As a0 and
b0 are odd we must have x + y ≥ 64λ. If H∞ is even we are done, so we can assume
x + y = 64λ and a0 b0 = H∞ without loss of generality.
   Notice now that A1 < 264µ as it fits in a zone of length µ. We therefore have x < 64µ
and a0 ≤ 264µ−x − 1, similarly y < 64ν and b0 ≤ 264ν−y − 1.
   Therefore,
                      H∞ = a0 b0 ≤ (264µ−x − 1)(264ν−y − 1)
                                  = β µ+ν−λ − 264ν−y − 264µ−x + 1
                                  ≤ β µ+ν−λ − 3,
and we can conclude that the propagation of v2 does not overflow.
  We did not expect this fact to require such a complex proof, especially because there
were no comments in the GMP source code to indicate that non-trivial reasoning was


                                                       Journal of Formalized Reasoning Vol.12, No.1, 2019
78       ·      Raphaël Rieu-Helft

needed. We discussed this with the developers, and they ended up changing the algorithm
to make it more clearly correct at no performance cost.4
5.2     Toom-2.5
Toom-2.5, also known as toom32, has the same signature as Toom-2. The algorithm is sim-
ilar (Alg. 13), but splits the larger operand into three parts rather than two. It is optimally
used when the first operand is longer than the second by half.
   The algorithm is organized into phases that are similar to those of Toom-2. First we split
the larger operand into three parts and the smaller one into two parts (Sec. 5.2.1). Then we
evaluate c at four points (Sec. 5.2.2 and 5.2.4), and finally we recompose the coefficients
of c (Sec. 5.2.3 and 5.2.5).
   5.2.1 Splitting (Alg. 13, lines 2-9). We pose either λ = 1 + b m−1                 n−1
                                                                    3 c or λ = 1 + b 2 c,
whichever is longest. We also pose µ = m − 2λ and ν = n − λ. The preconditions ensure
that 0 < ν ≤ λ, 0 < µ ≤ λ, and λ ≤ µ + ν.
   Similarly to Toom-2, we split a into three subwords a0 , a1 and a2 of lengths λ, λ and µ
respectively. We also split b into two subwords b0 and b1 of lengths λ and ν. We denote
A0 := value(a0 , λ) and so on. We define the polynomials a(X) := A0 + A1 X + A2 X 2
and b(X) := B0 + B1 X. The goal is to compute c(β λ ), where c(X) = a(X)b(X) is a
degree-3 polynomial. We pose c(X) = C0 + C1 X + C2 X 2 + C3 X 3 .


                                               C 0 = A0 B 0
                                               C 1 = A1 B 0 + A0 B 1
                                               C 2 = A2 B 0 + A1 B 1
                                               C 3 = A2 B 1
   The lengths of Ai and Bi coefficients imply that C0 has length 2λ, C1 and C2 have
length 2λ + 1 and C3 has length µ + ν.

   a              A0                     A1                   A2       b     B0           B1
                   λ                       λ                  µ              λ             ν


   5.2.2 Evaluation in 1 and −1 (Alg. 13, lines 10-25, and full Alg. 14). As c(X) has de-
gree 3, this time we need to obtain four values for interpolation. GMP chooses to evaluate
c at 0, −1, 1 and +∞. We remark that a(−1) = A0 − A1 + A2 and a(1) = A0 + A1 + A2 .
Therefore, some evaluation steps can be saved by computing A0 + A2 only once and using
this result for both evaluations (Alg. 14). Just like in Toom-2, we actually compute |a(−1)|
and |b(−1)| and store the sign of the product separately in a variable . The carries of the
evaluation of a(1), b(1) and |a(−1)| are stored in a] , b] and a[ respectively. The evaluation
of |b(−1)| = |B0 − B1 | cannot overflow (and this is why we bother computing absolute
values and storing the sign separately).

4 https://gmplib.org/repo/gmp/rev/02a2ec6e1bce




Journal of Formalized Reasoning Vol.12, No.1, 2019
A Why3 proof of GMP algorithms                                                                79




Algorithm 13 Toom-2.5 multiplication
Require: 30 ≤ n ≤ m < 30 × 2k
Require: valid(r, m + n), valid(a, m), valid(b, n), valid(s, 2(m + k))
Require: n + 2 ≤ m ∧ m + 6 ≤ 3 · n
Ensure: value(r, m + n) = value(a, m) × value(b, n)
 1: function TOOM 32_ MUL (r, a, b, s, m, n, k)
 2:    if 2m ≥ 3n then
 3:        λ ← 1 + b m−1
                      3 c
 4:    else
 5:        λ ← 1 + b n−1
                      2 c
 6:       µ ← m − 2λ
 7:       ν ←n−λ
 8:       (a0 , a1 , a2 ) ← (a, a + λ, a + 2λ)
 9:       (b0 , b1 ) ← (b, b + λ)
            Compute a(1), b(1), |a(−1)|, |b(−1)| in r, sign in  (see Alg. 14).
            Carries of a(1) and b(1) are in a] and b] , carry of |a(−1)| in a[ .
10:       s0 ← s + 2λ + 1
11:       TOOM 22_ MUL _ REC (s, r, r + λ, s0 , λ, λ, k − 1)        . Compute c(1) recursively.
12:       if a] = 1 then                                          . Propagate the carry for c(1).
13:            v ← b] + ADD _ N(s + λ, s + λ, r + λ, λ)
14:       else
15:            if a] = 2 then
16:                 v ← 2b] + ADDMUL _1(s + λ, r + λ, λ, 2)
17:            else                                                                     . a] = 0
18:                 v←0
19:       if b] then
20:            v ← v + ADD _ N(s + λ, s + λ, r, λ)
21:       s[2λ] ← v
22:       TOOM 22_ MUL _ REC (r, r + 2λ, r + 3λ, s0 , λ, λ, k − 1)          . Compute |c(−1)|
      recursively.
23:       if a[ then                                          . Propagate the carry for |c(−1)|.
24:            a[ ← ADD _ N(r + λ, r + λ, r + 3λ, λ)
25:       r[2λ] ← a[
26:       if  = −1 then
27:            SUB _ N (s, s, r, 2λ + 1)
28:       else
29:            ADD _ N (s, s, r, 2λ + 1)
30:       RSHIFT(s, s, 2λ + 1, 1)                                 . s ← c(1)+c(−1)
                                                                              2      = C0 + C2
31:       v ← ADD _ N(r + 2λ, s, s + λ, λ)                               . Add L and M in D1 .
32:       INCR (s + λ, λ + 1, v + s[2λ])                             . Propagate v and h to D2 .
80       ·      Raphaël Rieu-Helft

33:       if  = −1 then
34:            v ← ADD _ N(s, s, r, λ)
35:            v 0 ← r[2λ] + ADD _ NC(r + 2λ, r + 2λ, r + λ, λ, v)
36:            INCR (s + λ, λ + 1, v 0 )
37:       else
38:            v ← SUB _ N(s, s, r, λ)
39:            v 0 ← r[2λ] + SUB _ NC(r + 2λ, r + 2λ, r + λ, λ, v)
40:            DECR (s + λ, λ + 1, v 0 )
41:       TOOM 22_ MUL _ REC (r, a0 , b0 , s0 , λ, λ, k − 1)     . Compute c(0) recursively.
42:       if µ > ν then                                       . Compute c(+∞) recursively.
43:            MUL (r + 3λ, a2 , µ, b1 , ν)
44:       else
45:            MUL (r + 3λ), b1 , ν, a2 , µ)
46:       v ← SUB _ N(r + λ, r + λ, r + 3λ, λ)
47:       v 0 ← s[2λ] + v
48:       v ← SUB _ NC(r + 2λ, r + 2λ, r, λ, v)
49:       v 0 ← v 0 − SUB _ NC(r + 3λ, s + λ, r + λ, λ, v)
50:       v 0 ← v 0 + ADD(r + λ, r + λ, 3λ, s, λ)
51:       if µ + ν > λ then                                                 . Propagate v 0 .
52:            v 0 ← v 0 − SUB(r + 2λ, r + 2λ, 2λ, r + 4λ, µ + ν − λ)
53:            if v 0 < 0 then
54:                  DECR (r + 4λ, µ + ν − λ, −v 0 )
55:            else
56:                  INCR (r + 4λ, µ + ν − λ, v 0 )
                       ]




                                                                 [
                                               ]
                       a




                                                                 a
                                               b
                  →




                                          →




                                                            →




   r            a(1)                     b(1)            |a(−1)|            |b(−1)|
                   λ                       λ                 λ                   λ            µ+ν−λ


  After this step, we have a(1) = value(r, λ) + β λ a] and so on. A simple case analysis
shows 0 ≤ a] ≤ 2, 0 ≤ b] ≤ 1, and 0 ≤ a[ ≤ 1. The next thing to do is to compute c(1).
             c(1) = a(1)b(1) = (value(r, λ) + β λ a] ) · (value(r + λ, λ) + β λ b] )
  The recursive call at line 11 of Alg. 13 computes value(r, λ) · value(r + λ, λ) in s. The
if statement at line 12 adds a] value(r + λ, λ), shifted by λ as this term must be multiplied
by β λ . It also accumulates the carry and a] b] into the variable v. After this, we have
                                 2λ
         value(s, 2λ) + β             v = (value(r, λ) + β λ a] ) · value(r + λ, λ) + β 2λ a] b]
The only missing term is β λ b] · value(r, λ), which is handled in the if statement at line
16 of Alg. 13. After this, we have value(s, 2λ) + β 2λ v = c(1), and we can set s[2λ] to v.

   s                          c(1)
                             2λ + 1                                        remaining scratch space


  The product |c(−1)| = |a(−1)b(−1)| is computed similarly at lines 22-29. As there is
no carry corresponding to |b(−1)|, the procedure is a bit simpler. The term β λ a[ · value(r+
Journal of Formalized Reasoning Vol.12, No.1, 2019
                                               A Why3 proof of GMP algorithms                 ·      81

Algorithm 14 Toom-2.5: evaluation in 1 and -1
Ensure: value(r, λ) + β λ a] = a(1)
Ensure: value(r + λ, λ) + β λ b] = b(1)
Ensure: value(r + 2λ, λ) + β λ a[ = |a(−1)|
Ensure: value(r + 3λ, λ) = |b(−1)|
Ensure:  · (value(r + 2λ, λ) + β λ a[ ) · value(r + 3λ, λ) = a(−1)b(−1)
  a0 ← ADD(r, a0 , λ, a2 , µ)
  if a0 = 0 ∧ CMP(r, a1 , λ) < 0 then                                        . A0 + A2 < A1
       SUB _ N (r + 2λ, a1 , r, λ)
      a[ ← 0
       ← −1
  else                                                                       . A 1 ≤ A0 + A2
      v ← SUB _ N(r + 2λ, r, a1 , λ)
      a[ ← a0 − v
      ←1
  a] ← a0 + ADD _ N(r, r, a1 , λ)                                   . Finish computing a(1).
  if λ = ν then
      b] ← ADD _ N(r + λ, b0 , b1 , λ)
      if CMP(b0 , b1 , λ) < 0 then                                                  . B0 < B1
           SUB _ N (r + 3λ, b1 , b0 , λ)
            ← −
      else
           SUB _ N (r + 3λ, b0 , b1 , λ)
  else
      b] ← ADD(r + λ, b0 , b1 , λ, ν)
      if IS _ ZERO(b0 + ν, λ − ν) ∧ CMP(b0 , b1 , ν) < 0 then                       . B0 < B1
           SUB _ N (r + 3λ, b1 , b0 , ν)                      . b0 also has length at most ν.
           ZERO (r + 3λ + ν, λ − ν)      . We still have to initialize the rest of (r + 3λ, λ).
            ← −
      else
           SUB (r + 3λ, b0 , b1 , λ, ν)


3λ, λ) is added in the if statement at line 23. The product is stored in r(0, 2λ + 1). This
overwrites a(1), b(1) and part of |a(−1)|, but these intermediate results are no longer
needed.

  r                   |c(−1)|
                       2λ + 1                               λ+µ+ν−1


  s                    c(1)
                       2λ + 1                                      remaining scratch space


   5.2.3 Recomposition (Alg. 13, lines 26-40). We use the intermediate variable d :=
C1 + C3 + β λ · (C0 + C2 ). The sizes of the Ci imply that d has length 3λ + 1, so we pose
D0 , D1 , D2 such that d = D0 + β λ D1 + β 2λ D2 with D0 and D1 of length λ and D2 of
length λ + 1.
                                                       Journal of Formalized Reasoning Vol.12, No.1, 2019
82       ·       Raphaël Rieu-Helft

   Note that C1 + C3 = C0 + C2 − (C0 − C1 + C2 − C3 ) = C0 + C2 − c(−1), so:
                                  d = C0 + C2 − c(−1) + β λ · (C0 + C2 ).
   The computation of d goes as follows.

                                                     C0 + C2
          +                        C0 + C2
   d
             -                         c(−1)
          =            D0                      D1                 D2
                         λ                       λ                λ+1


   We compute D0 at s, D1 at r + 2λ and D2 at s + λ. The first step is to write C0 + C2 in
s. Noticing that c(1) + c(−1) = 2 · (C0 + C2 ), this is easy to do with one long addition or
subtraction (depending on ) of r and s (lines 26-29), and then a logical shift on the result
to do the division by 2 (line 30).

   r                         |c(−1)|
                              2λ + 1                                   λ+µ+ν−1


   s                         C0 + C2
                              2λ + 1                                      remaining scratch space


  Let us note h the highest limb of C0 + C2 and split the 2λ remaining limbs into two
subwords L and M , such that C0 + C2 = L + β λ M + β 2λ h.
  We add L and M together and write the result in r + 2λ, λ, the location of D1 (line 31
Alg. 13). This overwrites the most significant limb of |c(−1)|, but it is still stored in a[ .
The carry v is propagated and added to D2 , or s + λ, which already contains the higher
half of C0 + C2 (line 32).
                                  a[




                                                              v
                              →




                                                          →




   r                      |c(−1)|                       L+M
                              2λ                          λ                       µ+ν


   s              L               M + βλh + v + h
                   λ                       λ+1                            remaining scratch space


   The only thing left to do in the computation of d is to subtract c(−1). The if statement
at line 33 does so (by adding or subtracting |c(−1)| depending on ). As D0 and D1 are not
stored contiguously, we have to subtract both halves of c(−1) separately. The ADD _ NC and
SUB _ NC are variants of addition and subtraction that take a carry input to add or subtract
to the result, they allow to propagate the carry from the subtraction of the lower half of
c(−1) to the upper half. The carry of the upper half subtraction and a[ (the high limb of
|c(−1)|) are propagated to D2 .



Journal of Formalized Reasoning Vol.12, No.1, 2019
                                               A Why3 proof of GMP algorithms                   ·      83

  At that point, we can check the following:
                             λ
             value(s, λ) + β value(r    + 2λ, λ) + β 2λ value(s + λ, λ + 1)
             = L + β λ (L + M − β λ v) + β 2λ (M + β λ h + v + h) − c(−1)
             = L + β λ M + β 2λ h + β λ (L + M + β 2λ h) − c(−1)
             = C0 + C2 + β λ (C0 + C2 ) − c(−1)
             = d.


  r                                                 D1
                        2λ                          λ                    µ+ν


  s          D0                  D2
              λ                  λ+1                                 remaining scratch space


   5.2.4 Evaluation in 0 and +∞ (Alg. 13, lines 41-45). We recursively compute C0 in
(r, 2λ) and C3 in (r +3λ, µ+ν). The upper half of s is still available and is used as scratch
space. As the operands can be very unbalanced in the case of C3 , we have used the generic
multiplication (see Section 5.3) instead of calling Toom-2 or Toom-2.5 directly.

  r                     C0                          D1                         C3
                        2λ                          λ                    µ+ν


  s          D0                  D2
              λ                  λ+1                                 remaining scratch space


  5.2.5 Recomposition (Alg. 13, lines 46-56). Let us note L0 and H0 the two halves of
C0 , and split C3 = c(∞) into L∞ and H∞ of length λ and µ + ν − λ respectively.
                                    C0 = L0 + β λ H0
                                    C3 = L∞ + β λ H∞
  The product c(β λ ) = a(β λ )b(β λ ) can be expressed only in terms of C0 , C3 and d:
          C0 + β λ d + β 3λ C3 − C0 β 2λ − β λ C3
           = C0 + β λ (C1 + C3 ) + β 2λ (C0 + C2 ) + β 3λ C3 − C0 β 2λ − β λ C3
           = C0 + β λ C1 + β 2λ C2 + β 3λ C3
           = c(β λ ).
  This implies the following decomposition for c(β λ ):
c(β λ ) = C0 + β λ d + β 3λ C3 − C0 β 2λ − β λ C3
      = L0 + β λ (D0 + (H0 − L∞ )) + β 2λ (D1 − L0 − H∞ ) + β 3λ (D2 − (H0 − L∞ )) + β 4λ H∞ .
   The first step is to subtract L∞ from H0 at r + λ (lines 46-47). The borrow is stored in
v. The variable v 0 contains the sum of v and the high limb of D2 .

                                                         Journal of Formalized Reasoning Vol.12, No.1, 2019
84       ·      Raphaël Rieu-Helft




                                               v
                                          →
   r              L0                 H0 − L ∞              D1              L∞              H∞
                   λ                       λ                λ                λ            µ+ν−λ


   s             D0                        D2
                   λ                      λ+1                          remaining scratch space


  The next step is to subtract L0 from r + 2λ (line 48). The borrom v from the previous
computation at r + λ is also propagated into r + 2λ. The final carry is stored in v.




                                                                v
                                                           →
   r              L0                 H0 − L∞         →   D1 − L0           L∞              H∞
                   λ                       λ                λ                λ            µ+ν−λ


   s             D0                        D2
                   λ                      λ+1                          remaining scratch space


  We then subtract H0 − L∞ from D2 at r + 3λ, propagating the previous carry (line 49).
The carry out is accumulated in v 0 .




                                                                                  v  0
                                                                                 +
   r              L0                 H0 − L∞         →   D1 − L0            →
                                                                    →D2 − H0 + L∞ H∞
                   λ                       λ                λ                λ            µ+ν−λ


   s             D0                        D2
                   λ                      λ+1                          remaining scratch space


  We add D0 at r + λ, and propagate the carry all the way to r + 4λ (as the parameter
3λ is passed in the ADD call at line 50 of Alg. 13). Similarly, we then subtract H∞ from
r + 2λ and propagate the carry all the way to r + 4λ.
                                                                                 0
                                                                                 v
                                                                            →




   r              L0            D0 + H0 − L∞→D1 − L0 − H∞→D2 − H0 + L∞ H∞
                   λ                       λ                λ                λ            µ+ν−λ


   s             D0                        D2
                   λ                      λ+1                          remaining scratch space


   We then only have to propagate v 0 to r + 4λ to finish the computation. The proof that
this propagation does not overflow is much more straightforward than for Toom-2. Indeed,
this time there is only one carry to propagate, so we can rely on the fact that we know for
sure that the final result fits in space 3λ + µ + ν.
5.3     General case
When the two operands have similar sizes, we can compute their product using either
Toom-2 or Toom-2.5. When the operand sizes are very unbalanced, none of these algo-

Journal of Formalized Reasoning Vol.12, No.1, 2019
                                              A Why3 proof of GMP algorithms                 ·      85

rithms can be used directly. In this case, we perform a block product by decomposing
the larger operand in blocks of sizes 3/2 times that of the smaller operand, and calling
Toom-2.5 repeatedly (Alg. 15).

Algorithm 15 General case multiplication
Require: valid(r, m + n), valid(a, m), valid(b, n)
Require: 0 < n ≤ m
Ensure: value(r, m + n) = value(a, m) · value(b, n)
 1: function MUL (r, a, b, m, n)
 2:    if y < 30 then                         . Small operands, use schoolbook algorithm.
 3:         MUL _ BASECASE(r, a, b, m, n)
 4:    else
 5:        k ← 64
 6:        s ← ALLOC(5n + 128)                   . Allocate sufficiently large scratch space.
 7:        if 2m ≥ 5n then                      . Unbalanced operands, use block product.
 8:            m0 ← 3n/2                                                        . Block size.
 9:            w ← ALLOC(4n)               . Allocate workspace to store Toom-2.5 result.
10:            u←m                                     . Remaining section of a to multiply.
11:            TOOM 32_ MUL (r, a, b, s, m0 , n, k)
12:            u ← u − m0
13:            a0 ← a + m0
14:            r 0 ← r + m0
15:            while u ≥ 2n do
16:                TOOM 32_ MUL (w, a0 , b, s, m0 , n, k)
17:                v ← ADD _ N _ IN _ PLACE(r0 , w, n)              . Add result to subtotal.
18:                COPY (r + n, w + n, m0 )                         . Continue the addition.
19:                INCR (r + n, v)                                    . Propagate the carry.
20:                u ← u − m0
21:                a0 ← a0 + m0
22:                r 0 ← r 0 + m0
23:            if n ≤ u then                                       . Multiply the last block.
24:                if 4m < 5n then
25:                    TOOM 22_ MUL (w, a0 , b, s, u, n, k)
26:                else
27:                    TOOM 32_ MUL (w, a0 , b, s, u, n, k)
28:            else                                            . Operand sizes are reversed
29:                MUL (w, b, a0 , n, u)
30:            v ← ADD _ N _ IN _ PLACE(r0 , w, n)
31:            COPY (r + n, w + n, u)
32:            INCR (r + n, v)
33:        else
34:            if 4m < 5n then
35:                TOOM 22_ MUL (r, a, b, s, m, n, k)
36:            else
37:                TOOM 32_ MUL (r, a, b, s, m, n, k)



                                                      Journal of Formalized Reasoning Vol.12, No.1, 2019
86       ·      Raphaël Rieu-Helft

   This function is meant to be exposed to the user, with previous multiplication routines
left internal. It calls other multiplication algorithms that depend on the operand sizes. If
one operand is very small, the schoolbook multiplication is the best one (line 3). Otherwise,
if the operands are of sufficiently similar sizes, then we simply call Toom-2 or Toom-2.5
(line 34 of Alg. 15). Finally, if the operands are large and very unbalanced, we need to
perform a block product.
   The main loop invariants are as follow:
                       value(r, m + n − u)               = value(a, m − u) · value(b, n)
                                                     0
                                                 a =a+m−u
                                                 r0 = r + m − u.


   At the beginning of a loop iteration, a0 points to the first limb of a that has not been
multiplied yet, and r0 points to the zone in r where the next subresult should be added.
Note that the first n limbs after r0 already contain a part of the subtotal.
   After multiplying (a0 , m0 ) by (b, n), we need to add the product (stored in w) to r0 .
This is done in lines 17-19. Instead of performing one addition of length m0 + n, we take
advantage of the fact that only the first n limbs of r0 are occupied. We perform an addition
of length n on that zone and simply copy over the rest of w to r0 + n, and then propagate
the carry of the addition.
   What happens after the loop is very similar to one last iteration of the loop. The only
change is that the last block of a that is left to multiply has length u, not exactly 3n/2.
Depending on the ratio between n and u, either Toom-2 or Toom-2.5 is called. If u < n,
the recursive call at line 29 will necessarily jump to line 34 and call one of these two
functions, but with the order of the operands reversed to fulfill the precondition that the
larger operand is first.
   This algorithm is not exactly identical to GMP’s. GMP implements a third Toom-Cook
algorithm, toom42, which splits its larger operand into four parts and the smaller into two.
The general case uses this algorithm and splits its large operand in blocks of length 2n.
However, as toom42 is never called recursively, the asymptotic complexity is similar (the
exponent is the same, and the constant is about 7% worse). We did not verify toom42, as it
did not seem to present many new challenges compared to the other Toom-Cook algorithms
and the performance difference is not egregious.

6.    SQUARE ROOT
The GMP square root algorithm consists in four functions: two base case algorithms, for
operands of length 1 and 2 (sqrtrem1 and sqrtrem2 respectively), a general case algorithm
(dc_sqrtrem) that assumes its operand is normalized, and a wrapper function sqrtrem that
handles normalization. It computes the square root s and the remainder r of its operand a
such that a = s2 + r and 0 ≤ r ≤ 2s, or equivalently, s2 ≤ a < (s + 1)2 .
   The notion of normalization used here is slightly different from the division. The square
root algorithms that expect a normalized operand expect its high limb to be greater than or
equal to β/4, rather than β/2. In other terms, one of its two most significant bits must be
set. Note that the square root of such a number is normalized in the sense of the division,
which the general case square root algorithm makes use of.

Journal of Formalized Reasoning Vol.12, No.1, 2019
                                                   A Why3 proof of GMP algorithms                 ·      87

 6.1      Square root, n = 1
 The sqrtrem1 function computes the square root of a normalized one-limb machine inte-
 ger. It is used as the base case of GMP’s divide-and-conquer square root algorithm. We
 formally verified its correctness in previous work [21] using Why3 and the Gappa tool [9].
 Its implementation relies on C-specific features, such as the way numbers are represented
 in memory and the semantics of logical shifts for signed and unsigned types. Thus, a
 pseudocode would not make much sense and would be less readable than the actual C
 code. For the sake of completeness, we include a slightly edited version of the GMP C
 implementation in Alg. 16 and outline the algorithm briefly.

 Algorithm 16 Square root of a limb.
 1   #define MAGIC 0x10000000000
 2   /* 0xffe7debbfc < MAGIC < 0x232b1850f410 */
 3
 4   static mp_limb_t mpn_sqrtrem1 (mp_ptr rp, mp_limb_t a0) {
 5    mp_limb_t a1, x0, x1, x2, c, t, t1, t2, s;
 6       unsigned abits = a0 >> (64 - 1 - 8);
 7       x0 = 0x100 | invsqrttab[abits - 0x80];
 8       // x0 is the first approximation of 1/sqrt(a0)
 9       a1 = a0 >> (64 - 1 - 32);
10       t1 = (mp_limb_signed_t) (0x2000000000000 - 0x30000 - a1 * x0 * x0) >> 16;
11       x1 = (x0 << 16) + ((mp_limb_signed_t) (x0 * t1) >> (16+2));
12       // x1 is the second approximation of 1/sqrt(a0)
13       t2 = x1 * (a0 >> (32-8));
14       t = t2 >> 25;
15       t = ((mp_limb_signed_t) ((a0 << 14) - t * t - MAGIC) >> (32-8));
16       x2 = t2 + ((mp_limb_signed_t) (x1 * t) >> 15);
17       c = x2 >> 32;
18       // c is a full limb approximation of sqrt(a0)
19       s = c * c;
20       if (s + 2*c <= a0 - 1) {
21          s += 2*c + 1;
22            c++;
23        }
24       *rp = a0 - s;
25       return c;
26   }




    The program is best understood as a fixed-point arithmetic algorithm that computes the
 square root of a real number a ∈ [0.25, 1] using Newton’s method. Assuming we want to
 compute the square root of the integer a0 ∈ [262 , 264 − 1], we define a = 2−64 a0 . We start
 by fetching an 8-bit approximation x0 of a−1/2 from a precomputed table (line 7). Let us
 pose ε0 such that x0 = a−1/2 (1 + ε0 ). We have |ε0 | ≤ 2−8 .
    The algorithm consists in performing two steps of Newton’s method to reach a 32-bit
 approximation of a−1/2 . Equivalently, we are looking for roots of the function f (x) =
 x−2 −a. Given xi = a−1/2 (1+εi ), we compute xi+1 = xi −f (xi )/f 0 (xi ) = xi (3−ax2i ).
 The relative error εi decreases quadratically. More precisely, we find |εi+1 | ≈ 23 · |ε2i |.
 Intuitively, we need about two steps to reach a precision of 32 bits. The reason why we


                                                           Journal of Formalized Reasoning Vol.12, No.1, 2019
88       ·      Raphaël Rieu-Helft

are computing a−1/2 rather than a1/2 is that the Newton iteration for a1/2 would require a
division at each step, so it would be much more costly than the one for a−1/2 .
   The algorithm computes x1 between lines 9 and 11, and then x2 between lines 13 and
17. As part of the computation of x2 , there is an extra multiplication by a in order to
compute the square root of a rather than its inverse. We also multiply x2 by 232 to get back
                                                                               √
to computing the square root of a0 rather than a. At line 17, we obtain c = a0 (1 + ε2 ).
                                       −32
   At this point, we have ε2 ∈ [−2 , 0]. In other terms, either c is the square root of
a0 , or it is too small by one. The fact that |ε2 | is smaller than 2−32 is far from obvious.
Indeed, the constant 3/2 in the growth of εi is problematic. Furthermore, we are computing
using fixed-precision arithmetic rather than real numbers, so we introduce rounding errors.
However, the initial approximation x0 is actually a bit tighter than 8 bits (|ε0 | is smaller
than about 2−8.5 for all values of a0 ), which offsets the precision loss from the constant
3/2. Furthermore, the computations of x1 and x2 involve magic constants that would
not be present in the theoretical steps of the Newton iteration (0x30000 at line 10 and
MAGIC at line 15). These extra constants offset the rounding errors and ensure that the final
approximation is always by default.
   The final step of the algorithm (lines 19-25) performs a final adjustment step (in case we
             √                   √
computed a0 − 1 rather than a0 ) and computes the remainder.
6.2     Square root, n = 2
The sqrtrem2 function (Alg. 17) computes the square root of a, which must be two limbs
long. It takes two destinations operands, s to store the square root of a and r to store the
remainder. It returns the high limb c of the remainder, such that value(a, 2) = s[0]2 +
βc + r[0]. The parameters r and s are pointers that only need space for 1 limb. We assume
a to be normalized such that its high limb is greater than or equal to β/4. The algorithm
computes an initial estimate of the square root by calling sqrtrem1 on the high limb of a,
and then adjusts toward the correct square root.
  6.2.1 Initial estimate. After the call to sqrtrem1 at line 2, we have s20 + r[0] = a1
and a = 264 a1 + a0 . The value 232 s0 is the initial estimate for the square root of a. We
compute an initial remainder r0 such that 233 r0 + (a[0] mod 233 ) = βr[0] + a[0].
  6.2.2 Adjustment. We compute the quotient q of r0 by s0 , such that at line 6, s0 q is
somewhat close to r0 and q < 232 . We pose s1 = 232 s0 + q as a candidate square root.
Indeed, 233 qs0 ≈ βr[0] + a[0], so at line 8: s21 = βs20 + 233 qs0 + q 2 ≈ βa[1] + a[0] + q 2 .
More precisely,
                    s21 = βa[1] + a[0] + q 2 − (233 (r0 − qs0 ) + a[0] mod 233 ).
   We pose u = r0 − qs0 , and r1 , c such that r1 + βc = 233 (r0 − qs0 ) + a[0] mod 233 ,
r1 ∈ [0, β − 1]. Therefore, at line 11 we have s21 + r1 + βc = a + q 2 . At lines 12-13,
we subtract q 2 from r1 + βc in a way that avoids arithmetic underflows. At that point,
s21 + r1 + βc = a.
  6.2.3 No underapproximation of the square root. Let  √ us show that s1 is not an under-
approximation of the square root, or in other terms, b ac ≤ s1 . This is equivalent to
a < (s1 + 1)2 , or r1 + βc ≤ 2s1 . We have r1 + βc = 233 u + al − q 2 .
  If q0 < 232 , then q = q0 and u = r0 mod s0 < s0 , and al < 233 , so r1 + βc ≤ 233 s0 ≤
2s1 .

Journal of Formalized Reasoning Vol.12, No.1, 2019
                                             A Why3 proof of GMP algorithms                  ·      89

  The only remaining case is q0 = 232 , q = 232 −1. In this case, we notice that r[0] ≤ 2s0
(postcondition of sqrtrem1). This implies r0 ≤ 232 s0 + ah , so u = r0 − (232 − 1)s0 ≤
s0 + ah . We also have β = (q + 1)2 , so:
                          233 u + al − q 2 ≤ 233 s0 + a[0] − q 2
                                          ≤ 233 s0 + β − 1 − q 2
                                          = 233 s0 + 2q
                                          ≤ 2s1 .
                                                                        √
   6.2.4 Avoiding overapproximation. Let us now show that s1 < b ac + 1, or equiva-
lently, s21 ≤ a (no overapproximation of the square root). If we had s21 > a, it would imply
r1 + βc < 0, or equivalently, c ≤ −1. The conditional at line 14 takes care of this case by
adding 2s1 − 1 to r1 + βc and removing 1 from s1 , leaving the sum s21 + r1 + βc unchanged
and avoiding arithmetic overflows.

Algorithm 17 Square root of a 2-limb number
Require: valid(a, 2), valid(s, 1), valid(r, 1)
Require: a[1] ≥ β/4
Ensure: value(a, 2) = s[0]2 + β × result + r[0]
Ensure: 0 ≤ result ≤ 1
Ensure: r[0] + β × result ≤ 2s[0]
 1: function SQRTREM 2(s, r, a)
 2:    s0 ← SQRTREM 1(r, a[1])
 3:    ah ← a[0]  33
 4:    r0 ← (r[0]  31) + ah
 5:    q0 ← br0 /s0 c                                                             . q0 ≤ 232 .
                                                                                32
 6:    q ← q0 − (q0  32)                                        . If q0 = 2 , reduce it by 1.
 7:    u ← r0 − qs0
 8:    s1 ← (s0  32) + q
 9:    c ← u  31
10:    al ← a[0] mod 233
11:    r1 ← (u  33) + al
12:    c ← c − (r1 < q 2 )                                                                 . −1 ≤ c
13:    r1 ← r1 − q 2 mod β
14:    if c < 0 then                                          . Square root too large, adjust.
15:         r1 ← r1 + s1 mod β
16:         c ← c + (r1 < s1 )                                             . Carry propagation.
17:         s1 ← s1 − 1
18:         r1 ← r1 + s1 mod β
19:         c ← c + (r1 < s1 )                                             . Carry propagation.
20:    r[0] ← r1
21:    s[0] ← s1
22:    return c




                                                      Journal of Formalized Reasoning Vol.12, No.1, 2019
90       ·      Raphaël Rieu-Helft

6.3     Square root, general case
Our Why3 proof of the general case divide-and-conquer square root algorithm is largely
lifted from Bertot et al. [3], with only minor adjustments to account for small changes in the
GMP implementation since the publication of their article. For the sake of completeness,
the signature and specification of the algorithm are given in Alg. 18.

Algorithm 18 Square root of a normalized integer (specification)
Require: valid(a, 2n), valid(s, n), valid(w, bn/2c + 1)
Require: 1 ≤ n
Require: a[2n − 1] ≥ β/4
Ensure: value(s, n)2 + value(a, n) + β n result = value(old a, 2n)
Ensure: value(a, n) + β n result ≤ 2 × value(s, n)
Ensure: 0 ≤ result ≤ 1
  function DC _ SQRTREM(s, a, w, n)


6.4     Square root, normalizing wrapper
The previous “general case” algorithm still requires its input to have even length and to be
normalized (the highest limb of a must be greater than or equal to β/4). The final algorithm
(Alg. 19) is essentially a wrapper around the two previous ones that normalizes its operand,
calls the appropriate square root function, and denormalizes the result. It returns the size
of the remainder in limbs. If the result is 0, the operand is a perfect square.
   The key idea of this algorithm is the following lemma from Bertot et al. [3]:
  L EMMA NORMALIZATION . Let N, N1 , S1 , c such that S12 ≤ N1 < (S1 + 1)2 , N1 =
2 N . Let S, s0 such that S1 = 2c S + s0 , 0 ≤ s0 < 2c . Then S 2 ≤ N < (S + 1)2 .
 2c


This justifies that simply denormalizing the square root of the normalized operand gives
the correct square root.
   We first compute c, the floor of half the number of leading zeros of the high limb of a.
This means that shifting a to the left by 2c will multiply it by a power of 2 such that at
most one leading zero remains, which is exactly the precondition of the previous square
root algorithms.
   6.4.1 Special case n = 1. This case is a direct application of the lemma above. If the
operand is not already normalized, we shift it to the left by 2c and shift the result to the
right by c. The lemma ensures that this yields the correct square root, and we compute a
remainder straightforwardly.
   6.4.2 General case. We compute k = b(n + 1)/2c such that 2k = n + (n mod 2).
If the conditional at line 19 is true, then a is either not normalized or of odd length, so
we must normalize it before calling dc_sqrtrem. After line 27, this is done and we have
value(t, 2k) = 22c value(a, n). The value of c is incremented by 32 if n is odd, so that
we shift a by an extra limb (which is what the computation of t0 does) and end up with an
even-length number.
   At that point, we can call dc_sqrtrem. The normalization lemma implies that b value2(s,k)
                                                                                       c     c
is the correct square root (using the same variables as the lemma, S = bS1 2−c c). We still
need to compute the remainder.
Journal of Formalized Reasoning Vol.12, No.1, 2019
A Why3 proof of GMP algorithms                                  91


Algorithm 19 Square root of an integer
Require: valid(s, bn/2c + 1), valid(r, n), valid(a, n)
Require: 1 ≤ n
Require: a[n − 1] > 0
Ensure: value(a, n) = value(s, bn/2c + 1)2 + value(r, result)
Ensure: value(r, result) ≤ 2 × value(s, n)
Ensure: result > 0 ⇒ r[result − 1] > 0
 1: function SQRTREM (s, r, a, n)
 2:    h ← a[n − 1]
 3:    c = COUNT _ LEADING _ ZEROS(h)/2
 4:    if n = 1 then
 5:        if c = 0 then
 6:             s[0] = SQRT 1(r, h)
 7:        else
 8:             h0 ← h  2c
 9:             s0 ← SQRT 1(r, h0 )  c
10:             s[0] ← s0
11:             r[0] ← h − s02
12:        if r[0] = 0 then
13:             return 0
14:        else
15:             return 1
16:    k ← (n + 1)/2
17:    w ← ALLOC(k/2 + 1)
18:    t ← ALLOC(2k)
19:    if n = 1 mod 2 ∨ c 6= 0 then
20:        t[0] ← 0
21:        t0 ← t + (n mod 2)
22:        if c 6= 0 then
23:             LSHIFT (t0 , a, n, 2c)
24:        else
25:             COPY (t0 , a, n)
26:        if n mod 2 = 1 then
27:             c ← c + 32
28:        rl ← DC _ SQRTREM(s, t, w, k)
29:        s0 ← ALLOC(1)
30:        s0 [0] ← s[0] mod 2c
31:        rl ← rl + ADDMUL _1(t, s, k, 2s0 [0])
32:        b ← SUBMUL _1(t, s0 , s0 [0], 1)
33:        if k > 1 then
34:             b ← SUB _1_ IN _ PLACE(t + 1, b, k − 1)
35:        rl ← rl − b
36:        RSHIFT _ IN _ PLACE (s, k, c)
37:        t[k] ← rl
92       ·      Raphaël Rieu-Helft

38:           c2 ← 2c
39:           if c2 < 64 then
40:                k ←k+1
41:           else
42:                t←t+1
43:                c2 ← c2 − 64
44:           if c2 6= 0 then
45:                RSHIFT(r, t, k, c2 )
46:           else
47:                COPY (r, t, k)
48:           rn ← k
49:       else
50:           COPY (r, a, n)
51:           h ← DC _ SQRTREM(s, r, w, k)                                          .0≤h≤1
52:           r[k] ← h
53:           rn ← k + h
54:       while r[rn − 1] = 0 do
55:           rn ← rn − 1
56:           if rn = 0 then
57:                break
58:       return rn

   At line 29, we pose S1 = value(s, k), R1 = value(t, k) + β k rl , N = value(a, n) =,
N1 = 2c N . We have S12 + R1 = N1 . We pose s0 = S1 mod 2c . N1 can be written as
(S1 − s0 )2 + 2S1 s0 − s20 + R1 . We add 2S1 s0 to t at line 31, and then subtract s20 at line 32.
After line 32, value(t, k) + β k rl − βb = R1 + 2S1 s0 − s20 . Note that the borrow b has not
yet been propagated, instead the subtraction was only over a length of 1. The propagation
occurs at lines 33-34, and after line 35 we have value(t, k) + β k rl = R1 + 2S1 s0 − s20 .
At line 37 we write rl at t[k] such that value(t, k + 1) = R1 + 2S1 s0 − s20 .
   We pose S such that S1 = s0 + 2c S. The normalization lemma implies that S is the
square root of N . At line 36, we denormalize s such that value(s, k) = S. Furthermore,
S1 − s0 = 2c S, so t holds the remainder of the square root: value(t, k + 1) = R1 +
2S1 s0 − s20 = N1 − 22c S 2 = 22c (N − S 2 ). The only remaining thing to do is to shift t by
2c. As c can be up to 63, 2c can be up to 126 and rshift only accepts parameters smaller
than 64. The conditional at lines 39-43 takes care of this. If c2 ≥ 64, then the low limb of t
is all zeroes, so we can cheaply shift t by 64 by simply incrementing the pointer t. At line
44, value(t, k) = 2c2 (N − S 2 ) and 0 ≤ c2 ≤ 63, so we only have to shift t to the right to
get the correct remainder.
  6.4.3 Normalized case. If the operand is already normalized and of even length, we
can simply call dc_sqrtrem. This is done at lines 50-53.
   6.4.4 Normalizing the remainder. At line 54, rn contains an upper bound on the length
of the remainder. However, any number of high limbs may be zero. We simply iterate to
return rn such that rn = 0 or r[rn − 1] > 0.



Journal of Formalized Reasoning Vol.12, No.1, 2019
                                              A Why3 proof of GMP algorithms                  ·      93

7.   RELATED WORK
In this work, we have used the Why3 tool [4, 13, 11] to formally verify GMP’s integer
arithmetic layer. We obtain a verified and efficient C library. Previous work generally does
not deal with a large number of highly optimized algorithms. As far as we know, this work
is the first formal verification of a comprehensive arbitrary-precision integer library with
comparable performance to the state of the art. Let us now discuss a few examples of
existing verifications of multiprecision arithmetic functions or libraries.
   Bertot et al. verified GMP’s square-root general case algorithm [3] using Coq and the
Correctness tool, which translates an imperative program and its specifications into verifi-
cation conditions to be proved with Coq. Our Why3 proof of the same algorithm is directly
inspired from their article. They specify the memory as one large array of machine inte-
gers, so their specifications must include additional clauses to tell which memory zones
are left unchanged. Their formalization is otherwise quite similar to ours. However, their
proof is about 13,000 lines long, which is about as long as all our proofs combined. Indeed,
Why3 proofs are partially automatic, while Coq proofs are entirely interactive, so it is not
surprising that we enjoy a lower proof effort.
   Myreen and Curello verified an implementation of arbitrary-precision integer arithmetic
using the HOL4 theorem prover and separation logic [23]. Their library covers the four
basic arithmetic operations, but not the square root. They do not attempt to produce highly
efficient code. As a result, the algorithms they proved are simpler and less efficient than
the optimized ones that we proved. For example, their multiplication algorithm is the
schoolbook one. However, their verification goes all the way down to x86 machine code,
using formally verified proof-generating compilers and decompilers to do part of the proof
on a higher-level implementation. Using these tools, they also manage to avoid most proofs
involving pointer reasoning. Their proof effort per algorithm is similar to ours despite them
using an interactive tool.
   Affeldt used Coq to verify a binary extended GCD algorithm implemented in a variant
of MIPS assembly [1], as well as the functions it depends on, such as addition, subtraction
and halving. The work encompasses both signed and unsigned integer arithmetic. It uses
GMP’s number representation and a memory model based on separation logic. The author
verifies the algorithm in a pseudo-code language and proves a forward simulation relation
between the pseudo-code and the MIPS assembly code to prove the latter’s correctness. It
makes some simplifying assumptions, such as requiring that the operands share the same
length. The GCD algorithm that is verified is not trivial, but it is much less involved and
efficient than the one implemented in GMP (which we have not verified). The proof effort
is hard to quantify, as the development relies on preexisting frameworks by the same author
for pseudocode and assembly code, but it is rather high. The proofs of the algorithms
amount to about 15,000 lines of Coq.
   Fischer designed a modular exponentiation library [14] verified using Isabelle/HOL and
a framework for verifying imperative programs developed by Schirmer [25]. The verified
algorithms include multiplication, division and square-and-multiply modular exponentia-
tion. The library is not meant to be as efficient as GMP, as it represents arbitrary-precision
integers as garbage-collected doubly-linked lists of machine integers. The algorithms are
implemented in a restricted variant of the C language and are automatically transcribed into
Isabelle. The functional correctness of the algorithms and the pointer-level correctness of
the data structure are proved, but not termination or the absence of arithmetic overflows.

                                                       Journal of Formalized Reasoning Vol.12, No.1, 2019
94       ·      Raphaël Rieu-Helft

The author reports running into slowdown and memory issues inside the tool due to the
great number of invariants and conditions present in the logical context to keep track of
aliasing. By contrast, Why3 automatically keeps track of aliases inside its region-based
type system, rather than in the logic. This means that the user does not need to mention in
specifications and proofs that such and such pointers are not aliased, which would other-
wise cause large slowdown issues similar to those reported by Fischer.
   Berghofer developed a verified bignum library programmed in the S PARK fragment of
the Ada programming language, using a verification framework that sends goals to Is-
abelle/HOL [2]. The library provides modular exponentiation, as well as the primitives
required to implement it: modular multiplication and squaring, modular inverse, and basic
operations such as subtraction and doubling. A simple square-and-multiply algorithm is
used for modular exponentiation, without the sliding window optimization that GMP and
other state-of-the-art libraries implement. The author reports a 150% slowdown compared
to OpenSSL for their implementation of RSA using their library. However, OpenSSL uses
hand-written assembly code, which accounts for a large part of the discrepancy. The proof
effort for the library is only about 2,000 lines of Isabelle written over three weeks, which
is surprisingly low, even taking into account the low amount of verified algorithms.
   Unlike our semi-automatic Why3 proofs, most of the approaches described above use
interactive proof assistants. There have also been efforts to prove arithmetic libraries using
automated tools. Schoolderman used Why3 to verify hand-optimised Karatsuba multipli-
cation branch-free assembly routines for the AVR microcontroller architecture [26]. The
algorithms are not arbitrary-precision, instead there are many routines, each suited for a
particular operand size up to 96 × 96 bits. The fact that the size of the operands is fixed
and relatively small means the loops can be unrolled, which is why the algorithms are
branch-free. The size being known also makes the proof much easier for SMT solvers, and
the authors only needed to add a very small number of annotations to make the automatic
proof succeed.
   Finally, Zinzindohoué et al. developed an elliptic curve cryptography library written and
verified in F∗ that can be extracted to C [29]. It implements the full NaCl API. The integers
do not have an arbitrary large size, instead they have a small, fixed size that depends on the
choice of elliptic curve. A peculiarity is that only part of a machine word is used for each
limb. For instance, the most used curve uses 255-bit integers, which are not stored using
4 limbs of 64 bits, but 5 limbs of 51 bits. This means that arithmetic operations on limbs
do not overflow and carry bits do not have to be propagated. The performance of their
extracted C code is comparable to state-of-the-art C implementations. Part of the code is
now in use in several products such as the Mozilla Firefox web browser.

8.    CONCLUSION
We have used the Why3 platform to verify a library that implements state-of-the-art algo-
rithms from the mpn layer of GMP for comparison, addition, subtraction, multiplication,
division and square root. We have tackled two different algorithms for multiplication, one
suited for relatively small numbers (fewer than 1500 bits) and one suited for larger num-
bers. This verification effort covers functional correctnes and implementation aspects such
as memory safety and absence of arithmetic overflows.
   Due to the partially automatic nature of Why3 proofs, we are able to perform such a ver-
ification work with a high but reasonable amount of proof effort. The Why3 development
is about 13000 lines long, split between about 4500 lines of program code and 8500 lines
Journal of Formalized Reasoning Vol.12, No.1, 2019
                                             A Why3 proof of GMP algorithms                  ·      95

of annotations and proofs. The extracted library consists of about 2500 lines of C code.
The proof-to-code ratio is worse than for most existing Why3 proofs. Indeed, non-linear
arithmetic is difficult for automated provers, and the fact that the operands have unknown
lengths prevents the use of bit-blasting techniques. To make up for this, we had to write a
large amount of assertions that are practically transcriptions of paper proofs. As such, for
the mathematical parts of the proofs, the way we used Why3 was not much faster than the
way we would have used an interactive tool such as Coq.
   While the mathematical parts of the proof still require a lot of user guidance, most sub-
goals related to memory safety and arithmetic overflows can be discharged by automated
theorem provers. Our C memory model and Why3’s static treatment of aliasing make the
memory management part come almost for free. The GMP source code is rather well
documented, which made the specifications and invariants easy to write except for some
internal functions. All in all, the total proof effort is lower than for most related works
which use interactive proof assistants.
   Using Why3’s extraction mechanism, we obtain a verified C library that is compatible
with GMP. Our implementation contains all the algorithmic tricks that can be found in the
corresponding GMP functions, which makes the extracted code performance-competitive
with the non-assembly version of GMP [24]. What we have verified is not strictly GMP’s
implementation, but rather their WhyML transcription and, by extraction, a new C library
that uses the same algorithms. However, the fact that our implementation closely mirrors
GMP’s C code does give an additional measure of trust in GMP’s own correctness. While
we have not found any actual bug in GMP, we found something that can be seen as a bug
in the proof of GMP’s Toom-2 algorithm. Indeed, although it turned out that it was correct,
the complexity of our proof led the GMP developers to modify the algorithm to match the
correctness proof that they expected.
   Our library does not implement GMP’s full interface. We intend to verify the missing
primitives, e.g. side-channel resistant modular exponentiation and base conversions, in
future work. Another natural extension would be to tackle the mpz layer of GMP, which
is the higher-level interface for signed arbitrary-precision arithmetic. It is essentially a
wrapper around the mpn layer that we tackled in this work. It contains relatively little
arithmetic but poses many interesting verification challenges in its memory management.

Acknowledgements
We gratefully thank Guillaume Melquiond and the anonymous reviewers for their extensive
feedback on earlier versions of this article.

References
 [1] Reynald Affeldt. On construction of a library of formally verified low-level arithmetic
     functions. Innovations in Systems and Software Engineering, 9(2):59–77, 2013.
 [2] Stefan Berghofer. Verification of dependable software using SPARK and Isabelle. In
     Jörg Brauer, Marco Roveri, and Hendrik Tews, editors, 6th International Workshop
     on Systems Software Verification, volume 24 of OpenAccess Series in Informatics
     (OASIcs), pages 15–31, Dagstuhl, Germany, 2012.
 [3] Yves Bertot, Nicolas Magaud, and Paul Zimmermann. A proof of GMP square root.
     Journal of Automated Reasoning, 29(3-4):225–252, 2002.


                                                      Journal of Formalized Reasoning Vol.12, No.1, 2019
96       ·      Raphaël Rieu-Helft

 [4] François Bobot, Jean-Christophe Filliâtre, Claude Marché, and Andrei Paskevich.
     Why3: Shepherd your herd of provers. In Boogie 2011: First International Work-
     shop on Intermediate Verification Languages, pages 53–64, Wrocław, Poland, August
     2011. https://hal.inria.fr/hal-00790310.
 [5] Marco Bodrato and Alberto Zanoni. Integer and polynomial multiplication: Towards
     optimal Toom-Cook matrices. In 2007 International Symposium on Symbolic and
     Algebraic Computation, pages 17–24. ACM, 2007.
 [6] Richard P. Brent and Paul Zimmermann. Modern Computer Arithmetic. Cambridge
     University Press, 2010.
 [7] Martin Clochard. Preuves taillées en biseau. In vingt-huitièmes Journées Franco-
     phones des Langages Applicatifs (JFLA), Gourette, France, January 2017.
 [8] Stephen A. Cook. On the minimum computation time of functions. PhD thesis, De-
     partment of Mathematics, Harvard University, 1966.
 [9] Marc Daumas and Guillaume Melquiond. Certification of bounds on expressions
     involving rounded operators. Transactions on Mathematical Software, 37(1):1–20,
     2010.
[10] Edsger W. Dijkstra. A discipline of programming, volume 1. Prentice-Hall Engle-
     wood Cliffs, 1976.
[11] Jean-Christophe Filliâtre. One logic to use them all. In 24th International Confer-
     ence on Automated Deduction (CADE-24), volume 7898 of Lecture Notes in Artificial
     Intelligence, pages 1–20, Lake Placid, USA, June 2013.
[12] Jean-Christophe Filliâtre, Léon Gondelman, and Andrei Paskevich. The spirit of
     ghost code. Formal Methods in System Design, 48(3):152–174, 2016.
[13] Jean-Christophe Filliâtre and Andrei Paskevich. Why3 — where programs meet
     provers. In Matthias Felleisen and Philippa Gardner, editors, 22nd European Sym-
     posium on Programming, volume 7792 of Lecture Notes in Computer Science, pages
     125–128, Heidelberg, Germany, March 2013.
[14] Sabine Fischer. Formal verification of a big integer library. In DATE Workshop on
     Dependable Software Systems, 2008.
[15] Robert W. Floyd. Assigning meanings to programs. In Program Verification, pages
     65–81. Springer, 1993.
[16] Léon Gondelman. A Pragmatic Type System for Deductive Software Verification.
     PhD thesis, Université Paris-Saclay, December 2016.
[17] C. A. R. Hoare. An axiomatic basis for computer programming. Communications of
     the ACM, 12(10):576–580, 1969.
[18] Anatolii Karatsuba. Multiplication of multidigit numbers on automata. In Soviet
     Physics Doklady, volume 7, pages 595–596, 1963.
[19] Donald E. Knuth. The Art of Computer Programming, Volume 2 (3rd Ed.): Seminu-
     merical Algorithms. Addison-Wesley Longman Publishing Co., Inc., Boston, MA,
     USA, 1997.
[20] Guillaume Melquiond and Raphaël Rieu-Helft. A Why3 framework for reflection
     proofs and its application to GMP’s algorithms. In Didier Galmiche, Stephan Schulz,
     and Roberto Sebastiani, editors, 9th International Joint Conference on Automated
     Reasoning, number 10900 in Lecture Notes in Computer Science, pages 178–193,
     Oxford, United Kingdom, July 2018.
Journal of Formalized Reasoning Vol.12, No.1, 2019
                                             A Why3 proof of GMP algorithms                 ·      97

[21] Guillaume Melquiond and Raphaël Rieu-Helft. Formal verification of a state-of-the-
     art integer square root. In IEEE 26th Symposium on Computer Arithmetic (ARITH),
     Kyoto, Japan, June 2019.
[22] Niels Moller and Torbjörn Granlund. Improved division by invariant integers. IEEE
     Transactions on Computers, 60:165–175, 2011.
[23] Magnus O. Myreen and Gregorio Curello. Proof pearl: A verified bignum implemen-
     tation in x86-64 machine code. In Georges Gonthier and Michael Norrish, editors,
     3rd International Conference on Certified Programs and Proofs (CPP), volume 8307
     of Lecture Notes in Computer Science, pages 66–81, Melbourne, Australia, Decem-
     ber 2013.
[24] Raphaël Rieu-Helft, Claude Marché, and Guillaume Melquiond. How to get an ef-
     ficient yet verified arbitrary-precision integer library. In 9th Working Conference on
     Verified Software: Theories, Tools, and Experiments, volume 10712 of Lecture Notes
     in Computer Science, pages 84–101, Heidelberg, Germany, July 2017.
[25] Norbert Schirmer. A verification environment for sequential imperative programs
     in Isabelle/HOL. In International Conference on Logic for Programming Artificial
     Intelligence and Reasoning, pages 398–414, 2005.
[26] Marc Schoolderman. Verifying branch-free assembly code in Why3. In Working
     Conference on Verified Software: Theories, Tools, and Experiments, pages 66–83,
     2017.
[27] Andrei L. Toom. The complexity of a scheme of functional elements realizing the
     multiplication of integers. In Soviet Mathematics Doklady, volume 3, pages 714–716,
     1963.
[28] Henry S. Warren. Hacker’s Delight. Addison-Wesley Professional, 2nd edition, 2012.
[29] Jean Karim Zinzindohoué, Karthikeyan Bhargavan, Jonathan Protzenko, and Ben-
     jamin Beurdouche. HACL*: A verified modern cryptographic library. Cryptology
     ePrint Archive, Report 2017/536, 2017. https://eprint.iacr.org/2017/536.




                                                     Journal of Formalized Reasoning Vol.12, No.1, 2019
