jdr@sloth.mlb.semi.harris.com (Jim Ray) (04/09/91)
As the topic states, I am looking for some fortran program that one could use to compare the numeric (floating point) results with some "known" values? Specifically, I am looking for such a program that would exercise (extensively and repeditively ) most of the transcendental functions plus the normal run of the mill ones. This "program" would checking for roundoff and complience with the IEEE standard floating point specs. The reason I am looking for such a program, is the wide divergence of numerical results we have seen between several releases of fortran compilers on several workstations including SUN. Problem is that I can't figgure which one(s) are "correct". Any comments are welcome.... Thanks, -- Jim Ray Harris Semiconductor Internet: jdr@semi.harris.com PO Box 883 MS 62B-022 Phone: (407) 729-5059 Melbourne, FL 32901
alex@IASTATE.EDU (Alexander Roger K) (04/10/91)
In article <1991Apr9.023737.18140@mlb.semi.harris.com>, jdr@sloth.mlb.semi.harris.com (Jim Ray) writes: > I am looking for some fortran program that ... would exercise > transcendental functions ... checking for roundoff and complience with > the IEEE standard floating point specs. > > Jim Ray Harris Semiconductor > Internet: jdr@semi.harris.com PO Box 883 MS 62B-022 > Phone: (407) 729-5059 Melbourne, FL 32901 At least two such programs are available from netlib. To quote from the netlib index: ELEFUNT is a collection of transportable Fortran programs for testing the elementary function programs provided with Fortran compilers. The programs are described in detail in the book "Software Manual for the Elementary Functions" by W. J. Cody and W. Waite, Prentice Hall, 1980. PARANOIA is a rather large program, devised by Prof. Kahan of Berkeley, to explore the floating point system on your computer. "PARANOIA" comes in single and double precision versions for Fortran, and there is a BASIC ;-) version too. To see what is available send mail to netlib@ornl.gov -- send index from elefunt send index from paranoia -- *___________________________________________________________________________* | Roger Alexander, alex@iastate.edu | "...if nonreligion is a religion | | It's my *job* to express opinions |too, isn't *everybody* tax-exempt?"| | differing from those of my employer. | --Barbara Ehrenreich | *_______________________________________|___________________________________*
reilly@motcid.UUCP (Patrick L. Reilly) (04/13/91)
jdr@sloth.mlb.semi.harris.com (Jim Ray) writes:
:As the topic states, I am looking for some fortran program that one
:could use to compare the numeric (floating point) results with some
:"known" values? Specifically, I am looking for such a program that
:would exercise (extensively and repeditively ) most of the
:transcendental functions plus the normal run of the mill ones. This
:"program" would checking for roundoff and complience with the IEEE
:standard floating point specs.
You need to get a copy of "Testing Unconstrained Optimization Software"
by More, J. et. al, ACM Trans. on Mathematical Software, v. 7,n. 1,
March 1981, pages 17-41.
Although what you specifically seek is a bit different from
the referenced paper, you will find much insight and routines to
test with in the paper. Being up to date with Rosenbrock's,
Freudenstein and Roth's, Beale, Helical valley, etc. functions can
go a long way towards your problem.
As for roundoff, you might try the i1mach, r1mach, and d1mach
routines from the SLATEC library to get important parameters about your
target machine. Alternatively, give some of the code below a whirl.
Good luck,
========================== `` '' ===================================
Patrick Reilly, Ph.D. < @ @ > tel: 1+708+632-3109
Motorola CIG/Switch Dev. ( > ) fax: 1+708+632-2413
1501 W. Shure Dr. -BC569 \~/ UUNET: uunet!motcid!reilly
Arlington Hts., IL 60004
Any ill-formed premise may be
REilly's CAutioNary Truth : modeled and validated with
-- -- - - simulation.
====================================================================
Roundoff.f:
integer i
real h
double x
i = 0
x = 0.0
h = 0.1
10 continue
i = i + 1
x = x + h
print *, i, x
if(x.le.1.0)go to 10
end
Epsilon.f:
implicit double precision(a-h,o-z)
eps= 1.0
10 continue
eps= eps/2.0
t= 1.0 + eps
if(t.gt.1.0) go to 10
eps = 2.0*eps
print *,eps,'.le.machine epsilon.le.',2.0 * eps
end
****************************************************************
First results from the CrayX....
Output from epsilon.f:
2.5243548967072377773175314089E-29.le.machine
epsilon.le.5.0487097934144755546350628178E-29
(Epsilon.f output can be considered, roughly speaking, to be the fractional
accuracy to which floating point numbers are represented and
corresponds to a change of one in the least significant bit of the
mantissa. That is, machine epsilon is a positive floating-point
number for which fl(1+epsilon)>= 1, where fl(x) is the
floating-point representation of a real number x. This code yields this
value, epsilon, within a factor of two.)
Output from roundoff.f:
1, 1.0000000000000008881784197001E-1
2, 2.0000000000000017763568394002E-1
3, 3.0000000000000026645352591003E-1
4, 4.0000000000000035527136788005E-1
5, 5.0000000000000044408920985006E-1
6, 6.0000000000000053290705182007E-1
7, 7.0000000000000062172489379008E-1
8, 8.000000000000007105427357601E-1
9, 9.0000000000000079936057773011E-1
10, 1.0000000000000008881784197001E+0
Compare the above with the following after changing x = x + h to read
x = i * h in the roundoff.f program:
1, 1.0000000000000008881784197001E-1
2, 2.0000000000000017763568394002E-1
3, 3.000000000000007105427357601E-1
4, 4.0000000000000035527136788005E-1
5, 5.E-1
6, 6.000000000000014210854715202E-1
7, 6.9999999999999928945726423989E-1
8, 8.000000000000007105427357601E-1
9, 9.000000000000021316282072803E-1
10, 1.E+0
11, 1.1000000000000014210854715202E+0
If exact math were possible on the machine, both outputs should be 11 lines!
The problem with the 1st run is that the fraction 1/10 cannot be exactly
represented by the bits used in the machine word length. Note the direction
of the rounding errors...when i=10, x is greater than 1 and the loop stops.
Now for the Cray2:
Output from epsilon.f:
2.5243548967072377773175314D-29.le.machine
epsilon.le.5.0487097934144755546350628D-29
Output from roundoff.f:
1, 1.0000000000000008881784197D-1
2, 2.0000000000000017763568394D-1
3, 3.0000000000000026645352591D-1
4, 4.0000000000000035527136788D-1
5, 5.0000000000000044408920985D-1
6, 6.0000000000000053290705182D-1
7, 7.0000000000000062172489379D-1
8, 8.0000000000000071054273576D-1
9, 9.0000000000000079936057773D-1
10, 1.0000000000000008881784197D+0
Output from roundoff.f with same modification as above:
1, 1.0000000000000008881784197D-1
2, 2.0000000000000017763568394D-1
3, 3.0000000000000071054273576D-1
4, 4.0000000000000035527136788D-1
5, 5.0D-1
6, 6.0000000000000142108547152D-1
7, 6.9999999999999928945726423D-1
8, 8.0000000000000071054273576D-1
9, 9.0000000000000213162820728D-1
10, 1.0D+0
11, 1.1000000000000014210854715D+0
On a Sun 3/80, epsilon.f produces:
2.2204460492503D-16.le.machine
epsilon.le. 4.4408920985006D-16
roundoff.f output:
1 1.0000000149012D-01
2 0.20000000298023
3 0.30000000447035
4 0.40000000596046
5 0.50000000745058
6 0.60000000894070
7 0.70000001043081
8 0.80000001192093
9 0.90000001341105
10 1.0000000149012
Now, making the same modifications to the code as above:
1 1.0000000149012D-01
2 0.20000000298023
3 0.30000001192093
4 0.40000000596046
5 0.50000000000000
6 0.60000002384186
7 0.69999998807907
8 0.80000001192093
9 0.90000003576279
10 1.0000000000000
11 1.1000000238419
The epsilon value is useful for terminating iterations instead of waiting for
an output that will never come while the computer churns along seeking an exact
answer. Epsilon could be the stopping criteria.
The programs used were taken from Kalbasi, "Can You Trust Your Computer?",
IEEE Potentials, April 1990, pp15-18.