【GDOI2020模拟3.18】树上的数（树形、状压dp+min_25筛）

题目大意：

(nle100,mle10^{10})，树是随机的。

题解：

考虑设(f(x))表示乘积为(x)的权值和。

不难发现(f(x))是一个积性函数，且(f(p)=p*n)。

如果能够快速求出(f(p^k))，那么就可以min_25筛了。

设(cnt[x][y])表示：

在树上分配指数(a[i])，指数和(sum a[i]=x)，(sum_{p1,p2,…,pt是一条路径} min(a[p1],a[p2],…,a[pt])=y)，的方案数。

如果能求出(cnt[k][...])，那么就可以计算(f(p^k)=sum cnt[k][i]*p^i)。

对于(cnt[x][y])，先考虑一下(y)的取值范围。

对于每一个点(i)，出去的路径的(sum min)显然(<=sum min(a[i],a[j])<=x)，所以(y<=x*(x+1)/2)

那么设(f[i][j][S][v])表示：

以(i)为根子树中，选的指数和是(j)，每个点到(i)的路径(min)的可重集是(S)，子树内的路径(sum~min)是(v)。

(j<=log_3m)，(S)需要把(0)去掉。

指数和最多是(20)，此时(S)的方案数是(2714)种。

直接dp看似不可过，事实上把有值的位置提出来就可以过。

由([S1][v1]、[S2][v2])转移到([S'][v'])可以预处理后实现(O(1))。

Code：

#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, _b = y; i <= _b; i ++)
#define ff(i, x, y) for(int i = x, _b = y; i <  _b; i ++)
#define fd(i, x, y) for(int i = x, _b = y; i >= _b; i --)
#define ll long long
#define pp printf
#define hh pp("
")
using namespace std;

const int mo = 998244353;

ll ksm(ll x, ll y) {
	ll s = 1;
	for(; y; y /= 2, x = x * x % mo)
		if(y & 1) s = s * x % mo;
	return s;
}

ll ni2 = ksm(2, mo - 2);

int n, k; ll m;
int x, y;

#define V vector<int>
#define pb push_back
#define si size()

V e[105];

void Init() {
	scanf("%d %lld", &n, &m);
	fo(i, 1, n - 1) {
		scanf("%d %d", &x, &y);
		e[x].pb(y); e[y].pb(x);
	}
	k = log2(m) / log2(3);
}

int fa[105];

ll cnt[21][405];

namespace sub1 {
	const int N = 3005;
	
	struct P {
		int a[21];
		P() { memset(a, 0, sizeof a);}
	};
	
	int num(P a) {
		int s = 0;
		fd(i, k, 1) {
			s *= 2;
			fo(j, 1, a.a[i]) s = s * 2 + 1;
		}
		return s;
	}
	
	P d[N]; int d0;
	
	int id[1 << 21];
	
	void dg(P a, int s, int x) {
		d[++ d0] = a;
		id[num(a)] = d0;
		fo(i, x, s) fo(j, 1, s / i) {
			a.a[i] += j;
			dg(a, s - i * j, i + 1);
			a.a[i] -= j;
		}
	}
	
	int zid[N][21], mv[N][N], mer[N][N];
	
	int sa[21], siz[N];
	
	void init() {
		P a = P();
		dg(a, k, 1);
		fo(i, 1, d0) {
			fo(j, 0, k) {
				P a = d[i];
				fo(u, j + 1, k) {
					if(j) a.a[j] += a.a[u];
					a.a[u] = 0;
				}
				zid[i][j] = id[num(a)];
			}
		}
		fo(i, 1, d0) {
			fo(u, 1, k) sa[u] = sa[u - 1] + d[i].a[u] * u;
			fo(j, 1, d0) {
				fo(u, 1, k) {
					mv[i][j] += d[j].a[u] * sa[u];
					mv[j][i] += d[j].a[u] * sa[u - 1];
				}
			}
		}
		fo(i, 1, d0) fo(j, 1, k) siz[i] += j * d[i].a[j];
		fo(i, 1, d0) fo(j, 1, d0) if(siz[i] + siz[j] <= k) {
			P a = d[i];
			fo(u, 1, k) a.a[u] += d[j].a[u];
			mer[i][j] = id[num(a)];
		}
	}
	
	struct Q {
		int p, q; ll v;
	};
	
	Q f[105][21][20005]; int f0[105][21];
	Q g[21][20005], h[21][20005]; int g0[21], h0[21];
	
	int bz[21][N][200], bx[21][N][200];
	
	void add(int j, Q a) {
		int &t = bx[j][a.p][a.q];
		if(!t) {
			t = ++ g0[j];
			g[j][t] = a;
		} else {
			g[j][t].v = (g[j][t].v + a.v) % mo;
		}
	}
	
	void dg(int x) {
		ff(_y, 0, e[x].si) {
			int y = e[x][_y];
			if(y == fa[x]) continue;
			fa[y] = x; dg(y);
		}
		
		fo(v, 0, k) {
			
			P a = P();
			if(v) a.a[v] = 1;
			Q b; b.p = id[num(a)]; b.q = v; b.v = 1;
			add(v, b);
			
			ff(_y, 0, e[x].si) {
				int y = e[x][_y];
				if(y == fa[x]) continue;
				
				fo(i, 0, k)	{
					h0[i] = g0[i];
					fo(j, 1, h0[i]) h[i][j] = g[i][j], g[i][j].v = 0;
				}
				fo(i, 0, k) fo(j, 1, h0[i]) {
					Q t1 = h[i][j];
					fo(p, 0, k - i) fo(q, 1, f0[y][p]) {
						Q t2 = f[y][p][q]; t2.p = zid[t2.p][v];
						add(i + p, (Q) {mer[t1.p][t2.p], t1.q + t2.q + mv[t1.p][t2.p], t1.v * t2.v % mo});
					}
				}
			}
			
			fo(i, 0, k) {
				fo(j, 1, g0[i]) {
					int p = g[i][j].p, q = g[i][j].q;
					int &t = bz[i][p][q];
					if(!t) {
						t = ++ f0[x][i];
						f[x][i][t] = g[i][j];
					} else {
						f[x][i][t].v = (f[x][i][t].v + g[i][j].v) % mo;
					}
					bx[i][p][q] = 0;
				}
				g0[i] = 0;
			}
		}
		fo(i, 0, k) {
			fo(j, 1, f0[x][i]) {
				int p = f[x][i][j].p, q = f[x][i][j].q;
				bz[i][p][q] = 0;
			}
		}
	}
	
	void build() {
		dg(1);
		fo(i, 0, k) fo(j, 1, f0[1][i]) {
			Q a = f[1][i][j];
			cnt[i][a.q] = (cnt[i][a.q] + a.v) % mo;
		}
	}
}


ll calc(ll p, ll k) {
	if(p == 2) return 0;
	ll s = 0, x = 1;
	fo(i, 1, k * (k + 1) / 2) {
		x = x * p % mo;
		s = (s + x * cnt[k][i]) % mo;
	}
	return s;
}

namespace sub2 {
	const int N = 2e5 + 5;
	
	int sq;
	int bz[N], p0; ll p[N];
	ll f[N], sf[N];
	
	void sieve(int n) {
		fo(i, 2, n) {
			if(!bz[i]) {
				p[++ p0] = i;
				f[p0] = i;
				sf[p0] = (sf[p0 - 1] + f[p0]) % mo;
			}
			for(int j = 1; i * p[j] <= n; j ++) {
				int k = i * p[j]; bz[k] = 1;
				if(i % p[j] == 0) break;
			}
		}
	}
	
	ll w[N]; int w0, i1[N], i2[N];
	ll g[N];
	
	ll c[N][35];
	
	ll dg(ll x, int y) {
		if(x <= 1 || p[y] > x) return 0;
		int k = x <= sq ? i1[x] : i2[m / x];
		ll s = g[k] - sf[y - 1] * n;
		for(int i = y; i <= p0 && p[i] * p[i] <= x; i ++) {
			ll p1 = p[i], p2 = p1 * p1;
			for(int j = 1; p2 <= x; p1 = p2, p2 *= p[i], j ++) {
				s = (s + dg(x / p1, i + 1) * c[i][j] + c[i][j + 1]) % mo;
			}
		}
		return s;
	}
	
	void work() {
		sq = sqrt(m);
		sieve(sq);
		for(ll i = 1, j; i <= m; i = j + 1) {
			j = m / (m / i);
			w[++ w0] = m / i;
			if(w[w0] <= sq) i1[w[w0]] = w0; else i2[m / w[w0]] = w0;
		}
		fo(i, 1, w0) {
			g[i] = ((w[i] + 2) % mo) * ((w[i] - 1) % mo) % mo * ni2 % mo;
		}
		for(int j = 1; j <= p0; j ++) {
			for(int i = 1; i <= w0 && p[j] * p[j] <= w[i]; i ++) {
				ll y = w[i] / p[j];
				int k = (y <= sq ? i1[y] : i2[m / y]);
				g[i] = (g[i] - (g[k] - sf[j - 1]) * f[j]) % mo;
			}
		}
		
		fo(i, 1, w0) {
			if(w[i] >= 2) g[i] -= 2;
			g[i] = g[i] * n % mo;
		}
		f[1] = 0; fo(i, 1, p0) sf[i] -= 2;
		
		fo(i, 1, p0) {
			int k = log2(m) / log2(p[i]);
			fo(j, 1, k) c[i][j] = calc(p[i], j);
		}
		
		ll ans = dg(m, 1);
		ans = (ans % mo + mo) % mo;
		pp("%lld
", (ans + 1) % mo);
	}
}

int main() {
	freopen("number.in", "r", stdin);
	freopen("number.out", "w", stdout);
	Init();
	sub1 :: init();
	sub1 :: build();
	sub2 :: work();
}