Fortran - Re: why system providing routies will provide different result in

Author

Re: why system providing routies will provide different result in

Rich Townsend

2005-09-12, 6:58 pm

Ian Chivers wrote:
> c/c++ built in maths functions will convert/promote any argument to double
> and
> return a double result.
>
> literal numeric constants (e.g. 2.0) in fortran are c++ float.
> in c++ they are double.
>
> trig function evaluation of sin/cos/tan of 90 degrees
> in fortran (of c++ float) returns a result with cos and tan in the wrong
> quadrant.

I'm not sure what you mean by 'the wrong quadrant'. The trig functions return a
number between -1 and 1 (SIN and COS) or -Inf and +Inf (TAN); there is no notion
of 'quadrant', since the return value is not an angular measure.

If you instead mean the inverse trig functions, then the Fortran standard is
quite explicit about which quandrant the returned angular measure falls within.
To get a 'correct' quadrant, you should be using the ATAN2() function (there is
no ACOS2 or ASIN2, since there is no sensible way to define two-argument
variants of ACOS and ASIN).

cheers,

Rich

Brooks Moses

2005-09-14, 7:00 pm

Ian Chivers wrote:
> Here is a test program
[...]
> r=pi*angles(i)/180
> print *,angles(i),' ',sin(r),' ',cos(r),' ',tan(r)
[...]
> 61 0.8746197 0.4848096 1.804048
> 89 0.9998477 1.7452383E-02 57.29004
> 90 1.000000 -4.3711388E-08 -2.2877332E+07
> 91 0.9998477 -1.7452471E-02 -57.28975
>
> note the sign at 90 degrees of the cos and tan.

Yes. What's wrong with it?

Clearly, r is an approximation to pi/4 in any limited-precision
calculation, and thus these values are values for some angle that's
nearly but not exactly 90 degrees. Thus, I can't see any reason to
claim that either sign would be wrong here, unless you have some
specific claim about what side of pi/4 r happens to lie on.

- Brooks

--
The "bmoses-nospam" address is valid; no unmunging needed.

Rich Townsend

2005-09-14, 7:00 pm

Ian Chivers wrote:
> Here is a test program
>
> program
> ---------
>
> program trigtest
> implicit none
> integer , dimension(15) :: angles= &
> (/- 1,0,1,29,30,31,44,45,46,59,60,61,89,90,9
1/)
> real :: pi
> real :: r
> integer :: i
> pi=4.0*atan(1.0)
> do i=1,15
> r=pi*angles(i)/180
> print *,angles(i),' ',sin(r),' ',cos(r),' ',tan(r)
> end do
> end program trigtest
>
> end program
> --------------
>
> Here is the intel output
>
> -1 -1.7452406E-02 0.9998477 -1.7455066E-02
> 0 0.0000000E+00 1.000000 0.0000000E+00
> 1 1.7452406E-02 0.9998477 1.7455066E-02
> 29 0.4848096 0.8746197 0.5543091
> 30 0.5000000 0.8660254 0.5773503
> 31 0.5150381 0.8571673 0.6008607
> 44 0.6946584 0.7193398 0.9656888
> 45 0.7071068 0.7071068 1.000000
> 46 0.7193398 0.6946583 1.035530
> 59 0.8571673 0.5150381 1.664279
> 60 0.8660254 0.5000000 1.732051
> 61 0.8746197 0.4848096 1.804048
> 89 0.9998477 1.7452383E-02 57.29004
> 90 1.000000 -4.3711388E-08 -2.2877332E+07
> 91 0.9998477 -1.7452471E-02 -57.28975
>
> note the sign at 90 degrees of the cos and tan.
>
> we've found we get similar results
> with all of the compilers we currently use.

I'm still not following your point. COS(90) should ideally return zero, in
reality it returns some number close to zero. This number could either be
positive or negative. Why in your eyes should the number be positive?

cheers,

Rich

glen herrmannsfeldt

2005-09-14, 7:00 pm

Ian Chivers wrote:

> my main occupation is teaching programming
> to people with little or no numeric background.

> the mathematical value of the cosine 90 degrees is 0 and
> the tangent of 90 degrees is infinity.

> to the complete beginner to programming the actual
> values returned don't make much sense.

> i don't claim that the values are wrong, but that they
> are not at all obvious to the beginner.

This, along with different ways of doing argument reduction has
been discussed here and in other newsgroups. (*)

I am not sure how often this problem comes up, but, for
comparison purposes only, PL/I and some other languages
supply sind(), cosd(), and tand() for trig functions with
arguments in degrees. That avoids the problems of using
an approximate pi.

Maybe, though, this is a good example to use in analyzing the
problems that can come from using floating point arithmetic.
It is simple enough to analyze even by beginners.

(*) Some of the discussion involves doing argument reduction
using a very high precision value for pi. There may be some
reason for doing that, but I am not yet convinced, and some
of it has to do with the meaning of floating point numbers.

In physics, number are just about always approximations.
If one writes 2.0 they mean a value somewhere between
1.95 and 2.05. It makes sense, then, not to make any
assumptions about bits past the least significant bit
in a floating point value. Argument reduction with high
precision pi assumes that those bits are zero, which might be
true, but it isn't so obvious to me. It also means that there
is no floating point value whose cosine is zero or whose
sine is one, thus more reason for the degree functions.

At least one IBM machine had a flag that would supply ones
for fraction bits during post normalization. The idea was
that one could run a program both ways and compare the results.
The main problem is that you can have accidental cancellation
of such bits if you subtract two values, though it seems better
than nothing.

-- glen

glen herrmannsfeldt

2005-09-14, 7:00 pm

andy2O wrote:
(snip)

> OK, even good trig functions can only return closest floating point
> approximations - and I'm all for teaching it to all who will listen...

> But on reading the sample code, my first reaction was that the much
> bigger and more straightforward error would be if someone assumed that
> when angles==90, that the (rational) floating point value of r could
> ever be identical to pi/2 (an irrational number)....

I thought once about a floating point format with a single bit that
meant "multiply this value by pi", so that it would allow various
multiples of pi. See what Mathematica does with functions that
should return pi, and so with corresponding functions with pi
arguments.

> Clearly r==pi/2 is impossible for floating point r, so a (hypothetical)
> perfect precision cos function *should never* return cos(r)==0 for a
> floating point value of r. Infact if it did return cos(r)==0, wouldn't
> that be certain evidence that the cos function in question was not of
> perfect precision...

Traditionally argument reduction was done by dividing by an approximate
pi/2 or pi/4. In that case, it is possible to have an exact zero.
If you believe that floating point values approximate real values then
it seems fine to me. It seems that many now use a much more accurate
pi, doing argument reduction assuming that the bits past the LSB are
zero. That is, that floating point is not an approximation to a real
number but a subset of exact values for real numbers.

Some libraries will refuse to calculate trig functions on floating
point values large enough that the LSB has a place value close to or
greater than one. This was explained to me early in my Fortran career
when I accidentally passed a large value to sin(). It made sense
to me, but it seems not to everyone.

-- glen

N. Shamsundar

2005-09-15, 7:57 am

Ian Chivers wrote:
> Here is a test program
>

><--- SNIP ---->

we've found we get similar results
> with all of the compilers we currently use.
>
> it is one of the examples that Jane Sleightholme
> and I have developed to help beginners
> understand some of the problems they will
> encounter with floating point maths.
>

Here is a simpler, easier to appreciate, pedagogical example, suitable
for use with calculators and computers:

(4/3-1)*3-1

will give a non-zero result on most calculators. In fact, this
expression can be used to probe for the "machine-Epsilon" of the device.

Matlab on a PC gives

» (4/3-1)*3-1
ans =
-2.2204e-016
»

A TI-65 calculator gives -1 X 10^{-12}.

N. Shamsundar
University of Houston

N. Shamsundar

2005-09-15, 7:57 am

John Harper wrote:
> In article <1126735520.733036.224700@g49g2000cwa.googlegroups.com>,
> andy2O <andy2O@hotmail.com> wrote:
>
>
>
>
> Good idea, but when I looked at it just now it mentioned c99 but its
> latest Fortran was f77.
>

The language alluded to in the report does not matter. The issues dealt
with are intrinsic to finite-precision real arithmetic. The statements
of Kahan apply to any procedural language that has real types.

> Being an applied mathematician, it often seems to me that computer
> science is everything you can do in a computer except numerical
> analysis and floating-point arithmetic. Maybe I should add Fortran
> to that list.

That's an improper generalization. Most commercial applications of
computers involve little or no floating point operations. Early
microprocessors had no floating point instructions; some needed "math
co-processors" to provide the missing instructions. Floating point
arithmetic is a special topic that is of little interest to *most*
computer scientists.

>
> John Harper, School of Mathematics, Statistics and Computer Science,
> Victoria University, PO Box 600, Wellington, New Zealand
> e-mail john.harper@vuw.ac.nz phone (+64)(4)463 5341 fax (+64)(4)463 5045

N. Shamsundar
University of Houston

Rich Townsend

2005-09-15, 6:59 pm

N. Shamsundar wrote:
> Ian Chivers wrote:
>
>
>
>
> we've found we get similar results
>
>
> Here is a simpler, easier to appreciate, pedagogical example, suitable
> for use with calculators and computers:
>
> (4/3-1)*3-1
>
> will give a non-zero result on most calculators. In fact, this
> expression can be used to probe for the "machine-Epsilon" of the device.
>
> Matlab on a PC gives
>
> » (4/3-1)*3-1
> ans =
> -2.2204e-016
> »
>
> A TI-65 calculator gives -1 X 10^{-12}.

All Fortran compilers should give -1. How about that!

cheers,

Rich

Rich Townsend

2005-09-15, 6:59 pm

John Harper

2005-09-15, 6:59 pm

In article <dgbnpk$3u3k$1@masala.cc.uh.edu>,
N. Shamsundar <shamsundar_at_@uh.edu> wrote:
>John Harper wrote:
>
>That's an improper generalization.

I did not offer my comment as a strict definition; I said "it often
seems". Maybe I should also have added ":-)". Goldberg's "What every
computer scientist should know about floating point arithmetic" makes
the same point where it says

"It is not uncommon for computer system designers to neglect the parts
of a system related to floating-point. This is probably due to the fact
that floating-point is given very little (if any) attention in the
computer science curriculum."

I wish to withdraw my accusation that the latest Fortran referred to
was f77: f95 is mentioned in one place in the supplement. It would be
good to see a further supplement mentioning f2003 IEEE754 support,
and if f2003 overcomes some or all of the report's objections to
languages current at the time the first supplement was written it
would also be good for that to be mentioned. I refer to

"Current programming languages make it difficult for a program to
specify the precision it expects."

"In a portable program that uses double precision as its nominal
working precision, there are five ways we might want to control the
use of a wider precision: ...
No current language supports all five of these options."

If f2003 does not overcome those objections, will f2007?

John Harper, School of Mathematics, Statistics and Computer Science,
Victoria University, PO Box 600, Wellington, New Zealand
e-mail john.harper@vuw.ac.nz phone (+64)(4)463 5341 fax (+64)(4)463 5045

glen herrmannsfeldt

2005-09-15, 6:59 pm

John Harper wrote:

(snip referencing Goldberg's "What every computer scientist
should know about floating point arithmetic", which says:)

> "Current programming languages make it difficult for a program to
> specify the precision it expects."

> "In a portable program that uses double precision as its nominal
> working precision, there are five ways we might want to control the
> use of a wider precision: ...
> No current language supports all five of these options."

Extra precision helps some algorithms, hurts others, and
confuses some people who don't expect it.

> If f2003 does not overcome those objections, will f2007?

Much of the existing hardware doesn't make it easy to
do it much different. The standard isn't likely to
require things that the hardware can't supply easily.

-- glen

N. Shamsundar

2005-09-15, 9:57 pm

John Harper wrote, quoting Goldberg:
> "Current programming languages make it difficult for a program to
> specify the precision it expects."

In the next paragraph, Goldberg himself says:

> Language standards are not entirely to blame for the vagaries of expression evaluation.

Indeed, while an expert user should be provided with means for
*controlling* precision, simply specifying desired precision won't do.
For example, in computing a derivative with a first order difference
approximation, we are only entitled to half the precision in the
approximate derivative compared to the precision in the function value.

Once we move beyond the four basic operations of finite-precision real
arithmetic (for example, if we s

to perform matrix decompositions, or
to simulate fluid flow) we are in the domain of numerical analysis; to
do a decent job of teaching this subject to a computer science student,
we'd need at least two semesters after teaching multivariable calculus.
Few CS programs have room for this.

As Glen Herrmannsfeldt observes in this thread, giving additional
control over floating point precision is not necessarily beneficial to
the average user, who may be unprepared to use these esoteric features.

The language designers have to deliver a compromise that, while easy to
learn subset for most users and most applications, has special features
that are there for numerical analysts to use when needed.

And, while we are discussing "what every ... should know about ...", let
us ask how many of our descendants will learn how to compute square
roots or trigonometric functions by hand, since they may need to be able
to program microprocessors that will not have a built-in SQRT or ATAN
instruction.

N. Shamsundar
University of Houston

Sponsored Links