Submission #1137844

---

Source Code Expand

#include <stdio.h>
#include <complex.h>

#define PI2 6.28318530717958647692

const int k = 18;
const int m = (1 << 18);

complex double C[1<<18];
complex double w[1<<17];

void read_int(int *x){
  char c=0;
  while(c<'0'||c>'9') c = getchar_unlocked();

  *x = c-'0';
  c = getchar_unlocked();

  while(c>='0'&&c<='9') {
    *x *= 10;
    *x += c-'0';
    c = getchar_unlocked();
  }
}

void write_int(int x){
  static char buf[12];
  int i;
  if(x==0){
    puts("0");
    return;
  }

  for(i=10; x; i--){
    buf[i] = '0' + x % 10;
    x /= 10;
  }
  puts(buf+1+i);
}

unsigned minus(unsigned x) {
  int y = 1 << (31-__builtin_clz(x));
  return (~x&(y-1))^y;
}


void fft(complex double *A, int k){
  int m = 1 << k;
  int u = 1;
  int v = m/2;
  for(int i=k;i>0;i--){
    for(int jh=0;jh<u;jh++){
      complex double wj = w[jh];
      for(int j = jh << i, je = j|v;j<je; j++){
        complex double Ajv = wj * A[j|v];
        A[j|v] = A[j] - Ajv;
        A[j] += Ajv;
      }
    }
    u <<= 1;
    v >>= 1;
  }
}

void ifft(complex double *A, int k){
  int m = 1 << k;
  int u = m/2;
  int v = 1;
  for(int i=1;i<=k;i++){
    for(int jh=0;jh<u;jh++){
      complex double wj = conj(w[jh]);
      for(int j = jh << i, je = j|v;j<je; j++){
        complex double Ajv = A[j] - A[j|v];
        A[j] += A[j|v];
        A[j|v] = wj * Ajv;
      }
    }
    u >>= 1;
    v <<= 1;
  }
}

void genw(int i, int b, complex double z){
  if(b == 0){
    w[i] = z;
  }
  else {
    genw(i, b>>1, z);
    genw(i|b, b>>1, z*w[b]);
  }
}

int main(){
  int n;
  const double arg = -PI2/m;

  for(int i=1, j=m/4; j; i<<=1, j>>=1){
    w[i] = cexp(I * (arg * j));
  }
  genw(0, m/4, 1);

  scanf("%d", &n);
  for(int i=1;i<=n;i++){
    int a, b;
    read_int(&a);
    read_int(&b);
    C[i] = a + b*I;
  }
//  k = 33 - __builtin_clz(n);
//  m = 1 << k;

  fft(C, k);

  C[0] = 4 * creal(C[0]) * cimag(C[0]) * I;
  C[1] = 4 * creal(C[1]) * cimag(C[1]) * I;

  for(int i = 2; i < m; i+=2){
    int j = minus(i);
    C[i] = (C[i] + conj(C[j]))*(C[i] - conj(C[j]));
    C[j] = -conj(C[i]);
  }

  for(int i = 0; i < m; i+=2){
    C[i/2] = C[i]+C[i^1] + (C[i]-C[i^1])*w[i/2]*I;
  }

  ifft(C, k-1);

  for(int i=1;i<=n;i++){
    write_int(-creal(C[i])/(4*m)+0.5);
    write_int(cimag(C[i])/(4*m)+0.5);
  }
}

Submission Info

Judge Result