POJ 3358 Period of an Infinite Binary Expansion ★ (数论好题：欧拉函数)

题目链接：http://poj.org/problem?id=3358 题目大意：给定一个真分数p/q，求出在此种表示下的循环起点和循环节长度：{x} = 0.a₁a₂...a_r(a_r+1a_r+2...a_r+s)^w   我们可以观察一下1/10这组数据，按照二进制转换法(乘二法)，我们可以得到： 1/10  2/10 4/10 8/10 16/10 32/10 ... 然后都分子都模10，得到： 1/10  2/10 4/10 8/10 6/10 2/10 ... 这时候，发现出现了重复，那么这个重复就是我们要求的最小循环。   规律一般化：对于给定的p/q，我们先把它化成最简真分数，即gcd(p,q)=1. 那么我们就是要找p*2^i = p*2^j (mod q)，这样就找到了循环节. 因为gcd(p,q)==1，所以可以化简模方程得：2^i*(2^(j-i)-1) = 0 (mod q) 因为2^(j-i)-1是个奇数，所以i = （q的因子中2的个数）. 得q' = q / 2^i . 那么此时就剩下2^(j-1) = 1 (mod q')，并且j-i就是循环节长度，我们记为len. 因为gcd(2, q') == 1，所以由费马小定理的欧拉推广可知，2^Φ(i) = 1 (mod q')，所以一定有解。而且由定理可知：若a,p互质，且a^x = 1 (mod p)， 那么必定有x | Φ(p). 所以最后枚举phi(i)的因子即可.  

#include 
#include 
#include 
using namespace std;
long long gcd(long long a, long long b){
    return b ? gcd(b, a%b) : a;
}
long long phi(long long n){
    long long res = n;
    for (int i = 2; i * i <= n; i ++){
        if (n % i == 0){
            res = res / i * (i - 1);
            while(n % i == 0)
                n /= i;
        }
    }
    if (n > 1){
        res = res / n * (n - 1);
    }
    return res;
}
unsigned long long quick_add_mod(unsigned long long a, unsigned long long b, unsigned long long m){
    unsigned long long res = 0, tmp = a % m;
    while(b){
        if (b & 1)
        {
            res = res + tmp;
            res = (res >= m ? res - m : res);			//用减法比用mod快
        }
        b >>= 1;
        tmp <<= 1;
        tmp = (tmp >= m ? tmp - m : tmp);
    }
    return res;
}

long long exp_mod(long long a, long long b, long long m){
    long long res = 1 % m, tmp = a % m;
    while(b){
        if (b & 1){
            res = quick_add_mod(res, tmp, m);
        }
        tmp = quick_add_mod(tmp, tmp, m);
        b >>= 1;
    }
    return res;
}
vector  factor;
void Factor(long long n){
    factor.clear();
    factor.push_back(1);
    factor.push_back(n);
    for (int i = 2; i * i <= n; i ++){
        if (n % i == 0){
            factor.push_back(i);
            factor.push_back(n / i);
        }
    }
}
int main(){
    long long p, q, caseo = 1;
    while(scanf("%I64d%*c%I64d", &p, &q) == 2){
        //化成最简分数
        long long tmp_g = gcd(p, q);
        p /= tmp_g;
        q /= tmp_g;

        long long first_bit = 1;
        long long tmp_q = q;
        while(tmp_q % 2 == 0){
            first_bit ++;
            tmp_q >>= 1;
        }
        long long length = 0;
        Factor(phi(tmp_q));
        vector  :: iterator it = factor.begin();
        sort(it, it+factor.size());
        for (size_t i = 0; i < factor.size(); i ++){
            //printf("%d\n", factor[i]);
            if (exp_mod(2, factor[i], tmp_q) == 1){
                length = factor[i];
                break;
            }
        }
        printf("Case #%I64d: %I64d,%I64d\n", caseo ++, first_bit, length);
    }
    return 0;
}