1 #if !defined HAVE_BITPOLMOD_ARITH_H__
2 #define HAVE_BITPOLMOD_ARITH_H__
3 // This file is part of the FXT library.
4 // Copyright (C) 2010, 2011, 2012 Joerg Arndt
5 // License: GNU General Public License version 3 or later,
6 // see the file COPYING.txt in the main directory.
7
8 #include "bpol/bitpol-gcd.h" // bitpol_egcd()
9
10 #include "fxttypes.h"
11 #include "bits/bitsperlong.h"
12 #include "bits/bithigh.h" // highest_one() for powering
13 //#include "cmath.h"
14
15 //<<
16 // Modular arithmetic with polynomials over Z/2Z.
17 // The polynomial W represented by the word w is
18 // W = pol(w) =: \sum_k{ [bit_k(w)] * x^k}
19 // The modulus is C and
20 // h needs to be a mask with one bit set:
21 // h == highest_one(c) >> 1 == 1UL << (degree(C)-1)
22 //>>
23
24
25 #define MULT_UNROLL // define to unroll loops 4-fold
26 #define SQUARE_BY_MULT // whether to square via multiplication
27
28
bitpolmod_times_x(ulong a,ulong c,ulong h)29 static inline ulong bitpolmod_times_x(ulong a, ulong c, ulong h)
30 // Return (A * x) mod C
31 //
32 // If C is a primitive polynomial of degree n
33 // successive calls will cycle through all 2**n-1
34 // n-bit words and the sequence of bits
35 // (any fixed position) of a constitutes
36 // a shift register sequence (SRS).
37 // Start with a=2 to get an SRS that starts with
38 // n-1 consecutive zeros (use bit 0 of a)
39 {
40 ulong s = a & h;
41 a <<= 1;
42 if ( s ) a ^= c;
43 return a;
44 }
45 // -------------------------
46
bitpolmod_times_x2(ulong a,ulong c,ulong h)47 static inline ulong bitpolmod_times_x2(ulong a, ulong c, ulong h)
48 // Return (A * x * x) mod C
49 // Used for squaring.
50 {
51 { ulong s=a&h; a<<=1; if (s) a^=c; }
52 { ulong s=a&h; a<<=1; if (s) a^=c; }
53 return a;
54 }
55 // -------------------------
56
bitpolmod_div_x(ulong a,ulong c,ulong h)57 static inline ulong bitpolmod_div_x(ulong a, ulong c, ulong h)
58 // Return (A / x) mod C
59 // C must have nonzero constant term: (c&1)==1
60 //
61 // If C is a primitive polynomial of degree n
62 // successive calls will cycle through all 2**n-1
63 // n-bit words and the sequence of bits
64 // (any fixed position) of a constitutes
65 // a shift register sequence (SRS).
66 {
67 ulong s = a & 1;
68 a >>= 1;
69 if ( s )
70 {
71 a ^= (c>>1);
72 a |= h; // so it also works for n == BITS_PER_LONG
73 }
74 return a;
75 }
76 // -------------------------
77
78
bitpolmod_inv_x(ulong c,ulong h)79 static inline ulong bitpolmod_inv_x(ulong c, ulong h)
80 // Return (1 / x) mod C
81 // C must have nonzero constant term: (c&1)==1
82 // c>>1 == (C-1)/x = (0-1)/x == -1/x == 1/x (mod C).
83 {
84 // commented out lines take care of cases when x has no inverse mod C:
85 // if ( 0==(c & 1) ) return 0;
86 // if ( 1==c ) return 0;
87
88 ulong a = (c>>1);
89 a |= h; // so it also works for n == BITS_PER_LONG
90 return a;
91 }
92 // -------------------------
93
94
bitpolmod_mult(ulong a,ulong b,ulong c,ulong h)95 static inline ulong bitpolmod_mult(ulong a, ulong b, ulong c, ulong h)
96 // Return (A * B) mod C
97 //
98 // Must have deg(A) < deg(C) and deg(B) < deg(C)
99 //.
100 // With b=2 (== 'x') the result is identical to
101 // bitpolmod_times_x(a, c, h)
102 {
103 ulong t = 0;
104 do
105 {
106 #ifndef MULT_UNROLL
107 if ( b & 1 ) t ^= a;
108 b >>= 1;
109
110 ulong s = a & h;
111 a <<= 1;
112 if ( s ) a ^= c;
113
114 #else // unrolled:
115 { if (b&1) t^=a; b>>=1; ulong s=a&h; a<<=1; if (s) a^=c; }
116 { if (b&1) t^=a; b>>=1; ulong s=a&h; a<<=1; if (s) a^=c; }
117 { if (b&1) t^=a; b>>=1; ulong s=a&h; a<<=1; if (s) a^=c; }
118 { if (b&1) t^=a; b>>=1; ulong s=a&h; a<<=1; if (s) a^=c; }
119 #endif
120 }
121 while ( b );
122 return t;
123 }
124 // -------------------------
125
bitpolmod_square(ulong a,ulong c,ulong h)126 static inline ulong bitpolmod_square(ulong a, ulong c, ulong h)
127 // Return A*A mod C
128 // == bitpolmod_mult(a, a, c, h);
129 // compute \sum_{k=0}^{d}{a_k\,x^k} as \sum_{k=0}^{d}{a_k\,x^{2k}}
130 {
131 #ifdef SQUARE_BY_MULT
132
133 return bitpolmod_mult(a, a, c, h);
134
135 #else // SQUARE_BY_MULT
136
137 ulong t = 0, s = 1;
138 do
139 {
140 #ifndef MULT_UNROLL
141 if ( a&1 ) t ^= s;
142 a >>= 1;
143 s = bitpolmod_times_x2(s, c, h);
144 #else // unrolled:
145 if (a&1) t^=s; a>>=1; s=bitpolmod_times_x2(s, c, h);
146 if (a&1) t^=s; a>>=1; s=bitpolmod_times_x2(s, c, h);
147 if (a&1) t^=s; a>>=1; s=bitpolmod_times_x2(s, c, h);
148 if (a&1) t^=s; a>>=1; s=bitpolmod_times_x2(s, c, h);
149 #endif // MULT_UNROLL
150 }
151 while ( a );
152 return t;
153 #endif // SQUARE_BY_MULT
154 }
155 // -------------------------
156
157
bitpolmod_power(ulong a,ulong e,ulong c,ulong h)158 inline ulong bitpolmod_power(ulong a, ulong e, ulong c, ulong h)
159 // Return (A ** e) mod C
160 {
161 #if 1 // left-to-right powering:
162 ulong s = a;
163 ulong b = highest_one(e);
164 while ( b>1 )
165 {
166 b >>= 1;
167 s = bitpolmod_square(s, c, h); // s *= s;
168 if ( e & b ) s = bitpolmod_mult(s, a, c, h); // s *= a;
169 }
170 return s;
171
172 #else // right-to-left powering:
173
174 if ( 0==e ) return 1; // avoid hang with e==0 in following while
175 ulong s = a;
176 while ( 0==(e&1) )
177 {
178 s = bitpolmod_square(s, c, h);
179 e >>= 1;
180 }
181
182 a = s;
183 while ( 0!=(e>>=1) )
184 {
185 s = bitpolmod_square(s, c, h);
186 if ( e & 1 ) a = bitpolmod_mult(a, s, c, h);
187 }
188 return a;
189
190 #endif
191 }
192 // -------------------------
193
bitpolmod_xpower(ulong e,ulong c,ulong h)194 inline ulong bitpolmod_xpower(ulong e, ulong c, ulong h)
195 // Return (x ** e) mod C
196 {
197 ulong s = 2; // 'x'
198 ulong b = highest_one(e);
199 while ( b>1 )
200 {
201 b >>= 1;
202 s = bitpolmod_square(s, c, h); // s *= s;
203 if ( e & b ) s = bitpolmod_times_x(s, c, h); // s *= x;
204 }
205 return s;
206 }
207 // -------------------------
208
bitpolmod_inverse_irred(ulong a,ulong c,ulong h)209 inline ulong bitpolmod_inverse_irred(ulong a, ulong c, ulong h)
210 // Return (A ** -1) mod C
211 // Must have: C irreducible.
212 //
213 // With irreducible C the inverse of A can be obtained via
214 // i = bitpolmod_power(a, r1, c, h)
215 // where r1 = (h<<1)-2 = max_order - 1 = 2^degree(C) - 2
216 // Then 1 == bitpolmod_mult(a, i, c, h)
217 {
218 ulong r1 = (h<<1) - 2; // max order minus one
219 ulong i = bitpolmod_power(a, r1, c, h);
220 return i;
221 }
222 // -------------------------
223
bitpolmod_inverse(ulong a,ulong c)224 inline ulong bitpolmod_inverse(ulong a, ulong c)
225 // Return the inverse of A modulo C if it exists, else zero.
226 // Must have deg(A) < deg(C)
227 {
228 ulong i, t; // t unused
229 ulong g = bitpol_egcd(a, c, i, t);
230 if ( g!=1 ) i = 0;
231 return i;
232 }
233 // -------------------------
234
bitpolmod_divide(ulong a,ulong b,ulong c,ulong h)235 inline ulong bitpolmod_divide(ulong a, ulong b, ulong c, ulong h)
236 // Return A/B modulo C.
237 // Must have: gcd(b,c)==1
238 {
239 ulong i = bitpolmod_inverse(b, c);
240 a = bitpolmod_mult(a, i, c, h);
241 return a;
242 }
243 // -------------------------
244
245
bitpolmod_sqrt(ulong a,ulong c,ulong h)246 static inline ulong bitpolmod_sqrt(ulong a, ulong c, ulong h)
247 // Return sqrt(A) mod C
248 // Using the identity sqrt(A) = A^(2^(n-1)) where n=deg(C)
249 {
250 ulong d = h;
251 while ( (d>>=1) ) a = bitpolmod_square(a, c, h);
252 return a;
253 }
254 // -------------------------
255
256
257 #undef MULT_UNROLL
258 #undef SQUARE_BY_MULT
259
260 #endif // !defined HAVE_BITPOLMOD_ARITH_H__
261