博客
关于我
神犇和蒟蒻
阅读量:306 次
发布时间:2019-03-04

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

题目

题目描述

给定 n ( 1 ≤ n ≤ 1 0 9 ) n(1\le n\le 10^9) n(1n109) ,求下列两个值取模 1 0 9 + 7 10^9+7 109+7 后的结果。

∑ i = 1 n μ ( i 2 ) ,    ∑ i = 1 n φ ( i 2 ) \sum_{i=1}^{n}\mu(i^2),\;\sum_{i=1}^{n}\varphi(i^2) i=1nμ(i2),i=1nφ(i2)

思路

对于第一个,直接输出 1 1 1 就好了……定义告诉我们,指数为 2 2 2 就会凉。只有 μ ( 1 2 ) = 1 \mu(1^2)=1 μ(12)=1 可以挽回一些局面。

对于第二个,不难发现 φ ( i 2 ) = i ⋅ φ ( i ) \varphi(i^2)=i\cdot\varphi(i) φ(i2)=iφ(i) 。证明方法很多,我给两个比较好理解的,否则这篇博客就太短啦

等价转化

φ ( x 2 ) = ∑ i = 1 x 2 [ gcd ⁡ ( i , x 2 ) = 1 ] = ∑ i = 1 x 2 [ gcd ⁡ ( i , x ) = 1 ] = ∑ i = 1 x 2 [ gcd ⁡ ( i   m o d   x , x ) = 1 ] = ∑ d = 1 x [ gcd ⁡ ( d , x ) = 1 ] ∑ i ≡ d ( m o d x ) 1 = ∑ d = 1 x [ gcd ⁡ ( d , x ) = 1 ] ⋅ x = x ⋅ φ ( x ) \begin{aligned} \varphi(x^2)&=\sum_{i=1}^{x^2}[\gcd(i,x^2)=1]\\ &=\sum_{i=1}^{x^2}[\gcd(i,x)=1]\\ &=\sum_{i=1}^{x^2}[\gcd(i\bmod x,x)=1]\\ &=\sum_{d=1}^{x}[\gcd(d,x)=1]\sum_{i\equiv d\pmod{x}}1\\ &=\sum_{d=1}^{x}[\gcd(d,x)=1]\cdot x\\ &=x\cdot\varphi(x) \end{aligned} φ(x2)=i=1x2[gcd(i,x2)=1]=i=1x2[gcd(i,x)=1]=i=1x2[gcd(imodx,x)=1]=d=1x[gcd(d,x)=1]id(modx)1=d=1x[gcd(d,x)=1]x=xφ(x)

狄利克雷

φ ( x 2 ) = ∑ d ∣ x 2 μ ( d ) ⋅ x 2 d = x ∑ d ∣ x μ ( d ) ⋅ x d = x ⋅ φ ( x ) \begin{aligned} \varphi(x^2)&=\sum_{d|x^2}\mu(d)\cdot {x^2\over d}\\ &=x\sum_{d|x}\mu(d)\cdot{x\over d}\\ &=x\cdot \varphi(x) \end{aligned} φ(x2)=dx2μ(d)dx2=xdxμ(d)dx=xφ(x)

用上 φ = μ ∗ I d \varphi=\mu*I_d φ=μId 就非常好推呀!并且这个方向也比较自然。

当然, d ∣ x 2 d|x^2 dx2 变为 d ∣ x d|x dx 也是比较重要的操作。

最后一步,直接杜教筛即可。用 I d I_d Id 去筛。

代码

#include 
#include
#include
#include
using namespace std;typedef long long int_;inline int readint(){ int a = 0; char c = getchar(), f = 1; for(; c<'0'||c>'9'; c=getchar()) if(c == '-') f = -f; for(; '0'<=c&&c<='9'; c=getchar()) a = (a<<3)+(a<<1)+(c^48); return a*f;}const int MaxN = 5e6;const int Mod = 1e9+7;const int inv2 = (Mod+1)>>1;const int inv3 = (Mod+1)/3;vector< int > primes;bool isPrime[MaxN];int phi[MaxN], f[MaxN];void sievePrime(){ memset(isPrime+2,1,MaxN-2); f[1] = phi[1] = 1; for(int i=2,len=0; i
= MaxN) break; // good try isPrime[x] = false; if(i%primes[j] == 0){ phi[x] = phi[i]*primes[j]; break; } phi[x] = phi[i]*(primes[j]-1); } } for(int i=2; i
= SqrtN) ? haxi[1][n/(x)] : haxi[0][x])int w[SqrtN<<1]; // f(xez[x])int_ xez[SqrtN<<1]; // from haxi to realint djs(int_ n){ // Dc.Du sieve int id = 0; // allocate index for(int_ i=1; i<=n; i=n/(n/i)+1){ index_(n/i) = ++ id; xez[id] = n/i; // real value } int_ x; // for convenience for(int i=id; i>=1; --i){ if((x = xez[i]) < MaxN){ w[i] = f[x]; continue; } w[i] = ((x<<1)+1)*x%Mod *(x+1)%Mod*inv2%Mod *inv3%Mod; // sum x^2 int lst = 1, now; // init 1 for(int_ l=2,r; l<=x; l=r+1){ r = x/(x/l); now = (r*(r+1)>>1)%Mod; w[i] -= 1ll*(now-lst) *w[index_(x/l)]%Mod; w[i] = (w[i]%Mod+Mod)%Mod; lst = now; // move on } }// for(int i=id; i>=1; --i)// printf("f[%lld] = %d\n",xez[i],w[i]); return w[1]; // maximum = n}int main(){ int_ n; scanf("%lld",&n); sievePrime(); printf("1\n%d\n",djs(n)); return 0;}

转载地址:http://ovtq.baihongyu.com/

你可能感兴趣的文章
Java基础之反射
查看>>
线程池之SingleThreadPool学习
查看>>
对象的创建、内存布局和访问定位
查看>>
FreeRTOS学习笔记(9)——内存管理
查看>>
FreeRTOS学习笔记(10)——中断管理
查看>>
CC2640R2F学习笔记(1)——搭建环境、编译烧写
查看>>
ESP8266学习笔记(10)——官方WebServer
查看>>
CC2640R2F学习笔记(6)——UART串口使用
查看>>
SHELL命令
查看>>
redis命令学习
查看>>
自然划分的3-4-5规则
查看>>
剑指offer Leetcode 37.序列化二叉树
查看>>
剑指offer Leetcode 39.数组中出现次数超过一半的数字
查看>>
Latex中cases环境引入报错
查看>>
Latex排版的时候把图片放在指定位置
查看>>
用 Python 把你的朋友变成表情包(鼠标事件提取 ROI 版)
查看>>
Tensorflow2.0:基于循环卷积网络预测剩余寿命
查看>>
bzoj3879: SvT 后缀自动机
查看>>
4084: [Sdoi2015]双旋转字符串
查看>>
bzoj3439: Kpm的MC密码(四种做法)
查看>>