Hatena::Grouptopcoder

naoya_t@topcoder RSSフィード

2008-12-25

SRM400 Div1 Easy: StrongPrimePower

| 18:04 | SRM400 Div1 Easy: StrongPrimePower - naoya_t@topcoder を含むブックマーク はてなブックマーク - SRM400 Div1 Easy: StrongPrimePower - naoya_t@topcoder SRM400 Div1 Easy: StrongPrimePower - naoya_t@topcoder のブックマークコメント

ある数が素数の累乗(2乗以上)かを調べる。pow()とか出てこなくて12分ぐらいかかった。(けれど一発で通った)

素数判定にはMiller-Rabinを使っている。


template <typename T> T expt(T x, int n) // x^n
{
  T y=1; for(int i=0;i<n;i++) y*=x; return y;
}

////
long long random(long long n)
{
  return (long long)rand() % n;
}
long long check_nontrivial_sqrt_1(long long m, long long a)
{
  long long r = (a * a) % m;
  return (1LL < a && a < m-1 && r == 1)? 0LL : r;
}
long long expmod(long long base, long long exp, long long m)
{
  if (exp == 0)
    return 1LL;
  else if (!(exp & 1))
    return check_nontrivial_sqrt_1(m,expmod(base,exp/2,m));
  else
    return (base*expmod(base,exp-1,m)) % m;
}
bool miller_rabin_test(long long n)
{
  return expmod((1LL+random(n-1)),n-1,n) == 1LL;
}
bool fast_prime(long long n, int times)
{
  if (n == 1) return false;
  if (n == 2) return true;
  else if (!(n & 1)) return false;
  
  for (int i=times; i>0; i--)
        if (!miller_rabin_test(n)) return false;
  return true;
}
////

class StrongPrimePower {
public:
  vector<int> baseAndExponent(string n) {
    long long n_ = 0; // 2-10^18=1000^6<2^60
    for(int i=0;i<n.length();i++) {n_*=10; n_+=n[i]-'0';}

    for (int q=59;q>=2;q--) {
      double q_ = 1.0 / q;
      double d_ = pow(n_,q_); long long p = (long long)(d_ + 0.0001);
      if (p == 0) continue;
      if (!fast_prime(p,3)) continue;
      if (expt(p,q)==n_) {
        vector<int> ans(2);
        ans[0]=p; ans[1]=q; return ans;
      }
    }
    vector<int> ans; return ans;
  }
};