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

测试地址：Period of an Infinite Binary Expansion
题目大意：对于一个小于1的有理数$L$，将其写成二进制小数形式：$0.a_1a_2...$，小数部分无限延伸下去（如果有限就在后面填0）。若这个小数部分可以写成下列形式：$a_1a_2...a_r(a_{r+1}a_{r+2}...a_{r+s})^w$，其中$(s)^w$表示字符串$s$重复出现若干次，则称$a_1a_2...a_r$为一个长为$r$的前缀，$a_{r+1}a_{r+2}...a_{r+s}$称为一个长为$s$的循环。输入一个有理数$p/q(0\le p<q\le 2*10^9)$，$p,q$为整数，请你求这个有理数小数部分的最小循环长度，并求出循环最早在哪一位开始出现（即最小前缀长度+1）。
做法：可恶啊！明明方程都推出来了，可就是做不出来，好气啊，没办法只能看了题解……
首先我们可以将$p/q$约分，即令$p=p/\gcd(p,q),q=q/\gcd(p,q)$，简化计算而又不影响结果。分析题目，一个有理数小数部分的二进制表示可以用二进制转换法（乘二法）得到，那么我们可以得到下面这个序列：$2^i\times p/q(i\ge0)$，将每一个数的分子模$q$，得到：$(2^i\times p \mod q)/q$。将几组数据代入观察，发现分子出现循环，而这个循环长度正好相当于二进制位中的循环长度，那么问题就转化成了求关于$i,j$的同余方程：$2^i\times p\equiv 2^{i+j}\times p(\mod q)$的一组解，使得$i,j$最小。可以发现$i$和$j$就是我们说的最小前缀长度和最小循环长度。
将式子转换一下变成：$(2^j-1)\times p\times 2^i\equiv 0(\mod q)$，也就是说$q|(2^j-1)\times p\times 2^i$，因为$\gcd(p,q)=1$，所以$q|(2^j-1)\times 2^i$，又由于$2^j-1$是奇数，所以$i$就是$q$中包含的素因子$2$的数量，原因后面说明。设$q'=q/2^i$，则$q'|2^j-1$，那么问题转化成求关于$j$的同余方程$2^j\equiv1(\mod q')$的最小正整数解。这时候就表明前面确定$i$的大小的方法正确的原因了，如果这里$\gcd(q',2)\ne 1$，这个同余方程就无解了，所以要去掉$q$的所有因子$2$来确保方程必然有解。由欧拉定理得$2^{\varphi(q')}\equiv 1(\mod q')$，可以证明最小解$j$一定是$\varphi(q')$的因子。使用反证法证明，假设最小解$j$不是$\varphi(q')$的因子，那么设$r=\varphi(q)\mod i$，所以$r<j$，又因为$2^j\equiv 1(\mod q')$且$2^{\varphi(q')}\equiv 1(\mod q')$，得到$2^r\equiv 1(\mod q)$，这和$j$是最小解相矛盾，所以$j$必然是$\varphi(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\n",t,i,j);
  }

  return 0;
}