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

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)
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

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

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
if( cond )then
write(OUTPUT_UNIT,"(a)") "passed"
else
write(OUTPUT_UNIT,"(a)") "FAILED"
endif
end subroutine assert

end program TestPostfixCalculator


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.

1. 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
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


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.

• Thanks for looking into it. Fair points. The getToken is indeed awkward, since the program iterates over all of them anyway. I should have payed more attention to that item, rather than recycle the code. Returning an array of tokens is certainly the best option there. Jan 23 at 19:30
• Thanks for reviving the Fortran Code Review. ;-) Jan 24 at 13:13
• I have appended an answer to your code Jan 26 at 12:00

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
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.

• > 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. That's cool. I am in Ali Alavi's group in Stuttgart and at the moment developing GAS in FCIQMC (which Jeppe Olsen did in LUCIA for conventional CI calculation). The world is small ;-) Jan 26 at 12:32
• > I am on the fence regarding the select case statement. My reasoning is not because of performance but to be more self-documenting. Lots of ifs should IMHO only be used, if you explicitly want to fallthrough all possibilities and if the order of the ifs matter. Possibly increased performance is just a plus. Jan 26 at 12:38
• Small indeed. Lucia is a very nice tool, and Jeppe has done a tremendous job in extending it for our purposes. We use Lucia to compute the CI parent ions and their properties. We then cloak the ions with a hybrid GTOs+numerical functions for the photoelectron, and go from there. Soon we'll get to the multichannel TDSE, for which I plan to use coarrays (I used PETSc for the propagator in some past projects, but I have lost faith, after seeing that, in a single node, multithreaded MKL runs five or six times faster). Jan 26 at 15:00