結果

問題 No.377 背景パターン
コンテスト
ユーザー mudbdb
提出日時 2016-06-08 04:19:37
言語 C90
(gcc 12.3.0)
結果
CE  
(最新)
AC  
(最初)
実行時間 -
コード長 3,860 bytes
コンパイル時間 480 ms
コンパイル使用メモリ 27,520 KB
最終ジャッジ日時 2025-02-22 00:13:58
合計ジャッジ時間 1,091 ms
ジャッジサーバーID
（参考情報） 		judge5 / judge1
このコードへのチャレンジ
（要ログイン）
コンパイルエラー時のメッセージ・ソースコードは、提出者また管理者しか表示できないようにしております。（リジャッジ後のコンパイルエラーは公開されます）
ただし、clay言語の場合は開発者のデバッグのため、公開されます。

コンパイルメッセージ
main.c:25:1: error: C++ style comments are not allowed in ISO C90
   25 | // K^e (mod M)の計算
      | ^
main.c:25:1: note: (this will be reported only once per input file)

ソースコード

diff # 

#include <stdio.h>
#include <stdlib.h>
typedef long long int ll;
int H,W,K;
ll gcd(ll a, ll b)
{
  int c;
  if (a < b) {
    c = a;
    a = b;
    b = c;
  }
  while (b != 0) {
    c = a%b;
    a = b;
    b = c;
  }
  return a;
}
ll lcm(ll a, ll b)
{
  return a*b/gcd(a,b);
}
int M = 1000000007;
// K^e (mod M)の計算
ll pow_mod_K(ll K, ll e)
{
  if (K == 1) {
    return 1;
  }

  // K^(2^0), K^(2^1), ... , K^(2^n)
  static ll *KE = 0;
  int sz = 0;
  int i;
  if (KE == 0) {
    i = H;
    while (i != 0) { sz++; i>>=1; }
    sz++;
    i = W;
    while (i != 0) { sz++; i>>=1; }
    sz++;

    KE = (ll *)malloc(sizeof(ll)*sz);
    KE[0] = K%M;
    for (i=1; i<sz; i++) KE[i] = (KE[i-1]*KE[i-1])%M;
  }

  ll K2 = 1;
  i = 0;
  while (e != 0) {
    if (e&1) {
      K2 *= KE[i];
      K2 %= M;
    }
    i++;
    e>>=1;
  }
  return K2;
}
// B^e (mod M)の計算
ll pow_mod(ll B, ll e)
{
  if (B == 1) {
    return 1;
  }
  ll B2 = 1;
  B %= M;
  while (e != 0) {
    if (e&1) {
      B2 *= B;
      B2 %= M;
    }
    B = B*B;
    B %= M;
    e>>=1;
  }
  return B2;
}
// (2^3)^9 < 10^9 < (2^4)^9
ll P[4*9];
int E[4*9];
int Pend = 0;
int factor(ll N)
{
  ll p;
  int e;
  p=2;
  if (N%p == 0) {
    P[Pend] = p;
    e=1;
    for (N/=p; N%p==0; N/=p) e++;
    E[Pend] = e;
    Pend++;
  }
  for (p=3; p*p<=N; p+=2) {
    if (N%p == 0) {
      P[Pend] = p;
      e=1;
      for (N/=p; N%p==0; N/=p) e++;
      E[Pend] = e;
      Pend++;
    }
  }
  if (N != 1) {
    P[Pend] = N;
    E[Pend] = 1;
    Pend++;
  }
  return 0;
}
// e >= 0
ll pwr(ll d, int e)
{
  ll dd = 1;
  if (e == 0) {
  } else {
    int i;
    for (i=1; i<=e; i++) dd*=d;
  }
  return dd;
}
int cmpr(const void *a, const void *b)
{
  if (*(ll *)a < *(ll *)b) {
    return -1;
  } else if (*(ll *)a > *(ll *)b) {
    return 1;
  } else {
    return 0;
  }
}
// Nの全約数生成
void initD(int N, ll **aD, int *aDsz)
{
  Pend = 0;
  factor(N);

  int i;
  ll *D;
  int Dsz = 1;
  for (i=0; i<Pend; i++) Dsz *= (E[i]+1);
  D = (ll *)malloc(sizeof(ll)*Dsz);
  for (i=0; i<Dsz; i++) D[i] = 1;

  int j,k,l;
  int rpt = 1;
  for (i=0; i<Pend; i++) {
    for (j=0, k=0; j<Dsz; ) {
      for (l=0; l<rpt; l++) {
        D[j] *= pwr(P[i],k);
        j++;
      }
      k++;
      k %= (E[i]+1);
    }
    rpt *= (E[i]+1);
  }
  qsort(D, Dsz, sizeof(ll), cmpr);
  (*aD) = D;
  (*aDsz) = Dsz;
}
// 約数xの位置
int searchD(ll *D, int Dsz, ll x)
{
  int b,m,e;
  b = 0;
  e = Dsz-1;
  while (1) {
    m = (b + e)/2;
    if (D[m] == x) {
      break;
    } else if (D[m] < x) {
      b = m+1;
    } else { // x < D[m]
      e = m-1;
    }
  }
  return m;
}
ll T_mod(ll x, int pos, ll *D, int Dsz, ll *DP)
{
  ll T;
  if (pos == -1) {
    pos = searchD(D, Dsz, x);
  }
  if (DP[pos] == 0) {
    T = D[Dsz-1]/x;
    int i;
    for (i=0; i<Dsz; i++) {
      if ((x<D[i]) && (D[i]%x == 0)) {
        T -= T_mod(D[i], i, D, Dsz, DP);
        T %= M;
      }
    }
    T += (T<0?M:0);
    DP[pos] = T;
  } else {
    T = DP[pos];
  }
  return T;
}
int main()
{
  scanf("%d %d %d", &H, &W, &K);

  ll *H_D;
  ll *W_D;
  int H_Dsz;
  int W_Dsz;
  initD(H, &H_D, &H_Dsz);
  initD(W, &W_D, &W_Dsz);

  int i;
  ll *H_DP;
  ll *W_DP;
  H_DP = (ll*)malloc(sizeof(ll)*H_Dsz);
  W_DP = (ll*)malloc(sizeof(ll)*W_Dsz);
  for (i=0; i<H_Dsz; i++) H_DP[i] = 0;
  for (i=0; i<W_Dsz; i++) W_DP[i] = 0;

  int j,k;
  ll e,ed,B,a,b,h,w;
  ll sum = 0;
  h = H;
  w = W;
  for (i=0; i<H_Dsz; i++) {
    a = H_D[i];
    for (j=0; j<W_Dsz; j++) {
      b = W_D[j];
      ed = lcm(H/a,W/b);
      e = w*h/ed;
      B = pow_mod_K(K, e);
      B *= T_mod(a,-1,H_D,H_Dsz,H_DP);
      B %= M;
      B *= T_mod(b,-1,W_D,W_Dsz,W_DP);
      B %= M;
      sum += B;
      sum %= M;
    }
  }
  // フェルマーの小定理より
  // (W*H)^(-1)==(W*H)^(M-2) (mod M)
  printf("%lld\n", (sum*pow_mod(w*h,M-2))%M);
  return 0;
}
0