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 longusing 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);}