hdu 2865 Birthday Toy 及我对polya的总结

最新推荐文章于 2025-02-14 14:07:05 发布

小悟空

最新推荐文章于 2025-02-14 14:07:05 发布

阅读量2.2k

点赞数 4

分类专栏：数学 ACMer 文章标签：优化扩展图形 fun c

本文链接：https://blog.youkuaiyun.com/wukonwukon/article/details/7215467

版权

ACMer 同时被 2 个专栏收录

70 篇文章

订阅专栏

数学

17 篇文章

订阅专栏

hdu 2865 Birthday Toy 及我对polya的总结

    一直想总结一下这两天学的东西，今天借这个题总结一下。
   正如上篇所说的：组合计数问题中经常遇到两种困难，第一找出问题通解的表达式，第二是区分讨论问题中哪些应该看成相同的。换句话说，我们就可以将polya 问题分成两部分来分析，从代码上来说，我们也可以分成两部分来分析不同的实现。
    从区分哪些是相同的问题上分析题目，也就是置换群循环节之类的东西。这里分析的时候就不提了。不过就是旋转嘛！
    这个题的重点是如何计数相同珠子不相邻的方案数。
   分析；首先，由于中间一个大圆与每个小圆都相连，所以大圆用去一种颜色之后，只剩下K-1种颜色。
           设K-1种颜色染N个珠子的不同方案数为M，最后就是求M×K mod 1000000007。
           方法跟pku 2888一样，但是这次矩阵的规模很大，所以不能用矩阵来存这个图形了。
           但由于此处规定相邻珠子颜色不同，则邻接阵为对角线上元素全为0,，其余元素全为1。
           该矩阵的幂的迹，可以推导出公式 ( p - 1 ) ^ n + ( -1 ) ^ n * ( p - 1 )其中p是矩阵的阶数，也就是K-1。
    这个公式是怎么求出来的呢？？？？？？
    几乎所有的日志中都是这个公式，但是没见到有解释怎么求的这个公式，我说说我的想法：
    假设A的n-1次幂为：

其中x_n是对角线上的值。乘以对角线上全0，其余为1的矩阵后。
                                                                                 则1)    x_n = y_n-1*(p-1);
2)    y_n = x_n-1+y_n-1*(p-2);
上面两个式子可以解出来x_n = (p-2)*x_n-1 + (p-1)*x_n-2;
事实上，到这一步就能解了，利用矩阵的乘法，然后快速求出x_n的值，进而求出矩阵幂的迹。
当然到这一步，并没有推导出前面的那个公式，上面的递推公式怎么解呢？
注意到：x_n+x_n-1 = (p-1)*(x_n-1+x_n-2) ;
这样的话就能解出x_n+x_n-1;
接下来解出x_n不是问题了。

至此，我们解决了计数问题的两个难题。

我还想总结一下polya上得几个小的知识点。
一：大数的乘法逆元
        大数的乘法逆元我看到了三种方法

暴力法：

 int i;
    for (i=1;;i++) { if (((long long int)(n)*i-an)%M==0) break; }

欧拉函数利用了欧拉函数

long long inv( long long n )
{
    return pow( n, M - 2 )%M;
}

扩展欧几里得解同于方程

//扩展欧几里德
void exp_gcd( LL a ,LL b,LL &x,LL &y) {
     if( b == 0 ) {
         x = 1;
         y = 0;
     }
     else {
          exp_gcd( b,a%b,x,y );
          LL t;
          t = x;
          x = y;
          y = t - a/b*y;
     }
}
//逆元
inline LL getNN(LL x) {
        LL now , t;
        exp_gcd( x, M,now,t );
        return (now%M+M)%M;
}

4.今天看到一个求逆元的简洁写法

int64 inv(int64 x) {  
    //简洁版求逆元  
    if(x == 1) return 1;  
    return  inv(MOD%x) * (MOD - MOD/x) % MOD;  
}

原理还没有弄明白，应该是扩展欧几里得吧，我猜

第一种最慢！

这个解题报告已经过去好久了，但是我还是要对这个地方进行补充，众所周知，求a的逆元的时候，a和m互质。但是有过不呢？？如卡特兰数要是取模呢？看一个链接，讲的很好：http://hi.baidu.com/spellbreaker/blog/item/3b04ec27923ed91f8b82a11e.html

对于polya实现时也存在两种方法

循环时间复杂度O(sqrt(n))

long long ans = 0;
    for (long long i = 1; i*i <= n; i++) if (n % i == 0) { //这个地方可以优化，i*i<=n
        ans = (ans + (geter(p,i) * euler(n / i))%M ) % M;
        if(i*i != n) ans = (ans + (geter(p,n/i) * euler(i))%M ) % M;
    }

第二种是质因数分解，dfs枚举质因数

void dfs(int step, int now, int n) {
    if (step >= cnt + 1) {
        ans = (ans + fun(n % mod, now - 1) * euler(n / now) % mod) % mod;
        return;
    }
    for (int i = 0, t = 1; i <= c[step]; t *= p[step], i++)
        dfs(step + 1, now * t, n);
}

上面的第二种方法时间复杂度低，但是实现起来也复杂。

最后给出我hdu 2865的代码

#include <iostream>
#include <stdio.h>
#include <string.h>
using namespace std;

const long long M = 1000000007;
long long n,k;

#define PP 10000100
long long prime[PP];
bool isPrime[PP+1];
int size;
void getPrime() {
    memset(isPrime,true,sizeof(isPrime));
    int i;
    for(i=2; i<=PP/2; i++) {
        if(isPrime[i])  //i是素数
            for(int j=i+i; j<=PP; j+=i)
                isPrime[j]=0;
    }
    for(i=2; i<=PP; i++) {
        if(isPrime[i])
            prime[size++]=i;
    }
}

long long euler(long long n)//求欧拉函数
{
    long long i,l,t;
    l=n;
    for (i=0;i<size;i++)
    {
        t=0;
        while (l%prime[i]==0) { t++; l=l/prime[i]; }
        if (t>0) n=n/prime[i]*(prime[i]-1);
        if (l==1) break;
        if (l/prime[i]<prime[i]) { n=n/l*(l-1); break; }
    }
    return n%M;
}

long long pow(long long a,long long n)
{
    long long ret=1;
    long long A=a;
    while(n)
    {
        if (n & 1)
            ret=(ret*A)%M;
        A=(A*A)%M;
        n>>=1;
    }
    return ret%M;
}

long long geter(long long p,long long n)
{

    long long ans = 0;
    ans = pow(p-1,n);
    if(n%2 == 0) ans = (ans + p-1)%M;
    else ans = (ans + M - (p-1) )%M;
    return ans;
}

long long inv( long long n )
{
    return pow( n, M - 2 )%M;
}

long long polya(long long p,long long n)
{
    long long ans = 0;
    for (long long i = 1; i*i <= n; i++) if (n % i == 0) { //这个地方可以优化，i*i<=n
        ans = (ans + (geter(p,i) * euler(n / i))%M ) % M;
        if(i*i != n) ans = (ans + (geter(p,n/i) * euler(i))%M ) % M;
    }
    //乘法逆元！！
    return (ans*inv(n))%M;

}

int main()
{
    getPrime();
    while(scanf("%I64d%I64d",&n,&k) != EOF)
        printf("%I64d\n",((long long)k*polya(k-1,n)%M) );
    return 0;
}