• #### Spacing of real numbers around zero

From Beliavsky@21:1/5 to All on Tue Apr 19 08:56:32 2022
Spacing(0.0_wp) is supposed to give the number adjacent to zero, but it looks like that is actually given by epsilon(x) * tiny(x) , which is smaller. The output of the program

program main
use iso_fortran_env, only: int32, real32
implicit none
integer, parameter :: ikind = int32, wp = real32
integer (kind=ikind) :: i, nbits
integer, parameter :: ni = 4
real(kind=wp) :: x
nbits = bit_size(i)
print*,"bit_size(i) =",nbits
print*,"epsilon(x) =",epsilon(x)
print*,"tiny(x) =",tiny(x)
print*,"epsilon(x) * tiny(x) =",epsilon(x)*tiny(x)
print*,"spacing(0.0_wp) =",spacing(0.0_wp)
print "(/,a12,a16)","i","transfer(i,x)"
do i=0,ni
print*,i,transfer(i,x)
end do
end program main

with gfortran and ifort is

bit_size(i) = 32
epsilon(x) = 1.19209290E-07
tiny(x) = 1.17549435E-38
epsilon(x) * tiny(x) = 1.40129846E-45
spacing(0.0_wp) = 1.17549435E-38

i transfer(i,x)
0 0.00000000
1 1.40129846E-45
2 2.80259693E-45
3 4.20389539E-45
4 5.60519386E-45

and the output of the 64-bit program

program main
use iso_fortran_env, only: int64, real64
implicit none
integer, parameter :: ikind = int64, wp = real64
integer (kind=ikind) :: i, nbits
integer, parameter :: ni = 4
real(kind=wp) :: x
nbits = bit_size(i)
print*,"bit_size(i) =",nbits
print*,"epsilon(x) =",epsilon(x)
print*,"tiny(x) =",tiny(x)
print*,"epsilon(x) * tiny(x) =",epsilon(x)*tiny(x)
print*,"spacing(0.0_wp) =",spacing(0.0_wp)
print "(/,a21,a16)","i","transfer(i,x)"
do i=0,ni
print*,i,transfer(i,x)
end do
end program main

is

bit_size(i) = 64
epsilon(x) = 2.2204460492503131E-016
tiny(x) = 2.2250738585072014E-308
epsilon(x) * tiny(x) = 4.9406564584124654E-324
spacing(0.0_wp) = 2.2250738585072014E-308

i transfer(i,x)
0 0.0000000000000000
1 4.9406564584124654E-324
2 9.8813129168249309E-324
3 1.4821969375237396E-323
4 1.9762625833649862E-323

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From steve kargl@21:1/5 to Beliavsky on Tue Apr 19 19:10:03 2022
Beliavsky wrote:

Spacing(0.0_wp) is supposed to give the number adjacent to zero,
but it looks like that is actually given by epsilon(x) * tiny(x) , which
is smaller. The output of the program

If your processor supports subnormal number, first nonzero positive
closest to 0 is 2**(-p-emin).

% cat z.f90
program foo
print *, 2.**(minexponent(1.) - digits(1.))
end

% gfcx -o z a.f90 && ./z
1.40129846E-45

The IEEE facilities of Fortran allow to query about subnormal numbers.

--
steve

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)
• From gah4@21:1/5 to steve kargl on Tue Apr 19 12:59:30 2022
On Tuesday, April 19, 2022 at 12:10:06 PM UTC-7, steve kargl wrote:
Beliavsky wrote:

Spacing(0.0_wp) is supposed to give the number adjacent to zero,
but it looks like that is actually given by epsilon(x) * tiny(x) , which
is smaller.

If your processor supports subnormal number, first nonzero positive
closest to 0 is 2**(-p-emin).

Which I never thought was a good idea. You get the equivalent
of about 0.04 extra exponent bits, for a lot of extra work
and/or hardware.

Note that the actual smallest number has zero significant bits.

It especially complicates parallel processors, like vector processors,
that try to do about the same amount of work for each one.

Even more, some require software help to get the right result.

If someone really sees the need, give it another whole exponent bit.

--- SoupGate-Win32 v1.05
* Origin: fsxNet Usenet Gateway (21:1/5)