经典模型
[f(a,b,c,n)=sum_{i=0}^n lfloor frac{ai+b}{c}
floor
]
(a,b,c,n le 10^9,c ot= 0)，要求在 (log n) 的时间内求出。
约定
为了方便表示，做以下规定：
(x/y) 表示 (lfloor frac{a}{b} floor)
(a' = a mod c)
情况一
[a,b ge c
]
这个时候我们可以将 (a,b) 中大于 (c) 的部分拿出来，以简化问题：
[egin{aligned}
f(a,b,c,n ) &= sum_{i=0}^n lfloor frac{ai+b}{c}
floor\
&= frac{n imes (n + 1) imes a/c}{2} + (n+1) * (b/c)+ sum_{i = 0}^n lfloor frac{a'i+b'}{c}
floor\
&= frac{n imes (n + 1) imes a/c}{2} + (n+1) * (b/c)+ f(a',b',c,n)
end{aligned}
]
情况二
[a,b < c
]
这时候我们需要一些 trick 来想方设法地交换求和号：
[egin{aligned}
f(a,b,c,n) &= sum_{i=0}^n lfloor frac{ai+b}{c}
floor\
&= sum_{i=0}^n sum_{j=0}^{(ai+b)/c-1}1\
&= sum_{j=0}^{(an+b)/c-1}sum_{i=0}^n[j<(ai+b)/c]
end{aligned}
]
那个 (j < (ai+b)/c) 不是很好搞，但是如果关于底和顶的基本功比较好的话应该是问题不大的。
先抛俩较为基础的结论:（(n) 为整数，(x) 为实数）
[n le x Leftrightarrow n le lfloor x
floor
]
[n > x Leftrightarrow n > lfloor x
floor
]
然后推：
[j < (ai+b)/c\
j+1 le (ai+b)/c\
j+1 le frac{ai+b}{c}\
cj+c le ai+b\
cj+c-b le ai\
cj+c-b-1 < ai\
i > (cj+c-b-1)/a
]
然后原式子就好看很多了。定义 (m = (an+b)/c)，然后：
[egin{aligned}
f(a,b,c,n) &= sum_{j=0}^{(an+b)/c-1}sum_{i=0}^n[j<(ai+b)/c]\
&= sum_{j=0}^{m-1}n-lfloor frac{cj+c-b-1}{a}
floor\
&= nm-f(c,c-b-1,a,n)
end{aligned}
]
感觉没啥用？其实是有用的，观察到 (a,c) 换了位置，而情况一可以将 (a) 模 (c)，因此复杂度和辗转相除法的复杂度一样，都是 (log n)。
情况三(边界)
[a=0
]
[f(a,b,c,n)=(n+1)lfloor frac{b}{c}
floor
]
拓展
现在又跑出来俩类似 (f) 的函数 (g) 和 (h)，(h) 是二次的 (f)，(g) 是 (h) 的辅助数组：
[g(a,b,c,n) = sum_{i = 0}^n i lfloor frac{ai+b}{c}
floor\
h(a,b,c,n) = sum_{i=0}^n lfloor frac{ai+b}c
floor ^2
]
这俩式子的推导与 (f) 的思想类似，都是分为三种情况，第一种取模，第二种先办法搞出求和号然后交换以实现 (a,c) 互换位置的效果，第三种边界。第三种太水了，就不写了。
g函数
情况一
[egin{aligned}
g(a,b,c,n) &= sum_{i=0}^n i lfloor frac{ai+b}c
floor\
&= sum_{i=0}^n i(i imes a/c + b/c) lfloor frac{a'i+b'}c
floor \
&= a/c imes n(n+1)(2n+1)/6 + b/c imes n(n+1)/2 + g(a',b',c,n)
end{aligned}
]
情况二
[egin{aligned}
g(a,b,c,n) &= sum_{i=0}^n i sum_{j=0}^{(ai+b)/c-1}1\
&= sum_{j=0}^{m-1}sum_{i=0}^ni[j < (ai+b)/c]\
&= sum_{j=0}^{m-1}n(n+1)/2 - (jc+c-b-1)/a imes ((jc+c-b-1)/a-1)/2\
&= frac{1}{2}(mn(n+1)-h(c,c-b-1,a,m-1)-f(c,c-b-1,a,m-1))
end{aligned}
]
h函数
情况一
公式太难打了，就直接给结论吧。反正推导全是暴力拆式子
[h(a,b,c,n) = \
small (a/c)^2n(n+1)(2n+1)/6 + (a/c)(b/c)n(n+1)+(n+1)(b/c)^2+2(a/c)g(a',b',c,n)+2(b/c)f(a',b',c,n)+h(a',b',c,n)
]
情况二
这里还需要一个小 trick：（普通幂表示成自然数幂和）
[n^2 = (2sum_{i=0}^{n-1}i) + n
]
然后就可以推了：
[egin{aligned}
h(a,b,c,n) &= sum_{i=0}^n 2 sum_{j=0}^{(ai+b)/c-1}j + f(a,b,c,n)\
&= 2 sum_{j=0}^{m-1}j sum_{i=0}^n[(ai+ b)/c > j] + f(a,b,c,n)\
&= 2sum_{j=0}^{m-1}j(n-lfloor frac{cj+c-b-1}a
floor)+f(a,b,c,n)\
&= nm(m-1) - 2g(c,c-b-1,a,m-1)+f(a,b,c,n)
end{aligned}
]
由于 (f)， (g) 和 (h) 需要相互调用，我们可以每次打包返回三个值。
题目与模板
类欧模板（调试用）
struct node {
ll f, g, h;
node() {}
node(ll a, ll b, ll c) { f = a, g = b, h = c; }
};
node sol(int a, int b, int c, int n) {
ll nwf = 0, nwg = 0, nwh = 0;
if (!a) {
nwf = (b / c) * (n + 1) % P;
nwg = (b / c) * n % P * (n + 1) % P * inv2 % P;
nwh = (b / c) * (b / c) % P * (n + 1) % P;
return node(nwf, nwg, nwh);
}
if (a >= c || b >= c) {
node nd = sol(a % c, b % c, c, n);
ll tf = nd.f, tg = nd.g, th = nd.h;
nwf = ((a / c) * n % P * (n + 1) % P * inv2 + (b / c) * (n + 1) + tf) % P;//bug
nwg = ((a / c) * n % P * (n + 1) % P * (2ll * n + 1) % P * inv6 + (b / c) * n % P * (n + 1) % P * inv2 + tg) % P;
nwh = (a / c) * (a / c) % P * n % P * (n + 1) % P * (2ll * n + 1) % P * inv6;
nwh = (nwh + (a / c) * (b / c) % P * n % P * (n + 1) + (n + 1) * (b / c) % P * (b / c)) % P;
nwh = (nwh + 2ll * (b / c) * tf + 2ll * (a / c) * tg + th) % P;
return node(nwf, nwg, nwh);
}
ll m = (a * n + b) / c;//bug
node nd = sol(c, c - b - 1, a, m - 1);
m %= P;
ll tf = nd.f, tg = nd.g, th = nd.h;
nwf = (m * n - tf) % P;
nwg = (m * n % P * (n + 1) - th - tf) % P * inv2 % P;
nwh = (n * m % P * (m - 1) - 2ll * tg + nwf) % P;
return node(nwf, nwg, nwh);
}
注意一下取模的问题，计算答案要取模，但是向下递归的 (m) 不能扔个 (m mod P) 进去。