[Ada] Prove double precision integer arithmetic unit

gcc/ada/

	* libgnat/a-nbnbig.ads: Mark the unit as Pure.
	* libgnat/s-aridou.adb: Add contracts and ghost code for proof.
	(Scaled_Divide): Reorder operations and use of temporaries in
	two places to facilitate proof.
	* libgnat/s-aridou.ads: Add full functional contracts.
	* libgnat/s-arit64.adb: Mark in SPARK.
	* libgnat/s-arit64.ads: Add contracts similar to those from
	s-aridou.ads.
	* rtsfind.ads: Document the limitation that runtime units
	loading does not work for private with-clauses.

6 files changed, 2515 insertions, 113 deletions

diff --git a/gcc/ada/libgnat/a-nbnbig.ads b/gcc/ada/libgnat/a-nbnbig.ads
index 237bca1..f574e78 100644
--- a/gcc/ada/libgnat/a-nbnbig.ads
+++ b/gcc/ada/libgnat/a-nbnbig.ads
@@ -30,7 +30,7 @@ pragma Assertion_Policy (Ghost => Ignore);
 package Ada.Numerics.Big_Numbers.Big_Integers_Ghost with
    SPARK_Mode,
    Ghost,
-   Preelaborate
+   Pure
 is
    type Big_Integer is private
      with Integer_Literal => From_Universal_Image;
diff --git a/gcc/ada/libgnat/s-aridou.adb b/gcc/ada/libgnat/s-aridou.adb
index b47c319..67f2440 100644
--- a/gcc/ada/libgnat/s-aridou.adb
+++ b/gcc/ada/libgnat/s-aridou.adb
@@ -31,7 +31,20 @@
 
 with Ada.Unchecked_Conversion;
 
-package body System.Arith_Double is
+package body System.Arith_Double
+  with SPARK_Mode
+is
+   --  Contracts, ghost code, loop invariants and assertions in this unit are
+   --  meant for analysis only, not for run-time checking, as it would be too
+   --  costly otherwise. This is enforced by setting the assertion policy to
+   --  Ignore.
+
+   pragma Assertion_Policy (Pre            => Ignore,
+                            Post           => Ignore,
+                            Contract_Cases => Ignore,
+                            Ghost          => Ignore,
+                            Loop_Invariant => Ignore,
+                            Assert         => Ignore);
 
    pragma Suppress (Overflow_Check);
    pragma Suppress (Range_Check);
@@ -42,6 +55,30 @@ package body System.Arith_Double is
    Double_Size : constant Natural := Double_Int'Size;
    Single_Size : constant Natural := Double_Int'Size / 2;
 
+   --  Power-of-two constants. Use the names Big_2xx32, Big_2xx63 and Big_2xx64
+   --  even if Single_Size might not be 32 and Double_Size might not be 64, as
+   --  this facilitates code and proof understanding, compared to more generic
+   --  names.
+
+   pragma Warnings
+     (Off, "non-preelaborable call not allowed in preelaborated unit",
+      Reason => "Ghost code is not compiled");
+   pragma Warnings
+     (Off, "non-static constant in preelaborated unit",
+      Reason => "Ghost code is not compiled");
+   Big_2xx32 : constant Big_Integer :=
+     Big (Double_Int'(2 ** Single_Size))
+   with Ghost;
+   Big_2xx63 : constant Big_Integer :=
+     Big (Double_Uns'(2 ** (Double_Size - 1)))
+   with Ghost;
+   Big_2xx64 : constant Big_Integer :=
+     Big (Double_Uns'(2 ** Double_Size - 1)) + 1
+   with Ghost;
+   pragma Warnings
+     (On, "non-preelaborable call not allowed in preelaborated unit");
+   pragma Warnings (On, "non-static constant in preelaborated unit");
+
    -----------------------
    -- Local Subprograms --
    -----------------------
@@ -57,7 +94,9 @@ package body System.Arith_Double is
    --  Length doubling multiplication
 
    function "/" (A : Double_Uns; B : Single_Uns) return Double_Uns is
-     (A / Double_Uns (B));
+     (A / Double_Uns (B))
+   with
+     Pre => B /= 0;
    --  Length doubling division
 
    function "&" (Hi, Lo : Single_Uns) return Double_Uns is
@@ -66,16 +105,34 @@ package body System.Arith_Double is
 
    function "abs" (X : Double_Int) return Double_Uns is
      (if X = Double_Int'First
-      then 2 ** (Double_Size - 1)
+      then Double_Uns'(2 ** (Double_Size - 1))
       else Double_Uns (Double_Int'(abs X)));
    --  Convert absolute value of X to unsigned. Note that we can't just use
    --  the expression of the Else since it overflows for X = Double_Int'First.
 
    function "rem" (A : Double_Uns; B : Single_Uns) return Double_Uns is
-     (A rem Double_Uns (B));
+     (A rem Double_Uns (B))
+   with
+     Pre => B /= 0;
    --  Length doubling remainder
 
-   function Le3 (X1, X2, X3, Y1, Y2, Y3 : Single_Uns) return Boolean;
+   function Big_2xx (N : Natural) return Big_Integer is
+     (Big (Double_Uns'(2 ** N)))
+   with
+     Ghost,
+     Pre  => N < Double_Size,
+     Post => Big_2xx'Result > 0;
+   --  2**N as a big integer
+
+   function Big3 (X1, X2, X3 : Single_Uns) return Big_Integer is
+     (Big_2xx32 * Big_2xx32 * Big (Double_Uns (X1))
+                + Big_2xx32 * Big (Double_Uns (X2)) + Big (Double_Uns (X3)))
+   with Ghost;
+   --  X1&X2&X3 as a big integer
+
+   function Le3 (X1, X2, X3, Y1, Y2, Y3 : Single_Uns) return Boolean
+   with
+     Post => Le3'Result = (Big3 (X1, X2, X3) <= Big3 (Y1, Y2, Y3));
    --  Determines if (3 * Single_Size)-bit value X1&X2&X3 <= Y1&Y2&Y3
 
    function Lo (A : Double_Uns) return Single_Uns is
@@ -86,16 +143,27 @@ package body System.Arith_Double is
      (Single_Uns (Shift_Right (A, Single_Size)));
    --  High order half of double value
 
-   procedure Sub3 (X1, X2, X3 : in out Single_Uns; Y1, Y2, Y3 : Single_Uns);
+   procedure Sub3 (X1, X2, X3 : in out Single_Uns; Y1, Y2, Y3 : Single_Uns)
+   with
+     Pre  => Big3 (X1, X2, X3) >= Big3 (Y1, Y2, Y3),
+     Post => Big3 (X1, X2, X3) = Big3 (X1, X2, X3)'Old - Big3 (Y1, Y2, Y3);
    --  Computes X1&X2&X3 := X1&X2&X3 - Y1&Y1&Y3 mod 2 ** (3 * Single_Size)
 
-   function To_Neg_Int (A : Double_Uns) return Double_Int;
+   function To_Neg_Int (A : Double_Uns) return Double_Int
+   with
+     Annotate => (GNATprove, Terminating),
+     Pre      => In_Double_Int_Range (-Big (A)),
+     Post     => Big (To_Neg_Int'Result) = -Big (A);
    --  Convert to negative integer equivalent. If the input is in the range
    --  0 .. 2 ** (Double_Size - 1), then the corresponding nonpositive signed
    --  integer (obtained by negating the given value) is returned, otherwise
    --  constraint error is raised.
 
-   function To_Pos_Int (A : Double_Uns) return Double_Int;
+   function To_Pos_Int (A : Double_Uns) return Double_Int
+   with
+     Annotate => (GNATprove, Terminating),
+     Pre      => In_Double_Int_Range (Big (A)),
+     Post     => Big (To_Pos_Int'Result) = Big (A);
    --  Convert to positive integer equivalent. If the input is in the range
    --  0 .. 2 ** (Double_Size - 1) - 1, then the corresponding non-negative
    --  signed integer is returned, otherwise constraint error is raised.
@@ -104,6 +172,396 @@ package body System.Arith_Double is
    pragma No_Return (Raise_Error);
    --  Raise constraint error with appropriate message
 
+   ------------------
+   -- Local Lemmas --
+   ------------------
+
+   procedure Inline_Le3 (X1, X2, X3, Y1, Y2, Y3 : Single_Uns)
+   with
+     Ghost,
+     Pre  => Le3 (X1, X2, X3, Y1, Y2, Y3),
+     Post => Big3 (X1, X2, X3) <= Big3 (Y1, Y2, Y3);
+
+   procedure Lemma_Abs_Commutation (X : Double_Int)
+   with
+     Ghost,
+     Post => abs (Big (X)) = Big (Double_Uns'(abs X));
+
+   procedure Lemma_Abs_Div_Commutation (X, Y : Big_Integer)
+   with
+     Ghost,
+     Pre  => Y /= 0,
+     Post => abs (X / Y) = abs X / abs Y;
+
+   procedure Lemma_Abs_Mult_Commutation (X, Y : Big_Integer)
+   with
+     Ghost,
+     Post => abs (X * Y) = abs X * abs Y;
+
+   procedure Lemma_Abs_Rem_Commutation (X, Y : Big_Integer)
+   with
+     Ghost,
+     Pre  => Y /= 0,
+     Post => abs (X rem Y) = (abs X) rem (abs Y);
+
+   procedure Lemma_Add_Commutation (X : Double_Uns; Y : Single_Uns)
+   with
+     Ghost,
+     Pre  => X <= 2 ** Double_Size - 2 ** Single_Size,
+     Post => Big (X) + Big (Double_Uns (Y)) = Big (X + Double_Uns (Y));
+
+   procedure Lemma_Add_One (X : Double_Uns)
+   with
+     Ghost,
+     Pre  => X /= Double_Uns'Last,
+     Post => Big (X + Double_Uns'(1)) = Big (X) + 1;
+
+   procedure Lemma_Bounded_Powers_Of_2_Increasing (M, N : Natural)
+   with
+     Ghost,
+     Pre  => M < N and then N < Double_Size,
+     Post => Double_Uns'(2)**M < Double_Uns'(2)**N;
+
+   procedure Lemma_Deep_Mult_Commutation
+     (Factor : Big_Integer;
+      X, Y   : Single_Uns)
+   with
+     Ghost,
+     Post =>
+       Factor * Big (Double_Uns (X)) * Big (Double_Uns (Y)) =
+         Factor * Big (Double_Uns (X) * Double_Uns (Y));
+
+   procedure Lemma_Div_Commutation (X, Y : Double_Uns)
+   with
+     Ghost,
+     Pre  => Y /= 0,
+     Post => Big (X) / Big (Y) = Big (X / Y);
+
+   procedure Lemma_Div_Definition
+     (A : Double_Uns;
+      B : Single_Uns;
+      Q : Double_Uns;
+      R : Double_Uns)
+   with
+     Ghost,
+     Pre  => B /= 0 and then Q = A / B and then R = A rem B,
+     Post => Big (A) = Big (Double_Uns (B)) * Big (Q) + Big (R);
+
+   procedure Lemma_Div_Ge (X, Y, Z : Big_Integer)
+   with
+     Ghost,
+     Pre  => Z > 0 and then X >= Y * Z,
+     Post => X / Z >= Y;
+
+   procedure Lemma_Div_Lt (X, Y, Z : Big_Natural)
+   with
+     Ghost,
+     Pre  => Z > 0 and then X < Y * Z,
+     Post => X / Z < Y;
+
+   procedure Lemma_Div_Eq (A, B, S, R : Big_Integer)
+   with
+     Ghost,
+     Pre  => A * S = B * S + R and then S /= 0,
+     Post => A = B + R / S;
+
+   procedure Lemma_Double_Shift (X : Double_Uns; S, S1 : Double_Uns)
+   with
+     Ghost,
+     Pre  => S <= Double_Uns (Double_Size)
+       and then S1 <= Double_Uns (Double_Size),
+     Post => Shift_Left (Shift_Left (X, Natural (S)), Natural (S1)) =
+             Shift_Left (X, Natural (S + S1));
+
+   procedure Lemma_Double_Shift (X : Single_Uns; S, S1 : Natural)
+   with
+     Ghost,
+     Pre  => S <= Single_Size - S1,
+     Post => Shift_Left (Shift_Left (X, S), S1) = Shift_Left (X, S + S1);
+
+   procedure Lemma_Double_Shift (X : Double_Uns; S, S1 : Natural)
+   with
+     Ghost,
+     Pre  => S <= Double_Size - S1,
+     Post => Shift_Left (Shift_Left (X, S), S1) = Shift_Left (X, S + S1);
+
+   procedure Lemma_Double_Shift_Right (X : Double_Uns; S, S1 : Double_Uns)
+   with
+     Ghost,
+     Pre  => S <= Double_Uns (Double_Size)
+       and then S1 <= Double_Uns (Double_Size),
+     Post => Shift_Right (Shift_Right (X, Natural (S)), Natural (S1)) =
+             Shift_Right (X, Natural (S + S1));
+
+   procedure Lemma_Double_Shift_Right (X : Double_Uns; S, S1 : Natural)
+   with
+     Ghost,
+     Pre  => S <= Double_Size - S1,
+     Post => Shift_Right (Shift_Right (X, S), S1) = Shift_Right (X, S + S1);
+
+   procedure Lemma_Ge_Commutation (A, B : Double_Uns)
+   with
+     Ghost,
+     Pre  => A >= B,
+     Post => Big (A) >= Big (B);
+
+   procedure Lemma_Ge_Mult (A, B, C, D : Big_Integer)
+   with
+     Ghost,
+     Pre  => A >= B and then B * C >= D and then C > 0,
+     Post => A * C >= D;
+
+   procedure Lemma_Gt_Commutation (A, B : Double_Uns)
+   with
+     Ghost,
+     Pre  => A > B,
+     Post => Big (A) > Big (B);
+
+   procedure Lemma_Gt_Mult (A, B, C, D : Big_Integer)
+   with
+     Ghost,
+     Pre  => A >= B and then B * C > D and then C > 0,
+     Post => A * C > D;
+
+   procedure Lemma_Hi_Lo (Xu : Double_Uns; Xhi, Xlo : Single_Uns)
+   with
+     Ghost,
+     Pre  => Xhi = Hi (Xu) and Xlo = Lo (Xu),
+     Post =>
+       Big (Xu) = Big_2xx32 * Big (Double_Uns (Xhi)) + Big (Double_Uns (Xlo));
+
+   procedure Lemma_Hi_Lo_3 (Xu : Double_Uns; Xhi, Xlo : Single_Uns)
+   with
+     Ghost,
+     Pre  => Xhi = Hi (Xu) and then Xlo = Lo (Xu),
+     Post => Big (Xu) = Big3 (0, Xhi, Xlo);
+
+   procedure Lemma_Lo_Is_Ident (X : Double_Uns)
+   with
+     Ghost,
+     Pre  => Big (X) < Big_2xx32,
+     Post => Double_Uns (Lo (X)) = X;
+
+   procedure Lemma_Lt_Commutation (A, B : Double_Uns)
+   with
+     Ghost,
+     Pre  => A < B,
+     Post => Big (A) < Big (B);
+
+   procedure Lemma_Lt_Mult (A, B, C, D : Big_Integer)
+   with
+     Ghost,
+     Pre  => A < B and then B * C <= D and then C > 0,
+     Post => A * C < D;
+
+   procedure Lemma_Mult_Commutation (X, Y : Single_Uns)
+   with
+     Ghost,
+     Post =>
+       Big (Double_Uns (X)) * Big (Double_Uns (Y)) =
+         Big (Double_Uns (X) * Double_Uns (Y));
+
+   procedure Lemma_Mult_Commutation (X, Y : Double_Int)
+   with
+     Ghost,
+     Pre  => In_Double_Int_Range (Big (X) * Big (Y)),
+     Post => Big (X) * Big (Y) = Big (X * Y);
+
+   procedure Lemma_Mult_Commutation (X, Y, Z : Double_Uns)
+   with
+     Ghost,
+     Pre  => Big (X) * Big (Y) < Big_2xx64 and then Z = X * Y,
+     Post => Big (X) * Big (Y) = Big (Z);
+
+   procedure Lemma_Mult_Decomposition
+     (Mult               : Big_Integer;
+      Xu, Yu             : Double_Uns;
+      Xhi, Xlo, Yhi, Ylo : Single_Uns)
+   with
+     Ghost,
+     Pre  => Mult = Big (Xu) * Big (Yu)
+       and then Xhi = Hi (Xu)
+       and then Xlo = Lo (Xu)
+       and then Yhi = Hi (Yu)
+       and then Ylo = Lo (Yu),
+     Post => Mult =
+               Big_2xx32 * Big_2xx32 * (Big (Double_Uns'(Xhi * Yhi)))
+                         + Big_2xx32 * (Big (Double_Uns'(Xhi * Ylo)))
+                         + Big_2xx32 * (Big (Double_Uns'(Xlo * Yhi)))
+                                     + (Big (Double_Uns'(Xlo * Ylo)));
+
+   procedure Lemma_Mult_Distribution (X, Y, Z : Big_Integer)
+   with
+     Ghost,
+     Post => X * (Y + Z) = X * Y + X * Z;
+
+   procedure Lemma_Neg_Div (X, Y : Big_Integer)
+   with
+     Ghost,
+     Pre  => Y /= 0,
+     Post => X / Y = (-X) / (-Y);
+
+   procedure Lemma_Neg_Rem (X, Y : Big_Integer)
+   with
+     Ghost,
+     Pre  => Y /= 0,
+     Post => X rem Y = X rem (-Y);
+
+   procedure Lemma_Powers_Of_2 (M, N : Natural)
+   with
+     Ghost,
+     Pre  => M < Double_Size
+       and then N < Double_Size
+       and then M + N <= Double_Size,
+     Post =>
+       Big_2xx (M) * Big_2xx (N) =
+         (if M + N = Double_Size then Big_2xx64 else Big_2xx (M + N));
+
+   procedure Lemma_Powers_Of_2_Commutation (M : Natural)
+   with
+     Ghost,
+     Subprogram_Variant => (Decreases => M),
+     Pre  => M <= Double_Size,
+     Post => Big (Double_Uns'(2))**M =
+              (if M < Double_Size then Big_2xx (M) else Big_2xx64);
+
+   procedure Lemma_Powers_Of_2_Increasing (M, N : Natural)
+   with
+     Ghost,
+     Subprogram_Variant => (Increases => M),
+     Pre  => M < N,
+     Post => Big (Double_Uns'(2))**M < Big (Double_Uns'(2))**N;
+
+   procedure Lemma_Rem_Abs (X, Y : Big_Integer)
+   with
+     Ghost,
+     Pre  => Y /= 0,
+     Post => X rem Y = X rem (abs Y);
+   pragma Unreferenced (Lemma_Rem_Abs);
+
+   procedure Lemma_Rem_Commutation (X, Y : Double_Uns)
+   with
+     Ghost,
+     Pre  => Y /= 0,
+     Post => Big (X) rem Big (Y) = Big (X rem Y);
+
+   procedure Lemma_Rem_Is_Ident (X, Y : Big_Integer)
+   with
+     Ghost,
+     Pre  => abs X < abs Y,
+     Post => X rem Y = X;
+   pragma Unreferenced (Lemma_Rem_Is_Ident);
+
+   procedure Lemma_Rem_Sign (X, Y : Big_Integer)
+   with
+     Ghost,
+     Pre  => Y /= 0,
+     Post => Same_Sign (X rem Y, X);
+   pragma Unreferenced (Lemma_Rem_Sign);
+
+   procedure Lemma_Rev_Div_Definition (A, B, Q, R : Big_Natural)
+   with
+     Ghost,
+     Pre  => A = B * Q + R and then R < B,
+     Post => Q = A / B and then R = A rem B;
+
+   procedure Lemma_Shift_Right (X : Double_Uns; Shift : Natural)
+   with
+     Ghost,
+     Pre  => Shift < Double_Size,
+     Post => Big (Shift_Right (X, Shift)) = Big (X) / Big_2xx (Shift);
+
+   procedure Lemma_Shift_Without_Drop
+     (X, Y  : Double_Uns;
+      Mask  : Single_Uns;
+      Shift : Natural)
+   with
+     Ghost,
+     Pre  => (Hi (X) and Mask) = 0  --  X has the first Shift bits off
+       and then Shift <= Single_Size
+       and then Mask = Shift_Left (Single_Uns'Last, Single_Size - Shift)
+       and then Y = Shift_Left (X, Shift),
+     Post => Big (Y) = Big_2xx (Shift) * Big (X);
+
+   procedure Lemma_Simplify (X, Y : Big_Integer)
+   with
+     Ghost,
+     Pre  => Y /= 0,
+     Post => X * Y / Y = X;
+
+   procedure Lemma_Substitution (A, B, C, C1, D : Big_Integer)
+   with
+     Ghost,
+     Pre  => C = C1 and then A = B * C + D,
+     Post => A = B * C1 + D;
+
+   procedure Lemma_Subtract_Commutation (X, Y : Double_Uns)
+   with
+     Ghost,
+     Pre  => X >= Y,
+     Post => Big (X) - Big (Y) = Big (X - Y);
+
+   procedure Lemma_Subtract_Double_Uns (X, Y : Double_Int)
+   with
+     Ghost,
+     Pre  => X >= 0 and then X <= Y,
+     Post => Double_Uns (Y - X) = Double_Uns (Y) - Double_Uns (X);
+
+   procedure Lemma_Word_Commutation (X : Single_Uns)
+   with
+     Ghost,
+     Post => Big_2xx32 * Big (Double_Uns (X)) = Big (2**32 * Double_Uns (X));
+
+   -----------------------------
+   -- Local lemma null bodies --
+   -----------------------------
+
+   procedure Inline_Le3 (X1, X2, X3, Y1, Y2, Y3 : Single_Uns) is null;
+   procedure Lemma_Abs_Commutation (X : Double_Int) is null;
+   procedure Lemma_Abs_Mult_Commutation (X, Y : Big_Integer) is null;
+   procedure Lemma_Add_Commutation (X : Double_Uns; Y : Single_Uns) is null;
+   procedure Lemma_Add_One (X : Double_Uns) is null;
+   procedure Lemma_Bounded_Powers_Of_2_Increasing (M, N : Natural) is null;
+   procedure Lemma_Deep_Mult_Commutation
+     (Factor : Big_Integer;
+      X, Y   : Single_Uns)
+   is null;
+   procedure Lemma_Div_Commutation (X, Y : Double_Uns) is null;
+   procedure Lemma_Div_Definition
+     (A : Double_Uns;
+      B : Single_Uns;
+      Q : Double_Uns;
+      R : Double_Uns)
+   is null;
+   procedure Lemma_Div_Ge (X, Y, Z : Big_Integer) is null;
+   procedure Lemma_Div_Lt (X, Y, Z : Big_Natural) is null;
+   procedure Lemma_Div_Eq (A, B, S, R : Big_Integer) is null;
+   procedure Lemma_Double_Shift (X : Double_Uns; S, S1 : Double_Uns) is null;
+   procedure Lemma_Double_Shift (X : Single_Uns; S, S1 : Natural) is null;
+   procedure Lemma_Double_Shift_Right (X : Double_Uns; S, S1 : Double_Uns)
+   is null;
+   procedure Lemma_Ge_Commutation (A, B : Double_Uns) is null;
+   procedure Lemma_Ge_Mult (A, B, C, D : Big_Integer) is null;
+   procedure Lemma_Gt_Commutation (A, B : Double_Uns) is null;
+   procedure Lemma_Gt_Mult (A, B, C, D : Big_Integer) is null;
+   procedure Lemma_Lo_Is_Ident (X : Double_Uns) is null;
+   procedure Lemma_Lt_Commutation (A, B : Double_Uns) is null;
+   procedure Lemma_Lt_Mult (A, B, C, D : Big_Integer) is null;
+   procedure Lemma_Mult_Commutation (X, Y : Single_Uns) is null;
+   procedure Lemma_Mult_Commutation (X, Y : Double_Int) is null;
+   procedure Lemma_Mult_Commutation (X, Y, Z : Double_Uns) is null;
+   procedure Lemma_Mult_Distribution (X, Y, Z : Big_Integer) is null;
+   procedure Lemma_Neg_Rem (X, Y : Big_Integer) is null;
+   procedure Lemma_Rem_Commutation (X, Y : Double_Uns) is null;
+   procedure Lemma_Rem_Is_Ident (X, Y : Big_Integer) is null;
+   procedure Lemma_Rem_Sign (X, Y : Big_Integer) is null;
+   procedure Lemma_Rev_Div_Definition (A, B, Q, R : Big_Natural) is null;
+   procedure Lemma_Simplify (X, Y : Big_Integer) is null;
+   procedure Lemma_Substitution (A, B, C, C1, D : Big_Integer) is null;
+   procedure Lemma_Subtract_Commutation (X, Y : Double_Uns) is null;
+   procedure Lemma_Subtract_Double_Uns (X, Y : Double_Int) is null;
+   procedure Lemma_Word_Commutation (X : Single_Uns) is null;
+
    --------------------------
    -- Add_With_Ovflo_Check --
    --------------------------
@@ -111,18 +569,124 @@ package body System.Arith_Double is
    function Add_With_Ovflo_Check (X, Y : Double_Int) return Double_Int is
       R : constant Double_Int := To_Int (To_Uns (X) + To_Uns (Y));
 
+      --  Local lemmas
+
+      procedure Prove_Negative_X
+      with
+        Ghost,
+        Pre  => X < 0 and then (Y > 0 or else R < 0),
+        Post => R = X + Y;
+
+      procedure Prove_Non_Negative_X
+      with
+        Ghost,
+        Pre  => X >= 0 and then (Y < 0 or else R >= 0),
+        Post => R = X + Y;
+
+      procedure Prove_Overflow_Case
+      with
+        Ghost,
+        Pre  =>
+          (if X >= 0 then Y >= 0 and then R < 0
+                     else Y <= 0 and then R >= 0),
+        Post => not In_Double_Int_Range (Big (X) + Big (Y));
+
+      ----------------------
+      -- Prove_Negative_X --
+      ----------------------
+
+      procedure Prove_Negative_X is
+      begin
+         if X = Double_Int'First then
+            if Y > 0 then
+               null;
+            else
+               pragma Assert
+                 (To_Uns (X) + To_Uns (Y) =
+                    2 ** (Double_Size - 1) - Double_Uns (-Y));
+               pragma Assert  --  as R < 0
+                 (To_Uns (X) + To_Uns (Y) >= 2 ** (Double_Size - 1));
+               pragma Assert (Y = 0);
+            end if;
+
+         elsif Y = Double_Int'First then
+            pragma Assert
+              (To_Uns (X) + To_Uns (Y) =
+                 2 ** (Double_Size - 1) - Double_Uns (-X));
+            pragma Assert (False);
+
+         elsif Y <= 0 then
+            pragma Assert
+              (To_Uns (X) + To_Uns (Y) = -Double_Uns (-X) - Double_Uns (-Y));
+
+         else  --  Y > 0, 0 > X > Double_Int'First
+            declare
+               Ru : constant Double_Uns := To_Uns (X) + To_Uns (Y);
+            begin
+               pragma Assert (Ru = -Double_Uns (-X) + Double_Uns (Y));
+               if Ru < 2 ** (Double_Size - 1) then  --  R >= 0
+                  Lemma_Subtract_Double_Uns (-X, Y);
+                  pragma Assert (Ru = Double_Uns (X + Y));
+
+               elsif Ru = 2 ** (Double_Size - 1) then
+                  pragma Assert (Double_Uns (Y) < 2 ** (Double_Size - 1));
+                  pragma Assert (Double_Uns (-X) < 2 ** (Double_Size - 1));
+                  pragma Assert (False);
+
+               else
+                  pragma Assert
+                    (R = -Double_Int (-(-Double_Uns (-X) + Double_Uns (Y))));
+                  pragma Assert
+                    (R = -Double_Int (-Double_Uns (Y) + Double_Uns (-X)));
+               end if;
+            end;
+         end if;
+      end Prove_Negative_X;
+
+      --------------------------
+      -- Prove_Non_Negative_X --
+      --------------------------
+
+      procedure Prove_Non_Negative_X is
+      begin
+         if Y >= 0 or else Y = Double_Int'First then
+            null;
+         else
+            pragma Assert
+              (To_Uns (X) + To_Uns (Y) = Double_Uns (X) - Double_Uns (-Y));
+         end if;
+      end Prove_Non_Negative_X;
+
+      -------------------------
+      -- Prove_Overflow_Case --
+      -------------------------
+
+      procedure Prove_Overflow_Case is
+      begin
+         if X < 0 and then X /= Double_Int'First and then Y /= Double_Int'First
+         then
+            pragma Assert
+              (To_Uns (X) + To_Uns (Y) = -Double_Uns (-X) - Double_Uns (-Y));
+         end if;
+      end Prove_Overflow_Case;
+
+   --  Start of processing for Add_With_Ovflo_Check
+
    begin
       if X >= 0 then
          if Y < 0 or else R >= 0 then
+            Prove_Non_Negative_X;
             return R;
          end if;
 
       else -- X < 0
          if Y > 0 or else R < 0 then
+            Prove_Negative_X;
             return R;
          end if;
       end if;
 
+      Prove_Overflow_Case;
       Raise_Error;
    end Add_With_Ovflo_Check;
 
@@ -147,24 +711,173 @@ package body System.Arith_Double is
 
       T1, T2     : Double_Uns;
       Du, Qu, Ru : Double_Uns;
-      Den_Pos    : Boolean;
+      Den_Pos    : constant Boolean := (Y < 0) = (Z < 0);
+
+      --  Local ghost variables
+
+      Mult  : constant Big_Integer := abs (Big (Y) * Big (Z)) with Ghost;
+      Quot  : Big_Integer with Ghost;
+      Big_R : Big_Integer with Ghost;
+      Big_Q : Big_Integer with Ghost;
+
+      --  Local lemmas
+
+      function Is_Division_By_Zero_Case return Boolean is
+        (Y = 0 or else Z = 0)
+      with Ghost;
+
+      function Is_Overflow_Case return Boolean is
+        (not In_Double_Int_Range (Big (X) / (Big (Y) * Big (Z))))
+      with
+        Ghost,
+        Pre => Y /= 0 and Z /= 0;
+
+      procedure Prove_Overflow_Case
+      with
+        Ghost,
+        Pre  => X = Double_Int'First and then Big (Y) * Big (Z) = -1,
+        Post => Is_Overflow_Case;
+      --  Proves the special case where -2**(Double_Size - 1) is divided by -1,
+      --  generating an overflow.
+
+      procedure Prove_Quotient_Zero
+      with
+        Ghost,
+        Pre  => Mult >= Big_2xx64
+          and then
+            not (Mult = Big_2xx64 and then X = Double_Int'First and then Round)
+          and then Q = 0
+          and then R = X,
+        Post => Big (R) = Big (X) rem (Big (Y) * Big (Z))
+          and then
+            (if Round then
+               Big (Q) = Round_Quotient (Big (X), Big (Y) * Big (Z),
+                                         Big (X) / (Big (Y) * Big (Z)),
+                                         Big (R))
+             else Big (Q) = Big (X) / (Big (Y) * Big (Z)));
+      --  Proves the general case where divisor doesn't fit in Double_Uns and
+      --  quotient is 0.
+
+      procedure Prove_Round_To_One
+      with
+        Ghost,
+        Pre  => Mult = Big_2xx64
+          and then X = Double_Int'First
+          and then Q = (if Den_Pos then -1 else 1)
+          and then R = X
+          and then Round,
+        Post => Big (R) = Big (X) rem (Big (Y) * Big (Z))
+          and then Big (Q) = Round_Quotient (Big (X), Big (Y) * Big (Z),
+                                             Big (X) / (Big (Y) * Big (Z)),
+                                             Big (R));
+      --  Proves the special case where the divisor doesn't fit in Double_Uns
+      --  but quotient is still 1 or -1 due to rounding
+      --  (abs (Y*Z) = 2**Double_Size and X = -2**(Double_Size - 1) and Round).
+
+      procedure Prove_Rounding_Case
+      with
+        Ghost,
+        Pre  => Mult /= 0
+          and then Quot = Big (X) / (Big (Y) * Big (Z))
+          and then Big_R = Big (X) rem (Big (Y) * Big (Z))
+          and then Big_Q =
+            Round_Quotient (Big (X), Big (Y) * Big (Z), Quot, Big_R)
+          and then Big (Ru) = abs Big_R
+          and then Big (Du) = Mult
+          and then Big (Qu) =
+            (if Ru > (Du - Double_Uns'(1)) / Double_Uns'(2)
+             then abs Quot + 1
+             else abs Quot),
+        Post => abs Big_Q = Big (Qu);
+      --  Proves correctness of the rounding of the unsigned quotient
+
+      procedure Prove_Signs
+      with
+        Ghost,
+        Pre  => Mult /= 0
+          and then Quot = Big (X) / (Big (Y) * Big (Z))
+          and then Big_R = Big (X) rem (Big (Y) * Big (Z))
+          and then Big_Q =
+            (if Round then
+               Round_Quotient (Big (X), Big (Y) * Big (Z), Quot, Big_R)
+             else Quot)
+          and then Big (Ru) = abs Big_R
+          and then Big (Qu) = abs Big_Q
+          and then R = (if X >= 0 then To_Int (Ru) else To_Int (-Ru))
+          and then
+            Q = (if (X >= 0) = Den_Pos then To_Int (Qu) else To_Int (-Qu))
+          and then not (X = Double_Int'First and then Big (Y) * Big (Z) = -1),
+        Post => Big (R) = Big_R and then Big (Q) = Big_Q;
+      --  Proves final signs match the intended result after the unsigned
+      --  division is done.
+
+      -----------------------------
+      -- Local lemma null bodies --
+      -----------------------------
+
+      procedure Prove_Overflow_Case is null;
+      procedure Prove_Quotient_Zero is null;
+      procedure Prove_Round_To_One is null;
+
+      -------------------------
+      -- Prove_Rounding_Case --
+      -------------------------
+
+      procedure Prove_Rounding_Case is
+      begin
+         if Same_Sign (Big (X), Big (Y) * Big (Z)) then
+            null;
+         end if;
+      end Prove_Rounding_Case;
+
+      -----------------
+      -- Prove_Signs --
+      -----------------
+
+      procedure Prove_Signs is
+      begin
+         if (X >= 0) = Den_Pos then
+            pragma Assert (Quot >= 0);
+            pragma Assert (Big_Q >= 0);
+            pragma Assert (Big (Q) = Big_Q);
+         else
+            pragma Assert ((X >= 0) /= (Big (Y) * Big (Z) >= 0));
+            pragma Assert (Quot <= 0);
+            pragma Assert (Big_Q <= 0);
+            pragma Assert (Big (R) = Big_R);
+         end if;
+      end Prove_Signs;
+
+   --  Start of processing for Double_Divide
 
    begin
       if Yu = 0 or else Zu = 0 then
+         pragma Assert (Is_Division_By_Zero_Case);
          Raise_Error;
+         pragma Annotate
+           (GNATprove, Intentional, "call to nonreturning subprogram",
+            "Constraint_Error is raised in case of division by zero");
       end if;
 
-      --  Set final signs (RM 4.5.5(27-30))
-
-      Den_Pos := (Y < 0) = (Z < 0);
+      pragma Assert (Mult /= 0);
+      pragma Assert (Den_Pos = (Big (Y) * Big (Z) > 0));
+      Quot := Big (X) / (Big (Y) * Big (Z));
+      Big_R := Big (X) rem (Big (Y) * Big (Z));
+      if Round then
+         Big_Q := Round_Quotient (Big (X), Big (Y) * Big (Z), Quot, Big_R);
+      else
+         Big_Q := Quot;
+      end if;
+      Lemma_Mult_Decomposition (Mult, Yu, Zu, Yhi, Ylo, Zhi, Zlo);
 
       --  Compute Y * Z. Note that if the result overflows Double_Uns, then
       --  the rounded result is zero, except for the very special case where
-      --  X = -2 ** (Double_Size - 1) and abs(Y*Z) = 2 ** Double_Size, when
+      --  X = -2 ** (Double_Size - 1) and abs (Y * Z) = 2 ** Double_Size, when
       --  Round is True.
 
       if Yhi /= 0 then
          if Zhi /= 0 then
+            R := X;
 
             --  Handle the special case when Round is True
 
@@ -176,11 +889,23 @@ package body System.Arith_Double is
               and then Round
             then
                Q := (if Den_Pos then -1 else 1);
+
+               Prove_Round_To_One;
+
             else
                Q := 0;
+
+               pragma Assert (Big (Double_Uns'(Yhi * Zhi)) >= 1);
+               if Yhi > 1 or else Zhi > 1 then
+                  pragma Assert (Big (Double_Uns'(Yhi * Zhi)) > 1);
+               elsif Zlo > 0 then
+                  pragma Assert (Big (Double_Uns'(Yhi * Zlo)) > 0);
+               elsif Ylo > 0 then
+                  pragma Assert (Big (Double_Uns'(Ylo * Zhi)) > 0);
+               end if;
+               Prove_Quotient_Zero;
             end if;
 
-            R := X;
             return;
          else
             T2 := Yhi * Zlo;
@@ -191,9 +916,34 @@ package body System.Arith_Double is
       end if;
 
       T1 := Ylo * Zlo;
+
+      pragma Assert (Big (T2) = Big (Double_Uns'(Yhi * Zlo))
+                              + Big (Double_Uns'(Ylo * Zhi)));
+      Lemma_Mult_Distribution (Big_2xx32, Big (Double_Uns'(Yhi * Zlo)),
+                                          Big (Double_Uns'(Ylo * Zhi)));
+      pragma Assert (Mult = Big_2xx32 * Big (T2) + Big (T1));
+      Lemma_Hi_Lo (T1, Hi (T1), Lo (T1));
+      pragma Assert
+        (Mult = Big_2xx32 * Big (T2)
+              + Big_2xx32 * Big (Double_Uns (Hi (T1)))
+                          + Big (Double_Uns (Lo (T1))));
+      Lemma_Mult_Distribution (Big_2xx32, Big (T2),
+                                          Big (Double_Uns (Hi (T1))));
+      Lemma_Add_Commutation (T2, Hi (T1));
+
       T2 := T2 + Hi (T1);
 
+      pragma Assert (Mult = Big_2xx32 * Big (T2) + Big (Double_Uns (Lo (T1))));
+      Lemma_Hi_Lo (T2, Hi (T2), Lo (T2));
+      Lemma_Mult_Distribution (Big_2xx32, Big (Double_Uns (Hi (T2))),
+                                          Big (Double_Uns (Lo (T2))));
+      pragma Assert
+        (Mult = Big_2xx64 * Big (Double_Uns (Hi (T2)))
+              + Big_2xx32 * Big (Double_Uns (Lo (T2)))
+                          + Big (Double_Uns (Lo (T1))));
+
       if Hi (T2) /= 0 then
+         R := X;
 
          --  Handle the special case when Round is True
 
@@ -204,20 +954,53 @@ package body System.Arith_Double is
            and then Round
          then
             Q := (if Den_Pos then -1 else 1);
+
+            Prove_Round_To_One;
+
          else
             Q := 0;
+
+            pragma Assert (Big (Double_Uns (Hi (T2))) >= 1);
+            pragma Assert (Big (Double_Uns (Lo (T2))) >= 0);
+            pragma Assert (Big (Double_Uns (Lo (T1))) >= 0);
+            pragma Assert (Mult >= Big_2xx64);
+            if Hi (T2) > 1 then
+               pragma Assert (Big (Double_Uns (Hi (T2))) > 1);
+            elsif Lo (T2) > 0 then
+               pragma Assert (Big (Double_Uns (Lo (T2))) > 0);
+            elsif Lo (T1) > 0 then
+               pragma Assert (Double_Uns (Lo (T1)) > 0);
+               Lemma_Gt_Commutation (Double_Uns (Lo (T1)), 0);
+               pragma Assert (Big (Double_Uns (Lo (T1))) > 0);
+            end if;
+            Prove_Quotient_Zero;
          end if;
 
-         R := X;
          return;
       end if;
 
       Du := Lo (T2) & Lo (T1);
 
+      Lemma_Hi_Lo (Du, Lo (T2), Lo (T1));
+      pragma Assert (Mult = Big (Du));
+      pragma Assert (Du /= 0);
+      --  Multiplication of 2-limb arguments Yu and Zu leads to 4-limb result
+      --  (where each limb is a single value). Cases where 4 limbs are needed
+      --  require Yhi /= 0 and Zhi /= 0 and lead to early exit. Remaining cases
+      --  where 3 limbs are needed correspond to Hi(T2) /= 0 and lead to early
+      --  exit. Thus, at this point, the result fits in 2 limbs which are
+      --  exactly Lo (T2) and Lo (T1), which corresponds to the value of Du.
+      --  As the case where one of Yu or Zu is null also led to early exit,
+      --  we have Du /= 0 here.
+
       --  Check overflow case of largest negative number divided by -1
 
       if X = Double_Int'First and then Du = 1 and then not Den_Pos then
+         Prove_Overflow_Case;
          Raise_Error;
+         pragma Annotate
+           (GNATprove, Intentional, "call to nonreturning subprogram",
+            "Constraint_Error is raised in case of overflow");
       end if;
 
       --  Perform the actual division
@@ -231,31 +1014,52 @@ package body System.Arith_Double is
       --  exactly Lo(T2) and Lo(T1), which corresponds to the value of Du.
       --  As the case where one of Yu or Zu is null also led to early exit,
       --  we have Du/=0 here.
+
       Qu := Xu / Du;
       Ru := Xu rem Du;
 
+      Lemma_Div_Commutation (Xu, Du);
+      Lemma_Abs_Div_Commutation (Big (X), Big (Y) * Big (Z));
+      Lemma_Abs_Commutation (X);
+      pragma Assert (abs Quot = Big (Qu));
+      Lemma_Rem_Commutation (Xu, Du);
+      Lemma_Abs_Rem_Commutation (Big (X), Big (Y) * Big (Z));
+      pragma Assert (abs Big_R = Big (Ru));
+
       --  Deal with rounding case
 
-      if Round and then Ru > (Du - Double_Uns'(1)) / Double_Uns'(2) then
-         Qu := Qu + Double_Uns'(1);
+      if Round then
+         if Ru > (Du - Double_Uns'(1)) / Double_Uns'(2) then
+            Lemma_Add_Commutation (Qu, 1);
+
+            Qu := Qu + Double_Uns'(1);
+         end if;
+
+         Prove_Rounding_Case;
       end if;
 
+      pragma Assert (abs Big_Q = Big (Qu));
+
+      --  Set final signs (RM 4.5.5(27-30))
+
       --  Case of dividend (X) sign positive
 
       if X >= 0 then
          R := To_Int (Ru);
-         Q := (if Den_Pos then To_Int (Qu) else -To_Int (Qu));
-
-      --  Case of dividend (X) sign negative
+         Q := (if Den_Pos then To_Int (Qu) else To_Int (-Qu));
 
       --  We perform the unary minus operation on the unsigned value
       --  before conversion to signed, to avoid a possible overflow
       --  for value -2 ** (Double_Size - 1), both for computing R and Q.
 
+      --  Case of dividend (X) sign negative
+
       else
          R := To_Int (-Ru);
          Q := (if Den_Pos then To_Int (-Qu) else To_Int (Qu));
       end if;
+
+      Prove_Signs;
    end Double_Divide;
 
    ---------
@@ -278,6 +1082,259 @@ package body System.Arith_Double is
    end Le3;
 
    -------------------------------
+   -- Lemma_Abs_Div_Commutation --
+   -------------------------------
+
+   procedure Lemma_Abs_Div_Commutation (X, Y : Big_Integer) is
+   begin
+      if Y < 0 then
+         if X < 0 then
+            pragma Assert (abs (X / Y) = abs (X / (-Y)));
+         else
+            Lemma_Neg_Div (X, Y);
+            pragma Assert (abs (X / Y) = abs ((-X) / (-Y)));
+         end if;
+      end if;
+   end Lemma_Abs_Div_Commutation;
+
+   -------------------------------
+   -- Lemma_Abs_Rem_Commutation --
+   -------------------------------
+
+   procedure Lemma_Abs_Rem_Commutation (X, Y : Big_Integer) is
+   begin
+      if Y < 0 then
+         Lemma_Neg_Rem (X, Y);
+         if X < 0 then
+            pragma Assert (X rem Y = -((-X) rem (-Y)));
+            pragma Assert (abs (X rem Y) = (abs X) rem (abs Y));
+         else
+            pragma Assert (abs (X rem Y) = (abs X) rem (abs Y));
+         end if;
+      end if;
+   end Lemma_Abs_Rem_Commutation;
+
+   ------------------------
+   -- Lemma_Double_Shift --
+   ------------------------
+
+   procedure Lemma_Double_Shift (X : Double_Uns; S, S1 : Natural) is
+   begin
+      Lemma_Double_Shift (X, Double_Uns (S), Double_Uns (S1));
+      pragma Assert (Shift_Left (Shift_Left (X, S), S1)
+        = Shift_Left (Shift_Left (X, S), Natural (Double_Uns (S1))));
+      pragma Assert (Shift_Left (X, S + S1)
+        = Shift_Left (X, Natural (Double_Uns (S + S1))));
+   end Lemma_Double_Shift;
+
+   ------------------------------
+   -- Lemma_Double_Shift_Right --
+   ------------------------------
+
+   procedure Lemma_Double_Shift_Right (X : Double_Uns; S, S1 : Natural) is
+   begin
+      Lemma_Double_Shift_Right (X, Double_Uns (S), Double_Uns (S1));
+      pragma Assert (Shift_Right (Shift_Right (X, S), S1)
+        = Shift_Right (Shift_Right (X, S), Natural (Double_Uns (S1))));
+      pragma Assert (Shift_Right (X, S + S1)
+        = Shift_Right (X, Natural (Double_Uns (S + S1))));
+   end Lemma_Double_Shift_Right;
+
+   -----------------
+   -- Lemma_Hi_Lo --
+   -----------------
+
+   procedure Lemma_Hi_Lo (Xu : Double_Uns; Xhi, Xlo : Single_Uns) is
+   begin
+      pragma Assert (Double_Uns (Xhi) = Xu / Double_Uns'(2 ** Single_Size));
+      pragma Assert (Double_Uns (Xlo) = Xu mod 2 ** Single_Size);
+   end Lemma_Hi_Lo;
+
+   -------------------
+   -- Lemma_Hi_Lo_3 --
+   -------------------
+
+   procedure Lemma_Hi_Lo_3 (Xu : Double_Uns; Xhi, Xlo : Single_Uns) is
+   begin
+      Lemma_Hi_Lo (Xu, Xhi, Xlo);
+   end Lemma_Hi_Lo_3;
+
+   ------------------------------
+   -- Lemma_Mult_Decomposition --
+   ------------------------------
+
+   procedure Lemma_Mult_Decomposition
+     (Mult               : Big_Integer;
+      Xu, Yu             : Double_Uns;
+      Xhi, Xlo, Yhi, Ylo : Single_Uns)
+   is
+   begin
+      Lemma_Hi_Lo (Xu, Xhi, Xlo);
+      Lemma_Hi_Lo (Yu, Yhi, Ylo);
+
+      pragma Assert
+        (Mult = (Big_2xx32 * Big (Double_Uns (Xhi)) + Big (Double_Uns (Xlo))) *
+                (Big_2xx32 * Big (Double_Uns (Yhi)) + Big (Double_Uns (Ylo))));
+      pragma Assert (Mult =
+        Big_2xx32 * Big_2xx32 * Big (Double_Uns (Xhi)) * Big (Double_Uns (Yhi))
+                  + Big_2xx32 * Big (Double_Uns (Xhi)) * Big (Double_Uns (Ylo))
+                  + Big_2xx32 * Big (Double_Uns (Xlo)) * Big (Double_Uns (Yhi))
+                            + Big (Double_Uns (Xlo)) * Big (Double_Uns (Ylo)));
+      Lemma_Deep_Mult_Commutation (Big_2xx32 * Big_2xx32, Xhi, Yhi);
+      Lemma_Deep_Mult_Commutation (Big_2xx32, Xhi, Ylo);
+      Lemma_Deep_Mult_Commutation (Big_2xx32, Xlo, Yhi);
+      Lemma_Mult_Commutation (Xlo, Ylo);
+      pragma Assert (Mult =
+                       Big_2xx32 * Big_2xx32 * Big (Double_Uns'(Xhi * Yhi))
+                                 + Big_2xx32 * Big (Double_Uns'(Xhi * Ylo))
+                                 + Big_2xx32 * Big (Double_Uns'(Xlo * Yhi))
+                                             + Big (Double_Uns'(Xlo * Ylo)));
+   end Lemma_Mult_Decomposition;
+
+   -------------------
+   -- Lemma_Neg_Div --
+   -------------------
+
+   procedure Lemma_Neg_Div (X, Y : Big_Integer) is
+   begin
+      pragma Assert ((-X) / (-Y) = -(X / (-Y)));
+      pragma Assert (X / (-Y) = -(X / Y));
+   end Lemma_Neg_Div;
+
+   -----------------------
+   -- Lemma_Powers_Of_2 --
+   -----------------------
+
+   procedure Lemma_Powers_Of_2 (M, N : Natural) is
+   begin
+      if M + N < Double_Size then
+         pragma Assert (Double_Uns'(2**M) * Double_Uns'(2**N)
+                        = Double_Uns'(2**(M + N)));
+      end if;
+
+      Lemma_Powers_Of_2_Commutation (M);
+      Lemma_Powers_Of_2_Commutation (N);
+      Lemma_Powers_Of_2_Commutation (M + N);
+
+      if M + N < Double_Size then
+         pragma Assert (Big (Double_Uns'(2))**M * Big (Double_Uns'(2))**N
+                        = Big (Double_Uns'(2))**(M + N));
+         Lemma_Powers_Of_2_Increasing (M + N, Double_Size);
+         Lemma_Mult_Commutation (2 ** M, 2 ** N, 2 ** (M + N));
+      else
+         pragma Assert (Big (Double_Uns'(2))**M * Big (Double_Uns'(2))**N
+                        = Big (Double_Uns'(2))**(M + N));
+      end if;
+   end Lemma_Powers_Of_2;
+
+   -----------------------------------
+   -- Lemma_Powers_Of_2_Commutation --
+   -----------------------------------
+
+   procedure Lemma_Powers_Of_2_Commutation (M : Natural) is
+   begin
+      if M > 0 then
+         Lemma_Powers_Of_2_Commutation (M - 1);
+         pragma Assert (Big (Double_Uns'(2))**(M - 1) = Big_2xx (M - 1));
+         pragma Assert (Big (Double_Uns'(2))**M = Big_2xx (M - 1) * 2);
+         if M < Double_Size then
+            Lemma_Powers_Of_2_Increasing (M - 1, Double_Size - 1);
+            Lemma_Bounded_Powers_Of_2_Increasing (M - 1, Double_Size - 1);
+            pragma Assert (Double_Uns'(2 ** (M - 1)) * 2 = Double_Uns'(2**M));
+            Lemma_Mult_Commutation
+              (Double_Uns'(2 ** (M - 1)), 2, Double_Uns'(2**M));
+            pragma Assert (Big (Double_Uns'(2))**M = Big_2xx (M));
+         end if;
+      end if;
+   end Lemma_Powers_Of_2_Commutation;
+
+   ----------------------------------
+   -- Lemma_Powers_Of_2_Increasing --
+   ----------------------------------
+
+   procedure Lemma_Powers_Of_2_Increasing (M, N : Natural) is
+   begin
+      if M + 1 < N then
+         Lemma_Powers_Of_2_Increasing (M + 1, N);
+      end if;
+   end Lemma_Powers_Of_2_Increasing;
+
+   -------------------
+   -- Lemma_Rem_Abs --
+   -------------------
+
+   procedure Lemma_Rem_Abs (X, Y : Big_Integer) is
+   begin
+      Lemma_Neg_Rem (X, Y);
+   end Lemma_Rem_Abs;
+
+   -----------------------
+   -- Lemma_Shift_Right --
+   -----------------------
+
+   procedure Lemma_Shift_Right (X : Double_Uns; Shift : Natural) is
+      XX : Double_Uns := X;
+   begin
+      for J in 1 .. Shift loop
+         XX := Shift_Right (XX, 1);
+         Lemma_Double_Shift_Right (X, J - 1, 1);
+         pragma Loop_Invariant (XX = Shift_Right (X, J));
+         pragma Loop_Invariant (XX = X / Double_Uns'(2) ** J);
+      end loop;
+   end Lemma_Shift_Right;
+
+   ------------------------------
+   -- Lemma_Shift_Without_Drop --
+   ------------------------------
+
+   procedure Lemma_Shift_Without_Drop
+     (X, Y  : Double_Uns;
+      Mask  : Single_Uns;
+      Shift : Natural)
+   is
+      pragma Unreferenced (Mask);
+
+      procedure Lemma_Bound
+      with
+        Pre  => Shift <= Single_Size
+          and then X <= 2**Single_Size
+            * Double_Uns'(2**(Single_Size - Shift) - 1)
+            + Single_Uns'(2**Single_Size - 1),
+        Post => X <= 2**(Double_Size - Shift) - 1;
+
+      procedure Lemma_Exp_Pos (N : Integer)
+      with
+        Pre  => N in 0 .. Double_Size - 1,
+        Post => Double_Uns'(2**N) > 0;
+
+      -----------------------------
+      -- Local lemma null bodies --
+      -----------------------------
+
+      procedure Lemma_Bound is null;
+      procedure Lemma_Exp_Pos (N : Integer) is null;
+
+   --  Start of processing for Lemma_Shift_Without_Drop
+
+   begin
+      if Shift = 0 then
+         pragma Assert (Big (Y) = Big_2xx (Shift) * Big (X));
+         return;
+      end if;
+
+      Lemma_Bound;
+      Lemma_Exp_Pos (Double_Size - Shift);
+      pragma Assert (X < 2**(Double_Size - Shift));
+      pragma Assert (Big (X) < Big_2xx (Double_Size - Shift));
+      pragma Assert (Y = 2**Shift * X);
+      pragma Assert (Big_2xx (Shift) * Big (X)
+                     < Big_2xx (Shift) * Big_2xx (Double_Size - Shift));
+      Lemma_Powers_Of_2 (Shift, Double_Size - Shift);
+      Lemma_Mult_Commutation (2**Shift, X, Y);
+      pragma Assert (Big (Y) = Big_2xx (Shift) * Big (X));
+   end Lemma_Shift_Without_Drop;
+
+   -------------------------------
    -- Multiply_With_Ovflo_Check --
    -------------------------------
 
@@ -292,9 +1349,148 @@ package body System.Arith_Double is
 
       T1, T2 : Double_Uns;
 
+      --  Local ghost variables
+
+      Mult : constant Big_Integer := abs (Big (X) * Big (Y)) with Ghost;
+
+      --  Local lemmas
+
+      procedure Prove_Both_Too_Large
+      with
+        Ghost,
+        Pre  => Xhi /= 0
+          and then Yhi /= 0
+          and then Mult =
+            Big_2xx32 * Big_2xx32 * (Big (Double_Uns'(Xhi * Yhi)))
+                      + Big_2xx32 * (Big (Double_Uns'(Xhi * Ylo)))
+                      + Big_2xx32 * (Big (Double_Uns'(Xlo * Yhi)))
+                                  + (Big (Double_Uns'(Xlo * Ylo))),
+        Post => not In_Double_Int_Range (Big (X) * Big (Y));
+
+      procedure Prove_Final_Decomposition
+      with
+        Ghost,
+        Pre  => In_Double_Int_Range (Big (X) * Big (Y))
+          and then Mult = Big_2xx32 * Big (T2) + Big (Double_Uns (Lo (T1)))
+          and then Hi (T2) = 0,
+        Post => Mult = Big (Lo (T2) & Lo (T1));
+
+      procedure Prove_Neg_Int
+      with
+        Ghost,
+        Pre  => In_Double_Int_Range (Big (X) * Big (Y))
+          and then Mult = Big (T2)
+          and then ((X >= 0 and then Y < 0) or else (X < 0 and then Y >= 0)),
+        Post => To_Neg_Int (T2) = X * Y;
+
+      procedure Prove_Pos_Int
+      with
+        Ghost,
+        Pre  => In_Double_Int_Range (Big (X) * Big (Y))
+          and then Mult = Big (T2)
+          and then ((X >= 0 and then Y >= 0) or else (X < 0 and then Y < 0)),
+        Post => To_Pos_Int (T2) = X * Y;
+
+      procedure Prove_Result_Too_Large
+      with
+        Ghost,
+        Pre  => Mult = Big_2xx32 * Big (T2) + Big (Double_Uns (Lo (T1)))
+          and then Hi (T2) /= 0,
+        Post => not In_Double_Int_Range (Big (X) * Big (Y));
+
+      procedure Prove_Too_Large
+      with
+        Ghost,
+        Pre  => abs (Big (X) * Big (Y)) >= Big_2xx64,
+        Post => not In_Double_Int_Range (Big (X) * Big (Y));
+
+      --------------------------
+      -- Prove_Both_Too_Large --
+      --------------------------
+
+      procedure Prove_Both_Too_Large is
+      begin
+         pragma Assert
+           (Mult >= Big_2xx32 * Big_2xx32 * Big (Double_Uns'(Xhi * Yhi)));
+         pragma Assert (Double_Uns (Xhi) * Double_Uns (Yhi) >= 1);
+         pragma Assert (Mult >= Big_2xx32 * Big_2xx32);
+         Prove_Too_Large;
+      end Prove_Both_Too_Large;
+
+      -------------------------------
+      -- Prove_Final_Decomposition --
+      -------------------------------
+
+      procedure Prove_Final_Decomposition is
+      begin
+         Lemma_Hi_Lo (T2, Hi (T2), Lo (T2));
+         pragma Assert (Mult = Big_2xx32 * Big (Double_Uns (Lo (T2)))
+                                         + Big (Double_Uns (Lo (T1))));
+         pragma Assert (Mult <= Big_2xx63);
+         Lemma_Mult_Commutation (X, Y);
+         pragma Assert (Mult = abs (Big (X * Y)));
+         Lemma_Word_Commutation (Lo (T2));
+         pragma Assert (Mult = Big (Double_Uns'(2 ** Single_Size)
+                          * Double_Uns (Lo (T2)))
+                          + Big (Double_Uns (Lo (T1))));
+         Lemma_Add_Commutation (Double_Uns'(2 ** Single_Size)
+                                  * Double_Uns (Lo (T2)),
+                                Lo (T1));
+         pragma Assert (Mult = Big (Double_Uns'(2 ** Single_Size)
+                          * Double_Uns (Lo (T2)) + Lo (T1)));
+         pragma Assert (Lo (T2) & Lo (T1) = Double_Uns'(2 ** Single_Size)
+                          * Double_Uns (Lo (T2)) + Lo (T1));
+      end Prove_Final_Decomposition;
+
+      -------------------
+      -- Prove_Neg_Int --
+      -------------------
+
+      procedure Prove_Neg_Int is
+      begin
+         pragma Assert (X * Y <= 0);
+         pragma Assert (Mult = -Big (X * Y));
+      end Prove_Neg_Int;
+
+      -------------------
+      -- Prove_Pos_Int --
+      -------------------
+
+      procedure Prove_Pos_Int is
+      begin
+         pragma Assert (X * Y >= 0);
+         pragma Assert (Mult = Big (X * Y));
+      end Prove_Pos_Int;
+
+      ----------------------------
+      -- Prove_Result_Too_Large --
+      ----------------------------
+
+      procedure Prove_Result_Too_Large is
+      begin
+         pragma Assert (Mult >= Big_2xx32 * Big (T2));
+         Lemma_Hi_Lo (T2, Hi (T2), Lo (T2));
+         pragma Assert
+           (Mult >= Big_2xx32 * Big_2xx32 * Big (Double_Uns (Hi (T2))));
+         pragma Assert (Double_Uns (Hi (T2)) >= 1);
+         pragma Assert (Mult >= Big_2xx32 * Big_2xx32);
+         Prove_Too_Large;
+      end Prove_Result_Too_Large;
+
+      ---------------------
+      -- Prove_Too_Large --
+      ---------------------
+
+      procedure Prove_Too_Large is null;
+
+   --  Start of processing for Multiply_With_Ovflo_Check
+
    begin
+      Lemma_Mult_Decomposition (Mult, Xu, Yu, Xhi, Xlo, Yhi, Ylo);
+
       if Xhi /= 0 then
          if Yhi /= 0 then
+            Prove_Both_Too_Large;
             Raise_Error;
          else
             T2 := Xhi * Ylo;
@@ -311,28 +1507,50 @@ package body System.Arith_Double is
       --  result from the upper halves of the input values.
 
       T1 := Xlo * Ylo;
+
+      pragma Assert (Big (T2) = Big (Double_Uns'(Xhi * Ylo))
+                              + Big (Double_Uns'(Xlo * Yhi)));
+      Lemma_Mult_Distribution (Big_2xx32, Big (Double_Uns'(Xhi * Ylo)),
+                                          Big (Double_Uns'(Xlo * Yhi)));
+      pragma Assert (Mult = Big_2xx32 * Big (T2) + Big (T1));
+      Lemma_Add_Commutation (T2, Hi (T1));
+      pragma Assert
+        (Big (T2 + Hi (T1)) = Big (T2) + Big (Double_Uns (Hi (T1))));
+
       T2 := T2 + Hi (T1);
 
+      Lemma_Hi_Lo (T1, Hi (T1), Lo (T1));
+      pragma Assert (Mult = Big_2xx32 * Big (T2) + Big (Double_Uns (Lo (T1))));
+
       if Hi (T2) /= 0 then
+         Prove_Result_Too_Large;
          Raise_Error;
       end if;
 
+      Prove_Final_Decomposition;
+
       T2 := Lo (T2) & Lo (T1);
 
+      pragma Assert (Mult = Big (T2));
+
       if X >= 0 then
          if Y >= 0 then
+            Prove_Pos_Int;
             return To_Pos_Int (T2);
             pragma Annotate (CodePeer, Intentional, "precondition",
                              "Intentional Unsigned->Signed conversion");
          else
+            Prove_Neg_Int;
             return To_Neg_Int (T2);
          end if;
       else -- X < 0
          if Y < 0 then
+            Prove_Pos_Int;
             return To_Pos_Int (T2);
             pragma Annotate (CodePeer, Intentional, "precondition",
                              "Intentional Unsigned->Signed conversion");
          else
+            Prove_Neg_Int;
             return To_Neg_Int (T2);
          end if;
       end if;
@@ -346,6 +1564,9 @@ package body System.Arith_Double is
    procedure Raise_Error is
    begin
       raise Constraint_Error with "Double arithmetic overflow";
+      pragma Annotate
+        (GNATprove, Intentional, "exception might be raised",
+         "Procedure Raise_Error is called to signal input errors");
    end Raise_Error;
 
    -------------------
@@ -369,10 +1590,10 @@ package body System.Arith_Double is
       Zhi : Single_Uns := Hi (Zu);
       Zlo : Single_Uns := Lo (Zu);
 
-      D : array (1 .. 4) of Single_Uns;
+      D : array (1 .. 4) of Single_Uns with Relaxed_Initialization;
       --  The dividend, four digits (D(1) is high order)
 
-      Qd : array (1 .. 2) of Single_Uns;
+      Qd : array (1 .. 2) of Single_Uns with Relaxed_Initialization;
       --  The quotient digits, two digits (Qd(1) is high order)
 
       S1, S2, S3 : Single_Uns;
@@ -397,56 +1618,527 @@ package body System.Arith_Double is
       T1, T2, T3 : Double_Uns;
       --  Temporary values
 
+      --  Local ghost variables
+
+      Mult  : constant Big_Integer := abs (Big (X) * Big (Y)) with Ghost;
+      Quot  : Big_Integer with Ghost;
+      Big_R : Big_Integer with Ghost;
+      Big_Q : Big_Integer with Ghost;
+
+      --  Local lemmas
+
+      function Is_Division_By_Zero_Case return Boolean is (Z = 0) with Ghost;
+
+      procedure Prove_Dividend_Scaling
+      with
+        Ghost,
+        Pre  => D'Initialized
+          and then Scale <= Single_Size
+          and then Mult =
+            Big_2xx32 * Big_2xx32 * Big_2xx32 * Big (Double_Uns (D (1)))
+                      + Big_2xx32 * Big_2xx32 * Big (Double_Uns (D (2)))
+                                  + Big_2xx32 * Big (Double_Uns (D (3)))
+                                              + Big (Double_Uns (D (4)))
+          and then Big (D (1) & D (2)) * Big_2xx (Scale) < Big_2xx64
+          and then T1 = Shift_Left (D (1) & D (2), Scale)
+          and then T2 = Shift_Left (Double_Uns (D (3)), Scale)
+          and then T3 = Shift_Left (Double_Uns (D (4)), Scale),
+        Post => Mult * Big_2xx (Scale) =
+          Big_2xx32 * Big_2xx32 * Big_2xx32 * Big (Double_Uns (Hi (T1)))
+                    + Big_2xx32 * Big_2xx32 * Big (Double_Uns (Lo (T1) or
+                                                                      Hi (T2)))
+                                + Big_2xx32 * Big (Double_Uns (Lo (T2) or
+                                                                      Hi (T3)))
+                                            + Big (Double_Uns (Lo (T3)));
+      --  Proves the scaling of the 4-digit dividend actually multiplies it by
+      --  2**Scale.
+
+      procedure Prove_Multiplication (Q : Single_Uns)
+      with
+        Ghost,
+        Pre  => T1 = Q * Lo (Zu)
+          and then T2 = Q * Hi (Zu)
+          and then S3 = Lo (T1)
+          and then T3 = Hi (T1) + Lo (T2)
+          and then S2 = Lo (T3)
+          and then S1 = Hi (T3) + Hi (T2),
+        Post => Big3 (S1, S2, S3) = Big (Double_Uns (Q)) * Big (Zu);
+      --  Proves correctness of the multiplication of divisor by quotient to
+      --  compute amount to subtract.
+
+      procedure Prove_Overflow
+      with
+        Ghost,
+        Pre  => Z /= 0 and then Mult >= Big_2xx64 * Big (Double_Uns'(abs Z)),
+        Post => not In_Double_Int_Range (Big (X) * Big (Y) / Big (Z));
+      --  Proves overflow case when the quotient has at least 3 digits
+
+      procedure Prove_Qd_Calculation_Part_1 (J : Integer)
+      with
+        Ghost,
+        Pre  => J in 1 .. 2
+          and then D'Initialized
+          and then D (J) < Zhi
+          and then Hi (Zu) = Zhi
+          and then Qd (J)'Initialized
+          and then Qd (J) = Lo ((D (J) & D (J + 1)) / Zhi),
+        Post => Big (Double_Uns (Qd (J))) >=
+          Big3 (D (J), D (J + 1), D (J + 2)) / Big (Zu);
+      --  When dividing 3 digits by 2 digits, proves the initial calculation
+      --  of the quotient given by dividing the first 2 digits of the dividend
+      --  by the first digit of the divisor is not an underestimate (so
+      --  readjusting down works).
+
+      procedure Prove_Rescaling
+      with
+        Ghost,
+        Pre  => Scale <= Single_Size
+          and then Z /= 0
+          and then Mult * Big_2xx (Scale) = Big (Zu) * Big (Qu) + Big (Ru)
+          and then Big (Ru) < Big (Zu)
+          and then Big (Zu) = Big (Double_Uns'(abs Z)) * Big_2xx (Scale)
+          and then Quot = Big (X) * Big (Y) / Big (Z)
+          and then Big_R = Big (X) * Big (Y) rem Big (Z),
+        Post => abs Quot = Big (Qu)
+          and then abs Big_R = Big (Shift_Right (Ru, Scale));
+      --  Proves scaling back only the remainder is the right thing to do after
+      --  computing the scaled division.
+
+      procedure Prove_Rounding_Case
+      with
+        Ghost,
+        Pre  => Z /= 0
+          and then Quot = Big (X) * Big (Y) / Big (Z)
+          and then Big_R = Big (X) * Big (Y) rem Big (Z)
+          and then Big_Q =
+            Round_Quotient (Big (X) * Big (Y), Big (Z), Quot, Big_R)
+          and then Big (Ru) = abs Big_R
+          and then Big (Zu) = Big (Double_Uns'(abs Z)),
+        Post => abs Big_Q =
+          (if Ru > (Zu - Double_Uns'(1)) / Double_Uns'(2)
+           then abs Quot + 1
+           else abs Quot);
+      --  Proves correctness of the rounding of the unsigned quotient
+
+      procedure Prove_Sign_R
+      with
+        Ghost,
+        Pre  => Z /= 0 and then Big_R = Big (X) * Big (Y) rem Big (Z),
+        Post => In_Double_Int_Range (Big_R);
+
+      procedure Prove_Signs
+      with
+        Ghost,
+        Pre  => Z /= 0
+          and then Quot = Big (X) * Big (Y) / Big (Z)
+          and then Big_R = Big (X) * Big (Y) rem Big (Z)
+          and then Big_Q =
+            (if Round then
+               Round_Quotient (Big (X) * Big (Y), Big (Z), Quot, Big_R)
+             else Quot)
+          and then Big (Ru) = abs Big_R
+          and then Big (Qu) = abs Big_Q
+          and then In_Double_Int_Range (Big_Q)
+          and then In_Double_Int_Range (Big_R)
+          and then R =
+            (if (X >= 0) = (Y >= 0) then To_Pos_Int (Ru) else To_Neg_Int (Ru))
+          and then Q =
+            (if ((X >= 0) = (Y >= 0)) = (Z >= 0) then To_Pos_Int (Qu)
+             else To_Neg_Int (Qu)),  --  need to ensure To_Pos_Int precondition
+        Post => Big (R) = Big_R and then Big (Q) = Big_Q;
+      --  Proves final signs match the intended result after the unsigned
+      --  division is done.
+
+      procedure Prove_Z_Low
+      with
+        Ghost,
+        Pre  => Z /= 0
+          and then D'Initialized
+          and then Hi (abs Z) = 0
+          and then Lo (abs Z) = Zlo
+          and then Mult = Big_2xx32 * Big_2xx32 * Big (Double_Uns (D (2)))
+                                    + Big_2xx32 * Big (Double_Uns (D (3)))
+                                                + Big (Double_Uns (D (4)))
+          and then D (2) < Zlo
+          and then Quot = (Big (X) * Big (Y)) / Big (Z)
+          and then Big_R = (Big (X) * Big (Y)) rem Big (Z)
+          and then T1 = D (2) & D (3)
+          and then T2 = Lo (T1 rem Zlo) & D (4)
+          and then Qu = Lo (T1 / Zlo) & Lo (T2 / Zlo)
+          and then Ru = T2 rem Zlo,
+        Post => Big (Qu) = abs Quot
+          and then Big (Ru) = abs Big_R;
+      --  Proves the case where the divisor is only one digit
+
+      ----------------------------
+      -- Prove_Dividend_Scaling --
+      ----------------------------
+
+      procedure Prove_Dividend_Scaling is
+      begin
+         Lemma_Hi_Lo (D (1) & D (2), D (1), D (2));
+         pragma Assert (Mult * Big_2xx (Scale) =
+           Big_2xx32 * Big_2xx32 * Big_2xx (Scale) * Big (D (1) & D (2))
+                     + Big_2xx32 * Big_2xx (Scale) * Big (Double_Uns (D (3)))
+                                 + Big_2xx (Scale) * Big (Double_Uns (D (4))));
+         pragma Assert (Big_2xx (Scale) > 0);
+         Lemma_Lt_Mult (Big (Double_Uns (D (3))), Big_2xx32,
+                        Big_2xx (Scale), Big_2xx64);
+         Lemma_Lt_Mult (Big (Double_Uns (D (4))), Big_2xx32,
+                        Big_2xx (Scale), Big_2xx64);
+         Lemma_Mult_Commutation (2 ** Scale, D (1) & D (2), T1);
+         declare
+            Big_D12 : constant Big_Integer :=
+              Big_2xx (Scale) * Big (D (1) & D (2));
+            Big_T1  : constant Big_Integer := Big (T1);
+         begin
+            pragma Assert (Big_D12 = Big_T1);
+            pragma Assert (Big_2xx32 * Big_2xx32 * Big_D12
+                           = Big_2xx32 * Big_2xx32 * Big_T1);
+         end;
+         Lemma_Mult_Commutation (2 ** Scale, Double_Uns (D (3)), T2);
+         declare
+            Big_D3 : constant Big_Integer :=
+              Big_2xx (Scale) * Big (Double_Uns (D (3)));
+            Big_T2 : constant Big_Integer := Big (T2);
+         begin
+            pragma Assert (Big_D3 = Big_T2);
+            pragma Assert (Big_2xx32 * Big_D3 = Big_2xx32 * Big_T2);
+         end;
+         Lemma_Mult_Commutation (2 ** Scale, Double_Uns (D (4)), T3);
+         declare
+            Big_D4 : constant Big_Integer :=
+              Big_2xx (Scale) * Big (Double_Uns (D (4)));
+            Big_T3 : constant Big_Integer := Big (T3);
+         begin
+            pragma Assert (Big_D4 = Big_T3);
+         end;
+         pragma Assert (Mult * Big_2xx (Scale) =
+           Big_2xx32 * Big_2xx32 * Big (T1) + Big_2xx32 * Big (T2) + Big (T3));
+         Lemma_Hi_Lo (T1, Hi (T1), Lo (T1));
+         Lemma_Hi_Lo (T2, Hi (T2), Lo (T2));
+         Lemma_Hi_Lo (T3, Hi (T3), Lo (T3));
+         pragma Assert (Double_Uns (Lo (T1) or Hi (T2)) =
+                          Double_Uns (Lo (T1)) + Double_Uns (Hi (T2)));
+         pragma Assert (Double_Uns (Lo (T2) or Hi (T3)) =
+                          Double_Uns (Lo (T2)) + Double_Uns (Hi (T3)));
+      end Prove_Dividend_Scaling;
+
+      --------------------------
+      -- Prove_Multiplication --
+      --------------------------
+
+      procedure Prove_Multiplication (Q : Single_Uns) is
+      begin
+         Lemma_Hi_Lo (Zu, Hi (Zu), Lo (Zu));
+         Lemma_Hi_Lo (T1, Hi (T1), S3);
+         Lemma_Hi_Lo (T2, Hi (T2), Lo (T2));
+         Lemma_Hi_Lo (T3, Hi (T3), S2);
+         pragma Assert (Big (Double_Uns (Q)) * Big (Zu) =
+                          Big_2xx32 * Big_2xx32 * Big (Double_Uns (Hi (T2)))
+                        + Big_2xx32 * Big_2xx32 * Big (Double_Uns (Hi (T3)))
+                                    + Big_2xx32 * Big (Double_Uns (S2))
+                                                + Big (Double_Uns (S3)));
+         pragma Assert (Double_Uns (Hi (T3)) + Hi (T2) = Double_Uns (S1));
+         Lemma_Add_Commutation (Double_Uns (Hi (T3)), Hi (T2));
+         pragma Assert
+           (Big (Double_Uns (Hi (T3))) + Big (Double_Uns (Hi (T2))) =
+              Big (Double_Uns (S1)));
+      end Prove_Multiplication;
+
+      --------------------
+      -- Prove_Overflow --
+      --------------------
+
+      procedure Prove_Overflow is
+      begin
+         Lemma_Div_Ge (Mult, Big_2xx64, Big (Double_Uns'(abs Z)));
+         Lemma_Abs_Commutation (Z);
+         Lemma_Abs_Div_Commutation (Big (X) * Big (Y), Big (Z));
+      end Prove_Overflow;
+
+      ---------------------------------
+      -- Prove_Qd_Calculation_Part_1 --
+      ---------------------------------
+
+      procedure Prove_Qd_Calculation_Part_1 (J : Integer) is
+      begin
+         Lemma_Hi_Lo (D (J) & D (J + 1), D (J), D (J + 1));
+         Lemma_Lt_Commutation (Double_Uns (D (J)), Double_Uns (Zhi));
+         Lemma_Gt_Mult (Big (Double_Uns (Zhi)),
+                        Big (Double_Uns (D (J))) + 1,
+                        Big_2xx32, Big (D (J) & D (J + 1)));
+         Lemma_Div_Lt
+           (Big (D (J) & D (J + 1)), Big_2xx32, Big (Double_Uns (Zhi)));
+         Lemma_Div_Commutation (D (J) & D (J + 1), Double_Uns (Zhi));
+         Lemma_Lo_Is_Ident ((D (J) & D (J + 1)) / Zhi);
+         Lemma_Div_Definition (D (J) & D (J + 1), Zhi, Double_Uns (Qd (J)),
+                               (D (J) & D (J + 1)) rem Zhi);
+         Lemma_Lt_Commutation
+           ((D (J) & D (J + 1)) rem Zhi, Double_Uns (Zhi));
+         Lemma_Gt_Mult
+           ((Big (Double_Uns (Qd (J))) + 1) * Big (Double_Uns (Zhi)),
+            Big (D (J) & D (J + 1)) + 1, Big_2xx32,
+            Big3 (D (J), D (J + 1), D (J + 2)));
+         Lemma_Hi_Lo (Zu, Zhi, Lo (Zu));
+         Lemma_Gt_Mult (Big (Zu), Big_2xx32 * Big (Double_Uns (Zhi)),
+                        Big (Double_Uns (Qd (J))) + 1,
+                        Big3 (D (J), D (J + 1), D (J + 2)));
+         Lemma_Div_Lt (Big3 (D (J), D (J + 1), D (J + 2)),
+                       Big (Double_Uns (Qd (J))) + 1, Big (Zu));
+      end Prove_Qd_Calculation_Part_1;
+
+      ---------------------
+      -- Prove_Rescaling --
+      ---------------------
+
+      procedure Prove_Rescaling is
+      begin
+         Lemma_Div_Lt (Big (Ru), Big (Double_Uns'(abs Z)), Big_2xx (Scale));
+         Lemma_Div_Eq (Mult, Big (Double_Uns'(abs Z)) * Big (Qu),
+                       Big_2xx (Scale), Big (Ru));
+         Lemma_Rev_Div_Definition (Mult, Big (Double_Uns'(abs Z)),
+                                   Big (Qu), Big (Ru) / Big_2xx (Scale));
+         Lemma_Abs_Div_Commutation (Big (X) * Big (Y), Big (Z));
+         Lemma_Abs_Rem_Commutation (Big (X) * Big (Y), Big (Z));
+         Lemma_Abs_Commutation (Z);
+         Lemma_Shift_Right (Ru, Scale);
+      end Prove_Rescaling;
+
+      -------------------------
+      -- Prove_Rounding_Case --
+      -------------------------
+
+      procedure Prove_Rounding_Case is
+      begin
+         if Same_Sign (Big (X) * Big (Y), Big (Z)) then
+            null;
+         end if;
+      end Prove_Rounding_Case;
+
+      ------------------
+      -- Prove_Sign_R --
+      ------------------
+
+      procedure Prove_Sign_R is
+      begin
+         pragma Assert (In_Double_Int_Range (Big (Z)));
+      end Prove_Sign_R;
+
+      -----------------
+      -- Prove_Signs --
+      -----------------
+
+      procedure Prove_Signs is null;
+
+      -----------------
+      -- Prove_Z_Low --
+      -----------------
+
+      procedure Prove_Z_Low is
+      begin
+         Lemma_Hi_Lo (T1, D (2), D (3));
+         Lemma_Add_Commutation (Double_Uns (D (2)), 1);
+         pragma Assert
+           (Big (Double_Uns (D (2))) + 1 <= Big (Double_Uns (Zlo)));
+         Lemma_Div_Definition (T1, Zlo, T1 / Zlo, T1 rem Zlo);
+         pragma Assert (Double_Uns (Lo (T1 rem Zlo)) = T1 rem Zlo);
+         Lemma_Hi_Lo (T2, Lo (T1 rem Zlo), D (4));
+         pragma Assert (T1 rem Zlo + Double_Uns'(1) <= Double_Uns (Zlo));
+         Lemma_Add_Commutation (T1 rem Zlo, 1);
+         pragma Assert (Big (T1 rem Zlo) + 1 <= Big (Double_Uns (Zlo)));
+         Lemma_Div_Definition (T2, Zlo, T2 / Zlo, Ru);
+         pragma Assert
+           (Mult = Big (Double_Uns (Zlo)) *
+                     (Big_2xx32 * Big (T1 / Zlo) + Big (T2 / Zlo)) + Big (Ru));
+         Lemma_Div_Lt (Big (T1), Big_2xx32, Big (Double_Uns (Zlo)));
+         Lemma_Div_Commutation (T1, Double_Uns (Zlo));
+         Lemma_Lo_Is_Ident (T1 / Zlo);
+         Lemma_Div_Lt (Big (T2), Big_2xx32, Big (Double_Uns (Zlo)));
+         Lemma_Div_Commutation (T2, Double_Uns (Zlo));
+         Lemma_Lo_Is_Ident (T2 / Zlo);
+         Lemma_Hi_Lo (Qu, Lo (T1 / Zlo), Lo (T2 / Zlo));
+         Lemma_Substitution (Mult, Big (Double_Uns (Zlo)),
+                             Big_2xx32 * Big (T1 / Zlo) + Big (T2 / Zlo),
+                             Big (Qu), Big (Ru));
+         Lemma_Lt_Commutation (Ru, Double_Uns (Zlo));
+         Lemma_Rev_Div_Definition
+           (Mult, Big (Double_Uns (Zlo)), Big (Qu), Big (Ru));
+         pragma Assert (Double_Uns (Zlo) = abs Z);
+         Lemma_Abs_Commutation (Z);
+         Lemma_Abs_Div_Commutation (Big (X) * Big (Y), Big (Z));
+         Lemma_Abs_Rem_Commutation (Big (X) * Big (Y), Big (Z));
+      end Prove_Z_Low;
+
+   --  Start of processing for Scaled_Divide
+
    begin
+      if Z = 0 then
+         pragma Assert (Is_Division_By_Zero_Case);
+         Raise_Error;
+         pragma Annotate
+           (GNATprove, Intentional, "call to nonreturning subprogram",
+            "Constraint_Error is raised in case of division by zero");
+      end if;
+
+      Quot := Big (X) * Big (Y) / Big (Z);
+      Big_R := Big (X) * Big (Y) rem Big (Z);
+      if Round then
+         Big_Q := Round_Quotient (Big (X) * Big (Y), Big (Z), Quot, Big_R);
+      else
+         Big_Q := Quot;
+      end if;
+
       --  First do the multiplication, giving the four digit dividend
 
+      Lemma_Abs_Mult_Commutation (Big (X), Big (Y));
+      Lemma_Abs_Commutation (X);
+      Lemma_Abs_Commutation (Y);
+      Lemma_Mult_Decomposition (Mult, Xu, Yu, Xhi, Xlo, Yhi, Ylo);
+
       T1 := Xlo * Ylo;
       D (4) := Lo (T1);
       D (3) := Hi (T1);
 
+      Lemma_Hi_Lo (T1, D (3), D (4));
+
       if Yhi /= 0 then
          T1 := Xlo * Yhi;
+
+         Lemma_Hi_Lo (T1, Hi (T1), Lo (T1));
+
          T2 := D (3) + Lo (T1);
+
+         Lemma_Mult_Distribution (Big_2xx32,
+                                  Big (Double_Uns (D (3))),
+                                  Big (Double_Uns (Lo (T1))));
+         Lemma_Hi_Lo (T2, Hi (T2), Lo (T2));
+
          D (3) := Lo (T2);
          D (2) := Hi (T1) + Hi (T2);
 
+         pragma Assert (Double_Uns (Hi (T1)) + Hi (T2) = Double_Uns (D (2)));
+         Lemma_Add_Commutation (Double_Uns (Hi (T1)), Hi (T2));
+         pragma Assert
+           (Big (Double_Uns (Hi (T1))) + Big (Double_Uns (Hi (T2))) =
+              Big (Double_Uns (D (2))));
+
          if Xhi /= 0 then
             T1 := Xhi * Ylo;
+
+            Lemma_Hi_Lo (T1, Hi (T1), Lo (T1));
+
             T2 := D (3) + Lo (T1);
+
+            Lemma_Hi_Lo (T2, Hi (T2), Lo (T2));
+
             D (3) := Lo (T2);
             T3 := D (2) + Hi (T1);
+
+            Lemma_Add_Commutation (T3, Hi (T2));
+
             T3 := T3 + Hi (T2);
-            D (2) := Lo (T3);
-            D (1) := Hi (T3);
+            T2 := Double_Uns'(Xhi * Yhi);
 
-            T1 := (D (1) & D (2)) + Double_Uns'(Xhi * Yhi);
-            D (1) := Hi (T1);
+            Lemma_Hi_Lo (T2, Hi (T2), Lo (T2));
+            Lemma_Add_Commutation (T3, Lo (T2));
+
+            T1 := T3 + Lo (T2);
             D (2) := Lo (T1);
 
+            Lemma_Hi_Lo (T1, Hi (T1), D (2));
+
+            D (1) := Hi (T2) + Hi (T1);
+
+            pragma Assert
+              (Double_Uns (Hi (T2)) + Hi (T1) = Double_Uns (D (1)));
+            Lemma_Add_Commutation (Double_Uns (Hi (T2)), Hi (T1));
+            pragma Assert
+              (Big (Double_Uns (Hi (T2))) + Big (Double_Uns (Hi (T1))) =
+                   Big (Double_Uns (D (1))));
+
+            pragma Assert (Mult =
+              Big_2xx32 * Big_2xx32 * Big_2xx32 * Big (Double_Uns (D (1)))
+                        + Big_2xx32 * Big_2xx32 * Big (Double_Uns (D (2)))
+                                    + Big_2xx32 * Big (Double_Uns (D (3)))
+                                                + Big (Double_Uns (D (4))));
+
          else
             D (1) := 0;
          end if;
 
+         pragma Assert (Mult =
+           Big_2xx32 * Big_2xx32 * Big_2xx32 * Big (Double_Uns (D (1)))
+                     + Big_2xx32 * Big_2xx32 * Big (Double_Uns (D (2)))
+                                 + Big_2xx32 * Big (Double_Uns (D (3)))
+                                             + Big (Double_Uns (D (4))));
+
       else
          if Xhi /= 0 then
             T1 := Xhi * Ylo;
+
+            Lemma_Hi_Lo (T1, Hi (T1), Lo (T1));
+
             T2 := D (3) + Lo (T1);
+
+            Lemma_Mult_Distribution (Big_2xx32,
+                                     Big (Double_Uns (D (3))),
+                                     Big (Double_Uns (Lo (T1))));
+            Lemma_Hi_Lo (T2, Hi (T2), Lo (T2));
+
             D (3) := Lo (T2);
             D (2) := Hi (T1) + Hi (T2);
 
+            pragma Assert
+              (Double_Uns (Hi (T1)) + Hi (T2) = Double_Uns (D (2)));
+            Lemma_Add_Commutation (Double_Uns (Hi (T1)), Hi (T2));
+            pragma Assert
+              (Big (Double_Uns (Hi (T1))) + Big (Double_Uns (Hi (T2))) =
+                 Big (Double_Uns (D (2))));
+            pragma Assert
+              (Mult = Big_2xx32 * Big_2xx32 * Big (Double_Uns (D (2)))
+                                + Big_2xx32 * Big (Double_Uns (D (3)))
+                                            + Big (Double_Uns (D (4))));
          else
             D (2) := 0;
+
+            pragma Assert
+              (Mult = Big_2xx32 * Big_2xx32 * Big (Double_Uns (D (2)))
+                                + Big_2xx32 * Big (Double_Uns (D (3)))
+                                            + Big (Double_Uns (D (4))));
          end if;
 
          D (1) := 0;
       end if;
 
+      pragma Assert
+        (Mult = Big_2xx32 * Big_2xx32 * Big_2xx32 * Big (Double_Uns (D (1)))
+                          + Big_2xx32 * Big_2xx32 * Big (Double_Uns (D (2)))
+                                      + Big_2xx32 * Big (Double_Uns (D (3)))
+                                                  + Big (Double_Uns (D (4))));
+
       --  Now it is time for the dreaded multiple precision division. First an
       --  easy case, check for the simple case of a one digit divisor.
 
       if Zhi = 0 then
          if D (1) /= 0 or else D (2) >= Zlo then
+            if D (1) > 0 then
+               pragma Assert
+                 (Mult >= Big_2xx32 * Big_2xx32 * Big_2xx32
+                          * Big (Double_Uns (D (1))));
+               pragma Assert (Mult >= Big_2xx64 * Big_2xx32);
+               Lemma_Ge_Commutation (2 ** Single_Size, Zu);
+               pragma Assert (Mult >= Big_2xx64 * Big (Zu));
+            else
+               Lemma_Ge_Commutation (Double_Uns (D (2)), Zu);
+               pragma Assert (Mult >= Big_2xx64 * Big (Zu));
+            end if;
+
+            Prove_Overflow;
             Raise_Error;
+            pragma Annotate
+              (GNATprove, Intentional, "call to nonreturning subprogram",
+               "Constraint_Error is raised in case of overflow");
 
          --  Here we are dividing at most three digits by one digit
 
@@ -456,12 +2148,20 @@ package body System.Arith_Double is
 
             Qu := Lo (T1 / Zlo) & Lo (T2 / Zlo);
             Ru := T2 rem Zlo;
+
+            Prove_Z_Low;
          end if;
 
       --  If divisor is double digit and dividend is too large, raise error
 
       elsif (D (1) & D (2)) >= Zu then
+         Lemma_Hi_Lo (D (1) & D (2), D (1), D (2));
+         Lemma_Ge_Commutation (D (1) & D (2), Zu);
+         Prove_Overflow;
          Raise_Error;
+         pragma Annotate
+           (GNATprove, Intentional, "call to nonreturning subprogram",
+            "Constraint_Error is raised in case of overflow");
 
       --  This is the complex case where we definitely have a double digit
       --  divisor and a dividend of at least three digits. We use the classical
@@ -473,88 +2173,293 @@ package body System.Arith_Double is
          --  First normalize the divisor so that it has the leading bit on.
          --  We do this by finding the appropriate left shift amount.
 
-         Shift := Single_Size / 2;
-         Mask  := Shift_Left (2 ** (Single_Size / 2) - 1, Shift);
+         Shift := Single_Size;
+         Mask  := Single_Uns'Last;
          Scale := 0;
 
-         while Shift /= 0 loop
-            if (Hi (Zu) and Mask) = 0 then
-               Scale := Scale + Shift;
-               Zu := Shift_Left (Zu, Shift);
-            end if;
-
-            Shift := Shift / 2;
-            Mask := Shift_Left (Mask, Shift);
+         pragma Assert (Big_2xx (Scale) = 1);
+
+         while Shift > 1 loop
+            pragma Loop_Invariant (Scale <= Single_Size - Shift);
+            pragma Loop_Invariant ((Hi (Zu) and Mask) /= 0);
+            pragma Loop_Invariant
+              (Mask = Shift_Left (Single_Uns'Last, Single_Size - Shift));
+            pragma Loop_Invariant (Zu = Shift_Left (abs Z, Scale));
+            pragma Loop_Invariant (Big (Zu) =
+              Big (Double_Uns'(abs Z)) * Big_2xx (Scale));
+            pragma Loop_Invariant (Shift mod 2 = 0);
+            pragma Annotate
+              (GNATprove, False_Positive, "loop invariant",
+               "Shift actually is a power of 2");
+            --  Note : this scaling algorithm only works when Single_Size is a
+            --  power of 2.
+
+            declare
+               Shift_Prev : constant Natural := Shift with Ghost;
+               Mask_Prev  : constant Single_Uns := Mask with Ghost;
+               Zu_Prev    : constant Double_Uns := Zu with Ghost;
+
+               procedure Prove_Shifting
+               with
+                 Ghost,
+                 Pre  => Shift <= Single_Size / 2
+                   and then Zu = Shift_Left (Zu_Prev, Shift)
+                   and then Mask_Prev =
+                     Shift_Left (Single_Uns'Last, Single_Size - 2 * Shift)
+                   and then Mask =
+                     Shift_Left (Single_Uns'Last, Single_Size - Shift)
+                   and then (Hi (Zu_Prev) and Mask_Prev and not Mask) /= 0,
+                 Post => (Hi (Zu) and Mask) /= 0;
+
+               --------------------
+               -- Prove_Shifting --
+               --------------------
+
+               procedure Prove_Shifting is null;
+
+            begin
+               Shift := Shift / 2;
+
+               pragma Assert (Shift_Prev = 2 * Shift);
+
+               Mask := Shift_Left (Mask, Shift);
+
+               Lemma_Double_Shift
+                 (Single_Uns'Last, Single_Size - Shift_Prev, Shift);
+
+               if (Hi (Zu) and Mask) = 0 then
+                  Zu := Shift_Left (Zu, Shift);
+
+                  Prove_Shifting;
+                  pragma Assert (Big (Zu_Prev) =
+                    Big (Double_Uns'(abs Z)) * Big_2xx (Scale));
+                  Lemma_Shift_Without_Drop (Zu_Prev, Zu, Mask, Shift);
+                  Lemma_Substitution
+                    (Big (Zu), Big_2xx (Shift),
+                     Big (Zu_Prev), Big (Double_Uns'(abs Z)) * Big_2xx (Scale),
+                     0);
+                  Lemma_Powers_Of_2 (Shift, Scale);
+                  Lemma_Substitution
+                    (Big (Zu), Big (Double_Uns'(abs Z)),
+                     Big_2xx (Shift) * Big_2xx (Scale),
+                     Big_2xx (Shift + Scale), 0);
+                  Lemma_Double_Shift (abs Z, Scale, Shift);
+
+                  Scale := Scale + Shift;
+
+                  pragma Assert
+                    (Big (Zu) = Big (Double_Uns'(abs Z)) * Big_2xx (Scale));
+               end if;
+            end;
          end loop;
 
          Zhi := Hi (Zu);
          Zlo := Lo (Zu);
 
-         pragma Assert (Zhi /= 0);
-         --  We have Hi(Zu)/=0 before normalization. The sequence of Shift_Left
-         --  operations results in the leading bit of Zu being 1 by moving the
-         --  leftmost 1-bit in Zu to leading position, thus Zhi=Hi(Zu)/=0 here.
+         pragma Assert (Shift = 1);
+         pragma Assert (Mask = Shift_Left (Single_Uns'Last, Single_Size - 1));
+         pragma Assert ((Zhi and Mask) /= 0);
+         pragma Assert (Zhi >= 2 ** (Single_Size - 1));
+         pragma Assert (Big (Zu) = Big (Double_Uns'(abs Z)) * Big_2xx (Scale));
+         --  We have Hi (Zu) /= 0 before normalization. The sequence of
+         --  Shift_Left operations results in the leading bit of Zu being 1 by
+         --  moving the leftmost 1-bit in Zu to leading position, thus
+         --  Zhi = Hi (Zu) >= 2 ** (Single_Size - 1) here.
 
          --  Note that when we scale up the dividend, it still fits in four
          --  digits, since we already tested for overflow, and scaling does
          --  not change the invariant that (D (1) & D (2)) < Zu.
 
+         Lemma_Lt_Commutation (D (1) & D (2), abs Z);
+         Lemma_Lt_Mult (Big (D (1) & D (2)),
+                        Big (Double_Uns'(abs Z)), Big_2xx (Scale),
+                        Big_2xx64);
+
          T1 := Shift_Left (D (1) & D (2), Scale);
+         T2 := Shift_Left (Double_Uns (D (3)), Scale);
+         T3 := Shift_Left (Double_Uns (D (4)), Scale);
+
+         Prove_Dividend_Scaling;
+
          D (1) := Hi (T1);
-         T2 := Shift_Left (0 & D (3), Scale);
          D (2) := Lo (T1) or Hi (T2);
-         T3 := Shift_Left (0 & D (4), Scale);
          D (3) := Lo (T2) or Hi (T3);
          D (4) := Lo (T3);
 
-         --  Loop to compute quotient digits, runs twice for Qd(1) and Qd(2)
-
-         for J in 0 .. 1 loop
-
-            --  Compute next quotient digit. We have to divide three digits by
-            --  two digits. We estimate the quotient by dividing the leading
-            --  two digits by the leading digit. Given the scaling we did above
-            --  which ensured the first bit of the divisor is set, this gives
-            --  an estimate of the quotient that is at most two too high.
-
-            Qd (J + 1) := (if D (J + 1) = Zhi
-                           then 2 ** Single_Size - 1
-                           else Lo ((D (J + 1) & D (J + 2)) / Zhi));
-
-            --  Compute amount to subtract
-
-            T1 := Qd (J + 1) * Zlo;
-            T2 := Qd (J + 1) * Zhi;
-            S3 := Lo (T1);
-            T1 := Hi (T1) + Lo (T2);
-            S2 := Lo (T1);
-            S1 := Hi (T1) + Hi (T2);
-
-            --  Adjust quotient digit if it was too high
-
-            --  We use the version of the algorithm in the 2nd Edition of
-            --  "The Art of Computer Programming". This had a bug not
-            --  discovered till 1995, see Vol 2 errata:
-            --     http://www-cs-faculty.stanford.edu/~uno/err2-2e.ps.gz.
-            --  Under rare circumstances the expression in the test could
-            --  overflow. This version was further corrected in 2005, see
-            --  Vol 2 errata:
-            --     http://www-cs-faculty.stanford.edu/~uno/all2-pre.ps.gz.
-            --  This implementation is not impacted by these bugs, due to the
-            --  use of a word-size comparison done in function Le3 instead of
-            --  a comparison on two-word integer quantities in the original
-            --  algorithm.
-
-            loop
-               exit when Le3 (S1, S2, S3, D (J + 1), D (J + 2), D (J + 3));
-               Qd (J + 1) := Qd (J + 1) - 1;
-               Sub3 (S1, S2, S3, 0, Zhi, Zlo);
+         pragma Assert (Mult * Big_2xx (Scale) =
+           Big_2xx32 * Big_2xx32 * Big_2xx32 * Big (Double_Uns (D (1)))
+                     + Big_2xx32 * Big_2xx32 * Big (Double_Uns (D (2)))
+                                 + Big_2xx32 * Big (Double_Uns (D (3)))
+                                             + Big (Double_Uns (D (4))));
+         Lemma_Substitution (Big_2xx64 * Big (Zu), Big_2xx64, Big (Zu),
+                             Big (Double_Uns'(abs Z)) * Big_2xx (Scale), 0);
+         Lemma_Lt_Mult (Mult, Big_2xx64 * Big (Double_Uns'(abs Z)),
+                        Big_2xx (Scale), Big_2xx64 * Big (Zu));
+         Lemma_Div_Lt (Mult * Big_2xx (Scale), Big (Zu), Big_2xx64);
+         Lemma_Substitution (Mult * Big_2xx (Scale), Big_2xx32,
+                             Big_2xx32 * Big_2xx32 * Big (Double_Uns (D (1)))
+                                       + Big_2xx32 * Big (Double_Uns (D (2)))
+                                                   + Big (Double_Uns (D (3))),
+                             Big3 (D (1), D (2), D (3)),
+                             Big (Double_Uns (D (4))));
+
+         --  Loop to compute quotient digits, runs twice for Qd (1) and Qd (2)
+
+         declare
+            Qd1 : Single_Uns := 0 with Ghost;
+         begin
+            for J in 1 .. 2 loop
+               Lemma_Hi_Lo (D (J) & D (J + 1), D (J), D (J + 1));
+               pragma Assert (Big (D (J) & D (J + 1)) < Big (Zu));
+
+               --  Compute next quotient digit. We have to divide three digits
+               --  by two digits. We estimate the quotient by dividing the
+               --  leading two digits by the leading digit. Given the scaling
+               --  we did above which ensured the first bit of the divisor is
+               --  set, this gives an estimate of the quotient that is at most
+               --  two too high.
+
+               if D (J) > Zhi then
+                  Lemma_Lt_Commutation (Zu, D (J) & D (J + 1));
+                  pragma Assert (False);
+
+               elsif D (J) = Zhi then
+                  Qd (J) := Single_Uns'Last;
+
+                  Lemma_Gt_Mult (Big (Zu), Big (D (J) & D (J + 1)) + 1,
+                                 Big_2xx32,
+                                 Big3 (D (J), D (J + 1), D (J + 2)));
+                  Lemma_Div_Lt
+                    (Big3 (D (J), D (J + 1), D (J + 2)), Big_2xx32, Big (Zu));
+
+               else
+                  Qd (J) := Lo ((D (J) & D (J + 1)) / Zhi);
+
+                  Prove_Qd_Calculation_Part_1 (J);
+               end if;
+
+               Lemma_Gt_Mult
+                 (Big (Double_Uns (Qd (J))),
+                  Big3 (D (J), D (J + 1), D (J + 2)) / Big (Zu),
+                  Big (Zu), Big3 (D (J), D (J + 1), D (J + 2)) - Big (Zu));
+
+               --  Compute amount to subtract
+
+               T1 := Qd (J) * Zlo;
+               T2 := Qd (J) * Zhi;
+               S3 := Lo (T1);
+               T3 := Hi (T1) + Lo (T2);
+               S2 := Lo (T3);
+               S1 := Hi (T3) + Hi (T2);
+
+               Prove_Multiplication (Qd (J));
+
+               --  Adjust quotient digit if it was too high
+
+               --  We use the version of the algorithm in the 2nd Edition
+               --  of "The Art of Computer Programming". This had a bug not
+               --  discovered till 1995, see Vol 2 errata:
+               --     http://www-cs-faculty.stanford.edu/~uno/err2-2e.ps.gz.
+               --  Under rare circumstances the expression in the test could
+               --  overflow. This version was further corrected in 2005, see
+               --  Vol 2 errata:
+               --     http://www-cs-faculty.stanford.edu/~uno/all2-pre.ps.gz.
+               --  This implementation is not impacted by these bugs, due
+               --  to the use of a word-size comparison done in function Le3
+               --  instead of a comparison on two-word integer quantities in
+               --  the original algorithm.
+
+               Lemma_Hi_Lo_3 (Zu, Zhi, Zlo);
+
+               while not Le3 (S1, S2, S3, D (J), D (J + 1), D (J + 2)) loop
+                  pragma Loop_Invariant (Qd (J)'Initialized);
+                  pragma Loop_Invariant
+                    (Big3 (S1, S2, S3) = Big (Double_Uns (Qd (J))) * Big (Zu));
+                  pragma Assert (Big3 (S1, S2, S3) > 0);
+                  if Qd (J) = 0 then
+                     pragma Assert (Big3 (S1, S2, S3) = 0);
+                     pragma Assert (False);
+                  end if;
+                  Lemma_Ge_Commutation (Double_Uns (Qd (J)), 1);
+                  Lemma_Ge_Mult
+                    (Big (Double_Uns (Qd (J))), 1, Big (Zu), Big (Zu));
+
+                  Sub3 (S1, S2, S3, 0, Zhi, Zlo);
+
+                  pragma Assert
+                    (Big3 (S1, S2, S3) >
+                       Big3 (D (J), D (J + 1), D (J + 2)) - Big (Zu));
+                  Lemma_Subtract_Commutation (Double_Uns (Qd (J)), 1);
+                  Lemma_Substitution (Big3 (S1, S2, S3), Big (Zu),
+                                      Big (Double_Uns (Qd (J))) - 1,
+                                      Big (Double_Uns (Qd (J) - 1)), 0);
+
+                  Qd (J) := Qd (J) - 1;
+
+                  pragma Assert
+                    (Big3 (S1, S2, S3) = Big (Double_Uns (Qd (J))) * Big (Zu));
+               end loop;
+
+               --  Now subtract S1&S2&S3 from D1&D2&D3 ready for next step
+
+               pragma Assert
+                 (Big3 (S1, S2, S3) = Big (Double_Uns (Qd (J))) * Big (Zu));
+               pragma Assert (Big3 (S1, S2, S3) >
+                                Big3 (D (J), D (J + 1), D (J + 2)) - Big (Zu));
+               Inline_Le3 (S1, S2, S3, D (J), D (J + 1), D (J + 2));
+
+               Sub3 (D (J), D (J + 1), D (J + 2), S1, S2, S3);
+
+               pragma Assert (Big3 (D (J), D (J + 1), D (J + 2)) < Big (Zu));
+               if D (J) > 0 then
+                  pragma Assert (Big_2xx32 * Big_2xx32 = Big_2xx64);
+                  pragma Assert (Big3 (D (J), D (J + 1), D (J + 2)) =
+                                   Big_2xx32 * Big_2xx32
+                                   * Big (Double_Uns (D (J)))
+                                 + Big_2xx32 * Big (Double_Uns (D (J + 1)))
+                                             + Big (Double_Uns (D (J + 2))));
+                  pragma Assert (Big3 (D (J), D (J + 1), D (J + 2)) =
+                                   Big_2xx64 * Big (Double_Uns (D (J)))
+                                 + Big_2xx32 * Big (Double_Uns (D (J + 1)))
+                                             + Big (Double_Uns (D (J + 2))));
+                  pragma Assert (Big3 (D (J), D (J + 1), D (J + 2)) >=
+                                   Big_2xx64 * Big (Double_Uns (D (J))));
+                  Lemma_Ge_Commutation (Double_Uns (D (J)), Double_Uns'(1));
+                  pragma Assert
+                    (Big3 (D (J), D (J + 1), D (J + 2)) >= Big_2xx64);
+                  pragma Assert (False);
+               end if;
+
+               if J = 1 then
+                  Qd1 := Qd (1);
+                  pragma Assert
+                    (Big_2xx32 * Big_2xx32 * Big (Double_Uns (D (1))) = 0);
+                  pragma Assert
+                    (Mult * Big_2xx (Scale) =
+                       Big_2xx32 * Big3 (S1, S2, S3)
+                     + Big3 (D (2), D (3), D (4)));
+                  Lemma_Substitution (Mult * Big_2xx (Scale), Big_2xx32,
+                                      Big3 (S1, S2, S3),
+                                      Big (Double_Uns (Qd1)) * Big (Zu),
+                                      Big3 (D (2), D (3), D (4)));
+               else
+                  pragma Assert (Qd1 = Qd (1));
+                  pragma Assert
+                    (Big_2xx32 * Big_2xx32 * Big (Double_Uns (D (2))) = 0);
+                  pragma Assert
+                    (Mult * Big_2xx (Scale) =
+                       Big_2xx32 * Big (Double_Uns (Qd (1))) * Big (Zu)
+                     + Big3 (S1, S2, S3)
+                     + Big3 (D (2), D (3), D (4)));
+                  pragma Assert
+                    (Mult * Big_2xx (Scale) =
+                       Big_2xx32 * Big (Double_Uns (Qd (1))) * Big (Zu)
+                                 + Big (Double_Uns (Qd (2))) * Big (Zu)
+                     + Big_2xx32 * Big (Double_Uns (D (3)))
+                                 + Big (Double_Uns (D (4))));
+               end if;
             end loop;
-
-            --  Now subtract S1&S2&S3 from D1&D2&D3 ready for next step
-
-            Sub3 (D (J + 1), D (J + 2), D (J + 3), S1, S2, S3);
-         end loop;
+         end;
 
          --  The two quotient digits are now set, and the remainder of the
          --  scaled division is in D3&D4. To get the remainder for the
@@ -564,38 +2469,93 @@ package body System.Arith_Double is
          --  for rounding below.
 
          Qu := Qd (1) & Qd (2);
-         Ru := Shift_Right (D (3) & D (4), Scale);
+         Ru := D (3) & D (4);
+
+         pragma Assert
+           (Mult * Big_2xx (Scale) =
+              Big_2xx32 * Big (Double_Uns (Qd (1))) * Big (Zu)
+                        + Big (Double_Uns (Qd (2))) * Big (Zu)
+            + Big_2xx32 * Big (Double_Uns (D (3)))
+                        + Big (Double_Uns (D (4))));
+         Lemma_Hi_Lo (Qu, Qd (1), Qd (2));
+         Lemma_Hi_Lo (Ru, D (3), D (4));
+         Lemma_Substitution
+           (Mult * Big_2xx (Scale), Big (Zu),
+            Big_2xx32 * Big (Double_Uns (Qd (1))) + Big (Double_Uns (Qd (2))),
+            Big (Qu), Big (Ru));
+         Prove_Rescaling;
+
+         Ru := Shift_Right (Ru, Scale);
+
+         Lemma_Shift_Right (Zu, Scale);
+
          Zu := Shift_Right (Zu, Scale);
+
+         Lemma_Simplify (Big (Double_Uns'(abs Z)), Big_2xx (Scale));
       end if;
 
+      pragma Assert (Big (Ru) = abs Big_R);
+      pragma Assert (Big (Qu) = abs Quot);
+      pragma Assert (Big (Zu) = Big (Double_Uns'(abs Z)));
+
       --  Deal with rounding case
 
-      if Round and then Ru > (Zu - Double_Uns'(1)) / Double_Uns'(2) then
+      if Round then
+         Prove_Rounding_Case;
 
-         --  Protect against wrapping around when rounding, by signaling
-         --  an overflow when the quotient is too large.
+         if Ru > (Zu - Double_Uns'(1)) / Double_Uns'(2) then
+            pragma Assert (abs Big_Q = Big (Qu) + 1);
 
-         if Qu = Double_Uns'Last then
-            Raise_Error;
-         end if;
+            --  Protect against wrapping around when rounding, by signaling
+            --  an overflow when the quotient is too large.
+
+            if Qu = Double_Uns'Last then
+               pragma Assert (abs Big_Q = Big_2xx64);
+               Raise_Error;
+               pragma Annotate
+                 (GNATprove, Intentional, "call to nonreturning subprogram",
+                  "Constraint_Error is raised in case of overflow");
+            end if;
+
+            Lemma_Add_One (Qu);
 
-         Qu := Qu + Double_Uns'(1);
+            Qu := Qu + Double_Uns'(1);
+         end if;
       end if;
 
+      pragma Assert (Big (Qu) = abs Big_Q);
+
       --  Set final signs (RM 4.5.5(27-30))
 
       --  Case of dividend (X * Y) sign positive
 
       if (X >= 0 and then Y >= 0) or else (X < 0 and then Y < 0) then
          R := To_Pos_Int (Ru);
+         pragma Annotate
+           (GNATprove, Intentional, "precondition",
+            "Constraint_Error is raised in case of overflow");
+
          Q := (if Z > 0 then To_Pos_Int (Qu) else To_Neg_Int (Qu));
+         pragma Annotate
+           (GNATprove, Intentional, "precondition",
+            "Constraint_Error is raised in case of overflow");
 
       --  Case of dividend (X * Y) sign negative
 
       else
          R := To_Neg_Int (Ru);
+         pragma Annotate
+           (GNATprove, Intentional, "precondition",
+            "Constraint_Error is raised in case of overflow");
+
          Q := (if Z > 0 then To_Neg_Int (Qu) else To_Pos_Int (Qu));
+         pragma Annotate
+           (GNATprove, Intentional, "precondition",
+            "Constraint_Error is raised in case of overflow");
       end if;
+
+      Prove_Sign_R;
+      Prove_Signs;
    end Scaled_Divide;
 
    ----------
@@ -603,23 +2563,178 @@ package body System.Arith_Double is
    ----------
 
    procedure Sub3 (X1, X2, X3 : in out Single_Uns; Y1, Y2, Y3 : Single_Uns) is
+
+      --  Local ghost variables
+
+      XX1 : constant Single_Uns := X1 with Ghost;
+      XX2 : constant Single_Uns := X2 with Ghost;
+      XX3 : constant Single_Uns := X3 with Ghost;
+
+      --  Local lemmas
+
+      procedure Lemma_Add3_No_Carry (X1, X2, X3, Y1, Y2, Y3 : Single_Uns)
+      with
+        Ghost,
+        Pre  => X1 <= Single_Uns'Last - Y1
+          and then X2 <= Single_Uns'Last - Y2
+          and then X3 <= Single_Uns'Last - Y3,
+        Post => Big3 (X1 + Y1, X2 + Y2, X3 + Y3)
+              = Big3 (X1, X2, X3) + Big3 (Y1, Y2, Y3);
+
+      procedure Lemma_Ge_Expand (X1, X2, X3, Y1, Y2, Y3 : Single_Uns)
+      with
+        Ghost,
+        Pre  => Big3 (X1, X2, X3) >= Big3 (Y1, Y2, Y3),
+        Post => X1 > Y1
+          or else (X1 = Y1 and then X2 > Y2)
+          or else (X1 = Y1 and then X2 = Y2 and then X3 >= Y3);
+
+      procedure Lemma_Sub3_No_Carry (X1, X2, X3, Y1, Y2, Y3 : Single_Uns)
+      with
+        Ghost,
+        Pre  => X1 >= Y1 and then X2 >= Y2 and then X3 >= Y3,
+        Post => Big3 (X1 - Y1, X2 - Y2, X3 - Y3)
+              = Big3 (X1, X2, X3) - Big3 (Y1, Y2, Y3);
+
+      procedure Lemma_Sub3_With_Carry2 (X1, X2, X3, Y2 : Single_Uns)
+      with
+        Ghost,
+        Pre  => X2 < Y2,
+        Post => Big3 (X1, X2 - Y2, X3)
+              = Big3 (X1, X2, X3) + Big3 (1, 0, 0) - Big3 (0, Y2, 0);
+
+      procedure Lemma_Sub3_With_Carry3 (X1, X2, X3, Y3 : Single_Uns)
+      with
+        Ghost,
+        Pre  => X3 < Y3,
+        Post => Big3 (X1, X2, X3 - Y3)
+              = Big3 (X1, X2, X3) + Big3 (0, 1, 0) - Big3 (0, 0, Y3);
+
+      -------------------------
+      -- Lemma_Add3_No_Carry --
+      -------------------------
+
+      procedure Lemma_Add3_No_Carry (X1, X2, X3, Y1, Y2, Y3 : Single_Uns) is
+      begin
+         Lemma_Add_Commutation (Double_Uns (X1), Y1);
+         Lemma_Add_Commutation (Double_Uns (X2), Y2);
+         Lemma_Add_Commutation (Double_Uns (X3), Y3);
+         pragma Assert (Double_Uns (Single_Uns'(X1 + Y1))
+                        = Double_Uns (X1) + Double_Uns (Y1));
+         pragma Assert (Double_Uns (Single_Uns'(X2 + Y2))
+                        = Double_Uns (X2) + Double_Uns (Y2));
+         pragma Assert (Double_Uns (Single_Uns'(X3 + Y3))
+                        = Double_Uns (X3) + Double_Uns (Y3));
+      end Lemma_Add3_No_Carry;
+
+      ---------------------
+      -- Lemma_Ge_Expand --
+      ---------------------
+
+      procedure Lemma_Ge_Expand (X1, X2, X3, Y1, Y2, Y3 : Single_Uns) is null;
+
+      -------------------------
+      -- Lemma_Sub3_No_Carry --
+      -------------------------
+
+      procedure Lemma_Sub3_No_Carry (X1, X2, X3, Y1, Y2, Y3 : Single_Uns) is
+      begin
+         Lemma_Subtract_Commutation (Double_Uns (X1), Double_Uns (Y1));
+         Lemma_Subtract_Commutation (Double_Uns (X2), Double_Uns (Y2));
+         Lemma_Subtract_Commutation (Double_Uns (X3), Double_Uns (Y3));
+      end Lemma_Sub3_No_Carry;
+
+      ----------------------------
+      -- Lemma_Sub3_With_Carry2 --
+      ----------------------------
+
+      procedure Lemma_Sub3_With_Carry2 (X1, X2, X3, Y2 : Single_Uns) is
+         pragma Unreferenced (X1, X3);
+      begin
+         Lemma_Add_Commutation
+           (Double_Uns'(2 ** Single_Size) - Double_Uns (Y2), X2);
+         Lemma_Subtract_Commutation
+           (Double_Uns'(2 ** Single_Size), Double_Uns (Y2));
+      end Lemma_Sub3_With_Carry2;
+
+      ----------------------------
+      -- Lemma_Sub3_With_Carry3 --
+      ----------------------------
+
+      procedure Lemma_Sub3_With_Carry3 (X1, X2, X3, Y3 : Single_Uns) is
+         pragma Unreferenced (X1, X2);
+      begin
+         Lemma_Add_Commutation
+           (Double_Uns'(2 ** Single_Size) - Double_Uns (Y3), X3);
+         Lemma_Subtract_Commutation
+           (Double_Uns'(2 ** Single_Size), Double_Uns (Y3));
+      end Lemma_Sub3_With_Carry3;
+
+   --  Start of processing for Sub3
+
    begin
+      Lemma_Ge_Expand (X1, X2, X3, Y1, Y2, Y3);
+
       if Y3 > X3 then
          if X2 = 0 then
+            pragma Assert (X1 >= 1);
+            Lemma_Sub3_No_Carry (X1, X2, X3, 1, 0, 0);
+
             X1 := X1 - 1;
+
+            pragma Assert
+              (Big3 (X1, X2, X3) = Big3 (XX1, XX2, XX3) - Big3 (1, 0, 0));
+            pragma Assert
+              (Big3 (X1, X2, X3) = Big3 (XX1, XX2, XX3)
+               - Big3 (0, Single_Uns'Last, 0) - Big3 (0, 1, 0));
+            Lemma_Add3_No_Carry (X1, X2, X3, 0, Single_Uns'Last, 0);
+         else
+            Lemma_Sub3_No_Carry (X1, X2, X3, 0, 1, 0);
          end if;
 
          X2 := X2 - 1;
+
+         pragma Assert
+           (Big3 (X1, X2, X3) = Big3 (XX1, XX2, XX3) - Big3 (0, 1, 0));
+         Lemma_Sub3_With_Carry3 (X1, X2, X3, Y3);
+      else
+         Lemma_Sub3_No_Carry (X1, X2, X3, 0, 0, Y3);
       end if;
 
       X3 := X3 - Y3;
 
+      pragma Assert
+        (Big3 (X1, X2, X3) = Big3 (XX1, XX2, XX3) - Big3 (0, 0, Y3));
+
       if Y2 > X2 then
+         pragma Assert (X1 >= 1);
+         Lemma_Sub3_No_Carry (X1, X2, X3, 1, 0, 0);
+
          X1 := X1 - 1;
+
+         pragma Assert
+           (Big3 (X1, X2, X3) = Big3 (XX1, XX2, XX3)
+            - Big3 (0, 0, Y3) - Big3 (1, 0, 0));
+         Lemma_Sub3_With_Carry2 (X1, X2, X3, Y2);
+      else
+         Lemma_Sub3_No_Carry (X1, X2, X3, 0, Y2, 0);
       end if;
 
       X2 := X2 - Y2;
+
+      pragma Assert
+        (Big3 (X1, X2, X3) = Big3 (XX1, XX2, XX3) - Big3 (0, Y2, Y3));
+      pragma Assert (X1 >= Y1);
+      Lemma_Sub3_No_Carry (X1, Y2, X3, Y1, 0, 0);
+
       X1 := X1 - Y1;
+
+      pragma Assert
+        (Big3 (X1, X2, X3) = Big3 (XX1, XX2, XX3)
+         - Big3 (0, Y2, Y3) - Big3 (Y1, 0, 0));
+      Lemma_Add3_No_Carry (0, Y2, Y3, Y1, 0, 0);
+      pragma Assert
+        (Big3 (X1, X2, X3) = Big3 (XX1, XX2, XX3) - Big3 (Y1, Y2, Y3));
    end Sub3;
 
    -------------------------------
@@ -629,18 +2744,126 @@ package body System.Arith_Double is
    function Subtract_With_Ovflo_Check (X, Y : Double_Int) return Double_Int is
       R : constant Double_Int := To_Int (To_Uns (X) - To_Uns (Y));
 
+      --  Local lemmas
+
+      procedure Prove_Negative_X
+      with
+        Ghost,
+        Pre  => X < 0 and then (Y <= 0 or else R < 0),
+        Post => R = X - Y;
+
+      procedure Prove_Non_Negative_X
+      with
+        Ghost,
+        Pre  => X >= 0 and then (Y > 0 or else R >= 0),
+        Post => R = X - Y;
+
+      procedure Prove_Overflow_Case
+      with
+        Ghost,
+        Pre  =>
+          (if X >= 0 then Y <= 0 and then R < 0
+                     else Y > 0 and then R >= 0),
+        Post => not In_Double_Int_Range (Big (X) - Big (Y));
+
+      ----------------------
+      -- Prove_Negative_X --
+      ----------------------
+
+      procedure Prove_Negative_X is
+      begin
+         if X = Double_Int'First then
+            if Y = Double_Int'First or else Y > 0 then
+               null;
+            else
+               pragma Assert
+                 (To_Uns (X) - To_Uns (Y) =
+                    2 ** (Double_Size - 1) + Double_Uns (-Y));
+            end if;
+
+         elsif Y >= 0 or else Y = Double_Int'First then
+            null;
+
+         else
+            pragma Assert
+              (To_Uns (X) - To_Uns (Y) = -Double_Uns (-X) + Double_Uns (-Y));
+         end if;
+      end Prove_Negative_X;
+
+      --------------------------
+      -- Prove_Non_Negative_X --
+      --------------------------
+
+      procedure Prove_Non_Negative_X is
+      begin
+         if Y > 0 then
+            declare
+               Ru : constant Double_Uns := To_Uns (X) - To_Uns (Y);
+            begin
+               pragma Assert (Ru = Double_Uns (X) - Double_Uns (Y));
+               if Ru < 2 ** (Double_Size - 1) then  --  R >= 0
+                  Lemma_Subtract_Double_Uns (X => Y, Y => X);
+                  pragma Assert (Ru = Double_Uns (X - Y));
+
+               elsif Ru = 2 ** (Double_Size - 1) then
+                  pragma Assert (Double_Uns (Y) < 2 ** (Double_Size - 1));
+                  pragma Assert (False);
+
+               else
+                  pragma Assert
+                    (R = -Double_Int (-(Double_Uns (X) - Double_Uns (Y))));
+                  pragma Assert
+                    (R = -Double_Int (-Double_Uns (X) + Double_Uns (Y)));
+                  pragma Assert
+                    (R = -Double_Int (Double_Uns (Y) - Double_Uns (X)));
+               end if;
+            end;
+
+         elsif Y = Double_Int'First then
+            pragma Assert
+              (To_Uns (X) - To_Uns (Y) =
+                 Double_Uns (X) - 2 ** (Double_Size - 1));
+            pragma Assert (False);
+
+         else
+            pragma Assert
+              (To_Uns (X) - To_Uns (Y) = Double_Uns (X) + Double_Uns (-Y));
+         end if;
+      end Prove_Non_Negative_X;
+
+      -------------------------
+      -- Prove_Overflow_Case --
+      -------------------------
+
+      procedure Prove_Overflow_Case is
+      begin
+         if X >= 0 and then Y /= Double_Int'First then
+            pragma Assert
+              (To_Uns (X) - To_Uns (Y) = Double_Uns (X) + Double_Uns (-Y));
+
+         elsif X < 0 and then X /= Double_Int'First then
+            pragma Assert
+              (To_Uns (X) - To_Uns (Y) = -Double_Uns (-X) - Double_Uns (Y));
+         end if;
+      end Prove_Overflow_Case;
+
+   --  Start of processing for Subtract_With_Ovflo_Check
+
    begin
       if X >= 0 then
          if Y > 0 or else R >= 0 then
+            Prove_Non_Negative_X;
             return R;
          end if;
 
       else -- X < 0
          if Y <= 0 or else R < 0 then
+            Prove_Negative_X;
             return R;
          end if;
       end if;
 
+      Prove_Overflow_Case;
       Raise_Error;
    end Subtract_With_Ovflo_Check;
 
diff --git a/gcc/ada/libgnat/s-aridou.ads b/gcc/ada/libgnat/s-aridou.ads
index 0df27ca..7c32f7f 100644
--- a/gcc/ada/libgnat/s-aridou.ads
+++ b/gcc/ada/libgnat/s-aridou.ads
@@ -33,6 +33,9 @@
 --  signed integer values in cases where either overflow checking is required,
 --  or intermediate results are longer than the result type.
 
+with Ada.Numerics.Big_Numbers.Big_Integers_Ghost;
+use Ada.Numerics.Big_Numbers.Big_Integers_Ghost;
+
 generic
 
    type Double_Int is range <>;
@@ -50,27 +53,93 @@ generic
    with function Shift_Left (A : Single_Uns; B : Natural) return Single_Uns
      is <>;
 
-package System.Arith_Double is
-   pragma Pure;
+package System.Arith_Double
+  with Pure, SPARK_Mode
+is
+   --  Preconditions in this unit are meant for analysis only, not for run-time
+   --  checking, so that the expected exceptions are raised. This is enforced
+   --  by setting the corresponding assertion policy to Ignore. Postconditions
+   --  and contract cases should not be executed at runtime as well, in order
+   --  not to slow down the execution of these functions.
+
+   pragma Assertion_Policy (Pre            => Ignore,
+                            Post           => Ignore,
+                            Contract_Cases => Ignore,
+                            Ghost          => Ignore);
+
+   package Signed_Conversion is new Signed_Conversions (Int => Double_Int);
+
+   function Big (Arg : Double_Int) return Big_Integer is
+     (Signed_Conversion.To_Big_Integer (Arg))
+   with Ghost;
+
+   package Unsigned_Conversion is new Unsigned_Conversions (Int => Double_Uns);
 
-   function Add_With_Ovflo_Check (X, Y : Double_Int) return Double_Int;
+   function Big (Arg : Double_Uns) return Big_Integer is
+     (Unsigned_Conversion.To_Big_Integer (Arg))
+   with Ghost;
+
+   function In_Double_Int_Range (Arg : Big_Integer) return Boolean is
+     (In_Range (Arg, Big (Double_Int'First), Big (Double_Int'Last)))
+   with Ghost;
+
+   function Add_With_Ovflo_Check (X, Y : Double_Int) return Double_Int
+   with
+     Pre  => In_Double_Int_Range (Big (X) + Big (Y)),
+     Post => Add_With_Ovflo_Check'Result = X + Y;
    --  Raises Constraint_Error if sum of operands overflows Double_Int,
    --  otherwise returns the signed integer sum.
 
-   function Subtract_With_Ovflo_Check (X, Y : Double_Int) return Double_Int;
+   function Subtract_With_Ovflo_Check (X, Y : Double_Int) return Double_Int
+   with
+     Pre  => In_Double_Int_Range (Big (X) - Big (Y)),
+     Post => Subtract_With_Ovflo_Check'Result = X - Y;
    --  Raises Constraint_Error if difference of operands overflows Double_Int,
    --  otherwise returns the signed integer difference.
 
-   function Multiply_With_Ovflo_Check (X, Y : Double_Int) return Double_Int;
+   function Multiply_With_Ovflo_Check (X, Y : Double_Int) return Double_Int
+   with
+     Pre  => In_Double_Int_Range (Big (X) * Big (Y)),
+     Post => Multiply_With_Ovflo_Check'Result = X * Y;
    pragma Convention (C, Multiply_With_Ovflo_Check);
    --  Raises Constraint_Error if product of operands overflows Double_Int,
    --  otherwise returns the signed integer product. Gigi may also call this
    --  routine directly.
 
+   function Same_Sign (X, Y : Big_Integer) return Boolean is
+     (X = Big (Double_Int'(0))
+        or else Y = Big (Double_Int'(0))
+        or else (X < Big (Double_Int'(0))) = (Y < Big (Double_Int'(0))))
+   with Ghost;
+
+   function Round_Quotient (X, Y, Q, R : Big_Integer) return Big_Integer is
+     (if abs R > (abs Y - Big (Double_Int'(1))) / Big (Double_Int'(2)) then
+       (if Same_Sign (X, Y) then Q + Big (Double_Int'(1))
+        else Q - Big (Double_Int'(1)))
+      else
+        Q)
+   with
+     Ghost,
+     Pre => Y /= 0 and then Q = X / Y and then R = X rem Y;
+
    procedure Scaled_Divide
      (X, Y, Z : Double_Int;
       Q, R    : out Double_Int;
-      Round   : Boolean);
+      Round   : Boolean)
+   with
+     Pre  => Z /= 0
+       and then In_Double_Int_Range
+         (if Round then Round_Quotient (Big (X) * Big (Y), Big (Z),
+                                        Big (X) * Big (Y) / Big (Z),
+                                        Big (X) * Big (Y) rem Big (Z))
+          else Big (X) * Big (Y) / Big (Z)),
+     Post => Big (R) = Big (X) * Big (Y) rem Big (Z)
+       and then
+         (if Round then
+            Big (Q) = Round_Quotient (Big (X) * Big (Y), Big (Z),
+                                      Big (X) * Big (Y) / Big (Z), Big (R))
+          else
+            Big (Q) = Big (X) * Big (Y) / Big (Z));
    --  Performs the division of (X * Y) / Z, storing the quotient in Q
    --  and the remainder in R. Constraint_Error is raised if Z is zero,
    --  or if the quotient does not fit in Double_Int. Round indicates if
@@ -82,7 +151,22 @@ package System.Arith_Double is
    procedure Double_Divide
      (X, Y, Z : Double_Int;
       Q, R    : out Double_Int;
-      Round   : Boolean);
+      Round   : Boolean)
+   with
+     Pre  => Y /= 0
+       and then Z /= 0
+       and then In_Double_Int_Range
+         (if Round then Round_Quotient (Big (X), Big (Y) * Big (Z),
+                                        Big (X) / (Big (Y) * Big (Z)),
+                                        Big (X) rem (Big (Y) * Big (Z)))
+          else Big (X) / (Big (Y) * Big (Z))),
+     Post => Big (R) = Big (X) rem (Big (Y) * Big (Z))
+       and then
+         (if Round then
+            Big (Q) = Round_Quotient (Big (X), Big (Y) * Big (Z),
+                                      Big (X) / (Big (Y) * Big (Z)), Big (R))
+          else
+            Big (Q) = Big (X) / (Big (Y) * Big (Z)));
    --  Performs the division X / (Y * Z), storing the quotient in Q and
    --  the remainder in R. Constraint_Error is raised if Y or Z is zero,
    --  or if the quotient does not fit in Double_Int. Round indicates if the
diff --git a/gcc/ada/libgnat/s-arit64.adb b/gcc/ada/libgnat/s-arit64.adb
index 2f24a70..23b76d9 100644
--- a/gcc/ada/libgnat/s-arit64.adb
+++ b/gcc/ada/libgnat/s-arit64.adb
@@ -31,7 +31,9 @@
 
 with System.Arith_Double;
 
-package body System.Arith_64 is
+package body System.Arith_64
+  with SPARK_Mode
+is
 
    subtype Uns64 is Interfaces.Unsigned_64;
    subtype Uns32 is Interfaces.Unsigned_32;
diff --git a/gcc/ada/libgnat/s-arit64.ads b/gcc/ada/libgnat/s-arit64.ads
index c9141f5..fbfd0f6 100644
--- a/gcc/ada/libgnat/s-arit64.ads
+++ b/gcc/ada/libgnat/s-arit64.ads
@@ -36,31 +36,103 @@
 pragma Restrictions (No_Elaboration_Code);
 --  Allow direct call from gigi generated code
 
+--  Preconditions in this unit are meant for analysis only, not for run-time
+--  checking, so that the expected exceptions are raised. This is enforced
+--  by setting the corresponding assertion policy to Ignore. Postconditions
+--  and contract cases should not be executed at runtime as well, in order
+--  not to slow down the execution of these functions.
+
+pragma Assertion_Policy (Pre            => Ignore,
+                         Post           => Ignore,
+                         Contract_Cases => Ignore,
+                         Ghost          => Ignore);
+
+with Ada.Numerics.Big_Numbers.Big_Integers_Ghost;
 with Interfaces;
 
-package System.Arith_64 is
-   pragma Pure;
+package System.Arith_64
+  with Pure, SPARK_Mode
+is
+   use type Ada.Numerics.Big_Numbers.Big_Integers_Ghost.Big_Integer;
+   use type Interfaces.Integer_64;
 
    subtype Int64 is Interfaces.Integer_64;
 
-   function Add_With_Ovflo_Check64 (X, Y : Int64) return Int64;
+   subtype Big_Integer is
+     Ada.Numerics.Big_Numbers.Big_Integers_Ghost.Big_Integer
+     with Ghost;
+
+   package Signed_Conversion is new
+     Ada.Numerics.Big_Numbers.Big_Integers_Ghost.Signed_Conversions
+     (Int => Int64);
+
+   function Big (Arg : Int64) return Big_Integer is
+     (Signed_Conversion.To_Big_Integer (Arg))
+   with Ghost;
+
+   function In_Int64_Range (Arg : Big_Integer) return Boolean is
+     (Ada.Numerics.Big_Numbers.Big_Integers_Ghost.In_Range
+       (Arg, Big (Int64'First), Big (Int64'Last)))
+   with Ghost;
+
+   function Add_With_Ovflo_Check64 (X, Y : Int64) return Int64
+   with
+     Pre  => In_Int64_Range (Big (X) + Big (Y)),
+     Post => Add_With_Ovflo_Check64'Result = X + Y;
    --  Raises Constraint_Error if sum of operands overflows 64 bits,
    --  otherwise returns the 64-bit signed integer sum.
 
-   function Subtract_With_Ovflo_Check64 (X, Y : Int64) return Int64;
+   function Subtract_With_Ovflo_Check64 (X, Y : Int64) return Int64
+   with
+     Pre  => In_Int64_Range (Big (X) - Big (Y)),
+     Post => Subtract_With_Ovflo_Check64'Result = X - Y;
    --  Raises Constraint_Error if difference of operands overflows 64
    --  bits, otherwise returns the 64-bit signed integer difference.
 
-   function Multiply_With_Ovflo_Check64 (X, Y : Int64) return Int64;
+   function Multiply_With_Ovflo_Check64 (X, Y : Int64) return Int64
+   with
+     Pure_Function,
+     Pre  => In_Int64_Range (Big (X) * Big (Y)),
+     Post => Multiply_With_Ovflo_Check64'Result = X * Y;
    pragma Export (C, Multiply_With_Ovflo_Check64, "__gnat_mulv64");
    --  Raises Constraint_Error if product of operands overflows 64
    --  bits, otherwise returns the 64-bit signed integer product.
    --  Gigi may also call this routine directly.
 
+   function Same_Sign (X, Y : Big_Integer) return Boolean is
+     (X = Big (Int64'(0))
+        or else Y = Big (Int64'(0))
+        or else (X < Big (Int64'(0))) = (Y < Big (Int64'(0))))
+   with Ghost;
+
+   function Round_Quotient (X, Y, Q, R : Big_Integer) return Big_Integer is
+     (if abs R > (abs Y - Big (Int64'(1))) / Big (Int64'(2)) then
+       (if Same_Sign (X, Y) then Q + Big (Int64'(1))
+        else Q - Big (Int64'(1)))
+      else
+        Q)
+   with
+     Ghost,
+     Pre => Y /= 0 and then Q = X / Y and then R = X rem Y;
+
    procedure Scaled_Divide64
      (X, Y, Z : Int64;
       Q, R    : out Int64;
-      Round   : Boolean);
+      Round   : Boolean)
+   with
+     Pre  => Z /= 0
+       and then In_Int64_Range
+         (if Round then Round_Quotient (Big (X) * Big (Y), Big (Z),
+                                        Big (X) * Big (Y) / Big (Z),
+                                        Big (X) * Big (Y) rem Big (Z))
+          else Big (X) * Big (Y) / Big (Z)),
+     Post => Big (R) = Big (X) * Big (Y) rem Big (Z)
+       and then
+         (if Round then
+            Big (Q) = Round_Quotient (Big (X) * Big (Y), Big (Z),
+                                      Big (X) * Big (Y) / Big (Z), Big (R))
+          else
+            Big (Q) = Big (X) * Big (Y) / Big (Z));
    --  Performs the division of (X * Y) / Z, storing the quotient in Q
    --  and the remainder in R. Constraint_Error is raised if Z is zero,
    --  or if the quotient does not fit in 64 bits. Round indicates if
@@ -78,7 +150,22 @@ package System.Arith_64 is
    procedure Double_Divide64
      (X, Y, Z : Int64;
       Q, R    : out Int64;
-      Round   : Boolean);
+      Round   : Boolean)
+   with
+     Pre  => Y /= 0
+       and then Z /= 0
+       and then In_Int64_Range
+         (if Round then Round_Quotient (Big (X), Big (Y) * Big (Z),
+                                        Big (X) / (Big (Y) * Big (Z)),
+                                        Big (X) rem (Big (Y) * Big (Z)))
+          else Big (X) / (Big (Y) * Big (Z))),
+     Post => Big (R) = Big (X) rem (Big (Y) * Big (Z))
+       and then
+         (if Round then
+            Big (Q) = Round_Quotient (Big (X), Big (Y) * Big (Z),
+                                      Big (X) / (Big (Y) * Big (Z)), Big (R))
+          else
+            Big (Q) = Big (X) / (Big (Y) * Big (Z)));
    --  Performs the division X / (Y * Z), storing the quotient in Q and
    --  the remainder in R. Constraint_Error is raised if Y or Z is zero,
    --  or if the quotient does not fit in 64 bits. Round indicates if the
diff --git a/gcc/ada/rtsfind.ads b/gcc/ada/rtsfind.ads
index 99f870a..bedea07 100644
--- a/gcc/ada/rtsfind.ads
+++ b/gcc/ada/rtsfind.ads
@@ -60,6 +60,12 @@ package Rtsfind is
    --  the compilation except in the presence of use clauses, which might
    --  result in unexpected ambiguities.
 
+   --  IMPORTANT NOTE: the specs of packages and procedures with'ed using
+   --  this mechanism must not contain private with clauses. This is because
+   --  the special context installation/removal for private with clauses
+   --  only works in a clean environment for compilation, and could lead
+   --  here to removing visibility over lib units in the calling context.
+
    --  NOTE: If RTU_Id is modified, the subtypes of RTU_Id in the package body
    --  might need to be modified. See Get_Unit_Name.