AES中SBox的生成实现

• AES算法中SBox的生成主要有两步

1.求GF(2^8)有限域内各元素的乘法逆元(扩展欧几里得算法)
2.用第1步的结果做仿射变换(Affine Transformation)

下面结合代码进一步说明

typedef unsigned char uc;
typedef unsigned int ui;
typedef long long int ll;

// 输出: 输入u的仿射变换
uc affineTransformation(uc u);
// 输入:a, b, x, y
// 输出:返回gcd(a, b), 且通过指针返回使得满足 ax + by = gcd(a, b)的x,y
ui exGcd(ui a, ui b, ui* x, ui* y);

因为GF(2^8)域中的元素均可以用8bits表示,所以存储类型我采用char,且由于特殊需要转换成其他类型,不希望内容被符号填充,于是进一步采用unsigned char

那么生成S-Box的整体结构非常清晰：

for (ui i = 0x00; i <= 0xFF; ++i) {
    if(i && i % 16 == 0) printf("\n");

    ui x = 0, y = 0;
    exGcd(0x11b, i, &x, &y);            // 求得元素i的逆y
    uc res = affineTransformation(y);   // 得到逆y的仿射变换
    printf("%02X ", res);
}

细看exGcd(0x11b, i, &x, &y); 原理如下:

ax + by = d = gcd(a, b)  
现在如果gcd(a, b) = 1, 有ax + by = 1, 等式两边同mod a得
[(ax mod a) + (by mod a)]mod a = 1 mod a
0 + by mod a = 1
既然by mod a = 1, 则在GF(2^8)中
y = b^(-1)
只需令a=m(x)=x^8 + x^4 + x^3 + x + 1=0x11b
即可求得元素b的乘法逆元y

我们先看比较简单的仿射变换环节

仿射变换
观察这个矩阵,假设等号右边的三个矩阵分别编号1, 2, 3
留意,这个矩阵2待处理元素的位从上到下是b0..b7,即从低位到高位,我希望能将所有操作都以字节为单位
有些实现在这个地方的处理是把每一位放入bool数组再处理,显得多余。为了能每次操作单位都是一行/一列
比如图中的第一行10001111和b0b1b2b3b4b5b6b7，根据矩阵运算需要相乘(与），然后每位加（异或）起来得一项结果
由于输入u的排列是b7b6...b0和b0b1..b7相反，所以需要对矩阵1的每一行进行一次掉转,比如将10001111掉转成11110001,即0xF1
下面的就很好理解了：

/*
 * 对于某个元素u执行的一次仿射变换
 * 输入: u是某个元素
 * 输出: u的仿射变换
 */
uc affineTransformation(uc u) {
    /*
    * 对书上的mtx 和列矩阵C进行了掉转,从而使得低位对低位,高位对高位
    */
    // 矩阵1: mtx 和矩阵3: c
    const uc mtx[8] = {0xF1, 0xE3, 0xC7, 0x8F, 0x1F, 0x3E, 0x7C, 0xF8}; 
    const uc c = 0x63;  

    /*
     * mul是矩阵1某一行每一位对应 * {b0,...b7} 生成的一列结果
     * curBit从mul低位开始异或,即列结果相加
     * 结果就是矩阵1×矩阵2的中间矩阵tmp的某一位的值 0/1
     * 最后curBit左移i位,和tmp相或,就是置该位值
     */
    uc tmp = 0;
    for (int i = 0; i < 8; ++i) {
        uc mul = mtx[i] & u;
        uc curBit = 0;
        for (int j = 0; j < 8; ++j) {
            curBit ^= (mul & 1);
            mul >>= 1;
        }
        curBit <<= i;
        tmp |= curBit;
    }
    return tmp ^ c;
}

求乘法逆元环节

先分析exGcd算法, 思路如下:

将a和b看作GF(2^8)中的多项式
 ax + by = d = gcd(a, b) 
 显然b = 0时,gcd(a, b) = a,此时x = 1, y = 0
// y = 0 ? 上述时候,其实y为其它也无妨 但是为了
// 兼顾我们在GF(2^8)中自定义的0 的”逆元“取0,才令y=0

下面的目的是探求x2与x1之间的关系
a > b > 0时,设ax1 + by1 = gcd(a, b)
              bx2 + (a % b)y2 = gcd(b, a % b)
即有 ax1 + by1 = bx2 + (a % b)y2
即有 ax1 + by1 = bx2 + (a - [a / b] * b)y2
               = ay2 + bx2 - [a / b] * by2
               = ay2 + b(x2 - [a / b] * y2)
左右相等，故有 x1 = y2, y1 = x2 - [a / b] * y2
于是，得求解x1, y1的方法: x1, y1值基于x2, y2
于是可以递归，gcd不断递归一定会使得某时b = 0,递归可以结束

有了上面x1 y1(递归浅层)与下一层x2 y2(递归深层)的关系,可以写出exGcd算法

E exGcd(E a, E b, E& x, E& y) {
    if (b == 0) {
        x = 1; y = 0;
        return a;
    }
    E q = exGcd(b, a % b, x, y);
    int tmp = x;
    x = y;
    y = tmp - a / b * y;
    return q;
}

可以利用& 引用的性质简化代码

E exGcd(E a, E b, E& x, E& y) {
    if (b == 0) {
        x = 1; y = 0;
        return a;
    }
    E q = exGcd(b, a % b, y, x);
    y = y - a / b * x;
    return q;
}

所以我们可以得知我们的任务:

• 实现多项式的 u % v操作 u / v 操作和 u * v操作
  注意事项
  我们操作的单元一直都是uc,即unsigned char
  知GF(2^8)中,我们要求的是元素 u(x)以m(x)为模的乘法逆元,而m(x) = x^8 + x^4 + x^3 + x + 1
  即0X11B 即 0B100011011, m(x)如果用位来表示,需9bits，故而uc(unsigned char)已经无法胜任,所以我们多引入一种类型来兼顾,即unsigned int(ui)

/*
 * GF(2^8)中的除法, 思路是模拟手工
 * 输出: a / b的商 并且通过residue返回 余式
 * 留意,在求S-Box时,被除式可能是0x11b, 故而uc无法满足
 * 所以这里的被除式,除式都用ui,但是对于余式,不可能大于0xFF,所以用uc即可
 */
ui GFdiv(ui a, ui b, uc* residue) {
    ui res = 0;
    while (a >= b) {
        /*找到a和b的最高位ha hb*/
        ui ha = 8, hb = 8, t = a;
        while (!(t & 0x100)) {--ha; t <<= 1;}
        t = b;
        while (!(t & 0x100)) {--hb; t <<= 1;}

        ui d = ha - hb;
        ui curRes = 0x01 << d;
        res |= curRes;
        ui tmp = polMul(b, curRes);
        a ^= tmp;
    }
    *residue = a;
    return res;
} 

一开始我实现的乘法是GF(2^8)内的乘法
利用的是密码编码学与网络安全第七版 William Stallings中介绍的方法 思路简述如下:
首先基于这样一个事实:
在GF(2^n)中, x^n mod p(x) = [p(x) - x^n]
比如x^8 mod m(x) = [m(x)-x^8] = x^4 + x^3 + x + 1
我们可以利用这个性质对乘法进行拆分完成，比如f(x) * g(x)
g(x) = 11000001 = (10000000) ⊕ (01000000) ⊕（00000001）
预处理：对于f(x)本身,依次乘以1, 10, 100,...10000000,由于每次只左移一位,所以结果最高次只可能=8，此时需mod m(x),相当于异或x^4 + x^3 + x + 1，或者<8，此时无需进行额外操作。得到了这8个预知结果
f(x) * g(x)其实就是根据g(x)把上面预知的结果异或起来就ok.

/*
 * GF(2^8)中的乘法 
 * result = a * b
 */
uc GFmul(uc a, uc b) {
    uc res = 0, tmp;
    if ((b & 1) == 1) res = a;          /*因为0x01时,a不需移位,所以单独处理*/
    b >>= 1;
    for (int i = 0; i < 7; ++i) {       /*剩余处理7位*/
        tmp = a;
        a <<= 1;                        /*a先移位,根据原先的情况条件异或0x1b*/
        if (tmp > 127) a ^= 0x1b;
        if ((b & 1) == 1) res ^= a;
        b >>= 1;
    }
    return res;
}

但是留意这样一种情况：
回到我们所说的除法,模拟除法, m(x) / u(x)
如果u(x) = 1,那么商的第一项会是x^8 即需9位标识商,而上述乘法是GF(2^8)内的，最多8位,所以上述乘法不足以应付m(x) / u(x)，所以我多实现了个一般的多项式乘法供除法调用

/*
 * 一般意义的多项式乘法
 */
ui polMul(ui a, ui b) {
    ui res = 0;
    if ((b & 1) == 1) res = a;
    b >>= 1;
    for (int i = 0; i < 8; ++i) {
        a <<= 1;
        if ((b & 1) == 1) res ^= a;
        b >>= 1;
    }
    return res;
}

至此所有需要的函数实现完毕,最终操作一下,y即b的逆元

ui exGcd(ui a, ui b, ui* x, ui* y) {
    if (b == 0) {
        *x = 1; *y = 0;
        return a;
    }

    uc r;       /*余数r, 商q*/
    ui q = GFdiv(a, b, &r);

    ui res = exGcd(b, r, x, y);
    ui tmp = *x;
    *x = *y;
    *y = tmp ^ GFmul(q, *y);
    return res;
}

可得到S-Box：

63 7C 77 7B F2 6B 6F C5 30 01 67 2B FE D7 AB 76 
CA 82 C9 7D FA 59 47 F0 AD D4 A2 AF 9C A4 72 C0 
B7 FD 93 26 36 3F F7 CC 34 A5 E5 F1 71 D8 31 15 
04 C7 23 C3 18 96 05 9A 07 12 80 E2 EB 27 B2 75 
09 83 2C 1A 1B 6E 5A A0 52 3B D6 B3 29 E3 2F 84 
53 D1 00 ED 20 FC B1 5B 6A CB BE 39 4A 4C 58 CF 
D0 EF AA FB 43 4D 33 85 45 F9 02 7F 50 3C 9F A8 
51 A3 40 8F 92 9D 38 F5 BC B6 DA 21 10 FF F3 D2 
CD 0C 13 EC 5F 97 44 17 C4 A7 7E 3D 64 5D 19 73 
60 81 4F DC 22 2A 90 88 46 EE B8 14 DE 5E 0B DB 
E0 32 3A 0A 49 06 24 5C C2 D3 AC 62 91 95 E4 79 
E7 C8 37 6D 8D D5 4E A9 6C 56 F4 EA 65 7A AE 08 
BA 78 25 2E 1C A6 B4 C6 E8 DD 74 1F 4B BD 8B 8A 
70 3E B5 66 48 03 F6 0E 61 35 57 B9 86 C1 1D 9E 
E1 F8 98 11 69 D9 8E 94 9B 1E 87 E9 CE 55 28 DF 
8C A1 89 0D BF E6 42 68 41 99 2D 0F B0 54 BB 16