NOIP2016提高组 快速荷叶叶变换

Description

Data Constraint

对于 40% 的数据，1 ≤ N，M ≤ 1000

对于 60% 的数据，1 ≤ N，M ≤ 10^6

对于 100% 的数据，1 ≤ N，M ≤ 10^9

Solution

我们将N和M分开讨论。对于一个固定的i，∑Mj=1（M mod j）我们发现当⌊M/j⌋固定为一个值k时，（M mod j）会随着j的增大而减小，而且每次减小的值都为k，直到⌊M/j⌋改变。这样的话，我们就可以采用分块。每次将⌊M/j⌋相同的所有余数统计起来，这样就可以统计出所有的∑Mj=1（M mod j），最后在用相同的方法统计一下N的，两者相乘即可。因为一个数的约数不会超过2*N−−√，所以整个时间复杂度为O(N−−√)。

代码

#include<iostream>
#include<cmath>
#include<cstring>
#include<cstdio>
#include<algorithm>
#define ll long long
using namespace std;
const ll maxn=1000000007;
ll n,m,i,t,j,k,l,ans,x,y,p;
ll mi(int y){
if (y==1) return 2;
ll t=mi(y/2);
if (y%2) return t*t%maxn*2%maxn;return t*t%maxn;
}
int main(){
//  freopen("data.in","r",stdin);
scanf("%d%d",&n,&m);i=1;
while (i<m){
k=m%i;
j=m/(m/i);
l=m%j;
t=m/i;
p=(k-l)/t+1;
x=(x+(k+l)%maxn*p%maxn*mi(maxn-2)%maxn)%maxn;
i=j+1;
}
i=1;
while (i<n){
k=n%i;
j=n/(n/i);
l=n%j;
t=n/i;
p=(k-l)/t+1;
y=(y+(k+l)%maxn*p%maxn*mi(maxn-2)%maxn)%maxn;
i=j+1;
}
ans=x*y%maxn;
printf("%lld\n",ans);
}