[POJ 2888] | Eurce's Blog

数论。

本来这题是在分类表里靠后的位置的，但由于最近在啃具体数学，第四章讲了“莫比乌斯反演”，就顺势把题做一下。

题意是这样的，有一个长度为n的珠串，每个珠子的颜色有m种，但其中有k对颜色是不能相邻的，问在剔除因旋转导致的重复后总共有多少种可能的珠串。

先讲这题的简化版，是POJ 2154，区别在于没有珠子间不能连接的限制。

其基本的思想仍然是Polya定理。但是这里由于n值太大，无法直接枚举出所有的置换群。于是需要考虑将具有相同数量循环节的置换群合并处理。具体的合并方式是要利用到欧拉函数phi()的性质的。总共的置换群有n个，对n的每个因子d，有n/d个循环节的置换群其每个循环节长度为d，而这样的循环群总数是phi(d)个。又因为欧拉函数满足sum_{d\n} phi(d) == n，因此计数式子就可以转化为 (sum_{d\n} phi(d) * n^(n/d)) / n。计算phi(d)的时间复杂度是O(lgd)，而n的因子总数是不超过1000的（最多9个不同因子），n的n/d次方是可以在O(lg(n/d))时间内计算得到的，由此问题便能在合理时间内求解了。

另外就是计算的时候，要注意利用模运算的性质，使用int来处理（比long long快很多），同时尽量减少进行模运算的次数，否则还是会TLE的。

下面是2154的代码：

#include<cstdio>
#include<cstring>
#include<cmath>

#define pl 4000

int n,m,o;
int p[pl]={2,3},pn[pl][2],pnl;

void initp()
{
    int i,j,k,l;
    for(i=2;i<pl;i++)
    {
        for(j=p[i-1]+2;1;j+=2)
        {
            for(k=0,l=sqrt(j);p[k]<=l&&j%p[k];k++);
            if(p[k]>l)
            {
                p[i]=j;
                break;
            }
        }
    }
}

inline int phi(int d)
{
    int i,j,k;
    for(i=0,j=d,k=sqrt(d);p[i]<=k;i++)
        if(d%p[i]==0)
        {
            j-=j/p[i];
            while(d%p[i]==0)
                d/=p[i];
        }
    if(d>1)
        j-=j/d;
    return j%m;
}

inline int pm(int d)
{
    int i,j,k;
    k=n%m;
    for(i=1;i<=d;i<<=1);
    i>>=1;
    for(j=1;i;i>>=1)
    {
        j=j*j%m;
        if(i&d)
            j=j*k%m;
    }
    return j;
}

void init()
{
    int i,j,k,l;
    memset(pn,0,sizeof pn);
    for(i=pnl=0,j=n,k=sqrt(n);p[i]<=k&&j>1;i++)
    {
        if(j%p[i]==0)
        {
            pn[pnl][0]=p[i];
            while(j%p[i]==0)
            {
                j/=p[i];
                pn[pnl][1]++;
            }
            pnl++;
        }
    }
    if(j>1)
    {
        pn[pnl][0]=j;
        pn[pnl][1]=1;
        pnl++;
    }
}

void calc(int i,int d)
{
    if(i==pnl)
    {
        o=(o+phi(d)*pm(n/d-1))%m;
        return;
    }
    int j,k,l;
    for(j=0,k=1;j<=pn[i][1];j++,k*=pn[i][0])
        calc(i+1,d*k);
}

void solve()
{
    o=0;
    calc(0,1);
    printf("%d\n",o);
}

int main()
{
    initp();
    scanf("%*d");
    while(~scanf("%d%d",&n,&m))
    {
        init();
        solve();
    }
    return 0;
}

100

101

102

103

104

105

106

107

108

109

110

111

#include<cstdio>

#include<cstring>

#include<cmath>

#define pl 4000

int n,m,o;

int p[pl]={2,3},pn[pl][2],pnl;

void initp()

{

int i,j,k,l;

for(i=2;i<pl;i++)

{

for(j=p[i-1]+2;1;j+=2)

{

for(k=0,l=sqrt(j);p[k]<=l&&j%p[k];k++);

if(p[k]>l)

{

p[i]=j;

break;

}

inline int phi(int d)

{

int i,j,k;

for(i=0,j=d,k=sqrt(d);p[i]<=k;i++)

if(d%p[i]==0)

{

j-=j/p[i];

while(d%p[i]==0)

d/=p[i];

}

if(d>1)

j-=j/d;

return j%m;

}

inline int pm(int d)

{

int i,j,k;

k=n%m;

for(i=1;i<=d;i<<=1);

i>>=1;

for(j=1;i;i>>=1)

{

j=j*j%m;

if(i&d)

j=j*k%m;

}

return j;

}

void init()

{

int i,j,k,l;

memset(pn,0,sizeof pn);

for(i=pnl=0,j=n,k=sqrt(n);p[i]<=k&&j>1;i++)

{

if(j%p[i]==0)

{

pn[pnl][0]=p[i];

while(j%p[i]==0)

{

j/=p[i];

pn[pnl][1]++;

}

pnl++;

}

if(j>1)

{

pn[pnl][0]=j;

pn[pnl][1]=1;

pnl++;

}

void calc(int i,int d)

{

if(i==pnl)

{

o=(o+phi(d)*pm(n/d-1))%m;

return;

}

int j,k,l;

for(j=0,k=1;j<=pn[i][1];j++,k*=pn[i][0])

calc(i+1,d*k);

}

void solve()

{

o=0;

calc(0,1);

printf("%d\n",o);

}

int main()

{

initp();

scanf("%*d");

while(~scanf("%d%d",&n,&m))

{

init();

solve();

}

return 0;

}

回到2888。这题额外增加了珠子间的连接限制，对应的处理方式是这样的：初始化矩阵a[m][m]元素全为1，若i和j不能相邻，则a[i][j]=a[j][i]=0，然后计算b = a^d，则长度为d的珠串其总可能数是sum{b[i][i]}。这里当i<=j时，(a^d)[i][j]表示的是，以种类i开始长度为d，结尾能够与种类j的珠子相连的珠串总数。至于为什么是这样，推一下就明白了。于是这里计数式子就变成了(sum_{d\n} phi(d) * sum{(a^(n/d))[i][i]}) / n。落实到具体实现时，在对n做除法时，要转化成乘以模运算下n的乘法逆元。前面的每一步处理也要控制好在int的范围内以及限制模运算的使用次数。

#include<cstdio>
#include<cstring>
#include<cmath>

#define pl 4000
#define mod 9973

int n,m,o;
int p[pl]={2,3},pn[pl][2],pnl;
int a[10][10],b[10][10],t[10][10];

void initp()
{
    int i,j,k,l;
    for(i=2;i<pl;i++)
    {
        for(j=p[i-1]+2;1;j+=2)
        {
            for(k=0,l=sqrt(j);p[k]<=l&&j%p[k];k++);
            if(p[k]>l)
            {
                p[i]=j;
                break;
            }
        }
    }
}

int phi(int d)
{
    int i,j,k;
    for(i=0,j=d,k=sqrt(d);p[i]<=k;i++)
        if(d%p[i]==0)
        {
            j-=j/p[i];
            while(d%p[i]==0)
                d/=p[i];
        }
    if(d>1)
        j-=j/d;
    return j%mod;
}

void mul(int a[10][10],int b[10][10])
{
    int i,j,k,l;
    memset(t,0,sizeof t);
    for(i=0;i<m;i++)
        for(j=0;j<m;j++)
        {
            for(k=0;k<m;k++)
                t[i][j]+=a[i][k]*b[k][j];
            t[i][j]%=mod;
        }
    memcpy(a,t,sizeof t);
}

int pm(int d)
{
    int i,j,k,l;
    memset(b,0,sizeof b);
    for(i=0;i<m;i++)
        b[i][i]=1;
    for(i=1;i<=d;i<<=1);
    for(i>>=1;i;i>>=1)
    {
        mul(b,b);
        if(i&d)
            mul(b,a);
    }
    for(i=j=0;i<m;i++)
        j+=b[i][i];
    return j%mod;
}

void init()
{
    int i,j,k,l;
    for(i=0;i<m;i++)
        for(j=0;j<m;j++)
            a[i][j]=1;
    while(o--)
    {
        scanf("%d%d",&i,&j);
        a[i-1][j-1]=a[j-1][i-1]=0;
    }
    memset(pn,0,sizeof pn);
    for(i=pnl=0,j=n,k=sqrt(n);p[i]<=k&&j>1;i++)
    {
        if(j%p[i]==0)
        {
            pn[pnl][0]=p[i];
            while(j%p[i]==0)
            {
                j/=p[i];
                pn[pnl][1]++;
            }
            pnl++;
        }
    }
    if(j>1)
    {
        pn[pnl][0]=j;
        pn[pnl][1]=1;
        pnl++;
    }
}

void calc(int i,int d)
{
    if(i==pnl)
    {
        o=(o+phi(d)*pm(n/d))%mod;
        return;
    }
    int j,k,l;
    for(j=0,k=1;j<=pn[i][1];j++,k*=pn[i][0])
        calc(i+1,d*k);
}

void solve()
{
    o=0;
    calc(0,1);
    int i,j;
    for(i=1,j=n%mod;(i*j)%mod!=1;i++);
    printf("%d\n",o*i%mod);
}

int main()
{
    initp();
    scanf("%*d");
    while(~scanf("%d%d%d",&n,&m,&o))
    {
        init();
        solve();
    }
    return 0;
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

#include<cstdio>

#include<cstring>

#include<cmath>

#define pl 4000

#define mod 9973

int n,m,o;

int p[pl]={2,3},pn[pl][2],pnl;

int a[10][10],b[10][10],t[10][10];

void initp()

{

int i,j,k,l;

for(i=2;i<pl;i++)

{

for(j=p[i-1]+2;1;j+=2)

{

for(k=0,l=sqrt(j);p[k]<=l&&j%p[k];k++);

if(p[k]>l)

{

p[i]=j;

break;

}

int phi(int d)

{

int i,j,k;

for(i=0,j=d,k=sqrt(d);p[i]<=k;i++)

if(d%p[i]==0)

{

j-=j/p[i];

while(d%p[i]==0)

d/=p[i];

}

if(d>1)

j-=j/d;

return j%mod;

}

void mul(int a[10][10],int b[10][10])

{

int i,j,k,l;

memset(t,0,sizeof t);

for(i=0;i<m;i++)

for(j=0;j<m;j++)

{

for(k=0;k<m;k++)

t[i][j]+=a[i][k]*b[k][j];

t[i][j]%=mod;

}

memcpy(a,t,sizeof t);

}

int pm(int d)

{

int i,j,k,l;

memset(b,0,sizeof b);

for(i=0;i<m;i++)

b[i][i]=1;

for(i=1;i<=d;i<<=1);

for(i>>=1;i;i>>=1)

{

mul(b,b);

if(i&d)

mul(b,a);

}

for(i=j=0;i<m;i++)

j+=b[i][i];

return j%mod;

}

void init()

{

int i,j,k,l;

for(i=0;i<m;i++)

for(j=0;j<m;j++)

a[i][j]=1;

while(o--)

{

scanf("%d%d",&i,&j);

a[i-1][j-1]=a[j-1][i-1]=0;

}

memset(pn,0,sizeof pn);

for(i=pnl=0,j=n,k=sqrt(n);p[i]<=k&&j>1;i++)

{

if(j%p[i]==0)

{

pn[pnl][0]=p[i];

while(j%p[i]==0)

{

j/=p[i];

pn[pnl][1]++;

}

pnl++;

}

if(j>1)

{

pn[pnl][0]=j;

pn[pnl][1]=1;

pnl++;

}

void calc(int i,int d)

{

if(i==pnl)

{

o=(o+phi(d)*pm(n/d))%mod;

return;

}

int j,k,l;

for(j=0,k=1;j<=pn[i][1];j++,k*=pn[i][0])

calc(i+1,d*k);

}

void solve()

{

o=0;

calc(0,1);

int i,j;

for(i=1,j=n%mod;(i*j)%mod!=1;i++);

printf("%d\n",o*i%mod);

}

int main()

{

initp();

scanf("%*d");

while(~scanf("%d%d%d",&n,&m,&o))

{

init();

solve();

}

return 0;

}

虽然上面写了这么多但很多都是现学现卖的OTL。不过数论还真是很有意思的东西……