結果

問題 No.215 素数サイコロと合成数サイコロ (3-Hard)
ユーザー Ryuhei MoriRyuhei Mori
提出日時 2018-08-13 11:55:30
言語 C90
(gcc 12.3.0)
結果
AC  
実行時間 92 ms / 4,000 ms
コード長 6,267 bytes
コンパイル時間 404 ms
コンパイル使用メモリ 36,032 KB
実行使用メモリ 5,376 KB
最終ジャッジ日時 2024-06-26 03:15:02
合計ジャッジ時間 1,107 ms
ジャッジサーバーID
（参考情報） 		judge4 / judge5
このコードへのチャレンジ
（要ログイン）
ファイルパターン 結果
other AC * 2
権限があれば一括ダウンロードができます
コンパイルメッセージ
main.c: In function ‘main’:
main.c:186:3: warning: ignoring return value of ‘scanf’ declared with attribute ‘warn_unused_result’ [-Wunused-result]
  186 |   scanf("%lld%d%d", &n, &p, &c);
      |   ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~

ソースコード

diff # 

#pragma GCC optimize ("fast-math")
#include <stdio.h>
#include <complex.h>
#include <math.h>

#define PI2 6.28318530717958647692

typedef double complex cmplx;


#define K 14
#define M (1<<K)
#define B 83666LL

typedef unsigned Zp;
const Zp mod = 1000000007;

Zp add(Zp lhs, Zp rhs){
  Zp r = lhs + rhs;
  if(r >= mod) r -= mod;
  return r;
/*
  rhs = mod-rhs;
  if(lhs >= rhs) return lhs - rhs;
  else return mod + lhs - rhs;
*/
}

Zp sub(Zp lhs, Zp rhs){
  if(lhs >= rhs) return lhs - rhs;
  else return mod + lhs - rhs;
}

Zp mul(Zp lhs, Zp rhs){
  return ((unsigned long long) lhs * rhs) % mod;
}

Zp pos(int x){
  if(x < 0) x += mod;
  return x;
}

cmplx PL[M];
cmplx PH[M];
cmplx QL[M];
cmplx QH[M];

cmplx w[M>>1];

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

void init(int k, cmplx *w){
  int i, j;
  const int m = 1<<k;
  const double arg = -PI2/m;
  for(i=1, j=m/4; j; i<<=1, j>>=1){
    w[i] = cexp(I * (arg * j));
  }
  genw(0, m/4, 1, w);
}

void fft(cmplx *A, int k){
  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;
      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(cmplx *A, int k){
  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;
      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;
    }
  }
  for(j=0;j<m;j++) A[j] /= m;
}

void fft_real(cmplx *AL, cmplx *AH, int k){
  int i, y;
  fft(AL, k);
  AH[0] = cimag(AL[0]);
  AL[0] = creal(AL[0]);
  AH[1] = cimag(AL[1]);
  AL[1] = creal(AL[1]);
  for(i=y=2;y<(1<<k);y<<=1){
    for(;i<2*y;i+=2){
      int j = i^(y-1);
      AH[i] = -I*(AL[i] - conj(AL[j]))/2;
      AL[i] = (AL[i] + conj(AL[j]))/2;
      AL[j] = conj(AL[i]);
      AH[j] = conj(AH[i]);
    }
  }
}


const int P[] = { 2, 3, 5, 7, 11, 13 };
const int C[] = { 4, 6, 8, 9, 10, 12 };
const int S = sizeof(P)/sizeof(*P);

int dpP[7][13*300+1];
int dpC[7][12*300+1];
/*
dp[i][k][j] = dp[i][k-1][j-P[i]] + dp[i-1][k][j]
*/

int main(){
  int p, i, j, k, c;
  long long int n;
  scanf("%lld%d%d", &n, &p, &c);
  
  for(i=1;i<=S;i++) dpP[i][0] = 1;
  for(k=0;k<p;k++){
    for(i=1;i<=S;i++){
      for(j=P[S-1]*300;j>=P[i-1];j--){
        dpP[i][j] = add(dpP[i-1][j], dpP[i][j-P[i-1]]);
      }
      for(;j>=0;j--){
        dpP[i][j] = dpP[i-1][j];
      }
    }
  }
  for(i=1;i<=S;i++) dpC[i][0] = 1;
  for(k=0;k<c;k++){
    for(i=1;i<=S;i++){
      for(j=C[S-1]*300;j>=C[i-1];j--){
        dpC[i][j] = add(dpC[i-1][j], dpC[i][j-C[i-1]]);
      }
      for(;j>=0;j--){
        dpC[i][j] = dpC[i-1][j];
      }
    }
  }


  PL[0] = 0;
  for(i=1;i<=P[S-1]*300;i++){
    int x = dpP[S][i];
    PL[i] = x%B + x/B * I;
  }
  QL[0] = 0;
  for(i=1;i<=C[S-1]*300;i++){
    int x = dpC[S][i];
    QL[i] = x%B + x/B * I;
  }

  init(K, w);

  fft_real(PL, PH, K-1);
  fft_real(QL, QH, K-1);

  for(i=0;i<M/2;i++){
    cmplx tmp = PL[i] * QL[i] + PH[i] * QH[i] * I;
    QH[i] = PL[i] * QH[i] + PH[i] * QL[i];
    QL[i] = tmp;
  }

  ifft(QL, K-1);
  ifft(QH, K-1);

  int z = PL[0] = QL[0] = 1;
  for(i=1;i<M/2;i++){
    int x = add(add(pos((long long int)round(creal(QL[i]))%mod), mul(B, pos((long long int)round(creal(QH[i]))%mod))), mul(B*B%mod, pos((long long int)round(cimag(QL[i]))%mod)));
    x = sub(0, x);
    QL[i] = x%B + x/B * I;
    z = add(z, x);
    PL[i] = z%B + z/B * I;   
  }
  fft_real(PL, PH, K);
  fft_real(QL, QH, K);

  n += M/2-1;

  for(; ; n >>= 1){
    for(i=0;i<M;i++){
      cmplx tmp = PL[i] * QL[i^1] + PH[i] * QH[i^1] * (B*B%mod - mod);
      PH[i] = PL[i] * QH[i^1] + PH[i] * QL[i^1];
      PL[i] = tmp;
    }
    if(n&1){
      for(i=0;i<M/2;i++){
        PL[i] = ((PL[2*i] - PL[2*i+1]) + (PH[2*i] - PH[2*i+1]) * I) / 2 * conj(w[i]);
      }
    }
    else {
      for(i=0;i<M/2;i++){
        PL[i] = ((PL[2*i] + PL[2*i+1]) + (PH[2*i] + PH[2*i+1]) * I) / 2;
      }
    }
    if(n==1){
      for(i=1;i<M/2;i++){
        PL[0] += PL[i];
      }
      PL[0] /= M/2;
      break;
    }
    ifft(PL, K-1);
    for(i=0;i<M/2;i++){
      int x = add(pos((long long int)round(creal(PL[i]))%mod), mul(B, pos((long long int)round(cimag(PL[i]))%mod)));
      PL[i] = x%B + x/B * I;
      PL[M/2+i] = 0;
    }
    fft_real(PL, PH, K);

    for(i=0;i<M/2;i++){
      QL[i] = QL[2*i] * QL[2*i+1] + QH[2*i] * QH[2*i+1] * (B*B%mod - mod) + (QL[2*i] * QH[2*i+1] + QH[2*i] * QL[2*i+1]) * I;
    }
    ifft(QL, K-1);
    for(i=0;i<M/2;i++){
      int x = add(pos((long long int)round(creal(QL[i]))%mod), mul(B, pos((long long int)round(cimag(QL[i]))%mod)));
      QL[i] = x%B + x/B * I;
      QL[M/2+i] = 0;
    }

    fft_real(QL, QH, K);
  }

  printf("%d\n", add(pos((long long int)round(creal(PL[0]))%mod), mul(B, pos((long long int)round(cimag(PL[0]))%mod))));
  return 0;
}
0