source: git/ntl/src/GF2X.c @ 2cfffe

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