NCPC 2013: Dance Reconstruction

題目大意

    對一個初始矩陣進行置換操做，已知經K次置換後獲得的矩陣爲，求一組可能的。

 

樣例解釋

    這裏只選取第二組樣例進行解釋。

4 2 3 4 1 2

2 3 4 1

    初始矩陣爲，根據Sample Output得知置換操做爲，第一次置換後獲得矩陣，第二次置換後獲得矩陣（和給出的一致），所以是一組可能的。

 

解題思路

    考察一個特定的置換操做（咱們暫且稱這種矩陣爲「置換矩陣」），咱們將其寫成另一種形式，這樣寫的好處就是咱們能方便的看出有2個循環節，並且咱們會發現無論置換多少次，都是循環節內的數相互動了動位置而已。這啓發咱們，也許能夠從置換K次的結果出發，找到其中的循環節，而後反推出置換1次的結果（也就是）。

    並且咱們發現，其實K次置換也對應的一個置換矩陣。好比，若是K=2就至關於另外一個置換矩陣，簡寫爲。不難發現，隨着置換次數的增長，置換矩陣的循環節的個數有可能會發生變化，那麼有什麼變化的規律麼？研究一下就會發現，若是置換1次對應的置換矩陣中有一個長度爲m的循環節，當置換次數K知足gcd(m,K) = 1時，這個循環節在置換K次對應的置換矩陣中還是一個長度爲m的循環節，不然這個循環節將會變成gcd(m,K)個長度爲m/gcd(m,K)的循環節。反過來說，假設置換K次以後有一個循環節的長度爲n，若是gcd(n,K) > 1，那麼這個必定是某個長度爲m的較長的循環節「分裂」出來，而且必須知足m/gcd(m,K) = n，同時包含他在內至少應有gcd(m,K)個長度爲n的循環節（這裏只是給出了一個對長度爲n的循環節的個數給出了一個較爲寬泛的約束，嚴格的約束會在後面提到）；若是gcd(n,K) = 1，那麼既有多是「分裂」出來的，也有可能不是「分裂」出來的，對於此題而言，咱們能夠認爲其不是「分裂」出來的，能夠簡化問題的求解。

    通過上述分析以後，造成了一個大體的解題思路——分狀況處理置換K次對應的置換矩陣中的各個循環節。設循環節的長度爲n，若是：

1. gcd(n,K) = 1。那麼直接反推出置換1次時對應的置換矩陣便可。
2. gcd(n,K) > 1。那麼這個必定是某個長度爲m的較大的循環節「分裂」出來，並且知足m/gcd(m,K) = n。這時，咱們須要把gcd(m, K)個長爲n小循環節合併成一個長爲m的大循環節並反推出置換1次時對應的置換矩陣。

    合併循環節的關鍵：假設長度爲n的循環節一共有cnt個，咱們要找到一個最小的gcd(m,K)使得m=gcd(m,K)*n成立。若是這個最小的gcd(m, K)是cnt的約數，那麼必定有解，將gcd(m, K)個長度爲n的小循環節合成一個長度m的大循環便可；不然無解。

    如何找到這個最小的gcd(m,K)？不能單純地認爲gcd(m,K) = gcd(n,K)，好比n = 2, K = 4，那麼gcd(m,K) = gcd(n,K) = 2, m = gcd(m,K)*n = 4，這樣就又能推得gcd(m,K) = 4，矛盾。上面那個例子之因此出現矛盾就是由於n承擔了一部分公因子，所以咱們要讓n不承擔任何公因子：對於任意一個n裏面的素因子p，若是K裏面也有p，那麼p在K中是多少次方就應當在gcd(m,K)中是多少次方。不過實際寫代碼的時候不必按素因子逐個去檢查，只要不斷地求K和n的公因數d，並將K除以d，直到d = 1爲止，將前面全部的d累乘就是gcd(m,K)。

    至於怎麼根據置換K次時對應的置換矩陣反推出置換1次時對應的置換矩陣，在這裏就不細說了，這個不難。

    P.S. 實際寫代碼的時候就會發現其實並不用分狀況討論gcd(n,K) = 1仍是gcd(n,K) > 1，由於當gcd(n,K) = 1時求得的最小的gcd(m,K)就是1，並且1是任何數的約數，天然就會認爲這個循環節是由1個循環節「分裂」獲得的（也就是沒「分裂」）。

#include<cstdio>
#include<cstring>
#include<vector>
#include<algorithm>
#define MAXN 10010
typedef std::vector<int> VI;
int N, K, f[MAXN];
std::vector<VI> c; // cycles
bool cmp(VI v1, VI v2) {
    return v1.size() < v2.size();
}
int gcd(int x, int y) {
    return y == 0 ? x : gcd(y, x % y);
}
bool vis[MAXN];
// find the cycles
void divideCycle() {
    c.clear();
    memset(vis, 0, sizeof(vis[0]) * (N + 1));
    for(int i = 1; i <= N; i ++) {
        if(!vis[i]) {
            VI v;
            for(int j = i; !vis[j]; j = f[j]) {
                vis[j] = true;
                v.push_back(j);
            }
            c.push_back(v);
        }
    }
    std::sort(c.begin(), c.end(), cmp);
}
int ans[MAXN], a[MAXN];
void process() {
    int cn = c.size();
    for(int i = 0, j = 0; i < cn; i ++) {
        if(i == cn - 1 || c[i].size() != c[i + 1].size()) {
            int n = c[i].size(); // n是小循環節的長度
            int t = K, d, g = 1; // g是分裂出來的循環節個數
            while((d = gcd(n, t)) != 1) {
                g *= d, t /= d;
            }
            int m = g * n; // m是大循環節的長度

            if((i - j + 1) % g == 0) {
                for(; j <= i; j += g) { // j~j+g-1 合併成一個大循環節
                    for(int s = 0; s < g; s ++) {
                        for(int id = 0, loc = s; id < n; id ++, loc = (loc + K) % m) {
                            a[loc] = c[j + s][id];
                        }
                    }
                    a[m] = a[0];
                    for(int id = 0; id < m; id ++) {
                        ans[a[id]] = a[id + 1];
                    }
                }
            } else {
                printf("Impossible\n");
                return ;
            }
        }
    }

    printf("%d", ans[1]);
    for(int i = 2; i <= N; i ++) {
        printf(" %d", ans[i]);
    }
    printf("\n");
}
int main() {
    while(scanf("%d%d", &N, &K) == 2) {
        for(int i = 1; i <= N; i ++) {
            scanf("%d", &f[i]);
        }
        divideCycle();
        process();
    }
    return 0;
}