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

跟 51Nod 1356 代数数的次数 是一样的
不过这里都是质数 也就是就是 $2^n$
关键是输方案 这个不一定有二次剩余

感谢sxt
一个一个数加进答案 转化成 已知$F(x)$，求$F(x+\sqrt a)$和$F(x-\sqrt 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[N],B[N],C[N],D[N],E[N],F[N];

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

int m,num[N];
ll fac[N],inv[N];

int n;
ll f[N];

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