implicit none
INTEGER, PARAMETER :: SP=SELECTED_REAL_KIND(6,37)
INTEGER, PARAMETER :: DP=SELECTED_REAL_KIND(15,307)
REAL(KIND=SP) :: begin, B, C
REAL(KIND=DP) :: V, accum, summa, eps_fourth
integer, parameter :: r = 1000
real(KIND=SP), dimension(1:r) :: myarray
integer, dimension(1:r) :: ticks
integer i, mult
! main control
v = epsilon(begin)
print *, "epsilon is ", v
eps_fourth = v/ 10.0_dp
begin = .49999_sp
myarray(1) = begin
ticks(1) = 45
do i =2,r
! inner control
mult = 0
accum = 0.0_dp
b = myarray ( i - 1)
c = myarray ( i - 1)
summa =0.0_dp
do while ( b == c )
accum = mult * eps_fourth
summa = accum + real(b, kind=dp)
c = real(summa, kind=sp)
mult = mult + 1
end do
myarray(i) = c
ticks(i) = mult
end do
do i = 1, r
print *, "i= ", i, myarray(i), "mult is ", ticks(i)
end do
end program
! gfortran real3.f90 -Wall -o out
! ./out >text1.txt
Abridged output:
i= 331 0.49999982 mult is 3
i= 332 0.49999985 mult is 3
i= 333 0.49999988 mult is 3
i= 334 0.49999991 mult is 3
i= 335 0.49999994 mult is 3
i= 336 0.49999997 mult is 3
i= 337 0.50000000 mult is 3
i= 338 0.50000006 mult is 4
i= 339 0.50000012 mult is 4
i= 340 0.50000018 mult is 4
i= 341 0.50000024 mult is 4
i= 342 0.50000030 mult is 4
i= 343 0.50000036 mult is 4
This is suggestive of an asymmetry that Ron Shephard and steve mentioned.
I'm curious that these values seem to be spaced at epsilon over something
in the range of [2.5-3.33).
I was elected to church council today (They must need warm bodies).
Since I became an ecclesiastical authority, and I've been wrestling with
what title I would be known by. I always like Mencken's quip about the
Archbishop: A Christian ecclesiastic of a rank superior to that attained
by Christ. I've decided to become
Your Frankness,
--
frank
"Guns: yes, they are harmful."
program hj
integer(8) i
real x
x = 0
i = 0
! Count number of floats in [0,0.5)
do
i = i + 1
x = nearest(x, 1.)
if (x >= 0.5) exit
end do
print '(A,I0)', '[0,0.5) = ', i
x = 0.5
i = 0
! count number of floats in [0.5,1)
do
i = i + 1
x = nearest(x, 1.)
if (x >= 1.) exit
end do
print '(A,I0)', '[0.5,1) = ', i
x = 0.5
do i = 1, 10
x = nearest(x, -1.)
end do
do i = 1, 20
print *, x, spacing(x)
x = nearest(x, 1.)
end do
end program hj
REMOVE:kargl[209] gfc4x -o z -O hj.f90
REMOVE:kargl[210] ./z
[0,0.5) = 1056964608
[0.5,1) = 8388608
0.49999970 2.98023224E-08
0.49999973 2.98023224E-08
0.49999976 2.98023224E-08
0.49999979 2.98023224E-08
0.49999982 2.98023224E-08
0.49999985 2.98023224E-08
0.49999988 2.98023224E-08
0.49999991 2.98023224E-08
0.49999994 2.98023224E-08
0.49999997 2.98023224E-08
0.50000000 5.96046448E-08
0.50000006 5.96046448E-08
0.50000012 5.96046448E-08
0.50000018 5.96046448E-08
0.50000024 5.96046448E-08
0.50000030 5.96046448E-08
0.50000036 5.96046448E-08
0.50000042 5.96046448E-08
0.50000048 5.96046448E-08
0.50000054 5.96046448E-08
Note, the above counting does not include subnormal numbers.
--
steve
[code and output elided]
> Note, the above counting does not include subnormal numbers.
I don't know what those are. We seem to have the same machine parameters:
dan@dan-desktop:~/source$ gfortran real5.f90 -Wall -Wextra -o out
dan@dan-desktop:~/source$ ./out
epsilon is 1.19209290E-07
[0,0.5) = 1056964608
[0.5,1) = 8388608
0.49999970 2.98023224E-08 z is 0.25000000
0.49999973 2.98023224E-08 z is 0.25000000
0.49999976 2.98023224E-08 z is 0.25000000
0.49999979 2.98023224E-08 z is 0.25000000
0.49999982 2.98023224E-08 z is 0.25000000
0.49999985 2.98023224E-08 z is 0.25000000
0.49999988 2.98023224E-08 z is 0.25000000
0.49999991 2.98023224E-08 z is 0.25000000
0.49999994 2.98023224E-08 z is 0.25000000
0.49999997 2.98023224E-08 z is 0.25000000
0.50000000 5.96046448E-08 z is 0.50000000
0.50000006 5.96046448E-08 z is 0.50000000
0.50000012 5.96046448E-08 z is 0.50000000
0.50000018 5.96046448E-08 z is 0.50000000
0.50000024 5.96046448E-08 z is 0.50000000
0.50000030 5.96046448E-08 z is 0.50000000
0.50000036 5.96046448E-08 z is 0.50000000
0.50000042 5.96046448E-08 z is 0.50000000
0.50000048 5.96046448E-08 z is 0.50000000
0.50000054 5.96046448E-08 z is 0.50000000
dan@dan-desktop:~/source$ cat real5.c
cat: real5.c: No such file or directory
dan@dan-desktop:~/source$ cat real5.f90
program hj
implicit none
integer(8) i
real x, y, z
y = epsilon(x)
print *, "epsilon is ", y
x = 0
i = 0
! Count number of floats in [0,0.5)
do
i = i + 1
x = nearest(x, 1.)
if (x >= 0.5) exit
end do
print '(A,I0)', '[0,0.5) = ', i
x = 0.5
i = 0
! count number of floats in [0.5,1)
do
i = i + 1
x = nearest(x, 1.)
if (x >= 1.) exit
end do
print '(A,I0)', '[0.5,1) = ', i
x = 0.5
do i = 1, 10
x = nearest(x, -1.)
end do
do i = 1, 20
z = spacing(x)/y
print *, x, spacing(x), "z is ", z
x = nearest(x, 1.)
end do
end program hj
! gfortran real5.f90 -Wall -Wextra -o out
! ./out >text1.txt
dan@dan-desktop:~/source$
This program shows the virtues of the nearest function, whose second
argument is of the kind of the first and nonzero, which then returns the
closest value of that kind in the direction of the second argument.
Spacing is an intrinsic that does what I tried to calculate!
> On Sun, 22 Nov 2009 15:52:29 -0800, steve wrote:
>
> [code and output elided]
>> Note, the above counting does not include subnormal numbers.
>
> I don't know what those are.
Floating point numbers whose significand has less than full precision.
They used to be called denormals. Using them avoids sudden underflow.
> We seem to have the same machine parameters:
>
> dan@dan-desktop:~/source$ gfortran real5.f90 -Wall -Wextra -o out
> dan@dan-desktop:~/source$ ./out
> epsilon is 1.19209290E-07
> [0,0.5) = 1056964608
which is equal to 126 X 2^23. There are seven bits for the biased exponent;
2^7=128, but then your numbers don't include +Inf and NaN.
> [0.5,1) = 8388608
which is exactly 2^23. The significand takes all possible values, while the
sign and biased exponent are fixed (because the range is 2^-1 to 2^0).
<--CUT-->
-- mecej4
I was expecting it to split the difference:
dan@dan-desktop:~/source$ gfortran real6.f90 -Wall -Wextra -o out
dan@dan-desktop:~/source$ ./out
0.50000000 5.96046448E-08
dan@dan-desktop:~/source$ cat real6.f90
real x
x = .5
print *, x, spacing(x)
end program
! gfortran real6.f90 -Wall -Wextra -o out
! ./out >text1.txt
dan@dan-desktop:~/source$
[snip denormals]
>> [0,0.5) = 1056964608
>
> which is equal to 126 X 2^23. There are seven bits for the biased
> exponent; 2^7=128, but then your numbers don't include +Inf and NaN.
>
>> [0.5,1) = 8388608
> which is exactly 2^23. The significand takes all possible values, while
> the sign and biased exponent are fixed (because the range is 2^-1 to
> 2^0).
I read up on the model for real data when I last read fortran reference,
but I didn't do anything with it. Does what I'm doing here have anything
to do with the above numbers?
dan@dan-desktop:~/source$ gfortran real7.f90 -Wall -Wextra -o out
dan@dan-desktop:~/source$ ./out
radix is 2
range is 37
dan@dan-desktop:~/source$ cat real7.f90
real x
integer d, e
x = .5
d = radix(x)
e = range(x)
print *, "radix is ", d
print *, "range is ", e
end program
! gfortran real7.f90 -Wall -Wextra -o out
! ./out >text1.txt
dan@dan-desktop:~/source$
> I was expecting it to split the difference:
>
> dan@dan-desktop:~/source$ gfortran real6.f90 -Wall -Wextra -o out
> dan@dan-desktop:~/source$ ./out
> 0.50000000 5.96046448E-08
> dan@dan-desktop:~/source$ cat real6.f90
> real x
> x = .5
> print *, x, spacing(x)
> end program
Is that the right value to return for SPACING() for an argument of
0.5? From the specification of the function, I would have thought
it should return the smaller number, 2.98E-08.
I would not expect the average value to be returned, or split the
difference, or anything like that. I would expect the difference to
the nearest number to be returned, which is 2.98E-08, not 5.96E-08.
I can think of several reasons why I might prefer the larger of the
two differences to be returned, not the smaller, but that's not the
way the function is defined.
$.02 -Ron Shepard
In F95, spacing(x) = b**(e - p) with the radix b, precision p and e
exponent of x. For x = 0.5 and gfortran's default real kind, one has
b = 2, e = 0, p = 24, which gives 5.96xxxe-08. Note, F2003 has a
different but irrelevant definition of spacing().
--
steve
> In F95, spacing(x) = b**(e - p) with the radix b, precision p and e
> exponent of x. For x = 0.5 and gfortran's default real kind, one has b
> = 2, e = 0, p = 24, which gives 5.96xxxe-08. Note, F2003 has a
> different but irrelevant definition of spacing().
What do you mean with "precision" here?
dan@dan-desktop:~/source$ gfortran real8.f90 -Wall -Wextra -o out
dan@dan-desktop:~/source$ ./out
for sp :
radix is 2
range is 37
p is 24.000000
for dp :
radix is 2
range is 307
q is 52.999999854364731
dan@dan-desktop:~/source$ cat real8.f90
implicit integer(a-t)
INTEGER, PARAMETER :: SP=SELECTED_REAL_KIND(6,37)
INTEGER, PARAMETER :: DP=SELECTED_REAL_KIND(15,307)
REAL(KIND=SP) :: x, p
REAL(KIND=DP) :: y, q
x = .5_sp
y = .5_dp
e = 0
d = radix(x)
f = range(x)
print *, " for sp :"
print *, "radix is ", d
print *, "range is ", f
p = e - log(spacing(x))/(log(real(d)))
print *, "p is ", p
print *, " for dp :"
g = radix(y)
h = range(y)
print *, "radix is ", g
print *, "range is ", h
q = e - log(spacing(y))/(log(real(g)))
print *, "q is ", q
end program
! gfortran real8.f90 -Wall -Wextra -o out
dan@dan-desktop:~/source$
It looks like precision is 53 for dp on my machine.