source: git/ntl/src/GF2X.c @ 26e030

spielwiese
Last change on this file since 26e030 was 26e030, checked in by Hans Schönemann <hannes@…>, 15 years ago
*hannes: update to 5.5.1 git-svn-id: file:///usr/local/Singular/svn/trunk@11949 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 36.8 KB
Line 
1
2#include <NTL/GF2X.h>
3#include <NTL/vec_long.h>
4
5#include <NTL/new.h>
6#include <stdio.h>
7
8#if (defined(NTL_WIZARD_HACK) && defined(NTL_GF2X_LIB))
9#undef NTL_GF2X_LIB
10#endif
11
12#ifdef NTL_GF2X_LIB
13#include <gf2x.h>
14#endif
15
16NTL_START_IMPL
17
18long GF2X::HexOutput = 0;
19
20
21void GF2X::SetMaxLength(long n)
22{
23   if (n < 0) Error("GF2X::SetMaxLength: negative length");
24   if (NTL_OVERFLOW(n, 1, 0))
25      Error("GF2X::SetMaxLength: excessive length");
26   long w = (n + NTL_BITS_PER_LONG - 1)/NTL_BITS_PER_LONG;
27   xrep.SetMaxLength(w);
28}
29
30GF2X::GF2X(INIT_SIZE_TYPE, long n)
31{
32   SetMaxLength(n);
33}
34
35
36
37const GF2X& GF2X::zero()
38{
39   static GF2X z;
40   return z;
41}
42
43void GF2X::normalize()
44{
45   long n;
46   const _ntl_ulong *p;
47
48   n = xrep.length();
49   if (n == 0) return;
50   p = xrep.elts() + n;
51   while (n > 0 && (*--p) == 0) {
52      n--;
53   }
54   xrep.QuickSetLength(n);
55}
56
57long IsZero(const GF2X& a) 
58   { return a.xrep.length() == 0; }
59
60long IsOne(const GF2X& a)
61   { return a.xrep.length() == 1 && a.xrep[0] == 1; }
62
63
64
65
66
67
68
69long IsX(const GF2X& a)
70{
71   return a.xrep.length() == 1 && a.xrep[0] == 2;
72}
73
74GF2 coeff(const GF2X& a, long i)
75{
76   if (i < 0) return to_GF2(0);
77   long wi = i/NTL_BITS_PER_LONG;
78   if (wi >= a.xrep.length()) return to_GF2(0);
79   long bi = i - wi*NTL_BITS_PER_LONG;
80
81   return to_GF2((a.xrep[wi] & (1UL << bi)) != 0);
82}
83
84GF2 LeadCoeff(const GF2X& a)
85{
86   if (IsZero(a))
87      return to_GF2(0);
88   else
89      return to_GF2(1);
90}
91
92GF2 ConstTerm(const GF2X& a)
93{
94   if (IsZero(a))
95      return to_GF2(0);
96   else
97      return to_GF2((a.xrep[0] & 1) != 0);
98}
99
100
101void set(GF2X& x)
102{
103   x.xrep.SetLength(1);
104   x.xrep[0] = 1;
105}
106
107void SetX(GF2X& x)
108{
109   x.xrep.SetLength(1);
110   x.xrep[0] = 2;
111}
112
113void SetCoeff(GF2X& x, long i)
114{
115   if (i < 0) {
116      Error("SetCoeff: negative index");
117      return;
118   }
119
120   long n, j;
121
122   n = x.xrep.length();
123   long wi = i/NTL_BITS_PER_LONG;
124
125   if (wi >= n) {
126      x.xrep.SetLength(wi+1);
127      for (j = n; j <= wi; j++)
128         x.xrep[j] = 0;
129   }
130
131   long bi = i - wi*NTL_BITS_PER_LONG;
132
133   x.xrep[wi] |= (1UL << bi);
134}
135   
136
137
138void SetCoeff(GF2X& x, long i, long val)
139{
140   if (i < 0) {
141      Error("SetCoeff: negative index");
142      return;
143   }
144
145   val = val & 1;
146
147   if (val) {
148      SetCoeff(x, i);
149      return;
150   }
151
152   // we want to clear position i
153
154   long n;
155
156   n = x.xrep.length();
157   long wi = i/NTL_BITS_PER_LONG;
158
159   if (wi >= n) 
160      return;
161
162   long bi = i - wi*NTL_BITS_PER_LONG;
163
164   x.xrep[wi] &= ~(1UL << bi);
165   if (wi == n-1) x.normalize();
166}
167
168void SetCoeff(GF2X& x, long i, GF2 a)
169{
170   SetCoeff(x, i, rep(a));
171}
172
173
174void swap(GF2X& a, GF2X& b)
175{
176   swap(a.xrep, b.xrep);
177}
178
179
180
181long deg(const GF2X& aa)
182{
183   long n = aa.xrep.length();
184
185   if (n == 0)
186      return -1;
187
188   _ntl_ulong a = aa.xrep[n-1];
189   long i = 0;
190
191   if (a == 0) Error("GF2X: unnormalized polynomial detected in deg");
192
193   while (a>=256)
194      i += 8, a >>= 8;
195   if (a >=16)
196      i += 4, a >>= 4;
197   if (a >= 4)
198      i += 2, a >>= 2;
199   if (a >= 2)
200      i += 2;
201   else if (a >= 1)
202      i++;
203
204   return NTL_BITS_PER_LONG*(n-1) + i - 1;
205}
206   
207
208long operator==(const GF2X& a, const GF2X& b)
209{
210   return a.xrep == b.xrep;
211}
212
213long operator==(const GF2X& a, long b)
214{
215   if (b & 1) 
216      return IsOne(a);
217   else
218      return IsZero(a);
219}
220
221long operator==(const GF2X& a, GF2 b)
222{
223   if (b == 1) 
224      return IsOne(a);
225   else
226      return IsZero(a);
227}
228
229void random(GF2X& x, long n)
230{
231   if (n < 0) Error("GF2X random: negative length");
232
233   if (NTL_OVERFLOW(n, 1, 0))
234      Error("GF2X random: excessive length");
235
236   long wl = (n+NTL_BITS_PER_LONG-1)/NTL_BITS_PER_LONG;
237
238   x.xrep.SetLength(wl);
239
240   long i;
241   for (i = 0; i < wl-1; i++) {
242      x.xrep[i] = RandomWord();
243   }
244
245   if (n > 0) {
246      long pos = n % NTL_BITS_PER_LONG;
247      if (pos == 0) pos = NTL_BITS_PER_LONG;
248      x.xrep[wl-1] = RandomBits_ulong(pos);
249   }
250
251   x.normalize();
252}
253
254void add(GF2X& x, const GF2X& a, const GF2X& b)
255{
256   long sa = a.xrep.length();
257   long sb = b.xrep.length();
258
259   long i;
260
261   if (sa == sb) {
262      x.xrep.SetLength(sa);
263      if (sa == 0) return;
264
265      _ntl_ulong *xp = x.xrep.elts();
266      const _ntl_ulong *ap = a.xrep.elts();
267      const _ntl_ulong *bp = b.xrep.elts();
268
269      for (i = 0; i < sa; i++)
270         xp[i] = ap[i] ^ bp[i];
271
272      i = sa-1;
273      while (i >= 0 && !xp[i]) i--;
274      x.xrep.QuickSetLength(i+1);
275   }
276   
277   else if (sa < sb) {
278      x.xrep.SetLength(sb);
279      _ntl_ulong *xp = x.xrep.elts();
280      const _ntl_ulong *ap = a.xrep.elts();
281      const _ntl_ulong *bp = b.xrep.elts();
282
283      for (i = 0; i < sa; i++)
284         xp[i] = ap[i] ^ bp[i];
285
286      for (; i < sb; i++)
287         xp[i] = bp[i];
288   }
289   else { // sa > sb
290      x.xrep.SetLength(sa);
291      _ntl_ulong *xp = x.xrep.elts();
292      const _ntl_ulong *ap = a.xrep.elts();
293      const _ntl_ulong *bp = b.xrep.elts();
294
295      for (i = 0; i < sb; i++)
296         xp[i] = ap[i] ^ bp[i];
297
298      for (; i < sa; i++)
299         xp[i] = ap[i];
300   }
301}
302
303
304
305
306
307/*
308 * The bodies of mul1, Mul1, AddMul1, and mul_half
309 * are generated by the MakeDesc program, and the
310 * macros NTL_BB_MUL_CODE... are defined in mach_desc.h.
311 * Thanks to Paul Zimmermann for providing improvements
312 * to this approach.
313 */
314
315
316#if (defined(NTL_GF2X_ALTCODE1))
317
318#define NTL_EFF_BB_MUL_CODE0 NTL_ALT1_BB_MUL_CODE0
319#define NTL_EFF_BB_MUL_CODE1 NTL_ALT1_BB_MUL_CODE1
320#define NTL_EFF_BB_MUL_CODE2 NTL_ALT1_BB_MUL_CODE2
321#define NTL_EFF_SHORT_BB_MUL_CODE1 NTL_ALT1_SHORT_BB_MUL_CODE1
322#define NTL_EFF_HALF_BB_MUL_CODE0 NTL_ALT1_HALF_BB_MUL_CODE0
323
324#elif (defined(NTL_GF2X_ALTCODE))
325
326#define NTL_EFF_BB_MUL_CODE0 NTL_ALT_BB_MUL_CODE0
327#define NTL_EFF_BB_MUL_CODE1 NTL_ALT_BB_MUL_CODE1
328#define NTL_EFF_BB_MUL_CODE2 NTL_ALT_BB_MUL_CODE2
329#define NTL_EFF_SHORT_BB_MUL_CODE1 NTL_ALT_SHORT_BB_MUL_CODE1
330#define NTL_EFF_HALF_BB_MUL_CODE0 NTL_ALT_HALF_BB_MUL_CODE0
331
332#else
333
334#define NTL_EFF_BB_MUL_CODE0 NTL_BB_MUL_CODE0
335#define NTL_EFF_BB_MUL_CODE1 NTL_BB_MUL_CODE1
336#define NTL_EFF_BB_MUL_CODE2 NTL_BB_MUL_CODE2
337#define NTL_EFF_SHORT_BB_MUL_CODE1 NTL_SHORT_BB_MUL_CODE1
338#define NTL_EFF_HALF_BB_MUL_CODE0 NTL_HALF_BB_MUL_CODE0
339
340#endif
341
342
343
344static 
345void mul1(_ntl_ulong *c, _ntl_ulong a, _ntl_ulong b)
346{
347
348NTL_EFF_BB_MUL_CODE0
349
350
351}
352
353
354#ifdef NTL_GF2X_NOINLINE
355
356#define mul1_IL mul1
357
358#else
359
360static inline
361void mul1_inline(_ntl_ulong *c, _ntl_ulong a, _ntl_ulong b)
362{
363
364
365NTL_EFF_BB_MUL_CODE0
366
367
368}
369
370#define mul1_IL mul1_inline
371#endif
372
373
374static 
375void Mul1(_ntl_ulong *cp, const _ntl_ulong *bp, long sb, _ntl_ulong a)
376{
377 
378
379NTL_EFF_BB_MUL_CODE1
380
381
382}
383
384static 
385void AddMul1(_ntl_ulong *cp, const _ntl_ulong* bp, long sb, _ntl_ulong a)
386{
387
388
389NTL_EFF_BB_MUL_CODE2
390
391
392}
393
394
395static 
396void Mul1_short(_ntl_ulong *cp, const _ntl_ulong *bp, long sb, _ntl_ulong a)
397{
398 
399
400NTL_EFF_SHORT_BB_MUL_CODE1
401
402
403}
404
405
406
407
408
409static 
410void mul_half(_ntl_ulong *c, _ntl_ulong a, _ntl_ulong b)
411{
412
413
414NTL_EFF_HALF_BB_MUL_CODE0
415
416
417}
418
419
420// mul2...mul8 hard-code 2x2...8x8 word multiplies.
421// I adapted these routines from LiDIA (except mul3, see below).
422// Inlining these seems to hurt, not help.
423
424static
425void mul2(_ntl_ulong *c, const _ntl_ulong *a, const _ntl_ulong *b)
426{
427   _ntl_ulong hs0, hs1;
428   _ntl_ulong hl2[2];
429
430   hs0 = a[0] ^ a[1];
431   hs1 = b[0] ^ b[1];
432
433   mul1_IL(c, a[0], b[0]);
434   mul1_IL(c+2, a[1], b[1]);
435   mul1_IL(hl2, hs0, hs1);
436
437
438   hl2[0] = hl2[0] ^ c[0] ^ c[2];
439   hl2[1] = hl2[1] ^ c[1] ^ c[3];
440
441   c[1] ^= hl2[0];
442   c[2] ^= hl2[1];
443}
444
445
446/*
447 * This version of mul3 I got from Weimerskirch, Stebila,
448 * and Shantz, "Generic GF(2^m) arithmetic in software
449 * an its application to ECC" (ACISP 2003).
450 */
451
452static
453void mul3 (_ntl_ulong *c, const _ntl_ulong *a, const _ntl_ulong *b)
454{
455   _ntl_ulong d0[2], d1[2], d2[2], d01[2], d02[2], d12[2];
456
457   mul1_IL(d0, a[0], b[0]);
458   mul1_IL(d1, a[1], b[1]);
459   mul1_IL(d2, a[2], b[2]);
460   mul1_IL(d01, a[0]^a[1], b[0]^b[1]);
461   mul1_IL(d02, a[0]^a[2], b[0]^b[2]);
462   mul1_IL(d12, a[1]^a[2], b[1]^b[2]);
463
464   
465   c[0] = d0[0];
466   c[1] = d0[1] ^ d01[0] ^ d1[0] ^ d0[0];
467   c[2] = d01[1] ^ d1[1] ^ d0[1] ^ d02[0] ^ d2[0] ^ d0[0] ^ d1[0];
468   c[3] = d02[1] ^ d2[1] ^ d0[1] ^ d1[1] ^ d12[0] ^ d1[0] ^ d2[0];
469   c[4] = d12[1] ^ d1[1] ^ d2[1] ^ d2[0];
470   c[5] = d2[1];
471
472}
473
474static
475void mul4(_ntl_ulong *c, const _ntl_ulong *a, const _ntl_ulong *b)
476{
477   _ntl_ulong hs0[2], hs1[2];
478   _ntl_ulong hl2[4];
479
480   hs0[0] = a[0] ^ a[2];
481   hs0[1] = a[1] ^ a[3];
482   hs1[0] = b[0] ^ b[2];
483   hs1[1] = b[1] ^ b[3];
484
485   mul2(c, a, b);
486   mul2(c+4, a+2, b+2); 
487   mul2(hl2, hs0, hs1);
488
489   hl2[0] = hl2[0] ^ c[0] ^ c[4];
490   hl2[1] = hl2[1] ^ c[1] ^ c[5];
491   hl2[2] = hl2[2] ^ c[2] ^ c[6];
492   hl2[3] = hl2[3] ^ c[3] ^ c[7];
493   
494   c[2] ^= hl2[0];
495   c[3] ^= hl2[1];
496   c[4] ^= hl2[2];
497   c[5] ^= hl2[3];
498}
499
500static
501void mul5 (_ntl_ulong *c, const _ntl_ulong *a, const _ntl_ulong *b)
502{
503   _ntl_ulong hs0[3], hs1[3];
504   _ntl_ulong hl2[6];
505
506   hs0[0] = a[0] ^ a[3]; 
507   hs0[1] = a[1] ^ a[4];
508   hs0[2] = a[2];
509   hs1[0] = b[0] ^ b[3]; 
510   hs1[1] = b[1] ^ b[4];
511   hs1[2] = b[2];
512
513   mul3(c, a, b); 
514   mul3(hl2, hs0, hs1);
515   mul2(c+6, a+3, b+3);
516
517   hl2[0] = hl2[0] ^ c[0] ^ c[6]; 
518   hl2[1] = hl2[1] ^ c[1] ^ c[7];
519   hl2[2] = hl2[2] ^ c[2] ^ c[8];
520   hl2[3] = hl2[3] ^ c[3] ^ c[9];
521   hl2[4] = hl2[4] ^ c[4];
522   hl2[5] = hl2[5] ^ c[5];
523
524 
525   c[3] ^= hl2[0];
526   c[4] ^= hl2[1];
527   c[5] ^= hl2[2];
528   c[6] ^= hl2[3];
529   c[7] ^= hl2[4];
530   c[8] ^= hl2[5];
531}
532
533static
534void mul6(_ntl_ulong *c, const _ntl_ulong *a, const _ntl_ulong *b)
535{
536   _ntl_ulong hs0[3], hs1[3];
537   _ntl_ulong hl2[6];
538
539   hs0[0] = a[0] ^ a[3];   
540   hs0[1] = a[1] ^ a[4];
541   hs0[2] = a[2] ^ a[5];
542   hs1[0] = b[0] ^ b[3];   
543   hs1[1] = b[1] ^ b[4];
544   hs1[2] = b[2] ^ b[5];
545
546   mul3(c, a, b);   
547   mul3(c+6, a+3, b+3); 
548   mul3(hl2, hs0, hs1); 
549
550   hl2[0] = hl2[0] ^ c[0] ^ c[6];
551   hl2[1] = hl2[1] ^ c[1] ^ c[7];
552   hl2[2] = hl2[2] ^ c[2] ^ c[8];
553   hl2[3] = hl2[3] ^ c[3] ^ c[9];
554   hl2[4] = hl2[4] ^ c[4] ^ c[10];
555   hl2[5] = hl2[5] ^ c[5] ^ c[11];
556   
557   c[3] ^= hl2[0];
558   c[4] ^= hl2[1];
559   c[5] ^= hl2[2];
560   c[6] ^= hl2[3];
561   c[7] ^= hl2[4];
562   c[8] ^= hl2[5];
563}
564
565static
566void mul7(_ntl_ulong *c, const _ntl_ulong *a, const _ntl_ulong *b)
567{
568   _ntl_ulong hs0[4], hs1[4];
569   _ntl_ulong hl2[8];
570
571   hs0[0] = a[0] ^ a[4]; 
572   hs0[1] = a[1] ^ a[5];
573   hs0[2] = a[2] ^ a[6];
574   hs0[3] = a[3];
575   hs1[0] = b[0] ^ b[4]; 
576   hs1[1] = b[1] ^ b[5];
577   hs1[2] = b[2] ^ b[6];
578   hs1[3] = b[3];
579
580   mul4(c, a, b); 
581   mul4(hl2, hs0, hs1);
582   mul3(c+8, a+4, b+4); 
583
584   hl2[0] = hl2[0] ^ c[0] ^ c[8];
585   hl2[1] = hl2[1] ^ c[1] ^ c[9];
586   hl2[2] = hl2[2] ^ c[2] ^ c[10];
587   hl2[3] = hl2[3] ^ c[3] ^ c[11];
588   hl2[4] = hl2[4] ^ c[4] ^ c[12];
589   hl2[5] = hl2[5] ^ c[5] ^ c[13];
590   hl2[6] = hl2[6] ^ c[6];
591   hl2[7] = hl2[7] ^ c[7];
592   
593   c[4]  ^= hl2[0];
594   c[5]  ^= hl2[1];
595   c[6]  ^= hl2[2];
596   c[7]  ^= hl2[3];
597   c[8]  ^= hl2[4];
598   c[9]  ^= hl2[5];
599   c[10] ^= hl2[6];
600   c[11] ^= hl2[7];
601}
602
603static
604void mul8(_ntl_ulong *c, const _ntl_ulong *a, const _ntl_ulong *b)
605{
606   _ntl_ulong hs0[4], hs1[4];
607   _ntl_ulong hl2[8];
608
609   hs0[0] = a[0] ^ a[4]; 
610   hs0[1] = a[1] ^ a[5];
611   hs0[2] = a[2] ^ a[6];
612   hs0[3] = a[3] ^ a[7];
613   hs1[0] = b[0] ^ b[4]; 
614   hs1[1] = b[1] ^ b[5];
615   hs1[2] = b[2] ^ b[6];
616   hs1[3] = b[3] ^ b[7];
617
618   mul4(c, a, b); 
619   mul4(c+8, a+4, b+4);
620   mul4(hl2, hs0, hs1); 
621
622   hl2[0] = hl2[0] ^ c[0] ^ c[8];
623   hl2[1] = hl2[1] ^ c[1] ^ c[9];
624   hl2[2] = hl2[2] ^ c[2] ^ c[10];
625   hl2[3] = hl2[3] ^ c[3] ^ c[11];
626   hl2[4] = hl2[4] ^ c[4] ^ c[12];
627   hl2[5] = hl2[5] ^ c[5] ^ c[13];
628   hl2[6] = hl2[6] ^ c[6] ^ c[14];
629   hl2[7] = hl2[7] ^ c[7] ^ c[15];
630   
631   c[4]  ^= hl2[0];
632   c[5]  ^= hl2[1];
633   c[6]  ^= hl2[2];
634   c[7]  ^= hl2[3];
635   c[8]  ^= hl2[4];
636   c[9]  ^= hl2[5];
637   c[10] ^= hl2[6];
638   c[11] ^= hl2[7];
639}
640
641static
642void KarMul(_ntl_ulong *c, const _ntl_ulong *a, const _ntl_ulong *b, 
643            long len, _ntl_ulong *stk)
644{
645   if (len <= 8) {
646      switch (len) {
647         case 1: mul1(c, a[0], b[0]); break;
648         case 2: mul2(c, a, b); break;
649         case 3: mul3(c, a, b); break;
650         case 4: mul4(c, a, b); break;
651         case 5: mul5(c, a, b); break;
652         case 6: mul6(c, a, b); break;
653         case 7: mul7(c, a, b); break;
654         case 8: mul8(c, a, b); break;
655      }
656
657      return;
658   }
659
660   long ll, lh, i, ll2, lh2;
661   const _ntl_ulong *a0, *a1, *b0, *b1;
662   _ntl_ulong *a01, *b01, *h;
663
664
665   lh = len >> 1;
666   ll = (len+1) >> 1;
667
668   ll2 = ll << 1;
669   lh2 = lh << 1;
670
671   a01 = stk;  stk += ll+1;
672   b01 = stk;  stk += ll+1;
673   h = stk; stk += ll2+1;
674
675   a0 = a;
676   a1 = a+ll;
677   b0 = b;
678   b1 = b+ll;
679
680   KarMul(c, a0, b0, ll, stk);
681   KarMul(c+ll2, a1, b1, lh, stk);
682
683   for (i = 0; i < lh; i++) {
684      a01[i] = a[i] ^ a[i+ll];
685      b01[i] = b[i] ^ b[i+ll];
686   }
687
688   if (lh < ll) {
689      a01[lh] = a[lh];
690      b01[lh] = b[lh];
691   }
692
693   KarMul(h, a01, b01, ll, stk);
694
695   for (i = 0; i < ll2; i++)
696      h[i] ^= c[i];
697
698   for (i = 0; i < lh2; i++)
699      h[i] ^= c[i+ll2];
700
701   for (i = 0; i < ll2; i++)
702      c[i+ll] ^= h[i];
703}
704
705
706
707#ifdef NTL_GF2X_LIB
708
709
710void mul(GF2X& c, const GF2X& a, const GF2X& b)
711{
712   long sa = a.xrep.length();
713   long sb = b.xrep.length();
714
715   if (sa <= 0 || sb <= 0) {
716      clear(c);
717      return;
718   }
719 
720   _ntl_ulong a0 = a.xrep[0];
721   _ntl_ulong b0 = b.xrep[0];
722
723   if (sb == 1 && b0 == 1) {
724      c = a;
725      return;
726   }
727
728   if (sa == 1 && a0 == 1) {
729      c = b;
730      return;
731   }
732
733   if (&a == &b) {
734      sqr(c, a);
735      return;
736   }
737
738   if (sa == 1 && sb == 1) {
739      // special case...
740      _ntl_ulong v[2];
741      if (!(a0 >> NTL_BITS_PER_LONG/2))
742         mul_half(v, b0, a0);
743      else if (!(b0 >> NTL_BITS_PER_LONG/2))
744         mul_half(v, a0, b0);
745      else
746         mul1(v, a0, b0);
747
748      if (v[1]) {
749         c.xrep.SetLength(2);
750         _ntl_ulong *cp = &c.xrep[0];
751         cp[0] = v[0];
752         cp[1] = v[1];
753      }
754      else {
755         c.xrep.SetLength(1);
756         _ntl_ulong *cp = &c.xrep[0];
757         cp[0] = v[0];
758      }
759      return;
760   }
761
762   if (sa == 2 && sb == 2) {
763      // special case...
764      _ntl_ulong v[4];
765      mul2(v, &a.xrep[0], &b.xrep[0]);
766      if (v[3]) {
767         c.xrep.SetLength(4);
768         _ntl_ulong *cp = &c.xrep[0];
769         cp[0] = v[0];
770         cp[1] = v[1];
771         cp[2] = v[2];
772         cp[3] = v[3];
773      }
774      else {
775         c.xrep.SetLength(3);
776         _ntl_ulong *cp = &c.xrep[0];
777         cp[0] = v[0];
778         cp[1] = v[1];
779         cp[2] = v[2];
780      }
781      return;
782   }
783
784   // another special case:  one of the two inputs
785   // has length 1 (or less).
786
787   if (sa == 1) {
788      c.xrep.SetLength(sb + 1);
789      _ntl_ulong *cp = c.xrep.elts();
790      const _ntl_ulong *bp = b.xrep.elts();
791
792      if (a0 >> (NTL_BITS_PER_LONG-NTL_BB_MUL1_BITS+1))
793         Mul1(cp, bp, sb, a0);
794      else
795         Mul1_short(cp, bp, sb, a0);
796
797
798      c.normalize();
799      return;
800   }
801
802   if (sb == 1) {
803      c.xrep.SetLength(sa + 1);
804      _ntl_ulong *cp = c.xrep.elts();
805      const _ntl_ulong *ap = a.xrep.elts();
806
807
808      if (b0 >> (NTL_BITS_PER_LONG-NTL_BB_MUL1_BITS+1))
809         Mul1(cp, ap, sa, b0);
810      else
811         Mul1_short(cp, ap, sa, b0);
812
813      c.normalize();
814      return;
815   }
816
817   // finally: the general case
818
819
820   WordVector mem;
821
822   const _ntl_ulong *ap = a.xrep.elts(), *bp = b.xrep.elts();
823   _ntl_ulong *cp;
824
825
826   long sc = sa + sb;
827   long in_mem = 0;
828
829   if (&a == &c || &b == &c) {
830      mem.SetLength(sc);
831      cp = mem.elts();
832      in_mem = 1;
833   }
834   else {
835      c.xrep.SetLength(sc);
836      cp = c.xrep.elts();
837   }
838
839   gf2x_mul(cp, ap, sa, bp, sb);
840
841   if (in_mem) {
842      c.xrep = mem;
843   }
844
845   c.normalize();
846   mem.release();
847}
848#else
849void OldMul(GF2X& c, const GF2X& a, const GF2X& b)
850{
851   mul(c, a, b);
852}
853#endif
854
855
856
857#ifdef NTL_GF2X_LIB
858void OldMul(GF2X& c, const GF2X& a, const GF2X& b)
859#else
860void mul(GF2X& c, const GF2X& a, const GF2X& b)
861#endif
862{
863   long sa = a.xrep.length();
864   long sb = b.xrep.length();
865
866   if (sa <= 0 || sb <= 0) {
867      clear(c);
868      return;
869   }
870 
871   _ntl_ulong a0 = a.xrep[0];
872   _ntl_ulong b0 = b.xrep[0];
873
874   if (sb == 1 && b0 == 1) {
875      c = a;
876      return;
877   }
878
879   if (sa == 1 && a0 == 1) {
880      c = b;
881      return;
882   }
883
884   if (&a == &b) {
885      sqr(c, a);
886      return;
887   }
888
889   if (sa == sb && sa <= 8) {
890      // we treat these cases specially for efficiency reasons
891
892      switch (sa) {
893         case 1: {
894            _ntl_ulong v[2];
895            if (!(a0 >> NTL_BITS_PER_LONG/2))
896               mul_half(v, b0, a0);
897            else if (!(b0 >> NTL_BITS_PER_LONG/2))
898               mul_half(v, a0, b0);
899            else
900               mul1(v, a0, b0);
901
902            if (v[1]) {
903               c.xrep.SetLength(2);
904               _ntl_ulong *cp = &c.xrep[0];
905               cp[0] = v[0];
906               cp[1] = v[1];
907            }
908            else {
909               c.xrep.SetLength(1);
910               _ntl_ulong *cp = &c.xrep[0];
911               cp[0] = v[0];
912            }
913         }
914         return;
915
916         case 2: {
917            _ntl_ulong v[4];
918            mul2(v, &a.xrep[0], &b.xrep[0]);
919            if (v[3]) {
920               c.xrep.SetLength(4);
921               _ntl_ulong *cp = &c.xrep[0];
922               cp[0] = v[0];
923               cp[1] = v[1];
924               cp[2] = v[2];
925               cp[3] = v[3];
926            }
927            else {
928               c.xrep.SetLength(3);
929               _ntl_ulong *cp = &c.xrep[0];
930               cp[0] = v[0];
931               cp[1] = v[1];
932               cp[2] = v[2];
933            }
934         }
935         return;
936
937         case 3: {
938            _ntl_ulong v[6];
939            mul3(v, &a.xrep[0], &b.xrep[0]);
940            if (v[5]) {
941               c.xrep.SetLength(6);
942               _ntl_ulong *cp = &c.xrep[0];
943               cp[0] = v[0];
944               cp[1] = v[1];
945               cp[2] = v[2];
946               cp[3] = v[3];
947               cp[4] = v[4];
948               cp[5] = v[5];
949            }
950            else {
951               c.xrep.SetLength(5);
952               _ntl_ulong *cp = &c.xrep[0];
953               cp[0] = v[0];
954               cp[1] = v[1];
955               cp[2] = v[2];
956               cp[3] = v[3];
957               cp[4] = v[4];
958            }
959         }
960         return;
961
962         case 4: {
963            _ntl_ulong v[8];
964            mul4(v, &a.xrep[0], &b.xrep[0]);
965            if (v[7]) {
966               c.xrep.SetLength(8);
967               _ntl_ulong *cp = &c.xrep[0];
968               cp[0] = v[0];
969               cp[1] = v[1];
970               cp[2] = v[2];
971               cp[3] = v[3];
972               cp[4] = v[4];
973               cp[5] = v[5];
974               cp[6] = v[6];
975               cp[7] = v[7];
976            }
977            else {
978               c.xrep.SetLength(7);
979               _ntl_ulong *cp = &c.xrep[0];
980               cp[0] = v[0];
981               cp[1] = v[1];
982               cp[2] = v[2];
983               cp[3] = v[3];
984               cp[4] = v[4];
985               cp[5] = v[5];
986               cp[6] = v[6];
987            }
988         }
989         return;
990
991         case 5: {
992            _ntl_ulong v[10];
993            mul5(v, &a.xrep[0], &b.xrep[0]);
994            if (v[9]) {
995               c.xrep.SetLength(10);
996               _ntl_ulong *cp = &c.xrep[0];
997               cp[0] = v[0];
998               cp[1] = v[1];
999               cp[2] = v[2];
1000               cp[3] = v[3];
1001               cp[4] = v[4];
1002               cp[5] = v[5];
1003               cp[6] = v[6];
1004               cp[7] = v[7];
1005               cp[8] = v[8];
1006               cp[9] = v[9];
1007            }
1008            else {
1009               c.xrep.SetLength(9);
1010               _ntl_ulong *cp = &c.xrep[0];
1011               cp[0] = v[0];
1012               cp[1] = v[1];
1013               cp[2] = v[2];
1014               cp[3] = v[3];
1015               cp[4] = v[4];
1016               cp[5] = v[5];
1017               cp[6] = v[6];
1018               cp[7] = v[7];
1019               cp[8] = v[8];
1020            }
1021         }
1022         return;
1023
1024         case 6: {
1025            _ntl_ulong v[12];
1026            mul6(v, &a.xrep[0], &b.xrep[0]);
1027            if (v[11]) {
1028               c.xrep.SetLength(12);
1029               _ntl_ulong *cp = &c.xrep[0];
1030               cp[0] = v[0];
1031               cp[1] = v[1];
1032               cp[2] = v[2];
1033               cp[3] = v[3];
1034               cp[4] = v[4];
1035               cp[5] = v[5];
1036               cp[6] = v[6];
1037               cp[7] = v[7];
1038               cp[8] = v[8];
1039               cp[9] = v[9];
1040               cp[10] = v[10];
1041               cp[11] = v[11];
1042            }
1043            else {
1044               c.xrep.SetLength(11);
1045               _ntl_ulong *cp = &c.xrep[0];
1046               cp[0] = v[0];
1047               cp[1] = v[1];
1048               cp[2] = v[2];
1049               cp[3] = v[3];
1050               cp[4] = v[4];
1051               cp[5] = v[5];
1052               cp[6] = v[6];
1053               cp[7] = v[7];
1054               cp[8] = v[8];
1055               cp[9] = v[9];
1056               cp[10] = v[10];
1057            }
1058         }
1059         return; 
1060
1061         case 7: {
1062            _ntl_ulong v[14];
1063            mul7(v, &a.xrep[0], &b.xrep[0]);
1064            if (v[13]) {
1065               c.xrep.SetLength(14);
1066               _ntl_ulong *cp = &c.xrep[0];
1067               cp[0] = v[0];
1068               cp[1] = v[1];
1069               cp[2] = v[2];
1070               cp[3] = v[3];
1071               cp[4] = v[4];
1072               cp[5] = v[5];
1073               cp[6] = v[6];
1074               cp[7] = v[7];
1075               cp[8] = v[8];
1076               cp[9] = v[9];
1077               cp[10] = v[10];
1078               cp[11] = v[11];
1079               cp[12] = v[12];
1080               cp[13] = v[13];
1081            }
1082            else {
1083               c.xrep.SetLength(13);
1084               _ntl_ulong *cp = &c.xrep[0];
1085               cp[0] = v[0];
1086               cp[1] = v[1];
1087               cp[2] = v[2];
1088               cp[3] = v[3];
1089               cp[4] = v[4];
1090               cp[5] = v[5];
1091               cp[6] = v[6];
1092               cp[7] = v[7];
1093               cp[8] = v[8];
1094               cp[9] = v[9];
1095               cp[10] = v[10];
1096               cp[11] = v[11];
1097               cp[12] = v[12];
1098            }
1099         }
1100         return; 
1101
1102         case 8: {
1103            _ntl_ulong v[16];
1104            mul8(v, &a.xrep[0], &b.xrep[0]);
1105            if (v[15]) {
1106               c.xrep.SetLength(16);
1107               _ntl_ulong *cp = &c.xrep[0];
1108               cp[0] = v[0];
1109               cp[1] = v[1];
1110               cp[2] = v[2];
1111               cp[3] = v[3];
1112               cp[4] = v[4];
1113               cp[5] = v[5];
1114               cp[6] = v[6];
1115               cp[7] = v[7];
1116               cp[8] = v[8];
1117               cp[9] = v[9];
1118               cp[10] = v[10];
1119               cp[11] = v[11];
1120               cp[12] = v[12];
1121               cp[13] = v[13];
1122               cp[14] = v[14];
1123               cp[15] = v[15];
1124            }
1125            else {
1126               c.xrep.SetLength(15);
1127               _ntl_ulong *cp = &c.xrep[0];
1128               cp[0] = v[0];
1129               cp[1] = v[1];
1130               cp[2] = v[2];
1131               cp[3] = v[3];
1132               cp[4] = v[4];
1133               cp[5] = v[5];
1134               cp[6] = v[6];
1135               cp[7] = v[7];
1136               cp[8] = v[8];
1137               cp[9] = v[9];
1138               cp[10] = v[10];
1139               cp[11] = v[11];
1140               cp[12] = v[12];
1141               cp[13] = v[13];
1142               cp[14] = v[14];
1143            }
1144         }
1145         return; 
1146
1147      }
1148   }
1149
1150   // another special case:  one of the two inputs
1151   // has length 1 (or less).
1152
1153   if (sa == 1) {
1154      c.xrep.SetLength(sb + 1);
1155      _ntl_ulong *cp = c.xrep.elts();
1156      const _ntl_ulong *bp = b.xrep.elts();
1157
1158      if (a0 >> (NTL_BITS_PER_LONG-NTL_BB_MUL1_BITS+1))
1159         Mul1(cp, bp, sb, a0);
1160      else
1161         Mul1_short(cp, bp, sb, a0);
1162
1163
1164      c.normalize();
1165      return;
1166   }
1167
1168   if (sb == 1) {
1169      c.xrep.SetLength(sa + 1);
1170      _ntl_ulong *cp = c.xrep.elts();
1171      const _ntl_ulong *ap = a.xrep.elts();
1172
1173
1174      if (b0 >> (NTL_BITS_PER_LONG-NTL_BB_MUL1_BITS+1))
1175         Mul1(cp, ap, sa, b0);
1176      else
1177         Mul1_short(cp, ap, sa, b0);
1178
1179      c.normalize();
1180      return;
1181   }
1182
1183   // finally: the general case
1184
1185   
1186   static WordVector mem;
1187   static WordVector stk;
1188   static WordVector vec;
1189
1190   const _ntl_ulong *ap, *bp;
1191   _ntl_ulong *cp;
1192
1193
1194   long sc = sa + sb;
1195   long in_mem = 0;
1196
1197   if (&a == &c || &b == &c) {
1198      mem.SetLength(sc);
1199      cp = mem.elts();
1200      in_mem = 1;
1201   }
1202   else {
1203      c.xrep.SetLength(sc);
1204      cp = c.xrep.elts();
1205   }
1206
1207
1208   long n, hn, sp;
1209
1210   n = min(sa, sb);
1211   sp = 0;
1212   while (n > 8) {
1213      hn = (n+1) >> 1;
1214      sp += (hn << 2) + 3;
1215      n = hn;
1216   }
1217
1218   stk.SetLength(sp);
1219   _ntl_ulong *stk_p = stk.elts();
1220
1221   if (sa > sb) {
1222      { long t; t = sa; sa = sb; sb = t; }
1223      ap = b.xrep.elts();
1224      bp = a.xrep.elts();
1225   }
1226   else {
1227      ap = a.xrep.elts();
1228      bp = b.xrep.elts();
1229   }
1230
1231
1232   vec.SetLength(2*sa);
1233
1234   _ntl_ulong *v = vec.elts();
1235
1236   long i, j;
1237
1238   for (i = 0; i < sc; i++)
1239      cp[i] = 0;
1240
1241   do {
1242      if (sa == 0) break;
1243
1244      if (sa == 1) {
1245         AddMul1(cp, bp, sb, ap[0]);
1246         
1247         break;
1248      }
1249
1250      // general case
1251
1252      for (i = 0; i+sa <= sb; i += sa) {
1253         KarMul(v, ap, bp + i, sa, stk_p);
1254         for (j = 0; j < 2*sa; j++)
1255            cp[i+j] ^= v[j]; 
1256      }
1257
1258      { const _ntl_ulong *t; t = ap; ap = bp + i; bp = t; }
1259      { long t; t = sa; sa = sb - i; sb = t; }
1260      cp = cp + i;
1261   } while (1);
1262
1263   if (in_mem)
1264      c.xrep = mem;
1265
1266   c.normalize();
1267
1268   mem.release();
1269   stk.release();
1270   vec.release();
1271}
1272
1273
1274
1275void mul(GF2X& c, const GF2X& a, long b)
1276{
1277   if (b & 1)
1278      c = a;
1279   else
1280      clear(c);
1281}
1282
1283void mul(GF2X& c, const GF2X& a, GF2 b)
1284{
1285   if (b == 1)
1286      c = a;
1287   else
1288      clear(c);
1289}
1290
1291
1292void trunc(GF2X& x, const GF2X& a, long m)
1293{
1294   if (m < 0) Error("trunc: bad args");
1295
1296   long n = a.xrep.length();
1297   if (n == 0 || m == 0) {
1298      clear(x);
1299      return;
1300   }
1301
1302   if (&x == &a) {
1303      if (n*NTL_BITS_PER_LONG > m) {
1304         long wm = (m-1)/NTL_BITS_PER_LONG;
1305         long bm = m - NTL_BITS_PER_LONG*wm;
1306         _ntl_ulong msk;
1307         if (bm == NTL_BITS_PER_LONG)
1308            msk = ~(0UL);
1309         else
1310            msk = ((1UL << bm) - 1UL);
1311         x.xrep[wm] &= msk;
1312         x.xrep.QuickSetLength(wm+1);
1313         x.normalize();
1314      }
1315   }
1316   else if (n*NTL_BITS_PER_LONG <= m) 
1317      x = a;
1318   else {
1319      long wm = (m-1)/NTL_BITS_PER_LONG;
1320      long bm = m - NTL_BITS_PER_LONG*wm;
1321      x.xrep.SetLength(wm+1);
1322      _ntl_ulong *xp = &x.xrep[0];
1323      const _ntl_ulong *ap = &a.xrep[0];
1324      long i;
1325      for (i = 0; i < wm; i++)
1326         xp[i] = ap[i];
1327      _ntl_ulong msk;
1328      if (bm == NTL_BITS_PER_LONG)
1329         msk = ~(0UL);
1330      else
1331         msk = ((1UL << bm) - 1UL);
1332      xp[wm] = ap[wm] & msk;
1333      x.normalize();
1334   }
1335}
1336     
1337
1338/****** implementation of vec_GF2X ******/
1339
1340
1341NTL_vector_impl(GF2X,vec_GF2X)
1342
1343NTL_eq_vector_impl(GF2X,vec_GF2X)
1344
1345
1346
1347void MulByX(GF2X& x, const GF2X& a)
1348{
1349   long n = a.xrep.length();
1350   if (n == 0) {
1351      clear(x);
1352      return;
1353   }
1354
1355   if (a.xrep[n-1] & (1UL << (NTL_BITS_PER_LONG-1))) {
1356      x.xrep.SetLength(n+1);
1357      x.xrep[n] = 1;
1358   }
1359   else if (&x != &a)
1360      x.xrep.SetLength(n);
1361
1362   _ntl_ulong *xp = &x.xrep[0];
1363   const _ntl_ulong *ap = &a.xrep[0];
1364
1365   long i;
1366   for (i = n-1; i > 0; i--)
1367      xp[i] = (ap[i] << 1) | (ap[i-1] >> (NTL_BITS_PER_LONG-1));
1368
1369   xp[0] = ap[0] << 1;
1370
1371   // no need to normalize
1372}
1373
1374
1375
1376static _ntl_ulong sqrtab[256] = {
1377
13780UL, 1UL, 4UL, 5UL, 16UL, 17UL, 20UL, 21UL, 64UL, 
137965UL, 68UL, 69UL, 80UL, 81UL, 84UL, 85UL, 256UL, 
1380257UL, 260UL, 261UL, 272UL, 273UL, 276UL, 277UL, 320UL, 
1381321UL, 324UL, 325UL, 336UL, 337UL, 340UL, 341UL, 1024UL, 
13821025UL, 1028UL, 1029UL, 1040UL, 1041UL, 1044UL, 1045UL, 1088UL, 
13831089UL, 1092UL, 1093UL, 1104UL, 1105UL, 1108UL, 1109UL, 1280UL, 
13841281UL, 1284UL, 1285UL, 1296UL, 1297UL, 1300UL, 1301UL, 1344UL, 
13851345UL, 1348UL, 1349UL, 1360UL, 1361UL, 1364UL, 1365UL, 4096UL, 
13864097UL, 4100UL, 4101UL, 4112UL, 4113UL, 4116UL, 4117UL, 4160UL, 
13874161UL, 4164UL, 4165UL, 4176UL, 4177UL, 4180UL, 4181UL, 4352UL, 
13884353UL, 4356UL, 4357UL, 4368UL, 4369UL, 4372UL, 4373UL, 4416UL, 
13894417UL, 4420UL, 4421UL, 4432UL, 4433UL, 4436UL, 4437UL, 5120UL, 
13905121UL, 5124UL, 5125UL, 5136UL, 5137UL, 5140UL, 5141UL, 5184UL, 
13915185UL, 5188UL, 5189UL, 5200UL, 5201UL, 5204UL, 5205UL, 5376UL, 
13925377UL, 5380UL, 5381UL, 5392UL, 5393UL, 5396UL, 5397UL, 5440UL, 
13935441UL, 5444UL, 5445UL, 5456UL, 5457UL, 5460UL, 5461UL, 16384UL, 
139416385UL, 16388UL, 16389UL, 16400UL, 16401UL, 16404UL, 16405UL, 16448UL, 
139516449UL, 16452UL, 16453UL, 16464UL, 16465UL, 16468UL, 16469UL, 16640UL, 
139616641UL, 16644UL, 16645UL, 16656UL, 16657UL, 16660UL, 16661UL, 16704UL, 
139716705UL, 16708UL, 16709UL, 16720UL, 16721UL, 16724UL, 16725UL, 17408UL, 
139817409UL, 17412UL, 17413UL, 17424UL, 17425UL, 17428UL, 17429UL, 17472UL, 
139917473UL, 17476UL, 17477UL, 17488UL, 17489UL, 17492UL, 17493UL, 17664UL, 
140017665UL, 17668UL, 17669UL, 17680UL, 17681UL, 17684UL, 17685UL, 17728UL, 
140117729UL, 17732UL, 17733UL, 17744UL, 17745UL, 17748UL, 17749UL, 20480UL, 
140220481UL, 20484UL, 20485UL, 20496UL, 20497UL, 20500UL, 20501UL, 20544UL, 
140320545UL, 20548UL, 20549UL, 20560UL, 20561UL, 20564UL, 20565UL, 20736UL, 
140420737UL, 20740UL, 20741UL, 20752UL, 20753UL, 20756UL, 20757UL, 20800UL, 
140520801UL, 20804UL, 20805UL, 20816UL, 20817UL, 20820UL, 20821UL, 21504UL, 
140621505UL, 21508UL, 21509UL, 21520UL, 21521UL, 21524UL, 21525UL, 21568UL, 
140721569UL, 21572UL, 21573UL, 21584UL, 21585UL, 21588UL, 21589UL, 21760UL, 
140821761UL, 21764UL, 21765UL, 21776UL, 21777UL, 21780UL, 21781UL, 21824UL, 
140921825UL, 21828UL, 21829UL, 21840UL, 21841UL, 21844UL, 21845UL };
1410
1411
1412
1413
1414static inline 
1415void sqr1(_ntl_ulong *c, _ntl_ulong a)
1416{
1417   _ntl_ulong hi, lo;
1418
1419   NTL_BB_SQR_CODE
1420
1421   c[0] = lo;
1422   c[1] = hi;
1423}
1424
1425
1426
1427
1428void sqr(GF2X& c, const GF2X& a)
1429{
1430   long sa = a.xrep.length();
1431   if (sa <= 0) {
1432      clear(c);
1433      return;
1434   }
1435
1436   c.xrep.SetLength(sa << 1);
1437   _ntl_ulong *cp = c.xrep.elts();
1438   const _ntl_ulong *ap = a.xrep.elts();
1439   long i;
1440
1441   for (i = sa-1; i >= 0; i--)
1442      sqr1(cp + (i << 1), ap[i]);
1443
1444   c.normalize();
1445   return;
1446}
1447
1448
1449
1450void LeftShift(GF2X& c, const GF2X& a, long n)
1451{
1452   if (IsZero(a)) {
1453      clear(c);
1454      return;
1455   }
1456
1457   if (n == 1) {
1458      MulByX(c, a);
1459      return;
1460   }
1461
1462   if (n < 0) {
1463      if (n < -NTL_MAX_LONG) 
1464         clear(c);
1465      else
1466         RightShift(c, a, -n);
1467      return;
1468   }
1469
1470   if (NTL_OVERFLOW(n, 1, 0))
1471      Error("overflow in LeftShift");
1472
1473   if (n == 0) {
1474      c = a;
1475      return;
1476   }
1477
1478   long sa = a.xrep.length();
1479
1480   long wn = n/NTL_BITS_PER_LONG;
1481   long bn = n - wn*NTL_BITS_PER_LONG;
1482
1483   long sc = sa + wn;
1484   if (bn) sc++;
1485
1486   c.xrep.SetLength(sc);
1487
1488   _ntl_ulong *cp = c.xrep.elts();
1489   const _ntl_ulong *ap = a.xrep.elts();
1490
1491   long i;
1492
1493   if (bn == 0) {
1494      for (i = sa+wn-1; i >= wn; i--)
1495         cp[i] = ap[i-wn];
1496      for (i = wn-1; i >= 0; i--)
1497         cp[i] = 0;
1498   }
1499   else {
1500      cp[sa+wn] = ap[sa-1] >> (NTL_BITS_PER_LONG-bn);
1501      for (i = sa+wn-1; i >= wn+1; i--)
1502         cp[i] = (ap[i-wn] << bn) | (ap[i-wn-1] >> (NTL_BITS_PER_LONG-bn));
1503      cp[wn] = ap[0] << bn;
1504      for (i = wn-1; i >= 0; i--)
1505         cp[i] = 0;
1506   }
1507
1508   c.normalize();
1509}
1510
1511void ShiftAdd(GF2X& c, const GF2X& a, long n)
1512// c = c + a*X^n
1513{
1514   if (n < 0) Error("ShiftAdd: negative argument");
1515
1516   if (n == 0) {
1517      add(c, c, a);
1518      return;
1519   }
1520
1521   if (NTL_OVERFLOW(n, 1, 0))
1522      Error("overflow in ShiftAdd");
1523
1524   long sa = a.xrep.length();
1525   if (sa <= 0) {
1526      return;
1527   }
1528
1529   long sc = c.xrep.length();
1530
1531   long wn = n/NTL_BITS_PER_LONG;
1532   long bn = n - wn*NTL_BITS_PER_LONG;
1533
1534   long ss = sa + wn;
1535   if (bn) ss++;
1536
1537   if (ss > sc) 
1538      c.xrep.SetLength(ss);
1539
1540   _ntl_ulong *cp = c.xrep.elts();
1541   const _ntl_ulong *ap = a.xrep.elts();
1542
1543   long i;
1544
1545   for (i = sc; i < ss; i++)
1546      cp[i] = 0;
1547
1548   if (bn == 0) {
1549      for (i = sa+wn-1; i >= wn; i--)
1550         cp[i] ^= ap[i-wn];
1551   }
1552   else {
1553      cp[sa+wn] ^= ap[sa-1] >> (NTL_BITS_PER_LONG-bn);
1554      for (i = sa+wn-1; i >= wn+1; i--)
1555         cp[i] ^= (ap[i-wn] << bn) | (ap[i-wn-1] >> (NTL_BITS_PER_LONG-bn));
1556      cp[wn] ^= ap[0] << bn;
1557   }
1558
1559   c.normalize();
1560}
1561
1562
1563
1564
1565void RightShift(GF2X& c, const GF2X& a, long n)
1566{
1567   if (IsZero(a)) {
1568      clear(c);
1569      return;
1570   }
1571
1572   if (n < 0) {
1573      if (n < -NTL_MAX_LONG) Error("overflow in RightShift");
1574      LeftShift(c, a, -n);
1575      return;
1576   }
1577
1578   if (n == 0) {
1579      c = a;
1580      return;
1581   }
1582
1583   long sa = a.xrep.length();
1584   long wn = n/NTL_BITS_PER_LONG;
1585   long bn = n - wn*NTL_BITS_PER_LONG;
1586
1587   if (wn >= sa) {
1588      clear(c);
1589      return;
1590   }
1591
1592   c.xrep.SetLength(sa-wn);
1593   _ntl_ulong *cp = c.xrep.elts();
1594   const _ntl_ulong *ap = a.xrep.elts();
1595
1596   long i;
1597
1598   if (bn == 0) {
1599      for (i = 0; i < sa-wn; i++)
1600         cp[i] = ap[i+wn];
1601   }
1602   else {
1603      for (i = 0; i < sa-wn-1; i++)
1604         cp[i] = (ap[i+wn] >> bn) | (ap[i+wn+1] << (NTL_BITS_PER_LONG - bn));
1605
1606      cp[sa-wn-1] = ap[sa-1] >> bn;
1607   }
1608
1609   c.normalize();
1610}
1611
1612
1613
1614static _ntl_ulong revtab[256] = {
1615
16160UL, 128UL, 64UL, 192UL, 32UL, 160UL, 96UL, 224UL, 16UL, 144UL, 
161780UL, 208UL, 48UL, 176UL, 112UL, 240UL, 8UL, 136UL, 72UL, 200UL, 
161840UL, 168UL, 104UL, 232UL, 24UL, 152UL, 88UL, 216UL, 56UL, 184UL, 
1619120UL, 248UL, 4UL, 132UL, 68UL, 196UL, 36UL, 164UL, 100UL, 228UL, 
162020UL, 148UL, 84UL, 212UL, 52UL, 180UL, 116UL, 244UL, 12UL, 140UL, 
162176UL, 204UL, 44UL, 172UL, 108UL, 236UL, 28UL, 156UL, 92UL, 220UL, 
162260UL, 188UL, 124UL, 252UL, 2UL, 130UL, 66UL, 194UL, 34UL, 162UL, 
162398UL, 226UL, 18UL, 146UL, 82UL, 210UL, 50UL, 178UL, 114UL, 242UL, 
162410UL, 138UL, 74UL, 202UL, 42UL, 170UL, 106UL, 234UL, 26UL, 154UL, 
162590UL, 218UL, 58UL, 186UL, 122UL, 250UL, 6UL, 134UL, 70UL, 198UL, 
162638UL, 166UL, 102UL, 230UL, 22UL, 150UL, 86UL, 214UL, 54UL, 182UL, 
1627118UL, 246UL, 14UL, 142UL, 78UL, 206UL, 46UL, 174UL, 110UL, 238UL, 
162830UL, 158UL, 94UL, 222UL, 62UL, 190UL, 126UL, 254UL, 1UL, 129UL, 
162965UL, 193UL, 33UL, 161UL, 97UL, 225UL, 17UL, 145UL, 81UL, 209UL, 
163049UL, 177UL, 113UL, 241UL, 9UL, 137UL, 73UL, 201UL, 41UL, 169UL, 
1631105UL, 233UL, 25UL, 153UL, 89UL, 217UL, 57UL, 185UL, 121UL, 249UL, 
16325UL, 133UL, 69UL, 197UL, 37UL, 165UL, 101UL, 229UL, 21UL, 149UL, 
163385UL, 213UL, 53UL, 181UL, 117UL, 245UL, 13UL, 141UL, 77UL, 205UL, 
163445UL, 173UL, 109UL, 237UL, 29UL, 157UL, 93UL, 221UL, 61UL, 189UL, 
1635125UL, 253UL, 3UL, 131UL, 67UL, 195UL, 35UL, 163UL, 99UL, 227UL, 
163619UL, 147UL, 83UL, 211UL, 51UL, 179UL, 115UL, 243UL, 11UL, 139UL, 
163775UL, 203UL, 43UL, 171UL, 107UL, 235UL, 27UL, 155UL, 91UL, 219UL, 
163859UL, 187UL, 123UL, 251UL, 7UL, 135UL, 71UL, 199UL, 39UL, 167UL, 
1639103UL, 231UL, 23UL, 151UL, 87UL, 215UL, 55UL, 183UL, 119UL, 247UL, 
164015UL, 143UL, 79UL, 207UL, 47UL, 175UL, 111UL, 239UL, 31UL, 159UL, 
164195UL, 223UL, 63UL, 191UL, 127UL, 255UL  }; 
1642
1643static inline 
1644_ntl_ulong rev1(_ntl_ulong a)
1645{
1646   return NTL_BB_REV_CODE;
1647}
1648
1649
1650
1651void CopyReverse(GF2X& c, const GF2X& a, long hi)
1652// c[0..hi] = reverse(a[0..hi]), with zero fill as necessary
1653// input may alias output
1654
1655{
1656   if (hi < 0) { clear(c); return; }
1657
1658   if (NTL_OVERFLOW(hi, 1, 0))
1659      Error("overflow in CopyReverse");
1660
1661   long n = hi+1;
1662   long sa = a.xrep.length();
1663   if (n <= 0 || sa <= 0) {
1664      clear(c);
1665      return;
1666   }
1667
1668   long wn = n/NTL_BITS_PER_LONG;
1669   long bn = n - wn*NTL_BITS_PER_LONG;
1670
1671   if (bn != 0) {
1672      wn++;
1673      bn = NTL_BITS_PER_LONG - bn;
1674   }
1675
1676   c.xrep.SetLength(wn);
1677 
1678   _ntl_ulong *cp = c.xrep.elts();
1679   const _ntl_ulong *ap = a.xrep.elts();
1680
1681   long mm = min(sa, wn);
1682   long i;
1683
1684   for (i = 0; i < mm; i++)
1685      cp[i] = ap[i];
1686
1687   for (i = mm; i < wn; i++)
1688      cp[i] = 0;
1689
1690   if (bn != 0) {
1691      for (i = wn-1; i >= 1; i--)
1692         cp[i] = (cp[i] << bn) | (cp[i-1] >> (NTL_BITS_PER_LONG-bn));
1693      cp[0] = cp[0] << bn;
1694   }
1695
1696   for (i = 0; i < wn/2; i++) {
1697      _ntl_ulong t; t = cp[i]; cp[i] = cp[wn-1-i]; cp[wn-1-i] = t;
1698   }
1699
1700   for (i = 0; i < wn; i++)
1701      cp[i] = rev1(cp[i]);
1702
1703   c.normalize();
1704}
1705
1706
1707void div(GF2X& q, const GF2X& a, GF2 b)
1708{
1709   if (b == 0)
1710      Error("div: division by zero");
1711
1712   q = a;
1713}
1714
1715void div(GF2X& q, const GF2X& a, long b)
1716{
1717   if ((b & 1) == 0)
1718      Error("div: division by zero");
1719
1720   q = a;
1721}
1722
1723
1724
1725void GF2XFromBytes(GF2X& x, const unsigned char *p, long n)
1726{
1727   if (n <= 0) {
1728      x = 0;
1729      return;
1730   }
1731
1732   const long BytesPerLong = NTL_BITS_PER_LONG/8;
1733
1734   long lw, r, i, j;
1735
1736   lw = n/BytesPerLong;
1737   r = n - lw*BytesPerLong;
1738
1739   if (r != 0) 
1740      lw++;
1741   else
1742      r = BytesPerLong;
1743
1744   x.xrep.SetLength(lw);
1745   unsigned long *xp = x.xrep.elts();
1746
1747   for (i = 0; i < lw-1; i++) {
1748      unsigned long t = 0;
1749      for (j = 0; j < BytesPerLong; j++) {
1750         t >>= 8;
1751         t += (((unsigned long)(*p)) & 255UL) << ((BytesPerLong-1)*8);
1752         p++;
1753      }
1754      xp[i] = t;
1755   }
1756
1757   unsigned long t = 0;
1758   for (j = 0; j < r; j++) {
1759      t >>= 8;
1760      t += (((unsigned long)(*p)) & 255UL) << ((BytesPerLong-1)*8);
1761      p++;
1762   }
1763
1764   t >>= (BytesPerLong-r)*8;
1765   xp[lw-1] = t;
1766
1767   x.normalize();
1768}
1769
1770void BytesFromGF2X(unsigned char *p, const GF2X& a, long n)
1771{
1772   if (n < 0) n = 0;
1773
1774   const long BytesPerLong = NTL_BITS_PER_LONG/8;
1775
1776   long lbits = deg(a) + 1;
1777   long lbytes = (lbits+7)/8;
1778
1779   long min_bytes = min(lbytes, n);
1780
1781   long min_words = min_bytes/BytesPerLong;
1782   long r = min_bytes - min_words*BytesPerLong;
1783   if (r != 0)
1784      min_words++;
1785   else
1786      r = BytesPerLong;
1787
1788   const unsigned long *ap = a.xrep.elts();
1789
1790   long i, j;
1791
1792   for (i = 0; i < min_words-1; i++) {
1793      unsigned long t = ap[i];
1794      for (j = 0; j < BytesPerLong; j++) {
1795         *p = t & 255UL;
1796         t >>= 8;
1797         p++;
1798      }
1799   }
1800
1801   if (min_words > 0) {
1802      unsigned long t = ap[min_words-1];
1803      for (j = 0; j < r; j++) {
1804         *p = t & 255UL;
1805         t >>= 8;
1806         p++;
1807      }
1808   }
1809
1810   for (j = min_bytes; j < n; j++) {
1811      *p = 0;
1812      p++;
1813   }
1814}
1815
1816NTL_END_IMPL
Note: See TracBrowser for help on using the repository browser.