【POJ3358】Period of an Infinite Binary Expansion-欧拉定理+数论好题

测试地址：Period of an Infinite Binary Expansion
题目大意：对于一个小于1的有理数L，将其写成二进制小数形式：0.a1a2...，小数部分无限延伸下去（如果有限就在后面填0）。若这个小数部分可以写成下列形式：a1a2...ar(ar+1ar+2...ar+s)w，其中(s)w表示字符串s重复出现若干次，则称a1a2...ar为一个长为r的前缀，ar+1ar+2...ar+s称为一个长为s的循环。输入一个有理数p/q(0≤p<q≤2∗109)，p,q为整数，请你求这个有理数小数部分的最小循环长度，并求出循环最早在哪一位开始出现（即最小前缀长度+1）。
做法：可恶啊！明明方程都推出来了，可就是做不出来，好气啊，没办法只能看了题解……
首先我们可以将p/q约分，即令p=p/gcd(p,q),q=q/gcd(p,q)，简化计算而又不影响结果。分析题目，一个有理数小数部分的二进制表示可以用二进制转换法（乘二法）得到，那么我们可以得到下面这个序列：2i×p/q(i≥0)，将每一个数的分子模q，得到：(2i×pmodq)/q。将几组数据代入观察，发现分子出现循环，而这个循环长度正好相当于二进制位中的循环长度，那么问题就转化成了求关于i,j的同余方程：2i×p≡2i+j×p(modq)的一组解，使得i,j最小。可以发现i和j就是我们说的最小前缀长度和最小循环长度。
将式子转换一下变成：(2j−1)×p×2i≡0(modq)，也就是说q|(2j−1)×p×2i，因为gcd(p,q)=1，所以q|(2j−1)×2i，又由于2j−1是奇数，所以i就是q中包含的素因子2的数量，原因后面说明。设q′=q/2i，则q′|2j−1，那么问题转化成求关于j的同余方程2j≡1(modq′)的最小正整数解。这时候就表明前面确定i的大小的方法正确的原因了，如果这里gcd(q′,2)≠1，这个同余方程就无解了，所以要去掉q的所有因子2来确保方程必然有解。由欧拉定理得2φ(q′)≡1(modq′)，可以证明最小解j一定是φ(q′)的因子。使用反证法证明，假设最小解j不是φ(q′)的因子，那么设r=φ(q)modi，所以r<j，又因为2j≡1(modq′)且2φ(q′)≡1(modq′)，得到2r≡1(modq)，这和j是最小解相矛盾，所以j必然是φ(q′)的因子。那么我们就可以用和POJ3696一样的方法来做这道题，这道题我也写过题解，题解在此。这样我们就完美解决了这一题，最后的答案就是i+1,j。
以下是本人代码：

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>
#define ll long long
using namespace std;
ll p,q,fac[30];

ll gcd(ll a,ll b)
{
  return (b==0)?a:gcd(b,a%b);
}

ll mult(ll a,ll b,ll mod)
{
  a%=mod,b%=mod;
  ll s=a,sum=0;
  while(b)
  {
    if (b&1)
    {
      sum+=s;
      if (sum>=mod) sum-=mod;
    }
    b>>=1;s<<=1;
    if (s>=mod) s-=mod;
  }
  return sum;
}

ll power(ll a,ll b,ll mod)
{
  ll s=a,sum=1;
  while(b)
  {
    if (b&1) sum=mult(sum,s,mod);
    b>>=1;s=mult(s,s,mod);
  }
  return sum;
}

ll phi(ll x)
{
  ll p=x;
  for(int i=2;i*i<=x;i++)
    if (!(x%i))
    {
      p=p/i*(i-1);
      while(!(x%i)) x/=i;
    }
  if (x>1) p=p/x*(x-1);
  return p;
}

void find_factor(ll x)
{
  fac[0]=0;
  for(int i=2;i*i<=x;i++)
    if (!(x%i))
    {
      fac[++fac[0]]=i;
      while(!(x%i)) x/=i;
    }
  if (x>1) fac[++fac[0]]=x;
}

int main()
{
  int t=0;
  while(scanf("%lld/%lld",&p,&q)!=EOF)
  {
    t++;
    ll d=gcd(p,q),i=0,j;
    p/=d,q/=d;
    while(!(q%2)) {q>>=1;i++;}
    i++;
    j=phi(q);
    find_factor(j);
    for(int i=1;i<=fac[0];i++)
    {
      while(1)
      {
        j/=fac[i];
        if (power(2,j,q)!=1)
        {
          j*=fac[i];
          break;
        }
        else if (j%fac[i]) break;
      }
    }
    printf("Case #%d: %lld,%lld
",t,i,j);
  }

  return 0;
}