天天看點

【HDU】5197 DZY Loves Orzing 【FFT啟發式合并】

傳送門:【HDU】5197 DZY Loves Orzing

題目分析:

首先申明,我不會 dp 方程= =……這個東西給隊友找出來了,然後我就是套這個方程做題的Qrz……

對于這題,因為 n2 個數互不相同,是以每一列都可以單獨考慮。設 dpni 表示長度為 n 的排列,能恰好看見i個人的方案數,根據隊友的發現, dpni 就等于 |sni| ,其中 sni 是第一類 Stirling 數。(不要問我為什麼是這個,我不知道……)

由于對于第一類 Stirling 數我們有:

Sn=x(x−1)(x−2)⋯(x−n+1)

=∑ni=1(x+1−i)

用多項式表示則為:

Sn=∑ni=1sinxi

則我們的 dp 方程可以表示為:

dpn=∑ni=1|sin|xi

由 sin=(−1)n+i|sin| 得:

dpn=∑ni=1(−1)n+isinxi

是以我們隻要求的 Stirling 數就可以得到我們對應的dp值。

然而我們怎麼去求 Stirling 數呢?看到多個多項式相乘,很容易我們可以想到用 FFT 去優化,但是如果是簡單的從左往右做 n−1 次 FFT 的話,顯然複雜度爆表。

那我們該怎麼辦呢?

這時候一個自然的想法是對 FFT 之間的乘法分治!因為我們用高階多項式和低階多項式相乘,會導緻很多時間上的浪費,這時候我們就可以考慮用類似于啟發式合并的思想去進行多項式乘法:每次取出階數最小的兩個多項式進行 FFT ,這樣我們複雜度就保證了。(這個過程可以用優先隊列實作)

最後的答案即:

ans=(n2)!(n!)2∏ni=1dpain

其中 ai 為輸入資料。這裡還有個問題,就是 (n2)! 太大了怎麼辦?這個好說,打個表呗(對于給定的模數P,注意到 n2≥P 時 ans=0 )

時間複雜度即:

n2⋅2log2+n4⋅4log4+⋯+2⋅n2logn2+nlogn

=∑logni=1nlogi

=O(nlog2n)

然後由于題目給的模數P=999948289恰為費馬素數,是以我們可以将 FFT 換成 NTT 來保證答案精度。(算法中我用的原根為 g=13 )

時間複雜度: O(nlog2n)

空間複雜度: O(→_→ 看臉 )

my  code:

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <queue>
#include <algorithm>
using namespace std ;

typedef long long LL ;

#define clr( a , x ) memset ( a , x , sizeof a )
#define cpy( a , x ) memcpy ( a , x , sizeof a )

const int mod =  ;
const int MAXN =  ;
const int g =  ;

struct Node {
    int idx , n ;
    Node () {}
    Node ( int idx , int n ) : idx ( idx ) , n ( n ) {}
    bool operator < ( const Node& a ) const {
        if ( n != a.n ) return n > a.n ;
        return idx > a.idx ;
    }
} ;

vector < int > G[MAXN] ;
int x1[MAXN << ] , x2[MAXN << ] ;
int A[MAXN] ;
int f[MAXN] ;
int invf[MAXN] ;
int ff[MAXN] = {  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  } ;
int m , root ;

int power ( int a , int b ) {
    LL res =  , tmp = a ;
    while ( b ) {
        if ( b &  ) res = res * tmp % mod ;
        tmp = tmp * tmp % mod ;
        b >>=  ;
    }
    return ( int ) res ;
}

void NTT ( int y[] , int n , int rev ) {
    for ( int i =  , j , k , t ; i < n ; ++ i ) {
        for ( j =  , k = n >>  , t = i ; k ; k >>=  , t >>=  ) j = j <<  | t &  ;
        if ( i < j ) swap ( y[i] , y[j] ) ;
    }
    for ( int s =  , ds =  ; s <= n ; ds = s , s <<=  ) {
        int wn = power ( g , ( mod -  ) / s ) ;
        if ( rev <  ) wn = power ( wn , mod -  ) ;
        for ( int k =  , w =  ; k < ds ; ++ k , w = ( LL ) w * wn % mod ) {
            for ( int i = k , t ; i < n ; i += s ) {
                y[i + ds] = ( y[i] - ( t = ( LL ) w * y[i + ds] % mod ) + mod ) % mod ;
                y[i] = ( y[i] + t ) % mod ;
            }
        }
    }
    if ( rev <  ) {
        int inv = power ( n , mod -  ) ;
        for ( int i =  ; i < n ; ++ i ) y[i] = ( LL ) y[i] * inv % mod ;
    }
}

void solve () {
    priority_queue < Node > q ;
    for ( int i =  ; i <= m ; ++ i ) {
        scanf ( "%d" , &A[i] ) ;
        G[i].clear () ;
    }
    if ( ( LL ) m * m >= mod ) {
        printf ( "0\n" ) ;
        return ;
    }
    for ( int i =  ; i <= m ; ++ i ) {
        G[i].push_back ( (  - i + mod ) % mod ) ;
        G[i].push_back (  ) ;
        q.push ( Node ( i ,  ) ) ;
    }
    while ( !q.empty () ) {
        int x = q.top ().idx ;
        q.pop () ;
        if ( q.empty () ) {
            root = x ;
            break ;
        }
        int y = q.top ().idx ;
        q.pop () ;
        int n1 = G[x].size () , n2 = G[y].size () , nn = n1 + n2 -  ;
        int n =  ;
        while ( n < nn ) n <<=  ;
        for ( int i =  ; i < n ; ++ i ) x1[i] = i < n1 ? G[x][i] :  ;
        for ( int i =  ; i < n ; ++ i ) x2[i] = i < n2 ? G[y][i] :  ;
        NTT ( x1 , n ,  ) ;
        NTT ( x2 , n ,  ) ;
        for ( int i =  ; i < n ; ++ i ) x1[i] = ( LL ) x1[i] * x2[i] % mod ;
        NTT ( x1 , n , - ) ;
        G[x].clear () ;
        for ( int i =  ; i < nn ; ++ i ) G[x].push_back ( x1[i] ) ;
        q.push ( Node ( x , nn ) ) ;
    }
    int ans = power ( invf[m] , m ) ;
    int mm = m * m , x ;
    for ( x =  ; x +  <= mm ; x +=  ) ans = ( LL ) ans * ff[x / ] % mod ;
    for ( ++ x ; x <= mm ; ++ x ) ans = ( LL ) ans * x % mod ;
    for ( int i =  ; i <= m ; ++ i ) {
        if ( ( m + A[i] ) %  ==  ) ans = ( LL ) ans * G[root][A[i]] % mod ;
        else ans = ( LL ) ans * ( mod - G[root][A[i]] ) % mod ;
    }
    printf ( "%d\n" , ans ) ;
}

int main () {
    f[] =  ;
    for ( int i =  ; i < MAXN ; ++ i ) {
        f[i] = ( LL ) i * f[i - ] % mod ;
        invf[i] = power ( f[i] , mod -  ) ;
    }
    while ( ~scanf ( "%d" , &m ) ) solve () ;
    return  ;
}