注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

ydc的博客

 
 
 

日志

 
 

JZPKIL题解  

2013-09-12 22:14:37|  分类: bzoj |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
多组数据
输入n,x,y
输出
JZPKIL题解(未完待续) - ydc - ydc的博客
 

以下是若逼的题解时间:


JZPKIL题解(未完待续) - ydc - ydc的博客
 我们要解决的就是如何计算
JZPKIL题解(未完待续) - ydc - ydc的博客

 我们有
JZPKIL题解(未完待续) - ydc - ydc的博客
 于是我们继续推导
JZPKIL题解(未完待续) - ydc - ydc的博客
 
于是我们就需要求JZPKIL题解(未完待续) - ydc - ydc的博客 了:
JZPKIL题解 - ydc - ydc的博客
 
 
这是在省选集训时hzhwcmhf大神教我的:伯努利数,Bi就是传说中的伯努利数
大家可以去看维基百科
于是我们得到一个(m+1)次的通项公式。
JZPKIL题解(未完待续) - ydc - ydc的博客
 其中
JZPKIL题解(未完待续) - ydc - ydc的博客
 我们就可以继续推倒了
JZPKIL题解(未完待续) - ydc - ydc的博客
 
这个式子不是积性函数,我们把这么来看:
JZPKIL题解(未完待续) - ydc - ydc的博客
 
 无视掉t[k]后,就变成了若干个这种形式的多项式相加:
JZPKIL题解(未完待续) - ydc - ydc的博客
 
 他是三个函数的积性函数卷积形式,所以也是积性函数。既然如此,我们可以把n进行分解质因数,对每个p1^k的因子单独考虑,最后利用积性函数的性质可以将n代入这个多项式的结果,再将他*t[k]得到的多项式累加,就得到了答案了
分解质因数我们用Pollard_rho算法,可以在n^(0.25)时间里诡异的搞出来
由于这个时候的n=p^k的形式,JZPKIL题解(未完待续) - ydc - ydc的博客若不为0,d1要么为1,要么为p,然后就可以做了
于是算法流程是:
预处理sigma{i^k}的系数;对每组询问,将n分解质因数;对每个n=p^k,讨论d1=1以及d1=p,然后代入公式计算
关于a*b%c的问题,gyz有种神奇的O(1)的代码,但无数事实证明,那种算法是存在着精度误差的,那么这个精度误差会不会导致你的程序出bug,就看你写的如何了。vfk神啊,kzf神啊,都可以靠那东西过掉;但也的确存在用了那东西就错了的人(比如我),这根Pollard_rho的写法之类的东西是密切相关的(所以我考试打算还是用快速乘算了)
于是时间复杂度是O(T*(Pollard_rho的复杂度)+m^2)

本机我的程序开O2下是无压力的……但是bzoj不知道为什么,明明应该能跑过去,他偏偏就是没让我跑过去(不过bzoj只有50s时限,于是均摊变成了2.5s一个点)
最后无耻的加了个特判:如果x==y,就直接利用伯努利数计算答案,就成功降低了总耗时……卡过去了

update:有关于快速乘的精度误差已经解决
那个O(1)快速乘的精度误差在于,我要计算a*b/c取下整的时候,使用的是(LL)((long double)a*b/c)
这自然有可能有精度误差……所以gyz的处理方法是(LL)((long double)a*b/c+1e-3)
但是由于a,b,c可能很大,a*b/c=0.9999999991的之类的可能性是存在的……所以根据个人写法不同,那东西还是有可能会挂的
我这里有一个精度问题更小的写法,个人感觉几乎不会出错:
inline LL mul(LL x,LL y,LL z)
{
    return (x*y-(LL)(((long double)x*y+0.5)/(long double)z)*z+z)%z;
}
正确性与精度的优越性显然
改过的代码就不贴了
#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<ctime>
#define mod 1000000007
#define maxn 3010
#define maxm 70
#define maxp 1010
using namespace std;
typedef long long LL;
typedef unsigned long long ULL;
LL B[maxn];
int t[maxn][maxn],comb[maxn][maxn],prime[maxp],top;
pair<LL,LL> ex_gcd(LL a,LL b)
{
if(b==0)
return make_pair(1,0);
pair<LL,LL> temp=ex_gcd(b,a%b);
return make_pair(temp.second,temp.first-a/b*temp.second);
}
inline LL Anti(LL n)
{
LL ans=ex_gcd(n,mod).first;
if(ans<0)
ans+=mod;
return ans;
}
LL Fact[maxn],iFact[maxn];
void Prepare()
{
int n=3000;
comb[0][0]=1;
Fact[0]=1,iFact[0]=1;
for(int i=1;i<=n+1;++i)
Fact[i]=Fact[i-1]*i%mod,iFact[i]=Anti(Fact[i]);
for(int i=1;i<=n+1;++i)
{
comb[i][0]=1;
for(int j=1;j<=i;++j)
{
comb[i][j]=comb[i-1][j]+comb[i-1][j-1];
if(comb[i][j]>=mod)
comb[i][j]-=mod;
}
}
for(int i=0;i<=n;++i)
{
B[i]=i+1;
for(int j=0;j<i;++j)
{
B[i]=(B[i]-comb[i+1][j]*B[j])%mod;
if(B[i]<0)
B[i]+=mod;
}
B[i]=B[i]*Anti(comb[i+1][i])%mod;
}
for(int i=0;i<=n;++i)
{
LL a=Anti(i+1);
for(int j=0;j<=i;++j)
t[i][j]=B[j]*comb[i+1][j]%mod*a%mod;
}
static bool use[maxp];
int m=1000;
for(int i=2;i<=m;++i)
for(int j=2;j*i<=m;++j)
use[i*j]=true;
for(int i=2;i<=m;++i)
if(use[i]==false)
prime[++top]=i;
}
inline LL gcd(LL a,LL b)
{
static LL t;
while(b)
t=a,a=b,b=t%b;
return a;
}
inline LL mul(LL p,LL n,LL mo)
{
LL ans=0;
while(n)
{
if(n&1)
{
ans+=p;
if(ans>=mo)
ans-=mo;
}
n>>=1,p<<=1;
if(p>=mo)
p-=mo;
}
return ans;
}
inline LL power(LL p,LL n,LL mo)
{
LL ans=1;
for(;n;n>>=1,p=mul(p,p,mo))
if(n&1)
ans=mul(ans,p,mo);
return ans;
}
inline LL power(LL p,LL n)
{
if(p>=mod)
p%=mod;
LL ans=1;
while(n)
{
if(n&1)
{
ans=ans*p;
if(ans>=mod)
ans%=mod;
}
n>>=1,p*=p;
if(p>=mod)
p%=mod;
}
return ans;
}
inline bool Miller_Rabin(int p,LL n)
{
int s=0;
LL d=n-1;
while(!(d&1))
++s,d>>=1;
LL w=power(p,d,n);
if(w==1)
return true;
for(int i=0;i<=s-1;++i)
{
if(w==n-1)
return true;
w=mul(w,w,n);
}
return false;
}
inline bool isPrime(LL n)
{
static int prime[]={2,3,5,7,11,13,17,19,23};
if(n==1)
return false;
if(find(prime,prime+9,n)!=prime+9)
return true;
if(!(n&1))
return false;
for(int i=0;i<9;++i)
if(!Miller_Rabin(prime[i],n))
return false;
return true;
}
int x,y;
LL n;
LL d[maxm];
LL p[maxm],kp[maxm],pw[maxm][maxm];
int s[maxm],nd,sp;
void Pollard_rho(LL n)
{
if(isPrime(n))
{
d[++nd]=n;
return ;
}
int c=1;
while(1)
{
LL x1=1,x2=1;
int i=1,k=2;
while(1)
{
x1=mul(x1,x1,n)+c;
LL d=gcd(abs(x1-x2),n);
if(d!=1&&d!=n)
{
Pollard_rho(d),Pollard_rho(n/d);
return ;
}
if(x1==x2)
break;
if(++i==k)
k<<=1,x2=x1;
}
++c;
}
}
LL Kill(LL n,LL x,LL y)
{
LL ans=0;
nd=0;
for(int i=1;i<=top;++i)
while(n%prime[i]==0)
n/=prime[i],d[++nd]=prime[i];
if(n!=1)
Pollard_rho(n);
sort(d+1,d+nd+1);
sp=1,p[1]=d[1],s[1]=1,kp[1]=d[1];
for(int i=2;i<=nd;++i)
{
if(d[i]!=d[i-1])
++sp,p[sp]=d[i],s[sp]=1,kp[sp]=d[i];
else
++s[sp],kp[sp]*=d[i];
}
for(int i=1;i<=sp;++i)
{
pw[i][0]=1;
LL a=power(p[i],x);
for(int j=1;j<=s[i];++j)
{
pw[i][j]=pw[i][j-1]*a;
if(pw[i][j]>=mod)
pw[i][j]%=mod;
}
}
static LL pw1[maxm];
for(int i=0;i<=y;++i)
{
LL sum=1;
for(int j=1;j<=sp;++j)
{
LL S1=0,S2=0,a=power(p[j],y),b=power(p[j],y+1-i);
pw1[0]=1;
for(int k=1;k<=s[j];++k)
{
pw1[k]=pw1[k-1]*b;
if(pw1[k]>=mod)
pw1[k]%=mod;
}
for(int k=0;k<=s[j];++k)
{
S1=S1+pw[j][k]*pw1[s[j]-k];
if(S1>=mod)
S1%=mod;
}
for(int k=0;k<s[j];++k)
{
S2=S2+pw[j][k]*pw1[s[j]-k-1]%mod*a;
if(S2>=mod)
S2%=mod;
}
S1=(S1-S2)%mod;
if(S1<0)
S1+=mod;
sum=sum*power(kp[j],y);
if(sum>mod)
sum%=mod;
sum*=S1;
if(sum>mod)
sum%=mod;
}
ans=ans+sum*t[y][i];
if(ans>=mod)
ans%=mod;
}
return ans%mod;
}
LL calc(LL n)
{
LL ans=0,p=n%mod,d=p;
for(int i=y;i>=0;--i)
{
ans=ans+t[y][i]*d;
if(ans>=mod)
ans%=mod;
d=d*p;
if(d>=mod)
d%=mod;
}
return ans*power(n,y)%mod;
}
int main()
{
Prepare();
int test;
for(cin>>test;test;--test)
{
cin>>n>>x>>y;
if(x==y)
cout<<calc(n)<<endl;
else
cout<<Kill(n,x,y)<<endl;
}
return 0;
}



感谢cmg大神和yc5_yc大神教会我使用极其基础的Latex (为什么我用Latex打不出中文求解释)
  评论这张
 
阅读(1655)| 评论(4)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017