Why FORTRAN

[robert....@oracle.com]

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."

Robert Corbett

[Louis Krupp]

On Sun, 25 Mar 2012 18:41:52 -0700 (PDT), 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.

<snip>

> "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."

Do they offer a definition of "hi-fi numerical practice" anywhere?

Louis

[robert....@oracle.com]

On Mar 25, 6:41 pm, robert.corb...@oracle.com wrote:
> I recently purchased a copy of the book "Methods of Numerical
> Integration" by Ralston and Rabinowitz.

^^^^^^^
Davis

Oops. The book is by Davis and Rabinowitz. I purchased a book
by Ralston and Rabinowitz at the same time.

Robert Corbett

[robert....@oracle.com]

On Mar 25, 10:47 pm, Louis Krupp <lkr...@nospam.pssw.com.invalid>
wrote:
> On Sun, 25 Mar 2012 18:41:52 -0700 (PDT), robert.corb...@oracle.com

No, they do not. The phrase obviously is short for "high-fidelity,"
but the expanded phrase still requires interpretation.

I have seen enough use of array expressions in Fortran to think I know
what they mean. For example, I saw a bug report from someone who
replaced a DO-loop for doing a summation with a call of the
intrinsic function SUM, and then wondered why the results he got were
different. The call of SUM was software pipelined, which changed the
associativity of the summation.

Robert Corbett

[John Harper]

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), depending on exactly how they
had coded it. I hoped they would learn that conevrgence, and usefulness in
a computer, are quite different properties of series, and one can have the
first without the second. (When they did asymptotic series later they got
the second without the first.)

--
John Harper

[glen herrmannsfeldt]

John Harper <john....@vuw.ac.nz> wrote:

(snip)

> 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.)

Another one I have seen used is the quadratic formula.
For some specific values of A, B, and C, the results are
very different from the actual roots, though they don't look
at all unusual.

> Of course they got a wide variety of answers, none near
> exp(-x), depending on exactly how they had coded it.

More usual in actual problem solving are recursion relations
that are unstable, especially for functions defined by differential
equations.

-- glen

[Giorgio Pastore]

On 4/27/12 2:36 AM, glen herrmannsfeldt wrote:
....

> Another one I have seen used is the quadratic formula.
> For some specific values of A, B, and C, the results are
> very different from the actual roots, though they don't look
> at all unusual.

Your observation reminded me that I have a related question I planned to
ask here.

Exactly for the purpose of illustrating the dangers of subtractive
cancellations, I recently asked my students (introductory computational
physics course), to solve the quadratic equation x^2+8000x+1=0 by using
the usual formula.

In the past, I checked on different machines that default 32-bit
precision and double 64-bit precision nicely demonstrated a huge
relative error of default precision real arithmetic.

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.

Before considering such behavior as a bug, I am wondering if the
standard explicitely requires consistent results from real arithmetics
and complex arithmetics (with zero imaginary part) or if such
adiscrepancy is allowed by the standard.

Giorgio Pastore

PS the code I have used to check the behavior is the following:

program quadratic
implicit none
integer,parameter :: rk1=selected_real_kind(6)
integer,parameter :: rk2=selected_real_kind(15)
real(kind=rk1) :: a1,b1,c1
real(kind=rk2) :: a2,b2,c2

real (kind=rk1) :: discr1, xpr1
complex(kind=rk1) :: xpc1
real (kind=rk2) :: discr2, xpr2
complex(kind=rk2) :: xpc2

a1=1.0
b1=8000.0
c1=1.0

a2=a1
b2=b1
c2=c1


discr1 = b1**2 - 4*a1*c1
xpr1 = ( -b1 + sqrt(discr1) )/(2*a1)
xpc1 = ( -b1 + sqrt(cmplx(discr1,0.0_rk1) ) )/(2*a1)

print*,xpc1
print*,xpr1

discr2 = b2**2 - 4*a2*c2
xpr2 = ( -b2 + sqrt(discr2) )/(2*a2)
xpc2 = ( -b2 + sqrt(cmplx(discr2,0.0_rk2) ) )/(2*a2)

print*,xpc2
print*,xpr2
end program quadratic

On linux/ubuntu 10.04.4 and gfortran 4.4 I get
(-2.44140625E-04, 0.0000000 )
-1.25000006E-04
(-2.44140625000000000E-004, 0.0000000000000000 )
-1.25000001953035067E-004


while the same program compiled with gfortran 4.3.1 on my very old
powerPC G4 gives:

(-2.44140625E-04, 0.0000000 )
-2.44140625E-04
(-1.25000001844455255E-004, 0.0000000000000000 )
-1.25000001844455255E-004

with full equality between real and complex results.

[Richard Maine]

Giorgio Pastore <pas...@units.it> wrote:

> Before considering such behavior as a bug, I am wondering if the
> standard explicitely requires consistent results from real arithmetics
> and complex arithmetics (with zero imaginary part) or if such
> adiscrepancy is allowed by the standard.

The standard doesn't require much of *ANYTHING* about the accuracy of
floating-point arithmetic. Pretty much everything floating point is
defined just to be "a processor-dependent approximation." It has been
noted that a standard-conforming processor could return 42.0 as the
result of every floating point operation; that might not be a very good
approximation, but the standard doesn't require it to be very good.

It might be a slight bit of hyperbole in that I would not swear that one
couldn't argue 42.0 to be an invalid result in some cases, but it still
illustrates the general principle.

The IEEE stuff is probably one major exception, but that's an f2003
feature.

--
Richard Maine
email: last name at domain . net
domain: summer-triangle

[Steven G. Kargl]

On Sat, 28 Apr 2012 00:01:42 +0200, Giorgio Pastore wrote:

>
> I suspect something went wrong with the last gfortran distributions or
> the related libraries.

I suspect that there is a bug in your program.

--
steve

[Giorgio Pastore]

I do not see any bug in my program. If you find one, I would appreciate
very much to know where.

Giorgio Pastore

[Steven G. Kargl]

What is the kind type parameter for a value returned by CMPLX()?

xpc2 = ( -b2 + sqrt(cmplx(discr2,0.0_rk2) ) )/(2*a2)

--
steve

[glen herrmannsfeldt]

Richard Maine <nos...@see.signature> wrote:
> Giorgio Pastore <pas...@units.it> wrote:

>> Before considering such behavior as a bug, I am wondering if the
>> standard explicitely requires consistent results from real arithmetics
>> and complex arithmetics (with zero imaginary part) or if such
>> adiscrepancy is allowed by the standard.

> The standard doesn't require much of *ANYTHING* about the accuracy of
> floating-point arithmetic. Pretty much everything floating point is
> defined just to be "a processor-dependent approximation."

I was about to reply to the previous post, but then decided to
see if you had already replied.

> It has been
> noted that a standard-conforming processor could return 42.0 as the
> result of every floating point operation; that might not be a very good
> approximation, but the standard doesn't require it to be very good.

> It might be a slight bit of hyperbole in that I would not swear that one
> couldn't argue 42.0 to be an invalid result in some cases, but it still
> illustrates the general principle.

It seems to me that might not be allowed for sin() and cos() for real
arguments. It could be for complex sin() and cos(), though.

> The IEEE stuff is probably one major exception, but that's an f2003
> feature.

But isn't that only if the compiler claims to have IEEE support?

-- glen

[Giorgio Pastore]

On 4/28/12 1:35 AM, Steven G. Kargl wrote:
...

> What is the kind type parameter for a value returned by CMPLX()?
>
> xpc2 = ( -b2 + sqrt(cmplx(discr2,0.0_rk2) ) )/(2*a2)
>

You're right. As usual, afer a long set of checks, the final version of
the program is not the right one :-(.

Of course I did the chech with cmplx(discr2,0.0_rk2,rk2). But the result
I was carefully monitoring was the result with kind=rk1.

In any case, the following much simplified test should make clearer the
point.

program quadratic
implicit none
real :: a1,b1,x1
complex :: a2,b2,x2

a1=1.01
b1=2000.01

a2=cmplx(a1,0.0)
b2=cmplx(b1,0.0)

print*,a1,' is the real part of ',a2
print*,b1,' is the real part of ',b2
x1 = ( -b1**2 + (b1**2-4*a1) )
x2 = ( -b2**2 + (b2**2-4*a2) )
print*,' real math result ',x1
print*,' complex math result ',x2
end program quadratic

gfortran 4.4 result:

1.010000 is the real part of ( 1.010000 , 0.000000 )
2000.010 is the real part of ( 2000.010 , 0.000000 )
real math result -4.040000
complex math result ( -4.000000 , 0.000000 )


Can be considered an acceptable result ?


Giorgio Pastore

[glen herrmannsfeldt]

Steven G. Kargl <s...@removetroutmask.apl.washington.edu> wrote:

(snip, someone wrote)

>> I do not see any bug in my program. If you find one,
>> I would appreciate very much to know where.

> What is the kind type parameter for a value returned by CMPLX()?

> xpc2 = ( -b2 + sqrt(cmplx(discr2,0.0_rk2) ) )/(2*a2)

That is true, but it isn't the problem. Note that the single and
double complex give the same result, but significanlty different
from the single and double real result.

As to the OP, note that COMPLEX functions can easily give different
results from the corresponding REAL functions.

For one if, CABS(Z) = SQRT(REAL(Z)**2+AIMAG(Z)**2)

then it can easily overflow in cases where ABS(REAL(Z)) wouldn't.

(As a side note, I will complain about the use of REAL and CMPLX
functions for two different uses: One for converting to/from complex
values (the Fortran 77 sense), and for converting to the appropriate
Fortran data type and kind.)

According to the library with source nearby:

* COMPLEX SQUARE ROOT FUNCTION (SHORT)
* 1. THE PRINCIPLE BRANCH OF THE SQUARE ROOT IS TAKEN,
* I.E., THE REAL PART OF THE ANSWER IS POSITIVE.
* 2. WRITE SQRT(X+IY) = U+IV, WHERE U IS REAL, AND V IS
* IMAGINARY. IF X=Y=0, U=V=0.
* 3. IF X IS NON-NEGATIVE, U = SQRT((/X/+/X+IY/)/2) AND
* V = Y/(2*U).
* 4. IF X IS NEGATIVE, U = Y/(2*V) AND
* V = SIGN(Y)*SQRT((/X/+/X+IY/)/2).

Where the /X/ is the absolute value of X.

It then computes, conceptually, U=SQRT((X+SQRT(X**2+Y**2))/2).

Now, the actual program where I got that does a somewhat more
detailed calculation to avoid some intermediate underflow and
overflow, but you can see that, in the case of cancellation,
that the results could easily be different.

The actual calculation done in this case is:

A=MAX(X,Y), B=MIN(X,Y)
(If B is zero, or enough less than A, skip much of the computation.)

Otherwise, computes D=B/A, then F=SQRT(0.25+D*D/4) and /X+IY/=2*A*F.
Then more complication to avoid rounding problems computing
(/X/)/2+A*F, (/X/)+2*A*F, or (/X/)/4+A*F/2.

Note that the error source in computing CSQRT() are similar
to those in the quadratic formula, in the cases that can result
in cancellation from subtraction of nearly equal values.

In addition to all above, on systems with x87 floating point,
some intermediates might be computed in 80 bit registers,
while others are stored as 64 bits. That can easily result in
differences for the same expression in the same program.

-- glen

[Steven G. Kargl]

You failed to identify the hardware you used, but I suspect 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,

laptop:kargl[224] ./z
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.0400000
complex math result ( -4.0391626 , 0.0000000 )
laptop:kargl[225] gfc4x -o z -ffloat-store d.f90
laptop:kargl[226] ./z
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 )

--
steve

[Tim Prince]

Did you perhaps specify x87 extended precision (the default for -m32
mode), getting in effect at least double precision evaluation of real
expressions? complex frequently used less extended precision than
real. Options such as -fcx-limited-range could make a difference as
well. It's hard to sympathize with a blanket statement about gfortran
which may be valid only for a specific choice of options.
The standards never forbade use of extra precision, even though it might
produce inconsistencies.

[Richard Maine]

glen herrmannsfeldt <g...@ugcs.caltech.edu> wrote:

> Richard Maine <nos...@see.signature> wrote:

> > The standard doesn't require much of *ANYTHING* about the accuracy of
> > floating-point arithmetic. Pretty much everything floating point is
> > defined just to be "a processor-dependent approximation."

...

> > It has been
> > noted that a standard-conforming processor could return 42.0 as the
> > result of every floating point operation; that might not be a very good
> > approximation, but the standard doesn't require it to be very good.
>
> > It might be a slight bit of hyperbole in that I would not swear that one
> > couldn't argue 42.0 to be an invalid result in some cases, but it still
> > illustrates the general principle.
>
> It seems to me that might not be allowed for sin() and cos() for real
> arguments. It could be for complex sin() and cos(), though.

Those are exactly some of the examples that occured to me. I didn't
research enough to be confident of the answer either way, but those are
among the ones I'd look at if I actually cared enough to go to the work.
I could imagine people caring about it though. Not the silly 42.0 bit,
but whether sin() or cos() might be allowed to return something more
like 1.00000000001 (or therebouts), which might actually happen with
some algorithms and might upset some codes. Probably not with "decent"
algorithms, but then that's circular because that's probably one of the
criterion for one being "decent."

> > The IEEE stuff is probably one major exception, but that's an f2003
> > feature.
>
> But isn't that only if the compiler claims to have IEEE support?

Yes.

--
Richard Maine | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle | -- Mark Twain

[Giorgio Pastore]

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.

Giorgio

[Giorgio Pastore]

Moreover, I would be curious to know what is the default behavior of
other fortran compilers on i686 hardware.

Giorgio

[glen herrmannsfeldt]

Giorgio Pastore <pas...@units.it> wrote:

(snip)

> 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.

As previously noted, the inconsistent results are more general
than that, but that is just the case that you tested.

The standard doesn't require consistent results in any expression.

The x87 case has an interesting history. Part of the design
for the original 8087 included a virtual stack. When the eight
register stack on the processor overflowed, an interrupt routine
would copy some registers to memory. When it underflowed, they would
be copied back. However, no-one tried to write the interrupt
routine until after the chip was fabricated, and it turned out
not to be possible. So programs have to work within the eight
register stack.

The result, in combination with the way compilers use it,
is that some intermediate results have a 64 bit significand,
and others 53 bits. Without optmization, it would be usual for
the results of any expression to be stored in the given variable.

But with optimizations like common subexpression elimination,
it might keep an intermediate value at temporary real (the intel
name) precision. In that case, one can't be sure even by
storing into a variable.

And finally, as I previously noted, the computation of complex
expressions is often more, umm, complex, than real expressions.
Intermediate computations, such as squaring, can lose precision
in some cases when cancellation occurs.

-- glen

[Richard Maine]

Giorgio Pastore <pas...@units.it> wrote:

> 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.

Why in the world would the standard call out that particular case? Darn
near *everything* about floating-point numbers can give inconsistent
results per the standard. Some things actuyally do give surprising
results in actual implementations. The standard allows far worse
surprises than you'll actually see.

Note also that the standard is generally not a tutorial document. People
who need that level of warning would probably be better directed to read
other materials such as some of the books on the subject.

[robin....@gmail.com]

Before making a statement like that, you ought to correct
the program.

Before making a statement like that, you ought to correct
the program.

In the calculation of xpc2 in

xpc2 = ( -b2 + sqrt(cmplx(discr2,0.0_rk2) ) )/(2*a2)

CMPLX returns a single precision result, not double precision.
To obtain a double-precision result, you must specify a kind for CMPLX.

[robin....@gmail.com]

On Saturday, 28 April 2012 21:49:59 UTC+10, robin....@gmail.com wrote:

> Before making a statement like that, you ought to correct
> the program.
>
> In the calculation of xpc2 in
>
> xpc2 = ( -b2 + sqrt(cmplx(discr2,0.0_rk2) ) )/(2*a2)
>
> CMPLX returns a single precision result, not double precision.
> To obtain a double-precision result, you must specify a kind for CMPLX.

Results from program as posted:

(-2.441406E-04,0.0000)
-1.250000E-04
(-2.441406250000E-04,0.0000000000)
-1.250000019530E-04

Results from program after correction of CMPLX:

(-2.441406E-04,0.0000)
-1.250000E-04
(-1.250000018445E-04,0.0000000000)
-1.250000019530E-04

[Giorgio Pastore]

On 4/28/12 1:49 PM, robin....@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.

Giorgio

[Giorgio Pastore]

On 4/28/12 8:45 AM, Richard Maine wrote:
....

> Why in the world would the standard call out that particular case? Darn
> near *everything* about floating-point numbers can give inconsistent
> results per the standard. Some things actuyally do give surprising
> results in actual implementations. The standard allows far worse
> surprises than you'll actually see.
>
> Note also that the standard is generally not a tutorial document. People
> who need that level of warning would probably be better directed to read
> other materials such as some of the books on the subject.

I do not agree with you.

Your point of view that "People

> who need that level of warning would probably be better directed to

> read other materials such as some of the books on the subject." is
not fully consistent with the actual content of Fortran Standards.

For example, at section 7.1.8.3 of 2003 standard, I can read
a detailed explanation of the difference between mathematical and
computational equivalence. Then, the following note 7.17 provides an
explicit example.

As far as I can see, exactly like in note 7.17 the standard explains
that "Any difference between the values of the expressions (1./3.)*3.
and 1. is a computational difference,
not a mathematical difference. The difference between the values of the
expressions 5/2 and 5./2.
is a mathematical difference, not a computational difference."
it could have explained the difference between (1.0,0.0)/(3.0,0.0) and
1./3. And I think that such example could be even more useful than the
case integer vs real, since the mathematical properties of complex
numbers with zero imaginary part are definite (from mathematicians)
fully isomorph to those of real numbers.

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.


Giorgio

[Ron Shepard]

In article <4f9bfc52$0$1385$4faf...@reader1.news.tin.it>,

Giorgio Pastore <pas...@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:

discr1 = b1**2 - 4*a1*c1
xpr1 = ( -b1 + sqrt(discr1) )/(2*a1)
xpc1 = ( -b1 + sqrt(cmplx(discr1,0.0_rk1) ) )/(2*a1)

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.

$.02 -Ron Shepard

[glen herrmannsfeldt]

You have to get used to Robin. That always happens.

-- glen

[Giorgio Pastore]

On 4/28/12 5:35 PM, Ron Shepard wrote:
>.... 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.

Giorgio

[Thomas Koenig]

John Harper <john....@vuw.ac.nz> schrieb:

> 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

[Terence]

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).

[robin....@gmail.com]

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."
>
> Robert Corbett

[robin....@gmail.com]

On Sunday, 29 April 2012 02:03:36 UTC+10, glen herrmannsfeldt wrote:
> Giorgio Pastore

> wrote:
> > On 4/28/12 1:49 PM, robin....@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.

[robin....@gmail.com]

On Sunday, 29 April 2012 22:14:14 UTC+10, robin....@gmail.com

I thought that it wasn't possible for software to be worse than
Micro$oft's, but Google has managed it.

Google posted a response that I cancelled, and failed to post a response to
another article that I posted.

[robin....@gmail.com]

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.

[robert....@oracle.com]

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.

Robert Corbett

[Thomas Jahns]

On 04/28/2012 07:39 AM, Giorgio Pastore wrote:
> 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.

Thomas Jahns

[Qolin]

robin....@gmail.com wrote in message
news:6663432.360.1335714994206.JavaMail.geo-discussion-forums@pbbpz9...

...As does that.