C: replace AES FIPS-197 SubBytes table by constant-time code In FIPS-197 (the Advanced Encryption Standard, known as AES), it is made heavy use of SubBytes, which could be implemented as

unsigned char SubBytes(unsigned char x) {
static const unsigned char t[256] = {
  0x63,0x7C,0x77,0x7B,0xF2,0x6B,0x6F,0xC5,0x30,0x01,0x67,0x2B,0xFE,0xD7,0xAB,0x76,
  0xCA,0x82,0xC9,0x7D,0xFA,0x59,0x47,0xF0,0xAD,0xD4,0xA2,0xAF,0x9C,0xA4,0x72,0xC0,
  0xB7,0xFD,0x93,0x26,0x36,0x3F,0xF7,0xCC,0x34,0xA5,0xE5,0xF1,0x71,0xD8,0x31,0x15,
  0x04,0xC7,0x23,0xC3,0x18,0x96,0x05,0x9A,0x07,0x12,0x80,0xE2,0xEB,0x27,0xB2,0x75,
  0x09,0x83,0x2C,0x1A,0x1B,0x6E,0x5A,0xA0,0x52,0x3B,0xD6,0xB3,0x29,0xE3,0x2F,0x84,
  0x53,0xD1,0x00,0xED,0x20,0xFC,0xB1,0x5B,0x6A,0xCB,0xBE,0x39,0x4A,0x4C,0x58,0xCF,
  0xD0,0xEF,0xAA,0xFB,0x43,0x4D,0x33,0x85,0x45,0xF9,0x02,0x7F,0x50,0x3C,0x9F,0xA8,
  0x51,0xA3,0x40,0x8F,0x92,0x9D,0x38,0xF5,0xBC,0xB6,0xDA,0x21,0x10,0xFF,0xF3,0xD2,
  0xCD,0x0C,0x13,0xEC,0x5F,0x97,0x44,0x17,0xC4,0xA7,0x7E,0x3D,0x64,0x5D,0x19,0x73,
  0x60,0x81,0x4F,0xDC,0x22,0x2A,0x90,0x88,0x46,0xEE,0xB8,0x14,0xDE,0x5E,0x0B,0xDB,
  0xE0,0x32,0x3A,0x0A,0x49,0x06,0x24,0x5C,0xC2,0xD3,0xAC,0x62,0x91,0x95,0xE4,0x79,
  0xE7,0xC8,0x37,0x6D,0x8D,0xD5,0x4E,0xA9,0x6C,0x56,0xF4,0xEA,0x65,0x7A,0xAE,0x08,
  0xBA,0x78,0x25,0x2E,0x1C,0xA6,0xB4,0xC6,0xE8,0xDD,0x74,0x1F,0x4B,0xBD,0x8B,0x8A,
  0x70,0x3E,0xB5,0x66,0x48,0x03,0xF6,0x0E,0x61,0x35,0x57,0xB9,0x86,0xC1,0x1D,0x9E,
  0xE1,0xF8,0x98,0x11,0x69,0xD9,0x8E,0x94,0x9B,0x1E,0x87,0xE9,0xCE,0x55,0x28,0xDF,
  0x8C,0xA1,0x89,0x0D,0xBF,0xE6,0x42,0x68,0x41,0x99,0x2D,0x0F,0xB0,0x54,0xBB,0x16};
return t[x];}


This function is not arbitrary; it is a reversible mapping, consisting of an inversion in a Galois Field followed by an affine transformation. All the details are in FIPS-197 section 5.1.1 or here section 4.2.1 (under a slightly different name).

One problem with implementation as a table is that it opens to so-called cache-timing attacks.

Thus your mission is to devise an exact replacement for the above SubBytes() function, that exhibits constant-time behavior; we'll assume that's the case when nothing depending on the input x of SubBytes is used  either:


as an array index, 
as the control operand of if, while, for, case, or operator ?:;
as any operand of operators &amp;&amp;, ||, !, ==, !=, &lt;, &gt;, &lt;=, &gt;=, *, /, %;
as the right operand of operators &gt;&gt;, &lt;&lt;, *=, /=, %=, &lt;&lt;=, &gt;&gt;=.


The winning entry will be one with the lowest cost, obtained from the number of operators executed in the input-dependent data path, with a weight of 5 for unary operators - and ~ as well as for &lt;&lt;1, &gt;&gt;1, +1, -1; weight of 7 for all other operators, shifts with other counts, or adds/subs of other constants (type casts and promotions are free). By principle, that cost is unchanged by unrolling loops (if any), and independent of the input x. As a tie-breaker, the answer with the shortest code after removing whitespace and comments wins.

I plan to designate an entry as answer as early as I can in year 2013, UTC. I will consider answers in languages that I've some knowledge of, ranking them as a straight translation to C not optimized for size.

Apologies for the initial omission of +1 and -1 in the favored operators, of free casts and promotions, and ranking of size. Note that * is prohibited both when unary, and as multiplication.
Score: 940 933 926 910, field tower approach

public class SBox2
{
    public static void main(String[] args)
    {
        for (int i = 0; i &lt; 256; i++) {
            int s = SubBytes(i);
            System.out.format("%02x  ", s);
            if (i % 16 == 15) System.out.println();
        }
    }

    private static int SubBytes(int x) {
        int fwd;
        fwd  = 0x010001 &amp; -(x &amp; 1); x &gt;&gt;= 1; //   7+5+7+5+ | 24+
        fwd ^= 0x1d010f &amp; -(x &amp; 1); x &gt;&gt;= 1; // 7+7+5+7+5+ | 31+
        fwd ^= 0x4f020b &amp; -(x &amp; 1); x &gt;&gt;= 1; // 7+7+5+7+5+ | 31+
        fwd ^= 0x450201 &amp; -(x &amp; 1); x &gt;&gt;= 1; // 7+7+5+7+5+ | 31+
        fwd ^= 0xce080d &amp; -(x &amp; 1); x &gt;&gt;= 1; // 7+7+5+7+5+ | 31+
        fwd ^= 0xa20f0f &amp; -(x &amp; 1); x &gt;&gt;= 1; // 7+7+5+7+5+ | 31+
        fwd ^= 0xc60805 &amp; -(x &amp; 1); x &gt;&gt;= 1; // 7+7+5+7+5+ | 31+
        fwd ^= 0x60070e &amp; -x;                // 7+7+5+     | 19+

        // Running total so far: 229

        int p1;
        {
            int ma = fwd;
            int mb = fwd &gt;&gt; 16;         // 7+         | 7+
            p1  = ma &amp; -(mb&amp;1); ma&lt;&lt;=1; //   7+5+7+5+ | 24+
            p1 ^= ma &amp; -(mb&amp;2); ma&lt;&lt;=1; // 7+7+5+7+5+ | 31+
            p1 ^= ma &amp; -(mb&amp;4); ma&lt;&lt;=1; // 7+7+5+7+5+ | 31+
            p1 ^= ma &amp; -(mb&amp;8);         // 7+7+5+7+   | 26+
            int t = p1 &gt;&gt; 3;            // 7+         | 7+
            p1 ^= (t &gt;&gt; 1) ^ (t &amp; 0xe); // 7+5+7+7+   | 26+
        }

        // Running total so far: 229 + 152 = 381

        int y3, y2, y1, y0;
        {
            int Kinv = (fwd &gt;&gt; 20) ^ p1;     // 7+7+
            int w0 = Kinv &amp; 1; Kinv &gt;&gt;= 1;   // 7+5+
            int w1 = Kinv &amp; 1; Kinv &gt;&gt;= 1;   // 7+5+
            int w2 = Kinv &amp; 1; Kinv &gt;&gt;= 1;   // 7+5+
            int w3 = Kinv &amp; 1;               // 7+

            int t0 = w1 ^ w0 ^ (w2 &amp; w3);      // 7+7+7+
            int t1 = w2 ^ (w0 | w3);           // 7+7+
            int t2 = t0 ^ t1;                  // 7+

            y3 = t2 ^ (t1 &amp; (w1 &amp; w3));        // 7+7+7+
            y2 = t0 ^ (w0 | t2);               // 7+7+
            y1 = w0 ^ w3 ^ (t1 &amp; t0);          // 7+7+7+
            y0 = w3 ^ (t0 | (w1 ^ (w0 | w2))); // 7+7+7+7


        }

        // Running total so far: 381 + 24*7 + 3*5 = 564

        int p2;
        {
            int ma = fwd;
            p2  = ma &amp; -y0; ma&lt;&lt;=1;       //   7+5+5+ | 17+
            p2 ^= ma &amp; -y1; ma&lt;&lt;=1;       // 7+7+5+5+ | 24+
            p2 ^= ma &amp; -y2; ma&lt;&lt;=1;       // 7+7+5+5+ | 24+
            p2 ^= ma &amp; -y3;               // 7+7+5+   | 19+
            int t = p2 &gt;&gt; 3;              // 7+       | 7+
            p2 ^= (t &gt;&gt; 1) ^ (t &amp; 0xe0e); // 7+5+7+7+ | 26
        }

        // Running total so far: 564 + 117 = 681

        int inv8;
        inv8  =  31 &amp; -(p2 &amp; 1);           //   7+5+7+   | 19+
        inv8 ^= 178 &amp; -(p2 &amp; 2); p2 &gt;&gt;= 2; // 7+7+5+7+7+ | 33+
        inv8 ^= 171 &amp; -(p2 &amp; 1);           // 7+7+5+7+   | 26+
        inv8 ^=  54 &amp; -(p2 &amp; 2); p2 &gt;&gt;= 6; // 7+7+5+7+7+ | 33+
        inv8 ^= 188 &amp; -(p2 &amp; 1);           // 7+7+5+7+   | 26+
        inv8 ^=  76 &amp; -(p2 &amp; 2); p2 &gt;&gt;= 2; // 7+7+5+7+7+ | 33+
        inv8 ^= 127 &amp; -(p2 &amp; 1);           // 7+7+5+7+   | 26+
        inv8 ^= 222 &amp; -(p2 &amp; 2);           // 7+7+5+7    | 26+

        return inv8 ^ 0x63;                // 7+         | 7+

        // Grand total: 681 + 229 = 910
    }
}


The structure is essentially the same as Boyar and Peralta's implementation - reduce inversion in GF(2^8) to inversion in GF(2^4), break it down into a linear prologue, a non-linear body, and a linear epilogue, and then minimise them separately. I pay some penalties for bit extraction, but I compensate by being able to do operations in parallel (with some judicious padding of the bits of fwd). In more detail...

Background

As mentioned in the problem description, the S-box consists of an inversion in a particular implementation of the Galois field GF(2^8) followed by an affine transformation. If you know what both of those mean, skip this section.

An affine (or linear) transformation is a function f(x) which respects f(x + y) = f(x) + f(y) and f(a*x) = a*f(x).

A field is a set F of elements with two special elements, which we'll call 0 and 1, and two operators, which we'll call + and *, which respects various properties. In this section, assume that x, y, and z are elements of F.


The elements of F form an abelian group under + with 0 as the identity: i.e. x + y is an element of F; x + 0 = 0 + x = x; each x has a corresponding -x such that x + (-x) = (-x) + x = 0; x + (y + z) = (x + y) + z; and x + y = y + x.
The elements of F other than 0 form an abelian group under * with 1 as the identity.
Multiplication distributes over addition: x * (y + z) = (x * y) + (x * z).


It turns out that there are some quite severe limitations on finite fields:


They must have a number of elements which is a power of a prime.
They are isomorphic with all other finite fields of the same size (i.e. there's really only one finite field of a given size, and any others are just relabellings; call that field GF(p^k) where p is the prime and k is the power).
The multiplicative group F\{0} under * is cyclic; i.e. there's at least one element g such that every element is a power of g.
For powers greater than 1 there is a representation as univariate polynomials of order k of the field of the prime order. E.g. GF(2^8) has a representation in terms of polynomials of x over GF(2). In fact there's usually more than one representation. Consider x^7 * x in GF(2^8); it must be equivalent to some order-7 polynomial, but which? There are a lot of choices which give the right structure; AES chooses to make x^8 = x^4 + x^3 + x + 1 (the lexicographically smallest polynomial which works).


So how do we compute an inverse in that particular representation of GF(2^8)? It's too chunky a problem to tackle directly, so we need to break it down.

Field towers: representing GF(2^8) in terms of GF(2^4)

Instead of representing GF(2^8) with polynomials of 8 terms over GF(2) we can represent it with polynomials of 2 terms over GF(2^4). This time we need to choose a linear polynomial for x^2. Suppose we choose x^2 = px + q. Then (ax + b) * (cx + d) = (ad + bc + acp)x + (bd + acq).

Is it easier to find an inverse in this representation? If (cx + d) = (ax + b)^-1 we get the simultaneous equations


ad + (b + ap)c = 0
bd +     (aq)c = 1


Let D = [b(b+ap) + a^2 q] and set c = a * D^-1; d = (b + ap) * D^-1. So we can do an inverse in GF(2^8) for the cost of a conversion to GF(2^4), an inverse and a few additions and multiplications in GF(2^4), and a conversion back. Even if we do the inverse by means of a table, we've reduced the table size from 256 to 16.

Implementation details

To construct a representation of GF(4) we can choose between three polynomials to reduce x^4:


x^4 = x + 1
x^4 = x^3 + 1
x^4 = x^3 + x^2 + x + 1


The most important difference is in the implementation of multiplication. For any of the three (which correspond to poly of 3, 9, f) the following will work:

// 14x &amp;, 7x unary -, 3x &lt;&lt;1, 3x &gt;&gt;1, 3x &gt;&gt;3, 6x ^ gives score 226
int mul(int a, int b) {
    // Call current values a = a0, b = b0
    int p = a &amp; -(b &amp; 1);
    a = ((a &lt;&lt; 1) ^ (poly &amp; -(a &gt;&gt; 3))) &amp; 15;
    b &gt;&gt;= 1;
    // Now p = a0 * (b0 mod x); a = a0 x; b = b0 div x

    p ^= a &amp; -(b &amp; 1);
    a = ((a &lt;&lt; 1) ^ (poly &amp; -(a &gt;&gt; 3))) &amp; 15;
    b &gt;&gt;= 1;
    // Now p = a0 * (b0 mod x^2); a = a0 x^2; b = b0 div x^2

    p ^= a &amp; -(b &amp; 1);
    a = ((a &lt;&lt; 1) ^ (poly &amp; -(a &gt;&gt; 3))) &amp; 15;
    b &gt;&gt;= 1;
    // Now p = a0 * (b0 mod x^3); a = a0 x^3; b = b0 div x^3

    p ^= a &amp; -(b &amp; 1);
    // p = a0 * b0

    return p;
}


However, if we choose poly = 3 we can handle the overflow much more efficiently because it has a nice structure: there's no double-overflow, because the two inputs are both cubics and x^6 = x^2 (x + 1) is also cubic. In addition we can save the shifts of b: since we're leaving the overflow to last, a0 x^2 doesn't have any set bits corresponding to x or 1 and so we can mask it with -4 instead of -1. The result is

// 10x &amp;, 4x unary -, 3x &lt;&lt;1, 1x &gt;&gt;1, 1x &gt;&gt;3, 5x ^ gives score 152
int mul(int a, int b) {
    int p;
    p  = a &amp; -(b &amp; 1); a &lt;&lt;= 1;
    p ^= a &amp; -(b &amp; 2); a &lt;&lt;= 1;
    p ^= a &amp; -(b &amp; 4); a &lt;&lt;= 1;
    p ^= a &amp; -(b &amp; 8);
    // Here p = a0 * b0 but with overflow, which we need to bring back down.

    int t = p &gt;&gt; 3;
    p ^= (t &gt;&gt; 1) ^ (t &amp; 0xe);
    return p &amp; 15;
}


We still need to choose the values of p and q for the representation of GF(2^8) over GF(2^4), but we don't have many constraints on them. The one thing that matters is that we can get a linear function from the bits of our original representation to the bits of the working representation. This matters for two reasons: firstly, it's easy to do linear transformations, whereas a non-linear transformation would require optimisation equivalent in difficulty to just optimising the entire S-box; secondly, because we can get some side benefits. To recap the structure:

GF256 SubBytes(GF256 x) {
    GF16 a, b, t, D, Dinv, c, d;

    (a, b) = f(x); // f is linear

    t = b + a * p;
    D = b * t + a * a * q;
    Dinv = inverse_GF16(D);
    c = a * Dinv;
    d = t * Dinv;

    return finv(c, d); // finv is also linear
}


If the bits of x are x7 x6 ... x0 then since the transform is linear we get a = f({x7}0000000 + 0{x6}000000 + ... + 0000000{x0}) = f({x7}0000000) + f(0{x6}000000) + ... + f(0000000{x0}). Square it and we get a^2 = f({x7}0000000)^2 + f(0{x6}000000)^2 + ... + f(0000000{x0})^2 where the cross-terms cancel (because in GF(2), 1 + 1 = 0). So a^2 can also be computed as a linear function of x. We can augment the forward linear transform to get:

GF256 SubBytes(GF256 x) {
    GF16 a, b, t, a2q, D, Dinv, c, d;

    (a, b, t, a2q) = faug(x);

    D = b * t + a2q;
    Dinv = inverse_GF16(D);
    c = a * Dinv;
    d = t * Dinv;

    return finv(c, d);
}


and we're down to three multiplications and one addition. We can also extend the multiplication code to do the two multiplications by Dinv in parallel. So our total cost is a forward linear transformation, an addition, two multiplications, an inverse in GF(2^4), and a back linear transformation. We can roll in the post-inverse linear transformation of the S-box and get that essentially for free.

The computation of the coefficients for the linear transformations isn't very interesting, and nor is the micro-optimisation to save a mask here and a shift there. The remaining interesting part is the optimisation of inverse_GF16. There are 2^64 different functions from 4 bits to 4 bits, so a direct optimisation requires lots of memory and time. What I've done is to consider 4 functions from 4 bits to 1 bit, place a cap on the total cost permitted for any one function (with a maximum cost of 63 per function I can enumerate all suitable expressions in under a minute), and for each tuple of functions do common subexpression elimination. After 25 minutes of crunching, I find that the best inverse possible with that cap has a total cost of 133 (an average of 33.25 per bit of the output, which isn't bad considering that the cheapest expression for any individual bit is 35).

I'm still experimenting with other approaches to minimise the inversion in GF(2^4), and an approach which builds bottom-up rather than top-down has yielded an improvement from 133 to 126.