]> AND Private Git Repository - Cipher_code.git/blob - IDA_new/gf-complete/src/gf_w16.c
Logo AND Algorithmique Numérique Distribuée

Private GIT Repository
update
[Cipher_code.git] / 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 }