博客
关于我
强烈建议你试试无所不能的chatGPT,快点击我
【BZOJ3601】一个人的数论 - 莫比乌斯反演+高斯消元
阅读量:5088 次
发布时间:2019-06-13

本文共 2712 字,大约阅读时间需要 9 分钟。

咕了inf天终于来补这题了……

题意:

题解:

常规操作,先套莫比乌斯反演大力推式子:

$$f_d(n)=\sum\limits_{i=1}^{n}i^d[gcd(i,n)=1]$$

$$=\sum\limits_{i=1}^{n}i^d\sum\limits_{k|gcd(i,n)}\mu(k)$$

$$=\sum\limits_{k|n}\sum\limits_{i=1}^{\frac{n}{k}}(ik)^d\mu(k)$$

$$=\sum\limits_{k|n}k^d\mu(k)\sum\limits_{i=1}^{\frac{n}{k}}i^d$$

设$S(n)=\sum\limits_{i=1}^{n}i^d$,则:

$$ans=\sum\limits_{k|n}k^d\mu(k)S(\frac{n}{k})$$

经过一些微小的数学分析我们可知$S(n)$是一个关于$n$的$d+1$次多项式;

证明详见那题……

于是可以暴力高斯消元或者拉格朗日插值来求出这个多项式的每一项系数;

设$S(n)=\sum\limits_{i=1}^{d+1}s_in^i$,则:

$$ans=\sum\limits_{k|n}k^d\mu(k)\sum\limits_{i=1}^{d+1}s_i(\frac{n}{k})^i$$

$$=\sum\limits_{i=1}^{d+1}s_i\sum\limits_{k|n}\mu(k)k^d(\frac{n}{k})^i$$

$$=\sum\limits_{i=1}^{d+1}s_i\sum\limits_{k|n}\mu(k)k^{d-i}n^i$$

设$g_i(n)=\sum\limits_{k|n}\mu(k)k^{d-i}n^i$,易得这是个积性函数;

题目中也很良心的给出的是$n$的质因数分解式,那么可以分开来考虑$n$的所有质数的整数次幂的因数(即$p_i^{\alpha_i}$):

$$ans=\sum\limits_{i=1}^{d+1}s_i\prod\limits_{j=1}^{w}g_i(p_j^{\alpha_i})$$

只看其中的$g_i(p^\alpha)$一项:

$$g_i(p^\alpha)=\sum\limits_{j=0}^{\alpha}\mu(p^j)p^{j\times(d-i)}p^{\alpha i}$$

显然$j>1$时后面的都为0……因此只用考虑前两项即可:

$$g_i(p^\alpha)=p^{\alpha i}-p^{d-i}p^{\alpha i}$$

$$=p^{\alpha i}(1-p^{d-i})$$

于是就做完了。

代码:

1 #include
2 #include
3 #include
4 #include
5 #include
6 #include
7 #define inf 2147483647 8 #define eps 1e-9 9 #define mod 100000000710 using namespace std;11 typedef long long ll;12 typedef double db;13 int d,w,ans=0,tmp,p[1001],a[1001],sq[105][105],s[101];14 int fastpow(int x,int y){15 int ret=1;16 if(y<0)y+=mod-1;17 for(;y;y>>=1,x=(ll)x*x%mod){18 if(y&1)ret=(ll)ret*x%mod;19 }20 return ret;21 }22 void gauss(){23 int sum=0,nw,tmp;24 for(int i=0;i<=d+1;i++){25 sum=(sum+fastpow(i+1,d))%mod;26 sq[i][0]=1;27 sq[i][d+2]=sum;28 for(int j=1;j<=d+1;j++){29 sq[i][j]=(ll)sq[i][j-1]*(i+1)%mod;30 }31 }32 for(int i=0;i<=d+1;i++){33 nw=i;34 for(int j=i;j<=d+1;j++){35 if(sq[i][j]){36 nw=j;37 break;38 }39 }40 if(nw!=i){41 for(int j=0;j<=d+2;j++){42 swap(sq[i][j],sq[nw][j]);43 }44 }45 for(int j=0;j<=d+1;j++){46 if(i!=j&&sq[j][i]){47 tmp=(ll)sq[j][i]*fastpow(sq[i][i],mod-2)%mod;48 for(int k=0;k<=d+2;k++){49 sq[j][k]=(sq[j][k]-(ll)tmp*sq[i][k]%mod+mod)%mod;50 }51 }52 }53 }54 for(int i=0;i<=d+1;i++){55 s[i]=(ll)sq[i][d+2]*fastpow(sq[i][i],mod-2)%mod;56 }57 }58 int main(){59 scanf("%d%d",&d,&w);60 for(int i=1;i<=w;i++){61 scanf("%d%d",&p[i],&a[i]);62 }63 gauss();64 for(int i=1;i<=d+1;i++){65 tmp=1;66 for(int j=1;j<=w;j++){67 tmp=(ll)tmp*fastpow(p[j],(ll)a[j]*i%(mod-1))%mod;68 tmp=(ll)tmp*(mod-fastpow(p[j],d-i)+1)%mod;69 }70 ans=(ans+(ll)tmp*s[i]%mod)%mod;71 }72 printf("%d",ans);73 return 0;74 }

转载于:https://www.cnblogs.com/dcdcbigbig/p/10105039.html

你可能感兴趣的文章
Git随笔 -- 初始化远程仓库
查看>>
学习Unity -- 理解依赖注入(IOC)三种方式依赖注入
查看>>
高效法则 之 你还在用这么low的方法打开软件吗?
查看>>
k火车座位分布图
查看>>
C#集合之STACK
查看>>
容易混淆的 i++ 与 ++i
查看>>
bzoj 4480: [Jsoi2013]快乐的jyy
查看>>
Taxes
查看>>
算法导论(第三版)Exercises2.3(归并排序、二分查找、计算集合中是否有和为X的2个元素)...
查看>>
算法导论(第三版)Problems2(归并插入排序、数列逆序计算)
查看>>
d3.js:数据可视化利器之 交互行为:响应DOM事件
查看>>
微信小程序(18)-- 自定义头部导航栏
查看>>
CSS继承—深入剖析
查看>>
IOS开发中的分享到邮件
查看>>
Resharper插件的使用
查看>>
unity中UI的屏幕自适应代码
查看>>
lagou数据爬取
查看>>
井底飞天
查看>>
<a>标签实现锚点跳跃,<a>标签实现href不跳跃另外加事件(ref传参)
查看>>
C# async/await异步操作:异步执行方法封装
查看>>