Avoid signed left shift overflows in curve25519-donna: automated conv…

…ersion

Convert all << and <<= to SHL32/64 macros for safe unsigned bitwise
operations in generic curve25519-donna.c (c64 is not affected)
using the perl script below, a simplified version of
scripts/maint/safe_shl.pl

Behaviour and assembly code shouldn't change, as this patch makes explicit
the existing assumption that bits that are shifted out the left are ignored.

Fixes bug 13538.

#!/usr/bin/perl -p -w -i

BEGIN { %vartypes = (); }

if (/^[{}]/) {
    %vartypes = ();
}

# Discover the type of each variable
if (/(?:(?:^ *)|(?:\())(?:const +)?s(\d+) +([a-zA-Z_][_a-zA-Z0-9]*)(\[\d+\])?/) {
    # capture s\d\d declared in function arguments, and/or as arrays
    if (defined $3) {
        $vartypes{"*".$2} = $1;
    } else {
        $vartypes{$2} = $1;
    }
} elsif (/(?:(?:^ *)|(?:\())limb +((?:\*)?[a-zA-Z_][_a-zA-Z0-9]*)/) {
    # capture limb declared in function arguments, and/or as pointers
    $vartypes{$1} = '64';
}

# This fixes at most one shift per line. But that's all the code does.
if (/([a-zA-Z_][a-zA-Z_0-9]*)(?:\[\d+\])? *<< *(\d+)/) {
    $v = $1;
    if (exists $vartypes{$v}) {
        # replace a << n
        s/$v *<< *(\d+)/SHL$vartypes{$v}($v,$1)/;
    } elsif (exists $vartypes{"*".$v}) {
        # replace a[i] << n
        s/$v(\[\d+\]) *<< *(\d+)/SHL$vartypes{"*".$v}($v$1,$2)/;
    }
} elsif (/([a-zA-Z_][a-zA-Z_0-9]*)(?:\[\d+\])? *<<= *(\d+)/) {
    $v = $1;
    if (exists $vartypes{$v}) {
        # replace a <<= n
        s/$v *<<= *(\d+)/$v = SHL$vartypes{$v}($v,$1)/;
    } elsif (exists $vartypes{"*".$v}) {
        # replace a[i] <<= n
        s/$v(\[\d+\]) *<<= *(\d+)/$v$1 = SHL$vartypes{"*".$v}($v$1,$2)/;
    }
}

# Remove extra parenthesis
s/\(SHL64\((.*)\)\)/SHL64\($1\)/;
s/\(SHL32\((.*)\)\)/SHL32\($1\)/;

src/ext/curve25519_donna/curve25519-donna.c

@@ -218,32 +218,32 @@ static void freduce_degree(limb *output) {
    *
    * For output[0..8], the absolute entry value is < 14*2^54 and we add, at
    * most, 19*14*2^54 thus, on exit, |output[0..8]| < 280*2^54. */
-  output[8] += output[18] << 4;
-  output[8] += output[18] << 1;
+  output[8] += SHL64(output[18],4);
+  output[8] += SHL64(output[18],1);
   output[8] += output[18];
-  output[7] += output[17] << 4;
-  output[7] += output[17] << 1;
+  output[7] += SHL64(output[17],4);
+  output[7] += SHL64(output[17],1);
   output[7] += output[17];
-  output[6] += output[16] << 4;
-  output[6] += output[16] << 1;
+  output[6] += SHL64(output[16],4);
+  output[6] += SHL64(output[16],1);
   output[6] += output[16];
-  output[5] += output[15] << 4;
-  output[5] += output[15] << 1;
+  output[5] += SHL64(output[15],4);
+  output[5] += SHL64(output[15],1);
   output[5] += output[15];
-  output[4] += output[14] << 4;
-  output[4] += output[14] << 1;
+  output[4] += SHL64(output[14],4);
+  output[4] += SHL64(output[14],1);
   output[4] += output[14];
-  output[3] += output[13] << 4;
-  output[3] += output[13] << 1;
+  output[3] += SHL64(output[13],4);
+  output[3] += SHL64(output[13],1);
   output[3] += output[13];
-  output[2] += output[12] << 4;
-  output[2] += output[12] << 1;
+  output[2] += SHL64(output[12],4);
+  output[2] += SHL64(output[12],1);
   output[2] += output[12];
-  output[1] += output[11] << 4;
-  output[1] += output[11] << 1;
+  output[1] += SHL64(output[11],4);
+  output[1] += SHL64(output[11],1);
   output[1] += output[11];
-  output[0] += output[10] << 4;
-  output[0] += output[10] << 1;
+  output[0] += SHL64(output[10],4);
+  output[0] += SHL64(output[10],1);
   output[0] += output[10];
 }
 
@@ -309,7 +309,7 @@ static void freduce_coefficients(limb *output) {
      * most, 280*2^28 in the first iteration of this loop. This is added to the
      * next limb and we can approximate the resulting bound of that limb by
      * 281*2^54. */
-    output[i] -= over << 26;
+    output[i] -= SHL64(over,26);
     output[i+1] += over;
 
     /* For the first iteration, |output[i+1]| < 281*2^54, thus |over| <
@@ -319,12 +319,12 @@ static void freduce_coefficients(limb *output) {
      * For subsequent iterations of the loop, 281*2^54 remains a conservative
      * bound and no overflow occurs. */
     over = div_by_2_25(output[i+1]);
-    output[i+1] -= over << 25;
+    output[i+1] -= SHL64(over,25);
     output[i+2] += over;
   }
   /* Now |output[10]| < 281*2^29 and all other coefficients are reduced. */
-  output[0] += output[10] << 4;
-  output[0] += output[10] << 1;
+  output[0] += SHL64(output[10],4);
+  output[0] += SHL64(output[10],1);
   output[0] += output[10];
 
   output[10] = 0;
@@ -333,7 +333,7 @@ static void freduce_coefficients(limb *output) {
    * So |over| will be no more than 2^16. */
   {
     limb over = div_by_2_26(output[0]);
-    output[0] -= over << 26;
+    output[0] -= SHL64(over,26);
     output[1] += over;
   }
 
@@ -470,11 +470,11 @@ fexpand(limb *output, const u8 *input) {
 /* s32_eq returns 0xffffffff iff a == b and zero otherwise. */
 static s32 s32_eq(s32 a, s32 b) {
   a = ~(a ^ b);
-  a &= a << 16;
-  a &= a << 8;
-  a &= a << 4;
-  a &= a << 2;
-  a &= a << 1;
+  a &= SHL32(a,16);
+  a &= SHL32(a,8);
+  a &= SHL32(a,4);
+  a &= SHL32(a,2);
+  a &= SHL32(a,1);
   return a >> 31;
 }
 
@@ -509,12 +509,12 @@ fcontract(u8 *output, limb *input_limbs) {
          * non-negative by borrowing from the next-larger limb. */
         const s32 mask = input[i] >> 31;
         const s32 carry = -((input[i] & mask) >> 25);
-        input[i] = input[i] + (carry << 25);
+        input[i] = input[i] + SHL32(carry,25);
         input[i+1] = input[i+1] - carry;
       } else {
         const s32 mask = input[i] >> 31;
         const s32 carry = -((input[i] & mask) >> 26);
-        input[i] = input[i] + (carry << 26);
+        input[i] = input[i] + SHL32(carry,26);
         input[i+1] = input[i+1] - carry;
       }
     }
@@ -524,7 +524,7 @@ fcontract(u8 *output, limb *input_limbs) {
     {
       const s32 mask = input[9] >> 31;
       const s32 carry = -((input[9] & mask) >> 25);
-      input[9] = input[9] + (carry << 25);
+      input[9] = input[9] + SHL32(carry,25);
       input[0] = input[0] - (carry * 19);
     }
 
@@ -548,7 +548,7 @@ fcontract(u8 *output, limb *input_limbs) {
   {
     const s32 mask = input[0] >> 31;
     const s32 carry = -((input[0] & mask) >> 26);
-    input[0] = input[0] + (carry << 26);
+    input[0] = input[0] + SHL32(carry,26);
     input[1] = input[1] - carry;
   }
 
@@ -605,14 +605,14 @@ fcontract(u8 *output, limb *input_limbs) {
     }
   }
 
-  input[1] <<= 2;
-  input[2] <<= 3;
-  input[3] <<= 5;
-  input[4] <<= 6;
-  input[6] <<= 1;
-  input[7] <<= 3;
-  input[8] <<= 4;
-  input[9] <<= 6;
+  input[1] = SHL32(input[1],2);
+  input[2] = SHL32(input[2],3);
+  input[3] = SHL32(input[3],5);
+  input[4] = SHL32(input[4],6);
+  input[6] = SHL32(input[6],1);
+  input[7] = SHL32(input[7],3);
+  input[8] = SHL32(input[8],4);
+  input[9] = SHL32(input[9],6);
 #define F(i, s) \
   output[s+0] |=  input[i] & 0xff; \
   output[s+1]  = (input[i] >> 8) & 0xff; \