rwpenney / spfpm

Package for performing fixed-point, arbitrary-precision arithmetic in Python.
Other
64 stars 14 forks source link

Division Bug with Negative Divisor #14

Open WillGreen opened 1 year ago

WillGreen commented 1 year ago

Hello,

I've recently discovered your module and find it useful in hardware modelling. Thank you.

I think I've found a bug when dividing by a negative number. The answer seems to be off by one bit.

Dividing a negative number by a positive number produces the expected result.

I'm running SPFPM 1.6 on macOS 12.6.2 (arm64) with Python 3.10.9.

I've created a test script to show what I'm seeing:

from FixedPoint import FXfamily, FXnum

fp_family = FXfamily(n_bits=4, n_intbits=4)  # happens for other precisions too

x = fp_family(-1.0)
y = fp_family(1.0)
z = x / y
print("-1.0/1.0  ", z.toBinaryString(), " ", end="")
print(z.toDecimalString(precision=fp_family.fraction_bits))

x = fp_family(1.0)
y = fp_family(-1.0)
z = x / y
print("1.0/-1.0  ", z.toBinaryString(), " ", end="")
print(z.toDecimalString(precision=fp_family.fraction_bits))

x = fp_family(-1.0)
y = fp_family(-1.0)
z = x / y
print("-1.0/-1.0 ", z.toBinaryString(), " ", end="")
print(z.toDecimalString(precision=fp_family.fraction_bits))

x = fp_family(-3.0)
y = fp_family(2.0)
z = x / y
print("-3.0/2.0  ", z.toBinaryString(), " ", end="")
print(z.toDecimalString(precision=fp_family.fraction_bits))

x = fp_family(3.0)
y = fp_family(-2.0)
z = x / y
print("3.0/-2.0  ", z.toBinaryString(), " ", end="")
print(z.toDecimalString(precision=fp_family.fraction_bits))

x = fp_family(-3.0)
y = fp_family(-2.0)
z = x / y
print("-3.0/-2.0 ", z.toBinaryString(), " ", end="")
print(z.toDecimalString(precision=fp_family.fraction_bits))

x = fp_family(-7.5)
y = fp_family(2.5)
z = x / y
print("-7.5/2.5  ", z.toBinaryString(), " ", end="")
print(z.toDecimalString(precision=fp_family.fraction_bits))

x = fp_family(7.5)
y = fp_family(-2.5)
z = x / y
print("7.5/-2.5  ", z.toBinaryString(), " ", end="")
print(z.toDecimalString(precision=fp_family.fraction_bits))

x = fp_family(-7.5)
y = fp_family(-2.5)
z = x / y
print("-7.5/-2.5 ", z.toBinaryString(), " ", end="")
print(z.toDecimalString(precision=fp_family.fraction_bits))

Output

$ python3 spfpm-test.py
-1.0/1.0   1111.0000  -1
1.0/-1.0   1110.1111  -1.0625
-1.0/-1.0  0000.1111  0.9375
-3.0/2.0   1110.1000  -1.5
3.0/-2.0   1110.0111  -1.5625
-3.0/-2.0  0001.0111  1.4375
-7.5/2.5   1101.0000  -3
7.5/-2.5   1100.1111  -3.0625
-7.5/-2.5  0010.1111  2.9375
WillGreen commented 1 year ago

I've spotted another anomaly with division:

from FixedPoint import FXfamily, FXnum
fp_family = FXfamily(n_bits=4, n_intbits=5)  # five integer bits as answer is >8

x = fp_family(1.375)
y = fp_family(0.125)
z = x / y
print("1.375/0.125  ", z.toBinaryString(), " ", end="")
print(z.toDecimalString(precision=fp_family.fraction_bits))
1.375/0.125   01011.0100  11.25

If I do the division with Python numbers, then create a fixed-point number the answer is as expected:

z = fp_family(1.375 / 0.125)
1.375/0.125   01011.0000  11
rwpenney commented 1 year ago

Thanks, WIll, for these examples.

I'm not sure to what extent these are bugs rather than just symptoms of the fact that finite-precision arithmetic is always going to produce artefacts. Fixed-point arithmetic, especially when you have very limited numbers of fractional bits will make these artefacts much more significant. There is also a subtlety when constraining the whole-number bits (n_intbits) that spfpm loses the ability to preserve an implicit "sign bit" which may be the reason for the rounding differences between -1/1 versus 1/-1. This is a big difference from IEEE-754 floating-point arithmetic which besides being much more sophisticated overall, explicitly includes a sign-bit rather than working with twos-complement arithmetic.

There are points within spfpm where rounding necessarily has to occur, and a half-bit rounding term is included to try to bias things towards a better outcome on average. Perhaps that is something I need to take another look at to see if we can make some of the more common division operations slightly better behaved.

WillGreen commented 1 year ago

Thank you for the quick response. I'm using Python and your module to validate the fixed-point division logic in my FPGA design. I'm worked around this issue by doing the division in Python, then converting to fixed-point with spfpm:

model_q = fp_family(x / y)

I'd still be interested to know how rounding is occurring when spfpm does division. Is the internal representation two's complement? It looks like spfpm might be too keen to round towards negative infinity? But it's not causing me any problems, so feel free to close this issue if you like.

Verilog Design

In case it helps, this is what I've done in Verilog where I "round up half". I round up if the next bit (beyond the fixed-point range) is 1. For example, when I divide 13/7 with signed numbers and five integer and four fractional bits (a range of -16 to +15.9375):

z = x / y

x = 01101.0000  // 13
y = 00111.0000  // 7

Intermediate answer (with extra fractional bit calculated for rounding): 00001.11011

The final (fifth) fractional bit is 1, so in Verilog round up to:

z = 00001.1110 // 1.875

If I perform the same calculation with spfpm I get 13/7 00001.1101 1.8125

A more exact answer on a calculator is 1.857142857142857142857...

rwpenney commented 1 year ago

Thanks for the clarification. It seems like there's a choice to be made here between finding the closest 4-bit representation of 13/7 and what is the effect of trying the divide a 4-bit representation of 13 by a 4-bit representation of 7. Obviously, they aren't guaranteed to be the same, and even if they happen to be equal this doesn't mean that similar division operations will agree when implemented in fixed-point arithmetic.

It's interesting that your Verilog approach already uses an extra (fifth) fractional bit. Doing something similar in spfpm:

fam5 = FXfamily(5, 5)
fam4 = FXfamily(4, 5)
str(fam4(13) / fam4(7))
str(fam5(13) / fam5(7))

produces 1.812 (4 fractional bits) and 1.843 (5 fractional bits). If you rounded (fam5(13) / fam5(7)).toBinaryString() you'd get your Verilog-rounded 4-bit result.

It also looks like we're missing some convenient ways of doing this rounding within FXnum - I'll see if I can tweak the API accordingly.

WillGreen commented 1 year ago

I've settled on Gaussian rounding (round to even). It's slightly more complex to implement, but matches the default in Python and IEEE 745 (floating point standard).

You can see my implementation in Verilog for division here: https://github.com/projf/projf-explore/blob/main/lib/maths/div.sv - though this probably isn't particularly helpful if you don't know Verilog.

Thanks again for your module. I've used it in my multiplication maths test too: test/mul.py.