正六面体をn色で塗るパターン数

ポリアの数え上げを使って、正六面体をn色で塗るパターン（回転して同じになるものは同じとみなす）を数えてみました。

回転群の要素を列挙するクールな方法が思いつかなかったのですが、面の数が６個なので以下のようにしました。

(0,1,2,3,4,5)の順列をすべて挙げ、それが正六面体の回転群の元になっているかどうか判定する。
判定条件は以下のとおり。
『すべての面fについて、fを手前に向けてfに接する面を時計回りに挙げていったときに得られる列が写像前と写像後で（シフトすれば）一致していれば、回転群の要素である。』

上の判定で正しく回転群を抽出できていることは以下のように確認しました。

回転群の集合をA, 上記の判定条件を満たす集合をBとする。
要素xがAに含まれるなら、xはBに含まれる。が言えるため、AはBの部分集合になる。
Aの要素数は、上面にどれを持ってくるか * 側面の回転 = 6 * 4 = 24となる。
つまり、Bの要素数が24となれば、BとAは同じ集合であると言える。実際にBの要素数は24となったためA = B。

以上をふまえて書いたプログラムが以下です。
実は今回計算した数は多項式で計算できるようなので、それを使って検算しました。
 正しく計算できてました。

#include <algorithm>
#include <cassert>
#include <iostream>

using namespace std;

int adjacentFaces[][4] = {
    {1, 4, 3, 5}, 
    {2, 4, 0, 5}, 
    {3, 4, 1, 5}, 
    {0, 4, 2, 5}, 
    {1, 2, 3, 0}, 
    {3, 2, 1, 0}
};

bool isRotate(int a[]) {
    for (int i = 0; i < 6; i++) {
        int p = -1;
        for (int j = 0; j < 4; j++) {
            if (a[adjacentFaces[i][j]] == adjacentFaces[a[i]][0]) {
                p = j;
                break;
            }
        }
        if (p == -1)
            return false;

        for (int j = 0; j < 4; j++)
            if (adjacentFaces[a[i]][j] != a[adjacentFaces[i][(p+j)%4]])
                return false;
    }
    return true;
}

long long pow(long long x, long long p) {
    long long ret = 1;
    while (p > 0) {
        if (p & 1)
            ret *= x;
        x = x * x;
        p >>= 1;
    }
    return ret;
}

int cycleNum(int p[]) {
    bool vis[6] = {};
    
    int ret = 0;
    for (int i = 0; i < 6; i++) {
        if (vis[i])
            continue;
        ++ret;
        for (int j = i; !vis[j]; j = p[j])
            vis[j] = true;
    }

    return ret;
}

long long solve(int color) {
    int p[] = {0, 1, 2, 3, 4, 5};

    int cnt = 0;
    long long sum = 0;
    do {
        if (!isRotate(p))
            continue;
        ++cnt;
        sum += pow(color, cycleNum(p));

    } while (next_permutation(p, p+6));

    return sum / cnt;
}

int main(int argc, char **argv) {
    int c;
    cin >> c;
    
    long long ret = solve(c);
    assert(ret == (pow(c, 6) + 3*pow(c, 4) + 12*pow(c, 3) + 8*pow(c, 2)) / 24);
    cout << ret << endl;
    
    return 0;
}