2008-12-25
SRM400 Div1 Easy: StrongPrimePower
ある数が素数の累乗(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;
}
};
コメント
トラックバック - https://topcoder-g-hatena-ne-jp.jag-icpc.org/n4_t/20081225