pku3358（欧拉函数的应用）二进制循环小数

http://162.105.81.212/JudgeOnline/problem?id=3358

 

题目大意：

给出十进制数p和q，将p/q表示成二进制小数，并求出二进制小数的循环节起始位置和循环节的长度。

 

输入：

p/q   （ p >= 0, q > 0 , p,q < 2e10 ）

 

输出：

Case #n: s,l ( n第几个测试用例，s为循环节起始位置，l为循环节的长度)

思路来源：http://hi.baidu.com/lonewind/blog/item/2e7475dbb794bb64d1164e4f.html

 

解题思路：

求循环节最朴素的算法是通过迭代法求出第一个重复出现的余数，如果p与q互质

令X1= p, Xk = X(k-1) * 2 mod q, 直到Xn = p, 则循环节长度为n-1，但是，这道题给的数据范围很大，如采用这种算法，会执行很多次操作，可能多达数亿次，加之求余运算相对较慢，这种算法的时间效率就不能得到保证，而此题的时间限制是１Ｓ，因而需要寻求更优的算法。

         首先要将p/q化简，通过gcd(p,q)求得二者的最大公约数，并同时除去此数，可优化下面的计算，接着先求出循环节的起始位置，由于求的是二进制小数，所以需要对分母进行分解，p/q =p/(r*2^m)=(t+k*r)/(r*2^m)=(t/r+k)*1/(2^m)。由于p与q互质可知t与r也互质，而k为整数，对求小数部分没影响，又1/(2^m)只会使小数点向左偏移，只影响循环节的起始位置，对长度没影响，由于r与t,2都素质，可知t/r的循环节起始位置为１，又小数点向左偏移m位，因此p/q的循环节起始位置为m+1。下面的问题是求出t/r的循环节长度，而这个问题可转化为求满足t*2^k = t(mod r )的最小k值，上式可进一步化简，(t*2^k) mod r= (t mod r)*(2^k mod r ) mod r=t*(2^k mod r ) mod r = t，所以问题简化为求满足2^k = 1(mod r)的最小k值。

         上述问题如果只是简单枚举，当r为素数时可能需要计算r-1次，运算相当大，于是，这里引入了欧拉定理。

         任意一个整数n都可以表示为其素因子的乘积为：n=p1^q1*p2^q2*……*pm^qm,则它的欧拉函数为：Φ(n)=p1^(q1-1)*(p1-1)*p2^(q2-1)*(p2-1)*……*pm^(qm-1)*(pm-1)，进一步简为Φ(n)=n*(1-1/p1)*(1-1/p2)*……*(1-1/pm)，如果a与n互质，则a^Φ(n) = 1(mod n)。

#include<iostream>#include<algorithm>#include<cmath>using namespace std;const int maxn = 1000000;int a[maxn];__int64 mod(__int64 a, __int64 b, __int64 n)//a^b mod n;{__int64 res = 1;while(b){if(b & 1)res = res * a % n;a = (__int64)a * a % n;b >>= 1;}return res;}__int64 gcd(__int64 a, __int64 b){if(b == 0) return a;return gcd(b, a%b);}__int64 euler(__int64 x)//单个求欧拉函数，{__int64 i, res = x;for(i=2; i<(__int64)sqrt(x*1.0) + 1; i++)if(x % i == 0){res = res / i * (i - 1);while(x % i == 0)x /= i;}if(x > 1) res = res / x * (x - 1);return res;}int main(){__int64 p, q, oulasum, ans;int t = 1, len, i, j, x;char s;while(scanf("%I64d %c %I64d",&p,&s,&q) != EOF){if(p == 0){printf("Case #%d: 0,0/n",t++); continue;}//p /= gcd(p, q);//这一行加上就wa ???q /= gcd(p, q);len = 0;while(!(q & 1))//q%2==0,求非循环长度{q >>= 1;//q = q >> 1;len++;}if(q == 1){printf("Case #%d: 1,1/n",t++); continue;}oulasum = euler(q);//求q的欧拉函数x = 0;for(i=2; i*i<=oulasum; i++)//枚举欧拉函数的因子；要是从2到oulasum枚举就会超时。{if(oulasum % i == 0){a[x++] = i;a[x++] = oulasum / i;//加快速度，优化}}a[x++] = oulasum;sort(a, a+x);for(i=0; i<x; i++)//2^k = 1(mod r)if(mod(2, a[i], q) == 1){ans = a[i];break;}/*ans = oulasum;for(i=2; i*i<=oulasum; i++){if(oulasum % i==0 && mod(2, i, q)==1)if(i < ans) ans = i;if(oulasum%i==0 && mod(2, oulasum/i, q)==1)if(oulasum/i < ans) ans = oulasum/i;}*/printf("Case #%d: %d,%d/n",t++,len+1,ans);}return 0;} 

解析：我们可以观察一下1/10这组数据，按照二进制转换法(乘二法)，我们可以得到：

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 ...这时候，发现出现了重复，那么这个重复就是我们要求的最小循环。
抽象出模型如下：对p/q

首先p'=p/gcd(p,q)

q'=q/gcd(p,q);

然后我们就是求p'*2^i == p'*2^j (mod q')   (“==”表示同余,i<j)

经过变换得到：
p'*2^i*(2^(j-i)-1) ==0 (mod q')
也就是 q' | p'*2^i*(2^(j-i)-1)
由于gcd(p',q')=1,
得到： q' | 2^i*(2^(j-i)-1)
因为2^(j-i)-1为奇数，所以q'有多少个2的幂，i就是多少，而且i就是循环开始位置的前一位。
那么令q''为q'除去2的幂之后的数
此时 q'' | 2^(j-i)-1
也就是求出x,使得 2^x ==1 (mod q'')
由于q''和2互质，所以必定有解，而且这个是属于高次同余方程，有经典的方法可以做，大家可以自己去查阅资料解决。

 

来自J_Factory的学习空间

这个题目是求两个数相除p/q，结果的小数部分用二进制表示，当q不是2的幂时，这个二进制是个无线循环的01串。
下面是个模拟小数部分按二进制表示，可以发现二进制传一定会有循环，因为p=p%q,既然是循环，又是模运行，这和p模q的阶有关，p%q的阶一定是q的欧拉函数的因子。这样转换成一个模方程:p*2^n = x(mod q),当然p和q要互素，2和q互素。在计算前把q中的2去掉，p，q同除最大公因数。然后从1开始枚举，所以的欧拉数的因子。

#include <stdio.h>#define pr printfint main(){ int i,p,q; while(scanf("%d%d",&p,&q)==2){ pr("0."); for(i=1;i<=100;i++) { p*=2; pr("%d",p/q); p=p%q; } pr("/n"); }}/*5 1920.0000011010101010101010101010101010101010101010101010101010101010101010101010101010101010101010101010*/