Arndt - Algorithms for Programmers (523138), страница 45
Текст из файла (страница 45)
SHIFT REGISTER SEQUENCES{}}else264a ^= t;q ^= qb;do{t >>= 1;h >>= 1;qb >>= 1;}while ( h > a );}while ( t>=b );return q;return0;The polynomial greatest common divisor:inline ulong bitpol_gcd(ulong a, ulong b)// Return polynomial gcd(A, B){if ( 0==a ) return 0;if ( 0==b ) return 0;while ( 0!=b ){ulong c = bitpol_rem(a, b);a = b;b = c;}return a;}12.3.2Computations modulo a polynomialHere we consider arithmetic modulo a binary polynomial. The functions can be found in [FXT: fileauxbit/bitpolmodmult.h].Multiplication by x modulo (a polynomial) C is achieved by shifting left and subtracting (that is: XOR)C if the coefficient shifted out is one:static inline ulong bitpolmod_times_x(ulong a, ulong c, ulong h)// Return (A * x) mod C// where A and C represent polynomials over Z/2Z:// W = pol(w) =: \sum_k{ [bit_k(w)] * x^k}//// h needs to be a mask with one bit set:// h == highest_bit(c) >> 1 == 1UL << (degree(C)-1)//// If C is a primitive polynomial of degree n// successive calls will cycle through all 2**n-1// n-bit words and the sequence of bits// (any fixed position) of a constitutes// a shift register sequence (SRS).// Start with a=2 to get a SRS that starts with// n-1 consecutive zeroes (use bit 0 of a){ulong s = a & h;a <<= 1;if ( s ) a ^= c;return a;}Multiplication also needs to take care of the modulus:inline ulong bitpolmod_mult(ulong a, ulong b, ulong c, ulong h)// Return (A * B) mod C//// With b=2 (== ’x’) the result is identical to//bitpolmod_times_x(a, c){CHAPTER 12.
SHIFT REGISTER SEQUENCES}265ulong t = 0;while ( b ){if ( b & 1 ) t ^= a;b >>= 1;ulong s = a & h;a <<= 1;if ( s ) a ^= c;}return t;Now exponentiation is identical to bitpol_power() except for the call to the modulo-multiplicationfunctions:inline ulong bitpolmod_power(ulong a, ulong x, ulong c, ulong h)// Return (A ** x) mod C//// With primitive C the inverse of A can be obtained via// i = bitpolmod_power(a, c, r1, h)// where r1 = (h<<1)-2 = max_order - 1 = 2^degree(C) - 2// Then 1 == bitpolmod_mult(a, c, i, h){if ( 0==x ) return 1;ulong s = a;while ( 0==(x&1) ){s = bitpolmod_square(s, c, h);x >>= 1;}a = s;while ( 0!=(x>>=1) ){s = bitpolmod_square(s, c, h);if ( x & 1 ) a = bitpolmod_mult(a, s, c, h);}return a;}With primitive C the inverse of a polynomial A can be obtained by raising A to a power that equals themaximal order minus one:r1 = (h<<1)-2; // = max_order - 1 = 2^degree(C) - 2i = bitpolmod_power(a, c, r1, h);// here: 1==bitpolmod_mult(a, i, c, h);12.3.3Testing for irreducibilityThe following functions can be found in [FXT: file auxbit/bitpolirred.h].The derivative of a binary polynomial can be computed likeinline ulong bitpol_deriv(ulong x)// Return derived polynomial{#if BITS_PER_LONG >= 64x &= 0xaaaaaaaaaaaaaaaaUL;#elsex &= 0xaaaaaaaaUL;#endifreturn (x>>1);}The coefficients at the even powers have to be deleted because derivation multiplies them with an evenfactor (which is zero modulo two).Now extraction of squared factors can be achieved viainline ulong bitpol_test_squarefree(ulong x)CHAPTER 12.
SHIFT REGISTER SEQUENCES266// Return 0 if polynomial is square-free// else return square factor != 0{ulong d = bitpol_deriv(x);if ( 0==d ) return (1==x ? 0 : x);ulong g = bitpol_gcd(x, d);return (1==g ? 0 : g);}Testing for irreducibility (that is, whether the polynomial has no non-trivial factors) uses the fact thatkthe polynomial x2 + x has all irreducible polynomials of degree k as a factor. For example,5x2 + x ==x32 + xx · (x + 1) ·¡¢ ¡¢ ¡¢· x5 + x2 + 1 · x5 + x3 + 1 · x5 + x3 + x2 + x + 1 ·¡¢ ¡¢ ¡¢· x5 + x4 + x2 + x + 1 · x5 + x4 + x3 + x + 1 · x5 + x4 + x3 + x2 + 1(12.1)kFor a degree-d polynomial C one can compute uk = x2 (modulo C) for each k < d by successive squaringsand test whether gcd(C, uk + x) 6= 1 for all k.
But as a factor of degree f implies another one of degreed − f it suffices to do the first bd/2c of the tests.inline ulong bitpol_irreducible_q(ulong c, ulong h)// Return zero if C is reducible// else (i.e. C is irreducible) return value != 0// C must not be zero.//// h needs to be a mask with one bit set:// h == highest_bit(c) >> 1 == 1UL << (degree(C)-1){// if ( 0==(1&c) ) return (c==2 ? 1 : 0); // x is a factor// if ( 0==(c & 0xaa..a ) ) return 0; // at least one odd degree term// if ( 0==parity(c) ) return 0; // need odd number of nonzero coeff.// if ( 0!=bitpol_test_squarefree(c) ) return 0; // must be square freeulong d = c;ulong u = 2; // =^= xwhile ( 0 != (d>>=2) ) // floor( degree/2 ) times{// Square r-times for coefficients of c in GF(2^r).// We have r==1u = bitpolmod_mult(u, u, c, h);ulong upx = u ^ 2; // =^= u+xulong g = bitpol_gcd(upx, c);if ( 1!=g ) return 0;// reducible}return 1; // irreducible}Letzs=dx2 + x= 1+d−1X(12.2)x2k(12.3)k=0t=s−1= z/s(12.4)Then s is has all degree-d irreducible polynomials of trace 1 as factors and t those with trace 0.
As anexample consider the case d = 7. We use the notation [a, b, c, . . .] := xa + xb + xc + . . .:z = [128,1]s = [64,32,16,8,4,2,1,0] =[1,0] *[7,6,0] *[7,6,3,1,0] *CHAPTER 12. SHIFT REGISTER SEQUENCES267[7,6,4,1,0] *[7,6,4,2,0] *[7,6,5,2,0] *[7,6,5,3,2,1,0] *[7,6,5,4,0] *[7,6,5,4,2,1,0] *[7,6,5,4,3,2,0]t = [64,32,16,8,4,2,1] =[1] *[7,1,0] *[7,3,0] *[7,3,2,1,0] *[7,4,0] *[7,4,3,2,0] *[7,5,2,1,0] *[7,5,3,1,0] *[7,5,4,3,0] *[7,5,4,3,2,1,0]inline ulong bitpol_compose_xp1(ulong c)// Return C(x+1)// self-inverse{ulong z = 1;ulong r = 0;while ( c ){if ( c & 1 ) r ^= z;c >>= 1;z ^= (z<<1);}return r;}A version that avoids any branches and finishes in time log2 (b) (where b = bits per word) isinline ulong bitpol_compose_xp1(ulong c)// Return C(x+1){ulong s = BITS_PER_LONG >> 1;ulong m = ~0UL << s;while ( s ){c ^= ( (c&m) >> s );s >>= 1;m ^= (m>>s);}return c;}Which is exactly the blue_code() from section 8.23.When a polynomial is irreducible then the composition with x+1 is also irreducible.
Similar, the reversedword corresponds to another irreducible polynomial. The two statements remain true if ‘irreducible’ isreplaced by ‘primitive’ (see section 12.2.1).inline ulong bitpol_recip(ulong c)// Return x^deg(C) * C(1/x) (the reciprocal polynomial){ulong t = 0;while ( c ){t <<= 1;t |= (c & 1);c >>= 1;}return t;}In general the sequence of successive ‘compose’ and ‘reverse’ operations leads to 6 different polynomials:C= [11, 10, 4, 3, 0]CHAPTER 12.
SHIFT REGISTER SEQUENCES[11,[11,[11,[11,[11,[11,[11,10,8,10,10,9,9,10,4,7,9,7,7,4,4,3,1,7,6,2,2,3,2680] -- recip (C=bitpol_recip(C)) -->0]-- compose (C=bitpol_compose_xp1(C)) -->6, 5, 4, 1, 0] -- recip -->5, 4, 2, 1, 0] -- compose -->0] -- recip -->0] -- compose -->0] == initial valueTBD: factorization12.3.4Modulo multiplication with reversed polynomials *For the generation of shift register sequences (see section 12.2) a cycle-saving variant of the multiplicationby x can be found in [FXT: file auxbit/bitpolmodmultrev.h]:static inline ulong bitpolmod_times_x_rev(ulong a, ulong c)// Return (A * reverse(x)) mod C// where A and C represent polynomials over Z/2Z://W = pol(w) =: \sum_k{ [bit_k(w)] * x^k}//// If c is a primitive polynomial of degree n// successive calls will cycle through all 2**n-1// n-bit words and the sequence of bits// (any fixed position) of a constitutes// a shift register sequence (SRS).// Start with a=(1UL<<(n-1)) to get a SRS// that starts with n-1 consecutive zeroes{if ( a & 1 ) a ^= c;a >>= 1;return a;}For primitive polynomials see section 12.2.1.Similarly, a multiplication by the a reversed polynomial can be implemented as:inline ulong bitpolmod_mult_rev(ulong a, ulong b, ulong c)// Return (A * reverse(B)) mod C// where A, B and C represent polynomials over Z/2Z://W = pol(w) =: \sum_k{ [bit_k(w)] * x^k}//// With b=2 (== ’x’) the result is identical to//bitpolmod_times_x_rev(a, c){ulong t = 0;while ( b ){if ( b & 1 ) t ^= a;b >>= 1;if ( a & 1 ) a ^= c;a >>= 1;}return t;}12.4Feedback carry shift register (FCSR)There is a nice analogue of the LFSR in the modulo world, the feedback carry shift register (FCSR).
Withthe LFSR we needed an irreducible (‘prime’) polynomial C where x has maximal order. The powers ofx modulo C did run through all different (non-zero) words. Now take a prime c where 2 has maximalorder (that is, 2 is s primitive root modulo c). Then the powers of 2 modulo c run through all non-zerovalues less than c.The crucial part of the implementation isCHAPTER 12. SHIFT REGISTER SEQUENCES269ulong next(){a_ <<= 1;if ( a_ > c_ ) a_ -= c_;w_ <<= 1;ulong s = a_ & 1;w_ |= s;w_ &= mask_;a_ &= mask_;return w_;}[FXT: class fcsr in auxbit/fcsr.h] The routine is much simpler than the approach usually found inthe literature.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
