3512: DZY Loves Math IV
题意：求(sum_{i=1}^n sum_{j=1}^m varphi(ij))，(n le 10^5, m le 10^9)
n较小，考虑写成前缀和的形式，计算(S(n,m)=sum_{i=1}^m varphi(in))
一开始想出
[n= prod_i p_i, varphi(in) = varphi(i) cdot varphi(frac{n}{d})cdot d, d=(n,i)
]
比较好想，共有的质因子应该乘(p)而不是(p-1)
然后带进去枚举gcd用莫比乌斯反演的套路，中间的函数很奇怪不好算前缀和...
orz了题解，发现题解使用(varphi * 1 =id)来替换
[varphi(in) = varphi(i) cdot varphi(frac{n}{d})cdot sum_{emid d} varphi(e) = varphi(i) cdot sum_{emid d}varphi(frac{n}{e})
]
因为n是不同质因子的乘积，所以可以把两个(varphi)乘起来
这一步替换和用(mu * 1 = epsilon)替换有异曲同工之妙，都是将(gcd)等于的限制弱化了，变成了整除的关系
推倒后得到
[S(n,m) = sum_{dmid n}varphi(frac{n}{d})cdot S(d, lfloor frac{m}{d}
floor)
]
对于n不是不同质因子的乘积的，根据(varphi)的公式，多的质因子次数直接提出来乘上就行了
然后记忆化搜索，(n=1)就是(varphi)的前缀和
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <map>
using namespace std;
typedef long long ll;
const int N=1664512, U=1664510, mo = 1e9+7;
inline int read(){
char c=getchar(); int x=0,f=1;
while(c<'0' || c>'9') {if(c=='-')f=-1; c=getchar();}
while(c>='0' && c<='9') {x=x*10+c-'0'; c=getchar();}
return x*f;
}
int n, m;
inline void mod(int &x) {if(x>=mo) x-=mo; else if(x<0) x+=mo;}
bool notp[N]; int p[N/10], phi[N], pr[N];
void sieve(int n) {
phi[1]=1; pr[1] = 1;
for(int i=2; i<=n; i++) {
if(!notp[i]) p[++p[0]] = i, phi[i] = i-1, pr[i] = i;
for(int j=1; j <= p[0] && i*p[j] <= n; j++) {
int t = i*p[j];
notp[ t ] = 1;
if(i % p[j] == 0) {
phi[t] = phi[i] * p[j];
pr[t] = pr[i];
break;
}
phi[t] = phi[i] * (p[j] - 1);
pr[t] = pr[i] * p[j];
}
mod(phi[i] += phi[i-1]);
}
}
namespace ha {
const int p = 1001001;
struct meow{int ne, val, r;} e[3000];
int cnt=1, h[p];
inline void insert(int x, int val) {
int u = x % p;
for(int i=h[u];i;i=e[i].ne) if(e[i].r == x) return;
e[++cnt] = (meow){h[u], val, x}; h[u] = cnt;
}
inline int quer(int x) {
int u = x % p;
for(int i=h[u];i;i=e[i].ne) if(e[i].r == x) return e[i].val;
return -1;
}
} using ha::insert; using ha::quer;
int dj_s(int n) {
if(n <= U) return phi[n];
if(quer(n) != -1) return quer(n);
int ans = (ll) n * (n+1) / 2 %mo, r;
for(int i=2; i<=n; i=r+1) {
r = n/(n/i);
mod(ans -= (ll) dj_s(n/i) * (r-i+1) %mo);
}
insert(n, ans);
return ans;
}
inline int Pow(int a, int b) {
int ans=1;
for(; b; b>>=1, a=a*a)
if(b&1) ans=ans*a;
return ans;
}
inline ll Phi(int n) {
int ans = 1;
if(n <= U) {mod(ans = phi[n] - phi[n-1]); return ans;}
int m = sqrt(n);
for(int i=1; p[i] <= m; i++) if(n % p[i] == 0) {
int a = 0;
while(n % p[i] == 0) a++, n /= p[i];
ans *= Pow(p[i], a-1) * (p[i] - 1);
}
return ans;
}
map<int, int> Map[N];
int S(int n, int m) {
if(m == 0) return 0;
if(n == 1) return dj_s(m);
if(Map[n][m]) return Map[n][m];
//printf("S %d %d
", n, m);
int ans = 0;
for(int i=1; i*i <= n; i++) if(n%i == 0) {
int j = n/i;
mod(ans += Phi(j) * S(i, m/i) %mo);
if(i != j) mod(ans += Phi(i) * S(j, m/j) %mo);
}
Map[n][m]=ans;
return ans;
}
int main() {
freopen("in", "r", stdin);
sieve(U);
n=read(); m=read();
int ans = 0;
for(int i=1; i<=n; i++) mod(ans += (ll) i / pr[i] * S(pr[i], m) %mo);
//for(int i=1; i<=n; i++) printf("nnnnnnnn %d %d
", i, S(i, m));
printf("%d
", ans);
}