A postfix (a.k.a. Reverse-Polish Notation - RPN) calculator

Question

As an exercise, I put together a postfix calculator using modern Fortran. Language apart, I am interested in knowing your take on the algorithm. As far as I remember from my freshman year (chemistry - long ago), the problem has a standard solution in C, which I imagine is optimal in some sense. However, I did not look it up, and wrote something that is probably different in some respects. The program runs and passes the tests.

I am interested in knowing whether the present solution is acceptable, or if it has any major hidden flaws / inefficiencies. For folks not familiar with the simple syntax of modern Fortran, I suggest the following quick modern Fortran tutorial.

Thanks!

module mod_postfix

  implicit none

  private
  integer, parameter :: TOKEN_MAX_LEN = 50

  public :: EvalPostfix

contains

  real(kind(1d0)) function EvalPostfix( CmdStrn ) result(res)
    character(len=*), intent(in) :: CmdStrn
    integer :: iToken, nTokens, shift
    character(len=:)            , allocatable :: Token
    character(len=TOKEN_MAX_LEN), allocatable :: stack(:)
    nTokens = GetNTokens(CmdStrn)
    allocate(stack(nTokens))
    do iToken = 1, nTokens
       call GetToken(CmdStrn,iToken,Token)
       stack(iToken) = Token
    enddo
    shift=0
    call simplify_stack(nTokens,Stack,shift)
    read(Stack(nTokens),*)res
  end function EvalPostfix

  recursive subroutine simplify_stack(n,Stack,shift) 
    integer                     , intent(in)    :: n
    character(len=TOKEN_MAX_LEN), intent(inout) :: Stack(:)
    integer                     , intent(inout) :: shift
    character(len=:), allocatable :: sOp
    integer         :: i
    real(kind(1d0)) :: v1, v2, res
    logical :: IsBinary, IsUnary, IsNonary, IsOperator
    
    if(n==0)return

    sOp = trim(Stack(n))

    !.. Case Binary Operators
    IsBinary   = index( "+ - * / max min mod **", sOp ) > 0
    IsUnary    = index( " sin cos tan asin acos atan exp log int sqrt abs", sOp ) > 0
    IsNonary   = index( " random_number PI", sOp ) > 0
    IsOperator = IsBinary .or. IsUnary .or. IsNonary
    
    if( ( .not. IsOperator ) .and. n == shift + 1 )return

    call simplify_stack(n-1,stack,shift)

    if( IsBinary )then
       
       read(Stack(n-1),*)v2
       read(Stack(n-2),*)v1
       if( sOp == "+"    ) res = v1+v2
       if( sOp == "-"    ) res = v1-v2
       if( sOp == "*"    ) res = v1*v2
       if( sOp == "/"    ) res = v1/v2
       if( sOp == "max"  ) res = max(v1,v2)
       if( sOp == "min"  ) res = min(v1,v2)
       if( sOp == "mod"  ) res = mod(v1,v2)
       if( sOp == "**"   ) res = v1**v2
       write(Stack(n),"(e24.16)")res
       shift=shift+2
       do i=n-3,1,-1
          Stack(i+2)=Stack(i)
       enddo

    elseif( IsUnary )then

       read(Stack(n-1),*)v1
       if( sOp == "sin" ) res = sin (v1)
       if( sOp == "cos" ) res = cos (v1)
       if( sOp == "tan" ) res = tan (v1)
       if( sOp == "asin") res = asin(v1)
       if( sOp == "acos") res = acos(v1)
       if( sOp == "atan") res = atan(v1)
       if( sOp == "exp" ) res = exp (v1)
       if( sOp == "log" ) res = log (v1)
       if( sOp == "sqrt") res = sqrt(v1)
       if( sOp == "abs" ) res = abs (v1)
       if( sOp == "int" ) res = dble(int(v1))
       write(Stack(n),"(e24.16)")res
       shift=shift+1
       do i=n-2,1,-1
          Stack(i+1)=Stack(i)
       enddo

    elseif( IsNonary )then

       if( sOp == "random_number")call random_number(res)
       if( sOp == "PI"           )res=4.d0*atan(1.d0)
       write(Stack(n),"(e24.16)")res
       if(n == shift + 1)return
       call simplify_stack(n-1,stack,shift)

    end if

  end subroutine simplify_stack


  !> Counts the number of tokens 
  integer function GetNTokens( strn, separator_list_ ) result( n )
    implicit none
    character(len=*)          , intent(in) :: strn
    character(len=*), optional, intent(in) :: separator_list_
    !
    character       , parameter   :: SEPARATOR_LIST_DEFAULT = " "
    character(len=:), allocatable :: separator_list
    integer :: i,j
    n=0
    if(len_trim( strn ) == 0)return
    if(present(separator_list_))then
       allocate(separator_list,source=separator_list_)
    else
       allocate(separator_list,source=SEPARATOR_LIST_DEFAULT)
    endif
    i=1
    do
       j=verify(strn(i:),separator_list)
       if(j<=0)exit
       n=n+1
       j=i-1+j
       i=scan(strn(j:),separator_list)
       if(i<=0)exit
       i=j-1+i
       if(i>len(strn))exit
    enddo
    if(allocated(separator_list))deallocate(separator_list)
  end function GetNTokens

  subroutine GetToken( strn, iToken, token, separator_list_ )
    implicit none
    character(len=*),              intent(in) :: strn
    integer         ,              intent(in) :: iToken
    character(len=:), allocatable, intent(out):: token
    character(len=*), optional   , intent(in) :: separator_list_
    !
    character       , parameter   :: SEPARATOR_LIST_DEFAULT = " "
    character(len=:), allocatable :: separator_list
    integer :: i,j,n
    if(present(separator_list_))then
       allocate(separator_list,source=separator_list_)
    else
       allocate(separator_list,source=SEPARATOR_LIST_DEFAULT)
    endif
    if(iToken<1)return
    if(iToken>GetNTokens(strn,separator_list))return
    if(allocated(token))deallocate(token)
    i=1
    n=0
    do 
       j=verify(strn(i:),separator_list)
       if(j<=0)exit
       n=n+1
       j=i-1+j
       i=scan(strn(j:),separator_list)
       if(i<=0)then
          i=len_trim(strn)+1
       else
          i=j-1+i
       endif
       if(n==iToken)then
          allocate(token,source=strn(j:i-1))
          exit
       endif
    enddo
  end subroutine GetToken

end module Mod_Postfix


program TestPostfixCalculator
  use mod_postfix
  implicit none
  real(kind(1d0)) , parameter   :: THRESHOLD = 1.d-10
  real(kind(1d0))               :: res
  character(len=:), allocatable :: sPostfix
  
  call assert("+"  , abs( EvalPostfix(" 3 4 +")   -  7    ) < THRESHOLD )
  call assert("-"  , abs( EvalPostfix(" 3 4 -")   +  1    ) < THRESHOLD )
  call assert("*"  , abs( EvalPostfix(" 3 4 *")   - 12    ) < THRESHOLD )
  call assert("/"  , abs( EvalPostfix(" 3 4 /")   -  0.75 ) < THRESHOLD )
  call assert("max", abs( EvalPostfix(" 3 4 max") -  4    ) < THRESHOLD )
  call assert("min", abs( EvalPostfix(" 3 4 min") -  3    ) < THRESHOLD )
  call assert("mod", abs( EvalPostfix("13 5 mod") -  3    ) < THRESHOLD )
  call assert("**" , abs( EvalPostfix(" 2 5 **" ) - 32    ) < THRESHOLD )

  call assert("cos", abs( EvalPostfix(" PI 3 / cos" ) - 0.5 ) < THRESHOLD )

  res      = sqrt( (log(10.d0)-atan(2.d0))/max(cos(6.d0),exp(3.d0)) )
  sPostfix ="10 log 2 atan - 6 cos 3 exp max / sqrt" 
  call assert("expression1", abs( EvalPostfix(sPostfix) - res ) < THRESHOLD )
  !.. etc. etc.
  
contains
  
  subroutine assert(msg,cond)
    use, intrinsic :: iso_fortran_env, only : OUTPUT_UNIT
    character(len=*), intent(in) :: msg
    logical         , intent(in) :: cond
    write(OUTPUT_UNIT,"(a)",advance="no") "["//msg//"] "
    if( cond )then
       write(OUTPUT_UNIT,"(a)") "passed"
    else
       write(OUTPUT_UNIT,"(a)") "FAILED"
    endif
  end subroutine assert
  
end program TestPostfixCalculator

Answer 1

Things to improve in the current solution:

1. It is always better to use an integer, parameter for the desired kinds of types. You can still set integer, parameter :: wp = kind(1.d0) to achieve the same result as currently, but you can change it in one place, if you want to.

2. Some reused functionality should be encapsulated into functions. (For example the string to number conversion and back.)

3. intent(out), allocatable dummy arguments are automatically deallocated. if(allocated(token)) deallocate(token) can be ommited.

4. The check for specific operators is exclusive. (If it is a "+" you don't have to check anymore if it is a "-" etc.) It should either become if - else if - else if ... or a select case statement. It informs human readers of the code, that the cases are excluding each other. Enumerated ifs should be IMHO only used if you explicitly want to fall through all possibilities and if the order of the ifs matter. As in

if (use_mpi .and. .not. mpi_initialized) call MPI_Init(ierr)
! fancy library relies on MPI
if (use_fancy_library_to_speedup_fancy_algorithm) call init_fancy_library()

Possibly increased performance is a nice addition.

5. One implicit none per program and per module scope is sufficient.¹ If your compiler supports it implicit none(type, external) is preferred. If you forget to import names of subroutines it will fail at compile instead of linking time, which speeds up the trial-and-error loop in larger projects.

Architecture:

1. Using a string stack makes the code unnecessarily complicated and deteriorates precision (and possibly performance), because you convert back and forth between floating points and their string representaiton.. Basically only a real(wp) array of MAX_ARITY size is required as Stack.

2. The GetToken routine is a "Schlemiel the Painter's Algorithm". For the n-th token you have to loop through all previous tokens and you do this n times. It would be probably better to return an array of tokens, or to keep track of the current position in the string. This requires a bit more memory than in the current solution, but this memory demand only increases linearly with the length of the expression.

3. It would be perhaps better to separate parsing from evaluating. The parsing could return a function pointer which is then evaluated on the operands.

4. At the moment you test if something is an operator. If it is not, you assume that it can be converted to a real number. If an invalid operator is passed, the error will be something like Fortran runtime error: Bad real number in item 1 of list input depending on the IO functionality of the specific runtime library of your compiler. I would rather try to convert anything to a number and anything that cannot be converted might be a valid operator. Then you can check if that operator exists or not.

5. If the operators are not operating on individual arguments, but directly on the stack a lot of special casing code can go away. The - operator takes two values and appends one (The operation might look like this: [5, 1, 3] -> [5, 2]), the PI operator just appends one value, and so on. This makes the generalization to arbitrary aryness very easy. One can e.g. implement a mean function that consumes the whole stack and appends one element. If one looks to functional languages here that is the recommended way to go there as well.

Since I had a fixed stack class at hand the actual implementation became quite easy. (I did not implement all operators, but this should be straightforward.)

module constants_mod
    implicit none(type, external)
    public
    integer, parameter :: wp = kind(1.d0)
    real(wp), parameter :: PI = 4._wp * atan(1._wp)
end module


module stack_mod
    use constants_mod, only: wp
    implicit none(type, external)
    private
    public :: Stack_t

    integer, parameter :: STACK_SIZE = 50

    type :: Stack_t
        private
        real(wp) :: values(STACK_SIZE)
        integer :: pos = 0
    contains
        private
        procedure, public :: push_back
        procedure, public :: pop
        procedure, public :: size => my_size
        procedure, public :: capacity
    end type

contains

    !> @brief
    !> Push value onto stack. Aborts if stack size is exceeded.
    subroutine push_back(this, x)
        class(Stack_t), intent(inout) :: this
        real(wp), intent(in) :: x
        if (this%pos == size(this%values)) error stop 'push back would exceed stack size.'
        this%pos = this%pos + 1
        this%values(this%pos) = x
    end subroutine


    !> @brief
    !> Pop value from stack.  Aborts if stack is empty.
    real(wp) function pop(this)
        class(Stack_t), intent(inout) :: this
        if (this%pos == 0) error stop 'It is not possible to pop from empty stack.'
        pop = this%values(this%pos)
        this%pos = this%pos - 1
    end function

    !> @brief
    !> Return current size of stack.
    integer elemental function my_size(this)
        class(Stack_t), intent(in) :: this
        my_size = this%pos
    end function

    !> @brief
    !> Return the overall capacity (i.e. upper bound for size).
    integer elemental function capacity(this)
        class(Stack_t), intent(in) :: this
        capacity = size(this%values)
    end function
end module


module reverse_polish_calculator_mod
    use, intrinsic :: iso_fortran_env
    use, intrinsic :: ieee_arithmetic
    use constants_mod, only: wp, PI
    use stack_mod, only: Stack_t
    implicit none(type, external)

    private
    public :: RPN_eval

    type :: Token_t
        character(len=:), allocatable :: str
    end type

contains

    !> @brief
    !> Return true if `str` can be converted to floating point number.
    !>
    !> @details
    !> if true, the converted number is written to `x`.
    !> if false, `x` is set to NaN.
    logical function is_number(str, x)
        character(*), intent(in) :: str
        real(wp), intent(out) :: x
        integer :: ierr
        read(str, *, iostat=ierr) x
        is_number = ierr == 0
        if (.not. is_number) then
            x = ieee_value(x, ieee_quiet_nan)
        end if
    end function


    !> @brief
    !> Split string by whitespace.
    pure function tokenize(expr) result(res)
        character(*), intent(in) :: expr
        type(Token_t), allocatable :: res(:)
        character(len=1), parameter :: delimiter = ' '
        type(Token_t), allocatable :: tmp(:)

        integer :: n, low, high

        allocate(tmp(len(expr) / 2 + 1))
        low = 1; n = 0
        do while (low <= len(expr))
            do while (expr(low : low) == delimiter)
                low = low + 1
                if (low > len(expr)) exit
            end do
            if (low > len(expr)) exit

            high = low
            if (high < len(expr)) then
                do while (expr(high + 1 : high + 1) /= delimiter)
                    high = high + 1
                    if (high == len(expr)) exit
                end do
            end if
            n = n + 1
            tmp(n)%str = expr(low : high)
            low = high + 2
        end do
        res = tmp(: n)
    end function

    !> @brief
    !> Evaluate a string expression in reverse polish notation.
    function RPN_eval(expr) result(res)
        character(*), intent(in) :: expr
        real(wp) :: res

        type(Token_t), allocatable :: tokens(:)
        type(Stack_t) :: stack
        real(wp) :: x
        real(wp) :: A, B
        integer :: i

        tokens = tokenize(expr)
        do i = 1, size(tokens)
        associate(token => tokens(i)%str)
            if (is_number(token, x)) then
                call stack%push_back(x)
            else
                select case(token)
                ! 0-ary operators
                case("PI")
                    call stack%push_back(PI)
                case("random_number")
                    call random_number(A)
                    call stack%push_back(A)
                ! 1-ary operators
                case("exp")
                    A = stack%pop()
                    call stack%push_back(exp(A))
                ! 2-ary operators
                case("+")
                    A = stack%pop()
                    B = stack%pop()
                    call stack%push_back(A + B)
                case("-")
                    A = stack%pop()
                    B = stack%pop()
                    call stack%push_back(A - B)
                case("*")
                    A = stack%pop()
                    B = stack%pop()
                    call stack%push_back(A * B)
                case("/")
                    A = stack%pop()
                    B = stack%pop()
                    call stack%push_back(A / B)
                ! any-ary operators
                case("mean")
                    block
                        integer :: N
                        real(wp) :: acc
                        N = 0; acc = 0._wp
                        do while (stack%size() > 0)
                            acc = stack%pop() + acc
                            N = N + 1
                        end do
                        call stack%push_back(acc / real(N, wp))
                    end block
                case default
                    error stop "Operator "//token//" not known"
                end select
            end if
        end associate
        end do
        ! Here you could force that the Stack has to be reduced
        ! to one number using stack%size == 1.
        res = stack%pop()
    end function
end module


program reverse_polish_calculator_prog
    use reverse_polish_calculator_mod, only: RPN_eval
    implicit none(type, external)


    write(*, *) RPN_eval("7.2 0.8 +")
    write(*, *) RPN_eval("7.2 0.8 + 2 +")
    write(*, *) RPN_eval("PI PI - PI")
    write(*, *) RPN_eval("2 4 mean")

end program

Additional Notes:

1. If the stack does not have a fixed capacity, but reallocates and grows upon push_back (like C++'s std::vector::push_back). This implementation works on arbitrary large expressions. If the tokenize function does not return an array of tokens, but becomes something like a generator i.e. returns the next token upon request, the memory demand of the tokenizing step does not grow with the expression size.
2. It is tempting to write e.g. for the - operator:

call stack%push_back(stack%pop() - stack%pop())

This is unfortunately not valid. Which I had to clarify for myself here.

---

¹ Except if you are writing interfaces. There you have to repeat the implicit none.

Answer 2

Thank you for the code alternative. I like many of the changes you did: The non-quadratic tokenizer, of course, as well as the use of a stack of the minimum size needed, and the is_number check.

I am less fond of the allocation of a temporary array to the possible maximum number of tokens, though, even if the data in the unused elements is not allocated. As I see it, the function you wrote splits naturally in a token counter and in a token reader. Once the number of tokens is known, it is easy to fetch them from the input expression with the index function. Therefore, I would rather replace the tokenize function with something like the following

    !> @brief
    !> Count tokens in string.
    pure integer function countTokens(expr,delimiter) result(nTokens)
        character(*), intent(in)    :: expr
        character(*), intent(in)    :: delimiter
        !
        integer                     :: low, high

        low = 1; nTokens = 0
        do while (low <= len(expr))
            do while (expr(low : low) == delimiter)
                low = low + 1
                if (low > len(expr)) exit
            end do
            if (low > len(expr)) exit
            high = low
            if (high < len(expr)) then
                do while (expr(high + 1 : high + 1) /= delimiter)
                    high = high + 1
                    if (high == len(expr)) exit
                end do
            end if
            nTokens = nTokens + 1
            low = high + 2
        end do
    end function countTokens

    !> @brief
    !> Split string by whitespace.
    pure function tokenize(expr) result(res)
        character(*), intent(in)    :: expr
        type(Token_t), allocatable  :: res(:)
        !
        character(len=1), parameter :: delimiter     = ' '
        integer         , parameter :: TOKEN_MAX_LEN = 50
        !
        character(len=TOKEN_MAX_LEN):: sBuf
        integer                     :: iToken, nTokens, low

        nTokens = countTokens(expr,delimiter)
        allocate(res(nTokens))
        low=1
        do iToken = 1, nTokens
            read(expr(low:),*) sBuf
            res(iToken)%str = trim(adjustl(sBuf))
            if(iToken == nTokens)exit
            low = low + index(expr(low:),res(iToken)%str) - 1
            low = low + index(expr(low:)," ")
        enddo

    end function

I am on the fence regarding the select case statement. I originally used this same algorithm, but I did not like (and still do not like) repeating the boilerplate code A = stack%pop(); B = stack%pop() over and over. The only possible advantage I can see is if the compiler implements a binary search across the listed cases (?), which would be best for a large number of operators, of course. If, however, it goes linearly through the cases, then the gain over a list of ifs is just of about a factor of two.

I understand also the appeal of operators with arbitrary arity. I think a possible approach to avoid too much boilerplate would be to separately specify the kind (fixed number, number specified at run time, or to be determined from the stack) and, if applicable, the value of this operator attribute, and on this basis popping as many operands as needed.

Btw, as a MOLCAS contributor (if I understand correctly from your Chemistry stackExchange posts) you may know Jeppe Olsen, who used to be a MOLCAS contributor too. We are building a molecular photoionization code together, now.