IDA_new/gf-complete/src/gf_w16.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_w16.c
   7  *
   8  * Routines for 16-bit Galois fields
   9  */
  10
  11 #include "gf_int.h"
  12 #include <stdio.h>
  13 #include <stdlib.h>
  14 #include "gf_w16.h"
  15 #include "gf_cpu.h"
  16
  17 #define AB2(ip, am1 ,am2, b, t1, t2) {\
  18   t1 = (b << 1) & am1;\
  19   t2 = b & am2; \
  20   t2 = ((t2 << 1) - (t2 >> (GF_FIELD_WIDTH-1))); \
  21   b = (t1 ^ (t2 & ip));}
  22
  23 #define SSE_AB2(pp, m1 ,m2, va, t1, t2) {\
  24           t1 = _mm_and_si128(_mm_slli_epi64(va, 1), m1); \
  25           t2 = _mm_and_si128(va, m2); \
  26           t2 = _mm_sub_epi64 (_mm_slli_epi64(t2, 1), _mm_srli_epi64(t2, (GF_FIELD_WIDTH-1))); \
  27           va = _mm_xor_si128(t1, _mm_and_si128(t2, pp)); }
  28
  29 #define MM_PRINT(s, r) { uint8_t blah[16], ii; printf("%-12s", s); _mm_storeu_si128((__m128i *)blah, r); for (ii = 0; ii < 16; ii += 2) printf("  %02x %02x", blah[15-ii], blah[14-ii]); printf("\n"); }
  30
  31 #define GF_FIRST_BIT (1 << 15)
  32 #define GF_MULTBY_TWO(p) (((p) & GF_FIRST_BIT) ? (((p) << 1) ^ h->prim_poly) : (p) << 1)
  33
  34 static
  35 inline
  36 gf_val_32_t gf_w16_inverse_from_divide (gf_t *gf, gf_val_32_t a)
  37 {
  38   return gf->divide.w32(gf, 1, a);
  39 }
  40
  41 static
  42 inline
  43 gf_val_32_t gf_w16_divide_from_inverse (gf_t *gf, gf_val_32_t a, gf_val_32_t b)
  44 {
  45   b = gf->inverse.w32(gf, b);
  46   return gf->multiply.w32(gf, a, b);
  47 }
  48
  49 static
  50 void
  51 gf_w16_multiply_region_from_single(gf_t *gf, void *src, void *dest, gf_val_32_t val, int bytes, int xor)
  52 {
  53   gf_region_data rd;
  54   uint16_t *s16;
  55   uint16_t *d16;
  56
  57   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
  58   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
  59
  60   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 2);
  61   gf_do_initial_region_alignment(&rd);
  62
  63   s16 = (uint16_t *) rd.s_start;
  64   d16 = (uint16_t *) rd.d_start;
  65
  66   if (xor) {
  67     while (d16 < ((uint16_t *) rd.d_top)) {
  68       *d16 ^= gf->multiply.w32(gf, val, *s16);
  69       d16++;
  70       s16++;
  71     }
  72   } else {
  73     while (d16 < ((uint16_t *) rd.d_top)) {
  74       *d16 = gf->multiply.w32(gf, val, *s16);
  75       d16++;
  76       s16++;
  77     }
  78   }
  79   gf_do_final_region_alignment(&rd);
  80 }
  81
  82 #if defined(INTEL_SSE4_PCLMUL)
  83 static
  84 void
  85 gf_w16_clm_multiply_region_from_single_2(gf_t *gf, void *src, void *dest, gf_val_32_t val, int bytes, int xor)
  86 {
  87   gf_region_data rd;
  88   uint16_t *s16;
  89   uint16_t *d16;
  90   __m128i         a, b;
  91   __m128i         result;
  92   __m128i         prim_poly;
  93   __m128i         w;
  94   gf_internal_t * h = gf->scratch;
  95   prim_poly = _mm_set_epi32(0, 0, 0, (uint32_t)(h->prim_poly & 0x1ffffULL));
  96
  97   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
  98   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
  99
 100   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 2);
 101   gf_do_initial_region_alignment(&rd);
 102
 103   a = _mm_insert_epi32 (_mm_setzero_si128(), val, 0);
 104
 105   s16 = (uint16_t *) rd.s_start;
 106   d16 = (uint16_t *) rd.d_start;
 107
 108   if (xor) {
 109     while (d16 < ((uint16_t *) rd.d_top)) {
 110
 111       /* see gf_w16_clm_multiply() to see explanation of method */
 112
 113       b = _mm_insert_epi32 (a, (gf_val_32_t)(*s16), 0);
 114       result = _mm_clmulepi64_si128 (a, b, 0);
 115       w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 116       result = _mm_xor_si128 (result, w);
 117       w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 118       result = _mm_xor_si128 (result, w);
 119
 120       *d16 ^= ((gf_val_32_t)_mm_extract_epi32(result, 0));
 121       d16++;
 122       s16++;
 123     }
 124   } else {
 125     while (d16 < ((uint16_t *) rd.d_top)) {
 126
 127       /* see gf_w16_clm_multiply() to see explanation of method */
 128
 129       b = _mm_insert_epi32 (a, (gf_val_32_t)(*s16), 0);
 130       result = _mm_clmulepi64_si128 (a, b, 0);
 131       w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 132       result = _mm_xor_si128 (result, w);
 133       w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 134       result = _mm_xor_si128 (result, w);
 135
 136       *d16 = ((gf_val_32_t)_mm_extract_epi32(result, 0));
 137       d16++;
 138       s16++;
 139     }
 140   }
 141   gf_do_final_region_alignment(&rd);
 142 }
 143 #endif
 144
 145 #if defined(INTEL_SSE4_PCLMUL)
 146 static
 147 void
 148 gf_w16_clm_multiply_region_from_single_3(gf_t *gf, void *src, void *dest, gf_val_32_t val, int bytes, int xor)
 149 {
 150   gf_region_data rd;
 151   uint16_t *s16;
 152   uint16_t *d16;
 153
 154   __m128i         a, b;
 155   __m128i         result;
 156   __m128i         prim_poly;
 157   __m128i         w;
 158   gf_internal_t * h = gf->scratch;
 159   prim_poly = _mm_set_epi32(0, 0, 0, (uint32_t)(h->prim_poly & 0x1ffffULL));
 160
 161   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
 162   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
 163
 164   a = _mm_insert_epi32 (_mm_setzero_si128(), val, 0);
 165
 166   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 2);
 167   gf_do_initial_region_alignment(&rd);
 168
 169   s16 = (uint16_t *) rd.s_start;
 170   d16 = (uint16_t *) rd.d_start;
 171
 172   if (xor) {
 173     while (d16 < ((uint16_t *) rd.d_top)) {
 174
 175       /* see gf_w16_clm_multiply() to see explanation of method */
 176
 177       b = _mm_insert_epi32 (a, (gf_val_32_t)(*s16), 0);
 178       result = _mm_clmulepi64_si128 (a, b, 0);
 179       w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 180       result = _mm_xor_si128 (result, w);
 181       w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 182       result = _mm_xor_si128 (result, w);
 183       w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 184       result = _mm_xor_si128 (result, w);
 185
 186       *d16 ^= ((gf_val_32_t)_mm_extract_epi32(result, 0));
 187       d16++;
 188       s16++;
 189     }
 190   } else {
 191     while (d16 < ((uint16_t *) rd.d_top)) {
 192
 193       /* see gf_w16_clm_multiply() to see explanation of method */
 194
 195       b = _mm_insert_epi32 (a, (gf_val_32_t)(*s16), 0);
 196       result = _mm_clmulepi64_si128 (a, b, 0);
 197       w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 198       result = _mm_xor_si128 (result, w);
 199       w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 200       result = _mm_xor_si128 (result, w);
 201       w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 202       result = _mm_xor_si128 (result, w);
 203
 204       *d16 = ((gf_val_32_t)_mm_extract_epi32(result, 0));
 205       d16++;
 206       s16++;
 207     }
 208   }
 209   gf_do_final_region_alignment(&rd);
 210 }
 211 #endif
 212
 213 #if defined(INTEL_SSE4_PCLMUL)
 214 static
 215 void
 216 gf_w16_clm_multiply_region_from_single_4(gf_t *gf, void *src, void *dest, gf_val_32_t val, int bytes, int xor)
 217 {
 218   gf_region_data rd;
 219   uint16_t *s16;
 220   uint16_t *d16;
 221
 222   __m128i         a, b;
 223   __m128i         result;
 224   __m128i         prim_poly;
 225   __m128i         w;
 226   gf_internal_t * h = gf->scratch;
 227   prim_poly = _mm_set_epi32(0, 0, 0, (uint32_t)(h->prim_poly & 0x1ffffULL));
 228
 229   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
 230   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
 231
 232   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 2);
 233   gf_do_initial_region_alignment(&rd);
 234
 235   a = _mm_insert_epi32 (_mm_setzero_si128(), val, 0);
 236
 237   s16 = (uint16_t *) rd.s_start;
 238   d16 = (uint16_t *) rd.d_start;
 239
 240   if (xor) {
 241     while (d16 < ((uint16_t *) rd.d_top)) {
 242
 243       /* see gf_w16_clm_multiply() to see explanation of method */
 244
 245       b = _mm_insert_epi32 (a, (gf_val_32_t)(*s16), 0);
 246       result = _mm_clmulepi64_si128 (a, b, 0);
 247       w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 248       result = _mm_xor_si128 (result, w);
 249       w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 250       result = _mm_xor_si128 (result, w);
 251       w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 252       result = _mm_xor_si128 (result, w);
 253       w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 254       result = _mm_xor_si128 (result, w);
 255
 256       *d16 ^= ((gf_val_32_t)_mm_extract_epi32(result, 0));
 257       d16++;
 258       s16++;
 259     }
 260   } else {
 261     while (d16 < ((uint16_t *) rd.d_top)) {
 262
 263       /* see gf_w16_clm_multiply() to see explanation of method */
 264
 265       b = _mm_insert_epi32 (a, (gf_val_32_t)(*s16), 0);
 266       result = _mm_clmulepi64_si128 (a, b, 0);
 267       w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 268       result = _mm_xor_si128 (result, w);
 269       w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 270       result = _mm_xor_si128 (result, w);
 271       w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 272       result = _mm_xor_si128 (result, w);
 273       w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 274       result = _mm_xor_si128 (result, w);
 275
 276       *d16 = ((gf_val_32_t)_mm_extract_epi32(result, 0));
 277       d16++;
 278       s16++;
 279     }
 280   }
 281   gf_do_final_region_alignment(&rd);
 282 }
 283 #endif
 284
 285 static
 286 inline
 287 gf_val_32_t gf_w16_euclid (gf_t *gf, gf_val_32_t b)
 288 {
 289   gf_val_32_t e_i, e_im1, e_ip1;
 290   gf_val_32_t d_i, d_im1, d_ip1;
 291   gf_val_32_t y_i, y_im1, y_ip1;
 292   gf_val_32_t c_i;
 293
 294   if (b == 0) return -1;
 295   e_im1 = ((gf_internal_t *) (gf->scratch))->prim_poly;
 296   e_i = b;
 297   d_im1 = 16;
 298   for (d_i = d_im1; ((1 << d_i) & e_i) == 0; d_i--) ;
 299   y_i = 1;
 300   y_im1 = 0;
 301
 302   while (e_i != 1) {
 303
 304     e_ip1 = e_im1;
 305     d_ip1 = d_im1;
 306     c_i = 0;
 307
 308     while (d_ip1 >= d_i) {
 309       c_i ^= (1 << (d_ip1 - d_i));
 310       e_ip1 ^= (e_i << (d_ip1 - d_i));
 311       if (e_ip1 == 0) return 0;
 312       while ((e_ip1 & (1 << d_ip1)) == 0) d_ip1--;
 313     }
 314
 315     y_ip1 = y_im1 ^ gf->multiply.w32(gf, c_i, y_i);
 316     y_im1 = y_i;
 317     y_i = y_ip1;
 318
 319     e_im1 = e_i;
 320     d_im1 = d_i;
 321     e_i = e_ip1;
 322     d_i = d_ip1;
 323   }
 324
 325   return y_i;
 326 }
 327
 328 static
 329 gf_val_32_t gf_w16_extract_word(gf_t *gf, void *start, int bytes, int index)
 330 {
 331   uint16_t *r16, rv;
 332
 333   r16 = (uint16_t *) start;
 334   rv = r16[index];
 335   return rv;
 336 }
 337
 338 static
 339 gf_val_32_t gf_w16_composite_extract_word(gf_t *gf, void *start, int bytes, int index)
 340 {
 341   int sub_size;
 342   gf_internal_t *h;
 343   uint8_t *r8, *top;
 344   uint16_t a, b, *r16;
 345   gf_region_data rd;
 346
 347   h = (gf_internal_t *) gf->scratch;
 348   gf_set_region_data(&rd, gf, start, start, bytes, 0, 0, 32);
 349   r16 = (uint16_t *) start;
 350   if (r16 + index < (uint16_t *) rd.d_start) return r16[index];
 351   if (r16 + index >= (uint16_t *) rd.d_top) return r16[index];
 352   index -= (((uint16_t *) rd.d_start) - r16);
 353   r8 = (uint8_t *) rd.d_start;
 354   top = (uint8_t *) rd.d_top;
 355   sub_size = (top-r8)/2;
 356
 357   a = h->base_gf->extract_word.w32(h->base_gf, r8, sub_size, index);
 358   b = h->base_gf->extract_word.w32(h->base_gf, r8+sub_size, sub_size, index);
 359   return (a | (b << 8));
 360 }
 361
 362 static
 363 gf_val_32_t gf_w16_split_extract_word(gf_t *gf, void *start, int bytes, int index)
 364 {
 365   uint16_t *r16, rv;
 366   uint8_t *r8;
 367   gf_region_data rd;
 368
 369   gf_set_region_data(&rd, gf, start, start, bytes, 0, 0, 32);
 370   r16 = (uint16_t *) start;
 371   if (r16 + index < (uint16_t *) rd.d_start) return r16[index];
 372   if (r16 + index >= (uint16_t *) rd.d_top) return r16[index];
 373   index -= (((uint16_t *) rd.d_start) - r16);
 374   r8 = (uint8_t *) rd.d_start;
 375   r8 += ((index & 0xfffffff0)*2);
 376   r8 += (index & 0xf);
 377   rv = (*r8 << 8);
 378   r8 += 16;
 379   rv |= *r8;
 380   return rv;
 381 }
 382
 383 static
 384 inline
 385 gf_val_32_t gf_w16_matrix (gf_t *gf, gf_val_32_t b)
 386 {
 387   return gf_bitmatrix_inverse(b, 16, ((gf_internal_t *) (gf->scratch))->prim_poly);
 388 }
 389
 390 /* JSP: GF_MULT_SHIFT: The world's dumbest multiplication algorithm.  I only
 391    include it for completeness.  It does have the feature that it requires no
 392    extra memory.
 393  */
 394
 395 #if defined(INTEL_SSE4_PCLMUL)
 396 static
 397 inline
 398 gf_val_32_t
 399 gf_w16_clm_multiply_2 (gf_t *gf, gf_val_32_t a16, gf_val_32_t b16)
 400 {
 401   gf_val_32_t rv = 0;
 402
 403   __m128i         a, b;
 404   __m128i         result;
 405   __m128i         prim_poly;
 406   __m128i         w;
 407   gf_internal_t * h = gf->scratch;
 408
 409   a = _mm_insert_epi32 (_mm_setzero_si128(), a16, 0);
 410   b = _mm_insert_epi32 (a, b16, 0);
 411
 412   prim_poly = _mm_set_epi32(0, 0, 0, (uint32_t)(h->prim_poly & 0x1ffffULL));
 413
 414   /* Do the initial multiply */
 415
 416   result = _mm_clmulepi64_si128 (a, b, 0);
 417
 418   /* Ben: Do prim_poly reduction twice. We are guaranteed that we will only
 419      have to do the reduction at most twice, because (w-2)/z == 2. Where
 420      z is equal to the number of zeros after the leading 1
 421
 422      _mm_clmulepi64_si128 is the carryless multiply operation. Here
 423      _mm_srli_si128 shifts the result to the right by 2 bytes. This allows
 424      us to multiply the prim_poly by the leading bits of the result. We
 425      then xor the result of that operation back with the result.*/
 426
 427   w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 428   result = _mm_xor_si128 (result, w);
 429   w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 430   result = _mm_xor_si128 (result, w);
 431
 432   /* Extracts 32 bit value from result. */
 433
 434   rv = ((gf_val_32_t)_mm_extract_epi32(result, 0));
 435
 436   return rv;
 437 }
 438 #endif
 439
 440 #if defined(INTEL_SSE4_PCLMUL)
 441 static
 442 inline
 443 gf_val_32_t
 444 gf_w16_clm_multiply_3 (gf_t *gf, gf_val_32_t a16, gf_val_32_t b16)
 445 {
 446   gf_val_32_t rv = 0;
 447
 448   __m128i         a, b;
 449   __m128i         result;
 450   __m128i         prim_poly;
 451   __m128i         w;
 452   gf_internal_t * h = gf->scratch;
 453
 454   a = _mm_insert_epi32 (_mm_setzero_si128(), a16, 0);
 455   b = _mm_insert_epi32 (a, b16, 0);
 456
 457   prim_poly = _mm_set_epi32(0, 0, 0, (uint32_t)(h->prim_poly & 0x1ffffULL));
 458
 459   /* Do the initial multiply */
 460
 461   result = _mm_clmulepi64_si128 (a, b, 0);
 462
 463   w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 464   result = _mm_xor_si128 (result, w);
 465   w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 466   result = _mm_xor_si128 (result, w);
 467   w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 468   result = _mm_xor_si128 (result, w);
 469
 470   /* Extracts 32 bit value from result. */
 471
 472   rv = ((gf_val_32_t)_mm_extract_epi32(result, 0));
 473
 474   return rv;
 475 }
 476 #endif
 477
 478 #if defined(INTEL_SSE4_PCLMUL)
 479 static
 480 inline
 481 gf_val_32_t
 482 gf_w16_clm_multiply_4 (gf_t *gf, gf_val_32_t a16, gf_val_32_t b16)
 483 {
 484   gf_val_32_t rv = 0;
 485
 486   __m128i         a, b;
 487   __m128i         result;
 488   __m128i         prim_poly;
 489   __m128i         w;
 490   gf_internal_t * h = gf->scratch;
 491
 492   a = _mm_insert_epi32 (_mm_setzero_si128(), a16, 0);
 493   b = _mm_insert_epi32 (a, b16, 0);
 494
 495   prim_poly = _mm_set_epi32(0, 0, 0, (uint32_t)(h->prim_poly & 0x1ffffULL));
 496
 497   /* Do the initial multiply */
 498
 499   result = _mm_clmulepi64_si128 (a, b, 0);
 500
 501   w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 502   result = _mm_xor_si128 (result, w);
 503   w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 504   result = _mm_xor_si128 (result, w);
 505   w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 506   result = _mm_xor_si128 (result, w);
 507   w = _mm_clmulepi64_si128 (prim_poly, _mm_srli_si128 (result, 2), 0);
 508   result = _mm_xor_si128 (result, w);
 509
 510   /* Extracts 32 bit value from result. */
 511
 512   rv = ((gf_val_32_t)_mm_extract_epi32(result, 0));
 513
 514   return rv;
 515 }
 516 #endif
 517
 518
 519 static
 520 inline
 521  gf_val_32_t
 522 gf_w16_shift_multiply (gf_t *gf, gf_val_32_t a16, gf_val_32_t b16)
 523 {
 524   gf_val_32_t product, i, pp, a, b;
 525   gf_internal_t *h;
 526
 527   a = a16;
 528   b = b16;
 529   h = (gf_internal_t *) gf->scratch;
 530   pp = h->prim_poly;
 531
 532   product = 0;
 533
 534   for (i = 0; i < GF_FIELD_WIDTH; i++) {
 535     if (a & (1 << i)) product ^= (b << i);
 536   }
 537   for (i = (GF_FIELD_WIDTH*2-2); i >= GF_FIELD_WIDTH; i--) {
 538     if (product & (1 << i)) product ^= (pp << (i-GF_FIELD_WIDTH));
 539   }
 540   return product;
 541 }
 542
 543 static
 544 int gf_w16_shift_init(gf_t *gf)
 545 {
 546   SET_FUNCTION(gf,multiply,w32,gf_w16_shift_multiply)
 547   return 1;
 548 }
 549
 550 static
 551 int gf_w16_cfm_init(gf_t *gf)
 552 {
 553 #if defined(INTEL_SSE4_PCLMUL)
 554   if (gf_cpu_supports_intel_pclmul) {
 555     gf_internal_t *h;
 556
 557     h = (gf_internal_t *) gf->scratch;
 558
 559     /*Ben: Determining how many reductions to do */
 560
 561     if ((0xfe00 & h->prim_poly) == 0) {
 562       SET_FUNCTION(gf,multiply,w32,gf_w16_clm_multiply_2)
 563       SET_FUNCTION(gf,multiply_region,w32,gf_w16_clm_multiply_region_from_single_2)
 564     } else if((0xf000 & h->prim_poly) == 0) {
 565       SET_FUNCTION(gf,multiply,w32,gf_w16_clm_multiply_3)
 566       SET_FUNCTION(gf,multiply_region,w32,gf_w16_clm_multiply_region_from_single_3)
 567     } else if ((0xe000 & h->prim_poly) == 0) {
 568       SET_FUNCTION(gf,multiply,w32,gf_w16_clm_multiply_4)
 569       SET_FUNCTION(gf,multiply_region,w32,gf_w16_clm_multiply_region_from_single_4)
 570     } else {
 571       return 0;
 572     }
 573     return 1;
 574   }
 575 #endif
 576
 577   return 0;
 578 }
 579
 580 /* KMG: GF_MULT_LOGTABLE: */
 581
 582 static
 583 void
 584 gf_w16_log_multiply_region(gf_t *gf, void *src, void *dest, gf_val_32_t val, int bytes, int xor)
 585 {
 586   uint16_t *s16, *d16;
 587   int lv;
 588   struct gf_w16_logtable_data *ltd;
 589   gf_region_data rd;
 590
 591   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
 592   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
 593
 594   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 2);
 595   gf_do_initial_region_alignment(&rd);
 596
 597   ltd = (struct gf_w16_logtable_data *) ((gf_internal_t *) gf->scratch)->private;
 598   s16 = (uint16_t *) rd.s_start;
 599   d16 = (uint16_t *) rd.d_start;
 600
 601   lv = ltd->log_tbl[val];
 602
 603   if (xor) {
 604     while (d16 < (uint16_t *) rd.d_top) {
 605       *d16 ^= (*s16 == 0 ? 0 : ltd->antilog_tbl[lv + ltd->log_tbl[*s16]]);
 606       d16++;
 607       s16++;
 608     }
 609   } else {
 610     while (d16 < (uint16_t *) rd.d_top) {
 611       *d16 = (*s16 == 0 ? 0 : ltd->antilog_tbl[lv + ltd->log_tbl[*s16]]);
 612       d16++;
 613       s16++;
 614     }
 615   }
 616   gf_do_final_region_alignment(&rd);
 617 }
 618
 619 static
 620 inline
 621 gf_val_32_t
 622 gf_w16_log_multiply(gf_t *gf, gf_val_32_t a, gf_val_32_t b)
 623 {
 624   struct gf_w16_logtable_data *ltd;
 625
 626   ltd = (struct gf_w16_logtable_data *) ((gf_internal_t *) gf->scratch)->private;
 627   return (a == 0 || b == 0) ? 0 : ltd->antilog_tbl[(int) ltd->log_tbl[a] + (int) ltd->log_tbl[b]];
 628 }
 629
 630 static
 631 inline
 632 gf_val_32_t
 633 gf_w16_log_divide(gf_t *gf, gf_val_32_t a, gf_val_32_t b)
 634 {
 635   int log_sum = 0;
 636   struct gf_w16_logtable_data *ltd;
 637
 638   if (a == 0 || b == 0) return 0;
 639   ltd = (struct gf_w16_logtable_data *) ((gf_internal_t *) gf->scratch)->private;
 640
 641   log_sum = (int) ltd->log_tbl[a] - (int) ltd->log_tbl[b];
 642   return (ltd->d_antilog[log_sum]);
 643 }
 644
 645 static
 646 gf_val_32_t
 647 gf_w16_log_inverse(gf_t *gf, gf_val_32_t a)
 648 {
 649   struct gf_w16_logtable_data *ltd;
 650
 651   ltd = (struct gf_w16_logtable_data *) ((gf_internal_t *) gf->scratch)->private;
 652   return (ltd->inv_tbl[a]);
 653 }
 654
 655 static
 656 int gf_w16_log_init(gf_t *gf)
 657 {
 658   gf_internal_t *h;
 659   struct gf_w16_logtable_data *ltd;
 660   int i, b;
 661   int check = 0;
 662
 663   h = (gf_internal_t *) gf->scratch;
 664   ltd = h->private;
 665
 666   for (i = 0; i < GF_MULT_GROUP_SIZE+1; i++)
 667     ltd->log_tbl[i] = 0;
 668   ltd->d_antilog = ltd->antilog_tbl + GF_MULT_GROUP_SIZE;
 669
 670   b = 1;
 671   for (i = 0; i < GF_MULT_GROUP_SIZE; i++) {
 672       if (ltd->log_tbl[b] != 0) check = 1;
 673       ltd->log_tbl[b] = i;
 674       ltd->antilog_tbl[i] = b;
 675       ltd->antilog_tbl[i+GF_MULT_GROUP_SIZE] = b;
 676       b <<= 1;
 677       if (b & GF_FIELD_SIZE) {
 678           b = b ^ h->prim_poly;
 679       }
 680   }
 681
 682   /* If you can't construct the log table, there's a problem.  This code is used for
 683      some other implementations (e.g. in SPLIT), so if the log table doesn't work in
 684      that instance, use CARRY_FREE / SHIFT instead. */
 685
 686   if (check) {
 687     if (h->mult_type != GF_MULT_LOG_TABLE) {
 688       if (gf_cpu_supports_intel_pclmul) {
 689         return gf_w16_cfm_init(gf);
 690       }
 691       return gf_w16_shift_init(gf);
 692     } else {
 693       _gf_errno = GF_E_LOGPOLY;
 694       return 0;
 695     }
 696   }
 697
 698   ltd->inv_tbl[0] = 0;  /* Not really, but we need to fill it with something  */
 699   ltd->inv_tbl[1] = 1;
 700   for (i = 2; i < GF_FIELD_SIZE; i++) {
 701     ltd->inv_tbl[i] = ltd->antilog_tbl[GF_MULT_GROUP_SIZE-ltd->log_tbl[i]];
 702   }
 703
 704   SET_FUNCTION(gf,inverse,w32,gf_w16_log_inverse)
 705   SET_FUNCTION(gf,divide,w32,gf_w16_log_divide)
 706   SET_FUNCTION(gf,multiply,w32,gf_w16_log_multiply)
 707   SET_FUNCTION(gf,multiply_region,w32,gf_w16_log_multiply_region)
 708
 709   return 1;
 710 }
 711
 712 /* JSP: GF_MULT_SPLIT_TABLE: Using 8 multiplication tables to leverage SSE instructions.
 713 */
 714
 715
 716 /* Ben: Does alternate mapping multiplication using a split table in the
 717  lazy method without sse instructions*/
 718
 719 static
 720 void
 721 gf_w16_split_4_16_lazy_nosse_altmap_multiply_region(gf_t *gf, void *src, void *dest, gf_val_32_t val, int bytes, int xor)
 722 {
 723   uint64_t i, j, c, prod;
 724   uint8_t *s8, *d8, *top;
 725   uint16_t table[4][16];
 726   gf_region_data rd;
 727
 728   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
 729   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
 730
 731   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 32);
 732   gf_do_initial_region_alignment(&rd);
 733
 734   /*Ben: Constructs lazy multiplication table*/
 735
 736   for (j = 0; j < 16; j++) {
 737     for (i = 0; i < 4; i++) {
 738       c = (j << (i*4));
 739       table[i][j] = gf->multiply.w32(gf, c, val);
 740     }
 741   }
 742
 743   /*Ben: s8 is the start of source, d8 is the start of dest, top is end of dest region. */
 744
 745   s8 = (uint8_t *) rd.s_start;
 746   d8 = (uint8_t *) rd.d_start;
 747   top = (uint8_t *) rd.d_top;
 748
 749
 750   while (d8 < top) {
 751
 752     /*Ben: Multiplies across 16 two byte quantities using alternate mapping
 753        high bits are on the left, low bits are on the right. */
 754
 755     for (j=0;j<16;j++) {
 756
 757       /*Ben: If the xor flag is set, the product should include what is in dest */
 758       prod = (xor) ? ((uint16_t)(*d8)<<8) ^ *(d8+16) : 0;
 759
 760       /*Ben: xors all 4 table lookups into the product variable*/
 761
 762       prod ^= ((table[0][*(s8+16)&0xf]) ^
 763           (table[1][(*(s8+16)&0xf0)>>4]) ^
 764           (table[2][*(s8)&0xf]) ^
 765           (table[3][(*(s8)&0xf0)>>4]));
 766
 767       /*Ben: Stores product in the destination and moves on*/
 768
 769       *d8 = (uint8_t)(prod >> 8);
 770       *(d8+16) = (uint8_t)(prod & 0x00ff);
 771       s8++;
 772       d8++;
 773     }
 774     s8+=16;
 775     d8+=16;
 776   }
 777   gf_do_final_region_alignment(&rd);
 778 }
 779
 780 static
 781   void
 782 gf_w16_split_4_16_lazy_multiply_region(gf_t *gf, void *src, void *dest, gf_val_32_t val, int bytes, int xor)
 783 {
 784   uint64_t i, j, a, c, prod;
 785   uint16_t *s16, *d16, *top;
 786   uint16_t table[4][16];
 787   gf_region_data rd;
 788
 789   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
 790   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
 791
 792   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 2);
 793   gf_do_initial_region_alignment(&rd);
 794
 795   for (j = 0; j < 16; j++) {
 796     for (i = 0; i < 4; i++) {
 797       c = (j << (i*4));
 798       table[i][j] = gf->multiply.w32(gf, c, val);
 799     }
 800   }
 801
 802   s16 = (uint16_t *) rd.s_start;
 803   d16 = (uint16_t *) rd.d_start;
 804   top = (uint16_t *) rd.d_top;
 805
 806   while (d16 < top) {
 807     a = *s16;
 808     prod = (xor) ? *d16 : 0;
 809     for (i = 0; i < 4; i++) {
 810       prod ^= table[i][a&0xf];
 811       a >>= 4;
 812     }
 813     *d16 = prod;
 814     s16++;
 815     d16++;
 816   }
 817 }
 818
 819 static
 820 void
 821 gf_w16_split_8_16_lazy_multiply_region(gf_t *gf, void *src, void *dest, gf_val_32_t val, int bytes, int xor)
 822 {
 823   uint64_t j, k, v, a, prod, *s64, *d64, *top64;
 824   gf_internal_t *h;
 825   uint64_t htable[256], ltable[256];
 826   gf_region_data rd;
 827
 828   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
 829   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
 830
 831   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 8);
 832   gf_do_initial_region_alignment(&rd);
 833
 834   h = (gf_internal_t *) gf->scratch;
 835
 836   v = val;
 837   ltable[0] = 0;
 838   for (j = 1; j < 256; j <<= 1) {
 839     for (k = 0; k < j; k++) ltable[k^j] = (v ^ ltable[k]);
 840     v = GF_MULTBY_TWO(v);
 841   }
 842   htable[0] = 0;
 843   for (j = 1; j < 256; j <<= 1) {
 844     for (k = 0; k < j; k++) htable[k^j] = (v ^ htable[k]);
 845     v = GF_MULTBY_TWO(v);
 846   }
 847
 848   s64 = (uint64_t *) rd.s_start;
 849   d64 = (uint64_t *) rd.d_start;
 850   top64 = (uint64_t *) rd.d_top;
 851
 852 /* Does Unrolling Matter?  -- Doesn't seem to.
 853   while (d64 != top64) {
 854     a = *s64;
 855
 856     prod = htable[a >> 56];
 857     a <<= 8;
 858     prod ^= ltable[a >> 56];
 859     a <<= 8;
 860     prod <<= 16;
 861
 862     prod ^= htable[a >> 56];
 863     a <<= 8;
 864     prod ^= ltable[a >> 56];
 865     a <<= 8;
 866     prod <<= 16;
 867
 868     prod ^= htable[a >> 56];
 869     a <<= 8;
 870     prod ^= ltable[a >> 56];
 871     a <<= 8;
 872     prod <<= 16;
 873
 874     prod ^= htable[a >> 56];
 875     a <<= 8;
 876     prod ^= ltable[a >> 56];
 877     prod ^= ((xor) ? *d64 : 0);
 878     *d64 = prod;
 879     s64++;
 880     d64++;
 881   }
 882 */
 883
 884   while (d64 != top64) {
 885     a = *s64;
 886
 887     prod = 0;
 888     for (j = 0; j < 4; j++) {
 889       prod <<= 16;
 890       prod ^= htable[a >> 56];
 891       a <<= 8;
 892       prod ^= ltable[a >> 56];
 893       a <<= 8;
 894     }
 895
 896     //JSP: We can move the conditional outside the while loop, but we need to fully test it to understand which is better.
 897
 898     prod ^= ((xor) ? *d64 : 0);
 899     *d64 = prod;
 900     s64++;
 901     d64++;
 902   }
 903   gf_do_final_region_alignment(&rd);
 904 }
 905
 906 static void
 907 gf_w16_table_lazy_multiply_region(gf_t *gf, void *src, void *dest, gf_val_32_t val, int bytes, int xor)
 908 {
 909   uint64_t c;
 910   gf_internal_t *h;
 911   struct gf_w16_lazytable_data *ltd;
 912   gf_region_data rd;
 913
 914   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
 915   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
 916
 917   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 8);
 918   gf_do_initial_region_alignment(&rd);
 919
 920   h = (gf_internal_t *) gf->scratch;
 921   ltd = (struct gf_w16_lazytable_data *) h->private;
 922
 923   ltd->lazytable[0] = 0;
 924
 925   /*
 926   a = val;
 927   c = 1;
 928   pp = h->prim_poly;
 929
 930   do {
 931     ltd->lazytable[c] = a;
 932     c <<= 1;
 933     if (c & (1 << GF_FIELD_WIDTH)) c ^= pp;
 934     a <<= 1;
 935     if (a & (1 << GF_FIELD_WIDTH)) a ^= pp;
 936   } while (c != 1);
 937   */
 938
 939   for (c = 1; c < GF_FIELD_SIZE; c++) {
 940     ltd->lazytable[c] = gf_w16_shift_multiply(gf, c, val);
 941   }
 942
 943   gf_two_byte_region_table_multiply(&rd, ltd->lazytable);
 944   gf_do_final_region_alignment(&rd);
 945 }
 946
 947 #ifdef INTEL_SSSE3
 948 static
 949 void
 950 gf_w16_split_4_16_lazy_sse_multiply_region(gf_t *gf, void *src, void *dest, gf_val_32_t val, int bytes, int xor)
 951 {
 952   uint64_t i, j, *s64, *d64, *top64;;
 953   uint64_t c, prod;
 954   uint8_t low[4][16];
 955   uint8_t high[4][16];
 956   gf_region_data rd;
 957
 958   __m128i  mask, ta, tb, ti, tpl, tph, tlow[4], thigh[4], tta, ttb, lmask;
 959
 960   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
 961   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
 962
 963   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 32);
 964   gf_do_initial_region_alignment(&rd);
 965
 966   for (j = 0; j < 16; j++) {
 967     for (i = 0; i < 4; i++) {
 968       c = (j << (i*4));
 969       prod = gf->multiply.w32(gf, c, val);
 970       low[i][j] = (prod & 0xff);
 971       high[i][j] = (prod >> 8);
 972     }
 973   }
 974
 975   for (i = 0; i < 4; i++) {
 976     tlow[i] = _mm_loadu_si128((__m128i *)low[i]);
 977     thigh[i] = _mm_loadu_si128((__m128i *)high[i]);
 978   }
 979
 980   s64 = (uint64_t *) rd.s_start;
 981   d64 = (uint64_t *) rd.d_start;
 982   top64 = (uint64_t *) rd.d_top;
 983
 984   mask = _mm_set1_epi8 (0x0f);
 985   lmask = _mm_set1_epi16 (0xff);
 986
 987   if (xor) {
 988     while (d64 != top64) {
 989
 990       ta = _mm_load_si128((__m128i *) s64);
 991       tb = _mm_load_si128((__m128i *) (s64+2));
 992
 993       tta = _mm_srli_epi16(ta, 8);
 994       ttb = _mm_srli_epi16(tb, 8);
 995       tpl = _mm_and_si128(tb, lmask);
 996       tph = _mm_and_si128(ta, lmask);
 997
 998       tb = _mm_packus_epi16(tpl, tph);
 999       ta = _mm_packus_epi16(ttb, tta);
1000
1001       ti = _mm_and_si128 (mask, tb);
1002       tph = _mm_shuffle_epi8 (thigh[0], ti);
1003       tpl = _mm_shuffle_epi8 (tlow[0], ti);
1004
1005       tb = _mm_srli_epi16(tb, 4);
1006       ti = _mm_and_si128 (mask, tb);
1007       tpl = _mm_xor_si128(_mm_shuffle_epi8 (tlow[1], ti), tpl);
1008       tph = _mm_xor_si128(_mm_shuffle_epi8 (thigh[1], ti), tph);
1009
1010       ti = _mm_and_si128 (mask, ta);
1011       tpl = _mm_xor_si128(_mm_shuffle_epi8 (tlow[2], ti), tpl);
1012       tph = _mm_xor_si128(_mm_shuffle_epi8 (thigh[2], ti), tph);
1013
1014       ta = _mm_srli_epi16(ta, 4);
1015       ti = _mm_and_si128 (mask, ta);
1016       tpl = _mm_xor_si128(_mm_shuffle_epi8 (tlow[3], ti), tpl);
1017       tph = _mm_xor_si128(_mm_shuffle_epi8 (thigh[3], ti), tph);
1018
1019       ta = _mm_unpackhi_epi8(tpl, tph);
1020       tb = _mm_unpacklo_epi8(tpl, tph);
1021
1022       tta = _mm_load_si128((__m128i *) d64);
1023       ta = _mm_xor_si128(ta, tta);
1024       ttb = _mm_load_si128((__m128i *) (d64+2));
1025       tb = _mm_xor_si128(tb, ttb);
1026       _mm_store_si128 ((__m128i *)d64, ta);
1027       _mm_store_si128 ((__m128i *)(d64+2), tb);
1028
1029       d64 += 4;
1030       s64 += 4;
1031
1032     }
1033   } else {
1034     while (d64 != top64) {
1035
1036       ta = _mm_load_si128((__m128i *) s64);
1037       tb = _mm_load_si128((__m128i *) (s64+2));
1038
1039       tta = _mm_srli_epi16(ta, 8);
1040       ttb = _mm_srli_epi16(tb, 8);
1041       tpl = _mm_and_si128(tb, lmask);
1042       tph = _mm_and_si128(ta, lmask);
1043
1044       tb = _mm_packus_epi16(tpl, tph);
1045       ta = _mm_packus_epi16(ttb, tta);
1046
1047       ti = _mm_and_si128 (mask, tb);
1048       tph = _mm_shuffle_epi8 (thigh[0], ti);
1049       tpl = _mm_shuffle_epi8 (tlow[0], ti);
1050
1051       tb = _mm_srli_epi16(tb, 4);
1052       ti = _mm_and_si128 (mask, tb);
1053       tpl = _mm_xor_si128(_mm_shuffle_epi8 (tlow[1], ti), tpl);
1054       tph = _mm_xor_si128(_mm_shuffle_epi8 (thigh[1], ti), tph);
1055
1056       ti = _mm_and_si128 (mask, ta);
1057       tpl = _mm_xor_si128(_mm_shuffle_epi8 (tlow[2], ti), tpl);
1058       tph = _mm_xor_si128(_mm_shuffle_epi8 (thigh[2], ti), tph);
1059
1060       ta = _mm_srli_epi16(ta, 4);
1061       ti = _mm_and_si128 (mask, ta);
1062       tpl = _mm_xor_si128(_mm_shuffle_epi8 (tlow[3], ti), tpl);
1063       tph = _mm_xor_si128(_mm_shuffle_epi8 (thigh[3], ti), tph);
1064
1065       ta = _mm_unpackhi_epi8(tpl, tph);
1066       tb = _mm_unpacklo_epi8(tpl, tph);
1067
1068       _mm_store_si128 ((__m128i *)d64, ta);
1069       _mm_store_si128 ((__m128i *)(d64+2), tb);
1070
1071       d64 += 4;
1072       s64 += 4;
1073     }
1074   }
1075
1076   gf_do_final_region_alignment(&rd);
1077 }
1078 #endif
1079
1080 #ifdef INTEL_SSSE3
1081 static
1082 void
1083 gf_w16_split_4_16_lazy_sse_altmap_multiply_region(gf_t *gf, void *src, void *dest, gf_val_32_t val, int bytes, int xor)
1084 {
1085   uint64_t i, j, *s64, *d64, *top64;;
1086   uint64_t c, prod;
1087   uint8_t low[4][16];
1088   uint8_t high[4][16];
1089   gf_region_data rd;
1090   __m128i  mask, ta, tb, ti, tpl, tph, tlow[4], thigh[4];
1091
1092   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
1093   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
1094
1095   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 32);
1096   gf_do_initial_region_alignment(&rd);
1097
1098   for (j = 0; j < 16; j++) {
1099     for (i = 0; i < 4; i++) {
1100       c = (j << (i*4));
1101       prod = gf->multiply.w32(gf, c, val);
1102       low[i][j] = (prod & 0xff);
1103       high[i][j] = (prod >> 8);
1104     }
1105   }
1106
1107   for (i = 0; i < 4; i++) {
1108     tlow[i] = _mm_loadu_si128((__m128i *)low[i]);
1109     thigh[i] = _mm_loadu_si128((__m128i *)high[i]);
1110   }
1111
1112   s64 = (uint64_t *) rd.s_start;
1113   d64 = (uint64_t *) rd.d_start;
1114   top64 = (uint64_t *) rd.d_top;
1115
1116   mask = _mm_set1_epi8 (0x0f);
1117
1118   if (xor) {
1119     while (d64 != top64) {
1120
1121       ta = _mm_load_si128((__m128i *) s64);
1122       tb = _mm_load_si128((__m128i *) (s64+2));
1123
1124       ti = _mm_and_si128 (mask, tb);
1125       tph = _mm_shuffle_epi8 (thigh[0], ti);
1126       tpl = _mm_shuffle_epi8 (tlow[0], ti);
1127
1128       tb = _mm_srli_epi16(tb, 4);
1129       ti = _mm_and_si128 (mask, tb);
1130       tpl = _mm_xor_si128(_mm_shuffle_epi8 (tlow[1], ti), tpl);
1131       tph = _mm_xor_si128(_mm_shuffle_epi8 (thigh[1], ti), tph);
1132
1133       ti = _mm_and_si128 (mask, ta);
1134       tpl = _mm_xor_si128(_mm_shuffle_epi8 (tlow[2], ti), tpl);
1135       tph = _mm_xor_si128(_mm_shuffle_epi8 (thigh[2], ti), tph);
1136
1137       ta = _mm_srli_epi16(ta, 4);
1138       ti = _mm_and_si128 (mask, ta);
1139       tpl = _mm_xor_si128(_mm_shuffle_epi8 (tlow[3], ti), tpl);
1140       tph = _mm_xor_si128(_mm_shuffle_epi8 (thigh[3], ti), tph);
1141
1142       ta = _mm_load_si128((__m128i *) d64);
1143       tph = _mm_xor_si128(tph, ta);
1144       _mm_store_si128 ((__m128i *)d64, tph);
1145       tb = _mm_load_si128((__m128i *) (d64+2));
1146       tpl = _mm_xor_si128(tpl, tb);
1147       _mm_store_si128 ((__m128i *)(d64+2), tpl);
1148
1149       d64 += 4;
1150       s64 += 4;
1151     }
1152   } else {
1153     while (d64 != top64) {
1154
1155       ta = _mm_load_si128((__m128i *) s64);
1156       tb = _mm_load_si128((__m128i *) (s64+2));
1157
1158       ti = _mm_and_si128 (mask, tb);
1159       tph = _mm_shuffle_epi8 (thigh[0], ti);
1160       tpl = _mm_shuffle_epi8 (tlow[0], ti);
1161
1162       tb = _mm_srli_epi16(tb, 4);
1163       ti = _mm_and_si128 (mask, tb);
1164       tpl = _mm_xor_si128(_mm_shuffle_epi8 (tlow[1], ti), tpl);
1165       tph = _mm_xor_si128(_mm_shuffle_epi8 (thigh[1], ti), tph);
1166
1167       ti = _mm_and_si128 (mask, ta);
1168       tpl = _mm_xor_si128(_mm_shuffle_epi8 (tlow[2], ti), tpl);
1169       tph = _mm_xor_si128(_mm_shuffle_epi8 (thigh[2], ti), tph);
1170
1171       ta = _mm_srli_epi16(ta, 4);
1172       ti = _mm_and_si128 (mask, ta);
1173       tpl = _mm_xor_si128(_mm_shuffle_epi8 (tlow[3], ti), tpl);
1174       tph = _mm_xor_si128(_mm_shuffle_epi8 (thigh[3], ti), tph);
1175
1176       _mm_store_si128 ((__m128i *)d64, tph);
1177       _mm_store_si128 ((__m128i *)(d64+2), tpl);
1178
1179       d64 += 4;
1180       s64 += 4;
1181
1182     }
1183   }
1184   gf_do_final_region_alignment(&rd);
1185
1186 }
1187 #endif
1188
1189 uint32_t
1190 gf_w16_split_8_8_multiply(gf_t *gf, gf_val_32_t a, gf_val_32_t b)
1191 {
1192   uint32_t alow, blow;
1193   struct gf_w16_split_8_8_data *d8;
1194   gf_internal_t *h;
1195
1196   h = (gf_internal_t *) gf->scratch;
1197   d8 = (struct gf_w16_split_8_8_data *) h->private;
1198
1199   alow = a & 0xff;
1200   blow = b & 0xff;
1201   a >>= 8;
1202   b >>= 8;
1203
1204   return d8->tables[0][alow][blow] ^
1205          d8->tables[1][alow][b] ^
1206          d8->tables[1][a][blow] ^
1207          d8->tables[2][a][b];
1208 }
1209
1210 static
1211 int gf_w16_split_init(gf_t *gf)
1212 {
1213   gf_internal_t *h;
1214   struct gf_w16_split_8_8_data *d8;
1215   int i, j, exp;
1216   uint32_t p, basep, tmp;
1217
1218   h = (gf_internal_t *) gf->scratch;
1219
1220   if (h->arg1 == 8 && h->arg2 == 8) {
1221     d8 = (struct gf_w16_split_8_8_data *) h->private;
1222     basep = 1;
1223     for (exp = 0; exp < 3; exp++) {
1224       for (j = 0; j < 256; j++) d8->tables[exp][0][j] = 0;
1225       for (i = 0; i < 256; i++) d8->tables[exp][i][0] = 0;
1226       d8->tables[exp][1][1] = basep;
1227       for (i = 2; i < 256; i++) {
1228         if (i&1) {
1229           p = d8->tables[exp][i^1][1];
1230           d8->tables[exp][i][1] = p ^ basep;
1231         } else {
1232           p = d8->tables[exp][i>>1][1];
1233           d8->tables[exp][i][1] = GF_MULTBY_TWO(p);
1234         }
1235       }
1236       for (i = 1; i < 256; i++) {
1237         p = d8->tables[exp][i][1];
1238         for (j = 1; j < 256; j++) {
1239           if (j&1) {
1240             d8->tables[exp][i][j] = d8->tables[exp][i][j^1] ^ p;
1241           } else {
1242             tmp = d8->tables[exp][i][j>>1];
1243             d8->tables[exp][i][j] = GF_MULTBY_TWO(tmp);
1244           }
1245         }
1246       }
1247       for (i = 0; i < 8; i++) basep = GF_MULTBY_TWO(basep);
1248     }
1249     SET_FUNCTION(gf,multiply,w32,gf_w16_split_8_8_multiply)
1250     SET_FUNCTION(gf,multiply_region,w32,gf_w16_split_8_16_lazy_multiply_region)
1251     return 1;
1252
1253   }
1254
1255   /* We'll be using LOG for multiplication, unless the pp isn't primitive.
1256      In that case, we'll be using SHIFT. */
1257
1258   gf_w16_log_init(gf);
1259
1260   /* Defaults */
1261
1262 #ifdef INTEL_SSSE3
1263   if (gf_cpu_supports_intel_ssse3) {
1264     SET_FUNCTION(gf,multiply_region,w32,gf_w16_split_4_16_lazy_sse_multiply_region)
1265   } else {
1266 #elif ARM_NEON
1267   if (gf_cpu_supports_arm_neon) {
1268     gf_w16_neon_split_init(gf);
1269   } else {
1270 #endif
1271     SET_FUNCTION(gf,multiply_region,w32,gf_w16_split_8_16_lazy_multiply_region)
1272 #if defined(INTEL_SSSE3) || defined(ARM_NEON)
1273   }
1274 #endif
1275
1276   if ((h->arg1 == 8 && h->arg2 == 16) || (h->arg2 == 8 && h->arg1 == 16)) {
1277     SET_FUNCTION(gf,multiply_region,w32,gf_w16_split_8_16_lazy_multiply_region)
1278
1279   } else if ((h->arg1 == 4 && h->arg2 == 16) || (h->arg2 == 4 && h->arg1 == 16)) {
1280 #if defined(INTEL_SSSE3) || defined(ARM_NEON)
1281     if (gf_cpu_supports_intel_ssse3 || gf_cpu_supports_arm_neon) {
1282       if(h->region_type & GF_REGION_ALTMAP && h->region_type & GF_REGION_NOSIMD)
1283         SET_FUNCTION(gf,multiply_region,w32,gf_w16_split_4_16_lazy_nosse_altmap_multiply_region)
1284       else if(h->region_type & GF_REGION_NOSIMD)
1285         SET_FUNCTION(gf,multiply_region,w32,gf_w16_split_4_16_lazy_multiply_region)
1286 #if defined(INTEL_SSSE3)
1287       else if(h->region_type & GF_REGION_ALTMAP && gf_cpu_supports_intel_ssse3)
1288         SET_FUNCTION(gf,multiply_region,w32,gf_w16_split_4_16_lazy_sse_altmap_multiply_region)
1289 #endif
1290     } else {
1291 #endif
1292       if(h->region_type & GF_REGION_SIMD)
1293         return 0;
1294       else if(h->region_type & GF_REGION_ALTMAP)
1295         SET_FUNCTION(gf,multiply_region,w32,gf_w16_split_4_16_lazy_nosse_altmap_multiply_region)
1296       else
1297         SET_FUNCTION(gf,multiply_region,w32,gf_w16_split_4_16_lazy_multiply_region)
1298 #if defined(INTEL_SSSE3) || defined(ARM_NEON)
1299     }
1300 #endif
1301   }
1302
1303   return 1;
1304 }
1305
1306 static
1307 int gf_w16_table_init(gf_t *gf)
1308 {
1309   gf_w16_log_init(gf);
1310
1311   SET_FUNCTION(gf,multiply_region,w32,gf_w16_table_lazy_multiply_region)
1312   return 1;
1313 }
1314
1315 static
1316 void
1317 gf_w16_log_zero_multiply_region(gf_t *gf, void *src, void *dest, gf_val_32_t val, int bytes, int xor)
1318 {
1319   uint16_t lv;
1320   int i;
1321   uint16_t *s16, *d16, *top16;
1322   struct gf_w16_zero_logtable_data *ltd;
1323   gf_region_data rd;
1324
1325   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
1326   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
1327
1328   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 2);
1329   gf_do_initial_region_alignment(&rd);
1330
1331   ltd = (struct gf_w16_zero_logtable_data*) ((gf_internal_t *) gf->scratch)->private;
1332   s16 = (uint16_t *) rd.s_start;
1333   d16 = (uint16_t *) rd.d_start;
1334   top16 = (uint16_t *) rd.d_top;
1335   bytes = top16 - d16;
1336
1337   lv = ltd->log_tbl[val];
1338
1339   if (xor) {
1340     for (i = 0; i < bytes; i++) {
1341       d16[i] ^= (ltd->antilog_tbl[lv + ltd->log_tbl[s16[i]]]);
1342     }
1343   } else {
1344     for (i = 0; i < bytes; i++) {
1345       d16[i] = (ltd->antilog_tbl[lv + ltd->log_tbl[s16[i]]]);
1346     }
1347   }
1348
1349   /* This isn't necessary. */
1350
1351   gf_do_final_region_alignment(&rd);
1352 }
1353
1354 /* Here -- double-check Kevin */
1355
1356 static
1357 inline
1358 gf_val_32_t
1359 gf_w16_log_zero_multiply (gf_t *gf, gf_val_32_t a, gf_val_32_t b)
1360 {
1361   struct gf_w16_zero_logtable_data *ltd;
1362
1363   ltd = (struct gf_w16_zero_logtable_data *) ((gf_internal_t *) gf->scratch)->private;
1364   return ltd->antilog_tbl[ltd->log_tbl[a] + ltd->log_tbl[b]];
1365 }
1366
1367 static
1368 inline
1369 gf_val_32_t
1370 gf_w16_log_zero_divide (gf_t *gf, gf_val_32_t a, gf_val_32_t b)
1371 {
1372   int log_sum = 0;
1373   struct gf_w16_zero_logtable_data *ltd;
1374
1375   if (a == 0 || b == 0) return 0;
1376   ltd = (struct gf_w16_zero_logtable_data *) ((gf_internal_t *) gf->scratch)->private;
1377
1378   log_sum = ltd->log_tbl[a] - ltd->log_tbl[b] + (GF_MULT_GROUP_SIZE);
1379   return (ltd->antilog_tbl[log_sum]);
1380 }
1381
1382 static
1383 gf_val_32_t
1384 gf_w16_log_zero_inverse (gf_t *gf, gf_val_32_t a)
1385 {
1386   struct gf_w16_zero_logtable_data *ltd;
1387
1388   ltd = (struct gf_w16_zero_logtable_data *) ((gf_internal_t *) gf->scratch)->private;
1389   return (ltd->inv_tbl[a]);
1390 }
1391
1392 static
1393 inline
1394 gf_val_32_t
1395 gf_w16_bytwo_p_multiply (gf_t *gf, gf_val_32_t a, gf_val_32_t b)
1396 {
1397   uint32_t prod, pp, pmask, amask;
1398   gf_internal_t *h;
1399
1400   h = (gf_internal_t *) gf->scratch;
1401   pp = h->prim_poly;
1402
1403
1404   prod = 0;
1405   pmask = 0x8000;
1406   amask = 0x8000;
1407
1408   while (amask != 0) {
1409     if (prod & pmask) {
1410       prod = ((prod << 1) ^ pp);
1411     } else {
1412       prod <<= 1;
1413     }
1414     if (a & amask) prod ^= b;
1415     amask >>= 1;
1416   }
1417   return prod;
1418 }
1419
1420 static
1421 inline
1422 gf_val_32_t
1423 gf_w16_bytwo_b_multiply (gf_t *gf, gf_val_32_t a, gf_val_32_t b)
1424 {
1425   uint32_t prod, pp, bmask;
1426   gf_internal_t *h;
1427
1428   h = (gf_internal_t *) gf->scratch;
1429   pp = h->prim_poly;
1430
1431   prod = 0;
1432   bmask = 0x8000;
1433
1434   while (1) {
1435     if (a & 1) prod ^= b;
1436     a >>= 1;
1437     if (a == 0) return prod;
1438     if (b & bmask) {
1439       b = ((b << 1) ^ pp);
1440     } else {
1441       b <<= 1;
1442     }
1443   }
1444 }
1445
1446 static
1447 void
1448 gf_w16_bytwo_p_nosse_multiply_region(gf_t *gf, void *src, void *dest, gf_val_32_t val, int bytes, int xor)
1449 {
1450   uint64_t *s64, *d64, t1, t2, ta, prod, amask;
1451   gf_region_data rd;
1452   struct gf_w16_bytwo_data *btd;
1453
1454   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
1455   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
1456
1457   btd = (struct gf_w16_bytwo_data *) ((gf_internal_t *) (gf->scratch))->private;
1458
1459   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 8);
1460   gf_do_initial_region_alignment(&rd);
1461
1462   s64 = (uint64_t *) rd.s_start;
1463   d64 = (uint64_t *) rd.d_start;
1464
1465   if (xor) {
1466     while (s64 < (uint64_t *) rd.s_top) {
1467       prod = 0;
1468       amask = 0x8000;
1469       ta = *s64;
1470       while (amask != 0) {
1471         AB2(btd->prim_poly, btd->mask1, btd->mask2, prod, t1, t2);
1472         if (val & amask) prod ^= ta;
1473         amask >>= 1;
1474       }
1475       *d64 ^= prod;
1476       d64++;
1477       s64++;
1478     }
1479   } else {
1480     while (s64 < (uint64_t *) rd.s_top) {
1481       prod = 0;
1482       amask = 0x8000;
1483       ta = *s64;
1484       while (amask != 0) {
1485         AB2(btd->prim_poly, btd->mask1, btd->mask2, prod, t1, t2);
1486         if (val & amask) prod ^= ta;
1487         amask >>= 1;
1488       }
1489       *d64 = prod;
1490       d64++;
1491       s64++;
1492     }
1493   }
1494   gf_do_final_region_alignment(&rd);
1495 }
1496
1497 #define BYTWO_P_ONESTEP {\
1498       SSE_AB2(pp, m1 ,m2, prod, t1, t2); \
1499       t1 = _mm_and_si128(v, one); \
1500       t1 = _mm_sub_epi16(t1, one); \
1501       t1 = _mm_and_si128(t1, ta); \
1502       prod = _mm_xor_si128(prod, t1); \
1503       v = _mm_srli_epi64(v, 1); }
1504
1505 #ifdef INTEL_SSE2
1506 static
1507 void
1508 gf_w16_bytwo_p_sse_multiply_region(gf_t *gf, void *src, void *dest, gf_val_32_t val, int bytes, int xor)
1509 {
1510   int i;
1511   uint8_t *s8, *d8;
1512   uint32_t vrev;
1513   __m128i pp, m1, m2, ta, prod, t1, t2, tp, one, v;
1514   struct gf_w16_bytwo_data *btd;
1515   gf_region_data rd;
1516
1517   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
1518   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
1519
1520   btd = (struct gf_w16_bytwo_data *) ((gf_internal_t *) (gf->scratch))->private;
1521
1522   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 16);
1523   gf_do_initial_region_alignment(&rd);
1524
1525   vrev = 0;
1526   for (i = 0; i < 16; i++) {
1527     vrev <<= 1;
1528     if (!(val & (1 << i))) vrev |= 1;
1529   }
1530
1531   s8 = (uint8_t *) rd.s_start;
1532   d8 = (uint8_t *) rd.d_start;
1533
1534   pp = _mm_set1_epi16(btd->prim_poly&0xffff);
1535   m1 = _mm_set1_epi16((btd->mask1)&0xffff);
1536   m2 = _mm_set1_epi16((btd->mask2)&0xffff);
1537   one = _mm_set1_epi16(1);
1538
1539   while (d8 < (uint8_t *) rd.d_top) {
1540     prod = _mm_setzero_si128();
1541     v = _mm_set1_epi16(vrev);
1542     ta = _mm_load_si128((__m128i *) s8);
1543     tp = (!xor) ? _mm_setzero_si128() : _mm_load_si128((__m128i *) d8);
1544     BYTWO_P_ONESTEP;
1545     BYTWO_P_ONESTEP;
1546     BYTWO_P_ONESTEP;
1547     BYTWO_P_ONESTEP;
1548     BYTWO_P_ONESTEP;
1549     BYTWO_P_ONESTEP;
1550     BYTWO_P_ONESTEP;
1551     BYTWO_P_ONESTEP;
1552     BYTWO_P_ONESTEP;
1553     BYTWO_P_ONESTEP;
1554     BYTWO_P_ONESTEP;
1555     BYTWO_P_ONESTEP;
1556     BYTWO_P_ONESTEP;
1557     BYTWO_P_ONESTEP;
1558     BYTWO_P_ONESTEP;
1559     BYTWO_P_ONESTEP;
1560     _mm_store_si128((__m128i *) d8, _mm_xor_si128(prod, tp));
1561     d8 += 16;
1562     s8 += 16;
1563   }
1564   gf_do_final_region_alignment(&rd);
1565 }
1566 #endif
1567
1568 #ifdef INTEL_SSE2
1569 static
1570 void
1571 gf_w16_bytwo_b_sse_region_2_noxor(gf_region_data *rd, struct gf_w16_bytwo_data *btd)
1572 {
1573   uint8_t *d8, *s8;
1574   __m128i pp, m1, m2, t1, t2, va;
1575
1576   s8 = (uint8_t *) rd->s_start;
1577   d8 = (uint8_t *) rd->d_start;
1578
1579   pp = _mm_set1_epi16(btd->prim_poly&0xffff);
1580   m1 = _mm_set1_epi16((btd->mask1)&0xffff);
1581   m2 = _mm_set1_epi16((btd->mask2)&0xffff);
1582
1583   while (d8 < (uint8_t *) rd->d_top) {
1584     va = _mm_load_si128 ((__m128i *)(s8));
1585     SSE_AB2(pp, m1, m2, va, t1, t2);
1586     _mm_store_si128((__m128i *)d8, va);
1587     d8 += 16;
1588     s8 += 16;
1589   }
1590 }
1591 #endif
1592
1593 #ifdef INTEL_SSE2
1594 static
1595 void
1596 gf_w16_bytwo_b_sse_region_2_xor(gf_region_data *rd, struct gf_w16_bytwo_data *btd)
1597 {
1598   uint8_t *d8, *s8;
1599   __m128i pp, m1, m2, t1, t2, va, vb;
1600
1601   s8 = (uint8_t *) rd->s_start;
1602   d8 = (uint8_t *) rd->d_start;
1603
1604   pp = _mm_set1_epi16(btd->prim_poly&0xffff);
1605   m1 = _mm_set1_epi16((btd->mask1)&0xffff);
1606   m2 = _mm_set1_epi16((btd->mask2)&0xffff);
1607
1608   while (d8 < (uint8_t *) rd->d_top) {
1609     va = _mm_load_si128 ((__m128i *)(s8));
1610     SSE_AB2(pp, m1, m2, va, t1, t2);
1611     vb = _mm_load_si128 ((__m128i *)(d8));
1612     vb = _mm_xor_si128(vb, va);
1613     _mm_store_si128((__m128i *)d8, vb);
1614     d8 += 16;
1615     s8 += 16;
1616   }
1617 }
1618 #endif
1619
1620
1621 #ifdef INTEL_SSE2
1622 static
1623 void
1624 gf_w16_bytwo_b_sse_multiply_region(gf_t *gf, void *src, void *dest, gf_val_32_t val, int bytes, int xor)
1625 {
1626   int itb;
1627   uint8_t *d8, *s8;
1628   __m128i pp, m1, m2, t1, t2, va, vb;
1629   struct gf_w16_bytwo_data *btd;
1630   gf_region_data rd;
1631
1632   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
1633   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
1634
1635   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 16);
1636   gf_do_initial_region_alignment(&rd);
1637
1638   btd = (struct gf_w16_bytwo_data *) ((gf_internal_t *) (gf->scratch))->private;
1639
1640   if (val == 2) {
1641     if (xor) {
1642       gf_w16_bytwo_b_sse_region_2_xor(&rd, btd);
1643     } else {
1644       gf_w16_bytwo_b_sse_region_2_noxor(&rd, btd);
1645     }
1646     gf_do_final_region_alignment(&rd);
1647     return;
1648   }
1649
1650   s8 = (uint8_t *) rd.s_start;
1651   d8 = (uint8_t *) rd.d_start;
1652
1653   pp = _mm_set1_epi16(btd->prim_poly&0xffff);
1654   m1 = _mm_set1_epi16((btd->mask1)&0xffff);
1655   m2 = _mm_set1_epi16((btd->mask2)&0xffff);
1656
1657   while (d8 < (uint8_t *) rd.d_top) {
1658     va = _mm_load_si128 ((__m128i *)(s8));
1659     vb = (!xor) ? _mm_setzero_si128() : _mm_load_si128 ((__m128i *)(d8));
1660     itb = val;
1661     while (1) {
1662       if (itb & 1) vb = _mm_xor_si128(vb, va);
1663       itb >>= 1;
1664       if (itb == 0) break;
1665       SSE_AB2(pp, m1, m2, va, t1, t2);
1666     }
1667     _mm_store_si128((__m128i *)d8, vb);
1668     d8 += 16;
1669     s8 += 16;
1670   }
1671
1672   gf_do_final_region_alignment(&rd);
1673 }
1674 #endif
1675
1676 static
1677 void
1678 gf_w16_bytwo_b_nosse_multiply_region(gf_t *gf, void *src, void *dest, gf_val_32_t val, int bytes, int xor)
1679 {
1680   uint64_t *s64, *d64, t1, t2, ta, tb, prod;
1681   struct gf_w16_bytwo_data *btd;
1682   gf_region_data rd;
1683
1684   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
1685   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
1686
1687   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 16);
1688   gf_do_initial_region_alignment(&rd);
1689
1690   btd = (struct gf_w16_bytwo_data *) ((gf_internal_t *) (gf->scratch))->private;
1691   s64 = (uint64_t *) rd.s_start;
1692   d64 = (uint64_t *) rd.d_start;
1693
1694   switch (val) {
1695   case 2:
1696     if (xor) {
1697       while (d64 < (uint64_t *) rd.d_top) {
1698         ta = *s64;
1699         AB2(btd->prim_poly, btd->mask1, btd->mask2, ta, t1, t2);
1700         *d64 ^= ta;
1701         d64++;
1702         s64++;
1703       }
1704     } else {
1705       while (d64 < (uint64_t *) rd.d_top) {
1706         ta = *s64;
1707         AB2(btd->prim_poly, btd->mask1, btd->mask2, ta, t1, t2);
1708         *d64 = ta;
1709         d64++;
1710         s64++;
1711       }
1712     }
1713     break;
1714   case 3:
1715     if (xor) {
1716       while (d64 < (uint64_t *) rd.d_top) {
1717         ta = *s64;
1718         prod = ta;
1719         AB2(btd->prim_poly, btd->mask1, btd->mask2, ta, t1, t2);
1720         *d64 ^= (ta ^ prod);
1721         d64++;
1722         s64++;
1723       }
1724     } else {
1725       while (d64 < (uint64_t *) rd.d_top) {
1726         ta = *s64;
1727         prod = ta;
1728         AB2(btd->prim_poly, btd->mask1, btd->mask2, ta, t1, t2);
1729         *d64 = (ta ^ prod);
1730         d64++;
1731         s64++;
1732       }
1733     }
1734     break;
1735   case 4:
1736     if (xor) {
1737       while (d64 < (uint64_t *) rd.d_top) {
1738         ta = *s64;
1739         AB2(btd->prim_poly, btd->mask1, btd->mask2, ta, t1, t2);
1740         AB2(btd->prim_poly, btd->mask1, btd->mask2, ta, t1, t2);
1741         *d64 ^= ta;
1742         d64++;
1743         s64++;
1744       }
1745     } else {
1746       while (d64 < (uint64_t *) rd.d_top) {
1747         ta = *s64;
1748         AB2(btd->prim_poly, btd->mask1, btd->mask2, ta, t1, t2);
1749         AB2(btd->prim_poly, btd->mask1, btd->mask2, ta, t1, t2);
1750         *d64 = ta;
1751         d64++;
1752         s64++;
1753       }
1754     }
1755     break;
1756   case 5:
1757     if (xor) {
1758       while (d64 < (uint64_t *) rd.d_top) {
1759         ta = *s64;
1760         prod = ta;
1761         AB2(btd->prim_poly, btd->mask1, btd->mask2, ta, t1, t2);
1762         AB2(btd->prim_poly, btd->mask1, btd->mask2, ta, t1, t2);
1763         *d64 ^= (ta ^ prod);
1764         d64++;
1765         s64++;
1766       }
1767     } else {
1768       while (d64 < (uint64_t *) rd.d_top) {
1769         ta = *s64;
1770         prod = ta;
1771         AB2(btd->prim_poly, btd->mask1, btd->mask2, ta, t1, t2);
1772         AB2(btd->prim_poly, btd->mask1, btd->mask2, ta, t1, t2);
1773         *d64 = ta ^ prod;
1774         d64++;
1775         s64++;
1776       }
1777     }
1778     break;
1779   default:
1780     if (xor) {
1781       while (d64 < (uint64_t *) rd.d_top) {
1782         prod = *d64 ;
1783         ta = *s64;
1784         tb = val;
1785         while (1) {
1786           if (tb & 1) prod ^= ta;
1787           tb >>= 1;
1788           if (tb == 0) break;
1789           AB2(btd->prim_poly, btd->mask1, btd->mask2, ta, t1, t2);
1790         }
1791         *d64 = prod;
1792         d64++;
1793         s64++;
1794       }
1795     } else {
1796       while (d64 < (uint64_t *) rd.d_top) {
1797         prod = 0 ;
1798         ta = *s64;
1799         tb = val;
1800         while (1) {
1801           if (tb & 1) prod ^= ta;
1802           tb >>= 1;
1803           if (tb == 0) break;
1804           AB2(btd->prim_poly, btd->mask1, btd->mask2, ta, t1, t2);
1805         }
1806         *d64 = prod;
1807         d64++;
1808         s64++;
1809       }
1810     }
1811     break;
1812   }
1813   gf_do_final_region_alignment(&rd);
1814 }
1815
1816 static
1817 int gf_w16_bytwo_init(gf_t *gf)
1818 {
1819   gf_internal_t *h;
1820   uint64_t ip, m1, m2;
1821   struct gf_w16_bytwo_data *btd;
1822
1823   h = (gf_internal_t *) gf->scratch;
1824   btd = (struct gf_w16_bytwo_data *) (h->private);
1825   ip = h->prim_poly & 0xffff;
1826   m1 = 0xfffe;
1827   m2 = 0x8000;
1828   btd->prim_poly = 0;
1829   btd->mask1 = 0;
1830   btd->mask2 = 0;
1831
1832   while (ip != 0) {
1833     btd->prim_poly |= ip;
1834     btd->mask1 |= m1;
1835     btd->mask2 |= m2;
1836     ip <<= GF_FIELD_WIDTH;
1837     m1 <<= GF_FIELD_WIDTH;
1838     m2 <<= GF_FIELD_WIDTH;
1839   }
1840
1841   if (h->mult_type == GF_MULT_BYTWO_p) {
1842     SET_FUNCTION(gf,multiply,w32,gf_w16_bytwo_p_multiply)
1843     #ifdef INTEL_SSE2
1844     if (gf_cpu_supports_intel_sse2 && !(h->region_type & GF_REGION_NOSIMD)) {
1845       SET_FUNCTION(gf,multiply_region,w32,gf_w16_bytwo_p_sse_multiply_region)
1846     } else {
1847     #endif
1848       SET_FUNCTION(gf,multiply_region,w32,gf_w16_bytwo_p_nosse_multiply_region)
1849       if(h->region_type & GF_REGION_SIMD)
1850         return 0;
1851     #ifdef INTEL_SSE2
1852     }
1853     #endif
1854   } else {
1855     SET_FUNCTION(gf,multiply,w32,gf_w16_bytwo_b_multiply)
1856     #ifdef INTEL_SSE2
1857     if (gf_cpu_supports_intel_sse2 && !(h->region_type & GF_REGION_NOSIMD)) {
1858         SET_FUNCTION(gf,multiply_region,w32,gf_w16_bytwo_b_sse_multiply_region)
1859     } else {
1860     #endif
1861       SET_FUNCTION(gf,multiply_region,w32,gf_w16_bytwo_b_nosse_multiply_region)
1862       if(h->region_type & GF_REGION_SIMD)
1863         return 0;
1864     #ifdef INTEL_SSE2
1865     }
1866     #endif
1867   }
1868
1869   return 1;
1870 }
1871
1872 static
1873 int gf_w16_log_zero_init(gf_t *gf)
1874 {
1875   gf_internal_t *h;
1876   struct gf_w16_zero_logtable_data *ltd;
1877   int i, b;
1878
1879   h = (gf_internal_t *) gf->scratch;
1880   ltd = h->private;
1881
1882   ltd->log_tbl[0] = (-GF_MULT_GROUP_SIZE) + 1;
1883
1884   bzero(&(ltd->_antilog_tbl[0]), sizeof(ltd->_antilog_tbl));
1885
1886   ltd->antilog_tbl = &(ltd->_antilog_tbl[GF_FIELD_SIZE * 2]);
1887
1888   b = 1;
1889   for (i = 0; i < GF_MULT_GROUP_SIZE; i++) {
1890       ltd->log_tbl[b] = (uint16_t)i;
1891       ltd->antilog_tbl[i] = (uint16_t)b;
1892       ltd->antilog_tbl[i+GF_MULT_GROUP_SIZE] = (uint16_t)b;
1893       b <<= 1;
1894       if (b & GF_FIELD_SIZE) {
1895           b = b ^ h->prim_poly;
1896       }
1897   }
1898   ltd->inv_tbl[0] = 0;  /* Not really, but we need to fill it with something  */
1899   ltd->inv_tbl[1] = 1;
1900   for (i = 2; i < GF_FIELD_SIZE; i++) {
1901     ltd->inv_tbl[i] = ltd->antilog_tbl[GF_MULT_GROUP_SIZE-ltd->log_tbl[i]];
1902   }
1903
1904   SET_FUNCTION(gf,inverse,w32,gf_w16_log_zero_inverse)
1905   SET_FUNCTION(gf,divide,w32,gf_w16_log_zero_divide)
1906   SET_FUNCTION(gf,multiply,w32,gf_w16_log_zero_multiply)
1907   SET_FUNCTION(gf,multiply_region,w32,gf_w16_log_zero_multiply_region)
1908   return 1;
1909 }
1910
1911 static
1912 gf_val_32_t
1913 gf_w16_composite_multiply_recursive(gf_t *gf, gf_val_32_t a, gf_val_32_t b)
1914 {
1915   gf_internal_t *h = (gf_internal_t *) gf->scratch;
1916   gf_t *base_gf = h->base_gf;
1917   uint8_t b0 = b & 0x00ff;
1918   uint8_t b1 = (b & 0xff00) >> 8;
1919   uint8_t a0 = a & 0x00ff;
1920   uint8_t a1 = (a & 0xff00) >> 8;
1921   uint8_t a1b1;
1922   uint16_t rv;
1923
1924   a1b1 = base_gf->multiply.w32(base_gf, a1, b1);
1925
1926   rv = ((base_gf->multiply.w32(base_gf, a0, b0) ^ a1b1) | ((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)) << 8));
1927   return rv;
1928 }
1929
1930 static
1931 gf_val_32_t
1932 gf_w16_composite_multiply_inline(gf_t *gf, gf_val_32_t a, gf_val_32_t b)
1933 {
1934   gf_internal_t *h = (gf_internal_t *) gf->scratch;
1935   uint8_t b0 = b & 0x00ff;
1936   uint8_t b1 = (b & 0xff00) >> 8;
1937   uint8_t a0 = a & 0x00ff;
1938   uint8_t a1 = (a & 0xff00) >> 8;
1939   uint8_t a1b1, *mt;
1940   uint16_t rv;
1941   struct gf_w16_composite_data *cd;
1942
1943   cd = (struct gf_w16_composite_data *) h->private;
1944   mt = cd->mult_table;
1945
1946   a1b1 = GF_W8_INLINE_MULTDIV(mt, a1, b1);
1947
1948   rv = ((GF_W8_INLINE_MULTDIV(mt, a0, b0) ^ a1b1) | ((GF_W8_INLINE_MULTDIV(mt, a1, b0) ^ GF_W8_INLINE_MULTDIV(mt, a0, b1) ^ GF_W8_INLINE_MULTDIV(mt, a1b1, h->prim_poly)) << 8));
1949   return rv;
1950 }
1951
1952 /*
1953  * Composite field division trick (explained in 2007 tech report)
1954  *
1955  * Compute a / b = a*b^-1, where p(x) = x^2 + sx + 1
1956  *
1957  * let c = b^-1
1958  *
1959  * c*b = (s*b1c1+b1c0+b0c1)x+(b1c1+b0c0)
1960  *
1961  * want (s*b1c1+b1c0+b0c1) = 0 and (b1c1+b0c0) = 1
1962  *
1963  * let d = b1c1 and d+1 = b0c0
1964  *
1965  * solve s*b1c1+b1c0+b0c1 = 0
1966  *
1967  * solution: d = (b1b0^-1)(b1b0^-1+b0b1^-1+s)^-1
1968  *
1969  * c0 = (d+1)b0^-1
1970  * c1 = d*b1^-1
1971  *
1972  * a / b = a * c
1973  */
1974
1975 static
1976 gf_val_32_t
1977 gf_w16_composite_inverse(gf_t *gf, gf_val_32_t a)
1978 {
1979   gf_internal_t *h = (gf_internal_t *) gf->scratch;
1980   gf_t *base_gf = h->base_gf;
1981   uint8_t a0 = a & 0x00ff;
1982   uint8_t a1 = (a & 0xff00) >> 8;
1983   uint8_t c0, c1, d, tmp;
1984   uint16_t c;
1985   uint8_t a0inv, a1inv;
1986
1987   if (a0 == 0) {
1988     a1inv = base_gf->inverse.w32(base_gf, a1);
1989     c0 = base_gf->multiply.w32(base_gf, a1inv, h->prim_poly);
1990     c1 = a1inv;
1991   } else if (a1 == 0) {
1992     c0 = base_gf->inverse.w32(base_gf, a0);
1993     c1 = 0;
1994   } else {
1995     a1inv = base_gf->inverse.w32(base_gf, a1);
1996     a0inv = base_gf->inverse.w32(base_gf, a0);
1997
1998     d = base_gf->multiply.w32(base_gf, a1, a0inv);
1999
2000     tmp = (base_gf->multiply.w32(base_gf, a1, a0inv) ^ base_gf->multiply.w32(base_gf, a0, a1inv) ^ h->prim_poly);
2001     tmp = base_gf->inverse.w32(base_gf, tmp);
2002
2003     d = base_gf->multiply.w32(base_gf, d, tmp);
2004
2005     c0 = base_gf->multiply.w32(base_gf, (d^1), a0inv);
2006     c1 = base_gf->multiply.w32(base_gf, d, a1inv);
2007   }
2008
2009   c = c0 | (c1 << 8);
2010
2011   return c;
2012 }
2013
2014 static
2015 void
2016 gf_w16_composite_multiply_region(gf_t *gf, void *src, void *dest, gf_val_32_t val, int bytes, int xor)
2017 {
2018   gf_internal_t *h = (gf_internal_t *) gf->scratch;
2019   gf_t *base_gf = h->base_gf;
2020   uint8_t b0 = val & 0x00ff;
2021   uint8_t b1 = (val & 0xff00) >> 8;
2022   uint16_t *s16, *d16, *top;
2023   uint8_t a0, a1, a1b1, *mt;
2024   gf_region_data rd;
2025   struct gf_w16_composite_data *cd;
2026
2027   cd = (struct gf_w16_composite_data *) h->private;
2028   mt = cd->mult_table;
2029
2030   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
2031   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 2);
2032
2033   s16 = rd.s_start;
2034   d16 = rd.d_start;
2035   top = rd.d_top;
2036
2037   if (mt == NULL) {
2038     if (xor) {
2039       while (d16 < top) {
2040         a0 = (*s16) & 0x00ff;
2041         a1 = ((*s16) & 0xff00) >> 8;
2042         a1b1 = base_gf->multiply.w32(base_gf, a1, b1);
2043
2044         (*d16) ^= ((base_gf->multiply.w32(base_gf, a0, b0) ^ a1b1) |
2045                   ((base_gf->multiply.w32(base_gf, a1, b0) ^
2046                     base_gf->multiply.w32(base_gf, a0, b1) ^
2047                     base_gf->multiply.w32(base_gf, a1b1, h->prim_poly)) << 8));
2048         s16++;
2049         d16++;
2050       }
2051     } else {
2052       while (d16 < top) {
2053         a0 = (*s16) & 0x00ff;
2054         a1 = ((*s16) & 0xff00) >> 8;
2055         a1b1 = base_gf->multiply.w32(base_gf, a1, b1);
2056
2057         (*d16) = ((base_gf->multiply.w32(base_gf, a0, b0) ^ a1b1) |
2058                   ((base_gf->multiply.w32(base_gf, a1, b0) ^
2059                     base_gf->multiply.w32(base_gf, a0, b1) ^
2060                     base_gf->multiply.w32(base_gf, a1b1, h->prim_poly)) << 8));
2061         s16++;
2062         d16++;
2063       }
2064     }
2065   } else {
2066     if (xor) {
2067       while (d16 < top) {
2068         a0 = (*s16) & 0x00ff;
2069         a1 = ((*s16) & 0xff00) >> 8;
2070         a1b1 = GF_W8_INLINE_MULTDIV(mt, a1, b1);
2071
2072         (*d16) ^= ((GF_W8_INLINE_MULTDIV(mt, a0, b0) ^ a1b1) |
2073                   ((GF_W8_INLINE_MULTDIV(mt, a1, b0) ^
2074                     GF_W8_INLINE_MULTDIV(mt, a0, b1) ^
2075                     GF_W8_INLINE_MULTDIV(mt, a1b1, h->prim_poly)) << 8));
2076         s16++;
2077         d16++;
2078       }
2079     } else {
2080       while (d16 < top) {
2081         a0 = (*s16) & 0x00ff;
2082         a1 = ((*s16) & 0xff00) >> 8;
2083         a1b1 = GF_W8_INLINE_MULTDIV(mt, a1, b1);
2084
2085         (*d16) = ((GF_W8_INLINE_MULTDIV(mt, a0, b0) ^ a1b1) |
2086                   ((GF_W8_INLINE_MULTDIV(mt, a1, b0) ^
2087                     GF_W8_INLINE_MULTDIV(mt, a0, b1) ^
2088                     GF_W8_INLINE_MULTDIV(mt, a1b1, h->prim_poly)) << 8));
2089         s16++;
2090         d16++;
2091       }
2092     }
2093   }
2094 }
2095
2096 static
2097 void
2098 gf_w16_composite_multiply_region_alt(gf_t *gf, void *src, void *dest, gf_val_32_t val, int bytes, int xor)
2099 {
2100   gf_internal_t *h = (gf_internal_t *) gf->scratch;
2101   gf_t *base_gf = h->base_gf;
2102   uint8_t val0 = val & 0x00ff;
2103   uint8_t val1 = (val & 0xff00) >> 8;
2104   gf_region_data rd;
2105   int sub_reg_size;
2106   uint8_t *slow, *shigh;
2107   uint8_t *dlow, *dhigh, *top;;
2108
2109   /* JSP: I want the two pointers aligned wrt each other on 16 byte
2110      boundaries.  So I'm going to make sure that the area on
2111      which the two operate is a multiple of 32. Of course, that
2112      junks up the mapping, but so be it -- that's why we have extract_word.... */
2113
2114   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 32);
2115   gf_do_initial_region_alignment(&rd);
2116
2117   slow = (uint8_t *) rd.s_start;
2118   dlow = (uint8_t *) rd.d_start;
2119   top = (uint8_t *)  rd.d_top;
2120   sub_reg_size = (top - dlow)/2;
2121   shigh = slow + sub_reg_size;
2122   dhigh = dlow + sub_reg_size;
2123
2124   base_gf->multiply_region.w32(base_gf, slow, dlow, val0, sub_reg_size, xor);
2125   base_gf->multiply_region.w32(base_gf, shigh, dlow, val1, sub_reg_size, 1);
2126   base_gf->multiply_region.w32(base_gf, slow, dhigh, val1, sub_reg_size, xor);
2127   base_gf->multiply_region.w32(base_gf, shigh, dhigh, val0, sub_reg_size, 1);
2128   base_gf->multiply_region.w32(base_gf, shigh, dhigh, base_gf->multiply.w32(base_gf, h->prim_poly, val1), sub_reg_size, 1);
2129
2130   gf_do_final_region_alignment(&rd);
2131 }
2132
2133 static
2134 int gf_w16_composite_init(gf_t *gf)
2135 {
2136   gf_internal_t *h = (gf_internal_t *) gf->scratch;
2137   struct gf_w16_composite_data *cd;
2138
2139   if (h->base_gf == NULL) return 0;
2140
2141   cd = (struct gf_w16_composite_data *) h->private;
2142   cd->mult_table = gf_w8_get_mult_table(h->base_gf);
2143
2144   if (h->region_type & GF_REGION_ALTMAP) {
2145     SET_FUNCTION(gf,multiply_region,w32,gf_w16_composite_multiply_region_alt)
2146   } else {
2147     SET_FUNCTION(gf,multiply_region,w32,gf_w16_composite_multiply_region)
2148   }
2149
2150   if (cd->mult_table == NULL) {
2151     SET_FUNCTION(gf,multiply,w32,gf_w16_composite_multiply_recursive)
2152   } else {
2153     SET_FUNCTION(gf,multiply,w32,gf_w16_composite_multiply_inline)
2154   }
2155   SET_FUNCTION(gf,divide,w32,NULL)
2156   SET_FUNCTION(gf,inverse,w32,gf_w16_composite_inverse)
2157
2158   return 1;
2159 }
2160
2161 static
2162 void
2163 gf_w16_group_4_set_shift_tables(uint16_t *shift, uint16_t val, gf_internal_t *h)
2164 {
2165   int i, j;
2166
2167   shift[0] = 0;
2168   for (i = 0; i < 16; i += 2) {
2169     j = (shift[i>>1] << 1);
2170     if (j & (1 << 16)) j ^= h->prim_poly;
2171     shift[i] = j;
2172     shift[i^1] = j^val;
2173   }
2174 }
2175
2176 static
2177 inline
2178 gf_val_32_t
2179 gf_w16_group_4_4_multiply(gf_t *gf, gf_val_32_t a, gf_val_32_t b)
2180 {
2181   uint16_t p, l, ind, r, a16;
2182
2183   struct gf_w16_group_4_4_data *d44;
2184   gf_internal_t *h = (gf_internal_t *) gf->scratch;
2185
2186   d44 = (struct gf_w16_group_4_4_data *) h->private;
2187   gf_w16_group_4_set_shift_tables(d44->shift, b, h);
2188
2189   a16 = a;
2190   ind = a16 >> 12;
2191   a16 <<= 4;
2192   p = d44->shift[ind];
2193   r = p & 0xfff;
2194   l = p >> 12;
2195   ind = a16 >> 12;
2196   a16 <<= 4;
2197   p = (d44->shift[ind] ^ d44->reduce[l] ^ (r << 4));
2198   r = p & 0xfff;
2199   l = p >> 12;
2200   ind = a16 >> 12;
2201   a16 <<= 4;
2202   p = (d44->shift[ind] ^ d44->reduce[l] ^ (r << 4));
2203   r = p & 0xfff;
2204   l = p >> 12;
2205   ind = a16 >> 12;
2206   p = (d44->shift[ind] ^ d44->reduce[l] ^ (r << 4));
2207   return p;
2208 }
2209
2210 static
2211 void gf_w16_group_4_4_region_multiply(gf_t *gf, void *src, void *dest, gf_val_32_t val, int bytes, int xor)
2212 {
2213   uint16_t p, l, ind, r, a16, p16;
2214   struct gf_w16_group_4_4_data *d44;
2215   gf_region_data rd;
2216   uint16_t *s16, *d16, *top;
2217
2218   if (val == 0) { gf_multby_zero(dest, bytes, xor); return; }
2219   if (val == 1) { gf_multby_one(src, dest, bytes, xor); return; }
2220
2221   gf_internal_t *h = (gf_internal_t *) gf->scratch;
2222   d44 = (struct gf_w16_group_4_4_data *) h->private;
2223   gf_w16_group_4_set_shift_tables(d44->shift, val, h);
2224
2225   gf_set_region_data(&rd, gf, src, dest, bytes, val, xor, 2);
2226   gf_do_initial_region_alignment(&rd);
2227
2228   s16 = (uint16_t *) rd.s_start;
2229   d16 = (uint16_t *) rd.d_start;
2230   top = (uint16_t *) rd.d_top;
2231
2232   while (d16 < top) {
2233     a16 = *s16;
2234     p16 = (xor) ? *d16 : 0;
2235     ind = a16 >> 12;
2236     a16 <<= 4;
2237     p = d44->shift[ind];
2238     r = p & 0xfff;
2239     l = p >> 12;
2240     ind = a16 >> 12;
2241     a16 <<= 4;
2242     p = (d44->shift[ind] ^ d44->reduce[l] ^ (r << 4));
2243     r = p & 0xfff;
2244     l = p >> 12;
2245     ind = a16 >> 12;
2246     a16 <<= 4;
2247     p = (d44->shift[ind] ^ d44->reduce[l] ^ (r << 4));
2248     r = p & 0xfff;
2249     l = p >> 12;
2250     ind = a16 >> 12;
2251     p = (d44->shift[ind] ^ d44->reduce[l] ^ (r << 4));
2252     p ^= p16;
2253     *d16 = p;
2254     d16++;
2255     s16++;
2256   }
2257   gf_do_final_region_alignment(&rd);
2258 }
2259
2260 static
2261 int gf_w16_group_init(gf_t *gf)
2262 {
2263   int i, j, p;
2264   struct gf_w16_group_4_4_data *d44;
2265   gf_internal_t *h = (gf_internal_t *) gf->scratch;
2266
2267   d44 = (struct gf_w16_group_4_4_data *) h->private;
2268   d44->reduce[0] = 0;
2269   for (i = 0; i < 16; i++) {
2270     p = 0;
2271     for (j = 0; j < 4; j++) {
2272       if (i & (1 << j)) p ^= (h->prim_poly << j);
2273     }
2274     d44->reduce[p>>16] = (p&0xffff);
2275   }
2276
2277   SET_FUNCTION(gf,multiply,w32,gf_w16_group_4_4_multiply)
2278   SET_FUNCTION(gf,divide,w32,NULL)
2279   SET_FUNCTION(gf,inverse,w32,NULL)
2280   SET_FUNCTION(gf,multiply_region,w32,gf_w16_group_4_4_region_multiply)
2281
2282   return 1;
2283 }
2284
2285 int gf_w16_scratch_size(int mult_type, int region_type, int divide_type, int arg1, int arg2)
2286 {
2287   switch(mult_type)
2288   {
2289     case GF_MULT_TABLE:
2290       return sizeof(gf_internal_t) + sizeof(struct gf_w16_lazytable_data) + 64;
2291       break;
2292     case GF_MULT_BYTWO_p:
2293     case GF_MULT_BYTWO_b:
2294       return sizeof(gf_internal_t) + sizeof(struct gf_w16_bytwo_data);
2295       break;
2296     case GF_MULT_LOG_ZERO:
2297       return sizeof(gf_internal_t) + sizeof(struct gf_w16_zero_logtable_data) + 64;
2298       break;
2299     case GF_MULT_LOG_TABLE:
2300       return sizeof(gf_internal_t) + sizeof(struct gf_w16_logtable_data) + 64;
2301       break;
2302     case GF_MULT_DEFAULT:
2303     case GF_MULT_SPLIT_TABLE:
2304       if (arg1 == 8 && arg2 == 8) {
2305         return sizeof(gf_internal_t) + sizeof(struct gf_w16_split_8_8_data) + 64;
2306       } else if ((arg1 == 8 && arg2 == 16) || (arg2 == 8 && arg1 == 16)) {
2307         return sizeof(gf_internal_t) + sizeof(struct gf_w16_logtable_data) + 64;
2308       } else if (mult_type == GF_MULT_DEFAULT ||
2309                  (arg1 == 4 && arg2 == 16) || (arg2 == 4 && arg1 == 16)) {
2310         return sizeof(gf_internal_t) + sizeof(struct gf_w16_logtable_data) + 64;
2311       }
2312       return 0;
2313       break;
2314     case GF_MULT_GROUP:
2315       return sizeof(gf_internal_t) + sizeof(struct gf_w16_group_4_4_data) + 64;
2316       break;
2317     case GF_MULT_CARRY_FREE:
2318       return sizeof(gf_internal_t);
2319       break;
2320     case GF_MULT_SHIFT:
2321       return sizeof(gf_internal_t);
2322       break;
2323     case GF_MULT_COMPOSITE:
2324       return sizeof(gf_internal_t) + sizeof(struct gf_w16_composite_data) + 64;
2325       break;
2326
2327     default:
2328       return 0;
2329    }
2330    return 0;
2331 }
2332
2333 int gf_w16_init(gf_t *gf)
2334 {
2335   gf_internal_t *h;
2336
2337   h = (gf_internal_t *) gf->scratch;
2338
2339   /* Allen: set default primitive polynomial / irreducible polynomial if needed */
2340
2341   if (h->prim_poly == 0) {
2342     if (h->mult_type == GF_MULT_COMPOSITE) {
2343       h->prim_poly = gf_composite_get_default_poly(h->base_gf);
2344       if (h->prim_poly == 0) return 0;
2345     } else {
2346
2347      /* Allen: use the following primitive polynomial to make
2348                carryless multiply work more efficiently for GF(2^16).
2349
2350         h->prim_poly = 0x1002d;
2351
2352         The following is the traditional primitive polynomial for GF(2^16) */
2353
2354       h->prim_poly = 0x1100b;
2355     }
2356   }
2357
2358   if (h->mult_type != GF_MULT_COMPOSITE) h->prim_poly |= (1 << 16);
2359
2360   SET_FUNCTION(gf,multiply,w32,NULL)
2361   SET_FUNCTION(gf,divide,w32,NULL)
2362   SET_FUNCTION(gf,inverse,w32,NULL)
2363   SET_FUNCTION(gf,multiply_region,w32,NULL)
2364
2365   switch(h->mult_type) {
2366     case GF_MULT_LOG_ZERO:    if (gf_w16_log_zero_init(gf) == 0) return 0; break;
2367     case GF_MULT_LOG_TABLE:   if (gf_w16_log_init(gf) == 0) return 0; break;
2368     case GF_MULT_DEFAULT:
2369     case GF_MULT_SPLIT_TABLE: if (gf_w16_split_init(gf) == 0) return 0; break;
2370     case GF_MULT_TABLE:       if (gf_w16_table_init(gf) == 0) return 0; break;
2371     case GF_MULT_CARRY_FREE:  if (gf_w16_cfm_init(gf) == 0) return 0; break;
2372     case GF_MULT_SHIFT:       if (gf_w16_shift_init(gf) == 0) return 0; break;
2373     case GF_MULT_COMPOSITE:   if (gf_w16_composite_init(gf) == 0) return 0; break;
2374     case GF_MULT_BYTWO_p:
2375     case GF_MULT_BYTWO_b:     if (gf_w16_bytwo_init(gf) == 0) return 0; break;
2376     case GF_MULT_GROUP:       if (gf_w16_group_init(gf) == 0) return 0; break;
2377     default: return 0;
2378   }
2379   if (h->divide_type == GF_DIVIDE_EUCLID) {
2380     SET_FUNCTION(gf,divide,w32,gf_w16_divide_from_inverse)
2381     SET_FUNCTION(gf,inverse,w32,gf_w16_euclid)
2382   } else if (h->divide_type == GF_DIVIDE_MATRIX) {
2383     SET_FUNCTION(gf,divide,w32,gf_w16_divide_from_inverse)
2384     SET_FUNCTION(gf,inverse,w32,gf_w16_matrix)
2385   }
2386
2387   if (gf->divide.w32 == NULL) {
2388     SET_FUNCTION(gf,divide,w32,gf_w16_divide_from_inverse)
2389     if (gf->inverse.w32 == NULL) SET_FUNCTION(gf,inverse,w32,gf_w16_euclid)
2390   }
2391
2392   if (gf->inverse.w32 == NULL)  SET_FUNCTION(gf,inverse,w32,gf_w16_inverse_from_divide)
2393
2394   if (h->region_type & GF_REGION_ALTMAP) {
2395     if (h->mult_type == GF_MULT_COMPOSITE) {
2396       SET_FUNCTION(gf,extract_word,w32,gf_w16_composite_extract_word)
2397     } else {
2398       SET_FUNCTION(gf,extract_word,w32,gf_w16_split_extract_word)
2399     }
2400   } else if (h->region_type == GF_REGION_CAUCHY) {
2401     SET_FUNCTION(gf,multiply_region,w32,gf_wgen_cauchy_region)
2402     SET_FUNCTION(gf,extract_word,w32,gf_wgen_extract_word)
2403   } else {
2404     SET_FUNCTION(gf,extract_word,w32,gf_w16_extract_word)
2405   }
2406   if (gf->multiply_region.w32 == NULL) {
2407     SET_FUNCTION(gf,multiply_region,w32,gf_w16_multiply_region_from_single)
2408   }
2409   return 1;
2410 }
2411
2412 /* Inline setup functions */
2413
2414 uint16_t *gf_w16_get_log_table(gf_t *gf)
2415 {
2416   struct gf_w16_logtable_data *ltd;
2417
2418   if (gf->multiply.w32 == gf_w16_log_multiply) {
2419     ltd = (struct gf_w16_logtable_data *) ((gf_internal_t *) gf->scratch)->private;
2420     return (uint16_t *) ltd->log_tbl;
2421   }
2422   return NULL;
2423 }
2424
2425 uint16_t *gf_w16_get_mult_alog_table(gf_t *gf)
2426 {
2427   gf_internal_t *h;
2428   struct gf_w16_logtable_data *ltd;
2429
2430   h = (gf_internal_t *) gf->scratch;
2431   if (gf->multiply.w32 == gf_w16_log_multiply) {
2432     ltd = (struct gf_w16_logtable_data *) h->private;
2433     return (uint16_t *) ltd->antilog_tbl;
2434   }
2435   return NULL;
2436 }
2437
2438 uint16_t *gf_w16_get_div_alog_table(gf_t *gf)
2439 {
2440   gf_internal_t *h;
2441   struct gf_w16_logtable_data *ltd;
2442
2443   h = (gf_internal_t *) gf->scratch;
2444   if (gf->multiply.w32 == gf_w16_log_multiply) {
2445     ltd = (struct gf_w16_logtable_data *) h->private;
2446     return (uint16_t *) ltd->d_antilog;
2447   }
2448   return NULL;
2449 }