Error Handling for negative real numbers

Help
21

I have no idea of the internal functioning of a compiler. Perhaps @certik could tell us how the ** operator is (or will be) implemented in LFortran (or LLVM ?) ?

Does the compiler first have a look to the exponent type ? Which verifications are made ?
For real exponents, does the compiler (or its backend) use assembler instructions like F2XM1 and FYL2X ?
https://www.felixcloutier.com/x86/f2xm1
https://www.felixcloutier.com/x86/fyl2x

F2XM1: Computes the exponential value of 2 to the power of the source operand minus 1.
Values other than 2 can be exponentiated using the following formula: x^y ← 2^{y \times log_2{x}}

FYL2X: Computes ST(1) \times log_2 (ST(0))

Are these operations described in the IEEE 754 standard ?

22

Steven replied with great answers to the question. I would just add that use of English varies between computer language standards. Fortran has a concept of a variable being “undefined”, but the concept of a behavior being “undefined” is a terminology used by the C standard, not the Fortran standard. Fortran tries to stick to the ISO-recommended definitions for words like shall, may, might, can, … (sometimes followed by “not”).

3 Likes
23

Compilers typically call a library routine to do real exponentiation, and the routine can check for things such as a negative operand. They certainly know the type and kind of the operands, and can do some optimizations such as doing a SQRT for x**0.5. Integer powers are often converted into a series of multiplies, and if there is no special case then a library routine is called.

Library routines tend to be very careful in their arithmetic and often avoid using certain instructions if it can’t be proven that they would provide the mathematically-best result. I’ll also note that math library routines are often shared across languages - it’s not unusual for Fortran and C to use the same library routines for most operations. (There are some where language differences require distinct routines.)

3 Likes
24

I agree, and based on the development in this thread, that’s the only way to ensure the power operation is standard-conforming–test the input, rather than the result.

1 Like
25

With now 24 replies, the seemingly simple question of our friend @ENLS has initiated the longest discussion yet on this discourse!

4 Likes
26

your welcome :grinning:

2 Likes
27

I’ll just throw out here that back in 1999 I wrote Doctor Fortran and “The Dog that Did Not Bark”, about what makes a program standard-conforming or not, and that what the standard doesn’t say can bite you.

5 Likes
28

Thank you for this interesting read! It raised, however, a couple of questions, that never had I considered before:

  1. Would -std=f2018 enforce strict compliance with the Fortran 2018 standard?
  2. If no, what does this flag is for?
  3. If yes, do issues like the one you mention in your article still exist or have been addressed? I quote:

there is another class of non-conformity that, in general, can’t be detected at compile time

  1. I believe std=f2018 is gfortran’s flag, is there something similar for Intel?
1 Like
29

I can’t speak for gfortran - I think its handling of such things is different than ifort’s.

ifort has -stand f08. This turns on compile-time warnings of non-standard syntax, violations of constraints, and violation of many (but not all) of the un-numbered rules. If you want these to be errors instead of warnings, add -warn stderrors. The default if you just say -stand is f08. The current ifort has a f18 option.

If you turn on standards checking, things I mention in the article such as INTEGER*8 and use of logical operators for integers, will get flagged. It won’t detect assumptions about short-circuit evaluation nor, really, any run-time behavior that falls afoul of the standard (except that things that cause an error will be reported as such.)

I also want to repeat something I wrote earlier - just because you enabled standards checking and didn’t get any complaints, that doesn’t mean your program is standard-conforming. Certainly it has a better shot at being so than if the compiler DOES complain, but there are many things a compiler is not required to check or can’t check. It’s also possible (and I have seen cases of this in all the compilers I have worked on) that it fails to notice non-standard syntax or relationships that the standard said it should - those are compiler bugs.

4 Likes
30

A possible (naive) implementation of the cube root function could be:

module functions
    use iso_fortran_env, only: dp=>real64
    implicit none

    contains

    pure real(dp) function cbrt(x)
        real(dp), intent(in) :: x
        cbrt = sign(abs(x)**(1.0_dp / 3.0_dp), x)
    end function
end module functions

Perhaps an implementation with some if...else if... could be faster if speed matters. I don’t know.

1 Like
32

Are these procedures (cbrt,cbrtl) generally considered pure, meaning we can wrap them into an elemental function?

34

A good reference:

What Every Computer Scientist Should Know About Floating-Point Arithmetic, by David Goldberg, published in the March, 1991 issue of Computing Surveys.
https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

Search the “exponentiation” word in the page, there are some interesting paragraphs. For example:

The programming language Ada avoids this problem by only defining exponentiation for integer powers, while ANSI FORTRAN prohibits raising a negative number to a real power.
In fact, the FORTRAN standard says that
Any arithmetic operation whose result is not mathematically defined is prohibited…

1 Like
35

There’s no need for a wrapper - just define a PURE interface for the functions.

37

Here is an implementation of cbrt in openlibm: https://github.com/JuliaMath/openlibm/blob/97de1a46a9fedec7749e10da6770d189ff8012e9/src/s_cbrt.c#L38

They use a polynomial approximation and a Newton iteration to obtain a highly accurate answer and a high performance.

1 Like
39

@kargl thanks for the information!