Arndt - Algorithms for Programmers (523138), страница 21
Текст из файла (страница 21)
One can do better by solvinga slightly different problem.7.1.2A fast versionThe key idea is to update the value x̃ from the value x]− 1. As x is one added to x − 1, x̃ is one ‘reversed’added to x]− 1. If one finds a routine for that ‘reversed add’ update much of the computation can besaved.A routine to update r, that must be the same as the the result of revbin(x-1, n) to what would be theresult of revbin(x, n)CHAPTER 7.
PERMUTATIONS122function revbin_update(r, n){do{n := n >> 1r := r^n // bitwise XOR} while ((r&n) == 0)return r}In C this can be cryptified to an efficient piece of code:inline unsigned revbin_update(unsigned r, unsigned n){for (unsigned m=n>>1; (!((r^=m)&m)); m>>=1);return r;}[FXT: revbin update in auxbit/revbin.h]Now we are ready for a fast revbin-permute routine:procedure revbin_permute(a[], n)// a[0..n-1] input,result{if n<=2 returnr := 0 // the reversed 0for x:=1 to n-1{r := revbin_update(r, n) // inline meif r>x then swap(a[x],a[r])}}This routine is several times faster than the naive version.
revbin_update() needs for half of the callsjust one iteration because in half of the updates just the leftmost bit changes1 , in half of the remainingupdates it needs two iterations, in half of the still remaining updates it needs three and so on. The total4number of operations done by revbin_update() is therefore proportional to n ( 12 + 42 + 38 + 16+· · ·+ logn2 (n) )Plog2 (n) j= n j=1 2j . For n large this sum is close to 2n.
Thereby the asymptotics of revbin_permute() isimproved from proportional n log(n) to proportional n.7.1.3How many swaps?√How many swap()-statements will be executed in total for different n? About n − n, as there are onlyfew numbers with symmetric bit patterns: for even√ log2 (n) =: 2 b the left half of the bit pattern must bebthe reversed of the right half. There are 2 = 22b such numbers. For odd log2 (n) =: 2 b + 1 there aretwice as much symmetric patterns: the bit in the middle does not matter and can be 0 or 1.n248163264210220∞2 # swaps0241224569920.999 · √220n− n# symm. pairs22448832210√nSummarizing: almost all ‘revbin-pairs’ will be swapped by revbin_permute().1 correspondingto the change in only the rightmost bit if one is added to an even numberCHAPTER 7. PERMUTATIONS7.1.4123A still faster versionThe following table lists indices versus their revbin-counterpart.
The subscript 2 indicates printing inbase 2, ∆ := xe − x]− 1 and an ‘y’ in the last column marks index pairs where revbin_permute() willswap elements.x012345678910111213141516171819202122232425262728293031Observation one: ∆ =n2x20000000001000100001100100001010011000111010000100101010010110110001101011100111110000100011001010011101001010110110101111100011001110101101111100111011111011111x̃20000010000010001100000100101000110011100000101001001010110100011010110011101111000001100010100111001001011010101101111010001110011010111101100111101110111111111x̃016824420122821810266221430117925521132931911277231531∆-3116-816-2016-816-2616-816-2016-816-2916-816-2016-816-2616-816-2016-816x̃ > x?yyyyyyyyyyyyfor all odd x.Observation two: if for even x < n2 there is a swap (for the pair x, x̃) then there is also a swap for thepair n − 1 − x, n − 1 − x̃.
As x < n2 and x̃ < n2 one has n − 1 − x > n2 and n − 1 − x̃ > n2 , i.e. the swapsare independent.There should be no difficulties to cast these observations into a routine to put data into revbin order:procedure revbin_permute(a[], n){if n<=2 returnnh := n/2r := 0 // the reversed 0x := 1while x<nh{// x odd:r := r + nhswap(a[x], a[r])x := x + 1// x even:r := revbin_update(r,n) // inline meif r>x thenCHAPTER 7.
PERMUTATIONS{}124swap(a[x], a[r])swap(a[n-1-x], a[n-1-r])}x := x + 1}The revbin_update() would be in C, inlined and the first stage of the loop extractedr^=nh;for (unsigned m=(nh>>1); !((r^=m)&m); m>>=1){}The code above is an ideal candidate to derive an optimized version for zero padded data:procedure revbin_permute0(a[], n){if n<=2 returnnh := n/2r := 0 // the reversed 0x := 1while x<nh{// x odd:r := r + nha[r] := a[x]a[x] := 0x := x + 1// x even:r := revbin_update(r, n) // inline meif r>x then swap(a[x], a[r])// both a[n-1-x] and a[n-1-r] are zerox := x + 1}}One could carry the scheme that lead to the ‘faster’ revbin permute procedures further, e.g.
using 3hard-coded constants ∆1 , ∆2 , ∆3 depending on whether x mod 4 = 1, 2, 3 only calling revbin_update()for x mod 4 = 0. However, the code quickly gets quite complicated and there seems to be no measurablegain in speed, even for very large sequences.If, for complex data, one works with separate arrays for real and imaginary part2 one might be tempted todo away with half of the bookkeeping as follows: write a special procedure revbin_permute(a[],b[],n)that shall replace the two successive calls revbin_permute(a[],n) and revbin_permute(b[],n) andafter each statement swap(a[x],a[r]) has inserted a swap(b[x],b[r]). If you do so, be prepared fordisaster! Very likely the real and imaginary element for the same index lie apart in memory by a powerof two, leading to one hundred percent cache miss for the typical computer.
Even in the most favorablecase the cache miss rate will be increased. Do expect to hardly ever win anything noticeable but in mostcases to lose big. Think about it, whisper “direct mapped cache” and forget it.7.1.5The real world versionFinally we remark that the revbin_update can be optimized by usage of a small (length BITS_PER_LONG)table containing the reflected bursts of ones that change on the lower end with incrementing. A routinethat utilizes this idea, optionally uses the CPU-bitscan instruction(cf.
section 8.2) and further allows toselect the amount of symmetry optimizations looks like#include "inline.h" // swap()#include "fxttypes.h"#include "bitsperlong.h" // BITS_PER_LONG#include "revbin.h" // revbin(), revbin_update()#include "bitasm.h"#if defined BITS_USE_ASM#include "bitlow.h" // lowest_bit_idx()#define RBP_USE_ASM // use bitscan if available, comment out to disable2 asopposed to: using a data type ‘complex’ with real and imaginary part of each number in consecutive placesCHAPTER 7. PERMUTATIONS#endif // defined BITS_USE_ASM#define RBP_SYMM 4 // 1, 2, 4 (default is 4)#define idx_swap(f, k, r) { ulong kx=(k), rx=(r); swap(f[kx], f[rx]); }template <typename Type>void revbin_permute(Type *f, ulong n){if ( n<=8 ){if ( n==8 ){swap(f[1], f[4]);swap(f[3], f[6]);}else if ( n==4 ) swap(f[1], f[2]);return;}const ulong nh = (n>>1);ulong x[BITS_PER_LONG];x[0] = nh;{ // initialize xor-table:ulong i, m = nh;for (i=1; m!=0; ++i){m >>= 1;x[i] = x[i-1] ^ m;}}#if ( RBP_SYMM >= 2 )const ulong n1 = n - 1;// = 11111111#if ( RBP_SYMM >= 4 )const ulong nx1 = nh - 2;// = 01111110const ulong nx2 = n1 - nx1; // = 10111101#endif // ( RBP_SYMM >= 4 )#endif // ( RBP_SYMM >= 2 )ulong k=0, r=0;while ( k<n/RBP_SYMM ) // n>=16, n/2>=8, n/4>=4{// ----- k%4 == 0:if ( r>k ){swap(f[k], f[r]); // <nh, <nh 11#if ( RBP_SYMM >= 2 )idx_swap(f, n1^k, n1^r); // >nh, >nh 00#if ( RBP_SYMM >= 4 )idx_swap(f, nx1^k, nx1^r); // <nh, <nh 11idx_swap(f, nx2^k, nx2^r); // >nh, >nh 00#endif // ( RBP_SYMM >= 4 )#endif // ( RBP_SYMM >= 2 )}r ^= nh;++k;// ----- k%4 == 1:if ( r>k ){swap(f[k], f[r]); // <nh, >nh 10#if ( RBP_SYMM >= 4 )idx_swap(f, n1^k, n1^r); // >nh, <nh 01#endif // ( RBP_SYMM >= 4 )}{ // scan for lowest unset bit of k:#ifdef RBP_USE_ASMulong i = lowest_bit_idx(~k);#elseulong m = 2, i = 1;while ( m & k ) { m <<= 1; ++i; }#endif // RBP_USE_ASMr ^= x[i];}++k;// ----- k%4 == 2:if ( r>k ){125CHAPTER 7.
PERMUTATIONS126swap(f[k], f[r]); // <nh, <nh 11( RBP_SYMM >= 2 )idx_swap(f, n1^k, n1^r); // >nh, >nh 00#endif // ( RBP_SYMM >= 2 )}r ^= nh;++k;// ----- k%4 == 3:if ( r>k ){swap(f[k], f[r]);// <nh, >nh 10#if ( RBP_SYMM >= 4 )idx_swap(f, nx1^k, nx1^r);// <nh, >nh 10#endif // ( RBP_SYMM >= 4 )}{ // scan for lowest unset bit of k:#ifdef RBP_USE_ASMulong i = lowest_bit_idx(~k);#elseulong m = 4, i = 2;while ( m & k ) { m <<= 1; ++i; }#endif // RBP_USE_ASMr ^= x[i];}++k;}}#if. . . not the most readable piece of code but a nice example for a real-world optimized routine.This is [FXT: revbin permute in perm/revbinpermute.h], see [FXT: revbin permute0 inperm/revbinpermute0.h] for the respective version for zero padded data.7.2The radix permutationThe radix-permutation is the generalization of the revbin-permutation (corresponding to radix 2) toarbitrary radices.C++ code for the radix-r permutation of the array f[]:extern ulong nt[]; // nt[] = 9, 90, 900 for r=10, x=3extern ulong kt[]; // kt[] = 1, 10, 100 for r=10, x=3template <typename Type>void radix_permute(Type *f, ulong n, ulong r)//// swap elements with index pairs i, j were the// radix-r representation of i and j are mutually// digit-reversed (e.g.
436 <--> 634)//// This is a radix-r generalization of revbin_permute()// revbin_permute(f, n) =^= radix_permute(f, n, 2)//// must have:// n == p**x for some x>=1// r >= 2//{ulong x = 0;nt[0] = r-1;kt[0] = 1;while ( 1 ){ulong z = kt[x] * r;if ( z>n ) break;++x;kt[x] = z;nt[x] = nt[x-1] * r;}// here: n == p**xCHAPTER 7. PERMUTATIONS}127for (ulong i=0, j=0; i < n-1; i++){if ( i<j ) swap(f[i], f[j]);ulong t = x - 1;ulong k = nt[t]; // =^= k = (r-1) * n / r;while ( k<=j ){j -= k;k = nt[--t]; // =^= k /= r;}j += kt[t]; // =^= j += (k/(r-1));}[FXT: radix permute in perm/radixpermute.h]TBD: mixed-radix permute7.3In-place matrix transpositionTo transpose a nr × nc - matrix first identify the position i of then entry in row r and column c:i =r · nc + c(7.1)After the transposition the element will be at position i0 in the transposed n0r × n0c - matrixi0r0 · n0c + c0(7.2)c · nr + r(7.3)= c · nr · nc + r · nc(7.4)=Obviously, r0 = c, c0 = r, n0r = nc and n0c = nr , so:i0=Multiply the last equation by nci0 · ncWith n := nr · nc and r · nc = i − c we geti0 · nc =i =c·n+i−ci0 · nc + c · (n − 1)(7.5)(7.6)Take the equation modulo n − 1 to get3i ≡i0 · ncmod (n − 1)(7.7)That is, the transposition moves the element i = i0 · nc to position i0 .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
