当前位置:首页>编程日记>正文

【bzoj3601】一个人的数论 莫比乌斯反演+莫比乌斯函数性质+高斯消元

Description

【bzoj3601】一个人的数论 莫比乌斯反演+莫比乌斯函数性质+高斯消元 配图01

Sol

这题好难啊QAQ

反正不看题解我对自然数幂求和那里是一点思路都没有qwq

先推出一个可做一点的式子:

(f(n)=sum_{k=1}^{n}[(n,k)=1]k^d)

(=sum_{k=1}^{n}k^dsum_{e|n,e|k}mu(e))

(=sum_{e|n}sum_{k=1}^{n/e}(ek)^dmu(e))

(=sum_{e|n}e^dmu(e)sum_{k=1}^{n/e}k^d)

我们假装(反正就是可以但是我太弱了不会证)后面的式子是一个d+1次的关于n的多项式,因为d很小所以我们用高斯消元求出来系数ai,之后得到:

原式

(=sum_{e|n}e^dmu(e)sum_{k=1}^{d+1}a_k(n/e)^k)

(=sum_{k=1}^{d+1}a_ksum_{e|n}e^dmu(e)(n/e)^k)

(fk(n)=sum_{e|n}e^dmu(e)(n/e)^k)

因为(e^dmu(e))((n/e)^i)都是积性函数,所以他们的狄利克雷卷积(fk(n))也是积性函数

我们对于n分解质因数,对于每种质数的qi次幂单独计算,然后乘起来。

显然这样答案等于

(1^d*mu(1)*pi^{qi^{k}}+pi^d*mu(pi)*p^{(qi-1)^{k}})

题目都帮你分解好了就是一种暗示qwq

看到乘数里面有莫比乌斯函数,就要往次幂上想,然后只保留零次幂和一次幂(逃

这个可以直接暴力,ai还是已知的,那么直接枚举d即可。

Code

#include <bits/stdc++.h>
#define ll long long
using namespace std;
ll ans,v[105][105],pa[1005],pb[1005],P=1e9+7;int n,d;
ll ksm(ll a,ll b){ll res=1;for(;b;b>>=1,a=a*a%P) if(b&1) res=res*a%P;return res;}
void gauss(int n)
{for(int i=1,k,t;i<=n;i++){for(int j=i;j<=n;j++) if(v[j][i]) for(k=i;k<=n+1;k++) swap(v[j][k],v[i][k]);for(int j=i+1;j<=n;j++) if(i!=j) for(t=1ll*v[j][i]*ksm(v[i][i],P-2)%P,k=i;k<=n+1;k++) v[j][k]=(v[j][k]-1ll*t*v[i][k]%P+P)%P;}for(int i=n;i;i--) {for(int j=n;j>i;j--) v[i][n+1]=(v[i][n+1]-1ll*v[i][j]*v[j][n+1]%P+P)%P;v[i][n+1]=1ll*v[i][n+1]*ksm(v[i][i],P-2)%P;}
}
int main()
{scanf("%d%d",&d,&n);for(int i=1;i<=d+1;i++){for(int j=1;j<=d+1;j++) v[i][j]=ksm(i,j);for(int j=1;j<=i;j++) v[i][d+2]=(v[i][d+2]+ksm(j,d))%P;}gauss(d+1);for(int i=1;i<=n;i++) scanf("%lld%lld",&pa[i],&pb[i]);for(int i=1,tmp=1;i<=d+1;i++,tmp=1){for(int j=1;j<=n;j++) tmp=1ll*tmp*(ksm(pa[j],pb[j]*i)-ksm(pa[j],d+(pb[j]-1)*i)%P+P)%P;ans=(ans+1ll*tmp*v[i][d+2]%P)%P;}printf("%lld
",ans);
}

http://www.coolblog.cn/news/6fe1b9bcb135dd83.html

相关文章:

  • asp多表查询并显示_SpringBoot系列(五):SpringBoot整合Mybatis实现多表关联查询
  • s7day2学习记录
  • 【求锤得锤的故事】Redis锁从面试连环炮聊到神仙打架。
  • 矿Spring入门Demo
  • 拼音怎么写_老师:不会写的字用圈代替,看到孩子试卷,网友:人才
  • Linux 实时流量监测(iptraf中文图解)
  • Win10 + Python + GPU版MXNet + VS2015 + RTools + R配置
  • 美颜
  • shell访问php文件夹,Shell获取某目录下所有文件夹的名称
  • 如何优雅的实现 Spring Boot 接口参数加密解密?
  • LeCun亲授的深度学习入门课:从飞行器的发明到卷积神经网络
  • Mac原生Terminal快速登录ssh
  • java受保护的数据与_Javascript类定义语法,私有成员、受保护成员、静态成员等介绍...
  • mysql commit 机制_1024MySQL事物提交机制
  • 支撑微博千亿调用的轻量级RPC框架:Motan
  • jquery 使用小技巧
  • 2019-9
  • 法拉利虚拟学院2010 服务器,法拉利虚拟学院2010
  • vscode pylint 错误_将实际未错误的py库添加到pylint白名单
  • 科学计算工具NumPy(3):ndarray的元素处理
  • 工程师在工作电脑存 64G 不雅文件,被公司开除后索赔 41 万,结果…
  • linux批量创建用户和密码
  • newinsets用法java_Java XYPlot.setInsets方法代碼示例
  • js常用阻止冒泡事件
  • 气泡图在开源监控工具中的应用效果
  • 各类型土地利用图例_划重点!国土空间总体规划——土地利用
  • php 启动服务器监听
  • dubbo简单示例
  • 【设计模式】 模式PK:策略模式VS状态模式
  • [iptables]Redhat 7.2下使用iptables实现NAT
  • Ubuntu13.10:[3]如何开启SSH SERVER服务
  • CSS小技巧——CSS滚动条美化
  • JS实现-页面数据无限加载
  • 阿里巴巴分布式服务框架 Dubbo
  • 最新DOS大全
  • Django View(视图系统)
  • 阿里大鱼.net core 发送短信
  • 程序员入错行怎么办?
  • 两张超级大表join优化
  • 第九天函数
  • Linux软件安装-----apache安装
  • HDU 5988 最小费用流
  • Sorenson Capital:值得投资的 5 种 AI 技术
  • 《看透springmvc源码分析与实践》读书笔记一
  • 正式开课!如何学习相机模型与标定?(单目+双目+鱼眼+深度相机)
  • Arm芯片的新革命在缓缓上演
  • nagios自写插件—check_file
  • python3 错误 Max retries exceeded with url 解决方法
  • 行为模式之Template Method模式
  • 通过Spark进行ALS离线和Stream实时推荐