結果

問題 No.206 数の積集合を求めるクエリ
コンテスト
ユーザー Ryuhei Mori
提出日時 2018-07-23 09:02:34
言語 C
(gcc 13.3.0)
結果
AC  
実行時間 18 ms / 7,000 ms
コード長 4,875 bytes
コンパイル時間 521 ms
コンパイル使用メモリ 38,272 KB
実行使用メモリ 6,400 KB
最終ジャッジ日時 2024-12-26 05:35:27
合計ジャッジ時間 2,335 ms
ジャッジサーバーID
（参考情報） 		judge5 / judge2
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
sample AC * 3
other AC * 28
権限があれば一括ダウンロードができます

ソースコード

diff # 

#pragma GCC target ("avx")
#include <unistd.h>
#include <complex.h>

#define PI2 6.28318530717958647692528676655900577L

char ibuf[1200100];
char *ibufe = ibuf-1;

void readall(){
  int k, t = 0;
  while((k=read(STDIN_FILENO, ibuf+t, sizeof(ibuf)-t))>0) t += k;
}

int read_uint(){
  int x=0;
  while(*(++ibufe) <'0');
  do {
    x *= 10;
    x += *ibufe-'0';
  } while(*(++ibufe) >='0');

  return x;
}

char buf[700000];
char *bufe = buf;

void write_uintln(int x){
  int i;
  static char tmp[13];
  if(x==0){
    *bufe++ = '0';
    *bufe++ = '\n';
    return;
  }

  for(i=0; x; i++){
    tmp[i] = '0' + x % 10;
    x /= 10;
  }
  for(i--; i >= 0; i--){
    *bufe++ = tmp[i];
  }
  *bufe++ = '\n';
}

void writeall(){
  int k, t = 0;
  while((k=write(STDOUT_FILENO, buf+t, bufe-buf-t))>0) t += k;
}


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

typedef double complex cmplx;

cmplx C[1<<17];
cmplx w[1<<16];
long double complex v[16];

void fft(int k, cmplx *A, const cmplx *w) __attribute__((optimize("fast-math")));
void fft(int k, cmplx *A, const cmplx *w){
  const int m = 1 << k;
  int u = 1;
  int v = m/4;
  int i, j;
  if(k&1){
    for(j=0; j<m/2; j++){
      cmplx Ajv = A[j+(m/2)];
      A[j+(m/2)] = A[j] - Ajv;
      A[j] += Ajv;
    }
    u <<= 1;
    v >>= 1;
  }
  for(i=k&~1;i>0;i-=2){
    int jh;
    for(jh=0;jh<u;jh++){
      cmplx wj = w[jh<<1];
      cmplx wj2 = w[jh];
      cmplx wj3 = wj2 * wj;
//      cmplx wj3 = creal(wj) * (2 * creal(wj2) - 1) + (2 * creal(wj) * cimag(wj2) - cimag(wj)) * I;
      int je;
      for(j = jh << i, je = j+v;j<je; j++){
        cmplx tmp0 = A[j];
        cmplx tmp1 = wj * A[j+v];
        cmplx tmp2 = wj2 * A[j+2*v];
        cmplx tmp3 = wj3 * A[j+3*v];

        cmplx ttmp0 = tmp0 + tmp2;
        cmplx ttmp2 = tmp0 - tmp2;
        cmplx ttmp1 = tmp1 + tmp3;
        cmplx ttmp3 = -I * (tmp1 - tmp3);

        A[j] = ttmp0 + ttmp1;
        A[j+v] = ttmp0 - ttmp1;
        A[j+2*v] = ttmp2 + ttmp3;
        A[j+3*v] = ttmp2 - ttmp3;
      }
    }
    u <<= 2;
    v >>= 2;
  }
}


void ifft(int k, cmplx *A, const cmplx *w) __attribute__((optimize("fast-math")));
void ifft(int k, cmplx *A, const cmplx *w){
  const int m = 1 << k;
  int u = m/4;
  int v = 1;
  int i, j;
  for(i=2;i<=k;i+=2){
    int jh;
    for(jh=0;jh<u;jh++){
      cmplx wj = conj(w[jh<<1]);
      cmplx wj2 = conj(w[jh]);
      cmplx wj3 = wj2 * wj;
//      cmplx wj3 = creal(wj) * (2 * creal(wj2) - 1) + (2 * creal(wj) * cimag(wj2) - cimag(wj)) * I;
      int je;
      for(j = jh << i, je = j+v;j<je; j++){
        cmplx tmp0 = A[j];
        cmplx tmp1 = A[j+v];
        cmplx tmp2 = A[j+2*v];
        cmplx tmp3 = A[j+3*v];

        cmplx ttmp0 = tmp0 + tmp1;
        cmplx ttmp1 = tmp0 - tmp1;
        cmplx ttmp2 = tmp2 + tmp3;
        cmplx ttmp3 = I * (tmp2 - tmp3);

        A[j] = ttmp0 + ttmp2;
        A[j+v] = wj * (ttmp1 + ttmp3);
        A[j+2*v] = wj2 * (ttmp0 - ttmp2);
        A[j+3*v] = wj3 * (ttmp1 - ttmp3);
      }
    }
    u >>= 2;
    v <<= 2;
  }
  if(k&1){
    for(j = 0;j<m/2; j++){
      cmplx Ajv = A[j+(m/2)];
      A[j+(m/2)] = A[j] - Ajv;
      A[j] += Ajv;
    }
  }
}

void genw(int i, int b, long double complex z){
  if(b == -1){
    w[i] = z;
  }
  else {
    genw(i, b-1, z);
    genw(i|(1<<b), b-1, z*v[b]);
  }
}

void convolver(int k, cmplx *A, const cmplx *w) __attribute__((optimize("fast-math")));
void convolver(int k, cmplx *A, const cmplx *w){
  int i, y;
  const int m = 1 << k;

  fft(k, A, w);
  A[0] = 4 * creal(A[0]) * cimag(A[0]) * I;
  A[1] = 4 * creal(A[1]) * cimag(A[1]) * I;
  i = 2;
  for(y = 2; y < m; y <<= 1){
    for(; i < 2*y; i+=2){
      int j = i^(y-1);
      A[i] = (A[i] + conj(A[j]))*(A[i] - conj(A[j]));
      A[j] = -conj(A[i]);
    }
  }

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

  ifft(k-1, A, w);
}


int main() __attribute__((optimize("fast-math")));
int main(){
  int l, mm, q;
  int i, j;
  const long double arg = -PI2/m;
//  const int B = 100352;
  const int B = 100001;
//  const int B = 1<<17;

  for(i=0, j=m/4; j; i++, j>>=1){
    v[i] = cexpl(I * (arg * j));
  }
  genw(0, k-2, 1);

  readall();
  l = read_uint();
  mm = read_uint();
  read_uint();


  for(i=0;i<l;i++){
    int a;
    a = read_uint();
    if(a&1) C[a/2] += B;
    else    C[a/2] += 1;
  }
  for(i=0;i<mm;i++){
    int a;
    a = read_uint();
    if(a&1) C[((1<<18)-a)/2] += B*I;
    else    C[((1<<18)-a)/2] += I;
  }

  q = read_uint();

  convolver(k, C, w);

  double *CC = (double *)C;
  unsigned long long int tmp =(unsigned long long)(CC[0]+0.5)/B/B;
  for(i=1;i<=q/2;i++){
    tmp += CC[i]+0.5;
    write_uintln(tmp%B);
    tmp /= B;
    write_uintln(tmp%B);
    tmp /= B;
  }
  if(q&1){
    tmp += CC[i]+0.5;
    write_uintln(tmp%B);
  }
  writeall();

  return 0;
}
0