bzoj 4815 小Q的表格 —— 反演+分块

题目:https://www.lydsy.com/JudgeOnline/problem.php?id=4815

思路就和这里一样:https://blog.csdn.net/leolyun/article/details/70146612

不知为何乘逆元就错了,必须直接除...不过题目保证了是整数,所以直接除也没问题;

然后重新学习了一下分块的简洁写法,就能A了hhh

代码如下:

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
int const xn=4e6+5,mod=1e9+7,base=2000,xb=2005;
int n,cnt,pri[xn],phi[xn],g[xn],blk[xn],s[xn],pls[xb];
bool vis[xn];
ll rd()
{
  ll ret=0,f=1; char ch=getchar();
  while(ch<'0'||ch>'9'){if(ch=='-')f=0; ch=getchar();}
  while(ch>='0'&&ch<='9')ret=ret*10+ch-'0',ch=getchar();
  return f?ret:-ret;
}
ll pw(ll a,int b)
{
  ll ret=1;
  for(;b;b>>=1,a=(a*a)%mod)if(b&1)ret=(ret*a)%mod;
  return ret;
}
int upt(int x){while(x>=mod)x-=mod; while(x<0)x+=mod; return x;}
void init()
{
  phi[1]=1;
  for(int i=2;i<=n;i++)
    {
      if(!vis[i])pri[++cnt]=i,phi[i]=i-1;
      for(int j=1;j<=cnt&&(ll)i*pri[j]<=n;j++)
    {
      vis[i*pri[j]]=1;
      if(i%pri[j])phi[i*pri[j]]=(ll)phi[i]*(pri[j]-1)%mod;
      else {phi[i*pri[j]]=(ll)phi[i]*pri[j]%mod; break;}
    }
    }
  for(int i=1;i<=n;i++)g[i]=(ll)i*i%mod*phi[i]%mod;
  for(int i=2;i<=n;i++)g[i]=upt(g[i-1]+g[i]);
}
void init2()
{
  for(int i=1;i<=n;i++)blk[i]=(i-1)/base+1;
  for(int i=1;i<=n;i++)s[i]=(s[i-1]+(ll)i*i)%mod;
}
void change(int d,ll x)//add x
{
  for(int i=d;blk[i]==blk[d]&&i<=n;i++)s[i]=upt(s[i]+x);
  for(int i=blk[d]+1;i<=blk[n];i++)pls[i]=upt(pls[i]+x);
}
int cal(int ps)
{
  return upt(pls[blk[ps]]+s[ps]);
}
int gcd(int a,int b){return b?gcd(b,a%b):a;}
int main()
{
  int m=rd(); n=rd();
  init(); init2(); ll x;
  for(int i=1,a,b,k;i<=m;i++)
    {
      a=rd(); b=rd(); x=rd(); k=rd();
      int d=gcd(a,b); 
      //int t1=(ll)d*pw(a,mod-2)%mod,t2=(ll)d*pw(b,mod-2)%mod;
      ll t=((ll)a/d)*(b/d); x=(x/t)%mod;//
      change(d,upt(x-cal(d)+cal(d-1)));
      //change(d,x*t1%mod*t2%mod);
      int ans=0;
      for(int j=1,nxt;j<=k;j=nxt+1)
    {
      nxt=k/(k/j);
      ans=(ans+(ll)g[k/j]*upt(cal(nxt)-cal(j-1)))%mod;
    }
      printf("%d\n",ans);
    }
  return 0;
}

 

posted @ 2018-12-13 22:15  Zinn  阅读(134)  评论(0编辑  收藏  举报