Cube root code example 2

cbrt_2.f90

module cbrt_mod2

implicit none
private

public :: cbrt, sp

interface cbrt
  module procedure cbrt_sp_sp
  module procedure cbrt_sp_csp
  module procedure cbrt_csp_csp
end interface

integer, parameter :: sp = kind(1.)

complex(sp), parameter :: r = cmplx(-1,sqrt(3._sp))/2._sp

contains

  pure function cbrt_sp_sp(x) result(cbrt)
    real(sp), intent(in) :: x
    real(sp) :: cbrt

    if (x >= 0.0_sp) then
      cbrt = x**(1._sp/3)
    else 
      cbrt = -(-x)**(1._sp/3)
    end if
  end function

  pure function cbrt_csp_csp(z,k) result(cbrt)
    complex(sp), intent(in) :: z
    integer, intent(in), optional :: k
    complex(sp) :: cbrt
    integer :: k_
    
    cbrt = z**(1._sp/3)
    if (present(k)) then
        k_ = mod(k,3)
        select case(k_)
        case(1)
          cbrt = cbrt * r
        case(2)
          cbrt = cbrt * conjg(r)
        end select
    end if
  end function

  pure function cbrt_sp_csp(x,k) result(cbrt)
    real(sp), intent(in) :: x
    integer, intent(in) :: k
    complex(sp) :: cbrt
    cbrt = cbrt_csp_csp(cmplx(x),k)
  end function

end module

program test_cbrt

  use cbrt_mod2, only: cbrt,sp
  implicit none

  complex(sp) :: z
  character(40) :: fmtr, fmtc

  fmtr = '(A,F7.4)'
  fmtc = '(A,F7.4,SP,F7.4,"j")'

  write(*,'(A)')  "real to real"
  
  write(*,fmtr) "cbrt( 8._sp) = ", cbrt(8._sp)
  write(*,fmtr) "cbrt(-8._sp) = ", cbrt(-8._sp)

  write(*,'(/,A)')  "real to complex"

  write(*,fmtc) "cbrt(8._sp,k=0) = ", cbrt(8._sp,k=0)
  write(*,fmtc) "cbrt(8._sp,k=1) = ", cbrt(8._sp,k=1)
  write(*,fmtc) "cbrt(8._sp,k=2) = ", cbrt(8._sp,k=2)

  write(*,'(/,A)')  "complex to complex"

  z = cmplx(-8._sp)
  write(*,fmtc) "z = ", z
  write(*,fmtc) "cbrt(z) = ", cbrt(z)
  write(*,fmtc) "cbrt(z,k=0) = ",cbrt(z,k=0)
  write(*,fmtc) "cbrt(z,k=1) = ",cbrt(z,k=1)
  write(*,fmtc) "cbrt(z,k=2) = ",cbrt(z,k=2)
  write(*,fmtc) "cbrt(z,k=3) = ",cbrt(z,k=3)

end program