1 #if !defined HAVE_PERMAPPLY_H__
2 #define HAVE_PERMAPPLY_H__
3 // This file is part of the FXT library.
4 // Copyright (C) 2010, 2012, 2018 Joerg Arndt
5 // License: GNU General Public License version 3 or later,
6 // see the file COPYING.txt in the main directory.
7
8
9 #include "fxttypes.h"
10 #include "ds/bitarray.h"
11 #include "restrict.h"
12
13
14
15 template <typename Type>
apply_permutation(const ulong * x,const Type * f,Type * restrict g,ulong n)16 void apply_permutation(const ulong *x, const Type *f, Type * restrict g, ulong n)
17 // Apply the permutation x[] to the array f[],
18 // i.e. set g[x[k]] <-- f[k] for all k
19 {
20 for (ulong k=0; k<n; ++k) g[x[k]] = f[k];
21 }
22 // -------------------------
23
24
25 template <typename Type>
26 void apply_permutation(const ulong *x, Type * restrict f, ulong n, bitarray *bp=nullptr)
27 // Apply the permutation x[] to the array f[]
28 // i.e. set f[x[k]] <-- f[k] for all k
29 // In-place version.
30 {
31 bitarray *tp = bp;
32 if ( nullptr==bp ) tp = new bitarray(n); // tags
33 tp->clear_all();
34
35 for (ulong k=0; k<n; ++k)
36 {
37 if ( tp->test_clear(k) ) continue; // already processed
38 tp->set(k);
39
40 // --- do cycle: ---
41 ulong i = k; // start of cycle
42 Type t = f[i];
43 ulong g = x[i];
44 while ( 0==(tp->test_set(g)) ) // cf. gray_permute()
45 {
46 Type tt = f[g];
47 f[g] = t;
48 t = tt;
49 g = x[g];
50 }
51 f[g] = t;
52 // --- end (do cycle) ---
53 }
54
55 if ( nullptr==bp ) delete tp;
56 }
57 // -------------------------
58
59
60 template <typename Type>
apply_inverse_permutation(const ulong * x,const Type * f,Type * restrict g,ulong n)61 void apply_inverse_permutation(const ulong *x, const Type *f, Type * restrict g, ulong n)
62 // Apply the inverse permutation of x[] to the array f[]
63 // i.e. set g[k] <-- f[x[k]] for all k
64 //
65 // E.g. after
66 // idx_quick_sort(f, n, x); apply(x, f, g, n);
67 // g[] == sorted( f[] )
68 //.
69 // If f[] is a permutation, then after sort_idx()
70 // x[] is its inverse and
71 // idx_quick_sort(f, n, x); apply_inverse_permutation(f, x, g, n);
72 // (note that x, f are swapped in apply_inverse_permutation())
73 // will make x[] == sorted( f[] )
74 // ... i.e. x * f == f * x
75 //
76 {
77 for (ulong k=0; k<n; ++k) g[k] = f[x[k]];
78 }
79 // -------------------------
80
81
82 template <typename Type>
83 void apply_inverse_permutation(const ulong *x, Type * restrict f, ulong n, bitarray *bp= nullptr)
84 // Apply the inverse permutation of x[] to the array f[]
85 // i.e. set f[k] <-- f[x[k]] for all k
86 // In-place version.
87 {
88 bitarray *tp = bp;
89 if ( nullptr==bp ) tp = new bitarray(n); // tags
90 tp->clear_all();
91
92 for (ulong k=0; k<n; ++k)
93 {
94 if ( tp->test_clear(k) ) continue; // already processed
95 tp->set(k);
96
97 // --- do cycle: ---
98 ulong i = k; // start of cycle
99 Type t = f[i];
100 ulong g = x[i];
101 while ( 0==(tp->test_set(g)) ) // cf. inverse_gray_permute()
102 {
103 f[i] = f[g];
104 i = g;
105 g = x[i];
106 }
107 f[i] = t;
108 // --- end (do cycle) ---
109 }
110
111 if ( nullptr==bp ) delete tp;
112 }
113 // -------------------------
114
115
116
117
118 #endif // !defined HAVE_PERMAPPLY_H__
119