BZOJ2627 JZPKIL
一句話題意：求$\sum_{i=1}^ngcd(i,n)^xlcm(i,n)^y$，多組數據，模$1e9+7$。php
30%的數據，$x=y$
另30%的數據，$n<=10^9$, $x, y<=100$
100%的數據，$T<=100, 1<=n<=10^{18}, 0<=x, y<=3000$html
orz出題人顧宇宙！！！orz全網題解的源頭ydc！！！感受ydc的題解已經寫得很清楚了，可是爲了內部交流仍是寫一寫……下面的公式中$[x]$表示$x==1$。ios
推導部分算法
Part1 原式整理$$\sum_{i=1}^ngcd(i,n)^xlcm(i,n)^y$$$$=\sum_{i=1}^n(in)^ygcd(i,n)^{x-y}$$$$=n^y\sum_{i=1}^ni^ygcd(i,n)^{x-y}$$$$=n^y\sum_{d|n}d^{x-y}\sum_{k=1,[gcd(k,\frac{n}{d}))]}^{\frac{n}{d}}(kd)^y$$$$=n^y\sum_{d|n}d^x\sum_{k=1,[gcd(k,\frac{n}{d})]}^{\frac{n}{d}}k^y$$ide
Part2 套路反演函數
這一部分咱們來整理上式最後的$\sum_{k=1,[gcd(k,\frac{n}{d})]}^{\frac{n}{d}}k^y$spa
$$\sum_{i=1,[gcd(i,n)]}^ni^y$$$$=\sum_{i=1}^ni^y\sum_{d|gcd(i,n)}\mu(d)$$$$=\sum_{d=1}^n\mu(d)\sum_{i=1}^{\frac{n}{d}}(id)^y$$$$=\sum_{d=1}^n\mu(d)d^y\sum_{i=1}^{\frac{n}{d}}i^y$$code
Part3 伯努利數xml
能夠發現Part2的最後咱們把Part1的冪之和化爲了從1開始的連續一段……這有什麼用呢？咱們能夠套結論了！有一個叫伯努利數的東西（下文推導中用$B_k$表示第$k$個伯努利數），在網上搜一搜的話資料很多，百度百科裏基本的公式也給了一些。剛開始看證實說什麼從生成函數定義的角度能夠證實它和天然數冪和的數學關係，可是我後來看到它就是在求天然數冪和的過程當中被髮明的就感受證這個不必了……求伯努利數的方法放在實現部分裏，如今咱們繼續關注推導過程。htm
$$\sum_{i=1}^ni^m=\frac{1}{m+1}\sum_{k=0}^mC_{m+1,k}B_kn^{m-k+1}$$
大多數狀況下伯努利數都是這麼用的……能夠當作公式
你可能會發現伯努利數的公式很亂，有人用$n$的冪有人用$n+1$的冪還有人帶-1的冪……然而由於不會證我也不能確切地說哪種是對的。維基百科上大概是說伯努利數有兩類，兩類伯努利數的區別在於B1。第一類伯努利數，由美國國家標準技術研究所 (NIST)制定，在這標準下B1=-1/2;第二類伯努利數，又被稱爲是「原始的伯努利數」，在這標準下B1=1/2。$n+1$的式子是關於第一類伯努利數的，而$n$的式子是關於第二類伯努利數的;第一類伯努利數乘上-1的冪即爲第二類伯努利數……僞裝他們都是對的好了編套理論真難啊
設$t_k=\frac{1}{m+1}C_{m+1,k}B_k$，則$\sum_{i=1}^ni^m=\sum_{k=0}^mt_kn^{m-k+1}$
Part4 積性函數
即便化成這步代入原式它也很差處理……不過看一眼$10^{18}$的數據範圍，分解質因數來作好像是一種歷史必然？so如今咱們試一試把積性函數的部分整理到一塊兒會發生什麼……
$$n^y\sum_{d|n}d^x\sum_{d_1|\frac{n}{d}}\mu(d_1)d_1^y\sum_{k=0}^yt_k(\frac{n}{dd_1})^{y-k+1}$$
$$=\sum_{k=0}^yt_kn^y\sum_{d|n}d^x\sum_{d_1|\frac{n}{d}}\mu(d_1)d_1^y(\frac{n}{dd_1})^{y-k+1}$$
這時候除去$t_k$的部分就是一個積性函數了，$y$很小隻有$3000$，咱們能夠枚舉$k$而後分解質因數求後面的部分。那麼就變成了對每個質因子$p_i$的$k_i$冪$p_k$求$p_k^y\sum_{d|p_k}d^x\sum_{d_1|\frac{p_k}{d}}\mu(d_1)d_1^y(\frac{p_k}{dd_1})^{y-k+1}$，又由於$\mu(d_1)$只有在$d_1=1$時爲一、$d_1=p_i$時爲-1而其餘狀況下都爲0，式子能夠直接變成$p_k^y\sum_{d|p_k}d^x((\frac{p_k}{d})^{y-k+1}-p_i^y\frac{p_k}{dp_i}^{y-k+1})$，這就是推式子的終點了。
實現部分
Part1 伯努利數
你問我伯努利數怎麼求啊……
$\sum_{k=0}^nC_{n+1,k}B_k=n+1$
因此楊輝三角求一下組合數就能夠$O(k^2)$求出全部伯努利數了。
Part2 分解質因數
給$10^{18}$分解質因數？須要$n^{\frac{1}{4}}$的算法！Miller-Rabin和Pollard-Rho能夠解決，關於這兩個算法咱們的朋友Doggu的blog堪稱全網最良心。
Part3 快速乘
$10^{18}*10^{18}$應該怎麼作？快速乘顯然是必要的。這是ydc改良過的快速乘板子。
1 inline LL mul(LL x,LL y,LL z) 2 { 3 return (x*y-(LL)(((long double)x*y+0.5)/(long double)z)*z+z)%z; 4 }
好像少了點什麼東西
沒有錯，這題是有部分分的吧？對於$x=y$的$30$分，式子就是$\sum_{i=1}^n(in)^y$，若是你會用伯努利數就能夠拿到這一部分的分。對於$n<=10^9$的$30$分，伯努利數仍然是必要的，不過度解質因數能夠$\sqrt{n}$來作。看來出題人只打算給會伯努利數的人分數（顯然我不在此之列，因此部分分和正解同樣困難……）。想打部分分的同窗能夠去cogs上提交。
最後放代碼，完結撒花。
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 #include<cstdlib> 5 #include<algorithm> 6 #define ll long long 7 using namespace std; 8 const int mod=1e9+7; 9 const int sj=3010; 10 const int ss=110; 11 inline ll read() 12 { 13 ll jg=0;int jk=getchar()-'0'; 14 while(jk<0||jk>9) jk=getchar()-'0'; 15 while(jk>=0&&jk<=9) jg*=10,jg+=jk,jk=getchar()-'0'; 16 return jg; 17 } 18 int ca,ai,bi,cnt,tot,s[ss],sp,t[sj][sj],c[sj][sj]; 19 ll n,a1,a2,b[sj],pi[ss],pk[ss]; 20 ll p[sj],prime[10]={2,3,5,7,11,13,17,19,23,29},fac[ss]; 21 ll pw[ss][ss],pw1[ss]; 22 bool vi[sj]; 23 ll gcd(ll x,ll y) 24 { 25 if(!y) return x; 26 return gcd(y,x%y); 27 } 28 inline ll mul(ll x,ll y,ll z) 29 { 30 return (x*y-(ll)(((long double)x*y+0.5)/(long double)z)*z+z)%z; 31 } 32 ll qpow(ll x,ll y,ll z) 33 { 34 ll ret=1;x%=z; 35 while(y) 36 { 37 if(y&1) ret=mul(ret,x,z); 38 y>>=1,x=mul(x,x,z); 39 } 40 return ret; 41 } 42 ll ksm(ll x,int y) 43 { 44 ll ret=1;x%=mod; 45 while(y) 46 { 47 if(y&1) ret=ret*x%mod; 48 x=x*x%mod,y>>=1; 49 } 50 return ret; 51 } 52 void pre() 53 { 54 for(int i=0;i<=3005;i++) 55 { 56 c[i][0]=1; 57 for(int j=1;j<=i;j++) 58 { 59 c[i][j]=c[i-1][j]+c[i-1][j-1]; 60 if(c[i][j]>=mod) c[i][j]-=mod; 61 } 62 } 63 for(int i=0;i<=3005;i++) 64 { 65 b[i]=i+1; 66 for(int j=0;j<i;j++) 67 { 68 b[i]=(b[i]-c[i+1][j]*b[j])%mod; 69 if(b[i]<0) b[i]+=mod; 70 } 71 b[i]=b[i]*ksm(c[i+1][i],mod-2)%mod; 72 } 73 for(int i=0;i<=3005;i++) 74 { 75 a1=ksm(i+1,mod-2); 76 for(int j=0;j<=i;j++) 77 t[i][j]=c[i+1][j]*b[j]%mod*a1%mod; 78 } 79 for(int i=2;i<=1000;i++) 80 for(int j=2;j<=1000;j++) 81 if(i*j<=1000) 82 vi[i*j]=1; 83 for(int i=2;i<=1000;i++) 84 if(!vi[i]) p[++cnt]=i; 85 } 86 inline ll calc(ll x) 87 { 88 ll ret=0,q=x%mod; 89 ll nt=q; 90 for(int i=bi;i>=0;i--) 91 { 92 ret=(ret+t[bi][i]*nt)%mod; 93 nt=nt*q%mod; 94 } 95 return ret*ksm(q,bi)%mod; 96 } 97 bool Miller_Rabin(ll x) 98 { 99 if(x==1) return 0; 100 for(int j=0;j<=9;j++) 101 if(x==prime[j]) 102 return 1; 103 ll y=x-1; 104 int k=0; 105 while((y&1)==0) k++,y>>=1; 106 for(int i=0;i<10;i++) 107 { 108 ll z=rand()%(x-1)+1; 109 ll c=qpow(z,y,x),d; 110 for(int j=0;j<k;j++) 111 { 112 d=mul(c,c,x); 113 if(d==1&&c!=1&&c!=x-1) return 0; 114 c=d; 115 } 116 if(d!=1) return 0; 117 } 118 return 1; 119 } 120 ll pollard_rho(ll x,ll y) 121 { 122 ll i=1,k=2,c=rand()%(x-1)+1; 123 ll d=c; 124 while(1) 125 { 126 i++; 127 c=(mul(c,c,x)+y)%x; 128 ll g=gcd((d-c+x)%x,x); 129 if(g!=1&&g!=x) return g; 130 if(c==d) return x; 131 if(i==k) d=c,k<<=1; 132 } 133 } 134 void find(ll x,int tim) 135 { 136 if(x==1) return; 137 if(Miller_Rabin(x)) 138 { 139 fac[++tot]=x; 140 return; 141 } 142 ll yz=x,mk=tim; 143 while(yz>=x) yz=pollard_rho(yz,tim--); 144 find(yz,mk),find(x/yz,mk); 145 } 146 ll solve(ll x,ll y,ll z) 147 { 148 ll ret=0,qwq,s1,s2; 149 tot=sp=0; 150 for(int j=1;j<=cnt;j++) 151 while(x%p[j]==0) 152 fac[++tot]=p[j],x/=p[j]; 153 find(x,120); 154 sort(fac+1,fac+tot+1); 155 for(int i=1;i<=tot;i++) 156 { 157 if(fac[i]!=fac[i-1]) sp++,pi[sp]=pk[sp]=fac[i],s[sp]=1; 158 else pk[sp]*=fac[i],s[sp]++; 159 } 160 for(int i=1;i<=sp;i++) 161 { 162 pw[i][0]=1,a1=ksm(pi[i],y); 163 for(int j=1;j<=s[i];j++) 164 pw[i][j]=pw[i][j-1]*a1%mod; 165 } 166 for(int i=0;i<=z;i++) 167 { 168 qwq=1; 169 for(int j=1;j<=sp;j++) 170 { 171 s1=s2=0,pw1[0]=1,a1=ksm(pi[j],z),a2=ksm(pi[j],z-i+1); 172 for(int k=1;k<=s[j];k++) pw1[k]=pw1[k-1]*a2%mod; 173 for(int k=0;k<=s[j];k++) s1=(s1+pw[j][k]*pw1[s[j]-k])%mod; 174 for(int k=0;k<s[j];k++) s2=(s2+pw[j][k]*pw1[s[j]-k-1]%mod*a1)%mod; 175 qwq=qwq*(s1-s2+mod)%mod*ksm(pk[j],z)%mod; 176 } 177 ret=(ret+qwq*t[z][i])%mod; 178 } 179 return ret; 180 } 181 int main() 182 { 183 pre(); 184 ca=read(); 185 for(int i=1;i<=ca;i++) 186 { 187 n=read(),ai=read(),bi=read(); 188 if(ai==bi) printf("%lld\n",calc(n)); 189 else printf("%lld\n",solve(n,ai,bi)); 190 } 191 return 0; 192 }
-
每一个你不满意的现在，都有一个你没有努力的曾经。