Looking for a different kind of divmod function

Question

Python's divmod function works properly, and it is almost what I want. However, its behavior with non-integer numbers needs to be slightly different for an operation that needs to be performed. When running the following code, you might see what it is that is trying to be done.

>>> function = divmod
>>> from math import pi
>>> function(pi * pi, pi) == (pi, 0)
False
>>> 

How can function be defined above such that the final expression evaluates to True, not False? If anyone can figure out how to get (pi, 0) instead of (3.0, 0.4448...), that would be the answer.

Edit 1: Now for a more complicated example, the following code should yield [3, 2, 1, 3, 2, 1].

>>> x = 1 * pi ** 5 + \
        2 * pi ** 4 + \
        3 * pi ** 3 + \
        1 * pi ** 2 + \
        2 * pi ** 1 + \
        3 * pi ** 0
>>> digits = []
>>> while x:
        x, y = function(x, pi)
        digits.append(y)


>>> digits
[0.3989191524449005, 0.2212554774328268, 2.309739581793931, 0.1504440784612413,
2.858407346410207, 1.0]
>>> 

Edit 2: The following shows code that works fine except that it has unexpected but valid output.

import math

def convert_dec_to_pi(number):
    digits = get_pi_digits(number)
    digits, remainder = correct_pi_digits(digits)
    return make_pi_string(digits, remainder)

def get_pi_digits(number):
    digits = []
    while number:
        number, digit = divmod(number, math.pi)
        digits.append(digit)
    digits.reverse()
    return digits

def correct_pi_digits(digits):
    last = len(digits) - 1
    for index, digit in enumerate(digits):
        if index < last and digit % 1 != 0:
            a, b = get_digit_options(digit, digits[index + 1])
            digits[index:index+2] = a if 0 <= a[1] < math.pi else b
    digit, remainder = divmod(digits[-1], 1)
    digits[-1] = digit
    return digits, remainder

def get_digit_options(digit, next_digit):
    a, b = math.floor(digit), math.ceil(digit)
    if a not in range(4):
        return (b, (digit - b) * math.pi + next_digit), None
    if b not in range(4):
        return (a, (digit - a) * math.pi + next_digit), None
    c, d = ((a, (digit - a) * math.pi + next_digit),
            (b, (digit - b) * math.pi + next_digit))
    return (c, d) if digit - a < 0.5 else (d, c)

def make_pi_string(digits, remainder):
    return '{} base \u03C0 + {} base 10'.format(
        ''.join(str(int(d)) for d in digits), remainder)

The following function can be used to reverse the operation and check the results.

import re

def convert_pi_to_dec(string):
    match = re.search('^(\\d+) base \u03C0 \\+ (0\\.\\d+) base 10$', string)
    if not match:
        raise ValueError()
    digits, remainder = match.groups()
    return sum(int(x) * math.pi ** y for y, x in enumerate(reversed(digits))) \
           + float(remainder)

The following code does not raise an AssertionError, so it is evident that everything is working fine.

for n in range(1, 36):
    value = convert_dec_to_pi(n)
    print(value)
    assert convert_pi_to_dec(value) == n

So then this brings me to the following example. The output can be converted back without a problem, but one would expect something slightly different.

>>> convert_dec_to_pi(math.pi * math.pi)
'30 base π + 0.44482644031997864 base 10'
>>> convert_pi_to_dec(_) == math.pi * math.pi
True
>>> 

The string should have been 100 base π + 0.0 base 10. The output is accurate but not "proper" at this point.

Edit 3: The following example may provide some extra insight into what I am after. After running a loop with different powers of π, I would expect all outputs to be 10... base π + 0.0 base 10 in their form. The results are different from this as shown below.

>>> for power in range(20):
    print(convert_dec_to_pi(math.pi ** power))


1 base π + 0.0 base 10
10 base π + 0.0 base 10
30 base π + 0.44482644031997864 base 10
231 base π + 0.8422899173517213 base 10
2312 base π + 0.6461318165449161 base 10
23122 base π + 0.029882968108176033 base 10
231220 base π + 0.0938801130760924 base 10
2312130 base π + 0.7397595138779653 base 10
23121302 base π + 0.3240230542211062 base 10
231213021 base π + 0.017948446735832846 base 10
2312130210 base π + 0.05638670840988885 base 10
23121302100 base π + 0.17714406890720072 base 10
231213021000 base π + 0.5565145054551264 base 10
2312130133130 base π + 0.6366321966964654 base 10
23121301331302 base π + 3.9032618162071486e-05 base 10
231213013313020 base π + 0.00012262302157861615 base 10
2312130133123211 base π + 0.24905356925301847 base 10
23121301331232110 base π + 0.7824248909895828 base 10
231213013312321102 base π + 0.4580601707952492 base 10
2312130133123211021 base π + 0.4390387422112354 base 10
>>> convert_pi_to_dec('2312130133123211021 base π + 0.4390387422112354 base 10')
2791563949.5978436
>>> convert_pi_to_dec('10000000000000000000 base π + 0.0 base 10')
2791563949.5978436
>>> 

Also shown is how the last two strings are equivalent, yet the output should have been in the form of the second string. I find it fascinating that the difference between 10000000000000000000 base π and 2312130133123211021 base π is 0.4390387422112354 base 10, but that difference has a large influence over the representation. The output should have been as shown below.

1 base π + 0.0 base 10
10 base π + 0.0 base 10
100 base π + 0.0 base 10
1000 base π + 0.0 base 10
10000 base π + 0.0 base 10
100000 base π + 0.0 base 10
1000000 base π + 0.0 base 10
10000000 base π + 0.0 base 10
100000000 base π + 0.0 base 10
1000000000 base π + 0.0 base 10
10000000000 base π + 0.0 base 10
100000000000 base π + 0.0 base 10
1000000000000 base π + 0.0 base 10
10000000000000 base π + 0.0 base 10
100000000000000 base π + 0.0 base 10
1000000000000000 base π + 0.0 base 10
10000000000000000 base π + 0.0 base 10
100000000000000000 base π + 0.0 base 10
1000000000000000000 base π + 0.0 base 10
10000000000000000000 base π + 0.0 base 10

Is there something that I am missing, is there a solution to this problem, or should this be considered a fool's errand?

Answer 1

You're looking for an algorithm to determine a non-integer base representation of a floating-point number.

Wikipedia describes a greedy algorithm due to Rényi and Frougny; here's an attempt at an implementation:

from math import log, floor
def expansion(x, b):
    k = int(floor(log(x) / log(b)))
    d, r = divmod(x / float(b ** k), 1)
    digits = [int(d)]
    for _ in range(k):
        d, r = divmod(b * r, 1)
        digits.append(int(d))
    def rest(b, d, r):
        while r:
            d, r = divmod(b * r, 1)
            yield int(d)
    return digits, rest(b, d, r)

This gives the lexicographically initial expansion; you can get the lexicographically terminal expansion with a little fiddling:

def expansion(x, b, greedy=True):
    if not greedy:
        m = (floor(b) / (b - 1)) - 1
    k = int(floor(log(x) / log(b)))
    d, r = divmod(x / float(b ** k), 1)
    if not greedy and r < m:
        d, r = d - 1, r + 1
    digits = [int(d)]
    for _ in range(k):
        d, r = divmod(b * r, 1)
        if not greedy and r < m:
            d, r = d - 1, r + 1
        digits.append(int(d))
    def rest(d, r):
        while r:
            d, r = divmod(b * r, 1)
            if not greedy and r < m:
                d, r = d - 1, r + 1
            yield int(d)
    return digits, rest(d, r)

Unfortunately this still won't quite work, as OP's expansion is non-greedy in the first digit but greedy in the last digit.

Answer 2

Recognize that floating point arithmetic is by definition imprecise. Operations like pi*pi are not guaranteed to equal the mathematical constant π^2 (for that matter math.pi is only as accurate as "available precision" - meaning it's not the right value either). Therefore it's not actually possible to do operations on floating point numbers that treat them like real numbers.

The general solution is to check for a distance from some epsilon value, but this has obvious limitations. You'll do better to re-examine your underlying requirements (why do you need real-number precision?) and try to solve the problem from a different direction.

For the example you describe, why do you need to actually use a value for π? Could you just leave the actual calculation of π off til the very end, and do your operations only on the coefficients?

For instance, store the list [3, 2, 1, 3, 2, 1] directly, and do your operations and transformations with the implicit contract that they are coefficients, then define something like:

toFloat(ls,mult):
  pow = 0
  ret = 0
  for coef in ls:
    ret += coef * mult**pow
    pow += 1
  return ret

as the very last step before printing. Even better, you could wrap this behavior up in a class (and I'd be willing to bet someone has done this before) and make __str__() do toFloat()'s behavior, so that displaying your object gives you the most precise value you can get.