Periods of Reciprocals of Integers

[Macumazahn]

I'm tutoring a friend's son in math, and recently we've been practicing the long-division algorithm.
We started looking at the reciprocals of integers, which led to a discussion of their decimal representations.
Some of course "terminate" like 1/2 or 1/5. Others repeat forever, like 1/3 or 1/6.
For simplicity we considered the "terminating" decimal representations to be repeating decimals, with '0' repeating forever.
This made me wonder how to determine the period of the decimal representation of the reciprocal of an integer.
I'm sure this has been answered long ago for the general case, but I like to see for myself.
To this end I defined the "preamble" to be that part of the decimal that precedes the repeating part, and the "period" to be the length of one cycle of the repeating part.
The division algorithm shows both that every such reciprocal 1/n must eventually repeat, and that the maximum possible period is (n-1).
More to the point, the algorithm also reveals that the repetition of digits begins when a repeated remainder is first encountered.
So... I coded up a program, first for my workhorse 32Sii, then for my 41CX.
The challenge here is that the program must keep track of all the remainders it has seen, so that it can detect the first repetition.
That means a bitmap, which means that the range is limited by available memory.
That in turn means making the program as small as possible, while using as few registers as possible.
I was able to split the program into six routines, each stored in extended memory. These six routines load one another as required, using GETP.
I was able to squeeze the routines down to the point where the largest of the six requires 64 bytes. I just can't get it to 63 (which would improve range).
The silver lining in that cloud is that it means I can use a more capable version of one of the other programs ("Q"), as long as that version is no larger than 70 bytes.
Unlike the 32Sii which supports 36 bits per register, the 41CX supports only 30.
By squeezing the programs down to 70 bytes or less, I can execute SIZE 309 which (since I use two registers for other things) gives the program a maximum range of 30*(309-2)=9210.
The program actually runs through the division twice - first, to find the matching remainder, and second, to count off the preamble and period. I know it's slow, but I sacrificed time for space in order to maximize range.

The user interface is as follows:

If flag 00 is set, the program will display the decimal digits as it counts the preamble and period.
Flag 01 is set by the program when it completes the count of the preamble and begins counting the period, and cleared when the entire calculation is complete.
If flag 02 is set, the program will pause with TIME displayed, to allow a crude measure of elapsed time. Continue the program with R/S to see the actual results.

Place "P" in Alpha
XEQ GETP

Once the initial program "P" is loaded, functionality is as follows:

XEQ 01 = reset (reloads program "P")

Place positive integer argument n (1 <= n <= 9210) in stack-X
XEQ 02 = calculate preamble and period
Returns:
period-length in stack-X
preamble-length in stack-Y
n in stack-Z (also in Last-X)
matched-remainder in stack-T

If somehow we already know the matched-remainder, we can escape the 9210 limit as follows:

Place matched-remainder in stack-Y
Place n in stack-X
XEQ 03 = calculate preamble and period
Returns: the same stack of results as XEQ 02 above.

The program may be interrupted (and continued) as desired to change flag 00 or 02 settings.
However, the display of digits (SF 00) cannot be activated for the first time once the counting of the preamble has begun.
The program may be terminated at any time, and reset by XEQ 01.

If anyone could tell me how to squeeze program "R" down to 63 bytes without using an additional register or sacrificing functionality, I could extend the range to 9240...

Code:

1    LBL "P"
2    "P"
3    FIX 0
4    FC?C 01
5    GTO 00
6    FC? 02
7    GTO 08
8    TIME
9    TONE 9
10    FIX 4
11    STOP
12    FIX 0
13    RDN
14    LBL 08
15    LASTX
16    RCL 01
17    RCL 02
18    GTO 00
19    LBL 02
20    FIX 0
21    CLRG
22    1
23    ENTER↑
24    2
25    "R"
26    GETP
27    LBL 03
28    FIX 0
29    "T"
30    GETP
31    LBL 01
32    FIX 0
33    "P"
34    GETP
35    LBL 00
36    END

Code:

1    LBL "R"
2    "U"
3    LBL 10
4    STO 00
5    RDN
6    X<>Y
7    STO Z
8    MOD
9    ENTER↑
10    STO 01
11    30
12    ST/ Y
13    RDN
14    INT
15    ST+ 00
16    R↑
17    *
18    -
19    2
20    X<>Y
21    Y↑X
22    RCL IND 00
23    X<>Y
24    ST/ Y
25    X<>Y
26    INT
27    2
28    /
29    ENTER↑
30    INT
31    -
32    X≠0?
33    GETP
34    X<>Y
35    ST+ IND 00
36    R↑
37    RCL 01
38    10
39    *
40    2
41    GTO 10
42    LBL 01
43    "P"
44    GETP
45    END

Code:

1    LBL "U"
2    "T"
3    RCL 01
4    R↑
5    GETP
6    LBL 01
7    "P"
8    GETP
9    END

Code:

1    LBL "T"
2    "S"
3    FS? 00
4    "Q"
5    CLRG
6    CF 01
7    1
8    STO 00
9    GETP
10    LBL 01
11    CF 01
12    "P"
13    GETP
14    END

Code:

1    LBL "S"
2    "P"
3    LBL 11
4    X<>Y
5    MOD
6    LASTX
7    1
8    ST+ IND 00
9    RDN
10    X<> Z
11    X=Y?
12    GTO 00
13    LBL 13
14    STO T
15    CLX
16    10
17    *
18    GTO 11
19    LBL 00
20    FS? 01
21    GETP
22    SF 01
23    R↑
24    ST+ 00
25    RDN
26    GTO 13
27    LBL 01
28    CF 01
29    GETP
30    END

Code:

1    LBL "Q"
2    "P"
3    LBL 12
4    FS? 00
5    GTO 00
6    X<>Y
7    MOD
8    GTO 09
9    LBL 00
10    STO 03
11    X<>Y
12    ST/ Y
13    X<>Y
14    INT
15    PSE
16    X<>Y
17    *
18    ST- 03
19    CLX
20    RCL 03
21    LBL 09
22    LASTX
23    1
24    ST+ IND 00
25    RDN
26    X<> Z
27    X=Y?
28    GTO 00
29    LBL 14
30    STO T
31    CLX
32    10
33    *
34    GTO 12
35    LBL 00
36    FS? 01
37    GETP
38    SF 01
39    2
40    STO 00
41    RDN
42    GTO 14
43    LBL 01
44    CF 01
45    GETP
46    END

[Thomas Klemm]

(12-30-2017 02:19 PM)Macumazahn Wrote:  This made me wonder how to determine the period of the decimal representation of the reciprocal of an integer.

In case of \(n=13\) the period \(k=6\) and therefore we can write \(\frac{1}{13}=\frac{76,923}{999,999}\). This is a consequence of the formula to calculate the geometric series.

However this means that \(13\times76,923=999,999\) and thus in general: \(n\) divides \(10^k-1\).
For this to work properly we have to make sure that \(n\) doesn't share divisors with \(10\).

Here's a program for the HP-41:

Code:

LBL "PERIOD"
STO 00 ; n
LBL 00
RCL 00 ; n
2
MOD
X≠0?
GTO 01
LASTX
ST/ 00 ; n
GTO 00
LBL 01
RCL 00 ; n
5
MOD
X≠0?
GTO 02
LASTX
ST/ 00 ; n
GTO 01
LBL 02
CLX
STO 01 ; k
RCL 00 ; n
1
X=Y?
GTO 04
ENTER  ; 1 1
LBL 03 ; 1 p
ST+ 01 ; k'=k+1
X<>Y   ; p 1
10     ; 10 p 1
*      ; 10*p 1
RCL 00 ; n 10*p 1
MOD    ; 10*p%n 1
X<>Y   ; 1 p'
X≠Y?
GTO 03 ; while
LBL 04 ; done
RCL 01 ; k
END

This is a translation of the following Python program:

Code:

def period(n):
    while n % 2 == 0:
        n /= 2
    while n % 5 == 0:
        n /= 5
    k, p = 0, 1
    if n == 1:
        return k
    while True:
        k += 1
        p = p * 10 % n
        if p == 1:
            return k

There are probably more efficient ways to calculate the smallest integer \(k\) such that \(10^k\equiv1(\mod n)\) since it divides the Carmichael Function \(\lambda(n)\). But for this we'd need to know the prime-factors of \(n\).

Kind regards
Thomas

[Albert Chan]

(08-05-2018 08:49 PM)Thomas Klemm Wrote:  There are probably more efficient ways to calculate the smallest integer \(k\) such that \(10^k\equiv1(\mod n)\) since
it divides the Carmichael Function \(\lambda(n)\). But for this the we'd need to know the prime-factors of \(n\).

Hi, Thomas,

Your code is pretty good. Factoring number is hard.
A slight optimization is to check if p==n-1 too, and return k+k

The code with factoring (factor() not shown) look like this:

Code:

def order(a, n):             # assumed gcd(a, n) = 1
    m = n                    # m = phi(n) = k * period
    for (p,e) in factor(n): m -= m // p
    for (p,e) in factor(m):
        m //= p**e           # reduce k to 1, so m = period
        t = pow(a, m, n)
        while t != 1:        # m reduced too much
            m, e = m*p, e-1  # raise m, 1 p at a time
            if e==0: break
            t = pow(t, p, n)
    return m

Code:

def period(n):
    while n%2 == 0: n /= 2
    while n%5 == 0: n /= 5
    if n > 1: return order(10, n)

[Thomas Klemm]

Meanwhile I've extracted the redundant code into a "function" reduce:

Code:

LBL "PERIOD"
STO 00      ; n
2           ; 2
XEQ 00      ; reduce(2)
5           ; 5
XEQ 00      ; reduce(5)
CLX         ; 0
STO 01      ; k=0
RCL 00      ; n
1           ; 1 n
X=Y?        ; 1=n?
GTO 02      ; done
ENTER       ; 1 p=1
LBL 01      ; 1 p
ST+ 01      ; k'=k+1
X<>Y        ; p 1
10          ; 10 p 1
*           ; 10*p 1
RCL 00      ; n 10*p 1
MOD         ; 10*p%n 1
X<>Y        ; 1 p'
X≠Y?        ; 1=p?
GTO 01      ; while
LBL 02      ; done
RCL 01      ; k
RTN         ; return k
LBL 00      ; reduce(m)
RCL 00      ; n m
X<>Y        ; m n
MOD         ; n%m
X≠0?        ; m ∤ n
RTN         ; done
LASTX       ; m
ST/ 00      ; n'=n/m
GTO 00      ; reduce
END

If you want to also display the digits you can include the following lines:

Code:

RCL X       ; 10*p 10*p 1
RCL 00      ; n 10*p 10*p 1
/           ; 10*p/n 10*p 1 1
INT         ; ⌊10*p/n⌋ 10*p 1 1
VIEW X      ; digit
RDN         ; 10*p 1 …

However the result isn't correct for values divisible by 2 or 5.

(08-06-2018 12:59 AM)Albert Chan Wrote:  A slight optimization is to check if p==n-1 too, and return k+k

Fair point. But we'd still have to run to \(2k\) if we wanted to display the digits.

(12-30-2017 02:19 PM)Macumazahn Wrote:  The challenge here is that the program must keep track of all the remainders it has seen, so that it can detect the first repetition.

You could also use Floyd's cycle detection algorithm or a variant thereof. They use only a constant number of memory cells.

Cheers
Thomas

[Albert Chan]

(08-06-2018 03:37 AM)Thomas Klemm Wrote:  

(08-06-2018 12:59 AM)Albert Chan Wrote:  A slight optimization is to check if p==n-1 too, and return k+k

Fair point. But we'd still have to run to \(2k\) if we wanted to display the digits.

Hi, Thomas:

If p==n-1 exist, you can deduce the other half digits with 9-complements

To show why it work, i use a simpler fraction, say 1/13:
10^3 (mod 13) = 12, so period = 2*3 = 6

Firsrt 3 repeating digits of 1/13 = floor(1e3/13) = 076
To get the other half repeating digits d:

76000/999999 + d/999999 = 1/13 = 77/1001
76000 + d = 77*999 = (76+1)(1000-1) = 76000 + 1000 - 76 - 1
d = 999 - 76 = 923

-> 1/13 = 0.076 923

A harder example 3/17, 10^8 (mod 17) = 16, so period = 2*8 = 16

First 8 repeating digits of 3/17 = floor(3e8/17) = 17647058

--> 3/17 = 0.17647058 82352941

[Thomas Klemm]

(08-06-2018 02:05 PM)Albert Chan Wrote:  If p==n-1 exist, you can deduce the other half digits with 9-complements

Yes, I know. But where do you want to keep the digits once the period is longer than what the limited memory allows us to store? I was afraid that we don't gain much by this optimisation. Feel free to prove me wrong. I'd be interested to see an implementation.

Kind regards
Thomas

PS: I'd rather have my current solution display the digits as well in case of \(n\) having a common divisor with 10. That would be a requirement of Macumazahn's original post.

[Albert Chan]

I wrongly assumed the half digits had written down. Sorry.

For n != 3,6,9, here is an idea:

Instead of searching period and building digits 1 at a time, why not in pairs ?
Multiply by 100 instead of 10, stop if reminder = 1, or 10 (1 digit passed period)

Here is the search for period(n = 17):
1) 100/17 = 05 + 15/17
2) 1500/17 = 88 + 4/17
3) 400/17 = 23 + 9/17
4) 900/17 = 52 + 16/17
5) 1600/17 = 94 + 2/17
6) 200/17 = 11 + 13/17
7) 1300/17 = 76 + 8/17
8) 800/17 = 47 + 1/17 <-- stop, remainder = 1, thus period = 2*8 = 16

1/17 = 0.05 88 23 52 94 11 76 47


Unlike my previous optimization, this always cut search steps in half.
Also, it showed repeating decimals pairs without storing them.