HDU 6439(容斥原理+积性函数前缀和+打表）

Congruence equation

Time Limit: 8000/4000 MS (Java/Others)    Memory Limit: 65536/65536 K (Java/Others)
Total Submission(s): 153    Accepted Submission(s): 37

 

Problem Description

There is a sequence A which contains k integers.
Now we define f(m) is the number of different sequence C that satisfies for i from 1 to k:
1. If Ai=−1, Ci can be any integer in the range of [0,m).Otherwise Ci=Ai(modm)
2. ∑i=1kCixi≡1(modm)(xi are variables) have a solution in the range of integer.

Calculate the answer of ∑m=1nf(m)(mod109+7).

 

 

Input

The first line contains only one integer T(T≤100), which indicates the number of test cases. 
For each test case, the first line contains two integers k and n.(1≤k≤50,1≤n≤109)
The second line contains k integers: A1,A2…Ak(−1≤Ai≤109)
There are at most 10 test cases which satisfies n≥106

 

 

 

Output

For each test case, output one line "Case #x: y", where x is the case number (starting from 1) and y is the answer after mod 1000000007 (109+7).

 

 

Sample Input

2 5 10 -1 -1 8 -1 -1 3 20 -1 6 18

 

 

Sample Output

Case #1: 24354 Case #2: 140

 

 

Source

2018中国大学生程序设计竞赛 - 网络选拔赛

 

首先，的时候才能有解

我们记不是-1的a的gcd为g，并记录-1的数量num

如果g不为0，相当于问这num个变量有多少种取法，使得gcd为1

因为最后的gcd一定是g的因子，所以我们只需要对g的因子进行容斥

对于某个因子p,当前的答案就是

考虑到num最多50，可以采用伯努利数的方法来求自然数幂和

如果g=0 就是问这num个变量有多少种取法，使得gcd为1

考虑到如果num=1，答案显然为phi(n)

所以可以尝试打num=2的表，可以总结出规律：

这个函数与欧拉函数类似，显然，这也是个积性函数

之后就是用类似求欧拉函数前n项和的办法来就可以了

只需要改一改线性筛的部分以及将n*(n+1)/2改成1^k+2^k+...+n^k即可

中间需要记忆化，不然会T

#include<bits/stdc++.h>
#define mp make_pair
#define fir first
#define se second
#define ll long long
#define pb push_back
#define LL long long
#pragma comment(linker, "/STACK:1024000000,1024000000")
using namespace std;
const int maxn=2e6+10;
const ll mod=1e9+7;
const ll MOD=mod;
const int maxm=1e6+10;
const double eps=1e-7;
const int inf=0x3f3f3f3f;
const double INF = 1e20;
const double pi = acos (-1.0);
const int N=100;
int gcd(int a,int b){
    return (b==0?a:gcd(b,a%b));
}
ll qpow(ll a,ll b,ll p){
    ll res=1;
    while (b){
        if (b&1) res=res*a%p;
        b>>=1;
        a=a*a%p;
    }
    return res;
}
//伯努利数（自然数幂和的预处理部分）
LL C[N][N];
LL B[N],Inv[N];
LL Tmp[N];
void Init()
{
    //预处理组合数
    for(int i=0; i<N; i++)
    {
        C[i][0] = C[i][i] = 1;
        if(i == 0) continue;
        for(int j=1; j<i; j++)
            C[i][j] = (C[i-1][j] % MOD + C[i-1][j-1] % MOD) % MOD;
    }
    //预处理逆元
    Inv[1] = 1;
    for(int i=2; i<N; i++)
        Inv[i] = (MOD - MOD / i) * Inv[MOD % i] % MOD;
    //预处理伯努利数
    B[0] = 1;
    for(int i=1; i<N; i++)
    {
        LL ans = 0;
        if(i == N - 1) break;
        for(int j=0; j<i; j++)
        {
            ans += C[i+1][j] * B[j];
            ans %= MOD;
        }
        ans *= -Inv[i+1];
        ans = (ans % MOD + MOD) % MOD;
        B[i] = ans;
    }
}

ll k,n;
ll a[maxn];

ll sum(ll n,ll k){
    //用伯努利数算1^k+2^k+..+n^k;
    if (k==0) return n;
    ll res=0,now=n+1;
    for (int i=1;i<=k+1;i++){
        res=(res+C[k+1][i]*B[k+1-i]%mod*now%mod)%mod;
        now=now*(n+1)%mod;
    }
    res=res*Inv[k+1]%mod;
    return res;
}
bool vis[maxn];
ll phi[maxn], prime[maxn], cnt, prime1[maxn];
ll Sum[maxn];
void pre (int Maxn,int k)
{
    cnt = 0;
    for (int i=0;i<=Maxn;i++)
        vis[i]=0;
    phi[1] = 1;
    for ( int i = 2; i <= Maxn; i++)
    {
        if (!vis[i])
        {
            prime[cnt++] = i;
            prime1[i]=qpow(i,k,mod);
            phi[i] = qpow(i,k,mod)-1;
        }
        for ( int j = 0; j < cnt; j++)
        {
            if ( 1LL*i*prime[j] > Maxn)
                break;
            vis[i*prime[j]] = 1;
            if (i%prime[j] == 0)
            {
                phi[i*prime[j]] = phi[i] * prime1[prime[j]]%mod;
                //cout<<i*prime[j]<<" "<<phi[i]<<" "<<prime[j]<<" "<<phi[i*prime[j]]endl;
                break;
            }
            else
            {
                phi[i*prime[j]] = phi[i] * (prime1[prime[j]]-1)%mod;
            }
        }
    }
    Sum[0]=0;
    for (int i=1;i<=Maxn;i++)
        Sum[i]=(Sum[i-1]+phi[i])%mod;
}
map<long long ,long long> ma;
ll calc(long long n,int Maxn,int num){
    if(n<=Maxn) return Sum[n];
    if(ma.find(n)!=ma.end()) return ma[n];
    ll ans=sum(n,num);
    for(long long i=2,r;i<=n;i=r+1){
        r=n/(n/i);
        ans=(ans-1LL*calc(n/i,Maxn,num)*((r-i+1))%mod)%mod;
    }
   return ma[n]=(ans+mod)%mod;
}
int main(){
    int t;
    scanf("%d",&t);
    Init();
    int kase=0;
    while (t--){
        kase++;
        int num=0;
        scanf("%lld %lld",&k,&n);
        for (int i=1;i<=k;i++){
            scanf("%lld",&a[i]);
        }
        ll g=0;
        for (int i=1;i<=k;i++){
            if (a[i]==-1){
                num++;
            }
            else g=gcd(g,a[i]);
        }
        if (g!=0){
            //枚举g的因子 大力容斥
            vector<ll> v;
            ll temp=g;
            for (ll i=2;i*i<=g;i++){
                if (temp%i==0){
                    v.pb(i);
                    while(temp%i==0) temp/=i;
                }
            }
            if (temp>1) v.pb(temp);
            ll ans=0;
            for (int i=0;i<(1<<v.size());i++){
                ll d=1,flag=1;
                for (int  j=0;j<v.size();j++){
                    if (i&(1<<j)){
                        d*=v[j];
                        flag*=-1;
                    }
                }
                if (n>=d){
                    ans=(ans+flag*sum(n/d,num)%mod)%mod;
                }
            }
            ans=(ans+mod)%mod;
            printf("Case #%d: %lld\n",kase,ans);
        }
        else{
            int Maxn=pow(n,2.0/3.0);
            pre(Maxn,num);
            ma.clear();
            ll ans=calc(n,Maxn,num);
            ans=(ans+mod)%mod;
            printf("Case #%d: %lld\n",kase,ans);
        }
    }
    return 0;
}