Arndt - Algorithms for Programmers (523138), страница 24
Текст из файла (страница 24)
ThenEE−1Ē−1Ē−1======ĒĒRGRR G−1 RErr E −1R Ḡ RR Ḡ−1 R(7.20)(7.21)(7.22)(7.23)(7.24)(7.25)The relations between E and Ē areĒ −1 E=E −1 Ē = r = Xn−1(7.26)−1=Ē E −1 = X1(7.27)E Ēwhere Xx is the XOR permutation (with third argument equal to the subscript x, see section 7.7,page 132).Let Z̄ be the zip_rev-permutation described in section 7.6, page 131. ThenZ̄Z̄−1==G−1 EE −1 G(7.28)(7.29)CHAPTER 7. PERMUTATIONS140With an implementation of E and its inverse that is as fast as the Gray code permutation (TBD) thegiven relations constitute an algorithm that should be significantly faster than the triple revbin permutecode given on page 131.7.12Factorizing permutationsIn this section we will see some algorithms that use a certain type of decomposition (factorization interm of matrices) of some permutations we have studied so far.
The resulting algorithms involve proportional n · log(n) computations for length-n arrays. This might seem to render the schemes worthless(we generally can obtain a permutation with work proportional to n). There are, however, situationswhere on can use the algorithms advantageously. Firstly, with bit manipulations, where the whole binarywords are modified in one (or a few) statements. The corresponding algorithms are therefore only proportional log(n).
Secondly, when the algorithm can be used implicitly in order to integrate the permutationin a fast transform. The work can sometimes be reduced to zero in that case.The most simple example might the reversion viatemplate <typename Type>void perm1(Type *f, ulong ldn){ulong n = 1UL<<ldn;for (ulong ldk=1; ldk<=ldn; ++ldk) // starting with length 2{ulong k = 1UL<<ldk;for (ulong j=0; j<n; j+=k) func(f+j, ldk);}}where func swaps the upper and lower half of an array:template <typename Type>void func(Type *f, ulong ldn){ swap(f, f+n/2, n/2); }This idea has been exploited in section 8.7 in order to obtain a bit-reversal routine. The bit-reversalis self-inverse, therefore one can alternatively execute the steps in reverse order and still get the samepermutation:template <typename Type>void perm2(Type *f, ulong ldn){ulong n = 1UL<<ldn;for (ulong ldk=ldn; ldk>0; --ldk) // starting with full length{ulong k = 1UL<<ldk;for (ulong j=0; j<n; j+=k) func(f+j, ldk);}}Note that func in turn can be obtained viareverse(f, n);reverse(f, n/2);reverse(f+n/2, n/2);or the same statements in reversed order.Let us try a not so obvious example, use perm1 with func defined as:swap(f+n/2, f+n/2+n/4, n/4);The resulting permutation is the Gray-permutation.
The other way round (using perm2) one gets theinverse Gray-permutation. Using func defined asreverse(f+n/2, n/2);CHAPTER 7. PERMUTATIONS141one gets the Gray-permutation through perm2, its inverse through perm1. This idea has been used inthe core-routine for the sequency-ordered Walsh transform described in section 5.6. The work for theGray-permutation has been completely vaporized there. Note that the routine that swaps the halves ofthe upper half array could be obtained as either of//////inverse_gray_permute(f, n);gray_permute(f, n/2);gray_permute(f+n/2, n/2);=^=inverse_gray_permute(f, n/2); inverse_gray_permute(f+n/2, n/2);gray_permute(f, n);Similarly, the routine that reverses the upper half can be obtained as//////gray_permute(f, n/2);gray_permute(f+n/2, n/2);inverse_gray_permute(f, n);=^=gray_permute(f, n);inverse_gray_permute(f, n/2); inverse_gray_permute(f+n/2, n/2);Using func defined asswap(f+n/4, f+n/2, n/4);and perm2 one gets the zip-permutation, perm1 gives the inverse.Using func (lazily implementated as)Type g[n];copy(f, g, n/4);copy(f+n/4, g+n/2+n/4, n/4);copy(f+n/2, g+n/4, n/4);copy(f+n/2+n/4, g+n/2, n/4);copy(g, f, n);////////quarter: 0 -->1 --> 32 --> 13 --> 2which was obtained usingzip_rev(f, n);unzip_rev(f, n/2);unzip_rev(f+n/2, n/2);both the reversed zip-permutation and its inverse can be computed in the now hopefully obvious way.The revbin-permutation can be generated through the zip-permutation or its inverse.
However, the zippermutation is the more complicate one, so absorbing the revbin-permutation into fast transforms doesnot seem to be easy. The other way round it makes more sense:revbin_permute(f, n/2);revbin_permute(f, n);revbin_permute(f+n/2, n/2);Is a convenient (though not the most effective) way to compute the zip-permutation.Clearly, the idea presented here is in analogy with the decomposition of linear transforms. Finally thepermutations are (very simple forms of) linear transforms.7.13General permutationsSo far we treated special permutations that occurred as part of other algorithms. It is instructive tostudy permutations in general with the operations (as composition and inverse) on them.7.13.1Basic definitionsA straight forward way to describe a permutation is to consider the array of indices that for the original(un-permuted) data would be the length-n canonical sequence 0, 1, 2, . .
. , n − 1. The mentioned trivialCHAPTER 7. PERMUTATIONS142sequence describes the ‘do-nothing’ permutation or identity (wrt. composition of permutations). Theconcept is best described by the routine that applies a given permutation x on an array of data f : afterthe routine has finished the array g will contain the elements of f reordered according to xtemplate <typename Type>void apply_permutation(const ulong *x, const Type *f, Type * restrict g, ulong n)// apply x[] on f[]// i.e.
g[k] <-- f[x[k]] \forall k{for (ulong k=0; k<n; ++k) g[k] = f[x[k]];}[FXT: apply permutation in perm/permapply.h] An example using strings (arrays of characters): Thepermutation described by x = {7, 6, 3, 2, 5, 1, 0, 4} and the input dataf ="ABadCafe" would produceg ="efdaaBAC"All routines in this and the following section are declared in [FXT: file perm/permutation.h]Triviallybool is_identity(const ulong *f, ulong n)// check whether f[] is the identical permutation,// i.e.
whether f[k]==k for all k= 0...n-1{for (ulong k=0; k<n; ++k) if ( f[k] != k ) return false;return true;}A fixed point of a permutation is an index where the element is not moved:ulong count_fixed_points(const ulong *f, ulong n)// return number of fixed points in f[]{ulong ct = 0;for (ulong k=0; k<n; ++k) if ( f[k] == k ) ++ct;return ct;}A derangement is a permutation that has no fixed points (i.e. that moved every element to anotherposition so count_fixed_points() returns zero). To check whether a permutation is a derangement ofidentity use:bool is_derangement(const ulong *f, ulong n)//// check whether f[] is a derangement of identity// i.e.
whether f[k]!=k for all k//{for (ulong k=0; k<n; ++k) if ( f[k] == k ) return false;return true;}Whether two permutations are derangements of each other is obviously found by:bool is_derangement(const ulong *f, const ulong *g, ulong n)// check whether f[] is a derangement of g[],// i.e. whether f[k]!=g[k] for all k{for (ulong k=0; k<n; ++k) if ( f[k] == g[k] ) return false;return true;}To check whether a given array really describes a valid permutation one has to verify that each indexappears exactly once. The bitarray class described in 10.6 allows us to do the job without modificationof the input (like e.g.
sorting):CHAPTER 7. PERMUTATIONS143bool is_valid_permutation(const ulong *f, ulong n, bitarray *bp/*=0*/)// check whether all values 0...n-1 appear exactly once{// check whether any element is out of range:for (ulong k=0; k<n; ++k) if ( f[k]>=n ) return false;// check whether values are unique:bitarray *tp = bp;if ( 0==bp ) tp = new bitarray(n); // tagstp->clear_all();ulong k;for (k=0; k<n; ++k){if ( tp->test_set(f[k]) ) break;}if ( 0==bp ) delete tp;return (k==n);}7.13.2Compositions of permutationsOne can apply arbitrary many permutations to an array, one by one.
The resulting permutation is calledthe composition of the applied permutations. As an example, the check whether some permutation g isequal to f applied twice, or f · f , or f squared use:bool is_square(const ulong *f, const ulong *g, ulong n)// whether f * f == g as a permutation{for (ulong k=0; k<n; ++k) if ( g[k] != f[f[k]] ) return false;return true;}A permutation f is said to be the inverse of another permutation g if it undoes its effect, that is f · g = id(likewise g · f = id):bool is_inverse(const ulong *f, const ulong *g, ulong n)// check whether f[] is inverse of g[]{for (ulong k=0; k<n; ++k) if ( f[g[k]] != k ) return false;return true;}A permutation that is its own inverse (like the revbin-permutation) is called an involution. Checkingthat is easy:bool is_involution(const ulong *f, ulong n)// check whether max cycle length is <= 2{for (ulong k=0; k<n; ++k) if ( f[f[k]] != k )return true;}return false;Finding the inverse of a given permutation is trivial:void make_inverse(const ulong *f, ulong * restrict g, ulong n)// set g[] to the inverse of f[]{for (ulong k=0; k<n; ++k) g[f[k]] = k;}However, if one wants to do the operation in-place a little bit of thought is required.
The idea underlyingall subsequent routines working in-place is that every permutation entirely consists of disjoint cycles. Acycle (of a permutation) is a subset of the indices that is rotated (by one) by the permutation. The termdisjoint means that the cycles do not ‘cross’ each other. While this observation is pretty trivial it allowsus to do many operations by following the cycles of the permutation, one by one, and doing the necessaryoperation on each of them.