views:

660

answers:

5

Is there an algorithm that can calculate the digits of a repeating-decimal ratio without starting at the beginning?

I'm looking for a solution that doesn't use arbitrarily sized integers, since this should work for cases where the decimal expansion may be arbitrarily long.

For example, 33/59 expands to a repeating decimal with 58 digits. If I wanted to verify that, how could I calculate the digits starting at the 58th place?

Edited - with the ratio 2124679 / 2147483647, how to get the hundred digits in the 2147484600th through 2147484700th places.

A: 
Daniel Brückner
39 is not prime.
starblue
but it is coprime to 10 and that's what counts. (but you're right, the answer should be edited to reflect this)
Jason S
oops, no it's not sufficient, Fermat's Little Theorem only works on prime exponents.
Jason S
Uuuuppps! It was 04:00 AM when I answered the question and I was quite tired. This might explain why I believed that 39 is prime ... ;) I am going to fix that over the weekend.
Daniel Brückner
a=1 b=6 p=1, but k=(10^p-1)/b=9/6=1.5 not integral. will your solution only work for prime denominator?
Dan
+1  A: 

As a general technique, rational fractions have a non-repeating part followed by a repeating part, like this:

nnn.xxxxxxxxrrrrrr

xxxxxxxx is the nonrepeating part and rrrrrr is the repeating part.

  1. Determine the length of the nonrepeating part.
  2. If the digit in question is in the nonrepeating part, then calculate it directly using division.
  3. If the digit in question is in the repeating part, calculate its position within the repeating sequence (you now know the lengths of everything), and pick out the correct digit.

The above is a rough outline and would need more precision to implement in an actual algorithm, but it should get you started.

Greg Hewgill
Doesn't step 1 involve calculating the digits to 2x the repeat point?Step 2 would require division of large fractional integers?
caffiend
Step 1 can be done without actually calculating the digits themselves. To do that, find the smallest number that when multiplied by the denominator, gives a number that is of the form 999000, ie. some number of 9s followed by 0s. The number of 0s is the length of the nonrepeating part, and the number of 9s is the length of the repeating part.
Greg Hewgill
+4  A: 
Jason S
I'm still trying to wrap my head around the math here.It seems that doing the big number math is more complex than just doing the "long division" to get the digits and then finding repeats with cycle detection or duplicate remainders. The long division method gives one of 'denominator' remainders, in some sequence, I hope there's a clever way to shortcut and get straight to the desired remainder.
caffiend
Except that cycle detection is never going to be foolproof from looking at digits alone: Is 0.12345123451234512345 really enough to say that the repeat digits are 12345? If you keep going and you get 0.123451234512345123456123451234512345123456 then you get a different answer.
Jason S
But you do have a point about long division: there's less state involved. (the 33/57820 example is a good illustration of this!) I might try to combine long division with the Carmichael function which is guaranteed to give you a multiple of the cycle length.
Jason S
Still, you have to know where the nonrepeating part ends.
Jason S
A: 

AHA! caffiend: your comment to my other (longer) answer (specifically "duplicate remainders") leads me to a very simple solution that is O(n) where n = the sum of the lengths of the nonrepeating + repeating parts, and requires only integer math with numbers between 0 and 10*y where y is the denominator.

Here's a Javascript function to get the nth digit to the right of the decimal point for the rational number x/y:

function digit(x,y,n) 
{ 
   if (n == 0) 
      return Math.floor(x/y)%10; 
   return digit(10*(x%y),y,n-1);
}

It's recursive rather than iterative, and is not smart enough to detect cycles (the 10000th digit of 1/3 is obviously 3, but this keeps on going until it reaches the 10000th iteration), but it works at least until the stack runs out of memory.

Basically this works because of two facts:

  • the nth digit of x/y is the (n-1)th digit of 10x/y (example: the 6th digit of 1/7 is the 5th digit of 10/7 is the 4th digit of 100/7 etc.)
  • the nth digit of x/y is the nth digit of (x%y)/y (example: the 5th digit of 10/7 is also the 5th digit of 3/7)

We can tweak this to be an iterative routine and combine it with Floyd's cycle-finding algorithm (which I learned as the "rho" method from a Martin Gardner column) to get something that shortcuts this approach.

Here's a javascript function that computes a solution with this approach:

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;
}

And an example of running the nth digit of 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

Same thing for 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

And 122222/990000 (long nonrepeating part):

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

Here's another function that finds a stretch of digits:

// 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;
}

I've included the results for your answer (assuming that Javascript #'s didn't overflow):

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
Jason S
O(N) still seems brutal for a denominator like 214748367 or moreI can think of optimizations such as using a larger word size to reduce the constant factor, but then I think it's more difficult to track the intermediate remainders.
caffiend
It's not the O(N) where N is the denominator, but O(N) where N is the period.
Jason S
oops, I meant N = the period + the nonrepeating length.
Jason S
With both sides of the ratio prime, isn't the period the same as the denominator? I suppose that 2124679/214748367 won't run out of memory, just take lots of patience? For prime denominators, any fraction results in the same sequence, just rotated?
caffiend
if the denominator is prime, the period = denominator - 1. But your for your tougher example, 214748367 = 3*89*191*4211 and I get a period of 1759780 when I run my script.
Jason S
Yeah, omitted a digit - the thought was 2^31 - 1, so 2124679/2147483647. How long does the simple script take for 1759780?
caffiend
Several seconds, it's JSDB Javascript based on Mozilla's Spidermonkey (though JSDB is faster than Rhino). I would expect C++ to be quite a bit faster. It's O(N) and the loop operations here are pretty simple.
Jason S
But you'd need 64-bit math to handle 2^31-1, since you have to multiply by 10 before doing a remainder.
Jason S
Couldn't that multiply by 10 be done mod 2^31 - 1 ?
caffiend
True if you have a monolithic implementation of a*b mod c, or you use an algorithm which does the same w/o requiring larger intermediate products. Otherwise you need to do a*b first, then mod c.
Jason S
+3  A: 

OK, 3rd try's a charm :)

I can't believe I forgot about modular exponentiation.

So to steal/summarize from my 2nd answer, the nth digit of x/y is the 1st digit of (10n-1x mod y)/y = floor(10 * (10n-1x mod y) / y) mod 10.

The part that takes all the time is the 10n-1 mod y, but we can do that with fast (O(log n)) modular exponentiation. With this in place, it's not worth trying to do the cycle-finding algorithm.

However, you do need the ability to do (a * b mod y) where a and b are numbers that may be as large as y. (if y requires 32 bits, then you need to do 32x32 multiply and then 64-bit % 32-bit modulus, or you need an algorithm that circumvents this limitation. See my listing that follows, since I ran into this limitation with Javascript.)

So here's a new version.

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;
}

(You'll note that the structures of abmody() and the modular exponentiation are the same; both are based on Russian peasant multiplication.) And results:

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
Jason S