! Test: NEAREST; EPSILON; spacing.
! X+, X-: xH =nearest(x,+); xL =nearest(x,-) 

real(4) x4, x, x4H, x4L, ulp_4, ulp4
real(8) x8, x8H, x8L, ulp_8, ulp8, eps8, u82
real(16) x16, x16H, x16L, ulp_16, ulp16, eps16, u162

character(260) fname, name*20, w1*1
integer kind;
character(2) kind_char

!ulp_4(x4) =2.0**(exponent(x4)-24)
!ulp_8(x8) =2.D0**(exponent(x8)-53)
!ulp_16(x16) =2.Q0**(exponent(x16)-113)

write(*, '(a,/)') 'NEAREST(x)'

100 continue    ! input kind
write(*, '(a\)') 'Input real kind (4; 8; 16): '
read(*,*) kind

if(.NOT. (kind ==4 .or. kind ==8 .or. kind ==16) ) then
kind =4
write(*, '(a,/, a,/)') '** Invalid "kind".', '** "kind" set to default: 4'
endif

write(*,*) 'Output File Name: "name-kind.rez"'
write(*,'(a,\)') 'Input "name" (<= 20 chars): '
read(*,*) name

ik =1
if(kind ==16) ik =2
write(kind_char, '(I<ik>)') kind
fname =trim(name)//'-'//trim(kind_char)//'.rez'
open(1,file=fname, status='unknown')

write(1,*) 'SPACING = spatierea absoluta a numerelor in vecinatatea lui x;'

!ik =1
!if(kind ==16) ik =2
write(1,'(/,a, I<ik>, a)') 'REAL(', kind, ')'

200 continue    ! input x
write(*,'(a,\)') 'Input "x": '

select case(kind)
case(4)
read(*,*) x4

case(8)
read(*,*) x8
goto 800

case(16)
read(*,*) x16
goto 1600

case default
read(*,*) x4
end select

write(1,*) 
write(1,*) 'x, x4:', x4

x4H =nearest(x4, 1.)
x4L =nearest(x4, -1.)
write(1, '(/,a, 1pE15.6)') ' nearest(x4, +)=', x4H
write(1, '(a, 1pE15.6)') ' nearest(x4, -)=', x4L
write(1,*) 'nearest(x4, +)-x4=', x4H-x4
write(1,*) 'x4 -nearest(x4, -)=', x4-x4L

write(1,*)
ulp4 =ulp_4(x4)
write(1,*) 'ulp(x) =', ulp4

if(x4 +ulp4 ==x4) then
write(1,*) 'x +ulp(x) ==x !'
else
!write(1,*) 'x +ulp(x) =', x4 +ulp4, ' > x'
write(1, '(a, 1pe14.7, a)') 'x +ulp(x) =', x4 +ulp4, ' > x'
endif

if( x4+ulp4 > x4H ) then
write(1,*) '*** EXCEPTION: x+ulp(x) > nearest(x4, +)'
endif

if(x4 +ulp4/2 ==x4) then
write(1,*) 'x +ulp(x)/2 ==x'
else
write(1,'(a, 1pe14.7, a)') 'x +ulp(x)/2 =', x4 +ulp4/2., ' > x !'
! ulp/4:
    if(x4 +ulp4/4 ==x4) then
    write(1,*) 'x +ulp(x)/4 ==x'
    else
    write(1,*) 'x +ulp(x)/4 =', x4 +ulp4/4., ' > x !!'
    endif
endif


write(1,*)
write(1,*) 'spacing(x4)=', spacing(x4)
goto 900

800 continue
write(1,*) 
write(1,*) 'x, x8:', x8

x8H =nearest(x8, 1.d0)
x8L =nearest(x8, -1.d0)
write(1, '(/,a, 1pE23.15)') ' nearest(x8, +)=', x8H
write(1, '(a, 1pE23.15)') ' nearest(x8, -)=', x8L
write(1,*) 'nearest(x8, +)-x8=', x8H-x8
write(1,*) 'x8 -nearest(x8, -)=', x8-x8L

ulp8 =ulp_8(x8)
write(1,*)
write(1,*) 'ulp(x) =', ulp8

if(x8 +ulp8 ==x8) then
write(1,*) 'x +ulp(x) ==x !'
else
write(1,'(a, 1pe23.15, a)') 'x +ulp(x) =', x8+ulp8, ' > x'
endif


if(x8+ulp8/2 ==x8) then
write(1,*) 'x +ulp(x)/2 ==x'
else
write(1,'(a, 1pe23.15, a)') 'x +ulp(x)/2 =', x8+ulp8/2., ' > x !'
! ulp/4:
    if(x8+ulp8/4 ==x8) then
    write(1,*) 'x +ulp(x)/4 ==x'
    else
    write(1,*) 'x +ulp(x)/4 =', x8+ulp8/4., ' > x !!'
    endif
endif


write(1,*)
write(1,*) 'spacing(x8)=', spacing(x8)
goto 900

1600 continue
write(1,*) 
write(1,*) 'x; x16:', x16

x16H =nearest(x16, 1.q0)
x16L =nearest(x16, -1.q0)
write(1, '(/,a, 1pE43.33)') ' nearest(x16, +)=', x16H
write(1, '(a, 1pE43.33)') ' nearest(x16, -)=', x16L
write(1,*) 'nearest(x16, +)-x16=', x16H-x16
write(1,*) 'x16 -nearest(x16, -)=', x16-x16L

write(1,*)
ulp16 =ulp_16(x16)
write(1,*) 'ulp(x) =', ulp16

if(x16 +ulp16 ==x16) then
write(1,*) 'x +ulp(x) ==x !'
else
write(1,'(a, 1pE43.33,a)') 'x +ulp(x) =', x16+ulp16, ' > x'
endif


if(x16+ulp16/2 ==x16) then
write(1,*) 'x +ulp(x)/2 ==x'
else
write(1,'(a, 1pE43.33, a)') 'x +ulp(x)/2 =', x16+ulp16/2., ' > x !'
! ulp/4
    if(x16+ulp16/4 ==x16) then
    write(1,*) 'x+ulp(x)/4 ==x'
    else
    write(1,*) 'x +ulp(x)/4 =', x16+ulp16/4., ' > x !!'
    endif
endif


write(1,*)
write(1,*) 'spacing(x16)=', spacing(x16)

! New "x":
900 continue
print '(/,a,\)', 'Another x (y/char): '
read (*,'(a1)') w1
if(w1 =='y' .or. w1 =='Y') then
write(1,*) '---'
goto 200
else
endif

1000 continue
write(1, '(a)') '--------------------------------------------'
write(1, '(/,a)') 'EPSILON:'
write(1, '(/,a)') 'Single:'
EPS =epsilon(1.0)
write(1,*) 'EPS:', EPS
u =1.+epsilon(1.0)

if(u >1.) then
write(1,*) '1.0 +EPS > 1.0'
else
write(1,*) '1.0 +EPS = 1.0'
endif
write(1,*)

!macheps =2**(-24)
write(1,*) 'EPS/2:', epsilon(1.0)/2.
u2 =1. +epsilon(1.0)/2.
u_ =1. -epsilon(1.0)/2.

if(u2 >1.) then
write(1,*) '1.0 +EPS/2 > 1.0'
else
write(1,*) '1.0 +EPS/2 = 1.0'
endif

if(u_ <1.) then
write(1,*) '1.0 -EPS/2 = ', 1.0 -EPS/2, ' < 1.0'
else
write(1,*) '1.0 -EPS/2 = 1.0'
endif

write(1, '(/,a)') 'Double:'
eps8 =epsilon(1.0_8)
write(1,*) 'eps8:', eps8

uu =1._8 +epsilon(1.0_8)

if( uu >1._8) then
write(1,*) '1.0_8 +eps8 > 1.0_8'
else
write(1,*) '1.0_8 +eps8 = 1.0_8'
endif
write(1,*)

!macheps =2**(-53)
write(1,*) 'eps8/2:', epsilon(1.0_8)/2.
u82 =1.0_8 +epsilon(1.0_8)/2._8
u8_ =1.0_8 -epsilon(1.0_8)/2._8

if(u82 >1.) then
write(1,*) '1.0_8 +EPS/2 > 1.0_8'
else
write(1,*) '1.0_8 +EPS/2 = 1.0_8'
endif

if(u8_ <1._8) then
write(1,*) '1.0_8 -eps8/2 =', 1._8 -eps8/2, ' < 1.0_8'
else
write(1,*) '1.0_8 -eps8/2 = 1.0_8'
endif


write(1, '(/,a)') 'Quadruple:'
eps16 =epsilon(1.0_16)
write(1,*) 'eps16:', eps16

u16 =1._16 +epsilon(1.0_16)
if( u16 >1._16) then
write(1,*) '1.0_16 +eps16 > 1.0_16'
else
write(1,*) '1.0_16 +eps16 = 1.0_16'
endif
write(1,*)

!macheps =2**(-113)
write(1,*) 'eps16/2:', epsilon(1.0_16)/2.
u162 =1.0_16 +epsilon(1.0_16)/2._16
u16_ =1.0_16 -epsilon(1.0_16)/2._16

if(u162 >1.) then
write(1,*) '1.0_16 +eps16/2 > 1.0_16'
else
write(1,*) '1.0_16 +eps16/2 = 1.0_16'
endif

if(u16_ <1._16) then
write(1,*) '1.0_16 -eps16/2 =', 1._16 -eps16/2, ' < 1.0_16'
else
write(1,*) '1.0_16 -eps16/2 = 1.0_16'
endif


write(1,*)
write(1,*) '* END *'

write(*, '(/,a)') ' Fisier de rezultate: ', trim(fname)


END
