diff --git a/gcc/ada/libgnat/s-expmod.adb b/gcc/ada/libgnat/s-expmod.adb
--- a/gcc/ada/libgnat/s-expmod.adb
+++ b/gcc/ada/libgnat/s-expmod.adb
@@ -29,9 +29,168 @@
 --                                                                          --
 ------------------------------------------------------------------------------
 
-package body System.Exp_Mod is
+--  Preconditions, postconditions, 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,
+                         Ghost          => Ignore,
+                         Loop_Invariant => Ignore,
+                         Assert         => Ignore);
+
+with Ada.Numerics.Big_Numbers.Big_Integers_Ghost;
+use Ada.Numerics.Big_Numbers.Big_Integers_Ghost;
+
+package body System.Exp_Mod
+  with SPARK_Mode
+is
    use System.Unsigned_Types;
 
+   --  Local lemmas
+
+   procedure Lemma_Add_Mod (X, Y : Big_Natural; B : Big_Positive)
+   with
+     Ghost,
+     Post => (X + Y) mod B = ((X mod B) + (Y mod B)) mod B;
+
+   procedure Lemma_Exp_Expand (A : Big_Integer; Exp : Natural)
+   with
+     Ghost,
+     Post =>
+       (if Exp rem 2 = 0 then
+          A ** Exp = A ** (Exp / 2) * A ** (Exp / 2)
+        else
+          A ** Exp = A ** (Exp / 2) * A ** (Exp / 2) * A);
+
+   procedure Lemma_Exp_Mod (A : Big_Natural; Exp : Natural; B : Big_Positive)
+   with
+     Ghost,
+     Subprogram_Variant => (Decreases => Exp),
+     Post => ((A mod B) ** Exp) mod B = (A ** Exp) mod B;
+
+   procedure Lemma_Mod_Ident (A : Big_Natural; B : Big_Positive)
+   with
+     Ghost,
+     Pre => A < B,
+     Post => A mod B = A;
+
+   procedure Lemma_Mod_Mod (A : Big_Integer; B : Big_Positive)
+   with
+     Ghost,
+     Post => A mod B mod B = A mod B;
+
+   procedure Lemma_Mult_Div (X : Big_Natural; Y : Big_Positive)
+   with
+     Ghost,
+     Post => X * Y / Y = X;
+
+   procedure Lemma_Mult_Mod (X, Y : Big_Natural; B : Big_Positive)
+   with
+     Ghost,
+     --  The following subprogram variant can be added as soon as supported
+     --  Subprogram_Variant => (Decreases => Y),
+     Post => (X * Y) mod B = ((X mod B) * (Y mod B)) mod B;
+
+   -----------------------------
+   -- Local lemma null bodies --
+   -----------------------------
+
+   procedure Lemma_Mod_Ident (A : Big_Natural; B : Big_Positive) is null;
+   procedure Lemma_Mod_Mod (A : Big_Integer; B : Big_Positive) is null;
+   procedure Lemma_Mult_Div (X : Big_Natural; Y : Big_Positive) is null;
+
+   -------------------
+   -- Lemma_Add_Mod --
+   -------------------
+
+   procedure Lemma_Add_Mod (X, Y : Big_Natural; B : Big_Positive) is
+      Left  : constant Big_Natural := (X + Y) mod B;
+      Right : constant Big_Natural := ((X mod B) + (Y mod B)) mod B;
+      XQuot : constant Big_Natural := X / B;
+      YQuot : constant Big_Natural := Y / B;
+      AQuot : constant Big_Natural := (X mod B + Y mod B) / B;
+   begin
+      if Y /= 0 and B > 1 then
+         pragma Assert (X = XQuot * B + X mod B);
+         pragma Assert (Y = YQuot * B + Y mod B);
+         pragma Assert
+           (Left = ((XQuot + YQuot) * B + X mod B + Y mod B) mod B);
+         pragma Assert (X mod B + Y mod B = AQuot * B + Right);
+         pragma Assert (Left = ((XQuot + YQuot + AQuot) * B + Right) mod B);
+         pragma Assert (Left = Right);
+      end if;
+   end Lemma_Add_Mod;
+
+   ----------------------
+   -- Lemma_Exp_Expand --
+   ----------------------
+
+   procedure Lemma_Exp_Expand (A : Big_Integer; Exp : Natural) is
+   begin
+      if Exp rem 2 = 0 then
+         pragma Assert (Exp = Exp / 2 + Exp / 2);
+      else
+         pragma Assert (Exp = Exp / 2 + Exp / 2 + 1);
+         pragma Assert (A ** Exp = A ** (Exp / 2) * A ** (Exp / 2 + 1));
+         pragma Assert (A ** (Exp / 2 + 1) = A ** (Exp / 2) * A);
+         pragma Assert (A ** Exp = A ** (Exp / 2) * A ** (Exp / 2) * A);
+      end if;
+   end Lemma_Exp_Expand;
+
+   -------------------
+   -- Lemma_Exp_Mod --
+   -------------------
+
+   procedure Lemma_Exp_Mod (A : Big_Natural; Exp : Natural; B : Big_Positive)
+   is
+   begin
+      if Exp /= 0 then
+         declare
+            Left  : constant Big_Integer := ((A mod B) ** Exp) mod B;
+            Right : constant Big_Integer := (A ** Exp) mod B;
+         begin
+            Lemma_Mult_Mod (A mod B, (A mod B) ** (Exp - 1), B);
+            Lemma_Mod_Mod (A, B);
+            Lemma_Exp_Mod (A, Exp - 1, B);
+            Lemma_Mult_Mod (A, A ** (Exp - 1), B);
+            pragma Assert (Left = Right);
+         end;
+      end if;
+   end Lemma_Exp_Mod;
+
+   --------------------
+   -- Lemma_Mult_Mod --
+   --------------------
+
+   procedure Lemma_Mult_Mod (X, Y : Big_Natural; B : Big_Positive) is
+      Left : constant Big_Natural := (X * Y) mod B;
+      Right : constant Big_Natural := ((X mod B) * (Y mod B)) mod B;
+   begin
+      if Y /= 0 and B > 1 then
+         Lemma_Add_Mod (X * (Y - 1), X, B);
+         Lemma_Mult_Mod (X, Y - 1, B);
+         Lemma_Mod_Mod (X, B);
+         Lemma_Add_Mod ((X mod B) * ((Y - 1) mod B), X mod B, B);
+         Lemma_Add_Mod (Y - 1, 1, B);
+         pragma Assert (((Y - 1) mod B + 1) mod B = Y mod B);
+         if (Y - 1) mod B + 1 < B then
+            Lemma_Mod_Ident ((Y - 1) mod B + 1, B);
+            Lemma_Mod_Mod ((X mod B) * (Y mod B), B);
+            pragma Assert (Left = Right);
+         else
+            pragma Assert (Y mod B = 0);
+            pragma Assert ((X * Y) mod B = (X * Y) - (X * Y) / B * B);
+            pragma Assert
+              ((X * Y) mod B = (X * Y) - (X * (Y / B) * B) / B * B);
+            Lemma_Mult_Div (X * (Y / B), B);
+            pragma Assert (Left = 0);
+            pragma Assert (Left = Right);
+         end if;
+      end if;
+   end Lemma_Mult_Mod;
+
    -----------------
    -- Exp_Modular --
    -----------------
@@ -47,11 +206,36 @@ package body System.Exp_Mod is
 
       function Mult (X, Y : Unsigned) return Unsigned is
         (Unsigned (Long_Long_Unsigned (X) * Long_Long_Unsigned (Y)
-                    mod Long_Long_Unsigned (Modulus)));
+                    mod Long_Long_Unsigned (Modulus)))
+      with
+        Pre => Modulus /= 0;
       --  Modular multiplication. Note that we can't take advantage of the
       --  compiler's circuit, because the modulus is not known statically.
 
+      --  Local ghost variables, functions and lemmas
+
+      M : constant Big_Positive := Big (Modulus) with Ghost;
+
+      function Equal_Modulo (X, Y : Big_Integer) return Boolean is
+         (X mod M = Y mod M)
+      with
+        Ghost,
+        Pre => Modulus /= 0;
+
+      procedure Lemma_Mult (X, Y : Unsigned)
+      with
+        Ghost,
+        Post => Big (Mult (X, Y)) = (Big (X) * Big (Y)) mod M
+          and then Big (Mult (X, Y)) < M;
+
+      procedure Lemma_Mult (X, Y : Unsigned) is null;
+
+      Rest : Big_Integer with Ghost;
+      --  Ghost variable to hold Factor**Exp between Exp and Factor updates
+
    begin
+      pragma Assert (Modulus /= 1);
+
       --  We use the standard logarithmic approach, Exp gets shifted right
       --  testing successive low order bits and Factor is the value of the
       --  base raised to the next power of 2.
@@ -62,14 +246,66 @@ package body System.Exp_Mod is
 
       if Exp /= 0 then
          loop
+            pragma Loop_Invariant (Exp > 0);
+            pragma Loop_Invariant (Result < Modulus);
+            pragma Loop_Invariant (Equal_Modulo
+              (Big (Result) * Big (Factor) ** Exp, Big (Left) ** Right));
+            pragma Loop_Variant (Decreases => Exp);
+            pragma Annotate
+              (CodePeer, False_Positive,
+               "validity check", "confusion on generated code");
+
             if Exp rem 2 /= 0 then
+               pragma Assert
+                 (Big (Factor) ** Exp
+                  = Big (Factor) * Big (Factor) ** (Exp - 1));
+               pragma Assert (Equal_Modulo
+                 ((Big (Result) * Big (Factor)) * Big (Factor) ** (Exp - 1),
+                  Big (Left) ** Right));
+               Lemma_Mult_Mod (Big (Result) * Big (Factor),
+                                  Big (Factor) ** (Exp - 1),
+                                  Big (Modulus));
+               Lemma_Mult (Result, Factor);
+
                Result := Mult (Result, Factor);
+
+               Lemma_Mod_Ident (Big (Result), Big (Modulus));
+               Lemma_Mod_Mod (Big (Factor) ** (Exp - 1), Big (Modulus));
+               Lemma_Mult_Mod (Big (Result),
+                                  Big (Factor) ** (Exp - 1),
+                                  Big (Modulus));
+               pragma Assert (Equal_Modulo
+                 (Big (Result) * Big (Factor) ** (Exp - 1),
+                  Big (Left) ** Right));
+               Lemma_Exp_Expand (Big (Factor), Exp - 1);
             end if;
 
+            Lemma_Exp_Expand (Big (Factor), Exp);
+
             Exp := Exp / 2;
             exit when Exp = 0;
+
+            Rest := Big (Factor) ** Exp;
+            Lemma_Exp_Mod (Big (Factor) * Big (Factor), Exp, Big (Modulus));
+            pragma Assert
+              ((Big (Factor) * Big (Factor)) ** Exp = Rest * Rest);
+            pragma Assert (Equal_Modulo
+              ((Big (Factor) * Big (Factor)) ** Exp,
+               Rest * Rest));
+            Lemma_Mult (Factor, Factor);
+
             Factor := Mult (Factor, Factor);
+
+            Lemma_Mod_Mod (Rest * Rest, Big (Modulus));
+            Lemma_Mod_Ident (Big (Result), Big (Modulus));
+            Lemma_Mult_Mod (Big (Result), Rest * Rest, Big (Modulus));
+            Lemma_Mult_Mod (Big (Result), Big (Factor) ** Exp,
+                               Big (Modulus));
+            pragma Assert (Equal_Modulo
+              (Big (Result) * Big (Factor) ** Exp, Big (Left) ** Right));
          end loop;
+
+         pragma Assert (Big (Result) = Big (Left) ** Right mod Big (Modulus));
       end if;
 
       return Result;


diff --git a/gcc/ada/libgnat/s-expmod.ads b/gcc/ada/libgnat/s-expmod.ads
--- a/gcc/ada/libgnat/s-expmod.ads
+++ b/gcc/ada/libgnat/s-expmod.ads
@@ -35,22 +35,50 @@
 --  Note that 1 is a binary modulus (2**0), so the compiler should not (and
 --  will not) call this function with Modulus equal to 1.
 
+--  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 System.Unsigned_Types;
 
-package System.Exp_Mod is
-   pragma Pure;
+package System.Exp_Mod
+  with Pure, SPARK_Mode
+is
    use type System.Unsigned_Types.Unsigned;
+   subtype Unsigned is System.Unsigned_Types.Unsigned;
+
+   use type Ada.Numerics.Big_Numbers.Big_Integers_Ghost.Big_Integer;
+   subtype Big_Integer is
+     Ada.Numerics.Big_Numbers.Big_Integers_Ghost.Big_Integer
+   with Ghost;
+
+   package Unsigned_Conversion is
+     new Ada.Numerics.Big_Numbers.Big_Integers_Ghost.Unsigned_Conversions
+       (Int => Unsigned);
+
+   function Big (Arg : Unsigned) return Big_Integer is
+     (Unsigned_Conversion.To_Big_Integer (Arg))
+   with Ghost;
 
-   subtype Power_Of_2 is System.Unsigned_Types.Unsigned with
+   subtype Power_Of_2 is Unsigned with
      Dynamic_Predicate =>
         Power_Of_2 /= 0 and then (Power_Of_2 and (Power_Of_2 - 1)) = 0;
 
    function Exp_Modular
-     (Left    : System.Unsigned_Types.Unsigned;
-      Modulus : System.Unsigned_Types.Unsigned;
-      Right   : Natural) return System.Unsigned_Types.Unsigned
+     (Left    : Unsigned;
+      Modulus : Unsigned;
+      Right   : Natural) return Unsigned
    with
-       Pre  => Modulus /= 0 and then Modulus not in Power_Of_2,
-       Post => Exp_Modular'Result = Left ** Right mod Modulus;
+     Pre  => Modulus /= 0 and then Modulus not in Power_Of_2,
+     Post => Big (Exp_Modular'Result) = Big (Left) ** Right mod Big (Modulus);
 
 end System.Exp_Mod;


