Re: WHERE problems (longish) [message #35907] |
Tue, 22 July 2003 09:45  |
Paul Van Delst[1]
Messages: 1157 Registered: April 2002
|
Senior Member |
|
|
David Fanning wrote:
>
> Benjamin Panter writes:
>
>> The values which have -1 certainly exist - and were generated in exactly
>> the same way as the others. I've put the array online if anyone fancies
>> looking at it - http://www.roe.ac.uk/~bdp/where_problem.idl
>>
>> Am I being stupid again? What is special about 2980,3000 and 3020??
>
> There is nothing special about *those* numbers, but those
> are not the numbers you are using in your WHERE statement.
> You are using 2980.0, 3000.0, and 3020.0. While there isn't
> a big difference between integers and floats to you, there
> is a HUGE difference to a computer. Better read these articles:
>
> http://www.dfanning.com/math_tips/sky_is_falling.html
> http://www.dfanning.com/math_tips/razoredge.html
Hmm. Reading that razor's edge article made me dig out a little f90 program I
wrote to determine the real number spacing.
Given some number, A, one can do:
exponent=ceil(alog(abs(A))/alog(2.0d))
radix=2.0d ; (machar(/double)).ibeta ??
epsilon=(machar(/double)).eps
spacing=epsilon * (radix^(exponent-1))
So for, say, A = 1.234568d+16
EXPONENT LONG = 54
SPACING DOUBLE = 2.0000000
For A = 1.234568d+01
EXPONENT LONG = 4
SPACING DOUBLE = 1.7763568e-15
For A = 1.234568d-01
EXPONENT LONG = -3
SPACING DOUBLE = 1.3877788e-17
and for A = 1.234568d-16
EXPONENT LONG = -52
SPACING DOUBLE = 2.4651903e-32
which agree with the outputs of my f90 code using the intrinsic functions
EXPONENT and SPACING. The problem is, of course, what to do when A = 0.0. I
would just use spacing = epsilon/2.0. I think that, in this case, doing
something like
two=2.0d
radix=two
IF ( A .eq. 0.0 ) THEN $
spacing = epsilon/two $
ELSE BEGIN
exponent=ceil(alog(abs(A))/alog(two))
spacing=epsilon * (radix^(exponent-1))
ENDELSE
where you are counting on the equality check .EQ. *not* to work for numbers that
aren't represented as exactly zero (since exact zero can be represented). But
I'm not sure. You may also need a check to limit the calculated exponent to the
range allowed (the minexp and maxexp fields of the output from machar).
ANYway.... back to work....
paulv
--
Paul van Delst
CIMSS @ NOAA/NCEP/EMC
Ph: (301)763-8000 x7748
Fax:(301)763-8545
|
|
|