(mathcal{Description})
Link.
给定 ({f_1,f_2,cdots,f_n})，素数 (p)。求字典序最小的 ({a_1,a_2,cdots,a_n})，满足对于所有 (iin[1,n])，(a_iin{0,1}) 并且
[sum_{{k_{1..n}}}[(forall j)left((a_j=0land k_j=0)lor(a_j
ot=0land k_jge0
ight)]left[ sum_{j=1}^njk_j=i
ight]equiv f_ipmod p
]
(n<2^{18})。
(mathcal{Solution})
对于某个 (a_i=1)，其 OGF 为 (frac{1}{1-x^i})，所有 OGF 之积的 (1sim n) 次项系数在(mod p) 意义下就是 (f_{1..n})。记 (F(x)=sum_{i=0}^{+infty}f_ix^i)，默认运算在(mod p) 意义下进行，那么：
[F(x)=prod_{i=1}^nleft( frac{1}{1-x^i}
ight)^{a_i}
]
两边取 (-ln)：
[Rightarrow~~~~-ln F(x)=sum_{i=1}^na_iln(1-x^i)
]
Tayler 展开右式 (ln)：
[Rightarrow~~~~-ln F(x)=sum_{i=1}^na_isum_{j=1}^{+infty}-frac{x^{ij}}j
]
左右取负，枚举 (T=ij)：
[Rightarrow~~~~ln F(x)=sum_{T=1}^{+infty}x^Tsum_{i|T}a_ifrac{i}T
]
对已知的 (F(x)) 求 (ln)，就能取得 (Tin[1,n]) 内的所有 (frac{1}Tsum_{i|T}a_ii)，从 (T=1) 起刷表求解，每次枚举倍数消除贡献，可以 (mathcal O(nln n)) 解出所有 (a_{1..n})（没错，由于 (ln F(x)) 模意义下的多项式系数唯一，({a_n}) 唯一确定）。
综上，复杂度为多项式 (ln) 的 (mathcal O(nlog 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 ? ' ' : '
' ), flg = 1;
wint ( i );
}
}
putchar ( '
' );
return 0;
}