Re: Problem with averaging integers?

From: David Brown <dbrown_at_nyahnyahspammersnyahnyah>
Date: Tue Oct 18 2011 - 16:13:55 MDT

Hi Chuck,
This is ultimately a result of the fact that the integer 20060108 has too many digits to be represented precisely as a "float" value (32 bits in NCL). Only 6 digits of decimal precision can be guaranteed
to be correct with a 32 bit float. NCL's averaging function does its actual calculation in double precision, and gets the correct answer, but, as currently coded, unless the input is double, the output is converted to float before returning. You can get an accurate answer in this case only by converting the input to double, thereby ensuring a double return value.

ncl 5> tmp = new((/ 60 /), integer)
ncl 6> tmp = 20060109
ncl 7> print(avg(todouble(tmp)))

(0) 20060109

It might be argued that NCL should look at the size of the result before deciding whether to return a float or a double, or maybe it should return a double by default, but this result is in line with NCL's documented behavior.
 -dave

On Oct 17, 2011, at 8:53 AM, Charles Bardeen wrote:

> Hi Dennis,
>
> I am sorry I wasn't clear. I am not concerned about how the float is printed. I am concerned that the number ends in 108 rather than 109. It was set to be 20060109 and the average is 20060108. That appears to be true with your result too.
>
> Thanks,
>
> Chuck
>
> On Oct 17, 2011, at 8:18 AM, Dennis Shea <shea@ucar.edu> wrote:
>
>> There is nothing wrong. The 'print' value is
>> at 6 decimal values which is generally 'sufficient' for floats.
>> printing the full value ... gives the expected anser.
>>
>>
>> tmp = new((/ 60 /), integer)
>> tmp = 20060109 ; do not use tmp(:) [slower]
>> avg_tmp8 = avg(tmp)
>> print(avg_tmp8)
>> print(sprintf("%20.13f", avg_tmp8))
>>
>> ===
>> Variable: avg_tmp8
>> Type: float
>> Total Size: 4 bytes
>> 1 values
>> Number of Dimensions: 1
>> Dimensions and sizes: [1]
>> Coordinates:
>> Number Of Attributes: 1
>> _FillValue : -2.147484e+09
>> (0) 2.006011e+07
>>
>>
>> (0) 20060108.0000000000000 <=== sprintf
>>
>>
>> On 10/16/11 3:52 PM, Charles Bardeen wrote:
>>> The following code yields an average that is 1 from the true average …
>>>
>>>
>>> Copyright (C) 1995-2010 - All Rights Reserved
>>> University Corporation for Atmospheric Research
>>> NCAR Command Language Version 5.2.1
>>> The use of this software is governed by a License Agreement.
>>> See http://www.ncl.ucar.edu/ for more details.
>>> ncl 0> tmp = new((/ 60 /), integer)
>>> ncl 1> tmp(:) = 20060109
>>> ncl 2> print(avg(tmp))
>>> (0) 2.006011e+07
>>> ncl 3> print(doubletointeger(avg(tmp)))
>>> (0) 20060108
>>> ncl 4>
>>>
>>> Any idea why NCL is making this mistake. Is it an overflow or rounding
>>> problem? If I use smaller numbers, then it works correctly …
>>>
>>> ncl 8> tmp(:) = 2006
>>> ncl 9> print(avg(tmp))
>>> (0) 2006
>>> ncl 10> print(doubletointeger(avg(tmp)))
>>> (0) 2006
>>> ncl 11>
>>>
>>>
>>> Thanks,
>>>
>>> Chuck
>>>
>>>
>>> --------------------------------------------------------------
>>> Charles Bardeen
>>> National Center for Atmospheric Research
>>> P.O. Box 3000, Boulder, CO 80307-3000
>>> Phone: (303) 497-1752, Fax: (303) 497-1400
>>> bardeenc@ucar.edu <mailto:bardeenc@ucar.edu>
>>>
>>>
>>>
>>>
>>>
>>> _______________________________________________
>>> ncl-talk mailing list
>>> List instructions, subscriber options, unsubscribe:
>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Tue Oct 18 16:14:01 2011

This archive was generated by hypermail 2.1.8 : Wed Oct 19 2011 - 13:36:10 MDT