问题
给定整数 (n,p(1le nle 10^{12}))，求 (sumlimits_{i=1}^{n}lfloor frac{n}{i} floormod p)。
思路
暴力做肯定会超时。我们发现加数中有许多是相同的，并且这些加数单调不增（即相同加数必定在一起）。如果找出每段相同加数的长度，就能很快得到答案。
也就是说，如果对于一个 (i)，能找出最大的 (j) 使得 (lfloor frac{n}{i}
floor=lfloorfrac{n}{j}
floor)
，这一段的和即为 (lfloorfrac{n}{i}
floor(j-i+1))。这样就能大大提高效率。
根据向下取整的性质，可得 (lfloor frac{n}{i} floorlefrac{n}{j}<lfloor frac{n}{i} floor+1)
将前一个不等式变形，得 (jle frac{n}{lfloor frac{n}{i} floor})
由于 (j) 是正整数，可得 (jle iglfloorfrac{n}{lfloor frac{n}{i} floor}ig floor)
所以 (j) 最大为 (iglfloorfrac{n}{lfloor frac{n}{i} floor}ig floor)。
代码
for(i = 1; i <= n; i = j + 1) {
j = n / (n / i);
ans += (j - i + 1) * (n / i) % p;
}
复杂度证明
-
当 (1le ile sqrt n) 时：(i) 的取值有 (sqrt n) 种，所以此时 (lfloorfrac{n}{i} floor) 的取值最多有 (sqrt n) 种。
-
当 (sqrt n< ile n) 时：(1lelfloorfrac{n}{i} floorlesqrt n)，所以 (lfloorfrac{n}{i} floor) 的取值最多有 (sqrt n) 种。
也就是说，(lfloorfrac{n}{i} floor) 的取值最多有 (2sqrt n) 种，所以复杂度为 (mathcal{O}(sqrt n))。
进阶
注：下列各式中 (mod p) 已省略。
- 求 (sumlimits_{i=1}^{min(n,m)}lfloor frac{n}{i} floorlfloorfrac{m}{i} floor)
每次右边界取 (min(iglfloorfrac{n}{lfloor frac{n}{i} floor}ig floor,iglfloorfrac{m}{lfloor frac{m}{i} floor}ig floor))，保证了两者分别相等，所以每两项的乘积相等。复杂度并没有变化。
- 求 (sumlimits_{i=1}^{n}ilfloorfrac{n}{i} floor)
对于每个块： (sumlimits_{i=l}^{r}ilfloor frac{n}{i} floor=sumlimits_{i=l}^{r}ilfloor frac{n}{l} floor=lfloor frac{n}{l} floorsumlimits_{i=l}^{r}i=lfloor frac{n}{l} floor imesfrac{(r-l+1)(l+r)}{2})
- 求 (sumlimits_{i=1}^{n}lceilfrac{n}{i} ceil)
(lceilfrac{n}{i} ceil=egin{cases}lfloorfrac{n}{i} floor&imid n\lfloorfrac{n}{i} floor+1&i mid nend{cases})
一共有 (d(n)) （表示 (n) 的因数个数）个数满足 (imid n)，我们计算全部都不能整除时的结果，最后减去因数个数。
( ext{原式}=sumlimits_{i=1}^{n}(lfloorfrac{n}{i} floor+1)-d(n)=n-d(n)+sumlimits_{i=1}^{n}lfloorfrac{n}{i} floor)
- 求 (sumlimits_{i=1}^{n}lfloorfrac{n}{i^2} floor)
令 (l'=l^2)，(r'=r^2)，由那个最基本的柿子可以得到 (r'=iglfloorfrac{n}{lfloor frac{n}{l'} floor}ig floor)。
所以 (r=Biglfloorsqrt{iglfloorfrac{n}{lfloor frac{n}{l^2} floor}ig floor}Big floor)
- 已知 (n,a,b) 为正整数，求 (sumlimits_{i=1}^{n}lfloorfrac{n}{ai+b} floor)
令 (l'=al+b)，(r'=ar+b)，由那个最基本的柿子可以得到 (r'=iglfloorfrac{n}{lfloor frac{n}{l'} floor}ig floor)。
所以 (ar+b=iglfloorfrac{n}{lfloor frac{n}{l'} floor}ig floor)
(r=Biglfloordfrac{iglfloorfrac{n}{lfloor frac{n}{l'} floor}ig floor-b}{a}Big floor)
一般地，对于求 (sumlimits_{i=1}^{n}lfloorfrac{k}{f(i)} floormod p( ext{其中}n,p,k ext{是正整数},1le nle10^{12})) 时，块边界推导方法是：
(lfloor frac{k}{f(l)} floor=lfloor frac{k}{f(r)} floor)，由那个最基本的柿子可知 (f(r)=iglfloorfrac{n}{lfloorfrac{n}{f(l)} floor}ig floor)。
求出 (f(l))，然后再解这个关于 (r) 的方程，最后向下取整，就得到了该块的右边界。
例题
模板题，直接上代码。
q = read();
for(; q; --q) {
ans = 0;
n = read();
for(i = 1; i <= n; i = j + 1) {
j = n / (n / i);
ans += (j - i + 1) * (n / i);
}
printf("%lld
", ans);
}
先来推柿子：(sumlimits_{i=1}^nkmod i=sumlimits_{i=1}^n(k-ilfloorfrac{k}{i} floor)=nk-sumlimits_{i=1}^nilfloorfrac{k}{i} floor)
对于每个块：(sumlimits_{i=l}^rilfloorfrac{k}{i} floor=sumlimits_{i=l}^rilfloorfrac{k}{l} floor=lfloorfrac{k}{l} floorsumlimits_{i=l}^ri=frac{lfloorfrac{k}{l} floor(l+r)(r-l+1)}{2})
代码：
scanf("%lld%lld", &n, &k);
for(i = 1; i <= n && k / i; i = j + 1) {
j = min(n, k / (k / i));
ans += (i + j) * (j - i + 1) * (k / i);
}
printf("%lld
", n * k - (ans >> 1));
题目中要求 (i eq j)，不太好做，我们可以先算出没有这一限制条件时的答案，再减去当 (i=j) 时的答案。
先令 (n) 为较小的那个。
( ext{原式}=sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}(nmod i)(mmod j)-sumlimits_{i=1}^{n}(nmod i)(mmod i))
(=[sumlimits_{i=1}^{n}(n-ilfloorfrac{n}{i} floor)][sumlimits_{j=1}^{m}(m-jlfloorfrac{m}{j} floor)]-sumlimits_{i=1}^{n}(n-ilfloorfrac{n}{i} floor)(m-ilfloorfrac{m}{i} floor))
(=(n^2-sumlimits_{i=1}^{n}ilfloorfrac{n}{i} floor)(m^2-sumlimits_{j=1}^{m}jlfloorfrac{m}{j} floor)-sumlimits_{i=1}^{n}(nm-milfloorfrac{n}{i} floor-nilfloorfrac{m}{i} floor+i^2lfloorfrac{n}{i} floorlfloorfrac{m}{i} floor))
(=(n^2-sumlimits_{i=1}^{n}ilfloorfrac{n}{i} floor)(m^2-sumlimits_{j=1}^{m}jlfloorfrac{m}{j} floor)-n^2m+msumlimits_{i=1}^{n}ilfloorfrac{n}{i} floor+nsumlimits_{i=1}^{n}ilfloorfrac{m}{i} floor-sumlimits_{i=1}^{n}i^2lfloorfrac{n}{i} floorlfloorfrac{m}{i} floor)
每个和式都可以整除分块。
注意题目中模数不是质数，可用 exgcd 求逆元。
代码：
#include <cstdio>
#include <cstring>
using namespace std;
#define in inline
typedef long long ll;
in ll max(ll x, ll y) {return x > y ? x : y;}
in ll min(ll x, ll y) {return x < y ? x : y;}
in void swap(ll &x, ll &y) {x ^= y ^= x ^= y;}
#define rei register int
#define FOR(i, l, r) for(rei i = l, i##end = r; i <= i##end; ++i)
#define FOL(i, r, l) for(rei i = r, i##end = l; i >= i##end; --i)
const int Mod = 19940417;
ll ans, n, m, inv2, inv6, tmp;
in ll MOD(ll x) {return (x % Mod + Mod) % Mod;}
void exgcd(ll a, ll p, ll &x, ll &y) {//exgcd求逆元
if(!p) {x = 1; y = 0; return;}
exgcd(p, a % p, x, y);
ll tmp = x;
x = y;
y = tmp - y * (a / p);
}
in ll Sum2(ll n) {return n * (n + 1) % Mod * (n + n + 1) % Mod * inv6 % Mod;}//求1^2+2^2+3^2+...+n^2
ll Sumx(ll n, ll k) {//1*(k/1)+2*(k/2)+...+n*(k/n)
ll res = 0;
for(register ll i = 1, j; i <= n && k / i; i = j + 1) {
j = min(k / (k / i), n);
res = (res + (i + j) * (j - i + 1) % Mod * (k / i) % Mod) % Mod;
}
return res * inv2 % Mod;
}
inline ll SumMod(ll n, ll k) {return (n * k % Mod - Sumx(n, k) % Mod + Mod) % Mod;}//(kmod1)+(kmod2)+...+(kmodn)
ll Sumxy(ll n, ll m) {//最后一个和式
ll res = 0;
for(register ll i = 1, j; i <= n; i = j + 1) {
j = min(n / (n / i) ,m / (m / i));
res = (res + (n / i) * (m / i) % Mod * MOD(Sum2(j) - Sum2(i - 1))) % Mod;
}
return res % Mod;
}
signed main() {
scanf("%lld%lld", &n, &m);
if(n > m) swap(n, m);
exgcd(2, 19940417, inv2, tmp); inv2 = (inv2 % Mod + Mod) % Mod;
exgcd(6, 19940417, inv6, tmp); inv6 = (inv6 % Mod + Mod) % Mod;
ans = SumMod(n, n) * SumMod(m, m) % Mod;
ans = MOD(ans - n * n % Mod * m % Mod);
ans = (ans + m * Sumx(n, n) % Mod + n * Sumx(n, m) % Mod) % Mod;
ans = MOD(ans - Sumxy(n, m));
printf("%lld
", ans);
return 0;
}