POJ 3358 Period of an Infinite Binary Expansion【Euler関数】

タイトルリンク:
http://poj.org/problem?id=3358
テーマ:
分数形式p/qの有理数を入力し、{x}をその有理数バイナリ形式の小数部とし、{x}にループを持たせる
性,{x}= 0.A1A2A3…Ar(Ar+1Ar+2…Ar+s)^w.ループはr+1ビットから始まり,ループ節はsである.
現在、x 1=A 1 A 2 A 3…Arを{x}のループプレフィックスと呼び、x 2=Ar+1 Ar+2…Ar+sを{x}のループ部分と呼ぶ.
ループ接頭辞の長さとループ部分の長さをできるだけ小さくします.最小の循環部分の開始位置と最小を求めます
を行ないます.
たとえば、1/10=0.000110010011(0010010011)^w、0001100110011は1/10のループプレフィックスです.
0010011は1/10の循環部分である.でも一番小さいわけではありません.1/10=0.0(0011)^wなので、0は1/10です
最小ループ接頭辞、0011は1/10の最小ループ部分である.答えは:1/10の最小サイクル部分は2位から、最小
サイクル長は4です.
考え方:
小数を2進数に変換するのは毎回2を掛けて、1より大きい整数の部分を減らして(減らないで、引き続き2を掛けます)、小数の部分だけを残します.
このデータの1/10を観察します.小数表現は煩雑すぎて、ここでは点数でバイナリの変換を実現して、法則を観察します.
1/10を2乗法で得る.
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 …
2位からの繰返しが見られたが,繰返しのサイクル節は4であり,実は最小サイクル長である.
バイナリであることから,1*2^1=1*2^5(mod 10)が分かる.
今p/qを考えます.p 1=p/gcd(p,q),q 1=q/gcd(p,q)とする.バイナリなのでp 1*2^i≡p 1*2^j
(mod q1).変換により、p 1*2^i*(2^(j−i)−1)≡0(mod q 1)が得られる.
p 1とq 1は互いに質的であるため、2^i* (2^(j-i) - 1) % q1 = 0.2^(j-i)-1は奇数、2^iは偶数、2^i
すべてはq 1に由来し,q 1にどれだけの2のべき乗があるか,iはどれだけであり,ループ開始位置の上位でもある.
q 1中の2^iを簡略化して,(2^(j-i)-1)を得る. % q1 = 0.x=j-iとすると,2^x≡1(mod q 1)を求める.
のx最小はいくらです
q 1と2は互いに(すでに2が消えている)ので,必ず解がある.オラの定理による2^φ(q1) ≡ 1(mod q1).でも
φ(q 1)必ずしも最小の値ではないが、必ずφ(q 1)の約数.列挙φ(q 1)の約数で、最小の値が見つかります.
ACコード:

#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cstring>
using namespace std;

int GCD(int a,int b)
{
    if(b == 0)
        return a;
    return GCD(b,a%b);
}

int MultiPower(int a,int b,int m)   //a^b % m
{
    int ans = 1;
    while(b > 0)
    {
        if(b&1)
            ans = (__int64)ans * a % m; //ans = ans * a % m      。
        a = (__int64)a * a % m;
        b >>= 1;
    }
    return ans;
}

int Euler(int n)        //     φ(q1)
{
    int ret = n;
    for(int i = 2; i*i <= n; ++i)
    {
        if(n%i == 0)
        {
            n /= i;
            ret = ret - ret/i;
        }
        while(n % i == 0)
            n /= i;
    }
    if(n > 1)
        ret = ret - ret/n;
    return ret;
}

int factor[1010000];

int main()
{
    int kase = 0,p,q;
    while(~scanf("%d/%d",&p,&q))
    {
        int gcd = GCD(p,q);
        p /= gcd;
        q /= gcd;
        int time = 1;
        while(q%2 == 0)
        {
            q >>= 1;
            time++;
        }
        int phi = Euler(q);
        int ans = 0;
        if(phi == 1)
            ans = 1;
        else
        {
            int id = 0;
            for(int i = 1; i*i <= phi; ++i)     //  φ(q1)
            {
                if(phi % i == 0)
                {
                    factor[id++] = i;
                    factor[id++] = phi/i;
                }
            }
            sort(factor,factor+id);
            for(int i = 0; i < id; ++i)     //      
            {
                if(MultiPower(2,factor[i],q) == 1)
                {
                    ans = factor[i];
                    break;
                }
            }
        }

        printf("Case #%d: %d,%d
",++kase,time,ans);
    }

    return 0;
}