AtCoder-Sum of gcd of Tuples

AtCoder-Sum of gcd of Tuples (Hard) 题意 https://atcoder.jp/contests/abc162/tasks/abc162_e 求$\sum\gcd(a_1,a_2,\cdots,a_n)$,其中$a_i\in[1,K]$ Solution 直接计算肯定是不好计算的,可以考虑按$gcd$的值进行分类,问题就转化为一个计数问题 $\displaystyle \gcd(a,b)=d\Rightarrow\gcd(\frac{a}{d},\frac{b}{d})=1$ $\displaystyle {\gcd(a_1,\cdots, a_n)=d的数量}={\gcd(\frac{a_1}{d},\cdots,\frac{a_n}{d})=1的数量}$ 那么, $$ Ans=\sum_{d=1}^{K}dF(\lfloor\frac{K}{d}\rfloor, N) $$ 其中$F(K,N)$表示$\gcd(a_1,\cdots,a_N)=1$的个数,$a_i\in[1,K]$ 可以用容斥算出$\displaystyle F(K,N)=K^N-\sum_{i=1}^{K}F(\lfloor\frac{K}{i}\rfloor)$ #include <cstdio> #include <stack> #include <set> #include <cmath> #include <map> #include <time.h> #include <vector> #include <iostream> #include <string> #include <cstring> #include <algorithm> #include <memory.h> #include <cstdlib> #include <queue> #include <iomanip> // #include <unordered_map> #define P pair<int, int> #define LL long long #define LD long double #define PLL pair<LL, LL> #define mset(a, b) memset(a, b, sizeof(a)) #define rep(i, a, b) for (int i = a; i < b; i++) #define PI acos(-1.0) #define random(x) rand() % x #define debug(x) cout << #x << " " << x << "\n" using namespace std; const int inf = 0x3f3f3f3f; const LL __64inf = 0x3f3f3f3f3f3f3f3f; #ifdef DEBUG const int MAX = 1e6 + 50; #else const int MAX = 1e6 + 50; #endif const int mod = 1e9 + 7; LL N,K; LL f[MAX]; LL fact[MAX]; inline LL add(LL x, LL y){ LL res = x + y; return res >= mod ? res - mod : res; } inline LL qpow(LL x, LL n){ LL res = 1; while (n) { if(n &1) res = res * x % mod; x = x * x % mod; n >>= 1; } return res; } LL F(LL K, LL N){ if(f[K]) return f[K]; LL &res =f[K]; if(K==1){ return res = 1; } // res = qpow(K, N); res = fact[K]; for(LL i = 2, j; i <= K; i=j+1){ j = K/(K/i); // res = add(res, mod-F(K/i, N)); LL tmp = (j-i+1LL) * F(K/i, N) % mod; res = add(res, mod-tmp); } return res; } int main(){ #ifdef DEBUG freopen("in", "r", stdin); #endif scanf("%lld%lld", &N, &K); for(LL i = 1; i <= K; i++) fact[i] = qpow(i, N); LL ans = 0; f[1] = 1LL; for(LL k= 1; k <= K; k++){ F(k, N); } for(LL i = 1, j; i <= K; i++){ ans += f[K/i] % mod * i % mod; if(ans >= mod) ans -= mod; } printf("%lld\n", ans); }

📝 April 25, 2020&nbsp;·&nbsp;⌛ 2 min

CSAPP arch lab

arch lab Download archlab-handout 安装模拟器 解决undefined reference to ’matherr‘ 参考 Y86-64模拟器的安装与出现对’matherr’未定义引用问题的解决 Part A 在这部分要在sim/misc中完成,我们要编写和模拟三个Y86-64程序 sum.ys: 遍历链表求和 # begin at 0 .pos 0 irmovq stack,%rsp call main halt .align 8 ele1: .quad 0x00a .quad ele2 ele2: .quad 0x0b0 .quad ele3 ele3: .quad 0xc00 .quad 0 main: irmovq ele1,%rdi call rsum ret rsum: xorq %rax,%rax andq %rdi,%rdi je L1 pushq %rbx mrmovq (%rdi),%rbx mrmovq 8(%rdi),%rdi call rsum addq %rbx,%rax popq %rbx L1: ret .pos 0x200 stack: 用YAS编译后,再用YIS运行 ...

📝 March 16, 2020&nbsp;·&nbsp;⌛ 3 min

CodeChef - PRIMEDST

Prime Distance On Tree 题意 Prime Distance On Tree 给个树,从树上随机选取一对点$u,v$,求$\delta(u,v)$是素数的概率 Solution 可以从生成函数的角度考虑 假设rt是一个树的树根,而且rt的深度是d,将树中节点的深度统计出来,记为$f_{rt,d}$,如果$u,v…$是rt的子节点,那么rt对答案的贡献就是生成函数中素数项的系数,那么问题就是怎么计算生成函数了 $$ \sum_{u,v \in son(rt)} f_{u,1} * f_{v, 1} $$ 这个式子中的素数项系数和与下式是相等的,计算$f$是很简单的 $$ \Big(f_{rt,0}^2 - \sum_{u \in son(rt) }f_{u,1}^2\Big) / 2 $$ 为了确保复杂度不会太高,需要用点分治 #include <cstdio> #include <stack> #include <set> #include <cmath> #include <map> #include <time.h> #include <vector> #include <iostream> #include <string> #include <cstring> #include <algorithm> #include <memory.h> #include <cstdlib> #include <queue> #include <iomanip> // #include <unordered_map> #define P pair<int, int> #define LL long long #define LD long double #define PLL pair<LL, LL> #define mset(a, b) memset(a, b, sizeof(a)) #define rep(i, a, b) for (int i = a; i < b; i++) #define PI acos(-1.0) #define random(x) rand() % x #define debug(x) cout << #x << " " << x << "\n" using namespace std; const int inf = 0x3f3f3f3f; const LL __64inf = 0x3f3f3f3f3f3f3f3f; #ifdef DEBUG const int MAX = 2e3 + 50; #else const int MAX = 1e5 + 50; #endif const int mod = 1e9 + 7; #include <complex> namespace Prime { vector<int> prime; bool isprime[MAX<<1]; void init(int n){ mset(isprime, 1); isprime[0] = isprime[1] = 0; for(int i = 2; i < n; i++){ if(isprime[i]){ prime.push_back(i); } for(int j = 0; j < prime.size() and prime[j] * i < n; j++){ isprime[i * prime[j]] = 0; if(i * prime[j] == 0) break; } } } } // namespace namePrimPrime namespace FFT { complex<double> a[MAX<<2], b[MAX<<2]; int rev[MAX]; void fft(complex<double> *a, int n, int inv){ int bit = 0; while((1<<bit) < n) bit++; for(int i = 0; i < n; i++){ rev[i] = (rev[i>>1]>>1) | ((i & 1) << (bit-1)); if(i < rev[i]) swap(a[i], a[rev[i]]); } for(int mid = 1; mid < n; mid <<= 1){ complex<double> wn(cos(PI/mid), inv*sin(PI/mid)); for(int i = 0; i < n; i += (mid<<1)){ complex<double> w(1, 0); for(int j = 0; j < mid; j++, w *= wn){ complex<double> u = a[i+j], t = w * a[i+j+mid]; a[i+j] = u + t, a[i+j+mid] = u - t; } } } } void product(int n){ fft(a, n, 1); fft(b, n, 1); for(int i = 0; i <= n; i++) a[i] *= b[i]; fft(a, n, -1); } } // namespace nameFFTFFT vector<int> E[MAX]; bool vis[MAX]; // LL dep[MAX]; LL dep[MAX], dep_num[MAX]; int cnt; int siz[MAX]; int dfsG(int u, int p, int& Min, int &rt, int n){ siz[u] = 1; int Max = -1; for(int i = 0; i < E[u].size(); i++){ int v = E[u][i]; if(vis[v] or v == p) continue; siz[u] += dfsG(v, u, Min, rt, n); Max = max(Max, siz[v]); } Max = max(Max, n - siz[u]); if(Max < Min){ Min = Max, rt = u; } return siz[u]; } int dfsN(int u, int p){ int res = 1; for(int i = 0; i < E[u].size(); i++){ int v = E[u][i]; if(vis[v] or v == p) continue; res += dfsN(v, u); } return res; } int find_rt(int u){ int Min = inf; int rt = -1; int n = dfsN(u, -1); dfsG(u, -1, Min, rt, n); return rt; } LL dfs(int u, int p, int d){ // return max_dep dep[cnt++] = d; LL max_dep = d; for(int i = 0; i < E[u].size(); i++){ int v = E[u][i]; if(v == p or vis[v]) continue; max_dep = max(max_dep, dfs(v, u, d+1)); } return max_dep; } LL calcu(int u, int d){ // mset(dep, 0); cnt = 0; dfs(u, -1, d); LL Max = 0; for(int i = 0; i < cnt; i++){ dep_num[dep[i]]++; Max = max(Max, dep[i]); } LL res = 0; int N = 1; while(N <= Max * 2) N <<= 1; for(int i = 0; i < N; i++){ if(i <= Max) FFT::a[i] = FFT::b[i] = complex<double>(dep_num[i], 0); else FFT::a[i] = FFT::b[i] = complex<double>(0, 0); } FFT::product(N); for(int i = 0; i < Prime::prime.size() and Prime::prime[i] <= 2*Max; i++){ res += (LL)(FFT::a[Prime::prime[i]].real() / N + 0.5); } for(int i = 0; i < cnt; i++) dep_num[dep[i]]--; return res; } LL solve(int u){ // vis[u] = 1; int rt = find_rt(u); vis[rt] = 1; LL res = calcu(rt, 0); for(int i = 0; i < E[rt].size(); i++){ if(vis[E[rt][i]]) continue; res -= calcu(E[rt][i], 1); res += solve(E[rt][i]); } // for(int i = 0; i < E[rt].size(); i++){ // if(vis[E[rt][i]]) continue; // res += solve(E[rt][i]); // } return res; } int main(){ #ifdef DEBUG freopen("in", "r", stdin); #endif int n; scanf("%d", &n); Prime::init(MAX << 1); for(int i = 1; i < n; i++){ int u, v; scanf("%d%d", &u, &v); E[u].push_back(v); E[v].push_back(u); } LL up = solve(1); LL down = (LL)n * (n-1); printf("%.7lf\n", 1.0 * up / down); return 0; }

📝 March 1, 2020&nbsp;·&nbsp;⌛ 4 min

波里亚计数

Polya计数 将波利亚计数定理整理在这里,作为一个总结和介绍,也方便以后复习 为什么学习Polya计数定理 通过Polya计数定理,我们可以计算等价类的数量,比如下面这个问题: 用$m$种颜色给一个正方形染色,如果正方形可以自由转动,求染色方案数 让我们从一些概念开始 1 等价关系 1.1 等价关系的定义 假设$V$是一个集合,$S$是定义在$V$上的一个关系,若$S$有如下性质: 自反性 传递性 对称性 那么, $S$就是一个等价关系 $a$和$b$有关系$S$,可以记为$aSb$ 假设定义关系$S$,图形$a$可以旋转得到$b$ $\Leftrightarrow$ $aSb$ 例如图中的$方块_1$和$方块_2$具有关系$S$,即他们可以通过旋转得到彼此,而$方块_1$和$方块_2$则没有关系$S$ 1.2 等价类 通过上图,可以看出来:$方块_1$和$方块_2$是同一类的,而$方块_1$和$方块_2$则是另外两类,于是可以想到集合$V$上的等价关系$S$将集合的元素划分到不同的类中,我们把它称为等价类 而包含元素$a$的等价类则是由满足$aSb$的所有元素$b$组成的(当然也包含元素$a$),即$C(a)={b\in V | \ aSb}$ 仔细想一想,不难发现两个不同等价类是不相交的 2 置换群 2.1 置换群的定义 假设$A={1,2,\dotsc, n}$,通过置换,将$A$中的元素重新排列,得到另一个排列$a_1, a_2, \dotsc, a_n$,可以把这个过程写成 $$ \left ( \begin{matrix} 1 & 2 & \dotsc & n \newline a_1 & a_2 & \dotsc & a_n \end{matrix} \right ) $$ 所以,可以将置换看成一个双射函数$f:{1,2,\dotsc, n} \rightarrow {1,2,\dotsc, n}$ 置换之间也可进行合成运算 $$ \pi_1 = \left ( \begin{matrix} 1 & 2 & 3 & 4 \newline 4 & 2 & 1 & 3 \end{matrix} \right ) \ \ \pi_2 = \left ( \begin{matrix} 1 & 2 & 3 & 4 \newline 2 & 1 & 4 & 3 \end{matrix} \right ) $$ $$ \pi_1 \circ \pi_2 =\left ( \begin{matrix} 1 & 2 & 3 & 4 \newline 2 & 4 & 3 & 1 \end{matrix} \right ) $$ $\pi_1 \circ \pi_2 (1) = \pi_1(\pi_2(1))=2$ ...

📝 February 14, 2020&nbsp;·&nbsp;⌛ 5 min

Rock Paper Scissors Lizard Spock

Rock Paper Scissors Lizard Spock 题意: Rock Paper Scissors Lizard Spock 有五种手势,类似于石头剪刀布,有两个串$s, t$,由这五种手势组成,从某个位置开始匹配,如果$t_i$能赢$s_j$得一分,求一个$pos(0\le pos \le len(s)-len(t))$,使得得分最多 Solution: 将上图记为$G$,如果$op_1$可以赢$op_2$,则$G(op_1, op_2)=1$, 枚举可以得分的手势,假设当前手势为$op$, $$ \begin{aligned} t’_i &=(reverse\ t_i == op ? 1 : 0) \newline s’_i &= G[op][s_i] \end{aligned} $$ 将$t’,s’$做卷积,累加每次的结果,再遍历一遍匹配的起始位置,取最大值 #include <cstdio> #include <stack> #include <set> #include <cmath> #include <map> #include <time.h> #include <vector> #include <iostream> #include <string> #include <cstring> #include <algorithm> #include <memory.h> #include <cstdlib> #include <queue> #include <iomanip> #include <unordered_map> #define P pair<int, int> #define LL long long #define LD long double #define PLL pair<LL, LL> #define mset(a, b) memset(a, b, sizeof(a)) #define rep(i, a, b) for (int i = a; i < b; i++) #define PI acos(-1.0) #define random(x) rand() % x #define debug(x) cout << #x << " " << x << "\n" using namespace std; const int inf = 0x3f3f3f3f; const LL __64inf = 0x3f3f3f3f3f3f3f3f; #ifdef DEBUG const int MAX = 2e3 + 50; #else const int MAX = 1e6 + 50; #endif const int mod = 1e9 + 7; #include <complex> namespace FFT { int rev[MAX * 4]; int cnt[MAX * 4]; void fft(complex<double> *a, int n, int inv){ int bit = 0; while((1<<bit) < n) bit++; for(int i = 0; i < n; i++){ rev[i] = (rev[i>>1]>>1) | ((i & 1) << (bit-1)); if(i < rev[i]) swap(a[i], a[rev[i]]); } for(int mid = 1; mid < n; mid <<= 1){ complex<double> wn(cos(PI/mid), inv*sin(PI/mid)); for(int i = 0; i < n; i += (mid<<1)){ complex<double> w(1, 0); for(int j = 0; j < mid; j++, w*=wn){ complex<double> u = a[i+j], t = w * a[i+j+mid]; a[i+j] = u + t, a[i+j+mid] = u-t; } } } } }; int N; char s[MAX], t[MAX], rev_t[MAX]; complex<double> a[MAX*4], b[MAX*4]; int mark[100][100]; void init(){ mark['S']['P'] = 1; mark['P']['R'] = 1; mark['R']['L'] = 1; mark['L']['K'] = 1; mark['K']['S'] = 1; mark['S']['L'] = 1; mark['L']['P'] = 1; mark['P']['K'] = 1; mark['K']['R'] = 1; mark['R']['S'] = 1; } void solve(char op, int n, int m){ for(int i = 0; i < n; i++) a[i] = complex<double>(mark[op][s[i]], 0); for(int i = n; i < N; i++) a[i] = complex<double>(0, 0); for(int i = 0; i < m; i++) b[i] = complex<double>(rev_t[i]==op?1:0, 0); for(int i = m; i < N; i++) b[i] = complex<double>(0, 0); ; FFT::fft(a, N, 1); FFT::fft(b, N, 1); for(int i = 0; i < N;i++) a[i] *= b[i]; FFT::fft(a, N, -1); for(int i = m-1; i <= n-1; i++) FFT::cnt[i] += (int)(a[i].real() / N + 0.5); return; } int main(){ #ifdef DEBUG freopen("in", "r", stdin); #endif init(); scanf("%s%s", s, t); int n = strlen(s), m = strlen(t); N = 1; while(N <= n+m) N <<= 1; for(int i = 0; i < m; i++) rev_t[i] = t[m-i-1]; vector<char> ops({'R','P', 'S', 'K', 'L'}); for(auto op : ops){ solve(op, n, m); } int ans = 0; for(int i = m-1; i <= n-1; i++){ ans = max(ans, FFT::cnt[i]); // printf("%d ", FFT::cnt[i]); } // puts(""); printf("%d\n", ans); return 0; }

📝 February 11, 2020&nbsp;·&nbsp;⌛ 3 min