Proving properties of constant-time crypto code in SPARKNaCl

Over the last few months, I developed a SPARK version of the TweetNaCl cryptographic library. This was made public on GitHub here in February 2020, under the 2-clause BSD licence. TweetNaCl has achieved widespread popularity and adoption and is used to underpin many other projects. The JavaScript version of the library is used by many web sites, with the NPM package manager reporting over 14 Million downloads per week.

The README file at the top-level in that repository contains most of the background information that you'll need to build and use SPARKNaCl. The code has been proven "type safe" (aka "no run-time errors" or the "Silver" level of SPARK adoption) with the Community 2019 edition of SPARK.

This blog entry goes into a bit more technical detail on one particular aspect of the project: the challenge of re-writing and verifying "constant time" algorithms using SPARK. This is something of a big deal in crypto code. Some algorithms are known to "leak" information because their running time varies with respect to sensitive variables, like cryptographic keys, so a "Constant Time" programming style is adopted to prevent this kind of side-channel attack.

TweetNaCl is something of a masterclass in compact programming - it meets the unusual goal that the entire source code of the library can be transmitted in under 100 tweets. Naturally, this leads to a somewhat economical coding style, and a complete lack of comments in the code. A paper on the authors' website explains some of it, but leaves a lot of detail either unstated or assumed.

SPARKNaCl aims to reverse this situation, proving an alternative reference implementation of the NaCl API that is well-documented, provably free of runtime errors, and preserves the "constant time" nature of all the algorithms.

This blog entry concentrates on just one of the internal functions of TweetNaCl - a function called "Pack_25519".

Let's start with a look at the original C implementation, starting with a little background information on the types and data that we're working with.

TweetNaCl defines a few types and macros that we need to know about - namely:

#define FOR(i,n) for (i = 0; i < n; ++i)
#define sv static void
typedef unsigned char u8;
typedef long long i64;
typedef i64 gf[16];

In particular, the type "gf" is an array of 16 64-bit integers. This is actually used to represent a 256-bit "big integer" where each digit or "limb" is 16 bits, so in the range 0 .. 65535. (64 bits are used so these digits can "overflow" during intermediate computations...) A gf value where all the limbs are really in the range 0 .. 65535 we'll call a "Normal GF".

You also need to know that NaCl implements an Elliptic Curve called "Curve25519" where the "25519" is a reference to the modulus of the curve coordinates, which is actually 2**255 - 19 which (from now on) we will call P. The gf type represents such a value in little-endian format, so the limb at array index 0 represents the least-significant digit.

In the C code, there's a function called "pack25519" with declaration:

sv pack25519(u8 *o, const gf n);

A bit of reading and digging reveals that the parameter o is supposed to be pointing at an uninitialized array of 32 bytes, which is assigned in the function based on the value of the parameter n.

The TweetNaCl paper says that the purpose of Pack25519 is to "freeze integer mod 2**255 − 19 and store".

This gives a bit of a hint. Note that a Normal GF value stores sixteen 16-bit values, which is 256 bits (so representing integers in the range 0 .. 2**256 - 1) but that P is slightly less than half that upper bound. The point of Pack25519, then, is to take a Normal GF value (which might represent an integer greater than or equal to P) and reduce it modulo P, then "pack" the resulting value into 32 8-bit bytes.

This gives us enough to write a first specification for Pack_25519 in SPARK. In SPARK, it is perfectly normal for functions to return a composite type (such as an array) so it seems natural to use a function, not a procedure. We'll also need to introduce a few types as follows. For a longer explanation of the constants and how they are derived, see the SPARKNaCl sources in sparknacl.ads . For the complete source code for Pack_25519, see file sparknacl-utils.adb .

subtype I32 is Integer_32;
subtype N32 is I32 range 0 .. I32'Last;
subtype I64 is Integer_64;

subtype Index_32 is I32 range 0 .. 31;

type Byte_Seq is array (N32 range <>) of Byte;
subtype Bytes_32 is Byte_Seq (Index_32);

--  "LM"   = "Limb Modulus"
--  "LMM1" = "Limb Modulus Minus 1"
LM   : constant := 65536;
LMM1 : constant := 65535;
--  "R2256" = "Remainder of 2**256 (modulo 2**255-19)"
R2256 : constant := 38;

--  "Maximum GF Limb Coefficient"
MGFLC : constant := (R2256 * 15) + 1;

--  "Maximum GF Limb Product"
MGFLP : constant := LMM1 * LMM1;

subtype GF_Any_Limb is I64 range -LM .. (MGFLC * MGFLP);

type GF is array (Index_16) of GF_Any_Limb;

subtype GF_Normal_Limb is I64 range 0 .. LMM1;

subtype Normal_GF is GF
  with Dynamic_Predicate =>
     (for all I in Index_16 => Normal_GF (I) in GF_Normal_Limb);

Note that Normal_GF is a subtype of GF and uses Ada's "Dynamic_Predicate" aspect to constrain the value of each limb.

Finally, we can introduce a SPARK specification for Pack_25519:

--  Reduces N modulo (2**255 - 19) then packs the
--  value into 32 bytes little-endian.
function Pack_25519 (N : in Normal_GF) return Bytes_32
  with Global => null;

sv pack25519 (u8 *o, const gf n)
{
  int i, j, b;
  gf m, t;
  FOR(i,16) t[i]=n[i];
  car25519(t);
  car25519(t);
  car25519(t);
  FOR(j, 2) {
    m[0]=t[0]-0xffed;
    for(i=1;i<15;i++) {
      m[i]=t[i]-0xffff-((m[i-1]>>16)&1);
      m[i-1]&=0xffff;
    }
    m[15]=t[15]-0x7fff-((m[14]>>16)&1);
    b=(m[15]>>16)&1;
    m[14]&=0xffff;
    sel25519 (t, m, 1-b);
  }
  FOR(i, 16) {
    o[2*i]=t[i]&0xff;
    o[2*i+1]=t[i]>>8;
  }
}

sv sel25519 (gf p, gf q, int b);

which doesn't help much. Note that "gf" is an array type, so parameter p and q are by-reference (so probably "in out" in Ada terms). The parameter b is a "int" here but is only ever called with values 0 or 1, so probably being interpreted as a Boolean value. This gives is a reasonable chance of writing a SPARK specification for it. I have changed the name to "CSwap" which seems more readable:

--  Constant time conditional swap of P and Q.
procedure CSwap (P    : in out GF;
                 Q    : in out GF;
                 Swap : in     Boolean)
  with Global => null,
       Contract_Cases =>
         (Swap     => (P = Q'Old and Q = P'Old)
          not Swap => (P = P'Old and Q = Q'Old));

The "Contract_Cases" aspect here says exactly what the procedure does - swapping P and Q (or not) depending on the value of Swap. It might be tempting to write the body trivially:

if Swap then
   Temp := P;
   P := Q;
   Q := Temp;
end if;

but that won't do, since we require the running time to be constant, not even depending on the value of Swap. Let's see how the C code works:

sv sel25519 (gf p, gf q, int b)
{
  i64 t, i, c = ~(b-1);
  FOR(i, 16) {
    t= c&(p[i]^q[i]);
    p[i]^=t;
    q[i]^=t;
  }
}

Once again, this is less than obvious, but becomes clear when you consider the properties of the bitwise "xor" operator (^ in C), and the way the local variable c is initialised.

c is initialised to be "~(b-1)" where b can be only 0 or 1. If b is 0, then b-1 is -1. Bitwise "not" of -1 yields 0. Similarly, if b is 1, then c is initialised to -1 or (more usefully "all 1s" in binary).

This is best handled in SPARK with a constant look-up table:

type Bit_To_Swapmask_Table is array (Boolean) of U64;
Bit_To_Swapmask : constant Bit_To_Swapmask_Table :=
  (False => 16#0000_0000_0000_0000#,
   True  => 16#FFFF_FFFF_FFFF_FFFF#);

Translating and proving the body in SPARK reveals a slightly ugly side to the C version - it uses bitwise operators (like ~ and ^) on signed integers. This is implementation-defined in C since the semantics fall back on whatever representation of the integers happens to be in use. In SPARK, there's no equivalent, so we're forced to use instantiations of Unchecked_Conversion to switch from I64 to U64 and back again. We also need to introduce an axiom for proof:

pragma Assume
  (for all K in I64 => To_I64 (To_U64 (K)) = K);

Writing the body of CSwap in SPARK is then reasonably easy. For proof, we need a loop invariant - this states that after I loop iterations, the first I elements of P and Q have been swapped (or not...) It comes out like this:

procedure CSwap (P    : in out GF;
                 Q    : in out GF;
                 Swap : in     Boolean)
is
   T : U64;
   C : U64 := Bit_To_Swapmask (Swap);
begin
   for I in Index_16 loop
      T := C and (To_U64 (P (I)) xor To_U64 (Q (I)));
      P (I) := To_I64 (To_U64 (P (I)) xor T);
      Q (I) := To_I64 (To_U64 (Q (I)) xor T);

      pragma Loop_Invariant
        (if Swap then
           (for all J in Index_16 range 0 .. I =>
                (P (J) = Q'Loop_Entry (J) and
                 Q (J) = P'Loop_Entry (J)))
         else
           (for all J in Index_16 range 0 .. I =>
                (P (J) = P'Loop_Entry (J) and
                 Q (J) = Q'Loop_Entry (J)))
        );
   end loop;
end CSwap;

--  Subtracting P twice from a Normal_GF might result
--  in a GF where limb 15 can be negative with lower bound -65536
subtype Temp_GF_MSL is I64 range -LM .. LMM1;
subtype Temp_GF is GF
  with Dynamic_Predicate =>
    (Temp_GF (15) in Temp_GF_MSL and
      (for all K in Index_16 range 0 .. 14 =>
         Temp_GF (K) in GF_Normal_Limb));

procedure Subtract_P (T         : in     Temp_GF;
                      Result    :    out Temp_GF;
                      Underflow :    out Boolean)
  with Global => null,
       Pre    => T (15) >= -16#8000#,
       Post   => (Result (15) >= T (15) - 16#8000#);

The body of Subtract_P is basically a simple translation of the C code, but with a suitable loop invariant, which is sufficient to establish the absence of run-time errors. We also need one more subtype, thus:

subtype I64_Bit is I64 range 0 .. 1;

procedure Subtract_P (T         : in     Temp_GF;
                      Result    :    out Temp_GF;
                      Underflow :    out Boolean)
is
   Carry : I64_Bit;
   R     : GF;
begin
   R := (others => 0);

   --  Limb 0 - subtract LSL of P, which is 16#FFED#
   R (0) := T (0) - 16#FFED#;

   --  Limbs 1 .. 14 - subtract FFFF with carry
   for I in Index_16 range 1 .. 14 loop
      Carry     := ASR_16 (R (I - 1)) mod 2;
      R (I)     := T (I) - 16#FFFF# - Carry;
      R (I - 1) := R (I - 1) mod LM;

      pragma Loop_Invariant
        (for all J in Index_16 range 0 .. I - 1 =>
           R (J) in GF_Normal_Limb);
      pragma Loop_Invariant (T in Temp_GF);
   end loop;

   --  Limb 15 - Subtract MSL (Most Significant Limb)
   --  of P (16#7FFF#) with carry.
   --  Note that Limb 15 might become negative on underflow
   Carry  := ASR_16 (R (14)) mod 2;
   R (15) := (T (15) - 16#7FFF#) - Carry;
   R (14) := R (14) mod LM;

   --  Note that R (15) is not normalized here, so that the
   --  result of the first subtraction is numerically correct
   --  as the input to the second.
   Underflow := R (15) < 0;
   Result    := R;
end Subtract_P;

function Pack_25519 (N : in Normal_GF) return Bytes_32
is  
   L      : GF;
   R1, R2 : Temp_GF;
   First_Underflow  : Boolean;
   Second_Underflow : Boolean;
begin
   L := N;
   Subtract_P (L,  R1, First_Underflow);
   Subtract_P (R1, R2, Second_Underflow);
   CSwap (R1, R2, Second_Underflow);
   CSwap (L,  R2, First_Underflow);
   return To_Bytes_32 (R2);
end Pack_25519;

sparknacl-utils.adb:197:27: medium: predicate check might fail

Where line 197 is the return statement with the call to To_Bytes_32.

The problem is that R2 is a Temp_GF (which might have a negative value in the most significant limb), while the formal parameter of To_Bytes_32 is required to be a Normal_GF. Aha! We need to (somehow) keep track of the fact that the value in R2 is always a Normal_GF after the calls to Subtract_P and CSwap.

It suffices to establish that relationship with a post-condition on Subtract_P. In short, we need to know that a subtraction underflows if and only if the result of the Subtraction is not a Normal_GF. This means we extend the specification of Subtract_P as follows:

--  Result := T - P;
--  if     Underflow, then Result is not a Normal_GF
--  if not Underflow, then Result is     a Normal_GF
procedure Subtract_P (T         : in     Temp_GF;
                      Result    :    out Temp_GF;
                      Underflow :    out Boolean)
        with Global => null,
             Pre    => T (15) >= -16#8000#,
             Post   => (Result (15) >= T (15) - 16#8000#) and then
                       (Underflow /= (Result in Normal_GF));

R (14) := R (14) mod LM;

The buggy (pre-2014) version normalised limb 15 instead. To re-introduce this bug into SPARKNaCl, we would change it to:

R (15) := R (15) mod LM;

Running GNATprove on that buggy version yields:

sparknacl-utils.adb:139:23: medium: predicate check might fail