bzoj3601: 一个人的数论 莫比乌斯反演 高斯消元

bzoj3601: 一个人的数论

题目

在这里插入图片描述

分析

数论神仙题,话不多说开始推导
a n s = ∑ x n [ g c d ( x , n ) = = 1 ] x d ans=\sum_x^n[gcd(x,n)==1]x^d ans=xn[gcd(x,n)==1]xd
= ∑ x n ∑ c ∣ n , c ∣ x μ ( c ) x d =\sum_x^n\sum_{c|n,c|x}\mu(c)x^d =xncn,cxμ(c)xd
= ∑ c ∣ n n μ ( c ) ∑ x n c ( c x ) d =\sum_{c|n}^n\mu(c)\sum_{x}^{\frac{n}{c}}(cx)^d =cnnμ(c)xcn(cx)d
= ∑ c ∣ n n μ ( c ) c d ∑ x n c x d =\sum_{c|n}^n\mu(c)c^d\sum_{x}^{\frac{n}{c}}x^d =cnnμ(c)cdxcnxd

结论:等幂和

∑ i n i d = ∑ i d a i n i \sum_i^n i^d=\sum_i^d a_in^i inid=idaini

也就是 d d d次等幂和是关于 n n n d + 1 d+1 d+1次多项式,可以采用高斯消元搞出这个玩意儿。(推导?插值啊!!
a n s = ∑ c ∣ n n μ ( c ) c d ∑ i d + 1 a i ( n c ) i ans=\sum_{c|n}^n\mu(c)c^d\sum_i^{d+1}a_i(\frac{n}{c})^i ans=cnnμ(c)cdid+1ai(cn)i
= ∑ i d + 1 a i ∑ c ∣ n n μ ( c ) c d ( n c ) i =\sum_i^{d+1}a_i\sum_{c|n}^n\mu(c)c^d(\frac{n}{c})^i =id+1aicnnμ(c)cd(cn)i
g i ( n ) = ∑ c ∣ n n μ ( c ) c d ( n c ) i g_i(n)=\sum_{c|n}^n\mu(c)c^d(\frac{n}{c})^i gi(n)=cnnμ(c)cd(cn)i
显然这玩意儿是个积性函数,因为他本身是两个积性函数的狄利克雷卷积。
所以我们只需要研究 g i ( p q ) g_i(p^q) gi(pq),然后乘起来。
g i ( p q ) = ∑ j = 0 q μ ( p j ) p j d ( p q p j ) i g_i(p^q)=\sum_{j=0}^{q}\mu(p^j)p^{jd}(\frac{p^q}{p^j})^i gi(pq)=j=0qμ(pj)pjd(pjpq)i
当且仅当 j = 0 , 1 j=0,1 j=0,1 μ ( p j ) ≠ 0 \mu(p^j)\neq 0 μ(pj)̸=0
所以有
g i ( p q ) = p q i − c q i + d − i g_i(p^q)=p^{qi}-c^{qi+d-i} gi(pq)=pqicqi+di
呼呼,解决了!!

代码

#include<cstdio>
#include<algorithm>
const int P = 1e9 + 7, N = 1e3 + 10;
int ri() {
    char c = getchar(); int x = 0; for(;c < '0' || c > '9'; c = getchar()) ;
    for(;c >= '0' && c <= '9';  c = getchar()) x = (x << 1) + (x << 3) - '0' + c; return x;
}
int a[N][N], p[N], q[N], d, w;
int Pow(int x, int k) {
    k < 0 ? k += P - 1 : 0; int r = 1;
    for(;k; x = 1LL * x * x % P, k >>= 1) if(k & 1) r = 1LL * r * x % P;
    return r;
}
void Gas() {
    for(int i = 0;i <= d + 1; ++i) {
        a[i][0] = 1;
        for(int j = 1;j <= d + 1; ++j) a[i][j] = 1LL * a[i][j - 1] * (i + 1) % P;
        if(i) a[i][d + 2] = a[i - 1][d + 2];
        a[i][d + 2] = (a[i][d + 2] + a[i][d]) % P;
    }
    for(int i = 0, j; i <= d + 1; ++i) {
        for(j = i; j <= d + 1; ++j) if(a[j][i]) break;
        if(i != j) std::swap(a[i], a[j]);
        for(j = i + 1; j <= d + 1; ++j) {
            int t = 1LL * a[j][i] * Pow(a[i][i], P - 2) % P;
            for(int k = i; k <= d + 2; ++k)
                a[j][k] = (a[j][k] + P - 1LL * a[i][k] * t % P) % P;
        }
    }
    for(int i = d + 1; ~i; --i) {
        for(int j = i + 1; j <= d + 1; ++j)
            a[i][d + 2] = (a[i][d + 2] + P - 1LL * a[i][j] * a[j][d + 2] % P) % P;
        a[i][d + 2] = 1LL * a[i][d + 2] * Pow(a[i][i], P - 2) % P;
    }
}
int main() {
    d = ri(); w = ri(); Gas(); int ans = 0;
    for(int i = 1;i <= w; ++i) p[i] = ri(), q[i] = ri();
    for(int i = 0, r = 1; i <= d + 1; ans = (ans + 1LL * a[i++][d + 2] * r) % P, r = 1)
        for(int j = 1;j <= w; ++j)
            r = 1LL * r * Pow(p[j], 1LL * q[j] * i % (P - 1)) % P * (1 + P - Pow(p[j], d - i)) % P;
    printf("%d\n", ans);
    return 0;
}
  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值