IDA_new/gf-complete/src/gf_w64.c

   1 /*
   2  * GF-Complete: A Comprehensive Open Source Library for Galois Field Arithmetic
   3  * James S. Plank, Ethan L. Miller, Kevin M. Greenan,
   4  * Benjamin A. Arnold, John A. Burnum, Adam W. Disney, Allen C. McBride.
   5  *
   6  * gf_w64.c
   7  *
   8  * Routines for 64-bit Galois fields
   9  */
  10
  11 #include "gf_int.h"
  12 #include <stdio.h>
  13 #include <stdlib.h>
  14 #include "gf_w64.h"
  15 #include "gf_cpu.h"
  16
  17 static
  18 inline
  19 gf_val_64_t gf_w64_inverse_from_divide (gf_t *gf, gf_val_64_t a)
  20 {
  21   return gf->divide.w64(gf, 1, a);
  22 }
  23
  24 #define MM_PRINT8(s, r) { uint8_t blah[16], ii; printf("%-12s", s); _mm_storeu_si128((__m128i *)blah, r); for (ii = 0; ii < 16; ii += 1) printf("%s%02x", (ii%4==0) ? "   " : " ", blah[15-ii]); printf("\n"); }
  25
  26 static
  27 inline
  28 gf_val_64_t gf_w64_divide_from_inverse (gf_t *gf, gf_val_64_t a, gf_val_64_t b)
  29 {
  30   b = gf->inverse.w64(gf, b);
  31   return gf->multiply.w64(gf, a, b);
  32 }
  33
  34 static
  35 void
  36 gf_w64_multiply_region_from_single(gf_t *gf, void *src, void *dest, gf_val_64_t val, int bytes, int
  37 xor)
  38 {
  39   uint32_t i;
  40   gf_val_64_t *s64;
  41   gf_val_64_t *d64;
  42
  43   s64 = (gf_val_64_t *) src;
  44   d64 = (gf_val_64_t *) dest;
  45
  46   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
  47   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
  48
  49   if (xor) {
  50     for (i = 0; i < bytes/sizeof(gf_val_64_t); i++) {
  51       d64[i] ^= gf->multiply.w64(gf, val, s64[i]);
  52     }
  53   } else {
  54     for (i = 0; i < bytes/sizeof(gf_val_64_t); i++) {
  55       d64[i] = gf->multiply.w64(gf, val, s64[i]);
  56     }
  57   }
  58 }
  59
  60 #if defined(INTEL_SSE4_PCLMUL)
  61 static
  62 void
  63 gf_w64_clm_multiply_region_from_single_2(gf_t *gf, void *src, void *dest, gf_val_64_t val, int bytes, int
  64 xor)
  65 {
  66   gf_val_64_t *s64, *d64, *top;
  67   gf_region_data rd;
  68
  69   __m128i         a, b;
  70   __m128i         result, r1;
  71   __m128i         prim_poly;
  72   __m128i         w;
  73   __m128i         m1, m3, m4;
  74   gf_internal_t * h = gf->scratch;
  75
  76   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
  77   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
  78
  79   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 16);
  80   gf_do_initial_region_alignment(&rd);
  81
  82   prim_poly = _mm_set_epi32(0, 0, 0, (uint32_t)(h->prim_poly & 0xffffffffULL));
  83   b = _mm_insert_epi64 (_mm_setzero_si128(), val, 0);
  84   m1 = _mm_set_epi32(0, 0, 0, (uint32_t)0xffffffff);
  85   m3 = _mm_slli_si128(m1, 8);
  86   m4 = _mm_slli_si128(m3, 4);
  87
  88   s64 = (gf_val_64_t *) rd.s_start;
  89   d64 = (gf_val_64_t *) rd.d_start;
  90   top = (gf_val_64_t *) rd.d_top;
  91
  92   if (xor) {
  93     while (d64 != top) {
  94       a = _mm_load_si128((__m128i *) s64);
  95       result = _mm_clmulepi64_si128 (a, b, 1);
  96
  97       w = _mm_clmulepi64_si128 (_mm_and_si128(result, m4), prim_poly, 1);
  98       result = _mm_xor_si128 (result, w);
  99       w = _mm_clmulepi64_si128 (_mm_and_si128(result, m3), prim_poly, 1);
 100       r1 = _mm_xor_si128 (result, w);
 101
 102       result = _mm_clmulepi64_si128 (a, b, 0);
 103
 104       w = _mm_clmulepi64_si128 (_mm_and_si128(result, m4), prim_poly, 1);
 105       result = _mm_xor_si128 (result, w);
 106
 107       w = _mm_clmulepi64_si128 (_mm_and_si128(result, m3), prim_poly, 1);
 108       result = _mm_xor_si128 (result, w);
 109
 110       result = _mm_unpacklo_epi64(result, r1);
 111
 112       r1 = _mm_load_si128((__m128i *) d64);
 113       result = _mm_xor_si128(r1, result);
 114       _mm_store_si128((__m128i *) d64, result);
 115       d64 += 2;
 116       s64 += 2;
 117     }
 118   } else {
 119     while (d64 != top) {
 120
 121       a = _mm_load_si128((__m128i *) s64);
 122       result = _mm_clmulepi64_si128 (a, b, 1);
 123
 124       w = _mm_clmulepi64_si128 (_mm_and_si128(result, m4), prim_poly, 1);
 125       result = _mm_xor_si128 (result, w);
 126       w = _mm_clmulepi64_si128 (_mm_and_si128(result, m3), prim_poly, 1);
 127       r1 = _mm_xor_si128 (result, w);
 128
 129       result = _mm_clmulepi64_si128 (a, b, 0);
 130
 131       w = _mm_clmulepi64_si128 (_mm_and_si128(result, m4), prim_poly, 1);
 132       result = _mm_xor_si128 (result, w);
 133       w = _mm_clmulepi64_si128 (_mm_and_si128(result, m3), prim_poly, 1);
 134       result = _mm_xor_si128 (result, w);
 135
 136       result = _mm_unpacklo_epi64(result, r1);
 137
 138       _mm_store_si128((__m128i *) d64, result);
 139       d64 += 2;
 140       s64 += 2;
 141     }
 142   }
 143   gf_do_final_region_alignment(&rd);
 144 }
 145 #endif
 146
 147 #if defined(INTEL_SSE4_PCLMUL)
 148 static
 149 void
 150 gf_w64_clm_multiply_region_from_single_4(gf_t *gf, void *src, void *dest, gf_val_64_t val, int bytes, int
 151 xor)
 152 {
 153   gf_val_64_t *s64, *d64, *top;
 154   gf_region_data rd;
 155
 156   __m128i         a, b;
 157   __m128i         result, r1;
 158   __m128i         prim_poly;
 159   __m128i         w;
 160   __m128i         m1, m3, m4;
 161   gf_internal_t * h = gf->scratch;
 162
 163   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
 164   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
 165
 166   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 16);
 167   gf_do_initial_region_alignment(&rd);
 168
 169   prim_poly = _mm_set_epi32(0, 0, 0, (uint32_t)(h->prim_poly & 0xffffffffULL));
 170   b = _mm_insert_epi64 (_mm_setzero_si128(), val, 0);
 171   m1 = _mm_set_epi32(0, 0, 0, (uint32_t)0xffffffff);
 172   m3 = _mm_slli_si128(m1, 8);
 173   m4 = _mm_slli_si128(m3, 4);
 174
 175   s64 = (gf_val_64_t *) rd.s_start;
 176   d64 = (gf_val_64_t *) rd.d_start;
 177   top = (gf_val_64_t *) rd.d_top;
 178
 179   if (xor) {
 180     while (d64 != top) {
 181       a = _mm_load_si128((__m128i *) s64);
 182       result = _mm_clmulepi64_si128 (a, b, 1);
 183
 184       w = _mm_clmulepi64_si128 (_mm_and_si128(result, m4), prim_poly, 1);
 185       result = _mm_xor_si128 (result, w);
 186       w = _mm_clmulepi64_si128 (_mm_and_si128(result, m3), prim_poly, 1);
 187       r1 = _mm_xor_si128 (result, w);
 188
 189       result = _mm_clmulepi64_si128 (a, b, 0);
 190
 191       w = _mm_clmulepi64_si128 (_mm_and_si128(result, m4), prim_poly, 1);
 192       result = _mm_xor_si128 (result, w);
 193
 194       w = _mm_clmulepi64_si128 (_mm_and_si128(result, m3), prim_poly, 1);
 195       result = _mm_xor_si128 (result, w);
 196
 197       result = _mm_unpacklo_epi64(result, r1);
 198
 199       r1 = _mm_load_si128((__m128i *) d64);
 200       result = _mm_xor_si128(r1, result);
 201       _mm_store_si128((__m128i *) d64, result);
 202       d64 += 2;
 203       s64 += 2;
 204     }
 205   } else {
 206     while (d64 != top) {
 207       a = _mm_load_si128((__m128i *) s64);
 208       result = _mm_clmulepi64_si128 (a, b, 1);
 209
 210       w = _mm_clmulepi64_si128 (_mm_and_si128(result, m4), prim_poly, 1);
 211       result = _mm_xor_si128 (result, w);
 212       w = _mm_clmulepi64_si128 (_mm_and_si128(result, m3), prim_poly, 1);
 213       r1 = _mm_xor_si128 (result, w);
 214
 215       result = _mm_clmulepi64_si128 (a, b, 0);
 216
 217       w = _mm_clmulepi64_si128 (_mm_and_si128(result, m4), prim_poly, 1);
 218       result = _mm_xor_si128 (result, w);
 219       w = _mm_clmulepi64_si128 (_mm_and_si128(result, m3), prim_poly, 1);
 220       result = _mm_xor_si128 (result, w);
 221
 222       result = _mm_unpacklo_epi64(result, r1);
 223
 224       _mm_store_si128((__m128i *) d64, result);
 225       d64 += 2;
 226       s64 += 2;
 227     }
 228   }
 229   gf_do_final_region_alignment(&rd);
 230 }
 231 #endif
 232
 233 static
 234   inline
 235 gf_val_64_t gf_w64_euclid (gf_t *gf, gf_val_64_t b)
 236 {
 237   gf_val_64_t e_i, e_im1, e_ip1;
 238   gf_val_64_t d_i, d_im1, d_ip1;
 239   gf_val_64_t y_i, y_im1, y_ip1;
 240   gf_val_64_t c_i;
 241   gf_val_64_t one = 1;
 242
 243   if (b == 0) return -1;
 244   e_im1 = ((gf_internal_t *) (gf->scratch))->prim_poly;
 245   e_i = b;
 246   d_im1 = 64;
 247   for (d_i = d_im1-1; ((one << d_i) & e_i) == 0; d_i--) ;
 248   y_i = 1;
 249   y_im1 = 0;
 250
 251   while (e_i != 1) {
 252
 253     e_ip1 = e_im1;
 254     d_ip1 = d_im1;
 255     c_i = 0;
 256
 257     while (d_ip1 >= d_i) {
 258       c_i ^= (one << (d_ip1 - d_i));
 259       e_ip1 ^= (e_i << (d_ip1 - d_i));
 260       d_ip1--;
 261       if (e_ip1 == 0) return 0;
 262       while ((e_ip1 & (one << d_ip1)) == 0) d_ip1--;
 263     }
 264
 265     y_ip1 = y_im1 ^ gf->multiply.w64(gf, c_i, y_i);
 266     y_im1 = y_i;
 267     y_i = y_ip1;
 268
 269     e_im1 = e_i;
 270     d_im1 = d_i;
 271     e_i = e_ip1;
 272     d_i = d_ip1;
 273   }
 274
 275   return y_i;
 276 }
 277
 278 /* JSP: GF_MULT_SHIFT: The world's dumbest multiplication algorithm.  I only
 279    include it for completeness.  It does have the feature that it requires no
 280    extra memory.
 281 */
 282
 283 static
 284 inline
 285 gf_val_64_t
 286 gf_w64_shift_multiply (gf_t *gf, gf_val_64_t a64, gf_val_64_t b64)
 287 {
 288   uint64_t pl, pr, ppl, ppr, i, a, bl, br, one, lbit;
 289   gf_internal_t *h;
 290
 291   h = (gf_internal_t *) gf->scratch;
 292
 293   /* Allen: set leading one of primitive polynomial */
 294
 295   a = a64;
 296   bl = 0;
 297   br = b64;
 298   one = 1;
 299   lbit = (one << 63);
 300
 301   pl = 0; /* Allen: left side of product */
 302   pr = 0; /* Allen: right side of product */
 303
 304   /* Allen: unlike the corresponding functions for smaller word sizes,
 305    * this loop carries out the initial carryless multiply by
 306    * shifting b itself rather than simply looking at successively
 307    * higher shifts of b */
 308
 309   for (i = 0; i < GF_FIELD_WIDTH; i++) {
 310     if (a & (one << i)) {
 311       pl ^= bl;
 312       pr ^= br;
 313     }
 314
 315     bl <<= 1;
 316     if (br & lbit) bl ^= 1;
 317     br <<= 1;
 318   }
 319
 320   /* Allen: the name of the variable "one" is no longer descriptive at this point */
 321
 322   one = lbit >> 1;
 323   ppl = (h->prim_poly >> 2) | one;
 324   ppr = (h->prim_poly << (GF_FIELD_WIDTH-2));
 325   while (one != 0) {
 326     if (pl & one) {
 327       pl ^= ppl;
 328       pr ^= ppr;
 329     }
 330     one >>= 1;
 331     ppr >>= 1;
 332     if (ppl & 1) ppr ^= lbit;
 333     ppl >>= 1;
 334   }
 335   return pr;
 336 }
 337
 338 /*
 339  * ELM: Use the Intel carryless multiply instruction to do very fast 64x64 multiply.
 340  */
 341
 342 #if defined(INTEL_SSE4_PCLMUL)
 343
 344 static
 345 inline
 346 gf_val_64_t
 347 gf_w64_clm_multiply_2 (gf_t *gf, gf_val_64_t a64, gf_val_64_t b64)
 348 {
 349        gf_val_64_t rv = 0;
 350
 351         __m128i         a, b;
 352         __m128i         result;
 353         __m128i         prim_poly;
 354         __m128i         v, w;
 355         gf_internal_t * h = gf->scratch;
 356
 357         a = _mm_insert_epi64 (_mm_setzero_si128(), a64, 0);
 358         b = _mm_insert_epi64 (a, b64, 0);
 359         prim_poly = _mm_set_epi32(0, 0, 0, (uint32_t)(h->prim_poly & 0xffffffffULL));
 360         /* Do the initial multiply */
 361
 362         result = _mm_clmulepi64_si128 (a, b, 0);
 363
 364         /* Mask off the high order 32 bits using subtraction of the polynomial.
 365          * NOTE: this part requires that the polynomial have at least 32 leading 0 bits.
 366          */
 367
 368         /* Adam: We cant include the leading one in the 64 bit pclmul,
 369          so we need to split up the high 8 bytes of the result into two
 370          parts before we multiply them with the prim_poly.*/
 371
 372         v = _mm_insert_epi32 (_mm_srli_si128 (result, 8), 0, 0);
 373         w = _mm_clmulepi64_si128 (prim_poly, v, 0);
 374         result = _mm_xor_si128 (result, w);
 375         v = _mm_insert_epi32 (_mm_srli_si128 (result, 8), 0, 1);
 376         w = _mm_clmulepi64_si128 (prim_poly, v, 0);
 377         result = _mm_xor_si128 (result, w);
 378
 379         rv = ((gf_val_64_t)_mm_extract_epi64(result, 0));
 380         return rv;
 381 }
 382 #endif
 383
 384 #if defined(INTEL_SSE4_PCLMUL)
 385
 386 static
 387 inline
 388 gf_val_64_t
 389 gf_w64_clm_multiply_4 (gf_t *gf, gf_val_64_t a64, gf_val_64_t b64)
 390 {
 391   gf_val_64_t rv = 0;
 392
 393   __m128i         a, b;
 394   __m128i         result;
 395   __m128i         prim_poly;
 396   __m128i         v, w;
 397   gf_internal_t * h = gf->scratch;
 398
 399   a = _mm_insert_epi64 (_mm_setzero_si128(), a64, 0);
 400   b = _mm_insert_epi64 (a, b64, 0);
 401   prim_poly = _mm_set_epi32(0, 0, 0, (uint32_t)(h->prim_poly & 0xffffffffULL));
 402
 403   /* Do the initial multiply */
 404
 405   result = _mm_clmulepi64_si128 (a, b, 0);
 406
 407   v = _mm_insert_epi32 (_mm_srli_si128 (result, 8), 0, 0);
 408   w = _mm_clmulepi64_si128 (prim_poly, v, 0);
 409   result = _mm_xor_si128 (result, w);
 410   v = _mm_insert_epi32 (_mm_srli_si128 (result, 8), 0, 1);
 411   w = _mm_clmulepi64_si128 (prim_poly, v, 0);
 412   result = _mm_xor_si128 (result, w);
 413
 414   v = _mm_insert_epi32 (_mm_srli_si128 (result, 8), 0, 0);
 415   w = _mm_clmulepi64_si128 (prim_poly, v, 0);
 416   result = _mm_xor_si128 (result, w);
 417   v = _mm_insert_epi32 (_mm_srli_si128 (result, 8), 0, 1);
 418   w = _mm_clmulepi64_si128 (prim_poly, v, 0);
 419   result = _mm_xor_si128 (result, w);
 420
 421   rv = ((gf_val_64_t)_mm_extract_epi64(result, 0));
 422   return rv;
 423 }
 424 #endif
 425
 426
 427 #if defined(INTEL_SSE4_PCLMUL)
 428   void
 429 gf_w64_clm_multiply_region(gf_t *gf, void *src, void *dest, uint64_t val, int bytes, int xor)
 430 {
 431   gf_internal_t *h;
 432   uint8_t *s8, *d8, *dtop;
 433   gf_region_data rd;
 434   __m128i  v, b, m, prim_poly, c, fr, w, result;
 435
 436   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
 437   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
 438
 439   h = (gf_internal_t *) gf->scratch;
 440
 441   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 16);
 442   gf_do_initial_region_alignment(&rd);
 443
 444   s8 = (uint8_t *) rd.s_start;
 445   d8 = (uint8_t *) rd.d_start;
 446   dtop = (uint8_t *) rd.d_top;
 447
 448   v = _mm_insert_epi64(_mm_setzero_si128(), val, 0);
 449   m = _mm_set_epi32(0, 0, 0xffffffff, 0xffffffff);
 450   prim_poly = _mm_set_epi32(0, 0, 0, (uint32_t)(h->prim_poly & 0xffffffffULL));
 451
 452   if (xor) {
 453     while (d8 != dtop) {
 454       b = _mm_load_si128((__m128i *) s8);
 455       result = _mm_clmulepi64_si128 (b, v, 0);
 456       c = _mm_insert_epi32 (_mm_srli_si128 (result, 8), 0, 0);
 457       w = _mm_clmulepi64_si128 (prim_poly, c, 0);
 458       result = _mm_xor_si128 (result, w);
 459       c = _mm_insert_epi32 (_mm_srli_si128 (result, 8), 0, 1);
 460       w = _mm_clmulepi64_si128 (prim_poly, c, 0);
 461       fr = _mm_xor_si128 (result, w);
 462       fr = _mm_and_si128 (fr, m);
 463
 464       result = _mm_clmulepi64_si128 (b, v, 1);
 465       c = _mm_insert_epi32 (_mm_srli_si128 (result, 8), 0, 0);
 466       w = _mm_clmulepi64_si128 (prim_poly, c, 0);
 467       result = _mm_xor_si128 (result, w);
 468       c = _mm_insert_epi32 (_mm_srli_si128 (result, 8), 0, 1);
 469       w = _mm_clmulepi64_si128 (prim_poly, c, 0);
 470       result = _mm_xor_si128 (result, w);
 471       result = _mm_slli_si128 (result, 8);
 472       fr = _mm_xor_si128 (result, fr);
 473       result = _mm_load_si128((__m128i *) d8);
 474       fr = _mm_xor_si128 (result, fr);
 475
 476       _mm_store_si128((__m128i *) d8, fr);
 477       d8 += 16;
 478       s8 += 16;
 479     }
 480   } else {
 481     while (d8 < dtop) {
 482       b = _mm_load_si128((__m128i *) s8);
 483       result = _mm_clmulepi64_si128 (b, v, 0);
 484       c = _mm_insert_epi32 (_mm_srli_si128 (result, 8), 0, 0);
 485       w = _mm_clmulepi64_si128 (prim_poly, c, 0);
 486       result = _mm_xor_si128 (result, w);
 487       c = _mm_insert_epi32 (_mm_srli_si128 (result, 8), 0, 1);
 488       w = _mm_clmulepi64_si128 (prim_poly, c, 0);
 489       fr = _mm_xor_si128 (result, w);
 490       fr = _mm_and_si128 (fr, m);
 491
 492       result = _mm_clmulepi64_si128 (b, v, 1);
 493       c = _mm_insert_epi32 (_mm_srli_si128 (result, 8), 0, 0);
 494       w = _mm_clmulepi64_si128 (prim_poly, c, 0);
 495       result = _mm_xor_si128 (result, w);
 496       c = _mm_insert_epi32 (_mm_srli_si128 (result, 8), 0, 1);
 497       w = _mm_clmulepi64_si128 (prim_poly, c, 0);
 498       result = _mm_xor_si128 (result, w);
 499       result = _mm_slli_si128 (result, 8);
 500       fr = _mm_xor_si128 (result, fr);
 501
 502       _mm_store_si128((__m128i *) d8, fr);
 503       d8 += 16;
 504       s8 += 16;
 505     }
 506   }
 507   gf_do_final_region_alignment(&rd);
 508 }
 509 #endif
 510
 511 void
 512 gf_w64_split_4_64_lazy_multiply_region(gf_t *gf, void *src, void *dest, uint64_t val, int bytes, int xor)
 513 {
 514   gf_internal_t *h;
 515   struct gf_split_4_64_lazy_data *ld;
 516   int i, j, k;
 517   uint64_t pp, v, s, *s64, *d64, *top;
 518   gf_region_data rd;
 519
 520   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
 521   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
 522
 523   h = (gf_internal_t *) gf->scratch;
 524   pp = h->prim_poly;
 525
 526   ld = (struct gf_split_4_64_lazy_data *) h->private;
 527
 528   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 8);
 529   gf_do_initial_region_alignment(&rd);
 530
 531   if (ld->last_value != val) {
 532     v = val;
 533     for (i = 0; i < 16; i++) {
 534       ld->tables[i][0] = 0;
 535       for (j = 1; j < 16; j <<= 1) {
 536         for (k = 0; k < j; k++) {
 537           ld->tables[i][k^j] = (v ^ ld->tables[i][k]);
 538         }
 539         v = (v & GF_FIRST_BIT) ? ((v << 1) ^ pp) : (v << 1);
 540       }
 541     }
 542   }
 543   ld->last_value = val;
 544
 545   s64 = (uint64_t *) rd.s_start;
 546   d64 = (uint64_t *) rd.d_start;
 547   top = (uint64_t *) rd.d_top;
 548
 549   while (d64 != top) {
 550     v = (xor) ? *d64 : 0;
 551     s = *s64;
 552     i = 0;
 553     while (s != 0) {
 554       v ^= ld->tables[i][s&0xf];
 555       s >>= 4;
 556       i++;
 557     }
 558     *d64 = v;
 559     d64++;
 560     s64++;
 561   }
 562   gf_do_final_region_alignment(&rd);
 563 }
 564
 565 static
 566 inline
 567 uint64_t
 568 gf_w64_split_8_8_multiply (gf_t *gf, uint64_t a64, uint64_t b64)
 569 {
 570   uint64_t product, i, j, mask, tb;
 571   gf_internal_t *h;
 572   struct gf_split_8_8_data *d8;
 573
 574   h = (gf_internal_t *) gf->scratch;
 575   d8 = (struct gf_split_8_8_data *) h->private;
 576   product = 0;
 577   mask = 0xff;
 578
 579   for (i = 0; a64 != 0; i++) {
 580     tb = b64;
 581     for (j = 0; tb != 0; j++) {
 582       product ^= d8->tables[i+j][a64&mask][tb&mask];
 583       tb >>= 8;
 584     }
 585     a64 >>= 8;
 586   }
 587   return product;
 588 }
 589
 590 void
 591 gf_w64_split_8_64_lazy_multiply_region(gf_t *gf, void *src, void *dest, uint64_t val, int bytes, int xor)
 592 {
 593   gf_internal_t *h;
 594   struct gf_split_8_64_lazy_data *ld;
 595   int i, j, k;
 596   uint64_t pp, v, s, *s64, *d64, *top;
 597   gf_region_data rd;
 598
 599   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
 600   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
 601
 602   h = (gf_internal_t *) gf->scratch;
 603   pp = h->prim_poly;
 604
 605   ld = (struct gf_split_8_64_lazy_data *) h->private;
 606
 607   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 4);
 608   gf_do_initial_region_alignment(&rd);
 609
 610   if (ld->last_value != val) {
 611     v = val;
 612     for (i = 0; i < 8; i++) {
 613       ld->tables[i][0] = 0;
 614       for (j = 1; j < 256; j <<= 1) {
 615         for (k = 0; k < j; k++) {
 616           ld->tables[i][k^j] = (v ^ ld->tables[i][k]);
 617         }
 618         v = (v & GF_FIRST_BIT) ? ((v << 1) ^ pp) : (v << 1);
 619       }
 620     }
 621   }
 622   ld->last_value = val;
 623
 624   s64 = (uint64_t *) rd.s_start;
 625   d64 = (uint64_t *) rd.d_start;
 626   top = (uint64_t *) rd.d_top;
 627
 628   while (d64 != top) {
 629     v = (xor) ? *d64 : 0;
 630     s = *s64;
 631     i = 0;
 632     while (s != 0) {
 633       v ^= ld->tables[i][s&0xff];
 634       s >>= 8;
 635       i++;
 636     }
 637     *d64 = v;
 638     d64++;
 639     s64++;
 640   }
 641   gf_do_final_region_alignment(&rd);
 642 }
 643
 644 void
 645 gf_w64_split_16_64_lazy_multiply_region(gf_t *gf, void *src, void *dest, uint64_t val, int bytes, int xor)
 646 {
 647   gf_internal_t *h;
 648   struct gf_split_16_64_lazy_data *ld;
 649   int i, j, k;
 650   uint64_t pp, v, s, *s64, *d64, *top;
 651   gf_region_data rd;
 652
 653   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
 654   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
 655
 656   h = (gf_internal_t *) gf->scratch;
 657   pp = h->prim_poly;
 658
 659   ld = (struct gf_split_16_64_lazy_data *) h->private;
 660
 661   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 4);
 662   gf_do_initial_region_alignment(&rd);
 663
 664   if (ld->last_value != val) {
 665     v = val;
 666     for (i = 0; i < 4; i++) {
 667       ld->tables[i][0] = 0;
 668       for (j = 1; j < (1<<16); j <<= 1) {
 669         for (k = 0; k < j; k++) {
 670           ld->tables[i][k^j] = (v ^ ld->tables[i][k]);
 671         }
 672         v = (v & GF_FIRST_BIT) ? ((v << 1) ^ pp) : (v << 1);
 673       }
 674     }
 675   }
 676   ld->last_value = val;
 677
 678   s64 = (uint64_t *) rd.s_start;
 679   d64 = (uint64_t *) rd.d_start;
 680   top = (uint64_t *) rd.d_top;
 681
 682   while (d64 != top) {
 683     v = (xor) ? *d64 : 0;
 684     s = *s64;
 685     i = 0;
 686     while (s != 0) {
 687       v ^= ld->tables[i][s&0xffff];
 688       s >>= 16;
 689       i++;
 690     }
 691     *d64 = v;
 692     d64++;
 693     s64++;
 694   }
 695   gf_do_final_region_alignment(&rd);
 696 }
 697
 698 static
 699 int gf_w64_shift_init(gf_t *gf)
 700 {
 701   SET_FUNCTION(gf,multiply,w64,gf_w64_shift_multiply)
 702   SET_FUNCTION(gf,inverse,w64,gf_w64_euclid)
 703   SET_FUNCTION(gf,multiply_region,w64,gf_w64_multiply_region_from_single)
 704   return 1;
 705 }
 706
 707 static
 708 int gf_w64_cfm_init(gf_t *gf)
 709 {
 710   SET_FUNCTION(gf,inverse,w64,gf_w64_euclid)
 711   SET_FUNCTION(gf,multiply_region,w64,gf_w64_multiply_region_from_single)
 712
 713 #if defined(INTEL_SSE4_PCLMUL)
 714   if (gf_cpu_supports_intel_pclmul) {
 715     gf_internal_t *h;
 716
 717     h = (gf_internal_t *) gf->scratch;
 718
 719     if ((0xfffffffe00000000ULL & h->prim_poly) == 0){
 720       SET_FUNCTION(gf,multiply,w64,gf_w64_clm_multiply_2)
 721       SET_FUNCTION(gf,multiply_region,w64,gf_w64_clm_multiply_region_from_single_2)
 722     }else if((0xfffe000000000000ULL & h->prim_poly) == 0){
 723       SET_FUNCTION(gf,multiply,w64,gf_w64_clm_multiply_4)
 724       SET_FUNCTION(gf,multiply_region,w64,gf_w64_clm_multiply_region_from_single_4)
 725     } else {
 726       return 0;
 727     }
 728     return 1;
 729   }
 730 #endif
 731
 732   return 0;
 733 }
 734
 735 static
 736 void
 737 gf_w64_group_set_shift_tables(uint64_t *shift, uint64_t val, gf_internal_t *h)
 738 {
 739   uint64_t i;
 740   uint64_t j;
 741   uint64_t one = 1;
 742   int g_s;
 743
 744   g_s = h->arg1;
 745   shift[0] = 0;
 746
 747   for (i = 1; i < ((uint64_t)1 << g_s); i <<= 1) {
 748     for (j = 0; j < i; j++) shift[i|j] = shift[j]^val;
 749     if (val & (one << 63)) {
 750       val <<= 1;
 751       val ^= h->prim_poly;
 752     } else {
 753       val <<= 1;
 754     }
 755   }
 756 }
 757
 758 static
 759 inline
 760 gf_val_64_t
 761 gf_w64_group_multiply(gf_t *gf, gf_val_64_t a, gf_val_64_t b)
 762 {
 763   uint64_t top, bot, mask, tp;
 764   int g_s, g_r, lshift, rshift;
 765   struct gf_w64_group_data *gd;
 766
 767   gf_internal_t *h = (gf_internal_t *) gf->scratch;
 768   g_s = h->arg1;
 769   g_r = h->arg2;
 770   gd = (struct gf_w64_group_data *) h->private;
 771   gf_w64_group_set_shift_tables(gd->shift, b, h);
 772
 773   mask = (((uint64_t)1 << g_s) - 1);
 774   top = 0;
 775   bot = gd->shift[a&mask];
 776   a >>= g_s;
 777
 778   if (a == 0) return bot;
 779   lshift = 0;
 780   rshift = 64;
 781
 782   do {              /* Shifting out is straightfoward */
 783     lshift += g_s;
 784     rshift -= g_s;
 785     tp = gd->shift[a&mask];
 786     top ^= (tp >> rshift);
 787     bot ^= (tp << lshift);
 788     a >>= g_s;
 789   } while (a != 0);
 790
 791   /* Reducing is a bit gross, because I don't zero out the index bits of top.
 792      The reason is that we throw top away.  Even better, that last (tp >> rshift)
 793      is going to be ignored, so it doesn't matter how (tp >> 64) is implemented. */
 794
 795   lshift = ((lshift-1) / g_r) * g_r;
 796   rshift = 64 - lshift;
 797   mask = ((uint64_t)1 << g_r) - 1;
 798   while (lshift >= 0) {
 799     tp = gd->reduce[(top >> lshift) & mask];
 800     top ^= (tp >> rshift);
 801     bot ^= (tp << lshift);
 802     lshift -= g_r;
 803     rshift += g_r;
 804   }
 805
 806   return bot;
 807 }
 808
 809 static
 810 void gf_w64_group_multiply_region(gf_t *gf, void *src, void *dest, gf_val_64_t val, int bytes, int xor)
 811 {
 812   int i, fzb;
 813   uint64_t a64, smask, rmask, top, bot, tp;
 814   int lshift, rshift, g_s, g_r;
 815   gf_region_data rd;
 816   uint64_t *s64, *d64, *dtop;
 817   struct gf_w64_group_data *gd;
 818   gf_internal_t *h = (gf_internal_t *) gf->scratch;
 819
 820   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
 821   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
 822
 823   gd = (struct gf_w64_group_data *) h->private;
 824   g_s = h->arg1;
 825   g_r = h->arg2;
 826   gf_w64_group_set_shift_tables(gd->shift, val, h);
 827
 828   for (i = 63; !(val & (1ULL << i)); i--) ;
 829   i += g_s;
 830
 831   /* i is the bit position of the first zero bit in any element of
 832                            gd->shift[] */
 833
 834   if (i > 64) i = 64;
 835
 836   fzb = i;
 837
 838   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 4);
 839
 840   gf_do_initial_region_alignment(&rd);
 841
 842   s64 = (uint64_t *) rd.s_start;
 843   d64 = (uint64_t *) rd.d_start;
 844   dtop = (uint64_t *) rd.d_top;
 845
 846   smask = ((uint64_t)1 << g_s) - 1;
 847   rmask = ((uint64_t)1 << g_r) - 1;
 848
 849   while (d64 < dtop) {
 850     a64 = *s64;
 851
 852     top = 0;
 853     bot = gd->shift[a64&smask];
 854     a64 >>= g_s;
 855     i = fzb;
 856
 857     if (a64 != 0) {
 858       lshift = 0;
 859       rshift = 64;
 860
 861       do {
 862         lshift += g_s;
 863         rshift -= g_s;
 864         tp = gd->shift[a64&smask];
 865         top ^= (tp >> rshift);
 866         bot ^= (tp << lshift);
 867         a64 >>= g_s;
 868       } while (a64 != 0);
 869       i += lshift;
 870
 871       lshift = ((i-64-1) / g_r) * g_r;
 872       rshift = 64 - lshift;
 873       while (lshift >= 0) {
 874         tp = gd->reduce[(top >> lshift) & rmask];
 875         top ^= (tp >> rshift);
 876         bot ^= (tp << lshift);
 877         lshift -= g_r;
 878         rshift += g_r;
 879       }
 880     }
 881
 882     if (xor) bot ^= *d64;
 883     *d64 = bot;
 884     d64++;
 885     s64++;
 886   }
 887   gf_do_final_region_alignment(&rd);
 888 }
 889
 890 static
 891 inline
 892 gf_val_64_t
 893 gf_w64_group_s_equals_r_multiply(gf_t *gf, gf_val_64_t a, gf_val_64_t b)
 894 {
 895   int leftover, rs;
 896   uint64_t p, l, ind, a64;
 897   int bits_left;
 898   int g_s;
 899
 900   struct gf_w64_group_data *gd;
 901   gf_internal_t *h = (gf_internal_t *) gf->scratch;
 902   g_s = h->arg1;
 903
 904   gd = (struct gf_w64_group_data *) h->private;
 905   gf_w64_group_set_shift_tables(gd->shift, b, h);
 906
 907   leftover = 64 % g_s;
 908   if (leftover == 0) leftover = g_s;
 909
 910   rs = 64 - leftover;
 911   a64 = a;
 912   ind = a64 >> rs;
 913   a64 <<= leftover;
 914   p = gd->shift[ind];
 915
 916   bits_left = rs;
 917   rs = 64 - g_s;
 918
 919   while (bits_left > 0) {
 920     bits_left -= g_s;
 921     ind = a64 >> rs;
 922     a64 <<= g_s;
 923     l = p >> rs;
 924     p = (gd->shift[ind] ^ gd->reduce[l] ^ (p << g_s));
 925   }
 926   return p;
 927 }
 928
 929 static
 930 void gf_w64_group_s_equals_r_multiply_region(gf_t *gf, void *src, void *dest, gf_val_64_t val, int bytes, int xor)
 931 {
 932   int leftover, rs;
 933   uint64_t p, l, ind, a64;
 934   int bits_left;
 935   int g_s;
 936   gf_region_data rd;
 937   uint64_t *s64, *d64, *top;
 938   struct gf_w64_group_data *gd;
 939   gf_internal_t *h = (gf_internal_t *) gf->scratch;
 940
 941   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
 942   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
 943
 944   gd = (struct gf_w64_group_data *) h->private;
 945   g_s = h->arg1;
 946   gf_w64_group_set_shift_tables(gd->shift, val, h);
 947
 948   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 4);
 949   gf_do_initial_region_alignment(&rd);
 950
 951   s64 = (uint64_t *) rd.s_start;
 952   d64 = (uint64_t *) rd.d_start;
 953   top = (uint64_t *) rd.d_top;
 954
 955   leftover = 64 % g_s;
 956   if (leftover == 0) leftover = g_s;
 957
 958   while (d64 < top) {
 959     rs = 64 - leftover;
 960     a64 = *s64;
 961     ind = a64 >> rs;
 962     a64 <<= leftover;
 963     p = gd->shift[ind];
 964
 965     bits_left = rs;
 966     rs = 64 - g_s;
 967
 968     while (bits_left > 0) {
 969       bits_left -= g_s;
 970       ind = a64 >> rs;
 971       a64 <<= g_s;
 972       l = p >> rs;
 973       p = (gd->shift[ind] ^ gd->reduce[l] ^ (p << g_s));
 974     }
 975     if (xor) p ^= *d64;
 976     *d64 = p;
 977     d64++;
 978     s64++;
 979   }
 980   gf_do_final_region_alignment(&rd);
 981 }
 982
 983
 984 static
 985 int gf_w64_group_init(gf_t *gf)
 986 {
 987   uint64_t i, j, p, index;
 988   struct gf_w64_group_data *gd;
 989   gf_internal_t *h = (gf_internal_t *) gf->scratch;
 990   uint64_t g_r, g_s;
 991
 992   g_s = h->arg1;
 993   g_r = h->arg2;
 994
 995   gd = (struct gf_w64_group_data *) h->private;
 996   gd->shift = (uint64_t *) (&(gd->memory));
 997   gd->reduce = gd->shift + (1 << g_s);
 998
 999   gd->reduce[0] = 0;
1000   for (i = 0; i < ((uint64_t)1 << g_r); i++) {
1001     p = 0;
1002     index = 0;
1003     for (j = 0; j < g_r; j++) {
1004       if (i & (1 << j)) {
1005         p ^= (h->prim_poly << j);
1006         index ^= (1 << j);
1007         if (j > 0) index ^= (h->prim_poly >> (64-j));
1008       }
1009     }
1010     gd->reduce[index] = p;
1011   }
1012
1013   if (g_s == g_r) {
1014     SET_FUNCTION(gf,multiply,w64,gf_w64_group_s_equals_r_multiply)
1015     SET_FUNCTION(gf,multiply_region,w64,gf_w64_group_s_equals_r_multiply_region)
1016   } else {
1017     SET_FUNCTION(gf,multiply,w64,gf_w64_group_multiply)
1018     SET_FUNCTION(gf,multiply_region,w64,gf_w64_group_multiply_region)
1019   }
1020   SET_FUNCTION(gf,divide,w64,NULL)
1021   SET_FUNCTION(gf,inverse,w64,gf_w64_euclid)
1022
1023   return 1;
1024 }
1025
1026 static
1027 gf_val_64_t gf_w64_extract_word(gf_t *gf, void *start, int bytes, int index)
1028 {
1029   uint64_t *r64, rv;
1030
1031   r64 = (uint64_t *) start;
1032   rv = r64[index];
1033   return rv;
1034 }
1035
1036 static
1037 gf_val_64_t gf_w64_composite_extract_word(gf_t *gf, void *start, int bytes, int index)
1038 {
1039   int sub_size;
1040   gf_internal_t *h;
1041   uint8_t *r8, *top;
1042   uint64_t a, b, *r64;
1043   gf_region_data rd;
1044
1045   h = (gf_internal_t *) gf->scratch;
1046   gf_set_region_data(&rd, gf, start, start, bytes, 0, 0, 32);
1047   r64 = (uint64_t *) start;
1048   if (r64 + index < (uint64_t *) rd.d_start) return r64[index];
1049   if (r64 + index >= (uint64_t *) rd.d_top) return r64[index];
1050   index -= (((uint64_t *) rd.d_start) - r64);
1051   r8 = (uint8_t *) rd.d_start;
1052   top = (uint8_t *) rd.d_top;
1053   sub_size = (top-r8)/2;
1054
1055   a = h->base_gf->extract_word.w32(h->base_gf, r8, sub_size, index);
1056   b = h->base_gf->extract_word.w32(h->base_gf, r8+sub_size, sub_size, index);
1057   return (a | ((uint64_t)b << 32));
1058 }
1059
1060 static
1061 gf_val_64_t gf_w64_split_extract_word(gf_t *gf, void *start, int bytes, int index)
1062 {
1063   int i;
1064   uint64_t *r64, rv;
1065   uint8_t *r8;
1066   gf_region_data rd;
1067
1068   gf_set_region_data(&rd, gf, start, start, bytes, 0, 0, 128);
1069   r64 = (uint64_t *) start;
1070   if (r64 + index < (uint64_t *) rd.d_start) return r64[index];
1071   if (r64 + index >= (uint64_t *) rd.d_top) return r64[index];
1072   index -= (((uint64_t *) rd.d_start) - r64);
1073   r8 = (uint8_t *) rd.d_start;
1074   r8 += ((index & 0xfffffff0)*8);
1075   r8 += (index & 0xf);
1076   r8 += 112;
1077   rv =0;
1078   for (i = 0; i < 8; i++) {
1079     rv <<= 8;
1080     rv |= *r8;
1081     r8 -= 16;
1082   }
1083   return rv;
1084 }
1085
1086 static
1087 inline
1088 gf_val_64_t
1089 gf_w64_bytwo_b_multiply (gf_t *gf, gf_val_64_t a, gf_val_64_t b)
1090 {
1091   uint64_t prod, pp, bmask;
1092   gf_internal_t *h;
1093
1094   h = (gf_internal_t *) gf->scratch;
1095   pp = h->prim_poly;
1096
1097   prod = 0;
1098   bmask = 0x8000000000000000ULL;
1099
1100   while (1) {
1101     if (a & 1) prod ^= b;
1102     a >>= 1;
1103     if (a == 0) return prod;
1104     if (b & bmask) {
1105       b = ((b << 1) ^ pp);
1106     } else {
1107       b <<= 1;
1108     }
1109   }
1110 }
1111
1112 static
1113 inline
1114 gf_val_64_t
1115 gf_w64_bytwo_p_multiply (gf_t *gf, gf_val_64_t a, gf_val_64_t b)
1116 {
1117   uint64_t prod, pp, pmask, amask;
1118   gf_internal_t *h;
1119
1120   h = (gf_internal_t *) gf->scratch;
1121   pp = h->prim_poly;
1122
1123   prod = 0;
1124
1125   /* changed from declare then shift to just declare.*/
1126
1127   pmask = 0x8000000000000000ULL;
1128   amask = 0x8000000000000000ULL;
1129
1130   while (amask != 0) {
1131     if (prod & pmask) {
1132       prod = ((prod << 1) ^ pp);
1133     } else {
1134       prod <<= 1;
1135     }
1136     if (a & amask) prod ^= b;
1137     amask >>= 1;
1138   }
1139   return prod;
1140 }
1141
1142 static
1143 void
1144 gf_w64_bytwo_p_nosse_multiply_region(gf_t *gf, void *src, void *dest, gf_val_64_t val, int bytes, int xor)
1145 {
1146   uint64_t *s64, *d64, ta, prod, amask, pmask, pp;
1147   gf_region_data rd;
1148   gf_internal_t *h;
1149
1150   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
1151   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
1152
1153   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 8);
1154   gf_do_initial_region_alignment(&rd);
1155
1156   h = (gf_internal_t *) gf->scratch;
1157
1158   s64 = (uint64_t *) rd.s_start;
1159   d64 = (uint64_t *) rd.d_start;
1160   pmask = 0x80000000;
1161   pmask <<= 32;
1162   pp = h->prim_poly;
1163
1164   if (xor) {
1165     while (s64 < (uint64_t *) rd.s_top) {
1166       prod = 0;
1167       amask = pmask;
1168       ta = *s64;
1169       while (amask != 0) {
1170         prod = (prod & pmask) ? ((prod << 1) ^ pp) : (prod << 1);
1171         if (val & amask) prod ^= ta;
1172         amask >>= 1;
1173       }
1174       *d64 ^= prod;
1175       d64++;
1176       s64++;
1177     }
1178   } else {
1179     while (s64 < (uint64_t *) rd.s_top) {
1180       prod = 0;
1181       amask = pmask;
1182       ta = *s64;
1183       while (amask != 0) {
1184         prod = (prod & pmask) ? ((prod << 1) ^ pp) : (prod << 1);
1185         if (val & amask) prod ^= ta;
1186         amask >>= 1;
1187       }
1188       *d64 = prod;
1189       d64++;
1190       s64++;
1191     }
1192   }
1193   gf_do_final_region_alignment(&rd);
1194 }
1195
1196 static
1197 void
1198 gf_w64_bytwo_b_nosse_multiply_region(gf_t *gf, void *src, void *dest, gf_val_64_t val, int bytes, int xor)
1199 {
1200   uint64_t *s64, *d64, ta, tb, prod, bmask, pp;
1201   gf_region_data rd;
1202   gf_internal_t *h;
1203
1204   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
1205   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
1206
1207   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 8);
1208   gf_do_initial_region_alignment(&rd);
1209
1210   h = (gf_internal_t *) gf->scratch;
1211
1212   s64 = (uint64_t *) rd.s_start;
1213   d64 = (uint64_t *) rd.d_start;
1214   bmask = 0x80000000;
1215   bmask <<= 32;
1216   pp = h->prim_poly;
1217
1218   if (xor) {
1219     while (s64 < (uint64_t *) rd.s_top) {
1220       prod = 0;
1221       tb = val;
1222       ta = *s64;
1223       while (1) {
1224         if (tb & 1) prod ^= ta;
1225         tb >>= 1;
1226         if (tb == 0) break;
1227         ta = (ta & bmask) ? ((ta << 1) ^ pp) : (ta << 1);
1228       }
1229       *d64 ^= prod;
1230       d64++;
1231       s64++;
1232     }
1233   } else {
1234     while (s64 < (uint64_t *) rd.s_top) {
1235       prod = 0;
1236       tb = val;
1237       ta = *s64;
1238       while (1) {
1239         if (tb & 1) prod ^= ta;
1240         tb >>= 1;
1241         if (tb == 0) break;
1242         ta = (ta & bmask) ? ((ta << 1) ^ pp) : (ta << 1);
1243       }
1244       *d64 = prod;
1245       d64++;
1246       s64++;
1247     }
1248   }
1249   gf_do_final_region_alignment(&rd);
1250 }
1251
1252 #define SSE_AB2(pp, m1 ,m2, va, t1, t2) {\
1253           t1 = _mm_and_si128(_mm_slli_epi64(va, 1), m1); \
1254           t2 = _mm_and_si128(va, m2); \
1255           t2 = _mm_sub_epi64 (_mm_slli_epi64(t2, 1), _mm_srli_epi64(t2, (GF_FIELD_WIDTH-1))); \
1256           va = _mm_xor_si128(t1, _mm_and_si128(t2, pp)); }
1257
1258 #define BYTWO_P_ONESTEP {\
1259       SSE_AB2(pp, m1 ,m2, prod, t1, t2); \
1260       t1 = _mm_and_si128(v, one); \
1261       t1 = _mm_sub_epi64(t1, one); \
1262       t1 = _mm_and_si128(t1, ta); \
1263       prod = _mm_xor_si128(prod, t1); \
1264       v = _mm_srli_epi64(v, 1); }
1265
1266
1267 #ifdef INTEL_SSE2
1268 void gf_w64_bytwo_p_sse_multiply_region(gf_t *gf, void *src, void *dest, gf_val_64_t val, int bytes, int xor)
1269 {
1270   int i;
1271   uint8_t *s8, *d8;
1272   uint64_t vrev, one64;
1273   uint64_t amask;
1274   __m128i pp, m1, m2, ta, prod, t1, t2, tp, one, v;
1275   gf_region_data rd;
1276   gf_internal_t *h;
1277
1278   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
1279   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
1280
1281   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 16);
1282   gf_do_initial_region_alignment(&rd);
1283
1284   h = (gf_internal_t *) gf->scratch;
1285   one64 = 1;
1286   vrev = 0;
1287   for (i = 0; i < 64; i++) {
1288     vrev <<= 1;
1289     if (!(val & (one64 << i))) vrev |= 1;
1290   }
1291
1292   s8 = (uint8_t *) rd.s_start;
1293   d8 = (uint8_t *) rd.d_start;
1294
1295   amask = -1;
1296   amask ^= 1;
1297   pp = _mm_set1_epi64x(h->prim_poly);
1298   m1 = _mm_set1_epi64x(amask);
1299   m2 = _mm_set1_epi64x(one64 << 63);
1300   one = _mm_set1_epi64x(1);
1301
1302   while (d8 < (uint8_t *) rd.d_top) {
1303     prod = _mm_setzero_si128();
1304     v = _mm_set1_epi64x(vrev);
1305     ta = _mm_load_si128((__m128i *) s8);
1306     tp = (!xor) ? _mm_setzero_si128() : _mm_load_si128((__m128i *) d8);
1307     BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP;
1308     BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP;
1309     BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP;
1310     BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP;
1311     BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP;
1312     BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP;
1313     BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP;
1314     BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP;
1315     BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP;
1316     BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP;
1317     BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP;
1318     BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP;
1319     BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP;
1320     BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP;
1321     BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP;
1322     BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP; BYTWO_P_ONESTEP;
1323     _mm_store_si128((__m128i *) d8, _mm_xor_si128(prod, tp));
1324     d8 += 16;
1325     s8 += 16;
1326   }
1327   gf_do_final_region_alignment(&rd);
1328 }
1329 #endif
1330
1331 #ifdef INTEL_SSE2
1332 static
1333 void
1334 gf_w64_bytwo_b_sse_region_2_xor(gf_region_data *rd)
1335 {
1336   uint64_t one64, amask;
1337   uint8_t *d8, *s8;
1338   __m128i pp, m1, m2, t1, t2, va, vb;
1339   gf_internal_t *h;
1340
1341   s8 = (uint8_t *) rd->s_start;
1342   d8 = (uint8_t *) rd->d_start;
1343
1344   h = (gf_internal_t *) rd->gf->scratch;
1345   one64 = 1;
1346   amask = -1;
1347   amask ^= 1;
1348   pp = _mm_set1_epi64x(h->prim_poly);
1349   m1 = _mm_set1_epi64x(amask);
1350   m2 = _mm_set1_epi64x(one64 << 63);
1351
1352   while (d8 < (uint8_t *) rd->d_top) {
1353     va = _mm_load_si128 ((__m128i *)(s8));
1354     SSE_AB2(pp, m1, m2, va, t1, t2);
1355     vb = _mm_load_si128 ((__m128i *)(d8));
1356     vb = _mm_xor_si128(vb, va);
1357     _mm_store_si128((__m128i *)d8, vb);
1358     d8 += 16;
1359     s8 += 16;
1360   }
1361 }
1362 #endif
1363
1364 #ifdef INTEL_SSE2
1365 static
1366 void
1367 gf_w64_bytwo_b_sse_region_2_noxor(gf_region_data *rd)
1368 {
1369   uint64_t one64, amask;
1370   uint8_t *d8, *s8;
1371   __m128i pp, m1, m2, t1, t2, va;
1372   gf_internal_t *h;
1373
1374   s8 = (uint8_t *) rd->s_start;
1375   d8 = (uint8_t *) rd->d_start;
1376
1377   h = (gf_internal_t *) rd->gf->scratch;
1378   one64 = 1;
1379   amask = -1;
1380   amask ^= 1;
1381   pp = _mm_set1_epi64x(h->prim_poly);
1382   m1 = _mm_set1_epi64x(amask);
1383   m2 = _mm_set1_epi64x(one64 << 63);
1384
1385   while (d8 < (uint8_t *) rd->d_top) {
1386     va = _mm_load_si128 ((__m128i *)(s8));
1387     SSE_AB2(pp, m1, m2, va, t1, t2);
1388     _mm_store_si128((__m128i *)d8, va);
1389     d8 += 16;
1390     s8 += 16;
1391   }
1392 }
1393 #endif
1394
1395 #ifdef INTEL_SSE2
1396 static
1397 void
1398 gf_w64_bytwo_b_sse_multiply_region(gf_t *gf, void *src, void *dest, gf_val_64_t val, int bytes, int xor)
1399 {
1400   uint64_t itb, amask, one64;
1401   uint8_t *d8, *s8;
1402   __m128i pp, m1, m2, t1, t2, va, vb;
1403   gf_region_data rd;
1404   gf_internal_t *h;
1405
1406   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
1407   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
1408
1409   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 16);
1410   gf_do_initial_region_alignment(&rd);
1411
1412   if (val == 2) {
1413     if (xor) {
1414       gf_w64_bytwo_b_sse_region_2_xor(&rd);
1415     } else {
1416       gf_w64_bytwo_b_sse_region_2_noxor(&rd);
1417     }
1418     gf_do_final_region_alignment(&rd);
1419     return;
1420   }
1421
1422   s8 = (uint8_t *) rd.s_start;
1423   d8 = (uint8_t *) rd.d_start;
1424   h = (gf_internal_t *) gf->scratch;
1425
1426   one64 = 1;
1427   amask = -1;
1428   amask ^= 1;
1429   pp = _mm_set1_epi64x(h->prim_poly);
1430   m1 = _mm_set1_epi64x(amask);
1431   m2 = _mm_set1_epi64x(one64 << 63);
1432
1433   while (d8 < (uint8_t *) rd.d_top) {
1434     va = _mm_load_si128 ((__m128i *)(s8));
1435     vb = (!xor) ? _mm_setzero_si128() : _mm_load_si128 ((__m128i *)(d8));
1436     itb = val;
1437     while (1) {
1438       if (itb & 1) vb = _mm_xor_si128(vb, va);
1439       itb >>= 1;
1440       if (itb == 0) break;
1441       SSE_AB2(pp, m1, m2, va, t1, t2);
1442     }
1443     _mm_store_si128((__m128i *)d8, vb);
1444     d8 += 16;
1445     s8 += 16;
1446   }
1447
1448   gf_do_final_region_alignment(&rd);
1449 }
1450 #endif
1451
1452
1453 static
1454 int gf_w64_bytwo_init(gf_t *gf)
1455 {
1456   gf_internal_t *h;
1457
1458   h = (gf_internal_t *) gf->scratch;
1459
1460   if (h->mult_type == GF_MULT_BYTWO_p) {
1461     SET_FUNCTION(gf,multiply,w64,gf_w64_bytwo_p_multiply)
1462     #ifdef INTEL_SSE2
1463       if (gf_cpu_supports_intel_sse2 && !(h->region_type & GF_REGION_NOSIMD)) {
1464         SET_FUNCTION(gf,multiply_region,w64,gf_w64_bytwo_p_sse_multiply_region)
1465       } else {
1466     #endif
1467         SET_FUNCTION(gf,multiply_region,w64,gf_w64_bytwo_p_nosse_multiply_region)
1468         if(h->region_type & GF_REGION_SIMD)
1469           return 0;
1470     #ifdef INTEL_SSE2
1471       }
1472     #endif
1473   } else {
1474     SET_FUNCTION(gf,multiply,w64,gf_w64_bytwo_b_multiply)
1475     #ifdef INTEL_SSE2
1476       if (gf_cpu_supports_intel_sse2 && !(h->region_type & GF_REGION_NOSIMD)) {
1477         SET_FUNCTION(gf,multiply_region,w64,gf_w64_bytwo_b_sse_multiply_region)
1478       } else {
1479     #endif
1480       SET_FUNCTION(gf,multiply_region,w64,gf_w64_bytwo_b_nosse_multiply_region)
1481       if(h->region_type & GF_REGION_SIMD)
1482         return 0;
1483     #ifdef INTEL_SSE2
1484       }
1485     #endif
1486   }
1487   SET_FUNCTION(gf,inverse,w64,gf_w64_euclid)
1488   return 1;
1489 }
1490
1491
1492 static
1493 gf_val_64_t
1494 gf_w64_composite_multiply(gf_t *gf, gf_val_64_t a, gf_val_64_t b)
1495 {
1496   gf_internal_t *h = (gf_internal_t *) gf->scratch;
1497   gf_t *base_gf = h->base_gf;
1498   uint32_t b0 = b & 0x00000000ffffffff;
1499   uint32_t b1 = (b & 0xffffffff00000000) >> 32;
1500   uint32_t a0 = a & 0x00000000ffffffff;
1501   uint32_t a1 = (a & 0xffffffff00000000) >> 32;
1502   uint32_t a1b1;
1503
1504   a1b1 = base_gf->multiply.w32(base_gf, a1, b1);
1505
1506   return ((uint64_t)(base_gf->multiply.w32(base_gf, a0, b0) ^ a1b1) |
1507          ((uint64_t)(base_gf->multiply.w32(base_gf, a1, b0) ^ base_gf->multiply.w32(base_gf, a0, b1) ^ base_gf->multiply.w32(base_gf, a1b1, h->prim_poly)) << 32));
1508 }
1509
1510 /*
1511  * Composite field division trick (explained in 2007 tech report)
1512  *
1513  * Compute a / b = a*b^-1, where p(x) = x^2 + sx + 1
1514  *
1515  * let c = b^-1
1516  *
1517  * c*b = (s*b1c1+b1c0+b0c1)x+(b1c1+b0c0)
1518  *
1519  * want (s*b1c1+b1c0+b0c1) = 0 and (b1c1+b0c0) = 1
1520  *
1521  * let d = b1c1 and d+1 = b0c0
1522  *
1523  * solve s*b1c1+b1c0+b0c1 = 0
1524  *
1525  * solution: d = (b1b0^-1)(b1b0^-1+b0b1^-1+s)^-1
1526  *
1527  * c0 = (d+1)b0^-1
1528  * c1 = d*b1^-1
1529  *
1530  * a / b = a * c
1531  */
1532
1533 static
1534 gf_val_64_t
1535 gf_w64_composite_inverse(gf_t *gf, gf_val_64_t a)
1536 {
1537   gf_internal_t *h = (gf_internal_t *) gf->scratch;
1538   gf_t *base_gf = h->base_gf;
1539   uint32_t a0 = a & 0x00000000ffffffff;
1540   uint32_t a1 = (a & 0xffffffff00000000) >> 32;
1541   uint32_t c0, c1, d, tmp;
1542   uint64_t c;
1543   uint32_t a0inv, a1inv;
1544
1545   if (a0 == 0) {
1546     a1inv = base_gf->inverse.w32(base_gf, a1);
1547     c0 = base_gf->multiply.w32(base_gf, a1inv, h->prim_poly);
1548     c1 = a1inv;
1549   } else if (a1 == 0) {
1550     c0 = base_gf->inverse.w32(base_gf, a0);
1551     c1 = 0;
1552   } else {
1553     a1inv = base_gf->inverse.w32(base_gf, a1);
1554     a0inv = base_gf->inverse.w32(base_gf, a0);
1555
1556     d = base_gf->multiply.w32(base_gf, a1, a0inv);
1557
1558     tmp = (base_gf->multiply.w32(base_gf, a1, a0inv) ^ base_gf->multiply.w32(base_gf, a0, a1inv) ^ h->prim_poly);
1559     tmp = base_gf->inverse.w32(base_gf, tmp);
1560
1561     d = base_gf->multiply.w32(base_gf, d, tmp);
1562
1563     c0 = base_gf->multiply.w32(base_gf, (d^1), a0inv);
1564     c1 = base_gf->multiply.w32(base_gf, d, a1inv);
1565   }
1566
1567   c = c0 | ((uint64_t)c1 << 32);
1568
1569   return c;
1570 }
1571
1572 static
1573 void
1574 gf_w64_composite_multiply_region(gf_t *gf, void *src, void *dest, gf_val_64_t val, int bytes, int xor)
1575 {
1576   gf_internal_t *h = (gf_internal_t *) gf->scratch;
1577   gf_t *base_gf = h->base_gf;
1578   uint32_t b0 = val & 0x00000000ffffffff;
1579   uint32_t b1 = (val & 0xffffffff00000000) >> 32;
1580   uint64_t *s64, *d64;
1581   uint64_t *top;
1582   uint64_t a0, a1, a1b1;
1583   gf_region_data rd;
1584
1585   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
1586   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 8);
1587
1588   s64 = rd.s_start;
1589   d64 = rd.d_start;
1590   top = rd.d_top;
1591
1592   if (xor) {
1593     while (d64 < top) {
1594       a0 = *s64 & 0x00000000ffffffff;
1595       a1 = (*s64 & 0xffffffff00000000) >> 32;
1596       a1b1 = base_gf->multiply.w32(base_gf, a1, b1);
1597
1598       *d64 ^= ((uint64_t)(base_gf->multiply.w32(base_gf, a0, b0) ^ a1b1) |
1599                 ((uint64_t)(base_gf->multiply.w32(base_gf, a1, b0) ^ base_gf->multiply.w32(base_gf, a0, b1) ^ base_gf->multiply.w32(base_gf, a1b1, h->prim_poly)) << 32));
1600       s64++;
1601       d64++;
1602     }
1603   } else {
1604     while (d64 < top) {
1605       a0 = *s64 & 0x00000000ffffffff;
1606       a1 = (*s64 & 0xffffffff00000000) >> 32;
1607       a1b1 = base_gf->multiply.w32(base_gf, a1, b1);
1608
1609       *d64 = ((base_gf->multiply.w32(base_gf, a0, b0) ^ a1b1) |
1610                 ((uint64_t)(base_gf->multiply.w32(base_gf, a1, b0) ^ base_gf->multiply.w32(base_gf, a0, b1) ^ base_gf->multiply.w32(base_gf, a1b1, h->prim_poly)) << 32));
1611       s64++;
1612       d64++;
1613     }
1614   }
1615 }
1616
1617 static
1618 void
1619 gf_w64_composite_multiply_region_alt(gf_t *gf, void *src, void *dest, gf_val_64_t val, int bytes, int xor)
1620 {
1621   gf_internal_t *h = (gf_internal_t *) gf->scratch;
1622   gf_t *base_gf = h->base_gf;
1623   gf_val_32_t val0 = val & 0x00000000ffffffff;
1624   gf_val_32_t val1 = (val & 0xffffffff00000000) >> 32;
1625   uint8_t *slow, *shigh;
1626   uint8_t *dlow, *dhigh, *top;
1627   int sub_reg_size;
1628   gf_region_data rd;
1629
1630   if (!xor) {
1631     memset(dest, 0, bytes);
1632   }
1633
1634   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 32);
1635   gf_do_initial_region_alignment(&rd);
1636
1637   slow = (uint8_t *) rd.s_start;
1638   dlow = (uint8_t *) rd.d_start;
1639   top = (uint8_t*) rd.d_top;
1640   sub_reg_size = (top - dlow)/2;
1641   shigh = slow + sub_reg_size;
1642   dhigh = dlow + sub_reg_size;
1643
1644   base_gf->multiply_region.w32(base_gf, slow, dlow, val0, sub_reg_size, xor);
1645   base_gf->multiply_region.w32(base_gf, shigh, dlow, val1, sub_reg_size, 1);
1646   base_gf->multiply_region.w32(base_gf, slow, dhigh, val1, sub_reg_size, xor);
1647   base_gf->multiply_region.w32(base_gf, shigh, dhigh, val0, sub_reg_size, 1);
1648   base_gf->multiply_region.w32(base_gf, shigh, dhigh, base_gf->multiply.w32(base_gf, h->prim_poly, val1), sub_reg_size, 1);
1649
1650   gf_do_final_region_alignment(&rd);
1651 }
1652
1653
1654
1655 static
1656 int gf_w64_composite_init(gf_t *gf)
1657 {
1658   gf_internal_t *h = (gf_internal_t *) gf->scratch;
1659
1660   if (h->region_type & GF_REGION_ALTMAP) {
1661     SET_FUNCTION(gf,multiply_region,w64,gf_w64_composite_multiply_region_alt)
1662   } else {
1663     SET_FUNCTION(gf,multiply_region,w64,gf_w64_composite_multiply_region)
1664   }
1665
1666   SET_FUNCTION(gf,multiply,w64,gf_w64_composite_multiply)
1667   SET_FUNCTION(gf,divide,w64,NULL)
1668   SET_FUNCTION(gf,inverse,w64,gf_w64_composite_inverse)
1669
1670   return 1;
1671 }
1672
1673 #ifdef INTEL_SSSE3
1674 static
1675   void
1676 gf_w64_split_4_64_lazy_sse_altmap_multiply_region(gf_t *gf, void *src, void *dest, uint64_t val, int bytes, int xor)
1677 {
1678   gf_internal_t *h;
1679   int i, j, k;
1680   uint64_t pp, v, *s64, *d64, *top;
1681   __m128i si, tables[16][8], p[8], v0, mask1;
1682   struct gf_split_4_64_lazy_data *ld;
1683   uint8_t btable[16];
1684   gf_region_data rd;
1685
1686   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
1687   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
1688
1689   h = (gf_internal_t *) gf->scratch;
1690   pp = h->prim_poly;
1691
1692   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 128);
1693   gf_do_initial_region_alignment(&rd);
1694
1695   s64 = (uint64_t *) rd.s_start;
1696   d64 = (uint64_t *) rd.d_start;
1697   top = (uint64_t *) rd.d_top;
1698
1699   ld = (struct gf_split_4_64_lazy_data *) h->private;
1700
1701   v = val;
1702   for (i = 0; i < 16; i++) {
1703     ld->tables[i][0] = 0;
1704     for (j = 1; j < 16; j <<= 1) {
1705       for (k = 0; k < j; k++) {
1706         ld->tables[i][k^j] = (v ^ ld->tables[i][k]);
1707       }
1708       v = (v & GF_FIRST_BIT) ? ((v << 1) ^ pp) : (v << 1);
1709     }
1710     for (j = 0; j < 8; j++) {
1711       for (k = 0; k < 16; k++) {
1712         btable[k] = (uint8_t) ld->tables[i][k];
1713         ld->tables[i][k] >>= 8;
1714       }
1715       tables[i][j] = _mm_loadu_si128((__m128i *) btable);
1716     }
1717   }
1718
1719   mask1 = _mm_set1_epi8(0xf);
1720
1721   while (d64 != top) {
1722
1723     if (xor) {
1724       for (i = 0; i < 8; i++) p[i] = _mm_load_si128 ((__m128i *) (d64+i*2));
1725     } else {
1726       for (i = 0; i < 8; i++) p[i] = _mm_setzero_si128();
1727     }
1728     i = 0;
1729     for (k = 0; k < 8; k++) {
1730       v0 = _mm_load_si128((__m128i *) s64);
1731       /* MM_PRINT8("v", v0); */
1732       s64 += 2;
1733
1734       si = _mm_and_si128(v0, mask1);
1735
1736       for (j = 0; j < 8; j++) {
1737         p[j] = _mm_xor_si128(p[j], _mm_shuffle_epi8(tables[i][j], si));
1738       }
1739       i++;
1740       v0 = _mm_srli_epi32(v0, 4);
1741       si = _mm_and_si128(v0, mask1);
1742       for (j = 0; j < 8; j++) {
1743         p[j] = _mm_xor_si128(p[j], _mm_shuffle_epi8(tables[i][j], si));
1744       }
1745       i++;
1746     }
1747     for (i = 0; i < 8; i++) {
1748       /* MM_PRINT8("v", p[i]); */
1749       _mm_store_si128((__m128i *) d64, p[i]);
1750       d64 += 2;
1751     }
1752   }
1753   gf_do_final_region_alignment(&rd);
1754 }
1755 #endif
1756
1757 #ifdef INTEL_SSE4
1758 static
1759   void
1760 gf_w64_split_4_64_lazy_sse_multiply_region(gf_t *gf, void *src, void *dest, uint64_t val, int bytes, int xor)
1761 {
1762   gf_internal_t *h;
1763   int i, j, k;
1764   uint64_t pp, v, *s64, *d64, *top;
1765   __m128i si, tables[16][8], p[8], st[8], mask1, mask8, mask16, t1;
1766   struct gf_split_4_64_lazy_data *ld;
1767   uint8_t btable[16];
1768   gf_region_data rd;
1769
1770   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
1771   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
1772
1773   h = (gf_internal_t *) gf->scratch;
1774   pp = h->prim_poly;
1775
1776   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 128);
1777   gf_do_initial_region_alignment(&rd);
1778
1779   s64 = (uint64_t *) rd.s_start;
1780   d64 = (uint64_t *) rd.d_start;
1781   top = (uint64_t *) rd.d_top;
1782
1783   ld = (struct gf_split_4_64_lazy_data *) h->private;
1784
1785   v = val;
1786   for (i = 0; i < 16; i++) {
1787     ld->tables[i][0] = 0;
1788     for (j = 1; j < 16; j <<= 1) {
1789       for (k = 0; k < j; k++) {
1790         ld->tables[i][k^j] = (v ^ ld->tables[i][k]);
1791       }
1792       v = (v & GF_FIRST_BIT) ? ((v << 1) ^ pp) : (v << 1);
1793     }
1794     for (j = 0; j < 8; j++) {
1795       for (k = 0; k < 16; k++) {
1796         btable[k] = (uint8_t) ld->tables[i][k];
1797         ld->tables[i][k] >>= 8;
1798       }
1799       tables[i][j] = _mm_loadu_si128((__m128i *) btable);
1800     }
1801   }
1802
1803   mask1 = _mm_set1_epi8(0xf);
1804   mask8 = _mm_set1_epi16(0xff);
1805   mask16 = _mm_set1_epi32(0xffff);
1806
1807   while (d64 != top) {
1808
1809     for (i = 0; i < 8; i++) p[i] = _mm_setzero_si128();
1810
1811     for (k = 0; k < 8; k++) {
1812       st[k]  = _mm_load_si128((__m128i *) s64);
1813       s64 += 2;
1814     }
1815
1816     for (k = 0; k < 4; k ++) {
1817       st[k] = _mm_shuffle_epi32(st[k], _MM_SHUFFLE(3,1,2,0));
1818       st[k+4] = _mm_shuffle_epi32(st[k+4], _MM_SHUFFLE(2,0,3,1));
1819       t1 = _mm_blend_epi16(st[k], st[k+4], 0xf0);
1820       st[k] = _mm_srli_si128(st[k], 8);
1821       st[k+4] = _mm_slli_si128(st[k+4], 8);
1822       st[k+4] = _mm_blend_epi16(st[k], st[k+4], 0xf0);
1823       st[k] = t1;
1824     }
1825
1826 /*
1827     printf("After pack pass 1\n");
1828     for (k = 0; k < 8; k++) {
1829       MM_PRINT8("v", st[k]);
1830     }
1831     printf("\n");
1832  */
1833
1834     t1 = _mm_packus_epi32(_mm_and_si128(st[0], mask16), _mm_and_si128(st[2], mask16));
1835     st[2] = _mm_packus_epi32(_mm_srli_epi32(st[0], 16), _mm_srli_epi32(st[2], 16));
1836     st[0] = t1;
1837     t1 = _mm_packus_epi32(_mm_and_si128(st[1], mask16), _mm_and_si128(st[3], mask16));
1838     st[3] = _mm_packus_epi32(_mm_srli_epi32(st[1], 16), _mm_srli_epi32(st[3], 16));
1839     st[1] = t1;
1840     t1 = _mm_packus_epi32(_mm_and_si128(st[4], mask16), _mm_and_si128(st[6], mask16));
1841     st[6] = _mm_packus_epi32(_mm_srli_epi32(st[4], 16), _mm_srli_epi32(st[6], 16));
1842     st[4] = t1;
1843     t1 = _mm_packus_epi32(_mm_and_si128(st[5], mask16), _mm_and_si128(st[7], mask16));
1844     st[7] = _mm_packus_epi32(_mm_srli_epi32(st[5], 16), _mm_srli_epi32(st[7], 16));
1845     st[5] = t1;
1846
1847 /*
1848     printf("After pack pass 2\n");
1849     for (k = 0; k < 8; k++) {
1850       MM_PRINT8("v", st[k]);
1851     }
1852     printf("\n");
1853  */
1854     t1 = _mm_packus_epi16(_mm_and_si128(st[0], mask8), _mm_and_si128(st[1], mask8));
1855     st[1] = _mm_packus_epi16(_mm_srli_epi16(st[0], 8), _mm_srli_epi16(st[1], 8));
1856     st[0] = t1;
1857     t1 = _mm_packus_epi16(_mm_and_si128(st[2], mask8), _mm_and_si128(st[3], mask8));
1858     st[3] = _mm_packus_epi16(_mm_srli_epi16(st[2], 8), _mm_srli_epi16(st[3], 8));
1859     st[2] = t1;
1860     t1 = _mm_packus_epi16(_mm_and_si128(st[4], mask8), _mm_and_si128(st[5], mask8));
1861     st[5] = _mm_packus_epi16(_mm_srli_epi16(st[4], 8), _mm_srli_epi16(st[5], 8));
1862     st[4] = t1;
1863     t1 = _mm_packus_epi16(_mm_and_si128(st[6], mask8), _mm_and_si128(st[7], mask8));
1864     st[7] = _mm_packus_epi16(_mm_srli_epi16(st[6], 8), _mm_srli_epi16(st[7], 8));
1865     st[6] = t1;
1866
1867 /*
1868     printf("After final pack pass 2\n");
1869     for (k = 0; k < 8; k++) {
1870       MM_PRINT8("v", st[k]);
1871     }
1872  */
1873     i = 0;
1874     for (k = 0; k < 8; k++) {
1875       si = _mm_and_si128(st[k], mask1);
1876
1877       for (j = 0; j < 8; j++) {
1878         p[j] = _mm_xor_si128(p[j], _mm_shuffle_epi8(tables[i][j], si));
1879       }
1880       i++;
1881       st[k] = _mm_srli_epi32(st[k], 4);
1882       si = _mm_and_si128(st[k], mask1);
1883       for (j = 0; j < 8; j++) {
1884         p[j] = _mm_xor_si128(p[j], _mm_shuffle_epi8(tables[i][j], si));
1885       }
1886       i++;
1887     }
1888
1889     t1 = _mm_unpacklo_epi8(p[0], p[1]);
1890     p[1] = _mm_unpackhi_epi8(p[0], p[1]);
1891     p[0] = t1;
1892     t1 = _mm_unpacklo_epi8(p[2], p[3]);
1893     p[3] = _mm_unpackhi_epi8(p[2], p[3]);
1894     p[2] = t1;
1895     t1 = _mm_unpacklo_epi8(p[4], p[5]);
1896     p[5] = _mm_unpackhi_epi8(p[4], p[5]);
1897     p[4] = t1;
1898     t1 = _mm_unpacklo_epi8(p[6], p[7]);
1899     p[7] = _mm_unpackhi_epi8(p[6], p[7]);
1900     p[6] = t1;
1901
1902 /*
1903     printf("After unpack pass 1:\n");
1904     for (i = 0; i < 8; i++) {
1905       MM_PRINT8("v", p[i]);
1906     }
1907  */
1908
1909     t1 = _mm_unpacklo_epi16(p[0], p[2]);
1910     p[2] = _mm_unpackhi_epi16(p[0], p[2]);
1911     p[0] = t1;
1912     t1 = _mm_unpacklo_epi16(p[1], p[3]);
1913     p[3] = _mm_unpackhi_epi16(p[1], p[3]);
1914     p[1] = t1;
1915     t1 = _mm_unpacklo_epi16(p[4], p[6]);
1916     p[6] = _mm_unpackhi_epi16(p[4], p[6]);
1917     p[4] = t1;
1918     t1 = _mm_unpacklo_epi16(p[5], p[7]);
1919     p[7] = _mm_unpackhi_epi16(p[5], p[7]);
1920     p[5] = t1;
1921
1922 /*
1923     printf("After unpack pass 2:\n");
1924     for (i = 0; i < 8; i++) {
1925       MM_PRINT8("v", p[i]);
1926     }
1927  */
1928
1929     t1 = _mm_unpacklo_epi32(p[0], p[4]);
1930     p[4] = _mm_unpackhi_epi32(p[0], p[4]);
1931     p[0] = t1;
1932     t1 = _mm_unpacklo_epi32(p[1], p[5]);
1933     p[5] = _mm_unpackhi_epi32(p[1], p[5]);
1934     p[1] = t1;
1935     t1 = _mm_unpacklo_epi32(p[2], p[6]);
1936     p[6] = _mm_unpackhi_epi32(p[2], p[6]);
1937     p[2] = t1;
1938     t1 = _mm_unpacklo_epi32(p[3], p[7]);
1939     p[7] = _mm_unpackhi_epi32(p[3], p[7]);
1940     p[3] = t1;
1941
1942     if (xor) {
1943       for (i = 0; i < 8; i++) {
1944         t1 = _mm_load_si128((__m128i *) d64);
1945         _mm_store_si128((__m128i *) d64, _mm_xor_si128(p[i], t1));
1946         d64 += 2;
1947       }
1948     } else {
1949       for (i = 0; i < 8; i++) {
1950         _mm_store_si128((__m128i *) d64, p[i]);
1951         d64 += 2;
1952       }
1953     }
1954
1955   }
1956
1957   gf_do_final_region_alignment(&rd);
1958 }
1959 #endif
1960
1961 #define GF_MULTBY_TWO(p) (((p) & GF_FIRST_BIT) ? (((p) << 1) ^ h->prim_poly) : (p) << 1);
1962
1963 static
1964 int gf_w64_split_init(gf_t *gf)
1965 {
1966   gf_internal_t *h;
1967   struct gf_split_4_64_lazy_data *d4;
1968   struct gf_split_8_64_lazy_data *d8;
1969   struct gf_split_8_8_data *d88;
1970   struct gf_split_16_64_lazy_data *d16;
1971   uint64_t p, basep;
1972   int exp, i, j;
1973
1974   h = (gf_internal_t *) gf->scratch;
1975
1976   /* Defaults */
1977
1978   SET_FUNCTION(gf,multiply_region,w64,gf_w64_multiply_region_from_single)
1979
1980   SET_FUNCTION(gf,multiply,w64,gf_w64_bytwo_p_multiply)
1981
1982 #if defined(INTEL_SSE4_PCLMUL)
1983   if (gf_cpu_supports_intel_pclmul) {
1984     if ((!(h->region_type & GF_REGION_NOSIMD) &&
1985         (h->arg1 == 64 || h->arg2 == 64)) ||
1986         h->mult_type == GF_MULT_DEFAULT){
1987
1988       if ((0xfffffffe00000000ULL & h->prim_poly) == 0){
1989         SET_FUNCTION(gf,multiply,w64,gf_w64_clm_multiply_2)
1990         SET_FUNCTION(gf,multiply_region,w64,gf_w64_clm_multiply_region_from_single_2)
1991       }else if((0xfffe000000000000ULL & h->prim_poly) == 0){
1992         SET_FUNCTION(gf,multiply,w64,gf_w64_clm_multiply_4)
1993         SET_FUNCTION(gf,multiply_region,w64,gf_w64_clm_multiply_region_from_single_4)
1994       }else{
1995         return 0;
1996       }
1997     }
1998   }
1999 #endif
2000
2001   SET_FUNCTION(gf,inverse,w64,gf_w64_euclid)
2002
2003   /* Allen: set region pointers for default mult type. Single pointers are
2004    * taken care of above (explicitly for sse, implicitly for no sse). */
2005
2006   if (h->mult_type == GF_MULT_DEFAULT) {
2007 #if defined(INTEL_SSE4) || defined(ARCH_AARCH64)
2008     if (gf_cpu_supports_intel_sse4 || gf_cpu_supports_arm_neon) {
2009       d4 = (struct gf_split_4_64_lazy_data *) h->private;
2010       d4->last_value = 0;
2011 #if defined(INTEL_SSE4)
2012       if (gf_cpu_supports_intel_sse4)
2013         SET_FUNCTION(gf,multiply_region,w64,gf_w64_split_4_64_lazy_sse_multiply_region)
2014 #elif defined(ARCH_AARCH64)
2015       if (gf_cpu_supports_arm_neon)
2016         gf_w64_neon_split_init(gf);
2017 #endif
2018     } else {
2019 #endif
2020       d8 = (struct gf_split_8_64_lazy_data *) h->private;
2021       d8->last_value = 0;
2022       SET_FUNCTION(gf,multiply_region,w64,gf_w64_split_8_64_lazy_multiply_region)
2023 #if defined(INTEL_SSE4) || defined(ARCH_AARCH64)
2024     }
2025 #endif
2026   }
2027
2028   if ((h->arg1 == 4 && h->arg2 == 64) || (h->arg1 == 64 && h->arg2 == 4)) {
2029     d4 = (struct gf_split_4_64_lazy_data *) h->private;
2030     d4->last_value = 0;
2031
2032     if((h->region_type & GF_REGION_ALTMAP) && (h->region_type & GF_REGION_NOSIMD)) return 0;
2033     if(h->region_type & GF_REGION_ALTMAP)
2034     {
2035       #ifdef INTEL_SSSE3
2036         if (gf_cpu_supports_intel_ssse3) {
2037           SET_FUNCTION(gf,multiply_region,w64,gf_w64_split_4_64_lazy_sse_altmap_multiply_region)
2038         } else
2039       #elif defined(ARCH_AARCH64)
2040         if (gf_cpu_supports_arm_neon) {
2041           gf_w64_neon_split_init(gf);
2042         } else
2043       #endif
2044         return 0;
2045     }
2046     else //no altmap
2047     {
2048       #if defined(INTEL_SSE4) || defined(ARCH_AARCH64)
2049         if(gf_cpu_supports_intel_sse4 || gf_cpu_supports_arm_neon) {
2050           if (h->region_type & GF_REGION_NOSIMD) {
2051             SET_FUNCTION(gf,multiply_region,w64,gf_w64_split_4_64_lazy_multiply_region)
2052           } else
2053           #if defined(INTEL_SSE4)
2054             SET_FUNCTION(gf,multiply_region,w64,gf_w64_split_4_64_lazy_sse_multiply_region)
2055           #elif defined(ARCH_AARCH64)
2056             gf_w64_neon_split_init(gf);
2057           #endif
2058         } else {
2059       #endif
2060         SET_FUNCTION(gf,multiply_region,w64,gf_w64_split_4_64_lazy_multiply_region)
2061         if(h->region_type & GF_REGION_SIMD)
2062           return 0;
2063       #if defined(INTEL_SSE4) || defined(ARCH_AARCH64)
2064         }
2065       #endif
2066     }
2067   }
2068   if ((h->arg1 == 8 && h->arg2 == 64) || (h->arg1 == 64 && h->arg2 == 8)) {
2069     d8 = (struct gf_split_8_64_lazy_data *) h->private;
2070     d8->last_value = 0;
2071     SET_FUNCTION(gf,multiply_region,w64,gf_w64_split_8_64_lazy_multiply_region)
2072   }
2073   if ((h->arg1 == 16 && h->arg2 == 64) || (h->arg1 == 64 && h->arg2 == 16)) {
2074     d16 = (struct gf_split_16_64_lazy_data *) h->private;
2075     d16->last_value = 0;
2076     SET_FUNCTION(gf,multiply_region,w64,gf_w64_split_16_64_lazy_multiply_region)
2077   }
2078   if ((h->arg1 == 8 && h->arg2 == 8)) {
2079     d88 = (struct gf_split_8_8_data *) h->private;
2080     SET_FUNCTION(gf,multiply,w64,gf_w64_split_8_8_multiply)
2081
2082     /* The performance of this guy sucks, so don't bother with a region op */
2083
2084     basep = 1;
2085     for (exp = 0; exp < 15; exp++) {
2086       for (j = 0; j < 256; j++) d88->tables[exp][0][j] = 0;
2087       for (i = 0; i < 256; i++) d88->tables[exp][i][0] = 0;
2088       d88->tables[exp][1][1] = basep;
2089       for (i = 2; i < 256; i++) {
2090         if (i&1) {
2091           p = d88->tables[exp][i^1][1];
2092           d88->tables[exp][i][1] = p ^ basep;
2093         } else {
2094           p = d88->tables[exp][i>>1][1];
2095           d88->tables[exp][i][1] = GF_MULTBY_TWO(p);
2096         }
2097       }
2098       for (i = 1; i < 256; i++) {
2099         p = d88->tables[exp][i][1];
2100         for (j = 1; j < 256; j++) {
2101           if (j&1) {
2102             d88->tables[exp][i][j] = d88->tables[exp][i][j^1] ^ p;
2103           } else {
2104             d88->tables[exp][i][j] = GF_MULTBY_TWO(d88->tables[exp][i][j>>1]);
2105           }
2106         }
2107       }
2108       for (i = 0; i < 8; i++) basep = GF_MULTBY_TWO(basep);
2109     }
2110   }
2111   return 1;
2112 }
2113
2114 int gf_w64_scratch_size(int mult_type, int region_type, int divide_type, int arg1, int arg2)
2115 {
2116   switch(mult_type)
2117   {
2118     case GF_MULT_SHIFT:
2119       return sizeof(gf_internal_t);
2120       break;
2121     case GF_MULT_CARRY_FREE:
2122       return sizeof(gf_internal_t);
2123       break;
2124     case GF_MULT_BYTWO_p:
2125     case GF_MULT_BYTWO_b:
2126       return sizeof(gf_internal_t);
2127       break;
2128
2129     case GF_MULT_DEFAULT:
2130
2131       /* Allen: set the *local* arg1 and arg2, just for scratch size purposes,
2132        * then fall through to split table scratch size code. */
2133
2134 #if defined(INTEL_SSE4) || defined(ARCH_AARCH64)
2135     if (gf_cpu_supports_intel_sse4 || gf_cpu_supports_arm_neon) {
2136       arg1 = 64;
2137       arg2 = 4;
2138     } else {
2139 #endif
2140       arg1 = 64;
2141       arg2 = 8;
2142 #if defined(INTEL_SSE4) || defined(ARCH_AARCH64)
2143     }
2144 #endif
2145
2146     case GF_MULT_SPLIT_TABLE:
2147         if (arg1 == 8 && arg2 == 8) {
2148           return sizeof(gf_internal_t) + sizeof(struct gf_split_8_8_data) + 64;
2149         }
2150         if ((arg1 == 16 && arg2 == 64) || (arg2 == 16 && arg1 == 64)) {
2151           return sizeof(gf_internal_t) + sizeof(struct gf_split_16_64_lazy_data) + 64;
2152         }
2153         if ((arg1 == 8 && arg2 == 64) || (arg2 == 8 && arg1 == 64)) {
2154           return sizeof(gf_internal_t) + sizeof(struct gf_split_8_64_lazy_data) + 64;
2155         }
2156
2157         if ((arg1 == 64 && arg2 == 4) || (arg1 == 4 && arg2 == 64)) {
2158           return sizeof(gf_internal_t) + sizeof(struct gf_split_4_64_lazy_data) + 64;
2159         }
2160         return 0;
2161     case GF_MULT_GROUP:
2162       return sizeof(gf_internal_t) + sizeof(struct gf_w64_group_data) +
2163                sizeof(uint64_t) * (1 << arg1) +
2164                sizeof(uint64_t) * (1 << arg2) + 64;
2165       break;
2166     case GF_MULT_COMPOSITE:
2167       if (arg1 == 2) return sizeof(gf_internal_t) + 64;
2168       return 0;
2169       break;
2170     default:
2171       return 0;
2172    }
2173 }
2174
2175 int gf_w64_init(gf_t *gf)
2176 {
2177   gf_internal_t *h;
2178
2179   h = (gf_internal_t *) gf->scratch;
2180
2181   /* Allen: set default primitive polynomial / irreducible polynomial if needed */
2182
2183   /* Omitting the leftmost 1 as in w=32 */
2184
2185   if (h->prim_poly == 0) {
2186     if (h->mult_type == GF_MULT_COMPOSITE) {
2187       h->prim_poly = gf_composite_get_default_poly(h->base_gf);
2188       if (h->prim_poly == 0) return 0; /* This shouldn't happen */
2189     } else {
2190       h->prim_poly = 0x1b;
2191     }
2192   }
2193
2194   SET_FUNCTION(gf,multiply,w64,NULL)
2195   SET_FUNCTION(gf,divide,w64,NULL)
2196   SET_FUNCTION(gf,inverse,w64,NULL)
2197   SET_FUNCTION(gf,multiply_region,w64,NULL)
2198
2199   switch(h->mult_type) {
2200     case GF_MULT_CARRY_FREE:  if (gf_w64_cfm_init(gf) == 0) return 0; break;
2201     case GF_MULT_SHIFT:       if (gf_w64_shift_init(gf) == 0) return 0; break;
2202     case GF_MULT_COMPOSITE:   if (gf_w64_composite_init(gf) == 0) return 0; break;
2203     case GF_MULT_DEFAULT:
2204     case GF_MULT_SPLIT_TABLE: if (gf_w64_split_init(gf) == 0) return 0; break;
2205     case GF_MULT_GROUP:       if (gf_w64_group_init(gf) == 0) return 0; break;
2206     case GF_MULT_BYTWO_p:
2207     case GF_MULT_BYTWO_b:     if (gf_w64_bytwo_init(gf) == 0) return 0; break;
2208     default: return 0;
2209   }
2210   if (h->divide_type == GF_DIVIDE_EUCLID) {
2211     SET_FUNCTION(gf,divide,w64,gf_w64_divide_from_inverse)
2212     SET_FUNCTION(gf,inverse,w64,gf_w64_euclid)
2213   }
2214
2215   if (gf->inverse.w64 != NULL && gf->divide.w64 == NULL) {
2216     SET_FUNCTION(gf,divide,w64,gf_w64_divide_from_inverse)
2217   }
2218   if (gf->inverse.w64 == NULL && gf->divide.w64 != NULL) {
2219     SET_FUNCTION(gf,inverse,w64,gf_w64_inverse_from_divide)
2220   }
2221
2222   if (h->region_type == GF_REGION_CAUCHY) return 0;
2223
2224   if (h->region_type & GF_REGION_ALTMAP) {
2225     if (h->mult_type == GF_MULT_COMPOSITE) {
2226       SET_FUNCTION(gf,extract_word,w64,gf_w64_composite_extract_word)
2227     } else if (h->mult_type == GF_MULT_SPLIT_TABLE) {
2228       SET_FUNCTION(gf,extract_word,w64,gf_w64_split_extract_word)
2229     }
2230   } else {
2231     SET_FUNCTION(gf,extract_word,w64,gf_w64_extract_word)
2232   }
2233
2234   return 1;
2235 }