poj1286 polya计数法

#include <iostream>
#include <algorithm>
#include <vector>

using namespace std;
typedef long long LL;
int n;

LL pow(LL A, LL B)
{
    LL res = 1;
    while(B > 0)
    {
        if(B%2==1)
            res = res*A;
        A = A*A;
        B = B/2;
    }
    return res;
}

vector<int> divisor(int n)   //获取n的约数因子
{
    vector<int> res;
    for(int i=1; i*i<=n; i++)
    {
        if(n%i == 0)
        {
            res.push_back(i);
            if(i!=n/i)
                res.push_back(n/i);
        }
    }
    return res;
}

int pi[100], npi;
void prime_factor(int n)
{
    npi = 0;
    for(int i=2; i*i<=n; i++)
    {
        if(n%i==0)
        {
            pi[npi++] = i;
            while(n%i==0) n/=i;
        }
    }
    if(n!=1) pi[npi++] = n;
}

LL solve(int n)
{
    vector<int> div;
    div = divisor(n);
    prime_factor(n);
    LL res = 0;
    for(int i=0; i<div.size(); i++)
    {
        LL euler = div[i];
        for(int j=0; j<npi; j++)
            if(div[i]%pi[j]==0)
            euler = euler/pi[j]*(pi[j]-1);
        res += pow(3, n/div[i])*euler;
    }
    if(n&1) res += n*pow(3, n/2+1);
    else res += (pow(3, n/2)+pow(3, n/2+1))*(n/2);
    return res/(2*n);
}

int main()
{
    while(cin>>n)
    {
        if(n==-1) break;
        if(n==0)
        {
            cout<<0<<endl;
            continue;
        }
        LL res = solve(n);
        cout<<res<<endl;
    }
    return 0;
}