Thinking again, one of the best solutions would probably be to convert your numbers to
gmpy2.mpfr
instances on 128 bits. It would give you immediately all the arithmetic operations and more. The anyfloat module can be also used to help converting from/to
gmpy2.mpfr
.
Here is an example (not sufficiently tested, but it gives the idea)
from anyfloat import anyfloat
import gmpy2
def anyfloat_to_mpfr(a):
# XXX not sufficiently tested
g = gmpy2.mpfr(-a.significand if a.sign else a.significand)
p = 1 + a.exponent - a.significand.bit_length()
if p > 0:
g = gmpy2.mul_2exp(g, p)
elif p < 0:
g = gmpy2.div_2exp(g, -p)
return g
def anyfloat_from_mpfr(x):
# XXX not sufficiently tested
m, e = x.as_mantissa_exp()
m, e = abs(int(m)), int(e)
s = gmpy2.sign(x) == -1
return anyfloat(s, e + m.bit_length() - 1, m)
def main():
ctx = gmpy2.context(precision=113, subnormalize=True)
old = gmpy2.get_context()
gmpy2.set_context(ctx)
size = (15, 112) # the size of exponent and mantissa for 128 bits
encoded = '4000 921f b544 42d1 8469 898c c517 01b8'
a = anyfloat.from_ieee(int(encoded.replace(' ',''), 16), size)
g = anyfloat_to_mpfr(a)
print(repr(g))
aa = anyfloat_from_mpfr(g)
assert aa == a
print(hex(aa.to_ieee(size)))
gmpy2.set_context(old)
if __name__ == '__main__':
main()
Output:
mpfr('3.1415926535897932384626433832795028',113)
0x4000921fb54442d18469898cc51701b8
It seems to me that the mpfr precision must be set to 113 (it is normally 53 for double) for this to work in your quad floating point format. Also the use of subnormal=True in the gmpy2 context is an indication for the library to adhere to IEEE 754 rules. Now you're free to perform any calculation that you want with the numbers and convert them back later to IEEE.