[数学 FFT] Codechef July Challenge 2017 #APRPS Irrational Root

跟 51Nod 1356 代数数的次数 是一样的

不过这里都是质数 也就是就是 2n

关键是输方案 这个不一定有二次剩余

感谢sxt

一个一个数加进答案 转化成 已知F(x)，求F(x+a√)和F(x−a√)的系数

这个推一下就是一个FFT

拷了myy的板子 自己的太慢了QAQ

#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cstring>
#define cl(x) memset(x,0,sizeof(x))
using namespace std;
typedef long long ll;

inline void write(int x){
if (x>=10) write(x/10);
putchar('0'+x%10);
}

const int P=1e9+7;

inline ll Pow(ll a,int b){
ll ret=1;
for (;b;b>>=1,a=a*a%P)
if (b&1)
ret=ret*a%P;
return ret;
}

const int N=262144;

namespace MYY{
#define REP(i, a, b) for (int i = (a), _end_ = (b); i < _end_; ++i)
#define debug(...) fprintf(stderr, __VA_ARGS__)
#define mp make_pair
#define x first
#define y second
#define pb push_back
#define SZ(x) (int((x).size()))
#define ALL(x) (x).begin(), (x).end()

template<typename T> inline bool chkmin(T &a, const T &b) { return a > b ? a = b, 1 : 0; }
template<typename T> inline bool chkmax(T &a, const T &b) { return a < b ? a = b, 1 : 0; }

typedef long long LL;

const int oo = 0x3f3f3f3f;

const int Mod = 1e9 + 7;

const int max0 = 262144;

struct comp
{
double x, y;

comp(): x(0), y(0) { }
comp(const double &_x, const double &_y): x(_x), y(_y) { }

};

inline comp operator+(const comp &a, const comp &b) { return comp(a.x + b.x, a.y + b.y); }
inline comp operator-(const comp &a, const comp &b) { return comp(a.x - b.x, a.y - b.y); }
inline comp operator*(const comp &a, const comp &b) { return comp(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x); }
inline comp conj(const comp &a) { return comp(a.x, -a.y); }

const double PI = acos(-1);

int N, L;

comp w[max0 + 5];
int bitrev[max0 + 5];

void fft(comp *a, const int &n)
{
REP(i, 0, n) if (i < bitrev[i]) swap(a[i], a[bitrev[i]]);
for (int i = 2, lyc = n >> 1; i <= n; i <<= 1, lyc >>= 1)
for (int j = 0; j < n; j += i)
{
comp *l = a + j, *r = a + j + (i >> 1), *p = w;
REP(k, 0, i >> 1)
{
comp tmp = *r * *p;
*r = *l - tmp, *l = *l + tmp;
++l, ++r, p += lyc;
}
}
}

inline void fft_prepare()
{
REP(i, 0, N) bitrev[i] = bitrev[i >> 1] >> 1 | ((i & 1) << (L - 1));
REP(i, 0, N) w[i] = comp(cos(2 * PI * i / N), sin(2 * PI * i / N));
}

inline void conv(int *x, int *y, int *z)
{
REP(i, 0, N) (x[i] += Mod) %= Mod, (y[i] += Mod) %= Mod;
static comp a[max0 + 5], b[max0 + 5];
static comp dfta[max0 + 5], dftb[max0 + 5], dftc[max0 + 5], dftd[max0 + 5];

REP(i, 0, N) a[i] = comp(x[i] & 32767, x[i] >> 15);
REP(i, 0, N) b[i] = comp(y[i] & 32767, y[i] >> 15);
fft(a, N), fft(b, N);
REP(i, 0, N)
{
int j = (N - i) & (N - 1);
static comp da, db, dc, dd;
da = (a[i] + conj(a[j])) * comp(0.5, 0);
db = (a[i] - conj(a[j])) * comp(0, -0.5);
dc = (b[i] + conj(b[j])) * comp(0.5, 0);
dd = (b[i] - conj(b[j])) * comp(0, -0.5);
dfta[j] = da * dc;
dftb[j] = da * dd;
dftc[j] = db * dc;
dftd[j] = db * dd;
}
REP(i, 0, N) a[i] = dfta[i] + dftb[i] * comp(0, 1);
REP(i, 0, N) b[i] = dftc[i] + dftd[i] * comp(0, 1);
fft(a, N), fft(b, N);
REP(i, 0, N)
{
int da = (LL)(a[i].x / N + 0.5) % Mod;
int db = (LL)(a[i].y / N + 0.5) % Mod;
int dc = (LL)(b[i].x / N + 0.5) % Mod;
int dd = (LL)(b[i].y / N + 0.5) % Mod;
z[i] = (da + ((LL)(db + dc) << 15) + ((LL)dd << 30)) % Mod;
}
}

}

int A
,B
,C
,D
,E
,F
;

#define read(x) scanf("%d",&(x))

int m,num
;
ll fac
,inv
;

int n;
ll f
;

int main(){
int T;
freopen("t.in","r",stdin);
freopen("t.out","w",stdout);
read(T);
fac[0]=1; for (int i=1;i<N;i++) fac[i]=fac[i-1]*i%P;
inv[1]=1; for (int i=2;i<N;i++) inv[i]=(ll)(P-P/i)*inv[P%i]%P;
inv[0]=1; for (int i=1;i<N;i++) inv[i]=inv[i]*inv[i-1]%P;
while (T--){
read(m); for (int i=1;i<=m;i++) read(num[i]);
cl(f); cl(C); cl(D); cl(E); cl(F);
f[0]=0; f[1]=1; n=2;
for (int t=1;t<=m;t++){
MYY::L = 0;
for ( ; (1 << MYY::L) < (n<<2); ++MYY::L);
MYY::N = 1 << MYY::L;
MYY::fft_prepare();

for (int i=0;i<=n;i++){
A[i]=f[n-i]*fac[n-i]%P;
if (~i&1) B[i]=Pow(num[t],i>>1)*inv[i]%P;
}
MYY::conv(A,B,C);
for (int i=0;i<(n<<2);i++) A[i]=B[i]=0;
for (int i=0;i<=n;i++)
E[i]=C[n-i]*inv[i]%P;

for (int i=1;i<=n;i++){
A[i]=f[n-i]*fac[n-i]%P;
if (i&1) B[i]=Pow(num[t],i>>1)*inv[i]%P;
}
MYY::conv(A,B,D);
for (int i=0;i<(n<<2);i++) A[i]=B[i]=0;
for (int i=0;i<=n;i++)
F[i]=D[n-i]*inv[i]%P;

MYY::conv(E,E,E);
MYY::conv(F,F,F);

for (int i=0;i<=(n<<1);i++)
f[i]=(E[i]+P-(ll)F[i]*num[t]%P)%P;
n<<=1;
}
write(n>>1); putchar('\n');
for (int i=0;i<=(n>>1);i++)
write(f[i]),putchar(' ');
putchar('\n');
}
return 0;
}