In article <4f9bfc52$0$1385$4fafb...@reader1.news.tin.it>,
Giorgio Pastore <past...@units.it> wrote:
> In any case, while many dark corners of floating point arithmetic are > discussed in many books, I do not know a single book explicitly > mentioning the possibility that ((1.0,0.0)-(0.1,0.0))/(10.0,0.0) and > (1.0-0.1)/10.0 (or similar expressions) could be different.
Then they should, because most experienced fortran programmers know this already. Even the REAL expression that you give above can result in different values depending on, for example, optimization levels. Or in the context of a more complicated program, the compiler itself might decide in one case the keep the intermediate result in an extended precision register in one case, but flush to memory in order to use that register for another value in another case.
The complex expression is simply one example of this, the same real expression embedded within a different context in the program. The compiler might well want to use that register for another part of that complex expression evaluation, flushing the first result to memory (and truncating its mantissa in the process).
I don't think anyone mentioned this previously, but your original example:
demonstrates another kind of programming error. This is exactly the WRONG way to compute that root with b1>0 because of the large roundoff error. There are several expressions that can be used to compute that root without that roundoff error. I assume that was the point of the original exercise, but as I say, I don't think anyone mentioned this previously.
Giorgio Pastore <past...@units.it> wrote:
> On 4/28/12 1:49 PM, robin.vow...@gmail.com wrote:
> ...
>> Before making a statement like that, you ought to correct
>> the program.
> Robin, if you followed the thread, you realize that the mistaken > expression for kind rk2 does not affect the main issue I was > discussing.
You have to get used to Robin. That always happens.
>.... I assume that was
> the point of the original exercise, but as I say, I don't think
> anyone mentioned this previously.
Exactly. I started my post writing explicitly that the problem came from a program written to illustrate the dangers of subtractive cancellations for didactic purposes.
The whole discussion has been very interesting and instructive for me. The reasons for the difference were quite clear to me, although I did not think to use the -ffloat-store option. Still, it has been the first time, after years of numerical programming, that I realized that there is nothing in the standard requiring consistency between real and complex with zero imaginary part. I guess that, stated this way, most of my colleagues would be surprised as well. The fact that *after* thinking a moment the behavior is perfectly understandable, does not make it trivial at all, in my opinion and my point is that it should deserve more emphasize than presently. Then, if the emphasis should be put in books, in the examples of the standard or in compiler doc, it is probably matter of tastes.
> I used to show students the perils of misusing floating-point arithmetic by > getting them to calculate 1 - x + x**2/2! - x**3/3! + ... for x = 10.0 or > some such value. (It was many years ago and I don't now remember that > detail, nor what language and machine we then used.) Of course they got a > wide variety of answers, none near exp(-x),
For double precision, it is not too bad.
For single precision, the old trick of explicit summation over
alternating terms gives an improvement, but still a wrong answer
by a factor of five:
program main
implicit none
integer :: i
integer, parameter :: rp=selected_real_kind(6)
real(kind=rp) :: x,s,xp, fak, newterm
x = 7._rp
s = 0._rp
fak = 1._rp
xp = 1._rp
do i=0,100,2
fak = fak*(i+1)
newterm = xp * (i+1-x)/fak
if (abs(newterm) < abs(epsilon(s)*s)) exit
s = newterm + s
fak = fak*(i+2)
xp = xp * x**2
end do
print *,s
print *,exp(-x)
end program main
My five cents worth (about 4 UK pennies in 1945):
My first introduction to commercial use of Fortan was seeing runs of Shell's
UK acounting system being run on our IBM service bureau office 1401 computer
in 1960. (It had a Fortran compiler and ran in 16k. I was in the banking and
stockmarket sales area).Years later I joined Shell and once had a 3-year
stint as an auditor before actually writing accounting systems.
The point I'm getting to is that in all companies I worked in over 50 years,
the accounting (especially my own later comercial systems) was usually
written in Fortran (unless PL/1). The maths was always devoid of floating
point operations, being done entirely in integer arithmetic (based on
pennies in the UK and cents in USA, Caribbean and Central America, and, I
guess, every where else with decimal monetary systems. (I have a life-time
set of discovered-fraud stories about currency exchange accounting).
In survey software you HAVE to use integers to get the row, column and grand
totals to correlate.
The fun starts when percentage-based tables, so loved by Ad. sales
directors, also had to add up exactly. (There's a running-error propagation
and absorbance algorithm to handle that).
On Monday, 26 March 2012 12:41:52 UTC+11, robert....@oracle.com wrote:
> I recently purchased a copy of the book "Methods of Numerical
> Integration" by Ralston and Rabinowitz. The edition I purchased is
> the Dover Publications reprint dated 2007. It is a reprint of the
> second edition, and it is copyright 1984.
> The code examples are written in FORTRAN. The first chapter includes
> the following explanation for why FORTRAN is used.
> "ALGOL, or more accurately ALGOL 60, has virtually disappeared as
> a language for the dissemination of algorithms. Of the viable
> descendants of ALGOL, ALGOL 68 has not achieved wide acceptance in the
> computing community while PASCAL does not have the algorithmic power
> needed for mathematical software. Consequently, in spite of its many
> shortcomings, the venerable programming language FORTRAN has become
> the vehicle for almost all current mathematical software."
> Another statement I find interesting appears in Section 4.2.1
> "Alleviation of Roundoff Error."
> "Finally, it is our opinion that certain high-level languages
> such as APL offer the user seductively easy programming possibilities,
> and by concealing many numerical processes lead inevitably away from
> hi-fi numerical practice."
On Sunday, 29 April 2012 02:03:36 UTC+10, glen herrmannsfeldt wrote:
> Giorgio Pastore > wrote:
> > On 4/28/12 1:49 PM, robin.vow...@gmail.com wrote:
> > ...
> >> Before making a statement like that, you ought to correct
> >> the program.
> > Robin, if you followed the thread, you realize that the mistaken > > expression for kind rk2 does not affect the main issue I was > > discussing.
> You have to get used to Robin. That always happens.
Snide remarks are unbecoming.
Especially when your response is to a remark that is incorrect.
On Sunday, 29 April 2012 00:14:56 UTC+10, Giorgio Pastore wrote:
> On 4/28/12 1:49 PM, robin.vow...@gmail.com wrote:
> ...
> > Before making a statement like that, you ought to correct
> > the program.
> Robin, if you followed the thread, you realize that the mistaken > expression for kind rk2 does not affect the main issue I was discussing.
In you post, you were concerned with differences between
the result of a Real computation, and an identical computation
using Complex arithmetic.
You even suggested that there was a fault in library routine(s).
>However, with my big surprise, gfortran 4.4,4.5, 4.6 and 4.7 on ubuntu
>boxes give different results if one works with real or complex
>variables. In particular the real results is definitely more accurate
>than the complex one.
>I suspect something went wrong with the last gfortran distributions or
>the related libraries.
Since one of your computations - which contained an error - produced an
expected result, you concluded incorrectly that the program worked.
> On 4/28/12 2:48 AM, Steven G. Kargl wrote:
> ...
> > You failed to identify the hardware you used, but I suspect it is
> > i686 compatible.
> You are right. It is i686 compatible.
> > If it is, then the answer is (1) get better
> > hardware or (2) use -ffloat-store if you don't want intermediate
> > results stored in FPU 80-bit registers. On my i686-*-freebsd system,
> Now I can get a consistent result like you.
> > 1.0100000 is the real part of ( 1.0100000 , 0.0000000 )
> > 2000.0100 is the real part of ( 2000.0100 , 0.0000000 )
> > real math result -4.0000000
> > complex math result ( -4.0000000 , 0.0000000 )
> Ok, thank you and other people who helped me to understand the origin of
> such a behavior and the consistency with the standard.
> Still, I remain a little puzzled by the fact that a straightforward
> translation of the program in C99 (using float and complex C variables)
> and its compilation with gcc (same version as gfortran) provides a
> consistent
> 1.010000 is the real part of 1.010000 0.000000
> 2000.010010 is the real part of 2000.010010 0.000000
> real -4.040000
> complex -4.040000 0.000000
> *without* any special compilation option.
> Moreover, I think that the Fortran Standard should at least give a
> warning about this possibility of inconsistent results when using simple
> arithmetic operations between real or zero-imaginary part complex numbers.
The Fortran standard includes the mathematical equivalence rule, which
allows
any expression to be replaced with a mathematically equivalent
expression.
That rule effectively means you cannot trust the result of evaluating
an
expression. In the past, the power of optimizers was very limited,
which
protected users more than they knew. Over time, optimizers have
become
far more aggressive, and programs that used to work without problems
now
produce compromised results.
> On 4/28/12 2:48 AM, Steven G. Kargl wrote:
> ...
>> You failed to identify the hardware you used, but I suspect it is
>> i686 compatible.
> You are right. It is i686 compatible.
>> If it is, then the answer is (1) get better
>> hardware or (2) use -ffloat-store if you don't want intermediate
>> results stored in FPU 80-bit registers. On my i686-*-freebsd system,
One might also consider adjusting the default FPU precision or force the
compiler to use SSE operations (which will typically match the desired
precision). For gfortran, the -mpc64 option will usually result in the desired
precision without costly forcing register contents to memory.
If computations in multiple precisions are desired in a single program instead,
one needs to emit corresponding control register changes before issuing the
operation.
x86_64 aka amd64 does not exhibit this particular problem as severely, because
FP operations use the SSE registers by default.
> On Sunday, 29 April 2012 02:03:36 UTC+10, glen herrmannsfeldt wrote:
> > Giorgio Pastore
> > wrote:
> > > On 4/28/12 1:49 PM, robin.vow...@gmail.com wrote:
> > > ...
> > >> Before making a statement like that, you ought to correct
> > >> the program.
> > > Robin, if you followed the thread, you realize that the mistaken
> > > expression for kind rk2 does not affect the main issue I was
> > > discussing.
> > You have to get used to Robin. That always happens.
> Snide remarks are unbecoming.
> Especially when your response is to a remark that is incorrect.
