HDU 5628 Clarke and math
本文属于一个总结了一堆做法的玩意......
题目
简单的一个式子：给定$n,k,f(i)$，求
然后数据范围不重要，重要的是如何优化这个做法。
这个式子有$n$种问法，而且可以变式扩展，所以说这个式子也是比较重要的：
我们约定如果给定了$n,k$那么我们的$g$写作$g_k(n)$，如果给定了$n,k$中间的任意一个，枚举另一个，或者另一个是变化的，那么另一个数记为$i,j$
- 把$1~n$或$1~k$的$g_k(i)$或$g_i(n)$都求出来（狄利克雷卷积或者dp）
- 快速求单个$g_k(n)$（乘法同构，某考试的题（被我用dp水了））
- 所有$f$都是1或者某个特定函数（Luogu月赛题）
- 等等......
下面算法单词回答部分有些由于模数过大而要动态求逆元的原因需要多个$log$
为了区分我枚举的$k$和题目中的$k$，这里我们约定题目中的$k$写作大写的$K$
算法
狄利克雷卷积
比较大众的算法，就是根据狄利克雷卷积（$(f*g)(n)=sum_{d|n}{f(d)*g(frac{n}{d})}=sum_{a*b=n}{f(a)*g{b}}$（这个$f,g$与题目无关））。
我们首先令$k=1$，然后$g_1(i)=sum_{i_1|i}{f(i_1)}=(f*h)(n)$，这个时候我们推一下知道$h(i)=1$了
同时$k=2$时$g_2(i)=(f*h*h)(n)=(f*h^2)(n)$
答案就是$(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 ? ' ' : ' '); 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 ? ' ' : ' '); 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 ' ' 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 ? ' ' : ' '); 162 } 163 } 164 Main_Init(); 165 return 0; 166 }