Solution to 55e of Niuke challenge

I thought about it for nearly a whole day without any result , Asking big guys to be appointed is the nesting of routine questions ......

First, use Mobius function to deal with gcd, set up f（n） It's long for n Sequence gcd And , Obviously, the original problem can be deduced from this thing

Obviously the first 2 This item is difficult to handle , Consider fixed

Now use the second routine , Mobius function and Id The Dirichlet convolution of a function is an Euler function

At this point, the form is very concise , however m A very large , It is still difficult to do , Notice that it can be partitioned , and phi The sum of can be determined by min25 Sieve or Du Jiao sieve （ The third way ？XD） Quickly sift out the prefix and of the root sign positions （ I learned this step min25 The practice of sifting out the root number positions ）

And then convert it into

B about 2 Double root sign m, The question at this point is how to treat each a = 1 ~ n seek f(a), The fourth routine goes online , Generating function plus divide and conquer ntt Additive inversion

Introduce the generating function ,    ( For the convenience of the later convergent form , Subscripts start from zero )

be f(a) amount to  

take h(t) Write in convergent form , by

Equivalent to B A polynomial shaped like a fraction , Find the polynomials they add up to in mod x^(n + 1) The polynomials under , This thing can be divided into two parts ntt Add the inverse to do , Use divide and conquer ntt, We can work out ：

among

Then on F Find an inverse , let G Multiply F^(-1), You can get the answer polynomial

In the last step f Array push answer is very simple , It is sufficient to investigate the one bit prolongation to increase the contribution .

Code ：（ As if a、b Contrary , But it doesn't matter , There's another detail phi(1) = 1）

//#define LOCAL
#include<bits/stdc++.h>
#define pii pair<int,int>
#define fi first
#define sc second
#define pb push_back
#define ll long long
#define trav(v,x) for(auto v:x)
#define all(x) (x).begin(), (x).end()
#define VI vector<int>
#define VLL vector<ll>
#define pll pair<ll, ll>
#define double long double
//#define int long long
using namespace std;
const int N = 2e6 + 100;
const int inf = 1e9;
//const ll inf = 1e18
const ll mod = 23068673;

#ifdef LOCAL
void debug_out(){cerr << endl;}
template<typename Head, typename... Tail>
void debug_out(Head H, Tail... T)
{
	cerr << " " << to_string(H);
	debug_out(T...);
}
#define debug(...) cerr << "[" << #__VA_ARGS__ << "]:", debug_out(__VA_ARGS__)
#else
#define debug(...) 42
#endif

ll n, m;

namespace siever
{
	const int N = 3e5 + 100;
	ll n, B, num;
	ll val[N], g[2][N], sum[2][N], f[N], pre[N];
	int ps1[N], ps2[N];
	
	int pos(ll x)
	{
		if(x <= B)
			return ps1[x];
		return ps2[n / x];
	}
	
	int pnum;
	ll pri[N];
	bool np[N];

	void work(ll _n)
	{
		n = _n;
		B = sqrt(n);
		for(int i = 2; i <= B; i++)
		{
			if(!np[i])
			{
				pri[++pnum] = i;
				sum[0][pnum] = pnum;
				sum[1][pnum] = (sum[1][pnum - 1] + i) % mod;
			}
			for(int j = 1; j <= pnum && i * pri[j] <= B; j++)
			{
				np[i * pri[j]] = 1;
				if(i % pri[j] == 0)
					break;
			}
		}
		for(ll l = 1, r; l <= n; l = r + 1)
		{
			ll x = n / l;
			r = n / x;
			val[++num] = x;
			if(x <= B)
				ps1[x] = num;
			else ps2[n / x] = num;
			x %= mod;
			g[0][num] = (x - 1 + mod) % mod;
			g[1][num] = (x - 1 + mod) % mod * (x + 2) / 2 % mod;
		//	cerr << "!!" << ' ' << x << ' ' << g[0][num] << ' ' << g[1][num] << '\n';
		}
		for(int i = 1; i <= pnum; ++i)
		{
			for(int j = 1; j <= num && 1LL * pri[i] * pri[i] <= val[j]; ++j)
			{
				int bf = pos(val[j] / pri[i]);
				g[0][j] = (g[0][j] - g[0][bf] + sum[0][i - 1] + mod) % mod;
				g[1][j] = (g[1][j] - pri[i] * (g[1][bf] - sum[1][i - 1] + mod) % mod + mod) % mod;
			}
		}
		for(int i = 1; i <= num; i++)
			f[i] = (g[1][i] - g[0][i] + mod) % mod;
		for(int i = 1; i <= pnum; i++)
			pre[i] = (sum[1][i] - sum[0][i] + mod) % mod;
		for(int j = pnum; j >= 1; j--)
		{
			for(int i = 1; i <= num; i++)
			{
				if(1LL * pri[j] * pri[j] > val[i])
					break;
				ll tmp = pri[j];
				for(int e = 1; tmp * pri[j] <= val[i]; e++, tmp *= pri[j])
				{
					ll s = (tmp - tmp / pri[j]) % mod;
					ll t = (tmp * pri[j] - tmp) % mod;
					(f[i] += (s * (f[pos(val[i] / tmp)] - pre[j] + mod) + t)) %= mod;
				}
			}
		}
	}
}

ll qpow(ll x, ll y = mod - 2)
{
	ll res = 1;
	while(y)
	{
		if(y & 1)
			res = res * x % mod;
		x = x * x % mod;
		y >>= 1;
	}
	return res;
}

namespace Poly
{
	int wh[N], len, cc;
	void ntt(VLL &a, bool inv)
	{
		for(int i = 0; i < len; i++)
			if(i < wh[i])swap(a[i], a[wh[i]]);
		for(int l = 2; l <= len; l <<= 1)
		{
			int md = l >> 1;
			ll tp = qpow(3, (mod - 1) / l);
			for(int i = 0; i < len; i += l)
			{
				ll mo = 1;
				for(int j = 0; j < md; j++, mo = mo * tp % mod)
				{
					ll ha = mo * a[i + j + md] % mod;
					a[i + j + md] = (a[i + j] - ha + mod) % mod;
					a[i + j] = (a[i + j] + ha) % mod;
				}
			}
		}
		if(inv)
		{
			ll tmp = qpow(len);
			for(int i = 1; i < len / 2; i++)
				swap(a[i], a[len - i]);
			for(int i = 0; i < len; i++)
				a[i] = a[i] * tmp % mod;
		}
	}
	void pre_ntt(int l)
	{
		cc = 0, len = 1;
		while(len <= l)
			++cc, len <<= 1;
		for(int i = 1; i < len; i++)
		{
			wh[i] = (wh[i >> 1] >> 1) | ((i & 1) << (cc - 1));
		}
	}
	VLL operator *(VLL x, VLL y)
	{
		if(x.size() <= 50 && y.size() <= 50)
		{
			VLL z(x.size() + y.size() - 1);
			for(int i = 0; i < x.size(); i++)
				for(int j = 0; j < y.size(); j++)
					z[i + j] = (z[i + j] + x[i] * y[j]) % mod;
			return z;
		}
		pre_ntt(x.size() + y.size() - 2);
		int bf = x.size() + y.size() - 1;
		x.resize(len);
		y.resize(len);
		ntt(x, 0), ntt(y, 0);
		for(int i = 0; i < len; i++)
			x[i] = x[i] * y[i] % mod;
		ntt(x, 1);
		x.resize(bf);
		return x;
	}
	VLL operator +(VLL x, VLL y)
	{
		if(x.size() < y.size())
			swap(x, y);
		for(int i = 0; i < y.size(); i++)
			(x[i] += y[i]) %= mod;
		return x;
	}
	VLL A;
	void cal_inv(VLL &b, int n)
	{
		if(n == 1)
			return (void)(b[0] = qpow(A[0]));
		int mid = (n + 1) / 2;
		VLL a, c;
		pre_ntt(n + mid + mid);
		a.resize(len), c.resize(len);
		cal_inv(c, mid);
		pre_ntt(n + mid + mid);
		for(int i = 0; i < n; i++)
			a[i] = A[i];
		ntt(a, 0);
		ntt(c, 0);
		for(int i = 0; i < len; i++)
			a[i] = (2 * c[i] - c[i] * c[i] % mod * a[i] % mod + mod) % mod;
		ntt(a, 1);
		for(int i = 0; i < n; i++)
			b[i] = a[i];
	}
	
	VLL inv(VLL x, int n)
	{
		x.resize(n);
		A = x;
		x.clear();
		x.resize(n);
		cal_inv(x, n);
		return x;
	}
}
using namespace Poly;

#define pvv pair<VLL, VLL>

pvv cdq(int l, int r)
{
	if(l == r)
	{
		ll a = (siever::f[l] - siever::f[l + 1] + mod) % mod;
		if(l == siever::num)
			(a += 1) %= mod;
		ll b = (m / siever::val[l]) % mod;
		return pvv((VLL){a}, (VLL){1, mod - b});
	}
	int mid = (l + r) >> 1;
	pvv pl = cdq(l, mid);
	pvv pr = cdq(mid + 1, r);
	return pvv(pl.fi * pr.sc + pl.sc * pr.fi, pl.sc * pr.sc);
}
 
void sol()
{
	cin >> n >> m;
   // n = 3e5, m = 1e10;
	siever::work(m);
	pvv ans = cdq(1, siever::num);
//	trav(v, ans.sc)
//		cerr << v << ' ';
//	cerr << '\n';
	ans.sc = inv(ans.sc, n + 1);
	ans.fi = ans.fi * ans.sc;
//	trav(v, ans.fi)
//		cerr << v << ' ';
//	cerr << '\n';
	ll cur = 0, sum = 0;
	for(int i = 1; i <= n; i++)
	{
		cur = cur * m % mod;
		sum = (sum * m + ans.fi[i]) % mod;
		cur = (cur + sum) % mod;
		cout << cur << ' ';
	}
	cout << '\n';
}

signed main()
{
	ios::sync_with_stdio(0);
	cin.tie(0);
//	int tt;
//	cin >> tt;
//	while(tt--)
		sol();
}