B - bookshelf HDU - 6363 ——莫比乌斯反演

最新推荐文章于 2022-02-11 22:39:38 发布

天翼之城*

最新推荐文章于 2022-02-11 22:39:38 发布

阅读量216

点赞数

分类专栏：数学文章标签：莫比乌斯反演

本文链接：https://blog.csdn.net/tianyizhicheng/article/details/81673061

版权

数学专栏收录该内容

91 篇文章 1 订阅

订阅专栏

Patrick Star bought a bookshelf, he named it ZYG !!

Patrick Star has N book .

The ZYG has K layers (count from 1 to K) and there is no limit on the capacity of each layer !

Now Patrick want to put all N books on ZYG :

Assume that the i-th layer has cnti(0≤cnti≤N) books finally.
Assume that f[i] is the i-th fibonacci number (f[0]=0,f[1]=1,f[2]=1,f[i]=f[i−2]+f[i−1]).
Define the stable value of i-th layers stablei=f[cnti].
Define the beauty value of i-th layers beautyi=2stablei−1.
Define the whole beauty value of ZYG score=gcd(beauty1,beauty2,…,beautyk)(Note: gcd(0,x)=x).

Patrick Star wants to know the expected value of score if Patrick choose a distribute method randomly !
Input
The first line contain a integer T (no morn than 10), the following is T test case, for each test case :

Each line contains contains three integer n,k(0< n,k≤106).
Output
For each test case, output the answer as a value of a rational number modulo 109+7.

Formally, it is guaranteed that under given constraints the probability is always a rational number pq (p and q are integer and coprime, q is positive), such that q is not divisible by 109+7. Output such integer a between 0 and 109+6 that p−aq is divisible by 109+7.
Sample Input
1
6 8
Sample Output
797202805

第一次接触莫比乌斯反演，挺有趣的，也比较难，我不是负责数论的，但是因为有空就去学了一下
这篇博客主要讲的是反演的系数推导，
挺有帮助的 https://www.cnblogs.com/chenyang920/p/4811995.html

有两个需要知道的东西，我也是看了别人的博客才知道有这么个玩意
$gcd(x^a-1,x^b-1)==x^{gcd(a,b)}-1$
和
$g c d (f b [i], f b [j]) = f b [g c d (i, j)]$
之后根据要求的d求出它的倍数并构造函数，具体的看别的博客吧QAQ，因为2^f[n]很爆炸，所以要欧拉降幂
1e9+7是质数所以他的欧拉函数就是mod-1（我猜的）；

#include<bits/stdc++.h>
using namespace std;
#define ll long long
const ll mod=1e9+7;
const int maxn=2e6+5;
int n,k;
int nprime[maxn],prime[maxn],mu[maxn],cnt;
ll jc[maxn],ny[maxn],fb[maxn];
ll qpow(ll a,ll b)
{
    ll ret=a;
    ll ans=1;
    while(b)
    {
        if(b&1)
            ans=(ans*ret)%mod;
        ret=(ret*ret)%mod;
        b>>=1;
    }
    return ans;
}
void init()
{
    jc[0]=1;
    for(ll i=1;i<maxn;i++)
        jc[i]=(jc[i-1]*i)%mod;
    ny[maxn-1]=qpow(jc[maxn-1],mod-2);
    for(ll i=maxn-2;i>=0;i--)
        ny[i]=ny[i+1]*(i+1)%mod;


    cnt=0;
    mu[1]=1;
    for(ll i=2;i<maxn;i++)
    {
        if(!nprime[i])
        {
            prime[++cnt]=i;
            mu[i]=-1;
        }
        for(ll j=1;j<=cnt&&prime[j]*i<maxn;j++)
        {
            nprime[prime[j]*i]=1;
            if(i%prime[j]==0)
            {
                mu[prime[j]*i]=0;
                break;
            }
            mu[prime[j]*i]=-mu[i];
        }
    }

    fb[1]=fb[2]=1;
    for(int i=3;i<maxn;i++)
        fb[i]=(fb[i-1]+fb[i-2])%(mod-1);
}

ll c(ll n,ll m)
{
    if(n<0||m<0||m>n)
        return 0;
    if(n==m||m==0)
        return 1;
    return jc[n]*ny[n-m]%mod*ny[m]%mod;
}

int main()
{
    int t;
    init();
    scanf("%d",&t);
    while(t--)
    {
        ll ans=0;
        scanf("%d%d",&n,&k);
        for(ll i=1;i<=n;i++)
        {
            if(n%i)
                continue;
            ll ret=0;
            for(ll j=i;j<=n;j+=i)
            {
                if(n%j==0)
                {
                    ret=(ret+mu[j/i]*c(n/j+k-1,k-1)+mod)%mod;
                }
            }
            ans=(ans+ret*(qpow(2ll,fb[i])-1+mod))%mod;
            //这里可以+mod-1但是不加也行可能是因为fb取余的就是mod-1吧
        }
        ll a=c(n+k-1,k-1);
        a=qpow(a,mod-2);
        ans=ans*a%mod;
        printf("%lld\n",ans);
    }
    return 0;
}