[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
mpf: x * 1.0 != x
From: |
Peter Moulder |
Subject: |
mpf: x * 1.0 != x |
Date: |
Tue, 17 Jul 2001 00:37:19 +1000 |
Actually, gmp is behaving according to specification here; but x*1==x is
a desirable property, and easy to give (for mpf_mul, but not mpf_mul_ui
unfortunately), so I'm sending a patch.
A real-world situation where the problem arises is in performing a
"pivot" operation on matrix A at (i,j), to remove all elements from
column j other than in row i:
/* Scale row i to have 1 in column j. */
for(t = nColumns; t--;)
A[i][t] /= A[i][j];
/* Subtract multiples of row i from the other rows, so that A[k][j]
becomes zero for all k!=i. */
for(k = nRows; k--;)
if(k != i)
for(t = nColumns; t--;)
A[k][t] -= A[i][t] * A[k][j];
I hope you do apply this patch or similar for mpf_mul, since I'd guess
that this arises relatively often, and it works fine with floats
(so long as there are no NaN's or infinities), and the patch I'm sending
seems fairly clean.
I couldn't see a clean fix for mpf_mul_ui. But I can't off-hand think
of situations where one would rely on x - 1u * x == 0, so I'm going with
the clean code that still performs to spec.
Arguments relevant to whether to change mpf_mul and/or mpf_mul_ui:
- In both cases, the existing behavior conforms to documentation.
- For each case, making this change would also conform to
documentation.
- mpf_mul seems to have a clean fix that does not complicate the code.
- mpf_mul_ui appears not to have a no-readability-impact fix.
- It's desirable that x - 1 * x == 0, especially for the mpf_mul case.
Hardware float's and double's do have this property.
- A patch for mpf_mul is already written (attached).
- Making change to software is inherently dangerous; risk of
introducing bugs.
- It's ugly for mpf_mul and mpf_mul_ui to give different results
(different amount of precision) for numerically-equal arguments.
I attach a patch and some test case code. (Also a slightly cleaned-up
test case for mpf_div.)
pjm.
mpf-mul-patch
Description: Text document
mpf-mul-test.c
Description: Text document
mpf-div-test.c
Description: Text document
- mpf: x * 1.0 != x,
Peter Moulder <=