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

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

On Mar 25, 6:41=A0pm, 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

On Mar 25, 10:47=A0pm, Louis Krupp <lkr...@nospam.pssw.com.invalid>
wrote:
> On Sun, 25 Mar 2012 18:41:52 -0700 (PDT), robert.corb...@oracle.com
> wrote:
>
>
>
> >I recently purchased a copy of the book "Methods of Numerical
> >Integration" by Ralston and Rabinowitz. =A0The edition I purchased is
> >the Dover Publications reprint dated 2007. =A0It is a reprint of the
> >second edition, and it is copyright 1984.
> <snip>
> > =A0 =A0 "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?

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

robert.corbett@oracle.com wrote:

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

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

John Harper <john.harper@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

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.

Giorgio Pastore <pastgio@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

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

On 4/28/12 1:18 AM, Steven G. Kargl wrote:
> 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.
>
I do not see any bug in my program. If you find one, I would appreciate 
very much to know where.

Giorgio Pastore

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

> On 4/28/12 1:18 AM, Steven G. Kargl wrote:
>> 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.
>>
> 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)

-- 
steve

Richard Maine <nospam@see.signature> wrote:
> Giorgio Pastore <pastgio@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

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

Steven G. Kargl <sgk@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

On Sat, 28 Apr 2012 02:18:43 +0200, Giorgio Pastore wrote:

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

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

On 04/27/2012 06:01 PM, Giorgio Pastore wrote:
> 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.

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.

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

> Richard Maine <nospam@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

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

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

Giorgio

Giorgio Pastore <pastgio@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

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

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

On Saturday, 28 April 2012 08:01:42 UTC+10, Giorgio Pastore  wrote:

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

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




On Saturday, 28 April 2012 08:01:42 UTC+10, Giorgio Pastore  wrote:
> 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.

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.

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  

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

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

In article <4f9bfc52$0$1385$4fafbaef@reader1.news.tin.it>,
 Giorgio Pastore <pastgio@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

Giorgio Pastore <pastgio@units.it> wrote:
> On 4/28/12 1:49 PM, robin.vowels@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.

-- glen

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

John Harper <john.harper@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

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

On Sunday, 29 April 2012 00:14:56 UTC+10, Giorgio Pastore  wrote:
> On 4/28/12 1:49 PM, robin.vowels@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 Sunday, 29 April 2012 02:03:36 UTC+10, glen herrmannsfeldt  wrote:
> Giorgio Pastore 
>  wrote:
> > On 4/28/12 1:49 PM, robin.vowels@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 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.

On Apr 27, 10:39=A0pm, Giorgio Pastore <past...@units.it> 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. =A0On my i686-*-freebsd system,
>
> Now I can get a consistent result like you.
>
> > =A0 =A0 1.0100000 =A0 =A0 =A0is the real part of =A0( =A01.0100000 =A0 =
=A0, =A00.0000000 =A0 =A0)
> > =A0 =A0 2000.0100 =A0 =A0 =A0is the real part of =A0( =A02000.0100 =A0 =
=A0, =A00.0000000 =A0 =A0)
> > =A0 =A0real =A0 =A0math result =A0 -4.0000000
> > =A0 =A0complex math result =A0( -4.0000000 =A0 =A0, =A00.0000000 =A0 =
=A0)
>
> 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 =A0is the real part of 1.010000 0.000000
> 2000.010010 =A0is the real part of 2000.010010 0.000000
> =A0 real =A0-4.040000
> =A0 complex =A0-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.

Robert Corbett

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

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

> On Sunday, 29 April 2012 02:03:36 UTC+10, glen herrmannsfeldt  wrote:
> > Giorgio Pastore
> >  wrote:
> > > On 4/28/12 1:49 PM, robin.vowels@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.

....As does that.