poj 3358 Period of an Infinite Binary Expansion

由乘2取整得到分数的小数位，可以找到规律！！！

例如：1/10，2/10，4/10，8/10，16/10，32/10，64/10……

取整后：1/10，2/10，4/10，8/10，6/10，2/10，4/10……

这样我们就发现规律了！！！

也就是对于p/q而言，要满足2^x=2^y mod q (gcd(p,q)==1);

化简：2^x*(2^(x-y)-1) = 0 mod q;

q里面2的倍数有多少个，就是最小的循环起始位置。

继而化简：2^(x-y) = 1 mod q' (q'除以2的倍数之后的值)

也就是求2^t = 1 mod q'

由欧拉定理知道：t=phi(q');但是这求出的t不一定是最小的，所以应该枚举t的约数，继而得到答案……

链接http://poj.org/problem?id=3358

  1 #include<iostream>
  2 #include<stdio.h>
  3 #include<algorithm>
  4 #include<iomanip>
  5 #include<cmath>
  6 #include<string>
  7 using namespace std;
  8 int prime[5000],m;
  9 bool f[5001];
 10 void init()
 11 {
 12     __int64 i,j;
 13     m=0;
 14     for(i=2;i<=5000;i++)
 15     {
 16         if(f[i]==0)
 17         {
 18             prime[m++]=i;
 19             for(j=i*i;j<=5000;j+=i)
 20                 f[j]=1;
 21         }
 22     }
 23 }
 24 __int64 gcd(__int64 a,__int64 b)
 25 {
 26     __int64 t;
 27     if(a<b) swap(a,b);
 28     while(b)
 29     {
 30         t=a;
 31         a=b;
 32         b=t%b;
 33     }
 34     return a;
 35 }
 36 __int64 euler(__int64 n)
 37 {
 38     __int64 ans=1;
 39     int i;
 40     for(i=0;i<m&&prime[i]<=n;i++)
 41     {
 42         if(n%prime[i]==0)
 43         {
 44             n/=prime[i];
 45             ans*=prime[i]-1;
 46             while(n%prime[i]==0)
 47             {
 48                 ans*=prime[i];
 49                 n/=prime[i];
 50             }
 51         }
 52     }
 53     if(n!=1)
 54         ans*=n-1;
 55     return ans;
 56 }
 57 __int64 pows(__int64 a,int b,__int64 m)
 58 {
 59     __int64 ans=1;
 60     while(b)
 61     {
 62         if(b&1)
 63             ans=ans*a%m;
 64         b>>=1;
 65         a=a*a%m;
 66     }
 67     return ans;
 68 }
 69 int main()
 70 {
 71     init();
 72     __int64 a,b,q,g,mmin;
 73     int i,j,p,k=0;
 74     while(scanf("%I64d/%I64d",&a,&b)!=EOF)
 75     {
 76         if(a==0)
 77             printf("Case #%d: %d,%I64d
",++k,1,1);
 78         g=gcd(b,a);
 79         a=a/g;b=b/g;
 80         i=1;
 81         while(b%2==0)
 82         {
 83             b/=2;
 84             i++;
 85         }
 86         q=euler(b);j=i;
 87         mmin=q;
 88         for(i=1;i*i<=q;i++)
 89         {
 90             if(q%i==0)
 91             {
 92                 if(pows(2,i,b)==1)
 93                 {
 94                     mmin=i;
 95                     break;
 96                 }
 97                 p=q/i;
 98                 if(pows(2,p,b)==1&&p<mmin) mmin=p;
 99             }
100         }
101         printf("Case #%d: %d,%I64d
",++k,j,mmin);
102     }
103     return 0;
104 }

View Code