1 #if !defined HAVE_BITPOL_ARITH_H__
2 #define HAVE_BITPOL_ARITH_H__
3 // This file is part of the FXT library.
4 // Copyright (C) 2010, 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 "bits/bithigh.h" // highest_one()
9 //#include "bits/revgraycode.h" // inverse_rev_gray_code()
10 #include "bits/bitsperlong.h"
11 #include "fxttypes.h"
12
13
14 //<<
15 // Arithmetic with polynomials over Z/2Z
16 // The polynomial W represented by the word w is
17 // W = pol(w) =: \sum_k{ [bit_k(w)] * x^k}
18 //>>
19
20 #define MULT_UNROLL // define to unroll loops 4-fold
21
bitpol_mult(ulong a,ulong b)22 inline ulong bitpol_mult(ulong a, ulong b)
23 // Return A * B
24 // B=2 corresponds to multiplication with 'x'
25 // Note that the result will silently overflow
26 // if deg(A)+deg(B) >= BITS_PER_LONG
27 {
28 ulong t = 0;
29 do
30 {
31 if ( b & 1 ) t ^= a;
32 b >>= 1; a <<= 1;
33 #ifdef MULT_UNROLL
34 if ( b & 1 ) t ^= a;
35 b >>= 1; a <<= 1;
36 if ( b & 1 ) t ^= a;
37 b >>= 1; a <<= 1;
38 if ( b & 1 ) t ^= a;
39 b >>= 1; a <<= 1;
40 #endif
41 }
42 while ( b );
43 return t;
44 }
45 // -------------------------
46
bitpol_square(ulong a)47 inline ulong bitpol_square(ulong a)
48 // Return A * A
49 // == bit_zip(a)
50 // == bit_zip0(a)
51 // == bitpol_mult(a, a);
52 {
53 ulong t = 0, m = 1UL;
54 do
55 {
56 if ( a&1 ) t ^= m;
57 m <<= 2; a >>= 1;
58 #ifdef MULT_UNROLL
59 if ( a&1 ) t ^= m;
60 m <<= 2; a >>= 1;
61 if ( a&1 ) t ^= m;
62 m <<= 2; a >>= 1;
63 if ( a&1 ) t ^= m;
64 m <<= 2; a >>= 1;
65 #endif
66 }
67 while ( a );
68 return t;
69 }
70 // -------------------------
71
72
bitpol_power(ulong a,ulong e)73 inline ulong bitpol_power(ulong a, ulong e)
74 // Return A ** e
75 // Overflow will occur even for moderate exponents
76 {
77 if ( 0==a ) return 0;
78 if ( 0==e ) return 1;
79
80 ulong s = a;
81 while ( 0==(e&1) )
82 {
83 s = bitpol_square(s);
84 e >>= 1;
85 }
86
87 a = s;
88 while ( 0!=(e>>=1) )
89 {
90 s = bitpol_square(s);
91 if ( e & 1 ) a = bitpol_mult(a, s);
92 }
93 return a;
94 }
95 // -------------------------
96
97
bitpol_rem(ulong a,const ulong b)98 inline ulong bitpol_rem(ulong a, const ulong b)
99 // Return R = A % B = A - (A/B)*B
100 // Must have: B!=0
101 {
102 // if ( 0==b ) return a; // allow b==0
103 const ulong db = highest_one_idx(b);
104 ulong da;
105 while ( db <= (da=highest_one_idx(a)) )
106 {
107 if ( 0==a ) break; // needed because highest_one_idx(0)==highest_one_idx(1)
108 a ^= (b<<(da-db));
109 }
110 return a;
111
112 // version for polynomials of small degree:
113 // while ( b <= a )
114 // {
115 // ulong t = b;
116 // while ( (a^t) > t ) t <<= 1;
117 // // =^= while ( highest_one(a) > highest_one(t) ) t <<= 1;
118 // a ^= t;
119 // }
120 // return a;
121
122 // if ( b <= a )
123 // {
124 // ulong t = b;
125 // while ( (a^t) > t ) t <<= 1;
126 // ulong h = highest_one(t);
127 // do
128 // {
129 // a ^= t;
130 // do
131 // {
132 // t >>= 1;
133 // h >>= 1;
134 // }
135 // while ( h > a );
136 // }
137 // while ( t>=b );
138 // }
139 // return a;
140 }
141 // -------------------------
142
bitpol_div(ulong a,ulong b)143 inline ulong bitpol_div(ulong a, ulong b)
144 // Return Q = A / B
145 // Must have B!=0.
146 {
147 // if ( 0==b ) return 0; // allow b==0
148 const ulong db = highest_one_idx(b);
149 ulong q = 0; // quotient
150 ulong da;
151 while ( db <= (da=highest_one_idx(a)) )
152 {
153 if ( 0==a ) break; // needed because highest_one_idx(0)==highest_one_idx(1)
154 a ^= (b<<(da-db));
155 q ^= (1UL<<(da-db));
156 }
157 return q;
158
159 // if ( b > a ) return 0;
160 // ulong t = b;
161 // ulong qb = 1;
162 // while ( (a^t) > t ) { t <<= 1; qb <<= 1; }
163 // ulong h = highest_one(t);
164 // ulong q = 0;
165 // do
166 // {
167 // a ^= t;
168 // q ^= qb;
169 // do
170 // {
171 // t >>= 1;
172 // h >>= 1;
173 // qb >>= 1;
174 // }
175 // while ( h > a );
176 // }
177 // while ( t>=b );
178 // return q;
179 }
180 // -------------------------
181
182
bitpol_divrem(ulong a,ulong b,ulong & q,ulong & r)183 inline void bitpol_divrem(ulong a, ulong b, ulong &q, ulong &r)
184 // Set R, Q so that A == Q * B + R.
185 // Must have B!=0.
186 // Equivalent to: q=bitpol_div(a,b); r=bitpol_rem(a,b);
187 {
188 // if ( 0==b ) { r=a; q=0; return; } // allow b==0
189
190 const ulong db = highest_one_idx(b);
191 q = 0; // quotient
192 ulong da;
193 while ( db <= (da=highest_one_idx(a)) )
194 {
195 if ( 0==a ) break; // needed because highest_one_idx(0)==highest_one_idx(1)
196 a ^= (b<<(da-db));
197 q ^= (1UL<<(da-db));
198 }
199 r = a;
200 }
201 // -------------------------
202
bitpol_div_xp1(ulong a)203 inline ulong bitpol_div_xp1(ulong a)
204 // Return power series A / (x+1)
205 // If A is a multiple of x+1, then the returned value
206 // is the exact division by x+1
207 // Function identical to inverse_rev_gray_code(a)
208 {
209 a ^= a<<1; // rev_gray ** 1
210 a ^= a<<2; // rev_gray ** 2
211 a ^= a<<4; // rev_gray ** 4
212 a ^= a<<8; // rev_gray ** 8
213 a ^= a<<16; // rev_gray ** 16
214 #if BITS_PER_LONG >= 64
215 a ^= a<<32; // for 64bit words
216 #endif
217 return a;
218 }
219 // -------------------------
220
bitpol_div_x2p1(ulong a)221 inline ulong bitpol_div_x2p1(ulong a)
222 // Return power series A / (x^2+1)
223 // If A is a multiple of x^2+1, then the returned value
224 // is the exact division by x^2+1
225 {
226 a ^= a<<2; // rev_gray ** 2
227 a ^= a<<4; // rev_gray ** 4
228 a ^= a<<8; // rev_gray ** 8
229 a ^= a<<16; // rev_gray ** 16
230 #if BITS_PER_LONG >= 64
231 a ^= a<<32; // for 64bit words
232 #endif
233 return a;
234 }
235 // -------------------------
236
237
238
239 #undef MULT_UNROLL
240
241 #endif // !defined HAVE_BITPOL_ARITH_H__
242