1. 程式人生 > >BZOJ2137: submultiple(生成函數,二項式定理)

BZOJ2137: submultiple(生成函數,二項式定理)

mil for 數組 sam pro 由於 打出 led pri

Description

設函數g(N)表示N的約數個數。現在給出一個數M,求出所有M的約數x的g(x)的K次方和。

Input

第一行輸入N,K。N表示M由前N小的素數組成。接下來N行,第i+1行有一個正整數Pi,表示第Ai小的素數 有 Pi次。等式:

Output

輸出一個數,表示答案。只需輸出最後答案除以1000000007的余數。

Sample Input

2 3
1
3

Sample Output

900
【樣例說明】
M=2^1*3^3=54
M的約數有1,2,3,6,9,18,27,54.約數個數分別為1,2,2,4,3,6,4,8.
Answer=1^3+2^3+2^3+4^3+3^3+6^3+4^3+8^3=900


編號 N K Pi<=
1 50 3 10000
2 50 100 10000
3 50 20101125 10000
4 999 17651851 100000
5 5000 836954247 100000
6 4687 1073741823 100000
7 4321 123456789 100000
8 5216 368756432 100000
9 8080 2^31-1 100000
10 10086 3 2^63-1
11 64970 3 2^63-1
12 71321 3 2^63-1
13 350 5 2^31-1
14 250 6 2^31-1
15 110 7 2^31-1
16 99 8 2^31-1
17 80 9 2^31-1
18 70 10 2^31-1
19 60 11 2^31-1
20 50 12 2^31-1

解題思路:

拿到題感覺一臉不可做,反正不是反演的樣子。

先來考慮一下樣例:

$2^1*3^3=54$

考慮如何將答案分類。將貢獻重新累加起來。

首先$g(d)=\sum\limits_{i=1}^{n}{(1+c_i)}$()其中$c_i$為第$i$項的次數)

現在考慮如何將$M$的所有因數表示出來:

以54為例

其中每個因數可以按照2的次數分類,可以是0~1次共有2種選法。

對於每種選法,對應的數中3的次數有0~3共4種選法。

而對應3的每一種選法其答案都是$[(2的次數+1)*(3的次數+1)]^k$

那麽這種結果可以看做兩個多項式乘積的形式,這兩個多項式都是(1+2+3+...)的形式,而由於那個指數上的k是可以帶入括號的。

那麽答案就變成了$f(54)=(1^3+2^3)*(1^3+2^3+3^3+4^3)=9*100=900$

現在答案就可以解出來了,答案就變成了:

$f(M)=\prod\limits_{i=1}^{n}\sum\limits_{j=1}^{p_i}{j^k}$

n是可以枚舉的,對於前45分的點,可以暴力快速冪求解。

剩下的怎麽辦。

後面的那個高次求和$\sum\limits_{i=1}^{n}{i^k}$是非常公式化的,怎麽可能沒有公式。

高次求和公式會形如:

$\sum\limits_{i=1}^{n}{i^k}=\sum\limits_{i=1}^{k+1}{a_i*n^i}$

所以可以這樣求解(orz zwz):

$a_1x^1+a_2x^2+...+a_{k+1}x^{k+1}+(x+1)^k=a_1(x+1)+a_2(x+1)^2+...+a_{k+1}(x+1)^{k+1}$

二項式定理展開移項:

$\sum\limits_{i=1}^{k+1}{C_i^0x^0}+\sum\limits_{i=2}^{k+1}{C_i^1x^1}+...+\sum\limits_{i=k+1}^{k+1}{C_i^kx^k}=\sum\limits_{i=0}^{k}{C_k^ix^i}$

這個楊輝三角打出矩陣消元就可以解出系數了。

代碼:

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<algorithm>
 4 typedef long long lnt;
 5 const lnt mod=(lnt)(1e9+7);
 6 lnt n,k;
 7 lnt maxs;
 8 lnt C[100][100];
 9 lnt F[200001];
10 lnt p[200001];
11 lnt a[1001];
12 lnt b[100][100];
13 lnt c[100];
14 lnt ksm(lnt x,lnt y)
15 {
16     lnt ans=1;
17     while(y)
18     {
19         if(y&1)ans=ans*x%mod;
20         x=x*x%mod;
21         y=y/2;
22     }
23     return ans;
24 }
25 lnt squ(lnt x)
26 {
27     return x%mod*x%mod;
28 }
29 int main()
30 {
31     scanf("%lld%lld",&n,&k);
32     for(int i=1;i<=n;i++)
33         scanf("%lld",&p[i]),
34         maxs=std::max(p[i],maxs);
35     if(maxs<=100000)
36     {
37         lnt ans=1;
38         for(int i=1;i<=maxs+1;i++)F[i]=(F[i-1]+ksm(i,k))%mod;
39         for(int i=1;i<=n;i++)ans=(ans*F[p[i]+1])%mod;
40         printf("%lld\n",ans);
41         return 0;
42     }
43     if(k==3)
44     {
45         lnt ans=1;
46         lnt Inv=ksm(4,mod-2);
47         for(int i=1;i<=n;i++)
48             ans=ans*squ(p[i]%mod+1)%mod*squ(p[i]%mod+2)%mod*Inv%mod;
49         printf("%lld\n",ans);
50         return 0;
51     }
52     C[0][0]=1;
53     for(int i=1;i<=k+1;i++)
54     {
55         C[i][0]=1;
56         for(int j=1;j<=i;j++)
57             C[i][j]=(C[i-1][j-1]+C[i-1][j])%mod;
58     }
59     for(int i=1;i<=k+1;i++)
60     {
61         for(int j=i;j<=k+1;j++)
62         {
63             b[i][j]=C[j][i-1];
64         }
65     }
66     for(int i=1;i<=k+1;i++)a[i]=C[k][i-1];
67     for(int i=k+1;i>=1;i--)
68     {
69         a[i]=(a[i]*ksm(b[i][i],mod-2))%mod;
70         for(int j=i-1;j;j--)
71         {
72             (a[j]-=b[j][i]*a[i]%mod)%=mod;
73         }
74     }
75     lnt ans=1;
76     for(int i=1;i<=n;i++)
77     {
78         lnt tmp=0;
79         for(int j=1;j<=k+1;j++)
80             tmp=(tmp+(a[j]*ksm(p[i]%mod+1,j))%mod)%mod;
81         ans=(ans*tmp)%mod;
82     }
83     printf("%lld\n",(ans%mod+mod)%mod);
84     return 0;
85 }

BZOJ2137: submultiple(生成函數,二項式定理)