crypto/sha/keccak1600.c

   1 /*
   2  * Copyright 2016 The OpenSSL Project Authors. All Rights Reserved.
   3  *
   4  * Licensed under the OpenSSL license (the "License").  You may not use
   5  * this file except in compliance with the License.  You can obtain a copy
   6  * in the file LICENSE in the source distribution or at
   7  * https://www.openssl.org/source/license.html
   8  */
   9
  10 #include <openssl/e_os2.h>
  11 #include <string.h>
  12 #include <assert.h>
  13
  14 #ifndef KECCAK1600_ASM
  15
  16 #if defined(__x86_64__) || defined(__aarch64__) || \
  17     defined(__mips64) || defined(__ia64) || \
  18     (defined(__VMS) && !defined(__vax))
  19 /*
  20  * These are available even in ILP32 flavours, but even then they are
  21  * capable of performing 64-bit operations as efficiently as in *P64.
  22  * Since it's not given that we can use sizeof(void *), just shunt it.
  23  */
  24 # define BIT_INTERLEAVE (0)
  25 #else
  26 # define BIT_INTERLEAVE (sizeof(void *) < 8)
  27 #endif
  28
  29 #define ROL32(a, offset) (((a) << (offset)) | ((a) >> ((32 - (offset)) & 31)))
  30
  31 static uint64_t ROL64(uint64_t val, int offset)
  32 {
  33     if (offset == 0) {
  34         return val;
  35     } else if (!BIT_INTERLEAVE) {
  36         return (val << offset) | (val >> (64-offset));
  37     } else {
  38         uint32_t hi = (uint32_t)(val >> 32), lo = (uint32_t)val;
  39
  40         if (offset & 1) {
  41             uint32_t tmp = hi;
  42
  43             offset >>= 1;
  44             hi = ROL32(lo, offset);
  45             lo = ROL32(tmp, offset + 1);
  46         } else {
  47             offset >>= 1;
  48             lo = ROL32(lo, offset);
  49             hi = ROL32(hi, offset);
  50         }
  51
  52         return ((uint64_t)hi << 32) | lo;
  53     }
  54 }
  55
  56 static const unsigned char rhotates[5][5] = {
  57     {  0,  1, 62, 28, 27 },
  58     { 36, 44,  6, 55, 20 },
  59     {  3, 10, 43, 25, 39 },
  60     { 41, 45, 15, 21,  8 },
  61     { 18,  2, 61, 56, 14 }
  62 };
  63
  64 static const uint64_t iotas[] = {
  65     BIT_INTERLEAVE ? 0x0000000000000001U : 0x0000000000000001U,
  66     BIT_INTERLEAVE ? 0x0000008900000000U : 0x0000000000008082U,
  67     BIT_INTERLEAVE ? 0x8000008b00000000U : 0x800000000000808aU,
  68     BIT_INTERLEAVE ? 0x8000808000000000U : 0x8000000080008000U,
  69     BIT_INTERLEAVE ? 0x0000008b00000001U : 0x000000000000808bU,
  70     BIT_INTERLEAVE ? 0x0000800000000001U : 0x0000000080000001U,
  71     BIT_INTERLEAVE ? 0x8000808800000001U : 0x8000000080008081U,
  72     BIT_INTERLEAVE ? 0x8000008200000001U : 0x8000000000008009U,
  73     BIT_INTERLEAVE ? 0x0000000b00000000U : 0x000000000000008aU,
  74     BIT_INTERLEAVE ? 0x0000000a00000000U : 0x0000000000000088U,
  75     BIT_INTERLEAVE ? 0x0000808200000001U : 0x0000000080008009U,
  76     BIT_INTERLEAVE ? 0x0000800300000000U : 0x000000008000000aU,
  77     BIT_INTERLEAVE ? 0x0000808b00000001U : 0x000000008000808bU,
  78     BIT_INTERLEAVE ? 0x8000000b00000001U : 0x800000000000008bU,
  79     BIT_INTERLEAVE ? 0x8000008a00000001U : 0x8000000000008089U,
  80     BIT_INTERLEAVE ? 0x8000008100000001U : 0x8000000000008003U,
  81     BIT_INTERLEAVE ? 0x8000008100000000U : 0x8000000000008002U,
  82     BIT_INTERLEAVE ? 0x8000000800000000U : 0x8000000000000080U,
  83     BIT_INTERLEAVE ? 0x0000008300000000U : 0x000000000000800aU,
  84     BIT_INTERLEAVE ? 0x8000800300000000U : 0x800000008000000aU,
  85     BIT_INTERLEAVE ? 0x8000808800000001U : 0x8000000080008081U,
  86     BIT_INTERLEAVE ? 0x8000008800000000U : 0x8000000000008080U,
  87     BIT_INTERLEAVE ? 0x0000800000000001U : 0x0000000080000001U,
  88     BIT_INTERLEAVE ? 0x8000808200000000U : 0x8000000080008008U
  89 };
  90
  91 #if defined(KECCAK_REF)
  92 /*
  93  * This is straightforward or "maximum clarity" implementation aiming
  94  * to resemble section 3.2 of the FIPS PUB 202 "SHA-3 Standard:
  95  * Permutation-Based Hash and Extendible-Output Functions" as much as
  96  * possible. With one caveat. Because of the way C stores matrices,
  97  * references to A[x,y] in the specification are presented as A[y][x].
  98  * Implementation unrolls inner x-loops so that modulo 5 operations are
  99  * explicitly pre-computed.
 100  */
 101 static void Theta(uint64_t A[5][5])
 102 {
 103     uint64_t C[5], D[5];
 104     size_t y;
 105
 106     C[0] = A[0][0];
 107     C[1] = A[0][1];
 108     C[2] = A[0][2];
 109     C[3] = A[0][3];
 110     C[4] = A[0][4];
 111
 112     for (y = 1; y < 5; y++) {
 113         C[0] ^= A[y][0];
 114         C[1] ^= A[y][1];
 115         C[2] ^= A[y][2];
 116         C[3] ^= A[y][3];
 117         C[4] ^= A[y][4];
 118     }
 119
 120     D[0] = ROL64(C[1], 1) ^ C[4];
 121     D[1] = ROL64(C[2], 1) ^ C[0];
 122     D[2] = ROL64(C[3], 1) ^ C[1];
 123     D[3] = ROL64(C[4], 1) ^ C[2];
 124     D[4] = ROL64(C[0], 1) ^ C[3];
 125
 126     for (y = 0; y < 5; y++) {
 127         A[y][0] ^= D[0];
 128         A[y][1] ^= D[1];
 129         A[y][2] ^= D[2];
 130         A[y][3] ^= D[3];
 131         A[y][4] ^= D[4];
 132     }
 133 }
 134
 135 static void Rho(uint64_t A[5][5])
 136 {
 137     size_t y;
 138
 139     for (y = 0; y < 5; y++) {
 140         A[y][0] = ROL64(A[y][0], rhotates[y][0]);
 141         A[y][1] = ROL64(A[y][1], rhotates[y][1]);
 142         A[y][2] = ROL64(A[y][2], rhotates[y][2]);
 143         A[y][3] = ROL64(A[y][3], rhotates[y][3]);
 144         A[y][4] = ROL64(A[y][4], rhotates[y][4]);
 145     }
 146 }
 147
 148 static void Pi(uint64_t A[5][5])
 149 {
 150     uint64_t T[5][5];
 151
 152     /*
 153      * T = A
 154      * A[y][x] = T[x][(3*y+x)%5]
 155      */
 156     memcpy(T, A, sizeof(T));
 157
 158     A[0][0] = T[0][0];
 159     A[0][1] = T[1][1];
 160     A[0][2] = T[2][2];
 161     A[0][3] = T[3][3];
 162     A[0][4] = T[4][4];
 163
 164     A[1][0] = T[0][3];
 165     A[1][1] = T[1][4];
 166     A[1][2] = T[2][0];
 167     A[1][3] = T[3][1];
 168     A[1][4] = T[4][2];
 169
 170     A[2][0] = T[0][1];
 171     A[2][1] = T[1][2];
 172     A[2][2] = T[2][3];
 173     A[2][3] = T[3][4];
 174     A[2][4] = T[4][0];
 175
 176     A[3][0] = T[0][4];
 177     A[3][1] = T[1][0];
 178     A[3][2] = T[2][1];
 179     A[3][3] = T[3][2];
 180     A[3][4] = T[4][3];
 181
 182     A[4][0] = T[0][2];
 183     A[4][1] = T[1][3];
 184     A[4][2] = T[2][4];
 185     A[4][3] = T[3][0];
 186     A[4][4] = T[4][1];
 187 }
 188
 189 static void Chi(uint64_t A[5][5])
 190 {
 191     uint64_t C[5];
 192     size_t y;
 193
 194     for (y = 0; y < 5; y++) {
 195         C[0] = A[y][0] ^ (~A[y][1] & A[y][2]);
 196         C[1] = A[y][1] ^ (~A[y][2] & A[y][3]);
 197         C[2] = A[y][2] ^ (~A[y][3] & A[y][4]);
 198         C[3] = A[y][3] ^ (~A[y][4] & A[y][0]);
 199         C[4] = A[y][4] ^ (~A[y][0] & A[y][1]);
 200
 201         A[y][0] = C[0];
 202         A[y][1] = C[1];
 203         A[y][2] = C[2];
 204         A[y][3] = C[3];
 205         A[y][4] = C[4];
 206     }
 207 }
 208
 209 static void Iota(uint64_t A[5][5], size_t i)
 210 {
 211     assert(i < (sizeof(iotas) / sizeof(iotas[0])));
 212     A[0][0] ^= iotas[i];
 213 }
 214
 215 static void KeccakF1600(uint64_t A[5][5])
 216 {
 217     size_t i;
 218
 219     for (i = 0; i < 24; i++) {
 220         Theta(A);
 221         Rho(A);
 222         Pi(A);
 223         Chi(A);
 224         Iota(A, i);
 225     }
 226 }
 227
 228 #elif defined(KECCAK_1X)
 229 /*
 230  * This implementation is optimization of above code featuring unroll
 231  * of even y-loops, their fusion and code motion. It also minimizes
 232  * temporary storage. Compiler would normally do all these things for
 233  * you, purpose of manual optimization is to provide "unobscured"
 234  * reference for assembly implementation [in case this approach is
 235  * chosen for implementation on some platform]. In the nutshell it's
 236  * equivalent of "plane-per-plane processing" approach discussed in
 237  * section 2.4 of "Keccak implementation overview".
 238  */
 239 static void Round(uint64_t A[5][5], size_t i)
 240 {
 241     uint64_t C[5], E[2];        /* registers */
 242     uint64_t D[5], T[2][5];     /* memory    */
 243
 244     assert(i < (sizeof(iotas) / sizeof(iotas[0])));
 245
 246     C[0] = A[0][0] ^ A[1][0] ^ A[2][0] ^ A[3][0] ^ A[4][0];
 247     C[1] = A[0][1] ^ A[1][1] ^ A[2][1] ^ A[3][1] ^ A[4][1];
 248     C[2] = A[0][2] ^ A[1][2] ^ A[2][2] ^ A[3][2] ^ A[4][2];
 249     C[3] = A[0][3] ^ A[1][3] ^ A[2][3] ^ A[3][3] ^ A[4][3];
 250     C[4] = A[0][4] ^ A[1][4] ^ A[2][4] ^ A[3][4] ^ A[4][4];
 251
 252 #if defined(__arm__)
 253     D[1] = E[0] = ROL64(C[2], 1) ^ C[0];
 254     D[4] = E[1] = ROL64(C[0], 1) ^ C[3];
 255     D[0] = C[0] = ROL64(C[1], 1) ^ C[4];
 256     D[2] = C[1] = ROL64(C[3], 1) ^ C[1];
 257     D[3] = C[2] = ROL64(C[4], 1) ^ C[2];
 258
 259     T[0][0] = A[3][0] ^ C[0]; /* borrow T[0][0] */
 260     T[0][1] = A[0][1] ^ E[0]; /* D[1] */
 261     T[0][2] = A[0][2] ^ C[1]; /* D[2] */
 262     T[0][3] = A[0][3] ^ C[2]; /* D[3] */
 263     T[0][4] = A[0][4] ^ E[1]; /* D[4] */
 264
 265     C[3] = ROL64(A[3][3] ^ C[2], rhotates[3][3]);   /* D[3] */
 266     C[4] = ROL64(A[4][4] ^ E[1], rhotates[4][4]);   /* D[4] */
 267     C[0] =       A[0][0] ^ C[0]; /* rotate by 0 */  /* D[0] */
 268     C[2] = ROL64(A[2][2] ^ C[1], rhotates[2][2]);   /* D[2] */
 269     C[1] = ROL64(A[1][1] ^ E[0], rhotates[1][1]);   /* D[1] */
 270 #else
 271     D[0] = ROL64(C[1], 1) ^ C[4];
 272     D[1] = ROL64(C[2], 1) ^ C[0];
 273     D[2] = ROL64(C[3], 1) ^ C[1];
 274     D[3] = ROL64(C[4], 1) ^ C[2];
 275     D[4] = ROL64(C[0], 1) ^ C[3];
 276
 277     T[0][0] = A[3][0] ^ D[0]; /* borrow T[0][0] */
 278     T[0][1] = A[0][1] ^ D[1];
 279     T[0][2] = A[0][2] ^ D[2];
 280     T[0][3] = A[0][3] ^ D[3];
 281     T[0][4] = A[0][4] ^ D[4];
 282
 283     C[0] =       A[0][0] ^ D[0]; /* rotate by 0 */
 284     C[1] = ROL64(A[1][1] ^ D[1], rhotates[1][1]);
 285     C[2] = ROL64(A[2][2] ^ D[2], rhotates[2][2]);
 286     C[3] = ROL64(A[3][3] ^ D[3], rhotates[3][3]);
 287     C[4] = ROL64(A[4][4] ^ D[4], rhotates[4][4]);
 288 #endif
 289     A[0][0] = C[0] ^ (~C[1] & C[2]) ^ iotas[i];
 290     A[0][1] = C[1] ^ (~C[2] & C[3]);
 291     A[0][2] = C[2] ^ (~C[3] & C[4]);
 292     A[0][3] = C[3] ^ (~C[4] & C[0]);
 293     A[0][4] = C[4] ^ (~C[0] & C[1]);
 294
 295     T[1][0] = A[1][0] ^ (C[3] = D[0]);
 296     T[1][1] = A[2][1] ^ (C[4] = D[1]); /* borrow T[1][1] */
 297     T[1][2] = A[1][2] ^ (E[0] = D[2]);
 298     T[1][3] = A[1][3] ^ (E[1] = D[3]);
 299     T[1][4] = A[2][4] ^ (C[2] = D[4]); /* borrow T[1][4] */
 300
 301     C[0] = ROL64(T[0][3],        rhotates[0][3]);
 302     C[1] = ROL64(A[1][4] ^ C[2], rhotates[1][4]);   /* D[4] */
 303     C[2] = ROL64(A[2][0] ^ C[3], rhotates[2][0]);   /* D[0] */
 304     C[3] = ROL64(A[3][1] ^ C[4], rhotates[3][1]);   /* D[1] */
 305     C[4] = ROL64(A[4][2] ^ E[0], rhotates[4][2]);   /* D[2] */
 306
 307     A[1][0] = C[0] ^ (~C[1] & C[2]);
 308     A[1][1] = C[1] ^ (~C[2] & C[3]);
 309     A[1][2] = C[2] ^ (~C[3] & C[4]);
 310     A[1][3] = C[3] ^ (~C[4] & C[0]);
 311     A[1][4] = C[4] ^ (~C[0] & C[1]);
 312
 313     C[0] = ROL64(T[0][1],        rhotates[0][1]);
 314     C[1] = ROL64(T[1][2],        rhotates[1][2]);
 315     C[2] = ROL64(A[2][3] ^ D[3], rhotates[2][3]);
 316     C[3] = ROL64(A[3][4] ^ D[4], rhotates[3][4]);
 317     C[4] = ROL64(A[4][0] ^ D[0], rhotates[4][0]);
 318
 319     A[2][0] = C[0] ^ (~C[1] & C[2]);
 320     A[2][1] = C[1] ^ (~C[2] & C[3]);
 321     A[2][2] = C[2] ^ (~C[3] & C[4]);
 322     A[2][3] = C[3] ^ (~C[4] & C[0]);
 323     A[2][4] = C[4] ^ (~C[0] & C[1]);
 324
 325     C[0] = ROL64(T[0][4],        rhotates[0][4]);
 326     C[1] = ROL64(T[1][0],        rhotates[1][0]);
 327     C[2] = ROL64(T[1][1],        rhotates[2][1]); /* originally A[2][1] */
 328     C[3] = ROL64(A[3][2] ^ D[2], rhotates[3][2]);
 329     C[4] = ROL64(A[4][3] ^ D[3], rhotates[4][3]);
 330
 331     A[3][0] = C[0] ^ (~C[1] & C[2]);
 332     A[3][1] = C[1] ^ (~C[2] & C[3]);
 333     A[3][2] = C[2] ^ (~C[3] & C[4]);
 334     A[3][3] = C[3] ^ (~C[4] & C[0]);
 335     A[3][4] = C[4] ^ (~C[0] & C[1]);
 336
 337     C[0] = ROL64(T[0][2],        rhotates[0][2]);
 338     C[1] = ROL64(T[1][3],        rhotates[1][3]);
 339     C[2] = ROL64(T[1][4],        rhotates[2][4]); /* originally A[2][4] */
 340     C[3] = ROL64(T[0][0],        rhotates[3][0]); /* originally A[3][0] */
 341     C[4] = ROL64(A[4][1] ^ D[1], rhotates[4][1]);
 342
 343     A[4][0] = C[0] ^ (~C[1] & C[2]);
 344     A[4][1] = C[1] ^ (~C[2] & C[3]);
 345     A[4][2] = C[2] ^ (~C[3] & C[4]);
 346     A[4][3] = C[3] ^ (~C[4] & C[0]);
 347     A[4][4] = C[4] ^ (~C[0] & C[1]);
 348 }
 349
 350 static void KeccakF1600(uint64_t A[5][5])
 351 {
 352     size_t i;
 353
 354     for (i = 0; i < 24; i++) {
 355         Round(A, i);
 356     }
 357 }
 358
 359 #elif defined(KECCAK_1X_ALT)
 360 /*
 361  * This is variant of above KECCAK_1X that reduces requirement for
 362  * temporary storage even further, but at cost of more updates to A[][].
 363  * It's less suitable if A[][] is memory bound, but better if it's
 364  * register bound.
 365  */
 366
 367 static void Round(uint64_t A[5][5], size_t i)
 368 {
 369     uint64_t C[5], D[5];
 370
 371     assert(i < (sizeof(iotas) / sizeof(iotas[0])));
 372
 373     C[0] = A[0][0] ^ A[1][0] ^ A[2][0] ^ A[3][0] ^ A[4][0];
 374     C[1] = A[0][1] ^ A[1][1] ^ A[2][1] ^ A[3][1] ^ A[4][1];
 375     C[2] = A[0][2] ^ A[1][2] ^ A[2][2] ^ A[3][2] ^ A[4][2];
 376     C[3] = A[0][3] ^ A[1][3] ^ A[2][3] ^ A[3][3] ^ A[4][3];
 377     C[4] = A[0][4] ^ A[1][4] ^ A[2][4] ^ A[3][4] ^ A[4][4];
 378
 379     D[1] = C[0] ^  ROL64(C[2], 1);
 380     D[2] = C[1] ^  ROL64(C[3], 1);
 381     D[3] = C[2] ^= ROL64(C[4], 1);
 382     D[4] = C[3] ^= ROL64(C[0], 1);
 383     D[0] = C[4] ^= ROL64(C[1], 1);
 384
 385     A[0][1] ^= D[1];
 386     A[1][1] ^= D[1];
 387     A[2][1] ^= D[1];
 388     A[3][1] ^= D[1];
 389     A[4][1] ^= D[1];
 390
 391     A[0][2] ^= D[2];
 392     A[1][2] ^= D[2];
 393     A[2][2] ^= D[2];
 394     A[3][2] ^= D[2];
 395     A[4][2] ^= D[2];
 396
 397     A[0][3] ^= C[2];
 398     A[1][3] ^= C[2];
 399     A[2][3] ^= C[2];
 400     A[3][3] ^= C[2];
 401     A[4][3] ^= C[2];
 402
 403     A[0][4] ^= C[3];
 404     A[1][4] ^= C[3];
 405     A[2][4] ^= C[3];
 406     A[3][4] ^= C[3];
 407     A[4][4] ^= C[3];
 408
 409     A[0][0] ^= C[4];
 410     A[1][0] ^= C[4];
 411     A[2][0] ^= C[4];
 412     A[3][0] ^= C[4];
 413     A[4][0] ^= C[4];
 414
 415     C[1] = A[0][1];
 416     C[2] = A[0][2];
 417     C[3] = A[0][3];
 418     C[4] = A[0][4];
 419
 420     A[0][1] = ROL64(A[1][1], rhotates[1][1]);
 421     A[0][2] = ROL64(A[2][2], rhotates[2][2]);
 422     A[0][3] = ROL64(A[3][3], rhotates[3][3]);
 423     A[0][4] = ROL64(A[4][4], rhotates[4][4]);
 424
 425     A[1][1] = ROL64(A[1][4], rhotates[1][4]);
 426     A[2][2] = ROL64(A[2][3], rhotates[2][3]);
 427     A[3][3] = ROL64(A[3][2], rhotates[3][2]);
 428     A[4][4] = ROL64(A[4][1], rhotates[4][1]);
 429
 430     A[1][4] = ROL64(A[4][2], rhotates[4][2]);
 431     A[2][3] = ROL64(A[3][4], rhotates[3][4]);
 432     A[3][2] = ROL64(A[2][1], rhotates[2][1]);
 433     A[4][1] = ROL64(A[1][3], rhotates[1][3]);
 434
 435     A[4][2] = ROL64(A[2][4], rhotates[2][4]);
 436     A[3][4] = ROL64(A[4][3], rhotates[4][3]);
 437     A[2][1] = ROL64(A[1][2], rhotates[1][2]);
 438     A[1][3] = ROL64(A[3][1], rhotates[3][1]);
 439
 440     A[2][4] = ROL64(A[4][0], rhotates[4][0]);
 441     A[4][3] = ROL64(A[3][0], rhotates[3][0]);
 442     A[1][2] = ROL64(A[2][0], rhotates[2][0]);
 443     A[3][1] = ROL64(A[1][0], rhotates[1][0]);
 444
 445     A[1][0] = ROL64(C[3],    rhotates[0][3]);
 446     A[2][0] = ROL64(C[1],    rhotates[0][1]);
 447     A[3][0] = ROL64(C[4],    rhotates[0][4]);
 448     A[4][0] = ROL64(C[2],    rhotates[0][2]);
 449
 450     C[0] = A[0][0];
 451     C[1] = A[1][0];
 452     D[0] = A[0][1];
 453     D[1] = A[1][1];
 454
 455     A[0][0] ^= (~A[0][1] & A[0][2]);
 456     A[1][0] ^= (~A[1][1] & A[1][2]);
 457     A[0][1] ^= (~A[0][2] & A[0][3]);
 458     A[1][1] ^= (~A[1][2] & A[1][3]);
 459     A[0][2] ^= (~A[0][3] & A[0][4]);
 460     A[1][2] ^= (~A[1][3] & A[1][4]);
 461     A[0][3] ^= (~A[0][4] & C[0]);
 462     A[1][3] ^= (~A[1][4] & C[1]);
 463     A[0][4] ^= (~C[0]    & D[0]);
 464     A[1][4] ^= (~C[1]    & D[1]);
 465
 466     C[2] = A[2][0];
 467     C[3] = A[3][0];
 468     D[2] = A[2][1];
 469     D[3] = A[3][1];
 470
 471     A[2][0] ^= (~A[2][1] & A[2][2]);
 472     A[3][0] ^= (~A[3][1] & A[3][2]);
 473     A[2][1] ^= (~A[2][2] & A[2][3]);
 474     A[3][1] ^= (~A[3][2] & A[3][3]);
 475     A[2][2] ^= (~A[2][3] & A[2][4]);
 476     A[3][2] ^= (~A[3][3] & A[3][4]);
 477     A[2][3] ^= (~A[2][4] & C[2]);
 478     A[3][3] ^= (~A[3][4] & C[3]);
 479     A[2][4] ^= (~C[2]    & D[2]);
 480     A[3][4] ^= (~C[3]    & D[3]);
 481
 482     C[4] = A[4][0];
 483     D[4] = A[4][1];
 484
 485     A[4][0] ^= (~A[4][1] & A[4][2]);
 486     A[4][1] ^= (~A[4][2] & A[4][3]);
 487     A[4][2] ^= (~A[4][3] & A[4][4]);
 488     A[4][3] ^= (~A[4][4] & C[4]);
 489     A[4][4] ^= (~C[4]    & D[4]);
 490     A[0][0] ^= iotas[i];
 491 }
 492
 493 static void KeccakF1600(uint64_t A[5][5])
 494 {
 495     size_t i;
 496
 497     for (i = 0; i < 24; i++) {
 498         Round(A, i);
 499     }
 500 }
 501
 502 #elif defined(KECCAK_2X)
 503 /*
 504  * This implementation is variant of KECCAK_1X above with outer-most
 505  * round loop unrolled twice. This allows to take temporary storage
 506  * out of round procedure and simplify references to it by alternating
 507  * it with actual data (see round loop below). Just like original, it's
 508  * rather meant as reference for an assembly implementation. It's likely
 509  * to provide best instruction per processed byte ratio at minimal
 510  * round unroll factor...
 511  */
 512 static void Round(uint64_t R[5][5], uint64_t A[5][5], size_t i)
 513 {
 514     uint64_t C[5], D[5];
 515
 516     assert(i < (sizeof(iotas) / sizeof(iotas[0])));
 517
 518     C[0] = A[0][0] ^ A[1][0] ^ A[2][0] ^ A[3][0] ^ A[4][0];
 519     C[1] = A[0][1] ^ A[1][1] ^ A[2][1] ^ A[3][1] ^ A[4][1];
 520     C[2] = A[0][2] ^ A[1][2] ^ A[2][2] ^ A[3][2] ^ A[4][2];
 521     C[3] = A[0][3] ^ A[1][3] ^ A[2][3] ^ A[3][3] ^ A[4][3];
 522     C[4] = A[0][4] ^ A[1][4] ^ A[2][4] ^ A[3][4] ^ A[4][4];
 523
 524     D[0] = ROL64(C[1], 1) ^ C[4];
 525     D[1] = ROL64(C[2], 1) ^ C[0];
 526     D[2] = ROL64(C[3], 1) ^ C[1];
 527     D[3] = ROL64(C[4], 1) ^ C[2];
 528     D[4] = ROL64(C[0], 1) ^ C[3];
 529
 530     C[0] =       A[0][0] ^ D[0]; /* rotate by 0 */
 531     C[1] = ROL64(A[1][1] ^ D[1], rhotates[1][1]);
 532     C[2] = ROL64(A[2][2] ^ D[2], rhotates[2][2]);
 533     C[3] = ROL64(A[3][3] ^ D[3], rhotates[3][3]);
 534     C[4] = ROL64(A[4][4] ^ D[4], rhotates[4][4]);
 535
 536 #ifdef KECCAK_COMPLEMENTING_TRANSFORM
 537     R[0][0] = C[0] ^ ( C[1] | C[2]) ^ iotas[i];
 538     R[0][1] = C[1] ^ (~C[2] | C[3]);
 539     R[0][2] = C[2] ^ ( C[3] & C[4]);
 540     R[0][3] = C[3] ^ ( C[4] | C[0]);
 541     R[0][4] = C[4] ^ ( C[0] & C[1]);
 542 #else
 543     R[0][0] = C[0] ^ (~C[1] & C[2]) ^ iotas[i];
 544     R[0][1] = C[1] ^ (~C[2] & C[3]);
 545     R[0][2] = C[2] ^ (~C[3] & C[4]);
 546     R[0][3] = C[3] ^ (~C[4] & C[0]);
 547     R[0][4] = C[4] ^ (~C[0] & C[1]);
 548 #endif
 549
 550     C[0] = ROL64(A[0][3] ^ D[3], rhotates[0][3]);
 551     C[1] = ROL64(A[1][4] ^ D[4], rhotates[1][4]);
 552     C[2] = ROL64(A[2][0] ^ D[0], rhotates[2][0]);
 553     C[3] = ROL64(A[3][1] ^ D[1], rhotates[3][1]);
 554     C[4] = ROL64(A[4][2] ^ D[2], rhotates[4][2]);
 555
 556 #ifdef KECCAK_COMPLEMENTING_TRANSFORM
 557     R[1][0] = C[0] ^ (C[1] |  C[2]);
 558     R[1][1] = C[1] ^ (C[2] &  C[3]);
 559     R[1][2] = C[2] ^ (C[3] | ~C[4]);
 560     R[1][3] = C[3] ^ (C[4] |  C[0]);
 561     R[1][4] = C[4] ^ (C[0] &  C[1]);
 562 #else
 563     R[1][0] = C[0] ^ (~C[1] & C[2]);
 564     R[1][1] = C[1] ^ (~C[2] & C[3]);
 565     R[1][2] = C[2] ^ (~C[3] & C[4]);
 566     R[1][3] = C[3] ^ (~C[4] & C[0]);
 567     R[1][4] = C[4] ^ (~C[0] & C[1]);
 568 #endif
 569
 570     C[0] = ROL64(A[0][1] ^ D[1], rhotates[0][1]);
 571     C[1] = ROL64(A[1][2] ^ D[2], rhotates[1][2]);
 572     C[2] = ROL64(A[2][3] ^ D[3], rhotates[2][3]);
 573     C[3] = ROL64(A[3][4] ^ D[4], rhotates[3][4]);
 574     C[4] = ROL64(A[4][0] ^ D[0], rhotates[4][0]);
 575
 576 #ifdef KECCAK_COMPLEMENTING_TRANSFORM
 577     R[2][0] =  C[0] ^ ( C[1] | C[2]);
 578     R[2][1] =  C[1] ^ ( C[2] & C[3]);
 579     R[2][2] =  C[2] ^ (~C[3] & C[4]);
 580     R[2][3] = ~C[3] ^ ( C[4] | C[0]);
 581     R[2][4] =  C[4] ^ ( C[0] & C[1]);
 582 #else
 583     R[2][0] = C[0] ^ (~C[1] & C[2]);
 584     R[2][1] = C[1] ^ (~C[2] & C[3]);
 585     R[2][2] = C[2] ^ (~C[3] & C[4]);
 586     R[2][3] = C[3] ^ (~C[4] & C[0]);
 587     R[2][4] = C[4] ^ (~C[0] & C[1]);
 588 #endif
 589
 590     C[0] = ROL64(A[0][4] ^ D[4], rhotates[0][4]);
 591     C[1] = ROL64(A[1][0] ^ D[0], rhotates[1][0]);
 592     C[2] = ROL64(A[2][1] ^ D[1], rhotates[2][1]);
 593     C[3] = ROL64(A[3][2] ^ D[2], rhotates[3][2]);
 594     C[4] = ROL64(A[4][3] ^ D[3], rhotates[4][3]);
 595
 596 #ifdef KECCAK_COMPLEMENTING_TRANSFORM
 597     R[3][0] =  C[0] ^ ( C[1] & C[2]);
 598     R[3][1] =  C[1] ^ ( C[2] | C[3]);
 599     R[3][2] =  C[2] ^ (~C[3] | C[4]);
 600     R[3][3] = ~C[3] ^ ( C[4] & C[0]);
 601     R[3][4] =  C[4] ^ ( C[0] | C[1]);
 602 #else
 603     R[3][0] = C[0] ^ (~C[1] & C[2]);
 604     R[3][1] = C[1] ^ (~C[2] & C[3]);
 605     R[3][2] = C[2] ^ (~C[3] & C[4]);
 606     R[3][3] = C[3] ^ (~C[4] & C[0]);
 607     R[3][4] = C[4] ^ (~C[0] & C[1]);
 608 #endif
 609
 610     C[0] = ROL64(A[0][2] ^ D[2], rhotates[0][2]);
 611     C[1] = ROL64(A[1][3] ^ D[3], rhotates[1][3]);
 612     C[2] = ROL64(A[2][4] ^ D[4], rhotates[2][4]);
 613     C[3] = ROL64(A[3][0] ^ D[0], rhotates[3][0]);
 614     C[4] = ROL64(A[4][1] ^ D[1], rhotates[4][1]);
 615
 616 #ifdef KECCAK_COMPLEMENTING_TRANSFORM
 617     R[4][0] =  C[0] ^ (~C[1] & C[2]);
 618     R[4][1] = ~C[1] ^ ( C[2] | C[3]);
 619     R[4][2] =  C[2] ^ ( C[3] & C[4]);
 620     R[4][3] =  C[3] ^ ( C[4] | C[0]);
 621     R[4][4] =  C[4] ^ ( C[0] & C[1]);
 622 #else
 623     R[4][0] = C[0] ^ (~C[1] & C[2]);
 624     R[4][1] = C[1] ^ (~C[2] & C[3]);
 625     R[4][2] = C[2] ^ (~C[3] & C[4]);
 626     R[4][3] = C[3] ^ (~C[4] & C[0]);
 627     R[4][4] = C[4] ^ (~C[0] & C[1]);
 628 #endif
 629 }
 630
 631 static void KeccakF1600(uint64_t A[5][5])
 632 {
 633     uint64_t T[5][5];
 634     size_t i;
 635
 636 #ifdef KECCAK_COMPLEMENTING_TRANSFORM
 637     A[0][1] = ~A[0][1];
 638     A[0][2] = ~A[0][2];
 639     A[1][3] = ~A[1][3];
 640     A[2][2] = ~A[2][2];
 641     A[3][2] = ~A[3][2];
 642     A[4][0] = ~A[4][0];
 643 #endif
 644
 645     for (i = 0; i < 24; i += 2) {
 646         Round(T, A, i);
 647         Round(A, T, i + 1);
 648     }
 649
 650 #ifdef KECCAK_COMPLEMENTING_TRANSFORM
 651     A[0][1] = ~A[0][1];
 652     A[0][2] = ~A[0][2];
 653     A[1][3] = ~A[1][3];
 654     A[2][2] = ~A[2][2];
 655     A[3][2] = ~A[3][2];
 656     A[4][0] = ~A[4][0];
 657 #endif
 658 }
 659
 660 #else
 661 /*
 662  * This implementation is KECCAK_1X from above combined 4 times with
 663  * a twist that allows to omit temporary storage and perform in-place
 664  * processing. It's discussed in section 2.5 of "Keccak implementation
 665  * overview". It's likely to be best suited for processors with large
 666  * register bank...
 667  */
 668 static void FourRounds(uint64_t A[5][5], size_t i)
 669 {
 670     uint64_t B[5], C[5], D[5];
 671
 672     assert(i <= (sizeof(iotas) / sizeof(iotas[0]) - 4));
 673
 674     /* Round 4*n */
 675     C[0] = A[0][0] ^ A[1][0] ^ A[2][0] ^ A[3][0] ^ A[4][0];
 676     C[1] = A[0][1] ^ A[1][1] ^ A[2][1] ^ A[3][1] ^ A[4][1];
 677     C[2] = A[0][2] ^ A[1][2] ^ A[2][2] ^ A[3][2] ^ A[4][2];
 678     C[3] = A[0][3] ^ A[1][3] ^ A[2][3] ^ A[3][3] ^ A[4][3];
 679     C[4] = A[0][4] ^ A[1][4] ^ A[2][4] ^ A[3][4] ^ A[4][4];
 680
 681     D[0] = ROL64(C[1], 1) ^ C[4];
 682     D[1] = ROL64(C[2], 1) ^ C[0];
 683     D[2] = ROL64(C[3], 1) ^ C[1];
 684     D[3] = ROL64(C[4], 1) ^ C[2];
 685     D[4] = ROL64(C[0], 1) ^ C[3];
 686
 687     B[0] =       A[0][0] ^ D[0]; /* rotate by 0 */
 688     B[1] = ROL64(A[1][1] ^ D[1], rhotates[1][1]);
 689     B[2] = ROL64(A[2][2] ^ D[2], rhotates[2][2]);
 690     B[3] = ROL64(A[3][3] ^ D[3], rhotates[3][3]);
 691     B[4] = ROL64(A[4][4] ^ D[4], rhotates[4][4]);
 692
 693     C[0] = A[0][0] = B[0] ^ (~B[1] & B[2]) ^ iotas[i];
 694     C[1] = A[1][1] = B[1] ^ (~B[2] & B[3]);
 695     C[2] = A[2][2] = B[2] ^ (~B[3] & B[4]);
 696     C[3] = A[3][3] = B[3] ^ (~B[4] & B[0]);
 697     C[4] = A[4][4] = B[4] ^ (~B[0] & B[1]);
 698
 699     B[0] = ROL64(A[0][3] ^ D[3], rhotates[0][3]);
 700     B[1] = ROL64(A[1][4] ^ D[4], rhotates[1][4]);
 701     B[2] = ROL64(A[2][0] ^ D[0], rhotates[2][0]);
 702     B[3] = ROL64(A[3][1] ^ D[1], rhotates[3][1]);
 703     B[4] = ROL64(A[4][2] ^ D[2], rhotates[4][2]);
 704
 705     C[0] ^= A[2][0] = B[0] ^ (~B[1] & B[2]);
 706     C[1] ^= A[3][1] = B[1] ^ (~B[2] & B[3]);
 707     C[2] ^= A[4][2] = B[2] ^ (~B[3] & B[4]);
 708     C[3] ^= A[0][3] = B[3] ^ (~B[4] & B[0]);
 709     C[4] ^= A[1][4] = B[4] ^ (~B[0] & B[1]);
 710
 711     B[0] = ROL64(A[0][1] ^ D[1], rhotates[0][1]);
 712     B[1] = ROL64(A[1][2] ^ D[2], rhotates[1][2]);
 713     B[2] = ROL64(A[2][3] ^ D[3], rhotates[2][3]);
 714     B[3] = ROL64(A[3][4] ^ D[4], rhotates[3][4]);
 715     B[4] = ROL64(A[4][0] ^ D[0], rhotates[4][0]);
 716
 717     C[0] ^= A[4][0] = B[0] ^ (~B[1] & B[2]);
 718     C[1] ^= A[0][1] = B[1] ^ (~B[2] & B[3]);
 719     C[2] ^= A[1][2] = B[2] ^ (~B[3] & B[4]);
 720     C[3] ^= A[2][3] = B[3] ^ (~B[4] & B[0]);
 721     C[4] ^= A[3][4] = B[4] ^ (~B[0] & B[1]);
 722
 723     B[0] = ROL64(A[0][4] ^ D[4], rhotates[0][4]);
 724     B[1] = ROL64(A[1][0] ^ D[0], rhotates[1][0]);
 725     B[2] = ROL64(A[2][1] ^ D[1], rhotates[2][1]);
 726     B[3] = ROL64(A[3][2] ^ D[2], rhotates[3][2]);
 727     B[4] = ROL64(A[4][3] ^ D[3], rhotates[4][3]);
 728
 729     C[0] ^= A[1][0] = B[0] ^ (~B[1] & B[2]);
 730     C[1] ^= A[2][1] = B[1] ^ (~B[2] & B[3]);
 731     C[2] ^= A[3][2] = B[2] ^ (~B[3] & B[4]);
 732     C[3] ^= A[4][3] = B[3] ^ (~B[4] & B[0]);
 733     C[4] ^= A[0][4] = B[4] ^ (~B[0] & B[1]);
 734
 735     B[0] = ROL64(A[0][2] ^ D[2], rhotates[0][2]);
 736     B[1] = ROL64(A[1][3] ^ D[3], rhotates[1][3]);
 737     B[2] = ROL64(A[2][4] ^ D[4], rhotates[2][4]);
 738     B[3] = ROL64(A[3][0] ^ D[0], rhotates[3][0]);
 739     B[4] = ROL64(A[4][1] ^ D[1], rhotates[4][1]);
 740
 741     C[0] ^= A[3][0] = B[0] ^ (~B[1] & B[2]);
 742     C[1] ^= A[4][1] = B[1] ^ (~B[2] & B[3]);
 743     C[2] ^= A[0][2] = B[2] ^ (~B[3] & B[4]);
 744     C[3] ^= A[1][3] = B[3] ^ (~B[4] & B[0]);
 745     C[4] ^= A[2][4] = B[4] ^ (~B[0] & B[1]);
 746
 747     /* Round 4*n+1 */
 748     D[0] = ROL64(C[1], 1) ^ C[4];
 749     D[1] = ROL64(C[2], 1) ^ C[0];
 750     D[2] = ROL64(C[3], 1) ^ C[1];
 751     D[3] = ROL64(C[4], 1) ^ C[2];
 752     D[4] = ROL64(C[0], 1) ^ C[3];
 753
 754     B[0] =       A[0][0] ^ D[0]; /* rotate by 0 */
 755     B[1] = ROL64(A[3][1] ^ D[1], rhotates[1][1]);
 756     B[2] = ROL64(A[1][2] ^ D[2], rhotates[2][2]);
 757     B[3] = ROL64(A[4][3] ^ D[3], rhotates[3][3]);
 758     B[4] = ROL64(A[2][4] ^ D[4], rhotates[4][4]);
 759
 760     C[0] = A[0][0] = B[0] ^ (~B[1] & B[2]) ^ iotas[i + 1];
 761     C[1] = A[3][1] = B[1] ^ (~B[2] & B[3]);
 762     C[2] = A[1][2] = B[2] ^ (~B[3] & B[4]);
 763     C[3] = A[4][3] = B[3] ^ (~B[4] & B[0]);
 764     C[4] = A[2][4] = B[4] ^ (~B[0] & B[1]);
 765
 766     B[0] = ROL64(A[3][3] ^ D[3], rhotates[0][3]);
 767     B[1] = ROL64(A[1][4] ^ D[4], rhotates[1][4]);
 768     B[2] = ROL64(A[4][0] ^ D[0], rhotates[2][0]);
 769     B[3] = ROL64(A[2][1] ^ D[1], rhotates[3][1]);
 770     B[4] = ROL64(A[0][2] ^ D[2], rhotates[4][2]);
 771
 772     C[0] ^= A[4][0] = B[0] ^ (~B[1] & B[2]);
 773     C[1] ^= A[2][1] = B[1] ^ (~B[2] & B[3]);
 774     C[2] ^= A[0][2] = B[2] ^ (~B[3] & B[4]);
 775     C[3] ^= A[3][3] = B[3] ^ (~B[4] & B[0]);
 776     C[4] ^= A[1][4] = B[4] ^ (~B[0] & B[1]);
 777
 778     B[0] = ROL64(A[1][1] ^ D[1], rhotates[0][1]);
 779     B[1] = ROL64(A[4][2] ^ D[2], rhotates[1][2]);
 780     B[2] = ROL64(A[2][3] ^ D[3], rhotates[2][3]);
 781     B[3] = ROL64(A[0][4] ^ D[4], rhotates[3][4]);
 782     B[4] = ROL64(A[3][0] ^ D[0], rhotates[4][0]);
 783
 784     C[0] ^= A[3][0] = B[0] ^ (~B[1] & B[2]);
 785     C[1] ^= A[1][1] = B[1] ^ (~B[2] & B[3]);
 786     C[2] ^= A[4][2] = B[2] ^ (~B[3] & B[4]);
 787     C[3] ^= A[2][3] = B[3] ^ (~B[4] & B[0]);
 788     C[4] ^= A[0][4] = B[4] ^ (~B[0] & B[1]);
 789
 790     B[0] = ROL64(A[4][4] ^ D[4], rhotates[0][4]);
 791     B[1] = ROL64(A[2][0] ^ D[0], rhotates[1][0]);
 792     B[2] = ROL64(A[0][1] ^ D[1], rhotates[2][1]);
 793     B[3] = ROL64(A[3][2] ^ D[2], rhotates[3][2]);
 794     B[4] = ROL64(A[1][3] ^ D[3], rhotates[4][3]);
 795
 796     C[0] ^= A[2][0] = B[0] ^ (~B[1] & B[2]);
 797     C[1] ^= A[0][1] = B[1] ^ (~B[2] & B[3]);
 798     C[2] ^= A[3][2] = B[2] ^ (~B[3] & B[4]);
 799     C[3] ^= A[1][3] = B[3] ^ (~B[4] & B[0]);
 800     C[4] ^= A[4][4] = B[4] ^ (~B[0] & B[1]);
 801
 802     B[0] = ROL64(A[2][2] ^ D[2], rhotates[0][2]);
 803     B[1] = ROL64(A[0][3] ^ D[3], rhotates[1][3]);
 804     B[2] = ROL64(A[3][4] ^ D[4], rhotates[2][4]);
 805     B[3] = ROL64(A[1][0] ^ D[0], rhotates[3][0]);
 806     B[4] = ROL64(A[4][1] ^ D[1], rhotates[4][1]);
 807
 808     C[0] ^= A[1][0] = B[0] ^ (~B[1] & B[2]);
 809     C[1] ^= A[4][1] = B[1] ^ (~B[2] & B[3]);
 810     C[2] ^= A[2][2] = B[2] ^ (~B[3] & B[4]);
 811     C[3] ^= A[0][3] = B[3] ^ (~B[4] & B[0]);
 812     C[4] ^= A[3][4] = B[4] ^ (~B[0] & B[1]);
 813
 814     /* Round 4*n+2 */
 815     D[0] = ROL64(C[1], 1) ^ C[4];
 816     D[1] = ROL64(C[2], 1) ^ C[0];
 817     D[2] = ROL64(C[3], 1) ^ C[1];
 818     D[3] = ROL64(C[4], 1) ^ C[2];
 819     D[4] = ROL64(C[0], 1) ^ C[3];
 820
 821     B[0] =       A[0][0] ^ D[0]; /* rotate by 0 */
 822     B[1] = ROL64(A[2][1] ^ D[1], rhotates[1][1]);
 823     B[2] = ROL64(A[4][2] ^ D[2], rhotates[2][2]);
 824     B[3] = ROL64(A[1][3] ^ D[3], rhotates[3][3]);
 825     B[4] = ROL64(A[3][4] ^ D[4], rhotates[4][4]);
 826
 827     C[0] = A[0][0] = B[0] ^ (~B[1] & B[2]) ^ iotas[i + 2];
 828     C[1] = A[2][1] = B[1] ^ (~B[2] & B[3]);
 829     C[2] = A[4][2] = B[2] ^ (~B[3] & B[4]);
 830     C[3] = A[1][3] = B[3] ^ (~B[4] & B[0]);
 831     C[4] = A[3][4] = B[4] ^ (~B[0] & B[1]);
 832
 833     B[0] = ROL64(A[4][3] ^ D[3], rhotates[0][3]);
 834     B[1] = ROL64(A[1][4] ^ D[4], rhotates[1][4]);
 835     B[2] = ROL64(A[3][0] ^ D[0], rhotates[2][0]);
 836     B[3] = ROL64(A[0][1] ^ D[1], rhotates[3][1]);
 837     B[4] = ROL64(A[2][2] ^ D[2], rhotates[4][2]);
 838
 839     C[0] ^= A[3][0] = B[0] ^ (~B[1] & B[2]);
 840     C[1] ^= A[0][1] = B[1] ^ (~B[2] & B[3]);
 841     C[2] ^= A[2][2] = B[2] ^ (~B[3] & B[4]);
 842     C[3] ^= A[4][3] = B[3] ^ (~B[4] & B[0]);
 843     C[4] ^= A[1][4] = B[4] ^ (~B[0] & B[1]);
 844
 845     B[0] = ROL64(A[3][1] ^ D[1], rhotates[0][1]);
 846     B[1] = ROL64(A[0][2] ^ D[2], rhotates[1][2]);
 847     B[2] = ROL64(A[2][3] ^ D[3], rhotates[2][3]);
 848     B[3] = ROL64(A[4][4] ^ D[4], rhotates[3][4]);
 849     B[4] = ROL64(A[1][0] ^ D[0], rhotates[4][0]);
 850
 851     C[0] ^= A[1][0] = B[0] ^ (~B[1] & B[2]);
 852     C[1] ^= A[3][1] = B[1] ^ (~B[2] & B[3]);
 853     C[2] ^= A[0][2] = B[2] ^ (~B[3] & B[4]);
 854     C[3] ^= A[2][3] = B[3] ^ (~B[4] & B[0]);
 855     C[4] ^= A[4][4] = B[4] ^ (~B[0] & B[1]);
 856
 857     B[0] = ROL64(A[2][4] ^ D[4], rhotates[0][4]);
 858     B[1] = ROL64(A[4][0] ^ D[0], rhotates[1][0]);
 859     B[2] = ROL64(A[1][1] ^ D[1], rhotates[2][1]);
 860     B[3] = ROL64(A[3][2] ^ D[2], rhotates[3][2]);
 861     B[4] = ROL64(A[0][3] ^ D[3], rhotates[4][3]);
 862
 863     C[0] ^= A[4][0] = B[0] ^ (~B[1] & B[2]);
 864     C[1] ^= A[1][1] = B[1] ^ (~B[2] & B[3]);
 865     C[2] ^= A[3][2] = B[2] ^ (~B[3] & B[4]);
 866     C[3] ^= A[0][3] = B[3] ^ (~B[4] & B[0]);
 867     C[4] ^= A[2][4] = B[4] ^ (~B[0] & B[1]);
 868
 869     B[0] = ROL64(A[1][2] ^ D[2], rhotates[0][2]);
 870     B[1] = ROL64(A[3][3] ^ D[3], rhotates[1][3]);
 871     B[2] = ROL64(A[0][4] ^ D[4], rhotates[2][4]);
 872     B[3] = ROL64(A[2][0] ^ D[0], rhotates[3][0]);
 873     B[4] = ROL64(A[4][1] ^ D[1], rhotates[4][1]);
 874
 875     C[0] ^= A[2][0] = B[0] ^ (~B[1] & B[2]);
 876     C[1] ^= A[4][1] = B[1] ^ (~B[2] & B[3]);
 877     C[2] ^= A[1][2] = B[2] ^ (~B[3] & B[4]);
 878     C[3] ^= A[3][3] = B[3] ^ (~B[4] & B[0]);
 879     C[4] ^= A[0][4] = B[4] ^ (~B[0] & B[1]);
 880
 881     /* Round 4*n+3 */
 882     D[0] = ROL64(C[1], 1) ^ C[4];
 883     D[1] = ROL64(C[2], 1) ^ C[0];
 884     D[2] = ROL64(C[3], 1) ^ C[1];
 885     D[3] = ROL64(C[4], 1) ^ C[2];
 886     D[4] = ROL64(C[0], 1) ^ C[3];
 887
 888     B[0] =       A[0][0] ^ D[0]; /* rotate by 0 */
 889     B[1] = ROL64(A[0][1] ^ D[1], rhotates[1][1]);
 890     B[2] = ROL64(A[0][2] ^ D[2], rhotates[2][2]);
 891     B[3] = ROL64(A[0][3] ^ D[3], rhotates[3][3]);
 892     B[4] = ROL64(A[0][4] ^ D[4], rhotates[4][4]);
 893
 894     /* C[0] = */ A[0][0] = B[0] ^ (~B[1] & B[2]) ^ iotas[i + 3];
 895     /* C[1] = */ A[0][1] = B[1] ^ (~B[2] & B[3]);
 896     /* C[2] = */ A[0][2] = B[2] ^ (~B[3] & B[4]);
 897     /* C[3] = */ A[0][3] = B[3] ^ (~B[4] & B[0]);
 898     /* C[4] = */ A[0][4] = B[4] ^ (~B[0] & B[1]);
 899
 900     B[0] = ROL64(A[1][3] ^ D[3], rhotates[0][3]);
 901     B[1] = ROL64(A[1][4] ^ D[4], rhotates[1][4]);
 902     B[2] = ROL64(A[1][0] ^ D[0], rhotates[2][0]);
 903     B[3] = ROL64(A[1][1] ^ D[1], rhotates[3][1]);
 904     B[4] = ROL64(A[1][2] ^ D[2], rhotates[4][2]);
 905
 906     /* C[0] ^= */ A[1][0] = B[0] ^ (~B[1] & B[2]);
 907     /* C[1] ^= */ A[1][1] = B[1] ^ (~B[2] & B[3]);
 908     /* C[2] ^= */ A[1][2] = B[2] ^ (~B[3] & B[4]);
 909     /* C[3] ^= */ A[1][3] = B[3] ^ (~B[4] & B[0]);
 910     /* C[4] ^= */ A[1][4] = B[4] ^ (~B[0] & B[1]);
 911
 912     B[0] = ROL64(A[2][1] ^ D[1], rhotates[0][1]);
 913     B[1] = ROL64(A[2][2] ^ D[2], rhotates[1][2]);
 914     B[2] = ROL64(A[2][3] ^ D[3], rhotates[2][3]);
 915     B[3] = ROL64(A[2][4] ^ D[4], rhotates[3][4]);
 916     B[4] = ROL64(A[2][0] ^ D[0], rhotates[4][0]);
 917
 918     /* C[0] ^= */ A[2][0] = B[0] ^ (~B[1] & B[2]);
 919     /* C[1] ^= */ A[2][1] = B[1] ^ (~B[2] & B[3]);
 920     /* C[2] ^= */ A[2][2] = B[2] ^ (~B[3] & B[4]);
 921     /* C[3] ^= */ A[2][3] = B[3] ^ (~B[4] & B[0]);
 922     /* C[4] ^= */ A[2][4] = B[4] ^ (~B[0] & B[1]);
 923
 924     B[0] = ROL64(A[3][4] ^ D[4], rhotates[0][4]);
 925     B[1] = ROL64(A[3][0] ^ D[0], rhotates[1][0]);
 926     B[2] = ROL64(A[3][1] ^ D[1], rhotates[2][1]);
 927     B[3] = ROL64(A[3][2] ^ D[2], rhotates[3][2]);
 928     B[4] = ROL64(A[3][3] ^ D[3], rhotates[4][3]);
 929
 930     /* C[0] ^= */ A[3][0] = B[0] ^ (~B[1] & B[2]);
 931     /* C[1] ^= */ A[3][1] = B[1] ^ (~B[2] & B[3]);
 932     /* C[2] ^= */ A[3][2] = B[2] ^ (~B[3] & B[4]);
 933     /* C[3] ^= */ A[3][3] = B[3] ^ (~B[4] & B[0]);
 934     /* C[4] ^= */ A[3][4] = B[4] ^ (~B[0] & B[1]);
 935
 936     B[0] = ROL64(A[4][2] ^ D[2], rhotates[0][2]);
 937     B[1] = ROL64(A[4][3] ^ D[3], rhotates[1][3]);
 938     B[2] = ROL64(A[4][4] ^ D[4], rhotates[2][4]);
 939     B[3] = ROL64(A[4][0] ^ D[0], rhotates[3][0]);
 940     B[4] = ROL64(A[4][1] ^ D[1], rhotates[4][1]);
 941
 942     /* C[0] ^= */ A[4][0] = B[0] ^ (~B[1] & B[2]);
 943     /* C[1] ^= */ A[4][1] = B[1] ^ (~B[2] & B[3]);
 944     /* C[2] ^= */ A[4][2] = B[2] ^ (~B[3] & B[4]);
 945     /* C[3] ^= */ A[4][3] = B[3] ^ (~B[4] & B[0]);
 946     /* C[4] ^= */ A[4][4] = B[4] ^ (~B[0] & B[1]);
 947 }
 948
 949 static void KeccakF1600(uint64_t A[5][5])
 950 {
 951     size_t i;
 952
 953     for (i = 0; i < 24; i += 4) {
 954         FourRounds(A, i);
 955     }
 956 }
 957
 958 #endif
 959
 960 static uint64_t BitInterleave(uint64_t Ai)
 961 {
 962     if (BIT_INTERLEAVE) {
 963         uint32_t hi = (uint32_t)(Ai >> 32), lo = (uint32_t)Ai;
 964         uint32_t t0, t1;
 965
 966         t0 = lo & 0x55555555;
 967         t0 |= t0 >> 1;  t0 &= 0x33333333;
 968         t0 |= t0 >> 2;  t0 &= 0x0f0f0f0f;
 969         t0 |= t0 >> 4;  t0 &= 0x00ff00ff;
 970         t0 |= t0 >> 8;  t0 &= 0x0000ffff;
 971
 972         t1 = hi & 0x55555555;
 973         t1 |= t1 >> 1;  t1 &= 0x33333333;
 974         t1 |= t1 >> 2;  t1 &= 0x0f0f0f0f;
 975         t1 |= t1 >> 4;  t1 &= 0x00ff00ff;
 976         t1 |= t1 >> 8;  t1 <<= 16;
 977
 978         lo &= 0xaaaaaaaa;
 979         lo |= lo << 1;  lo &= 0xcccccccc;
 980         lo |= lo << 2;  lo &= 0xf0f0f0f0;
 981         lo |= lo << 4;  lo &= 0xff00ff00;
 982         lo |= lo << 8;  lo >>= 16;
 983
 984         hi &= 0xaaaaaaaa;
 985         hi |= hi << 1;  hi &= 0xcccccccc;
 986         hi |= hi << 2;  hi &= 0xf0f0f0f0;
 987         hi |= hi << 4;  hi &= 0xff00ff00;
 988         hi |= hi << 8;  hi &= 0xffff0000;
 989
 990         Ai = ((uint64_t)(hi | lo) << 32) | (t1 | t0);
 991     }
 992
 993     return Ai;
 994 }
 995
 996 static uint64_t BitDeinterleave(uint64_t Ai)
 997 {
 998     if (BIT_INTERLEAVE) {
 999         uint32_t hi = (uint32_t)(Ai >> 32), lo = (uint32_t)Ai;
1000         uint32_t t0, t1;
1001
1002         t0 = lo & 0x0000ffff;
1003         t0 |= t0 << 8;  t0 &= 0x00ff00ff;
1004         t0 |= t0 << 4;  t0 &= 0x0f0f0f0f;
1005         t0 |= t0 << 2;  t0 &= 0x33333333;
1006         t0 |= t0 << 1;  t0 &= 0x55555555;
1007
1008         t1 = hi << 16;
1009         t1 |= t1 >> 8;  t1 &= 0xff00ff00;
1010         t1 |= t1 >> 4;  t1 &= 0xf0f0f0f0;
1011         t1 |= t1 >> 2;  t1 &= 0xcccccccc;
1012         t1 |= t1 >> 1;  t1 &= 0xaaaaaaaa;
1013
1014         lo >>= 16;
1015         lo |= lo << 8;  lo &= 0x00ff00ff;
1016         lo |= lo << 4;  lo &= 0x0f0f0f0f;
1017         lo |= lo << 2;  lo &= 0x33333333;
1018         lo |= lo << 1;  lo &= 0x55555555;
1019
1020         hi &= 0xffff0000;
1021         hi |= hi >> 8;  hi &= 0xff00ff00;
1022         hi |= hi >> 4;  hi &= 0xf0f0f0f0;
1023         hi |= hi >> 2;  hi &= 0xcccccccc;
1024         hi |= hi >> 1;  hi &= 0xaaaaaaaa;
1025
1026         Ai = ((uint64_t)(hi | lo) << 32) | (t1 | t0);
1027     }
1028
1029     return Ai;
1030 }
1031
1032 /*
1033  * SHA3_absorb can be called multiple times, but at each invocation
1034  * largest multiple of |r| out of |len| bytes are processed. Then
1035  * remaining amount of bytes is returned. This is done to spare caller
1036  * trouble of calculating the largest multiple of |r|. |r| can be viewed
1037  * as blocksize. It is commonly (1600 - 256*n)/8, e.g. 168, 136, 104,
1038  * 72, but can also be (1600 - 448)/8 = 144. All this means that message
1039  * padding and intermediate sub-block buffering, byte- or bitwise, is
1040  * caller's reponsibility.
1041  */
1042 size_t SHA3_absorb(uint64_t A[5][5], const unsigned char *inp, size_t len,
1043                    size_t r)
1044 {
1045     uint64_t *A_flat = (uint64_t *)A;
1046     size_t i, w = r / 8;
1047
1048     assert(r < (25 * sizeof(A[0][0])) && (r % 8) == 0);
1049
1050     while (len >= r) {
1051         for (i = 0; i < w; i++) {
1052             uint64_t Ai = (uint64_t)inp[0]       | (uint64_t)inp[1] << 8  |
1053                           (uint64_t)inp[2] << 16 | (uint64_t)inp[3] << 24 |
1054                           (uint64_t)inp[4] << 32 | (uint64_t)inp[5] << 40 |
1055                           (uint64_t)inp[6] << 48 | (uint64_t)inp[7] << 56;
1056             inp += 8;
1057
1058             A_flat[i] ^= BitInterleave(Ai);
1059         }
1060         KeccakF1600(A);
1061         len -= r;
1062     }
1063
1064     return len;
1065 }
1066
1067 /*
1068  * SHA3_squeeze is called once at the end to generate |out| hash value
1069  * of |len| bytes.
1070  */
1071 void SHA3_squeeze(uint64_t A[5][5], unsigned char *out, size_t len, size_t r)
1072 {
1073     uint64_t *A_flat = (uint64_t *)A;
1074     size_t i, w = r / 8;
1075
1076     assert(r < (25 * sizeof(A[0][0])) && (r % 8) == 0);
1077
1078     while (len != 0) {
1079         for (i = 0; i < w && len != 0; i++) {
1080             uint64_t Ai = BitDeinterleave(A_flat[i]);
1081
1082             if (len < 8) {
1083                 for (i = 0; i < len; i++) {
1084                     *out++ = (unsigned char)Ai;
1085                     Ai >>= 8;
1086                 }
1087                 return;
1088             }
1089
1090             out[0] = (unsigned char)(Ai);
1091             out[1] = (unsigned char)(Ai >> 8);
1092             out[2] = (unsigned char)(Ai >> 16);
1093             out[3] = (unsigned char)(Ai >> 24);
1094             out[4] = (unsigned char)(Ai >> 32);
1095             out[5] = (unsigned char)(Ai >> 40);
1096             out[6] = (unsigned char)(Ai >> 48);
1097             out[7] = (unsigned char)(Ai >> 56);
1098             out += 8;
1099             len -= 8;
1100         }
1101         if (len)
1102             KeccakF1600(A);
1103     }
1104 }
1105 #else
1106 size_t SHA3_absorb(uint64_t A[5][5], const unsigned char *inp, size_t len,
1107                    size_t r);
1108 void SHA3_squeeze(uint64_t A[5][5], unsigned char *out, size_t len, size_t r);
1109 #endif
1110
1111 #ifdef SELFTEST
1112 /*
1113  * Post-padding one-shot implementations would look as following:
1114  *
1115  * SHA3_224     SHA3_sponge(inp, len, out, 224/8, (1600-448)/8);
1116  * SHA3_256     SHA3_sponge(inp, len, out, 256/8, (1600-512)/8);
1117  * SHA3_384     SHA3_sponge(inp, len, out, 384/8, (1600-768)/8);
1118  * SHA3_512     SHA3_sponge(inp, len, out, 512/8, (1600-1024)/8);
1119  * SHAKE_128    SHA3_sponge(inp, len, out, d, (1600-256)/8);
1120  * SHAKE_256    SHA3_sponge(inp, len, out, d, (1600-512)/8);
1121  */
1122
1123 void SHA3_sponge(const unsigned char *inp, size_t len,
1124                  unsigned char *out, size_t d, size_t r)
1125 {
1126     uint64_t A[5][5];
1127
1128     memset(A, 0, sizeof(A));
1129     SHA3_absorb(A, inp, len, r);
1130     SHA3_squeeze(A, out, d, r);
1131 }
1132
1133 # include <stdio.h>
1134
1135 int main()
1136 {
1137     /*
1138      * This is 5-bit SHAKE128 test from http://csrc.nist.gov/groups/ST/toolkit/examples.html#aHashing
1139      */
1140     unsigned char test[168] = { '\xf3', '\x3' };
1141     unsigned char out[512];
1142     size_t i;
1143     static const unsigned char result[512] = {
1144         0x2E, 0x0A, 0xBF, 0xBA, 0x83, 0xE6, 0x72, 0x0B,
1145         0xFB, 0xC2, 0x25, 0xFF, 0x6B, 0x7A, 0xB9, 0xFF,
1146         0xCE, 0x58, 0xBA, 0x02, 0x7E, 0xE3, 0xD8, 0x98,
1147         0x76, 0x4F, 0xEF, 0x28, 0x7D, 0xDE, 0xCC, 0xCA,
1148         0x3E, 0x6E, 0x59, 0x98, 0x41, 0x1E, 0x7D, 0xDB,
1149         0x32, 0xF6, 0x75, 0x38, 0xF5, 0x00, 0xB1, 0x8C,
1150         0x8C, 0x97, 0xC4, 0x52, 0xC3, 0x70, 0xEA, 0x2C,
1151         0xF0, 0xAF, 0xCA, 0x3E, 0x05, 0xDE, 0x7E, 0x4D,
1152         0xE2, 0x7F, 0xA4, 0x41, 0xA9, 0xCB, 0x34, 0xFD,
1153         0x17, 0xC9, 0x78, 0xB4, 0x2D, 0x5B, 0x7E, 0x7F,
1154         0x9A, 0xB1, 0x8F, 0xFE, 0xFF, 0xC3, 0xC5, 0xAC,
1155         0x2F, 0x3A, 0x45, 0x5E, 0xEB, 0xFD, 0xC7, 0x6C,
1156         0xEA, 0xEB, 0x0A, 0x2C, 0xCA, 0x22, 0xEE, 0xF6,
1157         0xE6, 0x37, 0xF4, 0xCA, 0xBE, 0x5C, 0x51, 0xDE,
1158         0xD2, 0xE3, 0xFA, 0xD8, 0xB9, 0x52, 0x70, 0xA3,
1159         0x21, 0x84, 0x56, 0x64, 0xF1, 0x07, 0xD1, 0x64,
1160         0x96, 0xBB, 0x7A, 0xBF, 0xBE, 0x75, 0x04, 0xB6,
1161         0xED, 0xE2, 0xE8, 0x9E, 0x4B, 0x99, 0x6F, 0xB5,
1162         0x8E, 0xFD, 0xC4, 0x18, 0x1F, 0x91, 0x63, 0x38,
1163         0x1C, 0xBE, 0x7B, 0xC0, 0x06, 0xA7, 0xA2, 0x05,
1164         0x98, 0x9C, 0x52, 0x6C, 0xD1, 0xBD, 0x68, 0x98,
1165         0x36, 0x93, 0xB4, 0xBD, 0xC5, 0x37, 0x28, 0xB2,
1166         0x41, 0xC1, 0xCF, 0xF4, 0x2B, 0xB6, 0x11, 0x50,
1167         0x2C, 0x35, 0x20, 0x5C, 0xAB, 0xB2, 0x88, 0x75,
1168         0x56, 0x55, 0xD6, 0x20, 0xC6, 0x79, 0x94, 0xF0,
1169         0x64, 0x51, 0x18, 0x7F, 0x6F, 0xD1, 0x7E, 0x04,
1170         0x66, 0x82, 0xBA, 0x12, 0x86, 0x06, 0x3F, 0xF8,
1171         0x8F, 0xE2, 0x50, 0x8D, 0x1F, 0xCA, 0xF9, 0x03,
1172         0x5A, 0x12, 0x31, 0xAD, 0x41, 0x50, 0xA9, 0xC9,
1173         0xB2, 0x4C, 0x9B, 0x2D, 0x66, 0xB2, 0xAD, 0x1B,
1174         0xDE, 0x0B, 0xD0, 0xBB, 0xCB, 0x8B, 0xE0, 0x5B,
1175         0x83, 0x52, 0x29, 0xEF, 0x79, 0x19, 0x73, 0x73,
1176         0x23, 0x42, 0x44, 0x01, 0xE1, 0xD8, 0x37, 0xB6,
1177         0x6E, 0xB4, 0xE6, 0x30, 0xFF, 0x1D, 0xE7, 0x0C,
1178         0xB3, 0x17, 0xC2, 0xBA, 0xCB, 0x08, 0x00, 0x1D,
1179         0x34, 0x77, 0xB7, 0xA7, 0x0A, 0x57, 0x6D, 0x20,
1180         0x86, 0x90, 0x33, 0x58, 0x9D, 0x85, 0xA0, 0x1D,
1181         0xDB, 0x2B, 0x66, 0x46, 0xC0, 0x43, 0xB5, 0x9F,
1182         0xC0, 0x11, 0x31, 0x1D, 0xA6, 0x66, 0xFA, 0x5A,
1183         0xD1, 0xD6, 0x38, 0x7F, 0xA9, 0xBC, 0x40, 0x15,
1184         0xA3, 0x8A, 0x51, 0xD1, 0xDA, 0x1E, 0xA6, 0x1D,
1185         0x64, 0x8D, 0xC8, 0xE3, 0x9A, 0x88, 0xB9, 0xD6,
1186         0x22, 0xBD, 0xE2, 0x07, 0xFD, 0xAB, 0xC6, 0xF2,
1187         0x82, 0x7A, 0x88, 0x0C, 0x33, 0x0B, 0xBF, 0x6D,
1188         0xF7, 0x33, 0x77, 0x4B, 0x65, 0x3E, 0x57, 0x30,
1189         0x5D, 0x78, 0xDC, 0xE1, 0x12, 0xF1, 0x0A, 0x2C,
1190         0x71, 0xF4, 0xCD, 0xAD, 0x92, 0xED, 0x11, 0x3E,
1191         0x1C, 0xEA, 0x63, 0xB9, 0x19, 0x25, 0xED, 0x28,
1192         0x19, 0x1E, 0x6D, 0xBB, 0xB5, 0xAA, 0x5A, 0x2A,
1193         0xFD, 0xA5, 0x1F, 0xC0, 0x5A, 0x3A, 0xF5, 0x25,
1194         0x8B, 0x87, 0x66, 0x52, 0x43, 0x55, 0x0F, 0x28,
1195         0x94, 0x8A, 0xE2, 0xB8, 0xBE, 0xB6, 0xBC, 0x9C,
1196         0x77, 0x0B, 0x35, 0xF0, 0x67, 0xEA, 0xA6, 0x41,
1197         0xEF, 0xE6, 0x5B, 0x1A, 0x44, 0x90, 0x9D, 0x1B,
1198         0x14, 0x9F, 0x97, 0xEE, 0xA6, 0x01, 0x39, 0x1C,
1199         0x60, 0x9E, 0xC8, 0x1D, 0x19, 0x30, 0xF5, 0x7C,
1200         0x18, 0xA4, 0xE0, 0xFA, 0xB4, 0x91, 0xD1, 0xCA,
1201         0xDF, 0xD5, 0x04, 0x83, 0x44, 0x9E, 0xDC, 0x0F,
1202         0x07, 0xFF, 0xB2, 0x4D, 0x2C, 0x6F, 0x9A, 0x9A,
1203         0x3B, 0xFF, 0x39, 0xAE, 0x3D, 0x57, 0xF5, 0x60,
1204         0x65, 0x4D, 0x7D, 0x75, 0xC9, 0x08, 0xAB, 0xE6,
1205         0x25, 0x64, 0x75, 0x3E, 0xAC, 0x39, 0xD7, 0x50,
1206         0x3D, 0xA6, 0xD3, 0x7C, 0x2E, 0x32, 0xE1, 0xAF,
1207         0x3B, 0x8A, 0xEC, 0x8A, 0xE3, 0x06, 0x9C, 0xD9
1208     };
1209
1210     test[167] = '\x80';
1211     SHA3_sponge(test, sizeof(test), out, sizeof(out), sizeof(test));
1212
1213     /*
1214      * Rationale behind keeping output [formatted as below] is that
1215      * one should be able to redirect it to a file, then copy-n-paste
1216      * final "output val" from official example to another file, and
1217      * compare the two with diff(1).
1218      */
1219     for (i = 0; i < sizeof(out);) {
1220         printf("%02X", out[i]);
1221         printf(++i % 16 && i != sizeof(out) ? " " : "\n");
1222     }
1223
1224     if (memcmp(out,result,sizeof(out))) {
1225         fprintf(stderr,"failure\n");
1226         return 1;
1227     } else {
1228         fprintf(stderr,"success\n");
1229         return 0;
1230     }
1231 }
1232 #endif