Hi all
How can we handle “Error: Raising a negative REAL at (1) to a REAL power is prohibited” in example below:
(-8) ** (1.0 / 3.0)
1 Like
Use complex numbers ?
With real numbers, you can not compute (-8)^{1/2} => square root of a negative number
Same (mathematical) problem with 1/3.
1 Like
@ELNS
See: https://math.stackexchange.com/questions/317528/how-do-you-compute-negative-numbers-to-fractional-powers
1 Like
@ELNS
There is some subtle things there, more subtle than it seems 
:
- you can compute \sqrt[3]{-27}: it’s -3, because -3 \times (-3) \times (-3)=-27
- But you can not compute (-27)^{1/3}
In fact, the mathematical function \sqrt[3]{x} = x^{1/3} when x \gt 0.
As far as I know, there is no cubic root function in Fortran
=> probably an interesting exercise…
1 Like
Thank you so much for helping me 
1 Like
I’m confused. Where is this prohibited? Not in Fortran, it’s not. You can do it but you won’t get a meaningful answer:
$ cat test_exception1.f90
program test_exception
implicit none
integer :: n
real :: x, y
do n = 1, 10
call random_number(x) ! random number in [0, 1]
x = x - 0.5 ! random number in [-0.5, 0.5]
y = x**(1./3.)
print *, 'x, y = ', x, y
end do
end program test_exception
$ gfortran -Wall -std=f2018 test_exception1.f90 && ./a.out
x, y = 8.94961357E-02 0.447302580
x, y = 0.190010607 0.574900389
x, y = 0.334387779 0.694091618
x, y = 6.04379177E-03 0.182153061
x, y = -0.412383974 NaN
x, y = 0.465656042 0.775095224
x, y = -0.412225127 NaN
x, y = 0.237244248 0.619058788
x, y = -0.436005294 NaN
x, y = -6.18785620E-02 NaN
This can happen in real-world applications–you need to work with real numbers and you can’t always predict that they’ll be positive. But you’re right, you should be able to handle it. Here’s how you can do it for this specific case:
$ cat test_exception2.f90
program test_exception
use ieee_arithmetic, only: ieee_is_nan
implicit none
integer :: n
real :: x, y
do n = 1, 10
call random_number(x) ! random number in [0, 1]
x = x - 0.5 ! random number in [-0.5, 0.5]
y = x**(1./3.)
if (ieee_is_nan(y)) then ! catch exception
print *, 'Exception caught in y = x**(1./3.), x, y = ', x, y
! handle exception here however you want
else
print *, 'x, y = ', x, y
end if
end do
end program test_exception
$ gfortran -Wall -std=f2018 test_exception2.f90 && ./a.out
x, y = 0.381638110 0.725354970
x, y = 1.73110366E-02 0.258686841
Exception caught in y = x**(1./3.), x, y = -0.201909781 NaN
Exception caught in y = x**(1./3.), x, y = -0.197003841 NaN
Exception caught in y = x**(1./3.), x, y = -0.303777575 NaN
x, y = 0.129012465 0.505293727
Exception caught in y = x**(1./3.), x, y = -0.205134809 NaN
x, y = 0.403641343 0.739035368
x, y = 0.304455400 0.672730684
x, y = 5.94465733E-02 0.390279382
1 Like
This is not a Fortran problem. But a mathematical one:
Real exponents are mathematically an extension of integer exponents. If a is a real, x^a can be defined by e^{a \cdot ln(x)}, and therefore it is not defined if x \le 0, unless you work in \mathbb{C}.
See:
1 Like
Okay, I just enabled a plugin for that. Here’s an example of inline and block math:
Inline: $\sqrt{x^2 + y^2}$
Inline: \sqrt{x^2 + y^2}
Block:
$$
\sqrt{x^2 + y^2}
$$
Block:
\sqrt{x^2 + y^2}
3 Likes
Perhaps, but I didn’t get that idea from the original post. 
Using complex variables seems to work here…
program main
use iso_fortran_env, only: dp => real64
implicit none
print *, cmplx( -8.0, 0.0 )**( 1.0 / 3.0 )
print *, ( cmplx( -8.0, 0.0 )**( 1.0 / 3.0 ) )**3
print *
print *, cmplx( -8.0, -0.0 )**( 1.0 / 3.0 )
print *, ( cmplx( -8.0, -0.0 )**( 1.0 / 3.0 ) )**3
print *
print *, cmplx( -8.0_dp, 0.0_dp, dp )**( 1.0_dp / 3.0_dp )
print *, ( cmplx( -8.0_dp, 0.0_dp, dp )**( 1.0_dp / 3.0_dp ) )**3
print *, cmplx( -8.0_dp, -0.0_dp, dp )**( 1.0_dp / 3.0_dp )
print *, ( cmplx( -8.0_dp, -0.0_dp, dp )**( 1.0_dp / 3.0_dp ) )**3
end
Result (gfortran-9.3):
(0.999999940,1.73205090)
(-8.00000095,-1.148161914E-06)
(0.999999940,-1.73205090)
(-8.00000095,1.148161914E-06)
(1.0000000000000000,1.7320508075688772)
(-7.9999999999999991,6.02105053308414143E-016)
(1.0000000000000000,-1.7320508075688772)
(-7.9999999999999991,-6.02105053308414143E-016)
Interestingly, Python gives a complex result from (-8)**(1.0 / 3.0) (w/o making -8 to complex),
>>> (-8)**(1.0 / 3.0)
(1.0000000000000002+1.7320508075688772j)
>>> (-8.0 + 0.0j)**(1.0 / 3.0)
(1.0000000000000002+1.7320508075688772j)
>>> (-8.0 - 0.0j)**(1.0 / 3.0)
(1.0000000000000002-1.7320508075688772j)
while Julia refuses it and requires a complex argument (so similar to Fortran):
> (-8.0)^(1.0 / 3.0)
**ERROR:** DomainError with -8.0:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
> (-8.0 + 0.0im)^(1.0 / 3.0)
1.0 + 1.7320508075688772im
> (-8.0 - 0.0im)^(1.0 / 3.0)
1.0 - 1.7320508075688772im
3 Likes
Perhaps, but I didn’t get that idea from the original post.
I agree with you, the question was:
so handling the exception is the exact answer.
I was obsessed by the mathematical aspect…
2 Likes
Fortran 2018, 10.1.5.2.4 Evaluation of numeric intrinsic operations, paragraph 1:
Raising a negative real value to a real power is prohibited.
6 Likes
I had no idea! Is this new in F2018?
Does it mean that a standard-conforming compiler must check at run-time that such operation doesn’t occur, like in my example above? gfortran-9.2.0 doesn’t do this, even with the -std=f2018 flag.
1 Like
No - it’s even in FORTRAN 66 (6.4):
No factor may be evaluated that requires a negative valued primary to be raised to a real or double precision exponent.
The standard NEVER requires a “processor” to do run-time checks, and this isn’t a numbered syntax rule or constraint that requires the ability to diagnose.
4 Likes
The thing to always keep in mind is that the Fortran standard describes the behavior of a standard-conforming program. If you color outside the lines, the standard does not specify the behavior. This gives a compiler great freedom to support extensions and behaviors such as giving you a NaN where the floating point model indicates is appropriate.
3 Likes
Thank you, Steve, this is very interesting. I think I got stuck on the word “prohibited”. So, this statement:
Raising a negative real value to a real power is prohibited.
can be interpreted as:
The result of raising a negative real value to a real power is undefined.
Do you agree?
Second, consider this program:
program schrodinger
character(10) :: arg
real :: x
call get_command_argument(1, arg)
read(arg, *) x
print *, x, x**(1./3.)
end program schrodinger
Is this program standard-conforming?
Based on my interpretation of the quoted rule, my conclusion is that the program is standard-conforming if the user inputs a positive value, and not if the user inputs a negative value. Until then, it’s both at the same time.
I bet there is a flaw somewhere in my thought process here, but I can’t find it. What do you think?
1 Like
No, I don’t agree with your first statement. The standard makes a distinction between behaviors that are prohibited (in that if you do so, your program is non-conforming), and behaviors whose effect is “processor-dependent” (the program is conforming, but the standard does not tell you what happens.)
Here’s a real-world example where the standard changed wording but not behavior. In Fortran 2008, 9.5.4 said:
A unit shall not be connected to more than one file at the same time, and a file shall not be connected to more than one unit at the same time.
Fortran 2018 changed that to (12.5.4):
A unit shall not be connected to more than one file at the same time. … It is processor dependent whether a file can be connected to more than one unit at the same time.
What is the practical difference here? A program that opened a file on more than one unit was non-conforming in Fortran 2008, but in Fortran 2018 that behavior was conforming. Compilers didn’t have to change - some already allowed it, some didn’t - both conformed to the standard in this regard (because it is describing what the program may or may not do), but programs that did this were no longer deemed nonconforming. This mattered to some users.
Your program is non-conforming if a negative value is input. It’s not “both at the same time”. There are many behaviors prohibited by the standard (deallocating an object that is not allocated, for example) which are entirely based on program flow. Whether a program is standard-conforming requires you to consider the totality of it, both in source form and execution behavior. Compilers are required to have the capability of diagnosing certain source form issues (numbered syntax rules and constraints), but the lack of a complaint doesn’t paste a gold star on your program’s forehead.
5 Likes
Okay, thank you for clarifying that, I think I understand now. I just learned something very important–you can’t determine standard-conformance of a program solely based on its source.
4 Likes
That is a fascinating way to think about standards compliance. I did not realize that a program’s compliance could be based on its input. Perhaps that is what the standards developers intended, but not at all the way I would have initially thought about it.
This leads me to an interesting question. Is it possible to prove that a program is standards compliant for all possible inputs?
Also an interesting implication, it would be possible to write esoteric programs that were only standards compliant on the third Thursday of the 4th month of every fifth leap year or some such strange case.
2 Likes
Hi,
This is a nice example how to catch exceptions in fortran.
Although, the question was about “Error Handling”, I think it is better to test if x is negative or positive or zero. It will be more clear about the “mathematical” aspect.
2 Likes