「学习笔记」类欧几里得算法
现有多组询问，要求在 (O(log n)) 求三个有趣式子的和。
设
[f(a,b,c,n)=sum_{i=0}^{n}lfloor frac {ai+b}{c}
floor
]
[g(a,b,c,n)=sum_{i=0}^{n}lfloor frac {ai+b}{c}
floor^2
]
[h(a,b,c,n)=sum_{i=0}^{n}ilfloor frac {ai+b}{c}
floor
]
(t_1=lfloor frac ac floor,t_2=lfloor frac bc floor)
(S_1(n)=sum_{i=0}^{n}i,S_2(n)=sum_{i=0}^{n}i^2)
(m=lfloor frac {an+b}{c} floor)
求 (f(a,b,c,n))
若 (ageq c) 或 (bgeq c)
(lfloor frac {ai+b}{c} floor=lfloor frac {(a ext{mod} c)i+(b ext{mod} c)}{c} floor+t_1i+t_2)
(Longrightarrow f(a,b,c,n)=f(a ext{mod} c,b ext{mod} c,c,n)+t_1S_1(n)+(n+1)t_2)
若 (a<c) 且 (b<c)
[f(a,b,c,n)=sum_{i=0}^{n}sum_{j=1}^{m}[lfloor frac {ai+b}{c}
floorgeq j]=sum_{i=0}^{n}sum_{j=0}^{m-1}[lfloor frac {ai+b}{c}
floorgeq j+1]
]
[=sum_{i=0}^{n}sum_{j=0}^{m-1}[ai+bgeq cj+c]=sum_{i=0}^{n}sum_{j=0}^{m-1}[aigeq cj+c-b]=sum_{i=0}^{n}sum_{j=0}^{m-1}[ai> cj+c-b-1]
]
[sum_{i=0}^{n}sum_{j=0}^{m-1}[i>lfloor frac {cj+c-b-1}{a}
floor]=sum_{j=0}^{m-1}sum_{i=0}^{n}[i>lfloor frac {cj+c-b-1}{a}
floor]=sum_{j=0}^{m-1}n-lfloor frac {cj+c-b-1}{a}
floor
]
[=mn-f(c,c-b-1,a,m-1)
]
求 (g(a,b,c,n)) 和 (h(a,b,c,n))
由于式子过长且与 (f(a,b,c,n)) 的推导过程类似，所以只有简化过程。
若 (ageq c) 或 (bgeq c)
(g(a,b,c,n)=g(a ext{mod} c,b ext{mod} c,c,n)+2t_1h(a ext{mod} c,b ext{mod} c,c,n)+2t_2f(a ext{mod} c,b ext{mod} c,c,n))
(+t_1^2S_2(n)+2t_1t_2S_1(n)+(n+1)t_2^2)
(h(a,b,c,n)=h(a ext{mod} c,b ext{mod} c,c,n)+t_1S_2(n)+t_2S_1(n))
若 (a<c) 且 (b<c)
[g(a,b,c,n)=sum_{j=0}^{m-1}sum_{k=0}^{m-1}n-max(lfloor frac {cj+c-b-1}{a}
floor,lfloor frac {ck+c-b-1}{a}
floor)
]
[=m^2n-2sum_{j=0}^{m-1}(j+1)lfloor frac {cj+c-b-1}{a}
floor+sum_{j=0}^{m-1}lfloor frac {cj+c-b-1}{a}
floor
]
[=m^2n-2h(c,c-b-1,a,m-1)-f(c,c-b-1,a,m-1)
]
[h(a,b,c,n)=sum_{j=0}^{m-1}sum_{i=0}^{n}i[i>lfloor frac {cj+c-b-1}{a}
floor]=sum_{j=0}^{m-1}S_1(n)-S_1(lfloor frac {cj+c-b-1}{a}
floor)
]
[=mS_1(n)-frac {g(c,c-b-1,a,m-1)}{2}-frac {f(c,c-b-1,a,m-1)}{2}
]
边界条件
若 (n=0)
[f(a,b,c,n)=t_2,g(a,b,c,n)=t_2^2,h(a,b,c,n)=0
]
若 (a=0)
[f(a,b,c,n)=(n+1)t_2,g(a,b,c,n)=(n+1)t_2^2,h(a,b,c,n)=t_2S_1(n)
]
(Code Below:)
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const ll mod=998244353;
const ll inv2=499122177;
const ll inv6=166374059;
ll n,a,b,c;
struct node{
ll f,g,h;
};
inline ll S1(ll n){
return n*(n+1)%mod*inv2%mod;
}
inline ll S2(ll n){
return n*(n+1)%mod*(2*n+1)%mod*inv6%mod;
}
inline node solve(ll a,ll b,ll c,ll n){
ll t1=a/c,t2=b/c,s1=S1(n),s2=S2(n),m=(a*n+b)/c;
node ans,now;ans.f=ans.g=ans.h=0;
if(!n){
ans.f=t2;
ans.g=t2*t2%mod;
return ans;
}
if(!a){
ans.f=(n+1)*t2%mod;
ans.g=(n+1)*t2%mod*t2%mod;
ans.h=t2*s1%mod;
return ans;
}
if(a>=c||b>=c){
now=solve(a%c,b%c,c,n);
ans.f=(now.f+t1*s1+(n+1)*t2)%mod;
ans.g=(now.g+2*t1*now.h+2*t2*now.f+t1*t1%mod*s2+2*t1*t2%mod*s1+(n+1)*t2%mod*t2)%mod;
ans.h=(now.h+t1*s2+t2*s1)%mod;
return ans;
}
now=solve(c,c-b-1,a,m-1);
ans.f=(m*n-now.f)%mod;ans.f=(ans.f+mod)%mod;
ans.g=(m*m%mod*n-2*now.h-now.f);ans.g=(ans.g+mod)%mod;
ans.h=(m*s1-now.g*inv2-now.f*inv2)%mod;ans.h=(ans.h+mod)%mod;
return ans;
}
int main()
{
ll T;
scanf("%lld",&T);
while(T--){
scanf("%lld%lld%lld%lld",&n,&a,&b,&c);
node ans=solve(a,b,c,n);
printf("%lld %lld %lld
",ans.f,ans.g,ans.h);
}
return 0;
}