前言
想感受一下敲公式的快感，于是有了这篇文章。
只有一些最简单的应用。
杜教筛
可以在低于线性的时间内筛积性函数前缀和，下文介绍原理：
假设我们要求 $S(n)=\sum_{i=1}^nf(i)$，其中 $f(x)$ 是个积性函数。
然后对于任意一个积性函数 $g(x)$，都有： $$ \begin{aligned} \sum\limits_{i=1}^n(f*g)(i) =&\sum\limits_{i=1}^n\sum\limits_{d\mid i}g(d)f(i/d)\\ =&\sum\limits_{d=1}^ng(d)\sum\limits_{d\mid i}^nf(i/d)\\ =&\sum\limits_{d=1}^ng(d)S(\lfloor n/d\rfloor)\\ =&\sum\limits_{i=1}^ng(i)S(\lfloor n/i\rfloor) \end{aligned} $$
接下来有：
$$ g(1)S(n)=\sum\limits_{i=1}^n(f*g)(i)-\sum_{i=2}^ng(i)S(\lfloor n/i\rfloor) $$ 需要构造合适的函数 $g(x)$ 来快速求解前一项，同时用整除分块计算后一项。
实现方式是先筛出来这个函数的前 $n^{\frac{2}{3}}$ 项，然后比这个数大的采用递归计算。用哈希表记忆化后的复杂度是 $O(n^{\frac{2}{3}})$，我不会证明。
luogu4213 sum
Description
求 $\sum_{i=1}^n\mu(i),\sum_{i=1}^n\varphi(i)$，$n<2^{31}$。
Solution
$\mu*1=\epsilon$，第一个问题做完了。
第二个虽然可以用 $id=\varphi*1$ 做，但是用莫比乌斯反演更快一些：先求出 $\sum_{i=1}^n\sum_{j=1}^n[\gcd(i,j)=1]$，简单莫反得到这个东西等于 $\sum_{i=1}^n\mu(i)\lfloor n/d\rfloor^2$，求出 $\mu$ 的前缀和之后可以直接 $O(\sqrt n)$ 求解。
然后要求的答案是 $\sum_{i=1}^n\sum_{j=1}^i[\gcd(i,j)=1]$，把刚刚的答案除以 $2$，再考虑一下 $i=j=1$ 就做完了。
Code
typedef long long ll;
typedef unsigned long long ull;
namespace solve
{
const int maxn = 1e7 + 10;
long long pri[maxn], phi[maxn], vis[maxn], mu[maxn], tot;
void init(const int n = 5e6)
{
mu[1] = phi[1] = 1;
for (int i = 2; i <= n; i++)
{
if (!vis[i])
pri[++tot] = i, mu[i] = -1, phi[i] = i - 1;
for (int j = 1; j <= tot && i * pri[j] <= n; j++)
{
vis[i * pri[j]] = 1;
if (i % pri[j] == 0)
{
mu[i * pri[j]] = 0;
phi[i * pri[j]] = phi[i] * pri[j];
break;
}
else
{
mu[i * pri[j]] = -mu[i];
phi[i * pri[j]] = phi[i] * (pri[j] - 1);
}
}
}
for (int i = 2; i <= n; i++)
mu[i] += mu[i - 1], phi[i] += phi[i - 1];
}
unordered_map<int, long long> mp;
int sum_mu(int n)
{
if (n <= 5e6)
return mu[n];
if (mp.count(n))
return mp[n];
long long res = 1;
for (long long i = 2, j; i <= n; i = j + 1)
{
j = n / (n / i);
res -= 1ll * (j - i + 1) * sum_mu(n / i);
}
return mp[n] = res;
}
long long sum_phi(int n)
{
if (n <= 5e6)
return phi[n];
long long res = 0;
for (long long i = 1, j; i <= n; i = j + 1)
{
j = n / (n / i);
res += (sum_mu(j) - sum_mu(i - 1)) * 1ll * (n / i) * (n / i);
}
return (res + 1) / 2;
}
void main()
{
init();
int T;
cin >> T;
while (T--)
{
int n;
cin >> n;
cout << sum_phi(n) << " " << sum_mu(n) << endl;
}
}
}
luogu3768 简单的数学题
Description
求： $$ \sum_{i=1}^n\sum_{j=1}^nij\gcd(i,j) $$ 对 $p$ 取模的结果。$n\leq10^{10}$。
Solution
不说啥了，直接开始推。 $$ \begin{aligned} ans =&\sum_{i=1}^n\sum_{j=1}^n\sum_{d=1,d|i,d|j}ijd[\gcd(i,j)=d]\\ =&\sum_{d=1}^nd\sum_{i=1,d|i}^ni\sum_{j=1,j|i}^nj[\gcd(i,j)=d]\\ =&\sum_{d=1}^nd^3\sum_{i=1}^{\lfloor n/d\rfloor}i\sum_{j=1}^{\lfloor n/d\rfloor}j[\gcd(i,j)=1]\\ =&\sum_{d=1}^nd^3\sum_{i=1}^{\lfloor n/d\rfloor}i\sum_{j=1}^{\lfloor n/d\rfloor}j\sum_{k|i,k|j}\mu(k)\\ =&\sum_{d=1}^nd^3\sum_{k=1}^{\lfloor n/d\rfloor}\mu(k)\sum_{i=1,k|i}^{\lfloor n/d\rfloor}i\sum_{j=1,k|j}^{\lfloor n/d\rfloor}j \end{aligned} $$ 设 $f(x)=\dfrac{x^2(x+1)^2}{4}$。 $$ \begin{aligned} ans=&\sum_{d=1}^nd^3\sum_{k=1}^{\lfloor n/d\rfloor}\mu(k)k^2f(\lfloor\dfrac{n}{dk}\rfloor) \end{aligned} $$ 设 $T=dk$。 $$ \begin{aligned} ans &=\sum_{d=1}^nd\sum_{T=1,d|T}^{\lfloor n/d\rfloor}\mu(\dfrac{T}{d})T^2f(\lfloor\dfrac{n}{T}\rfloor)\\ &=\sum_{T=1}^nf(\lfloor\dfrac{n}{T}\rfloor)T^2\sum_{d|T}d\times\mu(\dfrac{T}{d}) \end{aligned} $$ 由 $id*\mu=\varphi$，有： $$ ans=\sum_{T=1}^nf(\lfloor\dfrac{n}{T}\rfloor)T^2\varphi(T) $$ 设 $g(x)=x^2\varphi(x)$，我们需要求出它的前缀和，开始杜教筛。
设 $h(x)=x^2$，那么 $(h\ast g)(n)=n^2\sum_{d|n}\varphi(d)$，又由于 $\varphi\ast 1=id$，有 $(h\ast g)(x)=x^3$，这玩意的前缀和很好算。
然后就做完了。
Code
namespace solve
{
const int maxn = 1e7 + 10;
int mod;
int vis[maxn];
ll phi[maxn];
int pri[maxn], tot;
ll n;
void init(const int n = 1e7)
{
phi[1] = 1;
for (int i = 2; i <= n; i++)
{
if (!vis[i])
pri[++tot] = i, phi[i] = i - 1;
for (int j = 1; j <= tot && i * pri[j] <= n; j++)
{
vis[i * pri[j]] = 1;
if (i % pri[j])
phi[i * pri[j]] = phi[i] * phi[pri[j]];
else
{
phi[i * pri[j]] = phi[i] * pri[j];
break;
}
}
}
for (int i = 2; i <= n; i++)
phi[i] = 1ll * i * i % mod * phi[i] % mod, phi[i] += phi[i - 1], phi[i] %= mod;
}
ll qpow(ll a, ll x, ll p)
{
ll res = 1;
for (; x; x >>= 1, a = a * a % p)
if (x & 1)
res = res * a % p;
return res;
}
ll S(ll x) { return x * (x + 1) / 2 % mod; }
inline ll f(ll x) { return S(x % mod) * S(x % mod) % mod; }
ll inv6;
inline ll h(ll x)
{
x %= mod;
return 1ll * x * (x + 1) % mod * (2 * x + 1) % mod * inv6 % mod;
}
unordered_map<ll, ll> mp;
ll sum(ll n)
{
if (n <= 1e7)
return phi[n];
if (mp.count(n))
return mp[n];
ll res = f(n);
for (ll i = 2, j; i <= n; i = j + 1)
{
j = n / (n / i);
res -= (h(j) - h(i - 1) + mod) % mod * sum(n / i) % mod;
res %= mod;
}
return mp[n] = (res + mod) % mod;
}
void main()
{
cin >> mod >> n;
init();
inv6 = qpow(6, mod - 2, mod);
ll ans = 0;
for (ll i = 1, j; i <= n; i = j + 1)
{
j = n / (n / i);
ans += (sum(j) - sum(i - 1)) % mod * f(n / i) % mod;
ans %= mod;
}
ans += mod, ans %= mod;
cout << ans << endl;
}
}
NOI2016 循环之美
暂时不会推，先咕了。