Say I have two gmpy2.mpc objects x and y precise to p bits. When I compute x+y it may be that some parts of x and y cancel so I'm left with less precision.
For example
from gmpy2 import *
x = mpc('-0.55555')
y = mpc('0.5555500000001')
print(x+y)
The result is precise to only 4 significant figures even though x and y were precise to ~15.
I think I need to work out how many bits of cancellation occur when I do addition and subtraction and then take this away from the minimum of x or y's precision. For multiplication and division I think I will only lose 1 bit of precision at most.
So the question is quite general: how can I keep track of the precision of mpc objects, particularly when adding and subtracting them?
The following function will return the numbers of matching bits of two mpfr objects.
import gmpy2
def matching_bits(x, y):
'''Returns the number of bits that match between x and y. The
sign of x and y are ignored. x and y must be of type mpfr.'''
# Force both values to be positive, and x >= y.
x = abs(x)
y = abs(y)
if x < y:
x, y = y, x
if not isinstance(x, type(gmpy2.mpfr(0))) or not isinstance(y, type(gmpy2.mpfr(0))):
raise TypeError("Arguments must be of type 'mpfr'.")
x_bits, x_exp, x_prec = x.digits(2)
y_bits, y_exp, y_prec = y.digits(2)
# (x_exp - y_exp) is the number of zeros that must be prepended
# to x to align the mantissas. If that is greater than the precision
# y, then no bits in common.
if (x_exp - y_exp) > x_prec:
return 0
x_bits = "0" * (x_exp - y_exp) + x_bits
count = 0
while count < min(x_prec, y_prec) and x_bits[count] == y_bits[count]:
count += 1
return count
I haven't extensively tested this function but it should give you a start. You will need to check the real and imaginary components separately. You'll probably want to modify it to check for sign and whether you are performing addition vs. subtraction.