akaGM
Platinum Member | Редактировать | Профиль | Сообщение | Цитировать | Сообщить модератору
call ftest
call dtest
call qtest
end
subroutine ftest
real alog2, fstore
real x
real*16 fppow2
character*36 note
write (6,'(''single precision:'')')
x = 1.0e+00
10 if (fstore(1.0e+00 + x/2.0e+00) .ne. 1.0e+00) then
x = x / 2.0e+00
go to 10
end if
write (6,'(8x,''sizeof(REAL) = '',i2)') 4
write (6,'(8x,''FLT_MANT_DIG = '',i2)')
x (1 + nint(-alog2(x)))
if (x .eq. real(fppow2(-23))) then
note = '[IEEE 754 32-bit macheps]'
else
note = '[not IEEE 754 conformant]'
end if
write (6,'(8x,''machine epsilon = '', 1pe14.5e4, 1x, a)')
x x, note
x = 1.0e+00
20 if (fstore(x/2.0e+00) .ne. 0.0e+00) then
x = x / 2.0e+00
go to 20
end if
if (x .eq. real(fppow2(-149))) then
note = '[IEEE 754 smallest 32-bit subnormal]'
else if (x .eq. real(fppow2(-126))) then
note = '[IEEE 754 smallest 32-bit normal]'
else
note = '[not IEEE 754 conformant]'
end if
write (6,'(8x,''smallest positive number = '', 1pe14.5e4, 1x, a)')
x x, note
write (6,'()')
end
subroutine dtest
double precision dlog2, dstore
double precision x
real*16 fppow2
character*36 note
write (6,'(''double precision:'')')
x = 1.0d+00
10 if (dstore(1.0d+00 + x/2.0d+00) .ne. 1.0d+00) then
x = x / 2.0d+00
go to 10
end if
write (6,'(8x,''sizeof(DOUBLE PRECISION) = '',i2)') 8
write (6,'(8x,''DBL_MANT_DIG = '',i2)')
x (1 + nint(-dlog2(x)))
if (x .eq. dble(fppow2(-52))) then
note = '[IEEE 754 64-bit macheps]'
else
note = '[not IEEE 754 conformant]'
end if
write (6,'(8x,''machine epsilon = '', 1pe14.5e4, 1x, a)')
x x, note
x = 1.0d+00
20 if (dstore(x/2.0d+00) .ne. 0.0d+00) then
x = x / 2.0d+00
go to 20
end if
if (x .eq. dble(fppow2(-1074))) then
note = '[IEEE 754 smallest 64-bit subnormal]'
else if (x .eq. dble(fppow2(-1022))) then
note = '[IEEE 754 smallest 64-bit normal]'
else
note = '[not IEEE 754 conformant]'
end if
write (6,'(8x,''smallest positive number = '', 1pe14.5e4, 1x, a)')
x x, note
write (6,'()')
end
subroutine qtest
real*16 qlog2, qstore
real*16 x
real*16 fppow2
character*42 note
write (6,'(''quadruple precision:'')')
x = 1.0q+00
10 if (qstore(1.0q+00 + x/2.0q+00) .ne. 1.0q+00) then
x = x / 2.0q+00
go to 10
end if
write (6,'(8x,''sizeof(REAL*16) = '',i2)') 16
write (6,'(8x,''LDBL_MANT_DIG = '',i3)')
x (1 + nint(-qlog2(x)))
if (x .eq. fppow2(-52)) then
note = '[IEEE 754 64-bit macheps]'
else if (x .eq. fppow2(-63)) then
note = '[IEEE 754 80-bit macheps]'
else if (x .eq. fppow2(-112)) then
note = '[IEEE 754 128-bit macheps]'
else
note = '[not IEEE 754 conformant]'
end if
write (6,'(8x,''machine epsilon = '', 1pe14.5e4, 1x, a)')
x x, note
x = 1.0q+00
20 if (qstore(x/2.0q+00) .ne. 0.0q+00) then
x = x / 2.0q+00
go to 20
end if
if (x .eq. fppow2(-1074)) then
note = '[IEEE 754 smallest 64-bit subnormal]'
else if (x .eq. fppow2(-1022)) then
note = '[IEEE 754 smallest 64-bit normal]'
else if (x .eq. fppow2(-16446)) then
note = '[IEEE 754 smallest 80-bit subnormal]'
else if (x .eq. fppow2(-16382)) then
note = '[IEEE 754 smallest 80- and 128-bit normal]'
else if (x .eq. fppow2(-16494)) then
note = '[IEEE 754 smallest 128-bit subnormal]'
else if (x .eq. fppow2(-16382)) then
note = '[IEEE 754 smallest 128-bit normal]'
else
note = '[not IEEE 754 conformant]'
end if
write (6,'(8x,''smallest positive number = '', 1pe14.5e4, 1x, a)')
x x, note
write (6,'()')
end
real function fstore(x)
real x
fstore = x
end
double precision function dstore(x)
double precision x
dstore = x
end
real*16 function qstore(x)
real*16 x
qstore = x
end
real*16 function fppow2(n)
integer n
integer p
real*16 x
if (n .lt. 0) then
p = -n
x = 1.0q+00 / 2.0q+00
else
p = n
x = 2.0q+00
end if
fppow2 = 1.0q+00
10 if (p .gt. 0) then
p = p - 1
fppow2 = fppow2 * x
go to 10
end if
end
real function alog2(x)
real x
alog2 = alog(x)/alog(2.0e+00)
end
double precision function dlog2(x)
double precision x
dlog2 = dlog(x)/dlog(2.0d+00)
end
real*16 function qlog2(x)
real*16 x
qlog2 = qlog(x)/qlog(2.0q+00)
end
|
Всего записей: 24122 | Зарегистр. 06-12-2002 | Отправлено: 17:28 25-04-2008 | Исправлено: akaGM, 18:39 25-04-2008 |
|
