在任何进制下获取比率展开式中的特定数字(x/y的第n位数字)

8
有没有一种算法可以在不从开头开始的情况下计算出循环小数的数字?
我正在寻找一种解决方案,它不使用任意大小的整数,因为这应该适用于小数部分可能任意长的情况。
例如,33/59扩展为一个有58个数字的循环小数。如果我想要验证,如何计算从第58位开始的数字?
编辑 - 对于比率2124679 / 2147483647,如何获得在2147484600th到2147484700th位置之间的百位数字。

我正在寻找一种不使用任意大小整数的解决方案——我的猜测是,如果不使用任意精度整数算术,这是不可能的。 - Jason S
更正:您需要能够处理分母高达10倍的数字的整数算术来处理除法和余数。 - Jason S
5个回答

8

好的,第三次尝试成功了 :)

我简直不敢相信我竟然忘记了模幂运算。

所以为了窃取/总结我的第二个答案,x/y 的第 n 位数字是 (10n-1x mod y)/y 的第一位数字 = floor(10 * (10n-1x mod y) / y) mod 10。

花费所有时间的部分是 10n-1 mod y,但我们可以使用快速(O(log n))模幂运算来完成。有了这个,就不值得尝试进行循环查找算法。

但是,您需要能够执行(a * b mod y),其中 a 和 b 可能与 y 一样大。(如果 y 需要 32 位,则需要执行 32x32 乘法,然后进行 64 位 % 32 位模数,或者您需要绕过此限制的算法。请参见我的列表,因为我在 Javascript 中遇到了这个限制。)

所以这里是一个新版本。

function abmody(a,b,y)
{
  var x = 0;
  // binary fun here
  while (a > 0)
  {
    if (a & 1)
      x = (x + b) % y;
    b = (2 * b) % y;
    a >>>= 1;
  }
  return x;
}

function digits2(x,y,n1,n2)
{
  // the nth digit of x/y = floor(10 * (10^(n-1)*x mod y) / y) mod 10.
  var m = n1-1;
  var A = 1, B = 10;
  while (m > 0)
  {
    // loop invariant: 10^(n1-1) = A*(B^m) mod y

    if (m & 1)
    {
      // A = (A * B) % y    but javascript doesn't have enough sig. digits
      A = abmody(A,B,y);
    }
    // B = (B * B) % y    but javascript doesn't have enough sig. digits
    B = abmody(B,B,y);
    m >>>= 1;
  }

  x = x %  y;
  // A = (A * x) % y;
  A = abmody(A,x,y);

  var answer = "";
  for (var i = n1; i <= n2; ++i)
  {
    var digit = Math.floor(10*A/y)%10;
    answer += digit;
    A = (A * 10) % y;
  }
  return answer;
}

请注意,abmody() 和模幂运算的结构相同;两者都基于俄罗斯农民乘法

js>digits2(2124679,214748367,214748300,214748400)
20513882650385881630475914166090026658968726872786883636698387559799232373208220950057329190307649696
js>digits2(122222,990000,100,110)
65656565656
js>digits2(1,7,1,7)
1428571
js>digits2(1,7,601,607)
1428571
js>digits2(2124679,2147483647,2147484600,2147484700)
04837181235122113132440537741612893408915444001981729642479554583541841517920532039329657349423345806

4

编辑:(我将此帖子留在这里,以备后人参考。但请不要再点赞它:尽管理论上有用,但实际上并不实用。我已经发布了另一个答案,从实用的角度来看更加有用,不需要任何因数分解,并且不需要使用大整数。)


@Daniel Bruckner 的方法是正确的,我认为。(需要一些额外的扭曲)

也许有一种更简单的方法,但以下方法始终有效:

步骤1:因数分解分母(特别是因式分解2和5)

写成 q = x/y = x / (2a25a5z)。分母中的2和5不会导致重复的小数。因此,剩余的因子 z 与10互质。事实上,下一步需要对 z 进行因数分解,因此最好将整个分母进行因数分解。

计算 a10 = max(a2, a5),即因数2和5在 y 中的最小指数,其是10的倍数。

在我们的例子中,57820 = 2 * 2 * 5 * 7 * 7 * 59,因此 a2 = 2,a5 = 1,a10 = 2,z = 7 * 7 * 59 = 2891。

在我们的例子中,33/59,59是质数,不包含2或5的因子,因此 a2 = a5 = a10 = 0。

在我们的例子中,44/65,65 = 5*13,a2 = 0,a5 = a10 = 1。

仅供参考,我在这里找到了一个好的在线因数分解计算器(链接)。(甚至可以进行欧拉定理的计算,这对下一步很重要)

步骤2:使用欧拉定理Carmichael 定理

我们需要的是一个数字n,使得10的n次方减1能够被z整除,或者换句话说,10的n次方模z等于1。欧拉函数φ(z)和Carmichael函数λ(z)都会给出有效的n值,其中λ(z)会给出较小的数值,而φ(z)可能会更容易计算。这不太难,只需要分解z并进行一些计算即可。
φ(2891)=7*6*58=2436
λ(2891)=lcm(7*6,58)=1218
这意味着10的2436次方模2891,10的1218次方模2891都等于1。
对于简单的分数33/59,φ(59)=λ(59)=58,因此10的58次方模59等于1。
对于44/65=44/(5*13),φ(13)=λ(13)=12。
那又怎样呢?嗯,重复十进制的周期必须同时被φ(z)和λ(z)整除,因此它们实际上给出了重复十进制周期的上限。
步骤3:更多的数字计算
让我们使用n=λ(z)。如果我们从Q''=10的(a10+n)次方x/y中减去Q'=10的a10次方x/y,我们可以得到:
m=10的a10次方*(10的n次方-1)x/y
这是一个整数,因为10的a10次方是y的2和5的因子的倍数,而10的n次方-1是y的其余因子的倍数。
我们在这里所做的是将原始数字q向左移动a10个位置,以得到Q',并将q向左移动a10+n个位置,以得到Q'',这些都是重复的十进制数,但它们之间的差值是我们可以计算的整数。
然后,我们可以将x/y重写为m/10的a10次方/(10的n次方-1)。
考虑一下q=44/65=44/(5*13)的例子:
a10=1,λ(13)=12,所以Q'=10的1次方*q,Q''=10的(12+1)次方*q。

m = Q'' - Q' = (1012 - 1) * 101 * (44/65) = 153846153846*44 = 6769230769224

所以q = 6769230769224 / 10 / (1012 - 1)。

其他的分数33/57820和33/59会得到更大的分数。

步骤4:找出非循环和循环小数部分。

注意,对于1到9之间的k,k/9=0.kkkkkkkkkkkkk...

同样注意到2位数字kl在1到99之间,k/99=0.klklklklklkl...

这个是可以推广的:对于k位模式abc...ij,这个数字abc...ij/(10k-1)=0.abc...ijabc...ijabc...ij...

如果你按照这个模式进行计算,你会发现我们需要做的是将上一步得到的(潜在的)巨大整数m写成m=s*(10n-1)+r,其中1≤r<10n-1。

这就导致了最终答案:

  • s是非重复部分
  • r是重复部分(如果必要,在左侧填充零以确保它是n位数)
  • 对于a10=0,小数点位于非重复和重复部分之间;如果a10>0,则位于s和r交界处左侧a10个位置。

对于44/65,我们得到6769230769224=6*(1012-1)+769230769230

s=6,r=769230769230,44/65=0.6769230769230,这里的下划线表示重复部分。

你可以通过在步骤2中找到最小的n值来使数字更小,从Carmichael函数λ(z)开始,并查看它的因子是否导致n的值满足10n≡1(mod z)。

更新: 对于好奇的人,Python解释器似乎是用大数计算最简单的方法。 (pow(x,y)计算xy,//和%分别是整数除法和余数。)以下是一个示例:

>>> N = pow(10,12)-1
>>> m = N*pow(10,1)*44//65
>>> m
6769230769224
>>> r=m%N
>>> r
769230769230
>>> s=m//N
>>> s
6
>>> 44/65
0.67692307692307696

>>> N = pow(10,58)-1
>>> m=N*33//59
>>> m
5593220338983050847457627118644067796610169491525423728813
>>> r=m%N
>>> r
5593220338983050847457627118644067796610169491525423728813
>>> s=m//N
>>> s
0
>>> 33/59
0.55932203389830504

>>> N = pow(10,1218)-1
>>> m = N*pow(10,2)*33//57820
>>> m
57073676928398478035281909373919059149083362158422691110342442061570390868211691
45624351435489450017295053614666205465236942234520927014873746108612936700103770
32168799723279142165340712556208924247665167762020062262193012798339674852992044
27533725354548599100657212037357315807679003804911795226565202352127291594603943
27222414389484607402282947077135939121411276374956762365963334486336907644413697
68246281563472846765824974057419578000691802144586648218609477689380837080594949
84434451746800415081286751988931165686613628502248356969906606710480802490487720
51193358699411968177101349014181943964026288481494292632307160152196471809062608
09408509166378415773088896575579384296091317883085437564856451054998270494638533
37945347630577654790729851262538913870632998962296783120027672085783465928744379
10757523348322379799377378069872016603251470079557246627464545140089934278796264
26841923209961950882047734347976478727084053960567277758561051539259771705292286
40608785887236250432376340366655136630923555863023175371843652715323417502594258
04219993081978554133517813905223106191629194050501556554825319958491871324801106
88343133863714977516430300933932895191975095122794880664130058803182289865098581
80560359737115185
>>> r=m%N
>>> r
57073676928398478035281909373919059149083362158422691110342442061570390868211691
45624351435489450017295053614666205465236942234520927014873746108612936700103770
32168799723279142165340712556208924247665167762020062262193012798339674852992044
27533725354548599100657212037357315807679003804911795226565202352127291594603943
27222414389484607402282947077135939121411276374956762365963334486336907644413697
68246281563472846765824974057419578000691802144586648218609477689380837080594949
84434451746800415081286751988931165686613628502248356969906606710480802490487720
51193358699411968177101349014181943964026288481494292632307160152196471809062608
09408509166378415773088896575579384296091317883085437564856451054998270494638533
37945347630577654790729851262538913870632998962296783120027672085783465928744379
10757523348322379799377378069872016603251470079557246627464545140089934278796264
26841923209961950882047734347976478727084053960567277758561051539259771705292286
40608785887236250432376340366655136630923555863023175371843652715323417502594258
04219993081978554133517813905223106191629194050501556554825319958491871324801106
88343133863714977516430300933932895191975095122794880664130058803182289865098581
80560359737115185
>>> s=m//N
>>> s
0
>>> 33/57820
0.00057073676928398479

使用重载的Python %字符串运算符来进行零填充,可以查看完整的重复数字集:

>>> "%01218d" % r
'0570736769283984780352819093739190591490833621584226911103424420615703908682116
91456243514354894500172950536146662054652369422345209270148737461086129367001037
70321687997232791421653407125562089242476651677620200622621930127983396748529920
44275337253545485991006572120373573158076790038049117952265652023521272915946039
43272224143894846074022829470771359391214112763749567623659633344863369076444136
97682462815634728467658249740574195780006918021445866482186094776893808370805949
49844344517468004150812867519889311656866136285022483569699066067104808024904877
20511933586994119681771013490141819439640262884814942926323071601521964718090626
08094085091663784157730888965755793842960913178830854375648564510549982704946385
33379453476305776547907298512625389138706329989622967831200276720857834659287443
79107575233483223797993773780698720166032514700795572466274645451400899342787962
64268419232099619508820477343479764787270840539605672777585610515392597717052922
86406087858872362504323763403666551366309235558630231753718436527153234175025942
58042199930819785541335178139052231061916291940505015565548253199584918713248011
06883431338637149775164303009339328951919750951227948806641300588031822898650985
8180560359737115185'

我仍在试图理解这里的数学。似乎进行大数字计算比仅使用“长除法”获取数字并然后使用周期检测或重复余数来查找重复数字更为复杂。长除法方法给出了 '分母' 余数中的一个,以某个序列方式,希望有聪明的方法可以快速到达所需余数。 - caffiend
除非通过查看数字本身,否则循环检测永远不会是百分之百可靠的:0.12345123451234512345足以说明重复数字是12345吗?如果你继续计算并得到0.123451234512345123456123451234512345123456,那么你会得到另一个答案。 - Jason S
但是你对长除法确实有一点道理:它所涉及的状态较少。(33/57820的例子很好地说明了这一点!)我可能会尝试将长除法与卡迈克尔函数相结合,后者可以保证给出一个周期长度的倍数。 - Jason S
但是,您必须知道非重复部分何时结束。 - Jason S
@JasonS 当余数(不是数字!)形成一个循环时,你就找到了完整的循环小数部分。 - orlp
@orlp 噢,是的,没错。 - Jason S

2
AHA!caffiend:您对我另一个(更长的)答案的评论(具体来说是“重复余数”)引导我找到了一个非常简单的解决方案,其时间复杂度为O(n),其中n = 非重复部分和重复部分的长度之和,并且仅需要使用介于0和10*y之间的整数进行数学计算,其中y是分母。
以下是一个JavaScript函数,用于获取有理数x/y小数点右侧的第n位数字:
function digit(x,y,n) 
{ 
   if (n == 0) 
      return Math.floor(x/y)%10; 
   return digit(10*(x%y),y,n-1);
}

这是一个递归而不是迭代的算法,不能智能地检测循环(1/3的第10000位数字显然为3,但此过程会一直进行直到达到第10000次迭代),但至少可以工作直到堆栈耗尽内存。
基本上这个算法的原理有两点:
- x/y 的第n位数字是 10x/y 的 (n-1) 位数字(例如:1/7 的第6位数字是 10/7 的第5位数字,是 100/7 的第4位数字等); - x/y 的第n位数字是 (x%y)/y 的第n位数字(例如:10/7 的第5位数字也是 3/7 的第5位数字)。
我们可以将其调整为迭代例程,并与Floyd's cycle-finding algorithm(我从 Martin Gardner 的专栏中学到的“rho”方法)相结合,以获得一些可简化此方法的东西。
下面是一个使用这种方法计算解决方案的 JavaScript 函数:
function digit(x,y,n,returnstruct)
{
  function kernel(x,y) { return 10*(x%y); }

  var period = 0;
  var x1 = x;
  var x2 = x;
  var i = 0;
  while (n > 0)
  {
    n--;
    i++;
    x1 = kernel(x1,y); // iterate once
    x2 = kernel(x2,y);
    x2 = kernel(x2,y); // iterate twice  

    // have both 1x and 2x iterations reached the same state?
    if (x1 == x2)
    {
      period = i;
      n = n % period;
      i = 0; 
      // start again in case the nonrepeating part gave us a
      // multiple of the period rather than the period itself
    }
  }
  var answer=Math.floor(x1/y);
  if (returnstruct)
    return {period: period, digit: answer, 
      toString: function() 
      { 
        return 'period='+this.period+',digit='+this.digit;
      }};
  else
    return answer;
}

“计算1/700的第n位数字”的示例:”
js>1/700
0.0014285714285714286
js>n=10000000
10000000
js>rs=digit(1,700,n,true)
period=6,digit=4
js>n%6
4
js>rs=digit(1,700,4,true)
period=0,digit=4

33/59也是同样的事情:
js>33/59
0.559322033898305
js>rs=digit(33,59,3,true)
period=0,digit=9
js>rs=digit(33,59,61,true)
period=58,digit=9
js>rs=digit(33,59,61+58,true)
period=58,digit=9

"122222/990000(长的非循环部分):"
js>122222/990000
0.12345656565656565
js>digit(122222,990000,5,true)
period=0,digit=5
js>digit(122222,990000,7,true)
period=6,digit=5
js>digit(122222,990000,9,true)
period=2,digit=5
js>digit(122222,990000,9999,true)
period=2,digit=5
js>digit(122222,990000,10000,true)
period=2,digit=6

这里有另一个函数,用于查找连续的数字序列:
// find digits n1 through n2 of x/y
function digits(x,y,n1,n2,returnstruct)
{
  function kernel(x,y) { return 10*(x%y); }

  var period = 0;
  var x1 = x;
  var x2 = x;
  var i = 0;
  var answer='';
  while (n2 >= 0)
  {
    // time to print out digits?
    if (n1 <= 0) 
      answer = answer + Math.floor(x1/y);

    n1--,n2--;
    i++;
    x1 = kernel(x1,y); // iterate once
    x2 = kernel(x2,y);
    x2 = kernel(x2,y); // iterate twice  

    // have both 1x and 2x iterations reached the same state?
    if (x1 == x2)
    {
      period = i;
      if (n1 > period)
      {
        var jumpahead = n1 - (n1 % period);
        n1 -= jumpahead, n2 -= jumpahead;
      }
      i = 0; 
      // start again in case the nonrepeating part gave us a
      // multiple of the period rather than the period itself
    }    
  }
  if (returnstruct)
    return {period: period, digits: answer, 
      toString: function() 
      { 
        return 'period='+this.period+',digits='+this.digits;
      }};
  else
    return answer;
}

我已经包含了您的答案结果(假设Javascript #'s没有溢出):
js>digit(1,7,1,7,true)
period=6,digits=1428571
js>digit(1,7,601,607,true)
period=6,digits=1428571
js>1/7
0.14285714285714285
js>digit(2124679,214748367,214748300,214748400,true)
period=1759780,digits=20513882650385881630475914166090026658968726872786883636698387559799232373208220950057329190307649696
js>digit(122222,990000,100,110,true)
period=2,digits=65656565656

O(N) 对于如 214748367 或更大的分母仍然显得很高 我可以想到一些优化方法,例如使用更大的字长来减少常量因子,但是这样会更难跟踪中间的余数。 - caffiend
当比例的两侧都为质数时,周期不就是分母吗? 我认为2124679/214748367不会耗尽内存,只需要耐心等待即可。 对于质数分母,任何分数的结果序列都相同,只是旋转而已? - caffiend
如果分母是质数,周期= 分母 - 1。但对于更难的例子,214748367 = 389191*4211,当我运行我的脚本时,得到的周期为1759780。 - Jason S
是的,漏掉了一个数字 - 我想的是2^31 - 1,所以是2124679/2147483647。这个简单脚本对于1759780需要多长时间? - caffiend
但是你需要64位数学来处理2^31-1,因为你必须在执行余数之前乘以10。 - Jason S
显示剩余4条评论

2
作为一种常见的技术,有理分数由非重复部分和重复部分组成,如下所示:
nnn.xxxxxxxxrrrrrr 其中,xxxxxxxx是非重复部分,rrrrrr是重复部分。
1. 确定非重复部分的长度。 2. 如果要求的数字在非重复部分中,则直接使用除法进行计算。 3. 如果要求的数字在重复部分中,计算它在重复序列中的位置(您现在已经知道了所有内容的长度),并选择正确的数字。
以上是一个粗略的概述,在实际算法中需要更多的精确度来实现,但这应该可以让您开始工作。

第一步难道不是要计算出重复点的数字吗? 第二步需要对大分数进行除法运算? - caffiend
第一步可以在不实际计算数字本身的情况下完成。要做到这一点,找到最小的数字,使其乘以分母后得到一个形式为999000的数字,即一些由9组成的数字后跟着0。0的数量是非重复部分的长度,9的数量是重复部分的长度。 - Greg Hewgill

0

我暂时没有好主意。也许连分数可以帮助。我会再考虑一下...

更新

根据费马小定理,因为39是质数,所以以下内容成立。(=表示同余)

10^39 = 10 (39)

因为10与39互质。

10^(39 - 1) = 1 (39)

10^38 - 1 = 0 (39)

[to be continued tomorow]

我太累了,没能意识到39不是质数...^^ 我会在接下来的几天里更新答案并呈现整个想法。感谢指出39不是质数。

a/b的简短答案,其中a < b,假设周期长度为p...

  • 计算k = (10^p - 1) / b并验证它是否为整数,否则a/b没有一个周期为p
  • 计算c = k * a
  • c转换为其十进制表示,并在左侧填充零以达到总长度p
  • 小数点后第i位是填充的十进制表示的第(i mod p)位(i = 0是小数点后的第一位-我们是开发人员)

例子

a = 3
b = 7
p = 6

k = (10^6 - 1) / 7
  = 142,857

c = 142,857 * 3
  = 428,571

不需要填充,我们得出结论。

3     ______
- = 0.428571
7

但它与10互质,这才是最重要的。 (但你说得对,答案应该被编辑以反映这一点) - Jason S
哎呀,不够了,费马小定理只适用于质数指数。 - Jason S
哎呀!我回答这个问题的时候已经是凌晨4点了,而且我非常疲惫。这可能解释了为什么我认为39是质数... ;) 我打算在周末修复它。 - Daniel Brückner
a=1 b=6 p=1,但k=(10^p-1)/b=9/6=1.5不是整数。你的解决方案只适用于质数分母吗? - Dan

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接