编织幻梦

今天我在翻DP的时候，意外的发现了这个有趣的东西:-D。

其实很简单……只是居然可以用Dp。而且能做的事情很多。

#include <iostream>
#include <string.h>
#include <string>
#include <math.h>
#include <stdio.h>
#include <algorithm>
using namespace std;
const int maxh = 10,maxn = 1<<maxh,P = 12289;
const int g=10302,ig = 8974,ninv = 12277;
typedef long long ll;
struct fft{
    ll omg[maxn],iomg[maxn],inv[maxn];
    void initialize(){
        for(int i=0;i<maxn;i++) inv[i] = getinv(i);
        omg[0] = 1,iomg[0] = 1;
        for(int i=1;i<maxn;i++)
            omg[i] = omg[i-1]*g%P,iomg[i] = iomg[i-1]*ig%P;
    }
    ll getinv(ll x){
        ll ans = 0;
        for(int i=0;i<maxh;i++)
            ans |= (x>>i&1)<<maxh-i-1;
        return ans;
    }
    ll f[2][maxn];
    void ft(ll *o,ll *cof,ll *output){
        for(int i=0;i<maxn;i++) f[0][i] = cof[inv[i]];
        for(int h=1;h<=maxh;h++)
            for(int i=0;i<maxn;i++)
                f[h&1][i] = (f[h&1^1][~(1<<h-1)&i]+o[i<<maxh-h&maxn-1]*f[h&1^1][i|1<<h-1]%P)%P;
        for(int i=0;i<maxn;i++) output[i] = f[maxh&1][i];
    }void dft(ll *cof,ll *output){ft(omg,cof,output);}
    void idft(ll *cof,ll *output){ft(iomg,cof,output);for(int i=0;i<maxn;i++) output[i] = output[i]*ninv%P;}
}F;
ll cof1[maxn],cof2[maxn];
char input[maxn],input2[maxn];
int n;
void Init(){
    int len;
    gets(input2);len = strlen(input2);
    for(int i=0;i<len;i++) input[i] = input2[len-i-1];
    for(int i=0;i<len;i++) cof1[i] = input[i]-'0';
    gets(input2);len = strlen(input2);
    for(int i=0;i<len;i++) input[i] = input2[len-i-1];
    for(int i=0;i<len;i++) cof2[i] = input[i]-'0';
}
void solve(){
    F.initialize();
    F.dft(cof1,cof1);F.dft(cof2,cof2);
    for(int i=0;i<maxn;i++) cof1[i] = cof1[i]*cof2[i]%P;
    F.idft(cof1,cof1);
    int jw = 0;
    for(int i=0;i<maxn;i++){
        cof1[i] += jw%10,jw/=10;
        jw+=cof1[i]/10,cof1[i]%=10;
    }
    bool flag = false;
    for(int i=maxn-1;i>=0;i--){
        if(cof1[i]!=0) flag = true;
        if(flag) printf("%c",cof1[i]+'0');
    }
}
int main(){
    freopen("test.in","r",stdin);
    freopen("test.out","w",stdout);
 
    Init();
    solve();
 
    fclose(stdin);
    fclose(stdout);
    return 0;
}

考虑n! mod p

n! mod p

 = (1 * 2 * 3 * 4 * 5 * 6 * 7 * 8 * 9 * ... n) mod p

 = (1 * 2 * 3 * ... * p * p+1 * p+2 * ... * 2p * 2p+1 * 2p+2 * 2p+3 * ...) mod p

 = (1 * 2 * 3 * ... * p-1 * p+1 * p+2 * ... * 2p-1 * 2p+1 * ... *p^(n/p-1)* (2*3*...*n/p) ) mod p

其中

1 * 2 * 3 * ... p-1 mod p = p+1 * p+2 * p+3 * p+4 ... * 2p-1 mod p = ... = (n%p!) mod p

相当于每p段加上p，所以全部约掉，就是n%p! mod p。

所以，只需要计算C(n%p,m%p) * Lucas(n/p,m/p)

简化： 

C(N,M) mod p=C(Np0,Mp0)*C(Np1,Mp1) * C(Npx,Mpx)

NpX,MpX表示N和M在p进制下的第 X 位。

----------------------------------------------------------------------------------------------------------------------------------------

扩展，求C(N,M) MOD P^c