从任何基数的比率扩展中获取特定数字（x/y 的第 n 个数字）

Posted 2023-02-22

【中文标题】从任何基数的比率扩展中获取特定数字（x/y 的第 n 个数字）

【英文标题】：Getting a specific digit from a ratio expansion in any base (nth digit of x/y)

【发布时间】：2009-04-30 00:43:53

【问题描述】：

有没有一种算法可以不从头开始计算重复小数比的位数？

我正在寻找一种不使用任意大小的整数的解决方案，因为这应该适用于十进制扩展可能任意长的情况。

例如，33/59 扩展为具有 58 位数字的重复小数。如果我想验证这一点，我该如何计算从第 58 位开始的数字？

已编辑 - 比率为 2124679 / 2147483647，如何获得第 2147484600 到第 2147484700 位中的百位数字。

【问题讨论】：

“我正在寻找一种不使用任意大小整数的解决方案”——我的猜测是，如果不使用任意精度的整数运算，这是不可能的

更正：您需要整数运算来处理除法和余数，最多为分母的 10 倍。

【参考方案1】：

好的，第三次尝试很有趣:)

我不敢相信我忘记了模幂运算。

所以从我的第二个答案中窃取/总结，x/y 的第 n 个数字是 (10^n-1x mod y)/y = floor(10 * (10 ^n-1x mod y) / y) mod 10.

一直占用的部分是 10^n-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() 的结构和模幂运算是相同的；两者都是基于Russian peasant multiplication。） 结果：

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

【参考方案2】：

编辑：（我将把帖子留给后代。但请不要再投票了：它可能在理论上有用，但实际上并不实用。我已经发布了另一个从实际角度来看更有用的答案，不需要任何因式分解，也不需要使用 bignums。）

---

我认为@Daniel Bruckner 的方法是正确的。 （需要一些额外的曲折）

也许有更简单的方法，但以下方法总是有效的：

除了 33/59 之外，让我们使用示例 q = x/y = 33/57820 和 44/65，原因可能很快就会清楚。

第 1 步：分解分母（特别是分解出 2 和 5）

写 q = x/y = x/(2^a₂5^a₅z)。分母中的因数 2 和 5 不会导致小数重复。所以剩下的因子 z 是 10 的互质数。事实上，下一步需要对 z 进行因子分解，所以你不妨考虑整个事情。

计算 a₁₀ = max(a₂, a₅) 这是 10 的最小指数，它是因子的倍数y 中的 2 和 5。

在我们的例子中 57820 = 2 * 2 * 5 * 7 * 7 * 59，所以 a₂ = 2, a₅ = 1, a₁₀ = 2，z = 7 * 7 * 59 = 2891。

在我们的示例 33/59 中，59 是素数，不包含 2 或 5 的因数，因此 a₂ = a₅ = a_{10 = 0。}

在我们的示例中，44/65，65 = 5*13，a₂ = 0，a₅ = a₁₀ = 1 .

仅供参考，我找到了一个很好的在线保理计算器here。 （甚至对下一步很重要的totients）

第 2 步：使用 Euler's Theorem 或 Carmichael's Theorem。

我们想要的是一个数 n，使得 10ⁿ - 1 可以被 z 整除，或者换句话说，10ⁿ ≡ 1 mod z。 Euler 函数 φ(z) 和 Carmichael 函数 λ(z) 都将为您提供有效的 n 值，其中 λ(z) 为您提供较小的数字，而 φ(z) 可能更容易计算。这并不难，它只是意味着分解 z 并做一些数学运算。

φ(2891) = 7 * 6 * 58 = 2436

λ(2891) = lcm(7*6, 58) = 1218

这意味着 10²⁴³⁶ ≡ 10¹²¹⁸ ≡ 1 (mod 2891)。

对于更简单的分数 33/59，φ(59) = λ(59) = 58，所以 10⁵⁸ ≡ 1 (mod 59)。

对于 44/65 = 44/(5*13)，φ(13) = λ(13) = 12。

那又怎样？嗯，小数点的周期必须同时除以 φ(z) 和 λ(z)，因此它们有效地为您提供了小数点周期的上限。

第 3 步：更多数字运算

让我们使用 n = λ(z)。如果我们从 Q'' = 10^{(a₁₀ + n) 中减去 Q' = 10a₁₀x/y sup>x/y，我们得到：}

m = 10^a₁₀(10ⁿ - 1)x/y

这是一个整数，因为 10^a₁₀ 是 y 的 2 和 5 的因数的倍数，而 10 ⁿ-1 是 y 其余因数的倍数。

我们在这里所做的是将原始数字 q 左移 a₁₀ 位得到 Q'，然后将 q 左移 a₁₀ + n 位得到Q''，它们是重复的小数，但它们之间的差是一个我们可以计算的整数。

那么我们可以将 x/y 改写为 m / 10^a₁₀ / (10ⁿ - 1)。

考虑示例 q = 44/65 = 44/(5*13)

a₁₀ = 1，λ(13) = 12，所以 Q' = 10¹q 和 Q'' = 10^12+1问。

m = Q'' - Q' = (10¹² - 1) * 10¹ * (44/65) = 153846153846*44 = 6769230769224

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

其他分数 33/57820 和 33/59 导致更大的分数。

第 4 步：找出不重复和重复的小数部分。

请注意，对于介于 1 和 9 之间的 k，k/9 = 0.kkkkkkkkkkkkkk...

同样注意，1 到 99 之间的两位数 kl，k/99 = 0.klklklklklkl...

这概括了：对于 k 位模式 abc...ij，这个数字 abc...ij/(10^k-1) = 0.abc...ijabc...ijabc ...我...

如果你遵循这个模式，你会发现我们要做的就是把我们在上一步得到的这个（可能的）巨大整数 m，写成 m = s*(10ⁿ -1) + r，其中 1 ≤ r n-1.

这就引出了最终的答案：

s 是非重复部分 r 是重复部分（必要时在左边补零以确保它是 n 位） 带有₁₀ = 0，小数点在 非重复和重复部分；如果 a₁₀ > 0 则位于 a₁₀ 个位置在 s 和 r 之间的连接点。

对于 44/65，我们得到 6769230769224 = 6 * (10¹²-1) + 769230769230

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

您可以通过在步骤 2 中找到 n 的最小值来使数字更小，从 Carmichael 函数 λ(z) 开始并查看它的任何因子是否导致 n 值使得 10^{n ≡ 1 (mod z)。}

更新：对于好奇的人来说，Python interpeter 似乎是使用 bignums 进行计算的最简单方法。 （pow(x,y) 计算 x^y， // 和 % 分别是整数除法和余数。）这里有一个例子：

>>> 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'

【讨论】：

我仍然在努力思考这里的数学问题。似乎做大数数学比仅仅做“长除法”来得到数字然后用循环检测或重复余数找到重复更复杂。长除法给出“分母”余数之一，在某些序列中，我希望有一种聪明的方法可以捷径并直接得到所需的余数。

除了循环检测永远不会仅仅看数字是万无一失的：0.12345123451234512345真的足以说重复数字是12345吗？如果您继续前进并得到 0.123451234512345123456123451234512345123456，那么您会得到不同的答案。

但是您确实有关于长除法的观点：涉及的状态较少。 （33/57820 的例子很好地说明了这一点！）我可能会尝试将长除法与 Carmichael 函数结合起来，保证给你一个循环长度的倍数。

不过，您必须知道非重复部分的结束位置。

@JasonS 当余数（不是数字！）形成一个循环时，您已经找到了完整的循环小数部分。

【参考方案3】：

作为一种通用技术，有理分数有一个不重复的部分，后面跟着一个重复的部分，如下所示：

nnn.xxxxxxxxxrrrrrr

xxxxxxxx 是非重复部分，rrrrrr 是重复部分。

确定非重复部分的长度。 如果有问题的数字在非重复部分，则直接使用除法计算。 如果有问题的数字在重复部分中，请计算其在重复序列中的位置（您现在知道所有内容的长度），然后选出正确的数字。

以上是一个粗略的大纲，需要更精确才能在实际算法中实现，但它应该可以帮助您入门。

【讨论】：

第 1 步不涉及将数字计算为重复点的 2 倍吗？第 2 步需要大整数的除法？

步骤 1 可以在不实际计算数字本身的情况下完成。要做到这一点，找到最小的数字，当它乘以分母时，得到一个形式为 999000 的数字，即。一些 9 后跟 0。 0的个数是非重复部分的长度，9的个数是重复部分的长度。

【参考方案4】：

啊哈！ 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;

还有一个运行第n位1/700的例子：

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 # 没有溢出）：

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 或更多这样的分母来说似乎仍然很残酷中间余数。

不是以N为分母的O(N)，而是以N为周期的O(N)。

哎呀，我的意思是 N = 句点 + 非重复长度。

比素数两边，句号不是和分母一样吗？我想 2124679/214748367 不会耗尽内存，只是需要很多耐心？对于质数分母，任何分数都会产生相同的序列，只是旋转？

如果分母是素数，则句点 = 分母 - 1。但对于更严格的示例，214748367 = 3*89*191*4211 并且当我运行脚本时，我得到的句点为 1759780。

【参考方案5】：

特设我没有好主意。也许continued fractions 可以提供帮助。我会考虑一下...

更新

从Fermat's little theorem 开始，因为 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 互质，这才是最重要的。 （但你是对的，应该编辑答案以反映这一点）

哎呀，这还不够，费马小定理只适用于素数指数。

Uuuuppps！当我回答这个问题时，已经是凌晨 04:00，我已经很累了。这或许可以解释为什么我认为 39 是素数……；）我将在周末解决这个问题。

a=1 b=6 p=1，但 k=(10^p-1)/b=9/6=1.5 不是整数。您的解决方案是否仅适用于质数？