I implemented an algorithm which, for a given non-negative rational number r and a positive interger b, computes all of the digits of the expansion of r in base b. The algorithm itself actually outputs a function f(i: int) satisfying the equation n = ... f(-2)*b**-2 f(-1)*b**-1 f(0)*b**0 f(1)*b**1 f(2)*b**2 ..., and it computes the digits of the whole and fractional parts of r separately through two other auxiliary functions.

Below is my code in Python (3.10.2). It looks weird for Python code because I'm actually implementing the algorithm in MIT/GNU Scheme (15.3) and "sketching" it on Python. I'm showing the Python implementation instead of the Scheme one mostly because I believe it's easier to formulate this question (and actually have it answered) in Python.

# `base`: convert a non-negative fraction to any positive integer base

from fractions import Fraction

# a natural number is a integer greater than or equal to 0
def is_natural(n):
    return isinstance(n, int) and n >= 0

# a fractional is a positive rational number strictly between 0 and 1
def is_fractional(x):
    return isinstance(x, Fraction) and 0 < x < 1

# given a natural number `n` and a positive integer `b` as, compute the base
# `b` expansion of `n` as a unary procedure `f` over the integers such that
#   n = ...   f(-2)b^-2   f(-1)b^-1   f(0)b^0   f(1)b^1   f(2)b^2   ...
#
# since `n` is a  natural number, then f(i) = 0 for any integer `i` such that
# either i < 0 or i > log_b n. thus, there are only finitely many non-trivial
# values of `f` to compute
def natural_to_proc(n, b):
    if not is_natural(n):
        raise SyntaxError("%s is not a natural number." % n)
    if not (isinstance(b, int) and b > 0):
        raise SyntaxError("%s is not a positive integer." % b)

    def proc(i: int):
        if i < 0:
            return 0

        j = -1
        q = n
        while j != i:
            j  = 1
            r = q % b
            q = q // b

        return r

    return proc

# given a fractional number `x` and a positive integer `b` as, compute the base
# `b` expansion of `n` as a unary procedure `f` over the integers such that
#   n = ...   f(-2)b^-2   f(-1)b^-1   f(0)b^0   f(1)b^1   f(2)b^2   ...
#
# since `x` is a fractional number, then f(i) = 0 for any integer i >= 1 and,
# for some integer j < 0, f(i) is periodic for all integers i < j. thus, there
# are only finitely many non-trivial values of `f` to compute
def fractional_to_proc(x, b):
    if not is_fractional(x):
        raise SyntaxError("%s is not a fractional number." % x)
    if not (isinstance(b, int) and b > 0):
        raise SyntaxError("%s is not a positive integer." % b)

    n = x.numerator
    d = x.denominator

    def proc(i: int):
        if i > 0:
            return 0
        
        j = 1
        r = n
        while j != i:
            j -= 1
            q = r // d
            r = (r % d) * b

        return q

    return proc

# given a non-negative rational number `r` and a positive integer `b` as,
# compute the base `b` expansion of `n` as a unary procedure `f` over the
# integers such that
#   n = ...   f(-2)b^-2   f(-1)b^-1   f(0)b^0   f(1)b^1   f(2)b^2   ...
#
# since `r` is a rational number, then there are a natural number `n` and a
# fractional number `x` such that r = n   x. thus, the expansion of `r` is
# simply the concatenation of the expansions of `n` and `x`
def rational_to_proc(r, b):
    if not (isinstance(r, Fraction) and r >= 0):
        raise SyntaxError("%s is not a non-negative fraction." % r)
    if not (isinstance(b, int) and b > 0):
        raise SyntaxError("%s is not a positive integer." % b)

    n = r.numerator // r.denominator
    x = r - n

    def proc(i: int):
        if i < 0:
            return fractional_to_proc(x, b)(i)
        return natural_to_proc(n, b)(i)

    return proc

Now, I want to implement an auxiliary algorithm which, given such a function f(i), computes which digits are periodically repeating and discriminate them from the remaining digits. More especifically, I want it to output a list l such that l[0] is a (possibly empty) list of the digits of the whole part of r, l[1:-1] is a list of the non-periodically repeating digits of the fractional part of r if there are any (otherwise l[1:-1] is empty), and l[-1] is a list of the periodically repeating digits of the fractional part of r (if there are any, otherwise l[-1] = []), where l[0] is ordered from the least to most significant digit and the other two lists are ordered from the most to the least significant digit. For instance, the output for 2376661/33000 in base 10 should be [[2, 7], [0, 2, 0], [0, 3]] since its decimal expansion is 72.020[03].

My question is: how would I go about implementing this second algorithm? More especifically, how would I discriminate the periodically repeating digits from the others, provided I know ahead of time that the input is a non-negative rational number?

CodePudding user response：

Steps:

1. Split into whole and fractional part
2. Extract the whole digits by repeatedly dividing by b until there's nothing left.
3. Extract the fractional digits by repeatedly multiplying with b until we reach a fraction we've seen before.
4. Split the fractional digits at the place where we first saw the final fractional value.

def analyze(r, b):
    whole, fractional = divmod(r, 1)

    whole_digits = []
    while whole:
        whole, digit = divmod(whole, b)
        whole_digits.append(digit)

    fractional_digits = []
    seen = {}
    while fractional not in seen:
        seen[fractional] = len(fractional_digits)
        digit, fractional = divmod(fractional * b, 1)
        fractional_digits.append(digit)

    s = seen[fractional]
    return [whole_digits,
            fractional_digits[:s],
            fractional_digits[s:] if fractional else []]

Demo usage:

from fractions import Fraction

r = Fraction(2376661, 33000)
b = 10

print(float(r))
print(analyze(r, b))

Output (Try it online!):

72.0200303030303
[[2, 7], [0, 2, 0], [0, 3]]

Output for base 100 (I edited leading zeros in for clarity):

72.0200303030303
[[72], [02, 00], [30]]

Output for base 1000 (I edited leading zeros in for clarity):

72.0200303030303
[[072], [020], [030, 303]]