[[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^2\log n)\
그런데 전이의 형식이 권적의 모양이기 때문에.
따라서\(O (p\log p\log n)\를 수행할 수 있습니다.
그런데 출제자가 끊기지 않아서 안 썼는데...

#include
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;
}