Numerics question: reading floating point constants

> This post is one for the numerics fraction. I was just in the  
> process of transcribing some interesting numerical code from C to  
> Squeak while I was running into an issue with floating point  
> constants. The code that I have has constants like here:
>
>   pio2_2  =  6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
>
> but naively transcribing this into, e.g.,
>
>   pio22 := 6.07710050630396597660e-11.
>
> turned out to get some of the expected results wrong. Why? Because  
> the bit pattern is different, e.g., what we get is this:
>
>   6.07710050630396597660e-11 hex => '3DD0B4611A5FFFFF'
>
> Notice how the last hex digits are different.
>
> My question here is the following: While I can understand why  
> *printing* might give different results based on the algorithm  
> chosen (e.g., what amount of error do we allow when printing  
> floating point numbers etc) is it to be expected that *reading* has  
> similar properties? In particular considering that, e.g.,
>
>   (607710050630396597660 * 1.0e-31) hex => '3DD0B4611A600000'
>
> gives the "correct" bit pattern makes me think that there might be  
> something broken in our way of reading floating point numbers.
>
> Obviously, in a case like mine the best way to work around this  
> problem is to construct the floats from the bit patterns directly  
> but I can't stop but wonder what the impact of a subtle bug like  
> this might be on code which does *not* include the hex denotations  
> for the relevant bit patterns.
>
> So. Is this a bug?

A tricky question, indeed.

>> This post is one for the numerics fraction. I was just in the  
>> process of
>> transcribing some interesting numerical code from C to Squeak while I
>> was running into an issue with floating point constants. The code  
>> that I
>> have has constants like here:
>>
>>    pio2_2  =  6.07710050630396597660e-11, /* 0x3DD0B461,  
>> 0x1A600000 */
>>
>> but naively transcribing this into, e.g.,
>>
>>    pio22 := 6.07710050630396597660e-11.
>>
>> turned out to get some of the expected results wrong. Why? Because  
>> the
>> bit pattern is different, e.g., what we get is this:
>>
>>    6.07710050630396597660e-11 hex => '3DD0B4611A5FFFFF'
>>
>> Notice how the last hex digits are different.
>>
>> My question here is the following: While I can understand why  
>> *printing*
>> might give different results based on the algorithm chosen (e.g.,  
>> what
>> amount of error do we allow when printing floating point numbers  
>> etc) is
>> it to be expected that *reading* has similar properties? In  
>> particular
>> considering that, e.g.,
>>
>>    (607710050630396597660 * 1.0e-31) hex => '3DD0B4611A600000'
>>
>> gives the "correct" bit pattern makes me think that there might be
>> something broken in our way of reading floating point numbers.
>>
> A tricky question, indeed.
> I do not predend to have the answers to all problems that floats
> can introduce into computer science, but I think that method
> Number class>>readRemainderOf:from:base:withSign:
> has some questionable features.
> Here is what happens when we read a constant of the form
> <intPart>.<fractionalPart>e<expPart> :
>  1. the <intPart> is read. This is a nonnegative integer
>  2. the <fractionalPart> ist read. This is also a nonnegative integer.
>  3. The<fractionalPart> is converted into a float. A floating point
> division follows to scale that value. The floating point division
> requires a second #asFloat to bring the integer value of
>    (base raisedTo: aStream position - fracpos)
> into a float.
> 4. next, the <intPart> is converted into a float (this is the third  
> #asFloat
> and added to the float representation of the <fractionalPart>.
> 5. Finally, the <expPart> is read. The scale factor
>   (base raisedTo: <expPart> is multiplied into the result of
> step (4), this multiplication reqires again an asFloat, the forth one
> in our algorithm.
>
> Well, I really think that four calls of  #asFlaot are too much.
> We introduce approximated number representations far too
> early and we do this without need.
>
> The line
>
>>    (607710050630396597660 * 1.0e-31) hex => '3DD0B4611A600000'
>
> of your post clearly shows that it is possible to introduce floats
> much later in the computation
> Attached you find a proposal (for 3.9 #6713) to avoid three of the  
> four
> convertions that we use now.
>
>> Obviously, in a case like mine the best way to work around this  
>> problem
>> is to construct the floats from the bit patterns directly but I can't
>> stop but wonder what the impact of a subtle bug like this might be on
>> code which does *not* include the hex denotations for the relevant  
>> bit
>> patterns.
>>
>> So. Is this a bug?
>>
>> Cheers,
>>    - Andreas
>>
>
> So much for now, in the afternoon I will
>  * read your comments
>  * look into that sutff again
>  * write some test cases
>
> Cheers, Boris
>
> As always, comments are very welcome
> <ReadFloatProposal.1.cs>
>

> Hi Folks,
>
> This post is one for the numerics fraction. I was just in the process of
> transcribing some interesting numerical code from C to Squeak while I
> was running into an issue with floating point constants. The code that I
> have has constants like here:
>
>   pio2_2  =  6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
>
> but naively transcribing this into, e.g.,
>
>   pio22 := 6.07710050630396597660e-11.
>
> turned out to get some of the expected results wrong. Why? Because the
> bit pattern is different, e.g., what we get is this:
>
>   6.07710050630396597660e-11 hex => '3DD0B4611A5FFFFF'
>
> Notice how the last hex digits are different.
>
> My question here is the following: While I can understand why *printing*
> might give different results based on the algorithm chosen (e.g., what
> amount of error do we allow when printing floating point numbers etc) is
> it to be expected that *reading* has similar properties? In particular
> considering that, e.g.,
>
>   (607710050630396597660 * 1.0e-31) hex => '3DD0B4611A600000'
>
> gives the "correct" bit pattern makes me think that there might be
> something broken in our way of reading floating point numbers.
>
> Obviously, in a case like mine the best way to work around this problem
> is to construct the floats from the bit patterns directly but I can't
> stop but wonder what the impact of a subtle bug like this might be on
> code which does *not* include the hex denotations for the relevant bit
> patterns.
>
> So. Is this a bug?
>
> Cheers,
>   - Andreas
>
>

> 1 point for Andres, conversion from base 10 to base 2 is inexact.
>
> 1 point for Boris : successive inexact operations are not only  
> costly, they
> can degrade accuracy.
>
> I wonder if the conversion shouldn't be made with respect to IEEE  
> rounding
> mode, and if i remember, this rounding mode is a parameter of  
> underlying
> hardware.
>
> In the past, i tried several algo. and finally, i lazily decided to  
> use a call
> to sprintf/sscan from within VW.
> But with a good reference such as Tony's (1 point) we should make  
> it in
> Smalltalk.
>
> Nicolas
>
>

> Do you know a good book or lectures on the problems related to  
> number representation in our nice
> computers?
>
> Stef
>
> On 9 avr. 06, at 18:18, nicolas cellier wrote:
>
>> 1 point for Andres, conversion from base 10 to base 2 is inexact.
>>
>> 1 point for Boris : successive inexact operations are not only  
>> costly, they
>> can degrade accuracy.
>>
>> I wonder if the conversion shouldn't be made with respect to IEEE  
>> rounding
>> mode, and if i remember, this rounding mode is a parameter of  
>> underlying
>> hardware.
>>
>> In the past, i tried several algo. and finally, i lazily decided  
>> to use a call
>> to sprintf/sscan from within VW.
>> But with a good reference such as Tony's (1 point) we should make  
>> it in
>> Smalltalk.
>>
>> Nicolas
>>
>>
>
>

>
> Stef
>
> On 9 avr. 06, at 18:18, nicolas cellier wrote:
> > 1 point for Andres, conversion from base 10 to base 2 is inexact.
> >
> > 1 point for Boris : successive inexact operations are not only
> > costly, they
> > can degrade accuracy.
> >
> > I wonder if the conversion shouldn't be made with respect to IEEE
> > rounding
> > mode, and if i remember, this rounding mode is a parameter of
> > underlying
> > hardware.
> >
> > In the past, i tried several algo. and finally, i lazily decided to
> > use a call
> > to sprintf/sscan from within VW.
> > But with a good reference such as Tony's (1 point) we should make
> > it in
> > Smalltalk.
> >
> > Nicolas