POJ 3358 Period of an Infinite Binary Expansion 解题报告（欧拉函数+因式分解）

Period of an Infinite Binary Expansion

Time Limit: 1000MS Memory Limit: 65536KTotal Submissions: 2891 Accepted: 768

Description

Let {x} = 0.a₁a₂a₃... be the binary representation of the fractional part of a rational number z. Suppose that {x} is periodic then, we can write

{x} = 0.a₁a₂...a_r(a_r+1a_r+2...a_r+s)^w

for some integers r and s with r ≥ 0 and s > 0. Also, (a_r+1a_r+2...a_r+s)^wdenotes a nonterminating and repeating binary subsequence of {x}.

The subsequence x₁ = a₁a₂ ... a_ris called the preperiod of {x} and x₂ = a_r+1a_r+2 ... a_r+s is the period of {x}.

Suppose that |x₁| and |x₂| are chosen as small as possible then x₁ is called the least preperiod and x₂ is called the least period of {x}.

For example, x = 1/10 = 0.0001100110011(00110011)^w and 0001100110011 is a preperiod and 00110011 is a period of 1/10.

However, we can write 1/10 also as 1/10 = 0.0(0011)^w and 0 is the least preperiod and 0011 is the least period of 1/10.

The least period of 1/10 starts at the 2nd bit to the right of the binary point and the the length of the least period is 4.

Write a program that finds the position of the first bit of the least period and the length of the least period where the preperiod is also the minimum of a positive rational number less than 1.

Input

Each line is test case. It represents a rational number p/q where p and q are integers, p ≥ 0 and q > 0.

Output

Each line corresponds to a single test case. It represents a pair where the first number is the position of the first bit of the least period and the the second number is the length of the least period of the rational number.

Sample Input

1/10 1/5 101/120 121/1472

Sample Output

Case #1: 2,4 Case #2: 1,4 Case #3: 4,4 Case #4: 7,11

Source

Manila 2006

   

    解题报告： 一道非常好的数论题，考察了欧拉函数，因式分解， 费马小定理等，WA很多次，来细讲一下。

    首先，看清楚题目：给定的分数是一个小于1的正有理数，可以确定0<p<q。

    但是p和q可能存在公约数，即可以约分。我们首先分别对p，q进行因式分解，相同的素因子一起去掉。

    假设约分后的数字是1/10，那么小数的第一位是看2*1/10是否大于1，如果大于则分子减去分母，该位填1，如果小于则填0，并继续下去。

    如果b仍有素因子2，那么前几位都会被b的素因子抵消掉，这几位就是循环开始前的那几位。1/10就是0.0 00110011。

    将b的2全部约分掉后，得到a/b且b为奇数。如果b为1，小数后面都是0，循环节就是1。如果不是1，那么循环节的大小与a也无关了，而只与2^e有关。

    当2^e mod b==1时，e即为循环节。数论中称为2在模b系中的阶。

    如何求阶就是最后的问题。枚举肯定是可以的，但是如果b很大，时间上会有问题。大步小步算法亦然。

    这时可以用到数论的一个结论：一个群S的子群S'的大小 |S'| | |S|，即为S大小的约数。也就是Lagrange定理。

    模b系的群大小即为欧拉函数，我们可以对欧拉函数进行因式分解，测试pow(2, epi/fac[i], b)是否等于1，使最终的阶最小。

    表述还是有说不好的地方。数论的一些知识我是在刘汝佳老师的这本书上看的：http://pan.baidu.com/s/15tYTw

    代码如下：

#include<cstdio>#include<cstring>#include<ctime>#include<iostream>#include<algorithm>#include<map>using namespace std;typedef long long LL;const int S=40;LL mul_mod(LL a, LL b, LL mod){    a%=mod;    b%=mod;    LL res=0;    while(b)    {        if(b&1)            res+=a, res%=mod;        a<<=1, a%=mod;        b>>=1;    }    return res;}LL pow_mod(LL a, LL b, LL mod){    a%=mod;    LL res=1;    while(b)    {        if(b&1)            res=mul_mod(res, a, mod);        a=mul_mod(a, a, mod);        b>>=1;    }    return res;}bool check(LL a, LL n, LL x, LL t){    LL res = pow_mod(a, x, n);    if(res==1) return false;    LL last=res;    for(int i=1;i<=t;i++)    {        res=mul_mod(res, res, n);        if(res==1 && last!=1 && last!=n-1) return true;        last=res;    }    if(res!=1) return true;    return false;}bool Miller_Rabin(LL n){    if(n<2) return false;    if(n==2) return true;    if((n&1)==0) return false;    LL x=n-1;    LL t=0;    while((x&1)==0) x>>=1,t++;    for(int i=0;i<S;i++)        if(check(rand()%(n-1)+1, n, x, t))            return false;    return true;}long long factor[100];//质因数分解结果（刚返回时是无序的）int tol;//质因数的个数。数组小标从0开始LL gcd(LL a, LL b){    if(a==0) return 1;    if(a<0) return gcd(-a, b);    while(b)    {        LL t=a%b;        a=b;        b=t;    }    return a;}LL Pollard_rho(LL x, LL c){    LL i=1, k=2;    LL x0=rand()%x;    LL y=x0;    while(1)    {        i++;        x0=(mul_mod(x0, x0, x)+c)%x;        LL d=gcd(y-x0, x);        if(d!=1 && d!=x) return d;        if(y==x0) return x;        if(i==k) y=x0, k<<=1;    }}void findfac(LL n){    if(Miller_Rabin(n))    {        factor[tol++]=n;        return;    }    LL p=n;    while(p>=n) p=Pollard_rho(p, rand()%(n-1)+1);    findfac(p);    findfac(n/p);}LL calPhi(LL a){    tol=0;    findfac(a);    sort(factor, factor+tol);    tol=unique(factor, factor+tol)-factor;    for(int i=0;i<tol;i++)        a=a/factor[i]*(factor[i]-1);    return a;}int cas=1;void work(LL a, LL b){    printf("Case #%d: ", cas++);    map<LL, int> mm;    tol=0;    if(b>1) findfac(b);    for(int i=0;i<tol;i++)        mm[factor[i]]++;    tol=0;    if(a>1) findfac(a);    for(int i=0;i<tol;i++)        mm[factor[i]]--;    int t=mm[2]>0?mm[2]:0;    t++;    mm[2]=-1;printf("%d,", t);    b=1;    for(map<LL, int>::iterator it=mm.begin();it!=mm.end();it++) if(it->second>0)        for(int i=0;i<it->second;i++) b*=it->first;if(b==1)    {        puts("1");        return;    }    LL phi = calPhi(b);    tol=0;    if(phi>1) findfac(phi);    sort(factor, factor+tol);    LL ans=phi;    for(int i=0;i<tol;i++)        if(pow_mod(2, ans/factor[i], b)==1)            ans/=factor[i];    printf("%d\n", ans);}int main(){    LL a, b;    while(~scanf("%lld/%lld", &a, &b))        work(a, b);}