[Numeric.Interval_Lib] typecasting behaviour
Hi all, Background... I'm using the interval arithmetic library to evaluate the accuracy of different algorithms for my application (i.e. I'm using the library "offline"; once a winner is chosen, the interval arithmetic will be removed). One possibility that I am pursuing is using mixed precision arithmetic; i.e. do some parts in double precision, some parts in single precision - this is to explore the possibility of eventually utilising non-conventional processing hardware such as GPUs, which tend to have only single precision arithmetic hardware. The query... Before I started using the library, I had this assumption: let A be a float with an exact representation (i.e. width=0). If we cast A to a double, the resulting double will have a non-zero width - specifically the width will be equal to the extra precision afforded to doubles vs floats. Unless I've messed something up, I find that the library does not actually behave this way, and to be frank, I think it ought to! :-) Here's the simple testcase: typedef interval<float> I1; typedef interval<double> I2; I1 y = 1.0; I2 yy = (I1)y; cout << yy.upper() - yy.lower(); I expected yy.upper()-yy.lower() to be non-zero, but it was zero. I had a look at the policies but couldn't spot anything obviously relevant to this situation. Any thoughts on this behaviour per se, or any ideas on a neat and elegant way to implement my desired behaviour? FWIW, my system is PowerPC G4 with powerpc-apple-darwin8-g++-4.0.0 (GCC) 4.0.0 20041026 (Apple Computer, Inc. build 4061). -tirath
Le jeudi 19 octobre 2006 à 01:03 +0800, tirath a écrit :
Hi all,
The query... Before I started using the library, I had this assumption: let A be a float with an exact representation (i.e. width=0). If we cast A to a double, the resulting double will have a non-zero width - specifically the width will be equal to the extra precision afforded to doubles vs floats.
If you have a set of numbers that is exactly representable as a singleton interval<float>, then this set of numbers is also representable as a singleton interval<double>. There is no reason to use a bigger interval than necessary to represent this set of numbers.
Unless I've messed something up, I find that the library does not actually behave this way, and to be frank, I think it ought to! :-) Here's the simple testcase:
typedef interval<float> I1; typedef interval<double> I2; I1 y = 1.0; I2 yy = (I1)y; cout << yy.upper() - yy.lower();
I expected yy.upper()-yy.lower() to be non-zero, but it was zero.
This is the expected behavior. Interval "y" is a set with a single element 1 and you assign it to interval "yy". So interval "yy" has to contain 1 and it doesn't have to contain any other value. The best interval that satisfies this property is the interval [1,1], and this is the one the library computes.
I had a look at the policies but couldn't spot anything obviously relevant to this situation. Any thoughts on this behaviour per se, or any ideas on a neat and elegant way to implement my desired behaviour?
Could you detail why you desire this behavior? I can't think of a situation where you would like the intervals to get worse when you convert from interval<float> to interval<double>. Intervals may get worse when you convert from interval<double> to interval<float>, but not the other way around. Even for the application you described, I don't understand why you expect this behavior. A conversion from float to double is exact, so you want a conversion from interval<float> to interval<double> to be exact too, in order for your simulations to be useful. Best regards, Guillaume
On 19/10/2006, at 4:40 PM, Guillaume Melquiond wrote:
Le jeudi 19 octobre 2006 à 01:03 +0800, tirath a écrit :
Hi all,
The query... Before I started using the library, I had this assumption: let A be a float with an exact representation (i.e. width=0). If we cast A to a double, the resulting double will have a non-zero width - specifically the width will be equal to the extra precision afforded to doubles vs floats.
If you have a set of numbers that is exactly representable as a singleton interval<float>, then this set of numbers is also representable as a singleton interval<double>. There is no reason to use a bigger interval than necessary to represent this set of numbers.
Thanks for your response Guillaume. I think my explanation was a bit vague. Take two (with a more concrete example)... In general, when you have a 32-bit floating point representation of a real number, the precision is less than it would be with a 64-bit floating point representation of the number. Therefore when "promoting" a float to a double, say: interval<double> x; interval<float> y; ... y=x; shouldn't we (in the spirit of error propagation) acknowledge that the data assigned to y was of inferior precision? A range too small to be expressed in single precision may be large enough to be expressed in double precision, as single precision ulp ~ 1E-7 wheras double precision ulp ~ 1E-16. Consider the following... interval<double> x; interval<float> y; interval<double> z; x=1.00000000001; //the following loop forces error accumulation for(int i=0; i<MAX_ITER; i++) x*=x; print_width(x); y=x; print_width(y); z=y; print_width(z); with MAX_ITER=30, the error accumulated on x is pretty large, and consequently the output is: (variable: {lower} {upper} {width}) x: 1.0108 1.0108 2.40992e-07 y: 1.0108 1.0108 2.38419e-07 z: 1.0108 1.0108 2.40992e-07 with MAX_ITER=5, the error accumulated on x is small, and consequently: x: 1 1 6.88338e-15 y: 1 1 0 z: 1 1 0 I think this behaviour is a bit misleading because it seems to imply that assigning the double into the float then to another double actually zeros the interval width! The current situation implies that one can sanitise all double precision results by casting it to a float; this won't remove large intervals, but it will remove small intervals, and with iterative algorithms this can be very significant. On the contrary, shouldn't there be a precision cost associated with "promoting" a float to a double, since the source data is inferior to the precision of a double? Garbage in garbage out, no? Best regards, Tirath Ramdas
Le vendredi 20 octobre 2006 à 16:30 +0800, tirath a écrit :
In general, when you have a 32-bit floating point representation of a real number, the precision is less than it would be with a 64-bit floating point representation of the number. Therefore when "promoting" a float to a double, say: interval<double> x; interval<float> y; ... y=x;
shouldn't we (in the spirit of error propagation) acknowledge that the data assigned to y was of inferior precision? A range too small to be expressed in single precision may be large enough to be expressed in double precision, as single precision ulp ~ 1E-7 wheras double precision ulp ~ 1E-16.
This is not correct. If the range is too small to be represented in single precision, it will be inflated until it can be represented. So you don't have to inflate the range one more time when you convert it back to double precision, as it already was. But thanks to your testcase, I now understand where the confusion comes from. There must be a bug somewhere, as the width of the single precision interval should never have been zero. For now, I haven't found anything obvious. I'm not even sure whether the bug is in the interval library or not. I will have to investigate. Best regards, Guillaume
Le samedi 21 octobre 2006 à 16:06 +0200, Guillaume Melquiond a écrit :
But thanks to your testcase, I now understand where the confusion comes from. There must be a bug somewhere, as the width of the single precision interval should never have been zero. For now, I haven't found anything obvious. I'm not even sure whether the bug is in the interval library or not. I will have to investigate.
Following up on myself. There is a bug [*] in GCC 3.4 and later such that the compiler generates wrong code when converting from interval<double> to interval<float>. So, if your compiler is GCC, the behavior you observe is not surprising (although completely wrong). Best regards, Guillaume [*] http://gcc.gnu.org/bugzilla/show_bug.cgi?id=21032
On 21/10/2006, at 10:06 PM, Guillaume Melquiond wrote:
This is not correct. If the range is too small to be represented in single precision, it will be inflated until it can be represented. So you don't have to inflate the range one more time when you convert it back to double precision, as it already was.
Great, that works as I had hoped. On 22/10/2006, at 9:31 PM, Guillaume Melquiond wrote:
Following up on myself. There is a bug [*] in GCC 3.4 and later such that the compiler generates wrong code when converting from interval<double> to interval<float>. So, if your compiler is GCC, the behavior you observe is not surprising (although completely wrong).
You are correct; with gcc 4.0 the testcase fails as previously described, with gcc 3.3 it works as expected. Thanks for your help troubleshooting this. Sincerely, Tirath Ramdas
tirath wrote:
The query... Before I started using the library, I had this assumption: let A be a float with an exact representation (i.e. width=0). If we cast A to a double, the resulting double will have a non-zero width - specifically the width will be equal to the extra precision afforded to doubles vs floats.
I'm not familiar with the interval lib, but if a float is an exact value, then promoting it to a wider type will result in an exact value as well. I assume what you're looking for is a way to say "here is a float, that is holding the nearest representable value, but is inexact (has up to 0.5ulp error)" ? Then that error would propagate to the double as you expect. John.
participants (3)
-
Guillaume Melquiond -
John Maddock -
tirath