HDU 5628 Clarke and math——卷積，dp，組合

時間  2019-11-16
標籤 hdu clarke math 組合 欄目 應用數學 简体版
原文   原文鏈接

HDU 5628 Clarke and math

　　本文屬於一個總結了一堆作法的玩意......html

題目

　　簡單的一個式子：給定$n,k,f(i)$，求git

　　而後數據範圍不重要，重要的是如何優化這個作法。算法

　　這個式子有$n$種問法，並且能夠變式擴展，因此說這個式子也是比較重要的：數組

　　咱們約定若是給定了$n,k$那麼咱們的$g$寫做$g_k(n)$，若是給定了$n,k$中間的任意一個，枚舉另外一個，或者另外一個是變化的，那麼另外一個數記爲$i,j$函數

  1. 把$1~n$或$1~k$的$g_k(i)$或$g_i(n)$都求出來（狄利克雷卷積或者dp）
  2. 快速求單個$g_k(n)$（乘法同構，某考試的題（被我用dp水了））
  3. 全部$f$都是1或者某個特定函數（Luogu月賽題）
  4. 等等......

　　下面算法單詞回答部分有些因爲模數過大而要動態求逆元的緣由須要多個$log$優化

　　爲了區分我枚舉的$k$和題目中的$k$，這裏咱們約定題目中的$k$寫做大寫的$K$spa

算法

狄利克雷卷積

　　比較大衆的算法，就是根據狄利克雷卷積（$(f*g)(n)=\sum_{d|n}{f(d)*g(\frac{n}{d})}=\sum_{a*b=n}{f(a)*g{b}}$（這個$f,g$與題目無關））。code

　　咱們首先令$k=1$，而後$g_1(i)=\sum_{i_1|i}{f(i_1)}=(f*h)(n)$，這個時候咱們推一下知道$h(i)=1$了htm

　　同時$k=2$時$g_2(i)=(f*h*h)(n)=(f*h^2)(n)$blog

　　答案就是$(f*(h=1)^K)(i)$

　　因爲狄利克雷卷積知足結合律，因此能夠用快速冪優化，每次計算狄利克雷卷積複雜度$O(nlog_2n)$，快速冪$O(log_2n)$

　　代碼（早期）：

　　

 1 #include <cstdio>
 2 #include <cctype>
 3 #include <cstring>
 4 using namespace std;
 5 
 6 #define MAX_BUF 10000000
 7 #define MAXN 100005
 8 #define MODS 1000000007
 9 #define LL long long
10 
11 struct FR{
12     char buf[MAX_BUF + 1], *p1, *p2;
13     FR(){
14         p2 = (p1 = buf) + fread(buf, 1, MAX_BUF, stdin);
15     }
16     inline void flash(){
17         p2 = (p1 = buf) + fread(buf, 1, MAX_BUF, stdin);
18     }
19     inline char get_char(){
20         return p1 == p2 ? EOF : *p1++;
21     }
22 }fr;
23 
24 inline LL read(){
25     LL num = 0;
26     char c, sf = 1;
27     while(isspace(c = fr.get_char()));
28     if(c == '-') c = fr.get_char(), sf = -1;
29     while(num = num * 10 + c - 48, isdigit(c = fr.get_char()));
30     return num * sf;
31 }
32 
33 LL n, k, f[MAXN], g[MAXN], h[MAXN], tmp[MAXN]; 
34 
35 inline void Dirchlet(LL *scr, LL *dst){
36     memset(tmp, 0, sizeof(tmp));
37     for(int i = 1, j; i * i <= n; i++)
38         for((tmp[i * i] += scr[i] * dst[i] % MODS) %= MODS, j = i + 1; i * j <= n; j++)
39             (tmp[i * j] += scr[i] * dst[j] % MODS + scr[j] * dst[i] % MODS) %= MODS;
40     memcpy(dst, tmp, sizeof(tmp));
41 }
42 
43 int main(){
44     LL T = read();
45     while(T--){
46         n = read(), k = read();
47         for(LL i = 1; i <= n; i++)
48             f[i] = read(), h[i] = 0, g[i] = 1;
49         h[1] = 1;
50         while(k){
51             if(k & 1) Dirchlet(g, h);
52             Dirchlet(g, g);
53             k >>= 1;
54         }
55         Dirchlet(f, h);
56         for(LL i = 1; i <= n; i++) printf("%lld%c", h[i], i == n ? '\n' : ' ');
57     }
58     return 0;
59 }

 

dp+組合數學

　　之前看的：複雜度$O(nlog_2^2n)$（<--原文），之前拿這個水過題，例如快速求單個$g_k(n)$的時候效果比較好。

　　思路大概就是dp出因子的貢獻而後計算出組合，即某個因子在全部可能傳遞狀況下是如何傳遞的。

　　解：設$dp[i][j]$表示在$n=i$的時候，中間的$\Sigma$部分選擇了$j$個因子（原文是不一樣的因子，可是我並不理解）時對於第$i$個位置所產生的貢獻。

　　因此$dp[i][k] += dp[j][k - 1](j|i,k \leq log(n),dp[i][0] = f[i])$，舉個栗子：在最後一個$\Sigma$（第k層）的位置$i_{k}|i_{k-1}$，而第k-1層枚舉因子的時候就會有一個$x * i_{k} = i_{k - 1}$，這個貢獻最終是計入到了$k-1$上去的。可是從$f(j)$到$g(i)$中間的$\Sigma$在不一樣的位置消掉了$r$次因子不止一種方案，$f(j)$不止作了一次貢獻，這個方案數就是$C(K,r)$。

　　通俗地說：咱們求$g(i)$，可是這個$g(i)$是由多個$f(j)$湊出來的，如何將$g(i)$和$f(j)$扯上關係？就是中間$k$個$\Sigma$有些$\Sigma$會吃掉一些因子，假設$K$個$\Sigma$吃了$r$個因子，那麼也就是說有些$\Sigma$吃了個人因子，而有些沒有，這些狀況都會致使$f(j)$對$g(i)$作貢獻。所以若是如今我知道個人因子一共被$K$個$\Sigma$吃了$r$個的話那麼就至關於$K$個物品中選出$r$個物品，方案數$C(K,r)$。

　　也就是$g_K(n)=C(K,j)*dp[n][j](j \leq log(n))$

　　可是時間複雜度都會有瓶頸（dp是$O(nlog_2^2n)$，單詞詢問$O(log_2^2n)$），不過適用性很好。

　　代碼（早期）：

 1 #include <cstdio>
 2 #include <cctype>
 3 #include <cstring>
 4 using namespace std;
 5 
 6 #define MAX_BUF 10000000
 7 #define MAXK 25
 8 #define MAXN 100005
 9 #define MODS 1000000007
10 #define LL long long
11 
12 struct FR{
13     char buf[MAX_BUF + 1], *p1, *p2;
14     FR(){
15         p2 = (p1 = buf) + fread(buf, 1, MAX_BUF, stdin);
16     }
17     inline void flash(){
18         p2 = (p1 = buf) + fread(buf, 1, MAX_BUF, stdin);
19     }
20     inline char get_char(){
21         return p1 == p2 ? EOF : *p1++;
22     }
23 }fr;
24 
25 inline LL read(){
26     LL num = 0;
27     char c, sf = 1;
28     while(isspace(c = fr.get_char()));
29     if(c == '-') c = fr.get_char(), sf = -1;
30     while(num = num * 10 + c - 48, isdigit(c = fr.get_char()));
31     return num * sf;
32 }
33 
34 LL n, k, f[MAXN], c[MAXN][MAXK], dp[MAXN][MAXK], fac[MAXN];
35 
36 inline LL Get_Pow(LL base, LL t){
37     LL ret = 1;
38     while(t){
39         if(t & 1) (ret *= base) %= MODS;
40         (base *= base) %= MODS;
41         t >>= 1;
42     }
43     return ret;
44 }
45 
46 inline LL Get_C(LL a, LL b){
47     if(c[a][b]) return c[a][b];
48     LL u = fac[a], d = fac[a - b] * fac[b] % MODS;
49     return u * Get_Pow(d, MODS - 2) % MODS;
50 }
51 
52 const LL Up_K = 21;
53 
54 int main(){
55     LL T = read();
56     fac[0] = 1;
57     for(int i = 1; i < MAXN; i++) fac[i] = fac[i - 1] * i % MODS;
58     while(T--){
59         memset(dp, 0, sizeof(dp));
60         n = read(), k = read();
61         for(LL i = 1, j; i <= n; i++)
62             for(dp[i][0] = f[i] = read(), j = i << 1; j <= n; j += i)
63                 for(LL r = 0; r < Up_K; r++)
64                     (dp[j][r + 1] += dp[i][r]) %= MODS;
65         for(LL i = 1; i <= n; i++){
66             LL ans = 0;
67             for(LL r = 0; r < Up_K; r++)
68                 (ans += Get_C(k, r) * dp[i][r] % MODS) %= MODS;
69             printf("%lld%c", ans, i == n ? '\n' : ' ');
70         }
71     }
72     return 0;
73 }

乘法同構

　　這個我看不懂它的定義，可是乘法同構的思路好像和下面差很少，可是時間複雜度差了一些，可是這裏貼出原來的題和題解：對於上面式子在$n,k<=1e5$下詢問$1e3$次（能夠用dp+組合數學作）。

定義乘法同構：
$A=p[1]^{a[1]} * p[2]^{a[2]} * ... * p[n]^{a[n]}$
$B=q[1]^{b[1]} * q[2]^{b[2]} * ... * q[n]^{b[n]}$
其中p[i]與q[i]皆爲質數
將數組a與b降序排序後若是是徹底相同的，那麼稱A與B是乘法同構的
如 2*2*2*2*3*3*5 與 7*11*11*3*3*3*3 同構

咱們發現，10^5內在乘法同構下本質不一樣的數字只有165個
定義kind[x]:x在乘法同構下所屬於的種類
對着165個本質不一樣的數構造出矩陣A
咱們發現，$f_0[d]$對$f_k[x]$(d是x的因子)的貢獻爲$f_0[d] * A^k[kind[1]][kind[x/d]]$
因此咱們只須要A^k[kind[1]][]這一行就夠了
預處理$A,A^2,A^4,A^8,...O(165^3*log(n))$
對於每次詢問由於最終只須要一行，能夠優化到$O(165^2*log(n))$
總複雜度$O(165^3*log(n) + m*165^2*log(n))$

組合數學

　　這個最勁$O(nlog_2n)$，而後也能夠解決大部分問題，預處理$O(nlog_2n)$以及單次回答$O(log_2n)$，從HJWJBSR蒯來的.

　　並且目前在HDU上只有46ms（register優化居然失效了，速度還慢一些）。

　　思路：

　　咱們按照上一個dp+組合數學的思路，咱們設$h(j)$表示$f(i)$對$g(i*j)$的貢獻，若是算出了$h(1)~h(n)$的取值，那麼就能夠枚舉因數在$O(nlog_2n)$的時間內算出$f(1)~f(n)$。

　　如何計算$h(i)$：咱們首先發現$h(i)$只和$K$和$i*j$的因子有關，同時發現假設咱們咱們須要求$g(k_1^a)$和$g(k_2^a)$，在$K$同樣的狀況下，咱們傳遞因子的方案是相同的，也就是之前的$h(k^a)=C(a+K-1,a)$（上面一句話的$k,k_1,k_2$都是質因數）！而後因爲每一個因子傳遞的時候相互獨立！因此根據乘法原理，假設對要求的$h(n)$的$n$分解爲$n = k_1^{a_1}*k_2^{a_2}...k_x^{a_x}$，而後$h(n)=h(k_1^{a_1}*k_2^{a_2}...k_x^{a_x})=h(k_1^{a_1})*...*h(k_x^{a_x})$。咱們發現這是一個積性函數！能夠用線性篩求出來，因爲中間須要求某個數是某個質數的多少次冪，因此帶log。

　　代碼以下：

  1 #include <cstdio>
  2 #include <cctype>
  3 #include <cstring>
  4 
  5 //User's Lib
  6 
  7 using namespace std;
  8 
  9 char buf[11111111], *pc = buf;
 10 char outp[1111111], *op = outp;
 11 
 12 extern inline void Main_Init(){
 13     static bool INITED = false;
 14     if(INITED){
 15         fwrite(outp, 1, op - outp, stdout);
 16         fclose(stdin), fclose(stdout);
 17     } else {
 18         fread(buf, 1, 11111111, stdin); 
 19         INITED = true;           
 20     }
 21 }
 22 
 23 static inline int read(){
 24     int num = 0;
 25     char c;
 26     while((c = *pc ++) < 48);
 27     while(num = num * 10 + c - 48, (c = *pc ++) >= 48);
 28     return num;
 29 }
 30 
 31 inline void Write(const int &tar){//You Need To Write '-' and '\n' By Hand
 32     if(!tar) return ;
 33     Write(tar / 10);
 34     *op ++ = (tar % 10) ^ 48;
 35 }
 36 
 37 namespace LKF{
 38     template <typename T>
 39     extern inline T abs(T tar){
 40         return tar < 0 ? -tar : tar;
 41     }
 42     template <typename T>
 43     extern inline void swap(T &a, T &b){
 44         T t = a;
 45         a = b;
 46         b = t;
 47     }
 48     template <typename T>
 49     extern inline void upmax(T &x, const T &y){
 50         if(x < y) x = y;
 51     }
 52     template <typename T>
 53     extern inline void upmin(T &x, const T &y){
 54         if(x > y) x = y;
 55     }
 56     template <typename T>
 57     extern inline T max(T a, T b){
 58         return a > b ? a : b;
 59     }
 60     template <typename T>
 61     extern inline T min(T a, T b){
 62         return a < b ? a : b;
 63     }
 64 }
 65 
 66 //Source Code
 67 
 68 const int MAXN = 100023;
 69 const int MAXK = 18;
 70 const int MODS = 1000000007;
 71 
 72 int n, k;
 73 int fac[MAXN];
 74 int f[MAXN], h[MAXN], g[MAXN];
 75 
 76 inline int Get_Pow(int base, int times){
 77     int ret = 1;
 78     while(times){
 79         if(times & 1) ret = 1ll * ret * base % MODS;
 80         base = 1ll * base * base % MODS;
 81         times >>= 1;
 82     }
 83     return ret;
 84 }
 85 
 86 inline int Get_Inv(const int &tar){
 87     return Get_Pow(tar, MODS - 2);
 88 }
 89 
 90 int tot;
 91 int prime[MAXN];
 92 bool is_prime[MAXN];
 93 
 94 inline void Get_Prime(){
 95     for(int i = 2; i < MAXN; i++){
 96         if(!is_prime[i]) prime[++ tot] = i;
 97         for(int j = 1; j <= tot; j++){
 98             int tmp = i * prime[j];
 99             if(tmp >= MAXN) break;
100             is_prime[tmp] = true;
101             if(i % prime[j] == 0) break;
102         }
103     }
104 }
105 
106 inline void Get_Fac(){
107     fac[0] = 1;
108     for(int i = 1; i < MAXN; i++)
109         fac[i] = 1ll * fac[i - 1] * i % MODS;
110 }
111 
112 int c[MAXK];
113 
114 inline int Get_C(const int &a, const int &b){//a > b
115     int tmp = 1ll * fac[a - b] * fac[b] % MODS;
116     return 1ll * fac[a] * Get_Inv(tmp) % MODS;
117 }
118 
119 inline void Get_C(const int &k){
120     for(int i = 1; i < MAXK; i++)
121         c[i] = Get_C(i + k, i);
122 }
123 
124 inline void Get_H(){
125     int cnt = 0;
126     h[1] = 1;
127     for(int i = 2; i <= n; i++){
128         if(!is_prime[i]) h[i] = c[1], cnt ++;
129         for(int j = 1; j <= cnt; j++){
130             int tmp = i * prime[j];
131             if(tmp > n) break;
132             if(i % prime[j]) h[tmp] = 1ll * h[i] * h[prime[j]] % MODS;
133             else {
134                 int Cnt = 0;
135                 while(!(tmp % prime[j])) tmp /= prime[j], Cnt ++;//計算冪次
136                 h[i * prime[j]] = 1ll * h[tmp] * c[Cnt] % MODS;
137                 break;
138             }
139         }
140     }
141 }
142 
143 int main(){
144     Main_Init();
145     Get_Fac();
146     Get_Prime();
147     int T = read();
148     while(T --){
149         memset(g, 0, sizeof(int) * (n + 5));
150         n = read(), k = read() - 1;
151         for(int i = 1; i <= n; i++)//O(n)
152             f[i] = read();
153         Get_C(k);//O(log^2)
154         Get_H();//O(nlogn)
155         for(int i = 1; i <= n; i++)//O(n + n / 2 + ...) = O(nlogn)(?)
156             for(int j = 1; i * j <= n; j++)
157                 (g[i * j] += 1ll * f[i] * h[j] % MODS) %= MODS;
158         for(int i = 1; i <= n; i++){
159             if(g[i]) Write(g[i]);
160             else *op ++ = '0';
161             *op ++ = (i == n ? '\n' : ' ');
162         }
163     }
164     Main_Init();
165     return 0;
166 }
相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程