\(\mathcal{Description}\)
Link.
给定 \(\{f_1,f_2,\cdots,f_n\}\)，素数 \(p\)。求字典序最小的 \(\{a_1,a_2,\cdots,a_n\}\)，满足对于所有 \(i\in[1,n]\)，\(a_i\in\{0,1\}\) 并且
\[\sum_{\{k_{1..n}\}}[(\forall j)\left((a_j=0\land k_j=0)\lor(a_j\not=0\land k_j\ge0\right)]\left[ \sum_{j=1}^njk_j=i \right]\equiv f_i\pmod p
\]
\(n<2^{18}\)。
\(\mathcal{Solution}\)
对于某个 \(a_i=1\)，其 OGF 为 \(\frac{1}{1-x^i}\)，所有 OGF 之积的 \(1\sim n\) 次项系数在\(\bmod p\) 意义下就是 \(f_{1..n}\)。记 \(F(x)=\sum_{i=0}^{+\infty}f_ix^i\)，默认运算在\(\bmod p\) 意义下进行，那么：
\[F(x)=\prod_{i=1}^n\left( \frac{1}{1-x^i} \right)^{a_i}
\]
两边取 \(-\ln\)：
\[\Rightarrow~~~~-\ln F(x)=\sum_{i=1}^na_i\ln(1-x^i)
\]
Tayler 展开右式 \(\ln\)：
\[\Rightarrow~~~~-\ln F(x)=\sum_{i=1}^na_i\sum_{j=1}^{+\infty}-\frac{x^{ij}}j
\]
左右取负，枚举 \(T=ij\)：
\[\Rightarrow~~~~\ln F(x)=\sum_{T=1}^{+\infty}x^T\sum_{i|T}a_i\frac{i}T
\]
对已知的 \(F(x)\) 求 \(\ln\)，就能取得 \(T\in[1,n]\) 内的所有 \(\frac{1}T\sum_{i|T}a_ii\)，从 \(T=1\) 起刷表求解，每次枚举倍数消除贡献，可以 \(\mathcal O(n\ln n)\) 解出所有 \(a_{1..n}\)（没错，由于 \(\ln F(x)\) 模意义下的多项式系数唯一，\(\{a_n\}\) 唯一确定）。
综上，复杂度为多项式 \(\ln\) 的 \(\mathcal O(n\log n)\)。建议用毛爷爷的四次 FFT 科技做 MTT。
\(\mathcal{Code}\)
/* Clearink */#include <cmath>
#include <cstdio>
#include <cassert>#define rep( i, l, r ) for ( int i = l, rpbound##i = r; i <= rpbound##i; ++i )
#define per( i, r, l ) for ( int i = r, rpbound##i = l; i >= rpbound##i; --i )typedef long long LL;inline int rint () {
int x = 0; char s = getchar ();
for ( ; s < '0' || '9' < s; s = getchar () );
for ( ; '0' <= s && s <= '9'; s = getchar () ) x = x * 10 + ( s ^ '0' );
return x;
}inline void wint ( const int x ) {
if ( 9 < x ) wint ( x / 10 );
putchar ( x % 10 ^ '0' );
}const int MAXLEN = 1 << 19;
const double PI = acos ( -1 );
int n, M, f[MAXLEN + 5], a[MAXLEN + 5];inline int mul ( const long long a, const int b ) { return a * b % M; }
inline int sub ( int a, const int b ) { return ( a -= b ) < 0 ? a + M : a; }
inline int add ( int a, const int b ) { return ( a += b ) < M ? a : a - M; }
inline int mpow ( int a, int b ) {
int ret = 1;
for ( ; b; a = mul ( a, a ), b >>= 1 ) ret = mul ( ret, b & 1 ? a : 1 );
return ret;
}namespace PolyOper {int rev[MAXLEN + 5], inv[MAXLEN + 5];inline void initInv () {
inv[1] = 1;
rep ( i, 2, MAXLEN ) inv[i] = mul ( M - M / i, inv[M % i] );
}struct Complex {
double x, y;
Complex (): x ( 0 ), y ( 0 ) {}
Complex ( const double tx, const double ty ): x ( tx ), y ( ty ) {}
inline Complex operator + ( const Complex t ) const {
return Complex ( x + t.x, y + t.y );
}
inline Complex operator - ( const Complex t ) const {
return Complex ( x - t.x, y - t.y );
}
inline Complex operator * ( const Complex t ) const {
return Complex ( x * t.x - y * t.y, x * t.y + y * t.x );
}
inline Complex operator / ( const double t ) const {
return Complex ( x / t, y / t );
}
};
Complex omega[MAXLEN + 5], P[MAXLEN + 5], Q[MAXLEN + 5];
Complex C[MAXLEN + 5], D[MAXLEN + 5], E[MAXLEN + 5], F[MAXLEN + 5];inline void FFT ( const int n, Complex* A, const int tp ) {
rep ( i, 0, n - 1 ) if ( i < rev[i] ) std::swap ( A[i], A[rev[i]] );
for ( int i = 2, stp = 1; i <= n; i <<= 1, stp <<= 1 ) {
for ( int j = 0; j < n; j += i ) {
rep ( k, 0, stp - 1 ) {
Complex w ( omega[n / stp * k].x, tp * omega[n / stp * k].y );
Complex ev ( A[j + k] ), ov ( w * A[j + k + stp] );
A[j + k] = ev + ov, A[j + k + stp] = ev - ov;
}
}
}
if ( !~tp ) rep ( i, 0, n - 1 ) A[i] = A[i] / n;
}inline void initFFT ( const int lg ) {
int n = 1 << lg;
rep ( i, 0, n - 1 ) rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << lg >> 1 );
for ( int i = 1; i < n; i <<= 1 ) {
rep ( k, 0, i - 1 ) {
omega[n / i * k] = Complex ( cos ( PI * k / i ), sin ( PI * k / i ) );
}
}
}inline void polyConv ( const int n, const int m, const int* A, const int* B, int* R ) {
rep ( i, 0, n - 1 ) P[i] = Complex ( A[i] & 0x7fff, A[i] >> 15 );
rep ( i, 0, m - 1 ) Q[i] = Complex ( B[i] & 0x7fff, B[i] >> 15 );
int lg = 0, len = 1;
for ( ; len < n + m - 1; len <<= 1, ++ lg );
rep ( i, n, len ) P[i] = Complex ();
rep ( i, m, len ) Q[i] = Complex ();
initFFT ( lg );
FFT ( len, P, 1 ), FFT ( len, Q, 1 );
rep ( i, 0, len - 1 ) {
Complex t ( P[( len - i ) % len].x, -P[( len - i ) % len].y );
C[i] = ( P[i] + t ) / 2, D[i] = Complex ( 0, 1 ) * ( t - P[i] ) / 2;
}
rep ( i, 0, len - 1 ) {
Complex t ( Q[( len - i ) % len].x, -Q[( len - i ) % len].y );
E[i] = ( Q[i] + t ) / 2, F[i] = Complex ( 0, 1 ) * ( t - Q[i] ) / 2;
}
rep ( i, 0, len - 1 ) {
Complex c ( C[i] ), d ( D[i] ), e ( E[i] ), f ( F[i] );
C[i] = c * e, D[i] = c * f + d * e, E[i] = d * f;
P[i] = C[i] + Complex ( 0, 1 ) * E[i];
}
FFT ( len, D, -1 ), FFT ( len, P, -1 );
rep ( i, 0, n + m - 2 ) {
int c = ( LL ( P[i].x + 0.5 ) % M + M ) % M,
d = ( LL ( D[i].x + 0.5 ) % M + M ) % M,
e = ( LL ( P[i].y + 0.5 ) % M + M ) % M;
R[i] = ( c + ( 1ll << 15 ) % M * d % M + ( 1ll << 30 ) % M * e % M ) % M;
}
}inline void polyDer ( const int n, const int* A, int* R ) {
rep ( i, 1, n - 1 ) R[i - 1] = mul ( i, A[i] );
R[n - 1] = 0;
}inline void polyInt ( const int n, const int* A, int* R ) {
per ( i, n - 1, 0 )R[i + 1] = mul ( inv[i + 1], A[i] );
R[0] = 0;
}inline void polyInv ( const int n, const int* A, int* R ) {
if ( n == 1 ) return void ( R[0] = mpow ( A[0], M - 2 ) );
static int tmp[MAXLEN + 5];
polyInv ( n + 1 >> 1, A, R ); polyConv ( n, n, A, R, tmp );
tmp[0] = sub ( 2, tmp[0] );
rep ( i, 1, n - 1 ) tmp[i] = sub ( 0, tmp[i] );
polyConv ( n, n, tmp, R, R );
}inline void polyLn ( const int n, const int* A, int* R ) {
static int tmp[2][MAXLEN + 5];
polyDer ( n, A, tmp[0] ), polyInv ( n, A, tmp[1] );
polyConv ( n, n, tmp[0], tmp[1], tmp[0] );
polyInt ( n, tmp[0], R );
}} // namesapce PolyOper.int main () {
n = rint (), M = rint ();
PolyOper::initInv ();
f[0] = 1;
rep ( i, 1, n ) f[i] = rint ();
PolyOper::polyLn ( n + 1, f, f );
rep ( i, 1, n ) f[i] = mul ( f[i], i );
int cnt = 0;
rep ( i, 1, n ) {
assert ( !f[i] || f[i] == i );
cnt += a[i] = !!f[i];
if ( a[i] ) rep ( j, 2, n / i ) f[i * j] = sub ( f[i * j], i );
}
wint ( cnt );
for ( int i = 1, flg = 0; i <= n; ++i ) {
if ( a[i] ) {
putchar ( flg ? ' ' : '\n' ), flg = 1;
wint ( i );
}
}
putchar ( '\n' );
return 0;
}
