【[SDOI2017]序列计数】

这道题非常的迷幻

首先我们要容斥

考虑记(dp[i][j])表示前(i)位(\%p=j)的方案数

(g[i][j])表示前(i)位只用合数(\%p=j)的方案数

于是可以考虑最暴力的(dp)是(O(nm^*p))的

但是并没必要

我们可以提前处理(1-m)这些数(\%p)的值，用这些值来转移就好了

也就是额外记一个(cnt[i])表示(\%p)为(i)的数量和(comp[i])表示(\%p)为(i)的合数的数量

于是就没必要枚举(1-m)转移了

复杂度变成(O(np^2))

实际上这个序列的生成与位置无关

也就是说我们可以倍增的处理

(dp[n][j])可以由(dp[n/2][j])来直接转移

于是递归做就好了

复杂度(O(p^2log n))

然而因为转移的形式是一个卷积的模样

所以是可以做到(O(plog plog n))的

但是出题人又不卡，就没写...

#include<bits/stdc++.h>
using namespace std;
#define rep( i, s, t ) for( register int i = s; i <= t; ++ i )
#define re register
#define int long long
int read() {
	char cc = getchar(); int cn = 0, flus = 1;
	while(cc < '0' || cc > '9') {  if( cc == '-' ) flus = -flus;  cc = getchar();  }
	while(cc >= '0' && cc <= '9')  cn = cn * 10 + cc - '0', cc = getchar();
	return cn * flus;
}
const int P = 20170408 ; 
const int N = 2000000 + 5 ;
const int M = 100 + 5 ; 
const int R = 2 * 1e7 + 5 ; 
int n, m, p, dp[2][M], g[2][M], cnt[M], comp[M], Ed ; 
int pr[N], top;
bool isp[R];
void Init() {
	isp[1] = 1 ; 
	for( int i = 2; i <= m; ++ i ) {
		if( !isp[i] ) pr[++ top] = i, isp[i] = 0 ;
		for(int j = 1; j <= top; ++ j ) {
			if ( pr[j] * i > m ) break ;
			isp[pr[j] * i] = 1 ;
			if( i % pr[j] == 0 )  break;
		}
	}
}
void F( int x,  int w ) {
	if( x == 1 ) {
		rep( i, 0, ( p - 1 ) ) dp[w][i] = cnt[i], g[w][i] = comp[i] ; 
		return ; 
	}
	F( x / 2, w ^ 1 ) ;
	rep( j, 0, Ed ) dp[w][j] = 0, g[w][j] = 0 ;  
	rep( j, 0, Ed ) rep( k, 0, Ed ) 
	dp[w][j] = ( dp[w ^ 1][k] * dp[w ^ 1][(p + j - k) % p] % P + dp[w][j] ) % P ;
	rep( j, 0, Ed ) rep( k, 0, Ed ) 
	g[w][j] = ( g[w ^ 1][k] * g[w ^ 1][(p + j - k) % p] % P + g[w][j] ) % P ;
	if( x & 1 ) {
		rep( j, 0, Ed ) dp[w ^ 1][j] = dp[w][j], g[w ^ 1][j] = g[w][j] ;
		rep( j, 0, Ed ) dp[w][j] = 0, g[w][j] = 0 ;  
		rep( j, 0, Ed ) rep( k, 0, Ed )
		dp[w][j] = ( dp[w ^ 1][k] * cnt[(p + j - k) % p] % P + dp[w][j] ) % P ;
		rep( j, 0, Ed ) rep( k, 0, Ed ) 
		g[w][j] = ( g[w ^ 1][k] * comp[(p + j - k) % p] % P + g[w][j] ) % P ;
	}
}
signed main()
{
	n = read(), m = read(), p = read() ; isp[1] = 1, Init() ; 
	rep( i, 1, m ) cnt[i % p] ++, comp[i % p] += isp[i] ;
	 Ed = p - 1 ; F( n, 0 ) ;
	printf("%d
", ( dp[0][0] - g[0][0] + P ) % P ) ; 
	return 0;
}