Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.

Computing Arcsine of a complex do not compile

vincent_m_
Beginner
722 Views

Hello,

I have bought Intel(R) Visual Fortran Compiler XE 15.0.0.108 [IA-32] and Visual Studio 2012. But I have a problem with a code where I need to compute the Arcsine of a complex number, although it is in the Fortran 2008 standard (final draft) :

13.7.14 ASIN (X)
1 Description. Arcsine (inverse sine) function.
2 Class. Elemental function.
3 Argument. X shall be of type real with a value that satisfies the inequality |X| ≤ 1, or of type complex.

For example, the following code compiles and runs with gfortran (under Linux) but not with Visual Fortran Compiler XE 15.0.0.108 (Windows 8.1) :

complex(8) x, y
x = (0, 1)
y = asin(x)
print *, y
end

I obtain the following warning and error :
warning #7319: This argument's data type is incompatible with this intrinsic procedure; procedure assumed EXTERNAL.   [ASIN]
error #6404: This name does not have a type, and must have an explicit type.   [ASIN]

Is this a bug or is there an option to activate in Visual Studio ?

Thanks
Vincent

0 Kudos
16 Replies
Simon_Geard
New Contributor I
722 Views

Strange - presumably it's not yet implemented. I tried  Version 16.0.0.110 and got the same result. The following works though:

	complex function arcsin(z)
		complex, intent(in) :: z
		complex, parameter :: i = cmplx(0,1)
		arcsin = -i*log(i*z+sqrt(1-z**2))
	end function arcsin

 

0 Kudos
TimP
Honored Contributor III
722 Views

This feature isn't called out in the list of implemented new standard features

https://software.intel.com/en-us/articles/intel-fortran-compiler-support-for-fortran-language-standards

0 Kudos
mecej4
Honored Contributor III
722 Views

If you look up the description of the function ASIN in the IFort 16.0 manual, you see that its argument must satisfy "x   (Input)   Must be of type real. The | x | must be less than or equal to 1.".

There are many features of F2008 that are not yet implemented in IFort.

0 Kudos
vincent_m_
Beginner
722 Views

Thank you for your quick replies !

Although it is easy to compute the function (thank you sgeard), I am quite disappointed to see that such a simple feature of Fortran 2008 is not yet implemented in a commercial compiler. I can understand that implementing fully things like coarrays and other parallel stuff can take time, but basic mathematical functions like that...

Anyway, I have got the answer to my problem !

0 Kudos
Steven_L_Intel1
Employee
722 Views

This is on our list of F2008 features to be implemented - issue ID is DPD200378901. I note that we do support TAN here but not the other intrinsics for which COMPLEX support was added in F2008.

0 Kudos
vincent_m_
Beginner
722 Views

Steve Lionel (Intel) wrote:

This is on our list of F2008 features to be implemented - issue ID is DPD200378901. I note that we do support TAN here but not the other intrinsics for which COMPLEX support was added in F2008.

Thank you for this answer. It is the first time in my carreer I write a model needing ASIN(complex), but it can be useful.

0 Kudos
mecej4
Honored Contributor III
722 Views

Before the 1980s it was common for scientists and engineers to have within reach a handbook/pocketbook of mathematical formulae and function tables. See, for example:

     https://ia802604.us.archive.org/30/items/ashorttableinte00fostgoog/ashorttableinte00fostgoog.pdf

One relevant formula is No. 18. on p. 4.

0 Kudos
vincent_m_
Beginner
722 Views

mecej4 wrote:

Before the 1980s it was common for scientists and engineers to have within reach a handbook/pocketbook of mathematical formulae and function tables.

Of course we can compute it ourself, but Fortran being a language for intensive computations, an interesting question is: how a function like ASIN(complex) is implemented in a compiler ? Is it simply the computation of -i*log(i*z+sqrt(1-z**2)) which uses three multiplications, two additions/substractions, one log and one sqrt with complex values ? Or is it implemented in another way to be faster to compute and to limit rounding errors ? I am curious about that.

 

 

0 Kudos
Steven_L_Intel1
Employee
722 Views

The math library may have multiple algorithms it uses, depending on instruction set and choice of accuracy. Typically the algorithms are not the simple "canonical" ones but are indeed structured to be faster and limit rounding, out-of-range and cancellation errors.

The compiler itself rarely has such things - it calls into a math library. There may be several versions of a particular function depending if vectorization is in use or not.

In the case at hand, we believe that the underlying math library already has the necessary support, so what we need to do is have the compiler recognize these combinations and "hook up" to an appropriate math library call.

0 Kudos
JVanB
Valued Contributor II
722 Views

There is also the question of whether a replacement expression gets the right values both with positive and negative zeros on a branch cut. An example:

program P
   use ISO_FORTRAN_ENV
   implicit none
   real(REAL32), parameter :: two = real(z'40000000',REAL32)
   real(REAL32), parameter :: positive = real(z'00000000',REAL32)
   real(REAL32), parameter :: negative = real(z'80000000',REAL32)
   complex(REAL32) z

   z = cmplx(two,positive,REAL32)
   write(*,*) 'Above positive real axis:'
   write(*,*) asin(z)
   write(*,*) -(0,1)*log((0,1)*z+sqrt(1-z**2))
   z = cmplx(two,negative,REAL32)
   write(*,*) 'Below positive real axis:'
   write(*,*) asin(z)
   write(*,*) -(0,1)*log((0,1)*z+sqrt(1-z**2))
   z = cmplx(-two,positive,REAL32)
   write(*,*) 'Above negative real axis:'
   write(*,*) asin(z)
   write(*,*) -(0,1)*log((0,1)*z+sqrt(1-z**2))
   z = cmplx(-two,negative,REAL32)
   write(*,*) 'Below negative real axis:'
   write(*,*) asin(z)
   write(*,*) -(0,1)*log((0,1)*z+sqrt(1-z**2))
end program P

Output with gfortran:

 Above positive real axis:
 (  1.57079637    ,  1.31695783    )
 (  1.57079637    , -1.31695795    )
 Below positive real axis:
 (  1.57079637    , -1.31695795    )
 (  1.57079637    , -1.31695795    )
 Above negative real axis:
 ( -1.57079637    ,  1.31695783    )
 ( -1.57079637    ,  1.31695783    )
 Below negative real axis:
 ( -1.57079637    , -1.31695795    )
 ( -1.57079637    ,  1.31695783    )

So there are some inconsistencies observed in this case.

 

0 Kudos
TimP
Honored Contributor III
722 Views

While it was painful to edit this into compilable form, gfortran requires -O to produce other than

 Above positive real axis:
 (  1.57079637    , -1.31695795    )
 (  1.57079637    , -1.31695795    )
 Below positive real axis:
 (  1.57079637    , -1.31695795    )
 (  1.57079637    , -1.31695795    )
 Above negative real axis:
 ( -1.57079637    ,  1.31695783    )
 ( -1.57079637    ,  1.31695783    )
 Below negative real axis:
 ( -1.57079637    ,  1.31695783    )
 ( -1.57079637    ,  1.31695783    )

 

0 Kudos
JVanB
Valued Contributor II
722 Views

That's a bit surprising since what I posted compiled for me with gfortran 5.2.0 on Windows and always gave the results I showed with

gfortran asin.f90 -oasin

or

gfortran -O asin.f90 -oasin

both 32 and 64 bits. Are you compiling with some other version or in cygwin?

 

 

0 Kudos
JVanB
Valued Contributor II
722 Views

BTW, I tried coding up the recommendation of https://people.freebsd.org/~das/kahan86branch.pdf , but it didn't seem to help:

program P
   use ISO_FORTRAN_ENV
   implicit none
   real(REAL32), parameter :: two = real(z'40000000',REAL32)
   real(REAL32), parameter :: positive = real(z'00000000',REAL32)
   real(REAL32), parameter :: negative = real(z'80000000',REAL32)
   complex(REAL32) z

   z = cmplx(two,positive,REAL32)
   call S('Above positive real axis:',z)
   z = cmplx(two,negative,REAL32)
   call S('Below positive real axis:',z)
   z = cmplx(-two,positive,REAL32)
   call S('Above negative real axis:',z)
   z = cmplx(-two,negative,REAL32)
   call S('Below negative real axis:',z)
end program P

subroutine S(c,z)
   use ISO_FORTRAN_ENV
   implicit none
   character(*) c
   complex(REAL32) z

   write(*,*) c
   write(*,*) asin(z)
   write(*,*) cmplx(atan2(real(z),real(sqrt(1-z)*sqrt(1+z))), &
      asinh(imag(conjg(sqrt(1-z))*sqrt(1+z))),REAL32)
   write(*,*) -(0,1)*log((0,1)*z+sqrt(1-z**2))
end subroutine S

Results with gfortran:

 Above positive real axis:
 (  1.57079637    ,  1.31695783    )
 (  1.57079637    , -1.31695783    )
 (  1.57079637    , -1.31695795    )
 Below positive real axis:
 (  1.57079637    , -1.31695795    )
 (  1.57079637    , -1.31695783    )
 (  1.57079637    , -1.31695795    )
 Above negative real axis:
 ( -1.57079637    ,  1.31695783    )
 ( -1.57079637    ,  1.31695783    )
 ( -1.57079637    ,  1.31695783    )
 Below negative real axis:
 ( -1.57079637    , -1.31695795    )
 ( -1.57079637    ,  1.31695783    )
 ( -1.57079637    ,  1.31695783    )

Edit: fixed code and output.

0 Kudos
TimP
Honored Contributor III
722 Views

I got the same results as you with -O with cygwin64 gfortran 5.2 and 6.0.

0 Kudos
Steven_L_Intel1
Employee
722 Views

The missing complex-argument trigonometric functions from F2008 have been implemented for the next major release of the compiler (second half of 2017.)

0 Kudos
AShte
Beginner
723 Views
The correct answers, of course, are (floating point precision issue aside): z: ( 2.00000000 , 0.00000000 ) ASIN(z): ( 1.57079637 , 1.31695795 ) z: ( 2.00000000 , -0.00000000 ) ASIN(z): ( 1.57079637 , -1.31695795 ) z: ( -2.00000000 , 0.00000000 ) ASIN(z): ( -1.57079637 , 1.31695795 ) z: ( -2.00000000 , -0.00000000 ) ASIN(z): ( -1.57079637 , -1.31695795 ) program p use, intrinsic :: iso_fortran_env use, intrinsic :: ieee_arithmetic implicit none integer, parameter :: rk = real32 real( kind=rk ), parameter :: two = 2_rk, positive = 0.0_rk, & negative = -0.0_rk complex( kind=rk ) :: z z = cmplx(two,positive, kind=rk ) write(*,*) 'z:', z write (*,*) "ASIN(z):", asin(z) z = cmplx(two,negative, kind=rk ) write(*,*) 'z:', z write (*,*) "ASIN(z):", asin(z) z = cmplx(-two,positive, kind=rk ) write(*,*) 'z:', z write (*,*) "ASIN(z):", asin(z) z = cmplx(-two,negative, kind=rk ) write(*,*) 'z:', z write (*,*) "ASIN(z):", asin(z) end program p Compiled with gfortran6 -O3 -Wall -fsign-zero -Wl,-rpath="/usr/local/lib/gcc6" z.f90 on FreeBSD 11-current amd64. Note that -O3 optimisation preserves minus zero. Anton
0 Kudos
Reply