- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
This feature isn't called out in the list of implemented new standard features
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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 !
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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 )
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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?
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I got the same results as you with -O with cygwin64 gfortran 5.2 and 6.0.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
The missing complex-argument trigonometric functions from F2008 have been implemented for the next major release of the compiler (second half of 2017.)
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page