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 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 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 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 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 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, rem, w = r / 8;
1075
1076     assert(r < (25 * sizeof(A[0][0])) && (r % 8) == 0);
1077
1078     while (len >= r) {
1079         for (i = 0; i < w; i++) {
1080             uint64_t Ai = BitDeinterleave(A_flat[i]);
1081
1082             out[0] = (unsigned char)(Ai);
1083             out[1] = (unsigned char)(Ai >> 8);
1084             out[2] = (unsigned char)(Ai >> 16);
1085             out[3] = (unsigned char)(Ai >> 24);
1086             out[4] = (unsigned char)(Ai >> 32);
1087             out[5] = (unsigned char)(Ai >> 40);
1088             out[6] = (unsigned char)(Ai >> 48);
1089             out[7] = (unsigned char)(Ai >> 56);
1090             out += 8;
1091         }
1092         len -= r;
1093         if (len)
1094             KeccakF1600(A);
1095     }
1096
1097     rem = len % 8;
1098     len /= 8;
1099
1100     for (i = 0; i < len; i++) {
1101         uint64_t Ai = BitDeinterleave(A_flat[i]);
1102
1103         out[0] = (unsigned char)(Ai);
1104         out[1] = (unsigned char)(Ai >> 8);
1105         out[2] = (unsigned char)(Ai >> 16);
1106         out[3] = (unsigned char)(Ai >> 24);
1107         out[4] = (unsigned char)(Ai >> 32);
1108         out[5] = (unsigned char)(Ai >> 40);
1109         out[6] = (unsigned char)(Ai >> 48);
1110         out[7] = (unsigned char)(Ai >> 56);
1111         out += 8;
1112     }
1113
1114     if (rem) {
1115         uint64_t Ai = BitDeinterleave(A_flat[i]);
1116
1117         for (i = 0; i < rem; i++) {
1118             *out++ = (unsigned char)Ai;
1119             Ai >>= 8;
1120         }
1121     }
1122 }
1123 #else
1124 size_t SHA3_absorb(uint64_t A[5][5], const unsigned char *inp, size_t len,
1125                    size_t r);
1126 void SHA3_squeeze(uint64_t A[5][5], unsigned char *out, size_t len, size_t r);
1127 #endif
1128
1129 #ifdef SELFTEST
1130 /*
1131  * Post-padding one-shot implementations would look as following:
1132  *
1133  * SHA3_224     SHA3_sponge(inp, len, out, 224/8, (1600-448)/8);
1134  * SHA3_256     SHA3_sponge(inp, len, out, 256/8, (1600-512)/8);
1135  * SHA3_384     SHA3_sponge(inp, len, out, 384/8, (1600-768)/8);
1136  * SHA3_512     SHA3_sponge(inp, len, out, 512/8, (1600-1024)/8);
1137  * SHAKE_128    SHA3_sponge(inp, len, out, d, (1600-256)/8);
1138  * SHAKE_256    SHA3_sponge(inp, len, out, d, (1600-512)/8);
1139  */
1140
1141 void SHA3_sponge(const unsigned char *inp, size_t len,
1142                  unsigned char *out, size_t d, size_t r)
1143 {
1144     uint64_t A[5][5];
1145
1146     memset(A, 0, sizeof(A));
1147     SHA3_absorb(A, inp, len, r);
1148     SHA3_squeeze(A, out, d, r);
1149 }
1150
1151 # include <stdio.h>
1152
1153 int main()
1154 {
1155     /*
1156      * This is 5-bit SHAKE128 test from http://csrc.nist.gov/groups/ST/toolkit/examples.html#aHashing
1157      */
1158     unsigned char test[168] = { '\xf3', '\x3' };
1159     unsigned char out[512];
1160     size_t i;
1161     static const unsigned char result[512] = {
1162         0x2E, 0x0A, 0xBF, 0xBA, 0x83, 0xE6, 0x72, 0x0B,
1163         0xFB, 0xC2, 0x25, 0xFF, 0x6B, 0x7A, 0xB9, 0xFF,
1164         0xCE, 0x58, 0xBA, 0x02, 0x7E, 0xE3, 0xD8, 0x98,
1165         0x76, 0x4F, 0xEF, 0x28, 0x7D, 0xDE, 0xCC, 0xCA,
1166         0x3E, 0x6E, 0x59, 0x98, 0x41, 0x1E, 0x7D, 0xDB,
1167         0x32, 0xF6, 0x75, 0x38, 0xF5, 0x00, 0xB1, 0x8C,
1168         0x8C, 0x97, 0xC4, 0x52, 0xC3, 0x70, 0xEA, 0x2C,
1169         0xF0, 0xAF, 0xCA, 0x3E, 0x05, 0xDE, 0x7E, 0x4D,
1170         0xE2, 0x7F, 0xA4, 0x41, 0xA9, 0xCB, 0x34, 0xFD,
1171         0x17, 0xC9, 0x78, 0xB4, 0x2D, 0x5B, 0x7E, 0x7F,
1172         0x9A, 0xB1, 0x8F, 0xFE, 0xFF, 0xC3, 0xC5, 0xAC,
1173         0x2F, 0x3A, 0x45, 0x5E, 0xEB, 0xFD, 0xC7, 0x6C,
1174         0xEA, 0xEB, 0x0A, 0x2C, 0xCA, 0x22, 0xEE, 0xF6,
1175         0xE6, 0x37, 0xF4, 0xCA, 0xBE, 0x5C, 0x51, 0xDE,
1176         0xD2, 0xE3, 0xFA, 0xD8, 0xB9, 0x52, 0x70, 0xA3,
1177         0x21, 0x84, 0x56, 0x64, 0xF1, 0x07, 0xD1, 0x64,
1178         0x96, 0xBB, 0x7A, 0xBF, 0xBE, 0x75, 0x04, 0xB6,
1179         0xED, 0xE2, 0xE8, 0x9E, 0x4B, 0x99, 0x6F, 0xB5,
1180         0x8E, 0xFD, 0xC4, 0x18, 0x1F, 0x91, 0x63, 0x38,
1181         0x1C, 0xBE, 0x7B, 0xC0, 0x06, 0xA7, 0xA2, 0x05,
1182         0x98, 0x9C, 0x52, 0x6C, 0xD1, 0xBD, 0x68, 0x98,
1183         0x36, 0x93, 0xB4, 0xBD, 0xC5, 0x37, 0x28, 0xB2,
1184         0x41, 0xC1, 0xCF, 0xF4, 0x2B, 0xB6, 0x11, 0x50,
1185         0x2C, 0x35, 0x20, 0x5C, 0xAB, 0xB2, 0x88, 0x75,
1186         0x56, 0x55, 0xD6, 0x20, 0xC6, 0x79, 0x94, 0xF0,
1187         0x64, 0x51, 0x18, 0x7F, 0x6F, 0xD1, 0x7E, 0x04,
1188         0x66, 0x82, 0xBA, 0x12, 0x86, 0x06, 0x3F, 0xF8,
1189         0x8F, 0xE2, 0x50, 0x8D, 0x1F, 0xCA, 0xF9, 0x03,
1190         0x5A, 0x12, 0x31, 0xAD, 0x41, 0x50, 0xA9, 0xC9,
1191         0xB2, 0x4C, 0x9B, 0x2D, 0x66, 0xB2, 0xAD, 0x1B,
1192         0xDE, 0x0B, 0xD0, 0xBB, 0xCB, 0x8B, 0xE0, 0x5B,
1193         0x83, 0x52, 0x29, 0xEF, 0x79, 0x19, 0x73, 0x73,
1194         0x23, 0x42, 0x44, 0x01, 0xE1, 0xD8, 0x37, 0xB6,
1195         0x6E, 0xB4, 0xE6, 0x30, 0xFF, 0x1D, 0xE7, 0x0C,
1196         0xB3, 0x17, 0xC2, 0xBA, 0xCB, 0x08, 0x00, 0x1D,
1197         0x34, 0x77, 0xB7, 0xA7, 0x0A, 0x57, 0x6D, 0x20,
1198         0x86, 0x90, 0x33, 0x58, 0x9D, 0x85, 0xA0, 0x1D,
1199         0xDB, 0x2B, 0x66, 0x46, 0xC0, 0x43, 0xB5, 0x9F,
1200         0xC0, 0x11, 0x31, 0x1D, 0xA6, 0x66, 0xFA, 0x5A,
1201         0xD1, 0xD6, 0x38, 0x7F, 0xA9, 0xBC, 0x40, 0x15,
1202         0xA3, 0x8A, 0x51, 0xD1, 0xDA, 0x1E, 0xA6, 0x1D,
1203         0x64, 0x8D, 0xC8, 0xE3, 0x9A, 0x88, 0xB9, 0xD6,
1204         0x22, 0xBD, 0xE2, 0x07, 0xFD, 0xAB, 0xC6, 0xF2,
1205         0x82, 0x7A, 0x88, 0x0C, 0x33, 0x0B, 0xBF, 0x6D,
1206         0xF7, 0x33, 0x77, 0x4B, 0x65, 0x3E, 0x57, 0x30,
1207         0x5D, 0x78, 0xDC, 0xE1, 0x12, 0xF1, 0x0A, 0x2C,
1208         0x71, 0xF4, 0xCD, 0xAD, 0x92, 0xED, 0x11, 0x3E,
1209         0x1C, 0xEA, 0x63, 0xB9, 0x19, 0x25, 0xED, 0x28,
1210         0x19, 0x1E, 0x6D, 0xBB, 0xB5, 0xAA, 0x5A, 0x2A,
1211         0xFD, 0xA5, 0x1F, 0xC0, 0x5A, 0x3A, 0xF5, 0x25,
1212         0x8B, 0x87, 0x66, 0x52, 0x43, 0x55, 0x0F, 0x28,
1213         0x94, 0x8A, 0xE2, 0xB8, 0xBE, 0xB6, 0xBC, 0x9C,
1214         0x77, 0x0B, 0x35, 0xF0, 0x67, 0xEA, 0xA6, 0x41,
1215         0xEF, 0xE6, 0x5B, 0x1A, 0x44, 0x90, 0x9D, 0x1B,
1216         0x14, 0x9F, 0x97, 0xEE, 0xA6, 0x01, 0x39, 0x1C,
1217         0x60, 0x9E, 0xC8, 0x1D, 0x19, 0x30, 0xF5, 0x7C,
1218         0x18, 0xA4, 0xE0, 0xFA, 0xB4, 0x91, 0xD1, 0xCA,
1219         0xDF, 0xD5, 0x04, 0x83, 0x44, 0x9E, 0xDC, 0x0F,
1220         0x07, 0xFF, 0xB2, 0x4D, 0x2C, 0x6F, 0x9A, 0x9A,
1221         0x3B, 0xFF, 0x39, 0xAE, 0x3D, 0x57, 0xF5, 0x60,
1222         0x65, 0x4D, 0x7D, 0x75, 0xC9, 0x08, 0xAB, 0xE6,
1223         0x25, 0x64, 0x75, 0x3E, 0xAC, 0x39, 0xD7, 0x50,
1224         0x3D, 0xA6, 0xD3, 0x7C, 0x2E, 0x32, 0xE1, 0xAF,
1225         0x3B, 0x8A, 0xEC, 0x8A, 0xE3, 0x06, 0x9C, 0xD9
1226     };
1227
1228     test[167] = '\x80';
1229     SHA3_sponge(test, sizeof(test), out, sizeof(out), sizeof(test));
1230
1231     /*
1232      * Rationale behind keeping output [formatted as below] is that
1233      * one should be able to redirect it to a file, then copy-n-paste
1234      * final "output val" from official example to another file, and
1235      * compare the two with diff(1).
1236      */
1237     for (i = 0; i < sizeof(out);) {
1238         printf("%02X", out[i]);
1239         printf(++i % 16 && i != sizeof(out) ? " " : "\n");
1240     }
1241
1242     if (memcmp(out,result,sizeof(out))) {
1243         fprintf(stderr,"failure\n");
1244         return 1;
1245     } else {
1246         fprintf(stderr,"success\n");
1247         return 0;
1248     }
1249 }
1250 #endif