Loss of precision in calculations involving fractions and floating-point numbers in Python

Question

import sympy as sp
from fractions import Fraction

new_matrix_aux = [['X',   'B', 'X1', 'X2', 'X3', 'X4', 'X5', 'U1', 'U2'], 
                  ['X2', 10/3,    0,    1,  2/3,  1/3, -2/3, -1/3,  2/3],  
                  ['X1',  4/3,    1,    0,  2/3, -2/3,  1/3,  2/3, -1/3]]

z_function = "Z = 1 * X1 + 2 * X2 + 3 * X3 + 0 * X4 + 0 * X5 + M * U1 + M * U2"


lista_de_operaciones_zj = []

for j in range(1, len(new_matrix_aux[0])):
    z = z_function
    for i in range(1, len(new_matrix_aux)):
        z = z.replace(new_matrix_aux[i][0], str(new_matrix_aux[i][j]))
    for k in range(1, len(new_matrix_aux[0])):
        z = z.replace(new_matrix_aux[0][k], '0')
    z = z.replace(z[0], z[0] + '_' + new_matrix_aux[0][j])
    lista_de_operaciones_zj.append(z)
    print(f'"{z}"')

variables, resultados = [], []

for operacion in lista_de_operaciones_zj:
    variable = operacion.split("Z_")[1].split(" = ")[0]
    variables.append(variable)
    expression = operacion.split("= ")[1]
    resultado = sp.sympify(expression)
    resultados.append(resultado)

output = [variables, resultados]
print(output)

This is the incorrect output that I am getting, where for some reason throughout the mathematical processes of the code the precision of the calculated values is lost (especially in those operations that involved fractions as can be seen below):

[['B', 'X1', 'X2', 'X3', 'X4', 'X5', 'U1', 'U2'], [8.0000000000000003, 1, 2, 2.000000000000000, 0, -0.9999999999999999, 0, 0.9999999999999999]]

This should be the correct output matrix where there was no loss of precision during the calculations.

[['B', 'X1', 'X2', 'X3', 'X4', 'X5', 'U1', 'U2'], [8, 1, 2, 2, 0, -1, 0, 1]]

I assume this problem is due to the nature of floating-point calculations in Python. Floating point numbers are internally represented in binary and can introduce small rounding errors.

Answer 1

Your code begins with from fractions import Fraction, but then you never use Fraction at all.

To declare a Fraction in python, write the word explicitly:

from fractions import Fraction

x = Fraction(1, 3)  # fraction
y = 1/3             # not a fraction

print('x = {}  {}'.format(x, type(x)))
print('y = {}  {}'.format(y, type(y)))
# x = 1/3  <class 'fractions.Fraction'>
# y = 0.3333333333333333  <class 'float'>

Alternatively, since you're using sympy, you can use sympy.Rational:

import sympy

x = sympy.Rational(1,3)

If you only need fractions, then use fractions.Fraction. But if you need some more stuff from sympy, then I recommend sticking to sympy entirely and using sympy.Rational.

Answer 2

So as others have said, you're not using Fraction at all – also, if I understood correctly what your new_matrix_aux stuff should do, here's a version that

1. uses proper Fractions
2. substitutes the desired variable values into a series of expressions and sympy-evaluates them:

from fractions import Fraction

import sympy as sp


def generate_funcs(function_template, variables):
    variable_names, variable_values = zip(*variables.items())
    for this_values in zip(*variable_values, strict=True):
        replacements = dict(zip(variable_names, this_values))
        func = function_template
        for var, value in replacements.items():
            func = func.replace(var, str(value) if isinstance(value, int) else f"({value})")
        yield func


out_names = ["B", "X1", "X2", "X3", "X4", "X5", "U1", "U2"]
vals = {
    "U1": [0] * 8,
    "U2": [0] * 8,
    "X1": [Fraction(4, 3), 1, 0, Fraction(2, 3), Fraction(-2, 3), Fraction(1, 3), Fraction(2, 3), Fraction(-1, 3)],
    "X2": [Fraction(10, 3), 0, 1, Fraction(2, 3), Fraction(1, 3), Fraction(-2, 3), Fraction(-1, 3), Fraction(2, 3)],
    "X3": [0] * 8,
    "X4": [0] * 8,
    "X5": [0] * 8,
}

z_function = "1 * X1 + 2 * X2 + 3 * X3 + 0 * X4 + 0 * X5 + M * U1 + M * U2"

results = {}

for out_name, operation in zip(out_names, generate_funcs(z_function, vals)):
    print(out_name, "=", operation)
    results[out_name] = sp.sympify(operation)

print(results)

The output is e.g.

B = 1 * (4/3) + 2 * (10/3) + 3 * 0 + 0 * 0 + 0 * 0 + M * 0 + M * 0
X1 = 1 * 1 + 2 * 0 + 3 * 0 + 0 * 0 + 0 * 0 + M * 0 + M * 0
X2 = 1 * 0 + 2 * 1 + 3 * 0 + 0 * 0 + 0 * 0 + M * 0 + M * 0
X3 = 1 * (2/3) + 2 * (2/3) + 3 * 0 + 0 * 0 + 0 * 0 + M * 0 + M * 0
X4 = 1 * (-2/3) + 2 * (1/3) + 3 * 0 + 0 * 0 + 0 * 0 + M * 0 + M * 0
X5 = 1 * (1/3) + 2 * (-2/3) + 3 * 0 + 0 * 0 + 0 * 0 + M * 0 + M * 0
U1 = 1 * (2/3) + 2 * (-1/3) + 3 * 0 + 0 * 0 + 0 * 0 + M * 0 + M * 0
U2 = 1 * (-1/3) + 2 * (2/3) + 3 * 0 + 0 * 0 + 0 * 0 + M * 0 + M * 0
{'B': 8, 'X1': 1, 'X2': 2, 'X3': 2, 'X4': 0, 'X5': -1, 'U1': 0, 'U2': 1}