Arndt - Algorithms for Programmers (523138), страница 33
Текст из файла (страница 33)
The argument modulo 8 has to be equal to one.inline ulong invsqrt2adic(ulong d)// Return inverse square root modulo 2**BITS_PER_LONG// must have d==1 mod 8// number of correct bits are doubled with each step// ==> loop is executed prop. log_2(BITS_PER_LONG) times// precision is 4, 8, 16, 32, 64, ... bits (or better){if ( 1 != (d&7) ) return 0; // no inverse sqrt// start value: if d == ****10001 ==> x := ****1001ulong x = (d >> 1) | 1;ulong p, y;do{y = x;p = (3 - d * y * y);x = (y * p) >> 1;}while ( x!=y );return x;}The square root can be obtained by final multiplication with d.inline ulong sqrt2adic(ulong d)// Return square root modulo 2**BITS_PER_LONG// must have d==1 mod 8 or d==4 mod 32, d==16 mod 128// ... d==4**k mod 4**(k+3)// undefined return if condition does not hold{if ( 0==d ) return 0;ulong s = 0;while ( 0==(d&1) ) { d >>= 1; ++s; }d *= invsqrt2adic(d);d <<= (s>>1);returnd;}Note that the 2-adic square root is something completely different from the integer square root.The described functions can be found in [FXT: file auxbit/bit2adic.h].
A little demo is [FXT: filedemo/bit2adic-demo.cc], where no inverse or square root is given, it does not exit:CHAPTER 8. SOME BIT WIZARDRYxinvsqrtxinvxxxinvxinvxsqrtxxinvxinvxxxinvxinvsqrtxxxinvsqrt==============================191...............................1...............................1...............................11111111111111111111111111111111111111111111111111111111111111111..............................1.1111111111111111111111111111111...............................111.1.1.1.1.1.1.1.1.1.1.1.1.1.1.11111111111111111111111111111111.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.............................1................................1.111111111111111111111111111111...............................1.111..11..11..11..11..11..11..11.111111111111111111111111111111.11..11..11..11..11..11..11..11..11.............................11.11111111111111111111111111111.1..............................1111.11.11.11.11.11.11.11.11.11.11111111111111111111111111111111..1.1..1..1..1..1..1..1..1..1..1..11..111..1..11...11......1.11.1.1............................1...11111111111111111111111111111...............................1..1..111...111...111...111...111..1111111111111111111111111111111.1=1=-1===2-23=-3=4==-45=-5===6-67=-7===8-89For further information on 2-adic (more generally p-adic) numbers see [11], [15] and also [67].8.22Powers of the Gray codeConsider the Gray code as a transform of a binary word.
The following list shows the eight vectorsof the ‘canonical basis’ in the leftmost column. The columns to the right show the powers (repeatedapplications) of the Gray code for 8-bit words:.......1......1......1......1......1......1......1......1.......g**0=id.......1......11.....11.....11.....11.....11.....11.....11......g**1=g.......1......1......1.1....1.1....1.1....1.1....1.1....1.1.....g**2.......1......11.....111....1111...1111...1111...1111...1111....g**3.......1......1......1......1......1...1..1...1..1...1..1...1...g**4.......1......11.....11.....11.....11..1..11..11.11..11.11..11..g**5.......1......1......1.1....1.1....1.1.1..1.1.1..1.1.1.11.1.1.1.g**6This motivates theinline ulong gray_pow(ulong x, ulong e)// Return (gray_code**e)(x)// gray_pow(x, 1) == gray_code(x)// gray_pow(x, BITS_PER_LONG-1) == inverse_gray_code(x){e &= (BITS_PER_LONG-1); // modulo BITS_PER_LONGulong s = 1;while ( e ){if ( e & 1 ) x ^= x >> s; // gray ** ss <<= 1;e >>= 1;}return x;}The powers of the inverse Gray code cycle through the columns in the opposite direction, so.......1......11.....111....1111...11111..111111.111111111111111g**7=g**(-1)CHAPTER 8.
SOME BIT WIZARDRY192inline ulong inverse_gray_pow(ulong x, ulong e)// Return (inverse_gray_code**(e))(x)//== (gray_code**(-e))(x)// inverse_gray_pow(x, 1) == inverse_gray_code(x)// gray_pow(x, BITS_PER_LONG-1) == gray_code(x){return gray_pow(x, -e);}Same for the green code (this time using g for the green code):.......1......1......1......1......1......1......1......1.......g**0=id......11.....11.....11.....11.....11.....11.....11......1.......g**1=g.....1.1....1.1....1.1....1.1....1.1....1.1......1......1.......g**2....1111...1111...1111...1111...1111....111.....11......1.......g**3...1...1..1...1..1...1..1...1......1......1......1......1.......g**4..11..11.11..11.11..11..1..11.....11.....11.....11......1.......g**5.1.1.1.11.1.1.1..1.1.1..1.1.1....1.1....1.1......1......1.......g**6We just have to reverse the shift operator for green ** sinline ulong green_pow(ulong x, ulong e)// Return (green_code**e)(x)// green_pow(x, 1) == green_code(x)// green_pow(x, BITS_PER_LONG-1) == inverse_green_code(x){e &= (BITS_PER_LONG-1); // modulo BITS_PER_LONGulong s = 1;while ( e ){if ( e & 1 ) x ^= x << s; // green ** ss <<= 1;e >>= 1;}return x;}inline ulong inverse_green_pow(ulong x, ulong e)// Return (inverse_green_code**(e))(x)//== (green_code**(-e))(x)// inverse_green_pow(x, 1) == inverse_green_code(x)// green_pow(x, BITS_PER_LONG-1) == green_code(x){return green_pow(x, -e);}8.23Invertible transforms on wordsThe functions presented in this section are ‘transforms’ on binary words that are invertible.Considerinline ulong blue_code(ulong a){ulong s = BITS_PER_LONG >> 1;ulong m = ~0UL << s;while ( s ){a ^= ( (a&m) >> s );s >>= 1;m ^= (m>>s);}return a;}andinline ulong yellow_code(ulong a){ulong s = BITS_PER_LONG >> 1;111111111111111.111111..11111...1111....111.....11......1.......g**7=G**(-1)CHAPTER 8.
SOME BIT WIZARDRY}193ulong m = ~0UL >> s;while ( s ){a ^= ( (a&m) << s );s >>= 1;m ^= (m<<s);}return a;[FXT: file auxbit/bittransforms.h] Both involve a computational work ∼ log2 (b) where b is the number of bits per word (BITS_PER_LONG). The blue_code can be used as a fast implementation for thecomposition of a binary polynomial with x + 1, see page 267.[FXT: file demo/bittransforms-blue-demo.cc] gives (leftmost column gives Gray code for comparison,bit-counts at the right of each column):0:1:2:3:4:5:6:7:8:9:10:11:12:13:14:15:16:17:18:19:20:21:22:23:24:25:26:27:28:29:30:31:g=......g=.....1g=....11g=....1.g=...11.g=...111g=...1.1g=...1..g=..11..g=..11.1g=..1111g=..111.g=..1.1.g=..1.11g=..1..1g=..1...g=.11...g=.11..1g=.11.11g=.11.1.g=.1111.g=.11111g=.111.1g=.111..g=.1.1..g=.1.1.1g=.1.111g=.1.11.g=.1..1.g=.1..11g=.1...1g=.1....01212321234323212343454323432321b=......b=.....1b=....11b=....1.b=...1.1b=...1..b=...11.b=...111b=..1111b=..111.b=..11..b=..11.1b=..1.1.b=..1.11b=..1..1b=..1...b=.1...1b=.1....b=.1..1.b=.1..11b=.1.1..b=.1.1.1b=.1.111b=.1.11.b=.1111.b=.11111b=.111.1b=.111..b=.11.11b=.11.1.b=.11...b=.11..101212123432323212123234345434323y=................................y=11111111111111111111111111111111y=1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.y=.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1.1y=11..11..11..11..11..11..11..11..y=..11..11..11..11..11..11..11..11y=.11..11..11..11..11..11..11..11.y=1..11..11..11..11..11..11..11..1y=1...1...1...1...1...1...1...1...y=.111.111.111.111.111.111.111.111y=..1...1...1...1...1...1...1...1.y=11.111.111.111.111.111.111.111.1y=.1...1...1...1...1...1...1...1..y=1.111.111.111.111.111.111.111.11y=111.111.111.111.111.111.111.111.y=...1...1...1...1...1...1...1...1y=1111....1111....1111....1111....y=....1111....1111....1111....1111y=.1.11.1..1.11.1..1.11.1..1.11.1.y=1.1..1.11.1..1.11.1..1.11.1..1.1y=..1111....1111....1111....1111..y=11....1111....1111....1111....11y=1..1.11.1..1.11.1..1.11.1..1.11.y=.11.1..1.11.1..1.11.1..1.11.1..1y=.1111....1111....1111....1111...y=1....1111....1111....1111....111y=11.1..1.11.1..1.11.1..1.11.1..1.y=..1.11.1..1.11.1..1.11.1..1.11.1y=1.11.1..1.11.1..1.11.1..1.11.1..y=.1..1.11.1..1.11.1..1.11.1..1.11y=...1111....1111....1111....1111.y=111....1111....1111....1111....103216161616161682482482424816161616161616161616161616161616The parity of B(a) is equal to the lowest bit of a.
Up to the a = 47 the bit-count varies by ±1 betweensuccessive values of B(a), the transition B(47) → B(48) changes the bit-count by 3. The sequence of theindices a where the bit-count changes by more than one is 47, 51, 59, 67, 75, 79, 175, 179, 187, 195, 203,207, 291, 299, 339, 347, 419, 427, 467, 475, 531, 539, . . .The sequence of fixed points is 0, 1, 6, 7, 18, 19, 20, 21, 106, 107, 108, 109, 120, 121, 126, 127, 258, 259,260, 261, 272, 273, 278, 279, 360, 361, 366, 367, 378, 379, 380, 381, 1546, 1547, 1548, 1549, 1560, 1561,1566, 1567, 1632, 1633, . .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
