@gyrtle , it was mentioned upthread you may be looking for a āone line solutionā or some such to resolve this issue. And Julia and Python based program responses were shown as well.
Well, Iām unsure what kind of suggestions you seek but note the dynamic environments such as Julia and Python have minions and minions putting in countless hours pursuing essentially a āmonotheisticā path (e.g., have no worries about consensus in a ISO / IEC standard committee where every member can in principle come from a different ideology) in various library solutions.
If you want to match the program behavior of those solutions in Fortran, you can adopt a similar path and āreinvent the wheelā entirely with complex arithmetic and pursue something like the following, for after all, thatās what a language du jour brings, yet another reinvention!
PROGRAM sZ
use complex_m
REAL(DP),PARAMETER :: r = 0.2_dp
type(complex_t(WP=DP)) :: c0
c0 = (0.0,0.5) ; print *, sqrt(c0**2 - r)
c0 = (-0.0,0.5) ; print *, sqrt(c0**2 - r)
END PROGRAM sZ
to get the output you seek:
C:\temp>ifort /standard-semantics sz.f90
Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel(R) 64, Version 2021.5.0 Build 20211109_000000
Copyright (C) 1985-2021 Intel Corporation. All rights reserved.
Microsoft (R) Incremental Linker Version 14.30.30706.0
Copyright (C) Microsoft Corporation. All rights reserved.
-out:sz.exe
-subsystem:console
sz.obj
C:\temp>sz.exe
(0.00000000000000,0.670820393249937)
(0.00000000000000,0.670820393249937)
C:\temp>
You will note the above is served by this:
Click to see the silly "reinvention of the wheel"
module complex_m
! Add support for other precision to continue the silliness
integer, parameter :: SP = kind(1.0)
integer, parameter :: DP = kind(1D0)
real(DP), parameter :: ZERO = 0.0_dp
real(DP), parameter :: ONE = 1.0_dp
real(DP), parameter :: TWO = 2.0_dp
real(DP), parameter :: EPS = epsilon(ZERO)
type :: complex_t(WP)
integer, kind :: WP = DP
real(kind=WP) :: re = ZERO
real(kind=WP) :: im = ZERO
end type
interface assignment(=)
module procedure assign_complex_dp
end interface
interface operator(**)
module procedure exp_complex_dp
end interface
interface operator(-)
module procedure subtract_complex_real_dp
end interface
interface sqrt
module procedure sqrt_complex_dp
end interface
interface write(formatted)
module procedure write_complex_dp
end interface
private :: assign_complex_dp, exp_complex_dp, subtract_complex_real_dp, sqrt_complex_dp, &
write_complex_dp
contains
subroutine assign_complex_dp( lhs, rhs )
type(complex_t(WP=DP)), intent(out) :: lhs
complex(kind=SP), intent(in) :: rhs
lhs%re = rhs%re ; lhs%im = rhs%im
end subroutine
function exp_complex_dp( c, n ) result(r)
type(complex_t(WP=DP)), intent(in) :: c
integer, intent(in) :: n
type(complex_t(WP=DP)) :: r
select case ( n )
case ( 2 )
r%re = -( c%re**2 + c%im**2 )
if ( abs(c%re) > EPS ) then
r%im = -TWO*c%re*c%im
else
r%im = ZERO
end if
case default
error stop "exponent not yet supported."
end select
end function
function subtract_complex_real_dp( c, rhs ) result(r)
type(complex_t(WP=DP)), intent(in) :: c
real(kind=DP), intent(in) :: rhs
type(complex_t(WP=DP)) :: r
r%re = c%re - rhs
end function
function sqrt_complex_dp( c ) result(r)
type(complex_t(WP=DP)), intent(in) :: c
type(complex_t(WP=DP)) :: r
if ( abs(c%im) < EPS ) then
r%re = ZERO
r%im = sqrt(abs(c%re))
else
r%re = sqrt( (sqrt(c%re**2+c%im**2) + c%re)/TWO )
r%im = sign(ONE, c%im)*sqrt( (sqrt(c%re**2+c%im**2) - c%re)/TWO )
end if
end function
subroutine write_complex_dp( dtv, lun, iotype, vlist, istat, imsg )
! Argument list
class(complex_t(WP=DP)), intent(in) :: dtv
integer, intent(in) :: lun
character(len=*), intent(in) :: iotype
integer, intent(in) :: vlist(:)
integer, intent(out) :: istat
character (len=*), intent(inout) :: imsg
! local variable
character(len=20) :: pfmt
integer :: listv
complex(kind=DP) :: c
istat = 0
c%re = dtv%re ; c%im = dtv%im
if ( iotype == "LISTDIRECTED" ) then
write(lun, fmt=*, iostat=istat, iomsg=imsg) c
else if ( iotype(1:2) == "DT" ) then
! vlist(1) is to be used as the field widths of the
! component of the derived type variable. First set up the format to
! be used for output.
if ( size(vlist) > 0 ) then
if ( vlist(1) > 0 ) listv = vlist(1)
end if
write(pfmt,"(*(g0))" ) "(g", listv, ")"
write(lun, fmt=pfmt, advance="no", iostat=istat, iomsg=imsg) c
else
! Not supported
istat = 1
imsg = "Namelist option is not yet supported."
end if
return
end subroutine write_complex_dp
end module
P.S.> Apologies for any errors, it was slapped together rather hastily
