Added files

GNUmakefile

@@ -0,0 +1,49 @@
+.PHONY: veryclean clean force export help
+
+ifdef debug
+F_OPTS := -ggdb
+endif
+
+F:= ifort
+
+SRC := rpn_stack.f90 rpn_stack_sm.f90 linked_list.f90 amap.f90
+OBJ := $(SRC:%.f90=%.o)
+
+hp: hp.f90 $(OBJ)
+	$(F) -o $@ hp.f90 $(OBJ) $(F_OPTS)
+
+hp.exe: hp.f90 GNUmakefile
+	$(F) -o $@ $< $(F_OPTS)
+	
+rpn_stack_sm.o: rpn_stack_sm.f90
+	$(F) -c -o $@ $< $(F_OPTS)
+	
+rpn_stack.o: rpn_stack.f90
+	$(F) -c -o $@ $< $(F_OPTS)
+
+linked_list.o: linked_list.f90
+	$(F) -c -o $@ $< $(F_OPTS)
+
+amap.o: amap.f90
+	$(F) -c -o $@ $< $(F_OPTS)
+
+test_amap: test_amap.f90 amap.o
+	$(F) -o $@ test_amap.f90 amap.o
+
+clean:
+	@rm -vf *.o *.mod *.smod *~
+
+veryclean: clean
+	@rm -vf hp hp.exe hp.tar
+
+force: veryclean
+	$(MAKE)
+
+export: hp.tar
+
+hp.tar: GNUmakefile hp.f90 $(SRC)
+	tar cf $@ hp.f90 $(SRC) GNUmakefile
+
+help:
+	@echo "SRC = $(SRC)"
+	@echo "OBJ = $(OBJ)"

amap.f90

@@ -0,0 +1,195 @@
+! Associative map string -> real(8)
+module amap
+    implicit none
+    
+    ! The key
+    type key_t
+        private
+        character(len=16)  :: k = '-'
+    contains
+        procedure, private :: equals_key_t
+        generic, public    :: operator(==) => equals_key_t
+        procedure, private :: write_key_t
+        generic, public    :: write(formatted) => write_key_t
+        procedure, private :: set_to_key_t
+        generic, public    :: assignment(=) => set_to_key_t
+    end type key_t
+
+    ! The value
+    type value_t
+        private
+        real(8) :: v = huge(0.0d0) ! An out-of-band value
+    contains
+        procedure, private :: write_value_t
+        generic, public    :: write(formatted) => write_value_t
+        procedure, private :: set_to_value_t
+        generic, public    :: assignment(=) => set_to_value_t
+    end type value_t
+    
+    ! Map elements are (key,value) pairs
+    type pair_t
+        type(key_t)   :: k = key_t()
+        type(value_t) :: v = value_t()
+    end type pair_t
+    
+    ! The map
+    type amap_t
+        private
+        type(pair_t), allocatable :: pairs(:)
+        integer, private :: extent = 10
+        integer, private :: high_water = 0
+    contains
+        procedure, public  :: get => get_amap_t
+        procedure, public  :: get_value => get_value_amap_t
+        procedure, public  :: set => set_amap_t
+        procedure, public  :: find => find_amap_t
+        procedure, public  :: print => print_amap_t
+        procedure, public  :: clear => clear_amap_t
+        procedure, private :: is_key_kt
+        procedure, private :: is_key_kvt
+        generic, public    :: contains => is_key_kt, is_key_kvt
+    end type amap_t
+
+contains
+
+    subroutine clear_amap_t(this)
+        class(amap_t), intent(inout)  :: this
+        if (allocated(this%pairs)) then
+            deallocate(this%pairs)
+        end if
+        this%high_water = 0
+    end subroutine clear_amap_t
+    
+    subroutine set_to_key_t(this, k)
+        class(key_t), intent(inout)  :: this
+        character(len=*), intent(in) :: k
+        this%k = adjustl(k)
+    end subroutine set_to_key_t
+    
+    subroutine set_to_value_t(this, v)
+        class(value_t), intent(inout)  :: this
+        real(8), intent(in) :: v
+        this%v = v
+    end subroutine set_to_value_t
+
+    subroutine print_amap_t(this)
+        class(amap_t), intent(in) :: this
+        integer :: i
+        write(6,'(a,i0,a)') 'Map has ',this%high_water,' elements'
+        do i=1,this%high_water
+            write(6,'(4x,dt,a,dt)') this%pairs(i)%k,' -> ',this%pairs(i)%v
+        end do
+    end subroutine print_amap_t
+    
+    subroutine set_amap_t(this,kv,vv)
+        class(amap_t), intent(inout) :: this
+        character(len=*), intent(in) :: kv
+        real(8), intent(in)          :: vv
+        type(pair_t), allocatable :: tmp_pairs(:)
+        type(key_t)               :: k
+        type(value_t)             :: v
+        integer                   :: idx
+        
+        k = kv
+        v = vv
+        if (.not. allocated(this%pairs)) then
+            allocate(this%pairs(this%extent))
+        end if    
+        
+        idx = this%find(k)
+        if (idx > 0) then
+            this%pairs(idx) = pair_t(k,v)
+            return
+        end if
+        
+        if (this%high_water == size(this%pairs)) then
+            allocate(tmp_pairs(size(this%pairs)+this%extent))
+            tmp_pairs(1:this%high_water) = this%pairs
+            call move_alloc(tmp_pairs, this%pairs)
+        end if
+        this%high_water = this%high_water + 1
+        this%pairs(this%high_water) = pair_t(k,v)
+        
+    end subroutine set_amap_t
+    
+    function find_amap_t(this, k) result(r)
+        class(amap_t), intent(in) :: this
+        type(key_t), intent(in) :: k
+        integer :: r
+        do r = 1, this%high_water
+            if (this%pairs(r)%k == k) then
+                return
+            end if
+        end do
+        r = 0
+    end function find_amap_t
+    
+    function get_amap_t(this, kv) result(r)
+        class(amap_t), intent(in)    :: this
+        character(len=*), intent(in) :: kv
+        type(value_t) :: r
+        type(key_t) :: k
+        integer :: idx
+        k = kv
+        idx = this%find(k)
+        if (idx > 0) then
+            r = this%pairs(idx)%v
+        else
+            r = value_t()
+        end if
+    end function get_amap_t
+    
+    function get_value_amap_t(this, kv) result(r)
+        class(amap_t), intent(in)    :: this
+        character(len=*), intent(in) :: kv
+        real(8) :: r
+        type(value_t) :: s
+        s = this%get(kv)
+        r = s%v
+    end function get_value_amap_t
+    
+    function is_key_kt(this, k) result(r)
+        class(amap_t), intent(in) :: this
+        type(key_t), intent(in) :: k
+        logical :: r
+        r = this%find(k) > 0
+    end function is_key_kt
+    
+    function is_key_kvt(this, kv) result(r)
+        class(amap_t), intent(in) :: this
+        character(len=*), intent(in) :: kv
+        logical :: r
+        r = this%find(key_t(kv)) > 0
+    end function is_key_kvt
+    
+    function equals_key_t(this, k) result(r)
+        class(key_t), intent(in) :: this
+        type(key_t), intent(in) :: k
+        logical :: r
+        r = trim(adjustl(this%k)) == trim(adjustl(k%k))
+    end function equals_key_t
+    
+    subroutine write_key_t(key, unit, iotype, v_list, iostat, iomsg)
+        class(key_t), intent(in)    :: key
+        integer, intent(in)         :: unit
+        character(*), intent(in)    :: iotype
+        integer, intent(in)         :: v_list(:)
+        integer, intent(out)        :: iostat
+        character(*), intent(inout) :: iomsg
+        iostat = 0
+        iomsg = ""
+        write(6,'(a)', iostat=iostat, iomsg=iomsg) trim(adjustl(key%k))
+    end subroutine write_key_t
+
+    subroutine write_value_t(value, unit, iotype, v_list, iostat, iomsg)
+        class(value_t), intent(in)    :: value
+        integer, intent(in)         :: unit
+        character(*), intent(in)    :: iotype
+        integer, intent(in)         :: v_list(:)
+        integer, intent(out)        :: iostat
+        character(*), intent(inout) :: iomsg
+        iostat = 0
+        iomsg = ""
+        write(6,'(f0.6)', iostat=iostat, iomsg=iomsg) value%v
+    end subroutine write_value_t
+end module amap

hp.f90

@@ -0,0 +1,920 @@
+program hp15c
+    use rpn_stack
+    use linked_list, print_ll => print, clear_ll => clear, size_ll => size
+    use amap
+    
+    implicit none
+    
+    real(8)                   :: x
+    integer                   :: ios, i
+    integer                   :: verbosity = 0
+    character(100)            :: buff
+    integer                   :: blen
+    integer                   :: argl, argc
+    type(llist)               :: tokens
+    type(llist_node), pointer :: token
+
+    real(8), parameter :: ag = 9.80665d0
+    real(8), parameter :: g = 6.67430d-11
+    real(8), parameter :: e = exp(1.0d0)
+    real(8), parameter :: c = 2.99792458d8
+    type(amap_t)       :: constants
+    
+    type(amap_t)       :: stats
+    integer            :: in_sequence = 0
+    logical            :: seq_is_x
+    real(8), allocatable  :: x_seq(:), y_seq(:)
+    integer            :: n_seq = 0
+    
+    logical      :: veMode = .false.
+    logical      :: lang_en = .true.
+    logical      :: tmp_cmode
+    logical      :: ok
+    logical      :: getNext, numbers, have_expression
+    integer      :: stat
+    character(len=100) :: msg
+    character(5) :: lang
+
+    type(rpn_t)  :: mem(0:9) = rpn_t()
+    
+    ! Create a stack of size 4
+    type(stack_t(5)) :: stack
+
+    call stack%set_legend(['x:','y:','z:','s','t:'])
+    degrees_mode = .true.
+    complex_mode = .false.
+    eps = 1.0d-14
+    
+    ! Constants
+    call constants%set('g',ag)
+    call constants%set('G',g)
+    call constants%set('e',e)
+    call constants%set('c',c)
+    call constants%set('pi',pi)
+    call constants%set('two_pi',2*pi)
+    call constants%set('pi_over_2',pi/2)
+
+    ! Try to read the LANG environment variable
+    call get_environment_variable('LANG',lang,status=stat)
+    lang_en = stat /= 0
+    if (.not. lang_en) then
+        lang_en = merge(.true.,.false.,lang(1:3) == 'en_')
+    end if
+    lang = merge('POINT','COMMA',lang_en)
+
+    call init(lang)
+
+    ! Interrogate argument list
+    argc = command_argument_count()
+    have_expression = .false.
+    do i=1,argc
+        call get_command_argument(i, buff, argl)
+        if (buff(1:argl) == '-d') then
+            verbosity = 1
+            cycle
+            
+        else if (buff(1:argl) == '-c') then
+            complex_mode = .true.
+            cycle
+            
+        else if (buff(1:argl) == '-v') then
+            veMode = .true.
+            cycle
+        
+        else if (buff(1:argl) == '-h') then
+            call help
+            stop
+       
+        end if
+
+        have_expression = .true.
+        
+        ! Break the string up into a linked-list of tokens
+        call tokenize(buff(1:argl))
+        if (verbosity > 0) call print_ll(tokens)
+        
+        ! Interpret each token as a command and appky it
+        ok = tokens%iterate(apply_command)
+        
+        ! Do not print the stack at the end of a sequence -it's confusing
+        if (in_sequence /= 2) then
+            call stack%print(veMode)
+        end if
+        
+        ! Tidy
+        call clear_ll(tokens)
+
+        if (.not. ok) stop
+        
+    end do
+    
+    if (.not. have_expression) then
+        call stack%print(veMode)
+    end if
+    
+    ! Loop until quit
+    all :do
+        x = 0.0d0
+        buff = ''
+        write(6,'(a)',advance='no') ':: '
+        read(5,fmt='(a)',iostat=ios,iomsg=msg) buff
+        if (ios /= 0) then
+            write(6,'(/a)') 'Command:['//buff(1:blen)//']'//'; '//msg
+            cycle all
+        end if
+        buff = trim(adjustl(buff))
+        blen = len_trim(buff)
+        if (blen == 0) cycle all
+
+        ! Tokenize input string
+        call tokenize(buff(1:blen))
+        ok = tokens%iterate(apply_command)
+        if (.not. ok) exit all
+        
+        if (in_sequence == 1) then
+            write(6,'(i0)') n_seq
+        else if (in_sequence == 2) then
+            in_sequence = 0
+        else
+            call stack%print(veMode)
+        end if
+    end do all
+
+    call clear_ll(tokens)
+    stop
+
+contains
+       
+  subroutine tokenize(com)
+    character(*), intent(in)      :: com
+    integer                       :: start, end
+    character(len=:), allocatable :: command
+    
+    call clear_ll(tokens)
+    if (len_trim(com) == 0) then
+        return
+    end if
+    start = 1
+    ! Ensure there are no leading and trailing spaces
+    command = trim(adjustl(com))
+    end = index(command,' ')
+    end = merge(len(command),end-1,end==0)
+    do
+        call append(tokens,command(start:end))
+        if (end == len(command)) exit
+        start = end + nsp(command(end+1:))
+        end = index(command(start:),' ') - 1
+        end = merge(len(command),end+start-1,end == -1)
+     end do
+  end subroutine tokenize
+
+  subroutine help
+    write(6,'(/a)') 'Command Calculator'
+    write(6,'(a/)') '=================='
+    write(6,'(a)')  'Introduction'
+    write(6,'(a/)') '------------'
+    write(6,'(a)')  'This is a command-line calculator. It supports both real and complex modes, as well'
+    write(6,'(a)')  'as degrees/radians selection and precision control. It can be run interactively or as an'
+    write(6,'(a/)') 'expression parser. This help is deliberately terse to encourage exploration.'
+    write(6,'(a/)') '-------------------------------------------------------------------------------'
+    write(6,'(a)')  'Operators: + - * / ^ ^/x ^x ^2 ^/2 ^3 ^/3 ^*2 ^*10 || ! %'
+    write(6,'(a)')  'Constants: pi e g G c two_pi pi_over_2'
+    write(6,'(a)')  'Functions: sin cos tan asin acos atan sinh cosh tanh log2 log lg len sq sqrt cb cbrt'
+    write(6,'(a)')  '           alog2 alog alog10 gamma ncr npr rem int nint'
+    write(6,'(a)')  ' Controls: fix[0-9] clx cl cla '
+    write(6,'(a)')  '    Modes: real complex verbose terse degrees radians'
+    write(6,'(a)')  ' Memories: n=0...9  st<n> sw<n> rc<n> cl<n> m<n>+ m<n>- m<n>* m<n>/ msh'
+    write(6,'(a)')  '  Complex: ri _ || to_pol to_cart'
+    write(6,'(a)')  '  Actions: 1/ -- R r ? > < split drop'
+    write(6,'(a)')  '    Stats: { x1 x2 ... } { x1,y1 x2,y2 ... }'
+    write(6,'(a)')  '           n ux sx mx lqx uqx uy sy my lqy uqy a b cov corr'
+    write(6,'(a/)') '    Quits: q'
+    write(6,'(a/)') '-------------------------------------------------------------------------------'
+    write(6,'(a)')  'Examples'
+    write(6,'(a)')  '--------'
+    write(6,'(4x,a)')  'hp "fix2 18 2 - 8 2 / * ="                    -> 64.00'
+    write(6,'(4x,a)')  'hp "2 -- complex sqrt ="                      -> (0.00000,-1.414214)'
+    write(6,'(4x,a/)') 'hp -c "radians (1,pi_over_2)p ^ * degrees ="  -> (1.000000,180.000000) p'
+
+  end subroutine help
+
+  subroutine apply_command(command, ok)
+    use, intrinsic :: ieee_arithmetic
+    implicit none
+    character(*), intent(in) :: command
+    logical, intent(out)     :: ok
+
+    real(8)                 :: r, im, ang
+    complex(8)              :: u, z
+    real(8), allocatable    :: tmp_seq(:)
+    type(rpn_t)             :: us, zs
+    logical                 :: is_cart
+    integer                 :: m, idx
+    character(len=1)        :: comma
+    character(5), parameter :: lang(2) = ['POINT','COMMA']
+    
+    ok = .true.
+    if (len_trim(command) == 0) then
+        return
+    end if
+    
+    if (verbosity > 0) then
+        write(*,'(a)') 'Applying: '//command
+    end if
+    
+    if (in_sequence == 1) then
+        if (command == '}') then
+            in_sequence = 2
+            complex_mode = tmp_cmode
+            call calculate_stats
+            
+        else
+            ! All elements must be the same so either all x or all x,y
+            idx = index(command,',')
+            if (n_seq == 0) then
+                seq_is_x = (idx == 0)
+            end if
+            if (seq_is_x .neqv. (idx == 0)) then
+                goto 901
+            end if
+            if (seq_is_x) then
+                read(command,*,err=901) r
+                im = 0
+            else
+                read(command(1:idx-1),*,err=901) r
+                read(command(idx+1:len(command)),*,err=901) im
+            end if
+            ! Initial allocation
+            if (n_seq == 0 .and. .not. allocated(x_seq)) then
+                allocate(x_seq(10))
+                if (.not. seq_is_x) then
+                    allocate(y_seq(10))
+                end if
+            end if
+            
+            if (n_seq < size(x_seq)) then
+                n_seq = n_seq + 1
+            else
+                ! Expand array
+                allocate(tmp_seq(n_seq + 10))
+                tmp_seq(1:n_seq) = x_seq
+                call move_alloc(tmp_seq, x_seq)
+                if (.not. seq_is_x) then
+                    allocate(tmp_seq(n_seq + 10))
+                    tmp_seq(1:n_seq) = y_seq
+                    call move_alloc(tmp_seq, y_seq)
+               end if
+            end if
+            x_seq(n_seq) = r
+            if (.not. seq_is_x) then
+                y_seq(n_seq) = im
+            end if
+            if (verbosity > 0) then
+                print *,x_seq(1:n_seq)
+            end if
+        end if
+        return
+    end if
+    
+    select case(command)
+       
+    case('q','quit')
+       ok = .false.
+       return
+       
+    case('=')
+       ok = .false.
+       return
+       
+    case('{')
+        ! Start sequence
+        in_sequence = 1
+        n_seq = 0
+        tmp_cmode = complex_mode
+        complex_mode = .false.
+        call stats%clear()
+        
+    case('--')
+        call invoke_unary(chs_fr)
+       
+    case('^')
+        call stack%push(stack%peek(1))
+       
+    case('+')
+        call invoke_binary(add_fr)
+       
+    case('-')
+        call invoke_binary(subtract_fr)
+       
+    case('*')
+        call invoke_binary(multiply_fr)
+       
+    case('/')
+        call invoke_binary(divide_fr)
+       
+    case('^x')
+        call invoke_binary(power_fr)
+       
+    case('^/x')
+        ! Only raising to a real power is supported
+        zs = stack%peek(1)
+        if (zs%is_real()) then
+            call invoke_binary(root_fr)
+        else
+            goto 901
+        end if
+        
+    case('>')
+        call invoke_unary(next_root_fr)
+       
+    case('<')
+        call invoke_unary(previous_root_fr)
+        
+    case('%')
+        call invoke_binary(percent_fr)
+        
+    case('xy','XY')
+       call stack%swap
+       
+    case('R')
+       call stack%rotate_down
+       
+    case('r')
+       call stack%rotate_up
+      
+    case('CLA','cla')
+        call stack%clear
+        mem = rpn_t()
+        call stats%clear
+        
+    case('CL','cl')
+        call stack%clear
+
+        
+    case('CLX','clx')
+       call stack%set(rpn_t())
+       
+    case('_')
+        if (complex_mode) then
+            call invoke_unary(conj_fr)
+        endif
+
+    case('len','||')
+        if (complex_mode) then
+            call invoke_unary(len_fr)
+            ! Length is always reported as (x,0) and marked is_cartesian
+            zs = stack%peek(1)
+            call zs%set_value(is_cartesian=.true.)
+            call stack%set(zs,1)
+        else
+            zs = stack%peek(1)
+            z = zs%get_value()
+            us = stack%peek(2)
+            u = us%get_value()
+            call zs%set_value(cmplx(hypot(z%re,u%re),0,8))
+            call stack%set(zs)
+        end if
+        
+    case('split')
+        if (.not. complex_mode) then
+            zs = stack%pop()
+            x = zs%get_value()
+            if (x > 0) then
+                r = floor(x)
+            else
+                r = ceiling(x)
+            end if
+            im = x - r
+            call stack%push(im)
+            call stack%push(r)
+        end if
+         
+    case('int')
+        if (.not. complex_mode) then
+            zs = stack%peek(1)
+            x = zs%get_value()
+            if (x > 0) then
+                r = floor(x)
+            else
+                r = ceiling(x)
+            end if
+            call zs%set_value(cmplx(r,0,8))
+            call stack%set(zs,1)
+        end if
+         
+    case('nint')
+        if (.not. complex_mode) then
+            zs = stack%peek(1)
+            x = zs%get_value()
+            r = nint(x)
+            if (mod(r,2.0d0) == 1) then
+                r = r - 1
+            end if
+            call zs%set_value(cmplx(r,0,8))
+            call stack%set(zs,1)
+        end if
+         
+    case('rem')
+        if (.not. complex_mode) then
+            zs = stack%peek(1)
+            x = zs%get_value()
+            if (x > 0) then
+                r = floor(x)
+            else
+                r = ceiling(x)
+            end if
+            im = x - r
+            call zs%set_value(cmplx(im,0,8))
+            call stack%set(zs,1)
+        end if
+        
+    case('drop')
+       zs = stack%pop()
+       
+    case('ri')
+       ! Swap real and imaginary parts
+        if (complex_mode) then
+            call invoke_unary(swap_real_imaginary_fr)
+        end if
+
+    case('to_pol')
+        ! Convert x + iy to r + i theta
+        if (complex_mode) then
+            zs = stack%peek(1)
+            call stack%set(to_polar(zs))
+        end if
+
+    case('to_cart')
+        ! Convert (r,theta) to (x,y)
+        if (complex_mode) then
+            zs = stack%peek(1)
+            call stack%set(to_cartesian(zs))
+        end if
+       
+    case('1/')
+        call invoke_unary(reciprocal_fr)
+       
+    case('^2','sq')
+        call invoke_unary(power_2_fr)
+       
+    case('^/2','sqrt')
+        call invoke_unary(sqrt_fr)
+        
+    case('^3','cb')
+        call invoke_unary(power_3_fr)
+      
+    case('^/3','cbrt')
+        call invoke_unary(cbrt_fr)
+       
+    case('^*2','alog2')
+        call invoke_unary(exp_2_fr)
+    
+    case('^*10','alog10')
+        call invoke_unary(exp_10_fr)
+
+    case('exp','alog')
+        call invoke_unary(exp_e_fr)
+
+    case('ln')
+        call invoke_unary(ln_fr)
+
+    case('log2')
+        call invoke_unary(log2_fr)
+
+    case('lg')
+        call invoke_unary(lg_fr)
+
+    case('sinh')
+        call invoke_unary(hsine_fr)
+
+    case('cosh')
+        call invoke_unary(hcosine_fr)
+
+    case('tanh')
+        call invoke_unary(htangent_fr)
+
+    case('sin')
+        call invoke_unary(sine_fr)
+
+    case('cos')
+        call invoke_unary(cosine_fr)
+
+    case('tan')
+        call invoke_unary(tangent_fr)
+
+    case('asin')
+        call invoke_unary(asine_fr)
+
+    case('asinh')
+        call invoke_unary(ahsine_fr)
+
+    case('acos')
+        call invoke_unary(acosine_fr)
+
+    case('acosh')
+        call invoke_unary(ahcosine_fr)
+
+    case('atan')
+        call invoke_unary(atangent_fr)
+
+    case('atanh')
+        call invoke_unary(ahtangent_fr)
+
+    case('atan2')
+        call invoke_binary(atangent2_fr)
+
+    case('gamma')
+        call invoke_unary(gamma_fr)
+
+    case('!')
+        zs = stack%peek(1)
+        if (zs%is_positive_real()) then
+            call invoke_unary(fact_fr)
+        else
+            goto 901
+        end if
+
+    case('ncr')
+        zs = stack%peek(1)
+        us = stack%peek(2)
+        if (zs%is_positive_real() .and. us%is_positive_real()) then
+            call invoke_binary(ncr_fr)
+        else
+            goto 901
+        end if
+
+    case('npr')
+        zs = stack%peek(1)
+        us = stack%peek(2)
+        if (zs%is_positive_real() .and. us%is_positive_real()) then
+            call invoke_binary(npr_fr)
+        else
+            goto 901
+        end if
+       
+    case('m0+','m1+','m2+','m3+','m4+','m5+','m6+','m7+','m8+','m9+')
+       read(command(2:2),'(i1)',err=901) m
+       mem(m) = mem(m) + stack%peek(1)
+
+    case('m0-','m1-','m2-','m3-','m4-','m5-','m6-','m7-','m8-','m9-')
+       read(command(2:2),'(i1)',err=901) m
+       mem(m) = mem(m) - stack%peek(1)
+
+    case('m0*','m1*','m2*','m3*','m4*','m5*','m6*','m7*','m8*','m9*')
+       read(command(2:2),'(i1)',err=901) m
+       mem(m) = mem(m) * stack%peek(1)
+
+    case('m0/','m1/','m2/','m3/','m4/','m5/','m6/','m7/','m8/','m9/')
+       read(command(2:2),'(i1)',err=901) m
+       mem(m) = mem(m) / stack%peek(1)
+
+    case('st0','st1','st2','st3','st4','st5','st6','st7','st8','st9')
+       read(command(3:3),'(i1)',err=901) m
+       mem(m) = stack%peek(1)
+
+    case('sw0','sw1','sw2','sw3','sw4','sw5','sw6','sw7','sw8','sw9')
+       read(command(3:3),'(i1)',err=901) m
+       zs = stack%peek(1)
+       call stack%set(mem(m))
+       mem(m) = zs
+
+    case('rc0','rc1','rc2','rc3','rc4','rc5','rc6','rc7','rc8','rc9')
+       read(command(3:3),'(i1)',err=901) m
+       call stack%push(mem(m))
+
+    case('cl0','cl1','cl2','cl3','cl4','cl5','cl6','cl7','cl8','cl9')
+       read(command(3:3),'(i1)',err=901) m
+       mem(m) = rpn_t()
+
+    case('msh')
+        write(6,'(i3,a,dt)') (i,': ',mem(i),i=0,size(mem)-1)
+            
+    case('fix0','fix1','fix2','fix3','fix4','fix5','fix6','fix7','fix8','fix9')
+       read(command(4:4),'(i1)') m
+       call set_places(m)
+    
+    case('DEG','deg','DEGREES','degrees')
+        call toggle_degrees_mode(.true.)
+
+    case('RAD','rad','RADIANS','radians')
+        call toggle_degrees_mode(.false.)
+
+    case('mC','COMPLEX','complex')
+       complex_mode = .true.
+
+    case('mR','REAL','real')
+       complex_mode = .false.
+
+    case('mV','VERBOSE','verbose')
+       veMode = .true.
+
+    case('mT','TERSE','terse')
+       veMode = .false.
+
+    case('?')
+       write(6,advance='no',fmt='(a)') 'Status: '
+       write(6,advance='no',fmt='(2a)') merge('degrees','radians',degrees_mode),' ; '
+       write(6,advance='no',fmt='(a,i0)') 'dp = ',get_places()
+       if (complex_mode) then
+          write(6,advance='no',fmt='(a)') ' ; mode = complex'
+       else
+          write(6,advance='no',fmt='(a)') ' ; mode = real'
+       end if
+       write(6,advance='no',fmt='(a,i0)') ' ; stack size = ',stack%get_size()
+       write(6,'(a/)') ''
+
+    case('help','h')
+       call help
+
+    case default
+        ! Process constants first
+        block
+            integer :: lc, is_integer,split_idx,end_idx
+            character(len=:), allocatable :: re_comp, im_comp
+            lc = len_trim(command)
+            is_integer = (index(command,'.') == 0)
+            if (complex_mode) then
+                if (command(1:1) == '(') then
+                    split_idx = index(command,',')
+                    end_idx = index(command,')')
+                    re_comp = command(2:split_idx-1)
+                    im_comp = command(split_idx+1:end_idx-1)
+                    if (constants%contains(re_comp)) then
+                        z%re = constants%get_value(re_comp)
+                    else
+                        read(re_comp,*,err=901,end=901) z%re
+                    end if
+                    if (constants%contains(im_comp)) then
+                        z%im = constants%get_value(im_comp)
+                    else
+                        read(im_comp,*,err=901,end=901) z%im
+                    end if
+                    if (command(lc:lc) == 'p') then
+                        call stack%push(z,.false.)
+                    else
+                        call stack%push(z)
+                    end if
+                else
+                    if (constants%contains(command)) then
+                        x = constants%get_value(command)
+                    else
+                        read(command,*,err=901,end=901) x
+                    end if
+                    call stack%push(cmplx(x,0.0d0,8))
+                end if
+                
+            else
+                if (constants%contains(command)) then
+                    x = constants%get_value(command)
+                else if (stats%contains(command)) then
+                    x = stats%get_value(command)
+                else
+                    read(command,*,err=901,end=901) x
+                end if
+                call stack%push(cmplx(x,0.0d0,8))
+            end if
+        end block
+    end select
+    return
+
+901 continue
+    write(6,'(a)') command//' ???'
+    return
+            
+    end subroutine apply_command
+
+    subroutine calculate_stats
+        real(8) :: a, b, c, sxy
+        real(8) :: s(5,2)
+
+        call summary_stats(x_seq(1:n_seq),s(1,1),s(2,1),s(3,1),s(4,1),s(5,1))
+        call stats%set('n',real(n_seq,8))
+        call stats%set('ux',s(1,1))
+        call stats%set('mx',s(2,1))
+        call stats%set('sx',s(3,1))
+        call stats%set('lqx',s(4,1))
+        call stats%set('uqx',s(5,1))
+        if (seq_is_x) then
+                 write(6,10) '    count n -> ',n_seq
+            call print_value('    mean ux -> ',s(1,1))
+            call print_value('  stddev sx -> ',s(3,1))
+            call print_value('  median mx -> ',s(2,1))
+            call print_value('lower_q lqx -> ',s(4,1))
+            call print_value('upper_q uqx -> ',s(5,1))
+       else
+            call summary_stats(y_seq(1:n_seq),s(1,2),s(2,2),s(3,2),s(4,2),s(5,2))
+                 write(6,10) '    count n        -> ',n_seq
+            call print_value('    means ux , uy  -> ',s(1,1),s(1,2))
+            call print_value('  stddevs sx , xy  -> ',s(3,1),s(3,2))
+            call print_value('  medians mx , my  -> ',s(2,1),s(2,2))
+            call print_value('lower_qs lqx , lqy -> ',s(4,1),s(4,2))
+            call print_value('upper_qs uqx , uqy -> ',s(5,1),s(5,2))
+            call stats%set('uy',s(1,2))
+            call stats%set('my',s(2,2))
+            call stats%set('sy',s(3,2))
+            call stats%set('lqy',s(4,2))
+            call stats%set('uqy',s(5,2))
+
+            call calculate_regression(s(1,1),s(1,2),a,b,c,sxy)
+            call stats%set('a',a)
+            call stats%set('b',b)
+            call stats%set('corr',c)
+            call stats%set('cov',sxy)
+            write(6,'(/a)') 'Regression:  y = ax + b'
+            call print_value('   gradient a      ->',a)
+            call print_value('  intercept b      -> ',b)
+            call print_value(' covariance cov    -> ',sxy)
+            call print_value('correlation corr   -> ',c)
+        end if
+        10 format(a,i0)
+        
+    end subroutine calculate_stats
+    
+    subroutine calculate_regression(mean_x, mean_y, a, b, c, sxy)
+        real(8), intent(in) :: mean_x, mean_y
+        real(8), intent(out) :: a, b, c, sxy
+        integer :: i
+        real(8) :: sxx, syy
+        sxy = sum(x_seq(1:n_seq)*y_seq(1:n_seq))/n_seq - mean_x*mean_y
+        sxx = sum(x_seq(1:n_seq)*x_seq(1:n_seq))/n_seq - mean_x**2
+        syy = sum(y_seq(1:n_seq)*y_seq(1:n_seq))/n_seq - mean_y**2
+        a = sxy/sxx
+        b = mean_y - a*mean_x
+        c = sxy/sqrt(sxx*syy)
+    end subroutine calculate_regression
+    
+    subroutine print_value(name, x, y)
+        character(len=*), intent(in)  :: name
+        real(8), intent(in)           :: x
+        real(8), intent(in), optional :: y
+        character(len=:), allocatable :: fmt_x, fmt_y
+        call to_string(x, fmt_x)
+        if (present(y)) then
+            call to_string(y, fmt_y)
+            fmt_x = fmt_x//' , '//fmt_y
+        end if
+        write(6,'(a)') name//fmt_x
+    end subroutine print_value
+
+    subroutine summary_stats(a, mean, median, stddev, lower_q, upper_q)
+        real(8), intent(in)  :: a(:)
+        real(8), intent(out) :: mean, median, stddev, lower_q, upper_q
+        real(8) :: b(size(a))
+        real(8) :: s, s2
+        integer :: m, n
+        n = size(a)
+        b = a
+        s = sum(b)
+        s2 = sum(b**2)
+        mean = s/n
+        stddev = sqrt(s2/n - (s/n)**2)
+        call sort(b)
+        median = calc_median(b, m)
+        if (n < 3) then
+            lower_q = median
+            upper_q = median
+        else
+            lower_q = calc_median(b(1:m))
+            if (mod(n,2) == 0) then
+                upper_q = calc_median(b(m+1:n))
+            else
+                upper_q = calc_median(b(m:n))
+            end if
+        end if
+    end subroutine summary_stats
+    
+    function calc_median(a, mid) result(r)
+        real(8), intent(in)  :: a(:)
+        integer, intent(out), optional :: mid
+        real(8) :: r
+        integer :: m, n
+        n = size(a)
+        m = n/2
+        if (mod(n,2) == 0) then
+            r = (a(m) + a(m+1))/2.0d0
+        else
+            m = m + 1
+            r = a(m)
+        end if
+        if (present(mid)) then
+            mid = m
+        end if
+    end function calc_median
+        
+    ! 'a' won't be very big so a simple n**2 algorithm will do
+    subroutine sort(a)
+        real(8), intent(inout) :: a(:)
+        real(8) :: b(size(a))
+        integer :: i, j(size(a))
+        logical :: mask(size(a))
+        mask = .true.
+        b = a
+        do i=1,size(a)
+            j = minloc(b, mask)
+            associate (j1 => j(1))
+                a(i) = b(j1)
+                mask(j1) = .false.
+            end associate
+        end do
+    end subroutine
+    
+    subroutine toggle_degrees_mode(new_mode)
+        logical, intent(in) :: new_mode
+        integer    :: i
+        type(rpn_t) :: rz
+        
+        ! Only do something if the modes are different
+        if (new_mode .eqv. degrees_mode) return
+        
+        degrees_mode = .not. degrees_mode
+        
+        ! Convert all polar complex numbers
+        ! 1) In the stack
+        do i=1,stack%ssize
+            rz = stack%peek(i)
+            call update_angle_unit(rz)
+            call stack%set(rz,i)
+        end do
+        
+        ! 2) in memory
+        do i=lbound(mem,1),ubound(mem,1)
+            call update_angle_unit(mem(i))
+        end do
+        
+        ! 3) in multiple roots
+        do i=1,nroots
+            call update_angle_unit(roots(i))
+        end do
+          
+    end subroutine toggle_degrees_mode
+    
+    subroutine update_angle_unit(rz)
+        type(rpn_t), intent(inout) :: rz
+        complex(8) :: zs
+        logical    :: is_cart
+        zs = rz%get_value(is_cart)
+        if (is_cart) return
+        zs%im = zs%im*merge(to_deg, to_rad, degrees_mode)
+        call rz%set_value(zs,is_cart)
+    end subroutine update_angle_unit
+    
+    subroutine invoke_binary(action)
+        procedure(binary_f), pointer, intent(in) :: action
+        type(rpn_t) :: us, zs
+        logical :: is_cart
+
+        zs = stack%pop()
+        if (complex_mode) then
+            is_cart = zs%is_cartesian()
+            if (.not. is_cart) then
+                zs = to_cartesian(zs)
+            end if
+            us = stack%peek(1)
+            if (.not. us%is_cartesian()) then
+                us = to_polar(us)
+            end if
+            us = action(us,zs)
+            if (.not. is_cart) then
+                us = to_polar(us)
+            end if
+            call stack%set(us)
+        else
+            us = stack%peek(1)
+            call stack%set(action(us,zs))
+        end if
+    end subroutine invoke_binary
+        
+    subroutine invoke_unary(action)
+        procedure(unary_f), pointer, intent(in) :: action
+        logical :: is_cart
+        type(rpn_t) :: z
+        if (complex_mode) then
+            z = stack%peek(1)
+            is_cart = z%is_cartesian()
+            if (.not. is_cart) then
+                z = to_cartesian(z)
+            end if
+            z = action(z)
+            if (.not. is_cart) then
+                z = to_polar(z)
+            end if
+            call stack%set(z)
+        else
+            call stack%set(action(stack%peek(1)))
+        end if
+    end subroutine invoke_unary
+
+   integer function nsp(command)
+     implicit none
+     character(*), intent(in) :: command
+     integer :: i
+     do i=1,len(command)
+        if (command(i:i) /= ' ') then
+           nsp = i
+           return
+        end if
+     end do
+     nsp = 0
+   end function nsp
+
+end program hp15c

linked_list.f90

@@ -0,0 +1,97 @@
+module linked_list
+    public
+    
+    type llist
+        private
+        type(llist_node), pointer :: begin => null()
+        type(llist_node), pointer :: end => null()
+    contains
+        procedure                 :: iterate => iterate_ll
+    end type llist
+  
+    type llist_node
+        private
+        character(len=:), allocatable :: data
+        type(llist_node), pointer :: next => null()
+    end type llist_node
+
+    ! Interface for functions being applied to each list element in turn
+    ! when iterating
+    abstract interface
+        subroutine command_fun(command, ok)
+            character(*), intent(in) :: command
+            logical, intent(out)     :: ok ! Exit the loop if not true
+        end subroutine command_fun
+    end interface
+
+    
+contains
+
+    function iterate_ll(this, f) result(r)
+        class(llist), intent(inout), target :: this
+        procedure(command_fun)              :: f
+        logical :: r
+        type(llist_node), pointer :: token
+        token => this%begin
+        do
+            if (.not. associated(token)) exit
+            call f(trim(token%data), r)
+            if (.not. r) then
+                exit
+            end if
+            token => token%next
+        end do
+    end function iterate_ll
+       
+    subroutine append(lst, data)
+      type(llist), intent(inout) :: lst
+      character(*), intent(in)    :: data
+      if (.not. associated(lst%begin)) then
+         allocate(lst%begin)
+         lst%begin%data = data
+         lst%end => lst%begin
+      else
+         allocate(lst%end%next)
+         lst%end%next%data = data
+         lst%end => lst%end%next
+      end if
+    end subroutine append
+
+    subroutine print(lst)
+      type(llist), intent(in) :: lst
+      type(llist_node), pointer :: next
+      write(*,'(a)') 'Tokens:'
+      next => lst%begin
+      do
+         if (.not. associated(next)) exit
+         write(*,'(4x,a)') next%data
+         next => next%next
+      end do
+    end subroutine print
+
+    integer function size(lst)
+      type(llist), intent(inout) :: lst
+      type(llist_node), pointer  :: this
+      size = 0
+      this => lst%begin
+      do
+         if (.not. associated(this)) exit
+         size = size + 1
+         this => this%next
+      end do
+    end function size
+
+    subroutine clear(lst)
+      type(llist), intent(inout) :: lst
+      type(llist_node), pointer  :: this, next
+      this => lst%begin
+      do
+         if (.not. associated(this)) exit
+         next => this%next
+         deallocate(this)
+         this => next
+      end do
+      nullify(lst%end)
+      nullify(lst%begin)
+    end subroutine clear
+end module linked_list

rpn_stack.f90

@@ -0,0 +1,467 @@
+module rpn_stack
+    implicit none
+    
+    ! Type for the data that's going on to the stack
+    type rpn_t
+        private
+        complex(8), private   :: zdata   = 0
+        logical, private      :: is_cart = .true.
+    contains
+        procedure, private :: write_rpns
+        generic, public    :: write(formatted) => write_rpns
+        procedure :: get_value                 => get_value_rpns
+        procedure :: set_value                 => set_value_rpns
+        procedure :: is_integer                => is_integer_rpns
+        procedure :: is_real                   => is_real_rpns
+        procedure :: is_positive_real          => is_positive_real_rpns
+        procedure :: is_cartesian              => is_cartesian_rpns
+        procedure :: set_angle_unit            => set_angle_unit_rpns
+        procedure, private :: add_rpns
+        generic, public    :: operator(+) => add_rpns
+        procedure, private :: subtract_rpns
+        generic, public    :: operator(-) => subtract_rpns
+        procedure, private :: multiply_rpns
+        generic, public    :: operator(*) => multiply_rpns
+        procedure, private :: divide_rpns
+        generic, public    :: operator(/) => divide_rpns
+        procedure, private :: power_rpns
+        generic, public    :: operator(**) => power_rpns
+        procedure, private :: set_to_rpns
+        generic, public    :: assignment(=) => set_to_rpns
+    end type rpn_t
+    
+    ! Make the stack a parameterized derived type in case we want a different size
+    type stack_t(ssize)
+        integer, len :: ssize
+        private
+        type(rpn_t)  :: sdata(ssize)
+        character(2) :: legend(ssize)
+        integer      :: high_water = 0
+    contains
+        procedure, private :: push_stackt
+        procedure, private :: push_all_stackt
+        procedure, private :: push_r_stackt
+        generic, public    :: push => push_stackt, push_all_stackt, push_r_stackt
+        procedure          :: peek => peek_stackt
+        procedure          :: pop  => pop_stackt
+        procedure          :: set  => set_stackt
+        procedure          :: clear => clear_stackt
+        procedure          :: swap  => swap_stackt
+        procedure          :: rotate_up  => rotate_up_stackt
+        procedure          :: rotate_down  => rotate_down_stackt
+        procedure          :: print => print_stackt
+        procedure          :: get_size => get_size_stackt
+        procedure          :: set_legend => set_legend_stackt
+    end type stack_t
+    
+    
+    interface
+        
+        module subroutine set_legend_stackt(stk, legend)
+            class(stack_t(*)), intent(inout) :: stk
+            character(len=2), intent(in)     :: legend(:)
+        end subroutine set_legend_stackt
+        module function get_size_stackt(stk) result(r)
+            class(stack_t(*)), intent(in) :: stk
+            integer :: r
+        end function get_size_stackt
+        module subroutine print_stackt(stk, ve_mode)
+            class(stack_t(*)), intent(in) :: stk
+            logical, intent(in)           :: ve_mode
+        end subroutine print_stackt
+        module subroutine push_stackt(stk, z)
+            class(stack_t(*)), intent(inout) :: stk
+            type(rpn_t) :: z
+        end subroutine push_stackt
+        module subroutine push_r_stackt(stk, x)
+            class(stack_t(*)), intent(inout) :: stk
+            real(8) :: x
+        end subroutine push_r_stackt
+        module subroutine push_all_stackt(stk, z, is_cart)
+            class(stack_t(*)), intent(inout) :: stk
+            complex(8), intent(in) :: z
+            logical, intent(in), optional :: is_cart
+        end subroutine push_all_stackt
+        module subroutine set_stackt(stk, z, idx)
+            class(stack_t(*)), intent(inout) :: stk
+            type(rpn_t), intent(in) :: z
+            integer, optional, intent(in) :: idx
+        end subroutine set_stackt
+        module function peek_stackt(stk, idx) result(r)
+            class(stack_t(*)), intent(inout) :: stk
+            integer, intent(in) :: idx 
+            type(rpn_t) :: r
+        end function peek_stackt
+        module function pop_stackt(stk) result(r)
+            class(stack_t(*)), intent(inout) :: stk
+            type(rpn_t) :: r
+        end function pop_stackt
+        module subroutine clear_stackt(stk)
+            class(stack_t(*)), intent(inout) :: stk
+        end subroutine clear_stackt
+        module subroutine swap_stackt(stk)
+            class(stack_t(*)), intent(inout) :: stk
+        end subroutine swap_stackt
+        module subroutine rotate_up_stackt(stk)
+            class(stack_t(*)), intent(inout) :: stk
+        end subroutine rotate_up_stackt
+        module subroutine rotate_down_stackt(stk)
+            class(stack_t(*)), intent(inout) :: stk
+        end subroutine rotate_down_stackt
+    end interface
+    
+    interface
+        module subroutine write_rpns(se, unit, iotype, v_list, iostat, iomsg)
+            class(rpn_t), intent(in)    :: se
+            integer, intent(in)         :: unit
+            character(*), intent(in)    :: iotype
+            integer, intent(in)         :: v_list(:)
+            integer, intent(out)        :: iostat
+            character(*), intent(inout) :: iomsg
+        end subroutine write_rpns
+        module function is_integer_rpns(this) result(r)
+            class(rpn_t), intent(in) :: this
+            logical :: r
+        end function is_integer_rpns
+        module function is_cartesian_rpns(this) result(r)
+            class(rpn_t), intent(in) :: this
+            logical :: r
+        end function is_cartesian_rpns
+        module function is_real_rpns(this) result(r)
+            class(rpn_t), intent(in) :: this
+            logical :: r
+        end function is_real_rpns
+        module function is_positive_real_rpns(this) result(r)
+            class(rpn_t), intent(in) :: this
+            logical :: r
+        end function is_positive_real_rpns
+        module subroutine set_angle_unit_rpns(this, degrees)
+            class(rpn_t), intent(inout) :: this
+            logical, intent(in) :: degrees
+        end subroutine set_angle_unit_rpns
+        module function get_value_rpns(this, is_cartesian) result(r)
+            class(rpn_t), intent(in)       :: this
+            logical, optional, intent(out) :: is_cartesian
+            complex(8) :: r
+        end function get_value_rpns
+        module subroutine set_value_rpns(this, z, is_cartesian)
+            class(rpn_t), intent(inout)      :: this
+            complex(8), optional, intent(in) :: z
+            logical, optional, intent(in)    :: is_cartesian
+        end subroutine set_value_rpns
+        module subroutine set_to_rpns(this, z)
+            class(rpn_t), intent(inout) :: this
+            type(rpn_t), intent(in) :: z
+        end subroutine set_to_rpns
+        module function add_rpns(a, b) result(r)
+            class(rpn_t), intent(in) :: a
+            type(rpn_t), intent(in)  :: b
+            type(rpn_t) :: r
+        end function add_rpns
+        module function subtract_rpns(a, b) result(r)
+            class(rpn_t), intent(in) :: a
+            type(rpn_t), intent(in)  :: b
+            type(rpn_t) :: r
+        end function subtract_rpns
+        module function multiply_rpns(a, b) result(r)
+            class(rpn_t), intent(in) :: a
+            type(rpn_t), intent(in)  :: b
+            type(rpn_t) :: r
+        end function multiply_rpns
+        module function divide_rpns(a, b) result(r)
+            class(rpn_t), intent(in) :: a
+            type(rpn_t), intent(in)  :: b
+            type(rpn_t) :: r
+        end function divide_rpns
+        module function power_rpns(this, x) result(r)
+            class(rpn_t), intent(in) :: this
+            real(8), intent(in) :: x
+            type(rpn_t) :: r
+        end function power_rpns        
+    end interface
+    
+    real(8), parameter :: pi = 4*atan(1.0d0)
+    real(8), parameter :: to_rad = pi/180
+    real(8), parameter :: to_deg = 180/pi
+
+    character(5), private :: decimal = 'POINT'
+    
+    integer                  :: nroots = 0
+    type(rpn_t), allocatable :: roots(:)
+    integer, private         :: current_root = 0
+    
+    character(9), private :: f_large
+    character(9), private :: f_small
+    integer               :: dec_places = 6
+    logical               :: degrees_mode = .true.
+    logical               :: complex_mode = .false.
+    real(8)               :: eps = 1.0d-14
+    
+    ! Functions interface
+    interface
+        module function add_fr(a,b) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t), intent(in) :: b
+            type(rpn_t) :: r
+        end function add_fr
+        
+        module function subtract_fr(a,b) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t), intent(in) :: b
+            type(rpn_t) :: r
+        end function subtract_fr
+        
+        module function multiply_fr(a,b) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t), intent(in) :: b
+            type(rpn_t) :: r
+        end function multiply_fr
+        
+        module function divide_fr(a,b) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t), intent(in) :: b
+            type(rpn_t) :: r
+        end function divide_fr
+
+        module function power_fr(a, b) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t), intent(in) :: b
+            type(rpn_t) :: r
+        end function power_fr
+        
+        module function percent_fr(a, b) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+            type(rpn_t), intent(in) :: b
+        end function percent_fr
+        
+        module function power_2_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function power_2_fr
+
+        module function power_3_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function power_3_fr
+
+        module function sqrt_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function sqrt_fr
+
+        module function cbrt_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function cbrt_fr
+
+        module function reciprocal_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function reciprocal_fr
+
+        module function conj_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function conj_fr
+    
+        module function len_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function len_fr
+        
+        module function swap_real_imaginary_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function swap_real_imaginary_fr
+        
+        module function chs_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function chs_fr
+        
+        module function sine_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function sine_fr
+
+        module function cosine_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function cosine_fr
+
+        module function tangent_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function tangent_fr
+
+        module function hsine_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function hsine_fr
+
+        module function hcosine_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function hcosine_fr
+
+        module function htangent_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function htangent_fr
+
+        module function asine_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function asine_fr
+
+        module function acosine_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function acosine_fr
+
+        module function atangent_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function atangent_fr
+
+        module function ahsine_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function ahsine_fr
+
+        module function ahcosine_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function ahcosine_fr
+
+        module function ahtangent_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function ahtangent_fr
+    
+        module function exp_2_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function exp_2_fr
+
+        module function exp_e_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function exp_e_fr
+
+        module function exp_10_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function exp_10_fr
+
+        module function ln_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function ln_fr
+
+        module function log2_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function log2_fr
+
+        module function lg_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function lg_fr
+
+        module function gamma_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function gamma_fr
+
+        module function fact_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function fact_fr
+
+        module function ncr_fr(a, b) result(r)
+            type(rpn_t), intent(in) :: a, b
+            type(rpn_t) :: r
+        end function ncr_fr
+
+        module function npr_fr(a, b) result(r)
+            type(rpn_t), intent(in) :: a, b
+            type(rpn_t) :: r
+        end function npr_fr
+        
+        module function root_fr(a, b) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t), intent(in) :: b
+            type(rpn_t) :: r
+        end function root_fr
+        
+        module function next_root_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function next_root_fr
+        
+        module function previous_root_fr(a) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function previous_root_fr
+
+        module function atangent2_fr(a, b) result(r)
+            type(rpn_t), intent(in) :: a
+            type(rpn_t), intent(in) :: b
+            type(rpn_t) :: r
+        end function atangent2_fr
+        
+        module function round(x) result(r)
+            real(8), intent(in) :: x
+            real(8) ::r
+        end function round
+        
+        module subroutine init(lang)
+            character(5), intent(in), optional :: lang
+        end subroutine init
+
+        module subroutine set_places(n)
+            integer, intent(in) :: n
+        end subroutine set_places
+        
+        module function get_places() result(r)
+            integer :: r
+        end function get_places
+        
+        module subroutine to_string(x, str)
+            real(8), intent(in) :: x
+            character(len=:), allocatable, intent(out) :: str
+        end subroutine to_string
+    end interface
+         
+    public
+    
+    interface to_polar
+        module function to_polar_rpns(stk_z) result(r)
+            type(rpn_t), intent(in) :: stk_z
+            type(rpn_t) :: r
+        end function to_polar_rpns
+    end interface
+    
+    interface to_cartesian
+        module function to_cartesian_rpns(stk_z) result(r)
+            type(rpn_t), intent(in) :: stk_z
+            type(rpn_t) :: r
+        end function to_cartesian_rpns
+    end interface
+    
+    abstract interface
+        function binary_f(a,b) result(r)
+            import
+            type(rpn_t), intent(in) :: a, b
+            type(rpn_t) :: r
+        end function binary_f
+        function unary_f(a) result(r)
+            import
+            type(rpn_t), intent(in) :: a
+            type(rpn_t) :: r
+        end function unary_f
+    end interface
+        
+end module rpn_stack

rpn_stack_sm.f90

@@ -0,0 +1,759 @@
+! Implementation code for stack
+submodule (rpn_stack) stack_sm
+    implicit none
+
+contains
+
+    module subroutine set_legend_stackt(stk, legend)
+        class(stack_t(*)), intent(inout) :: stk
+        character(len=2), intent(in)     :: legend(:)
+        stk%legend = legend
+    end subroutine set_legend_stackt
+
+    module function get_size_stackt(stk) result(r)
+        class(stack_t(*)), intent(in) :: stk
+        integer :: r
+        r = stk%high_water
+    end function get_size_stackt
+
+    module subroutine print_stackt(stk, ve_mode)
+        class(stack_t(*)), intent(in) :: stk
+        logical, intent(in)           :: ve_mode
+        integer :: i, j
+        if (ve_mode) then
+            do i=stk%high_water,1,-1
+                write(6,fmt='(a)',advance='no') stk%legend(i)//' '
+                write(6,'(dt)') stk%sdata(i)
+            end do
+        else
+            write(6,fmt='(dt)') stk%sdata(1)
+        end if
+    end subroutine print_stackt
+    
+    module subroutine push_stackt(stk, z)
+        class(stack_t(*)), intent(inout) :: stk
+        type(rpn_t) :: z
+        integer :: i
+        do i=stk%ssize,2,-1
+            stk%sdata(i) = stk%sdata(i-1)
+        end do
+        stk%sdata(1) = z
+        if (stk%high_water < stk%ssize) &
+            stk%high_water = stk%high_water + 1
+    end subroutine push_stackt
+    
+    module subroutine push_r_stackt(stk, x)
+        class(stack_t(*)), intent(inout) :: stk
+        real(8) :: x
+        type(rpn_t) :: z
+        z = rpn_t(cmplx(x,0.0d0))
+        call stk%push_stackt(z)
+    end subroutine push_r_stackt
+    
+    module subroutine push_all_stackt(stk, z, is_cart)
+        class(stack_t(*)), intent(inout) :: stk
+        complex(8), intent(in) :: z
+        logical, intent(in), optional :: is_cart
+        integer :: i
+        do i=stk%ssize,2,-1
+            stk%sdata(i) = stk%sdata(i-1)
+        end do
+        if (present(is_cart)) then
+            call stk%set(rpn_t(z,is_cart))
+        else
+            call stk%set(rpn_t(z))
+        end if
+        if (stk%high_water < stk%ssize) &
+            stk%high_water = stk%high_water + 1
+    end subroutine push_all_stackt
+
+    module subroutine set_stackt(stk, z, idx)
+        class(stack_t(*)), intent(inout) :: stk
+        type(rpn_t), intent(in) :: z
+        integer, optional, intent(in) :: idx
+        if (present(idx)) then
+            stk%sdata(idx) = z
+        else
+            stk%sdata(1) = z
+        end if
+    end subroutine set_stackt
+    
+    module function peek_stackt(stk, idx) result(r)
+        class(stack_t(*)), intent(inout) :: stk
+        integer, intent(in) :: idx 
+        type(rpn_t) :: r
+        if (idx >= 1 .and. idx <= stk%ssize) then
+            r = stk%sdata(idx)
+        else
+            write(*,'(a,i0,a)') '***Invalid index (',idx,')'
+            r = rpn_t()
+        end if
+    end function peek_stackt
+    
+    module function pop_stackt(stk) result(r)
+        class(stack_t(*)), intent(inout) :: stk
+        type(rpn_t) :: r
+        integer :: i
+        r = stk%sdata(1)
+        do i=1,stk%ssize-1
+            stk%sdata(i) = stk%sdata(i+1)
+        end do
+        stk%sdata(stk%ssize) = rpn_t()
+        if (stk%high_water > 0) &
+            stk%high_water = stk%high_water - 1
+    end function pop_stackt
+    
+    module subroutine clear_stackt(stk)
+        class(stack_t(*)), intent(inout) :: stk
+        integer :: i
+        do i=1,stk%ssize
+            stk%sdata(i) = rpn_t()
+        end do
+        stk%high_water = 0
+    end subroutine clear_stackt
+    
+    module subroutine swap_stackt(stk)
+        class(stack_t(*)), intent(inout) :: stk
+        integer :: i
+        type(rpn_t) :: z
+        z = stk%sdata(1)
+        stk%sdata(1) = stk%sdata(2)
+        stk%sdata(2) = z
+    end subroutine swap_stackt
+    
+    module subroutine rotate_up_stackt(stk)
+        class(stack_t(*)), intent(inout) :: stk
+        type(rpn_t) :: z
+        z = stk%pop()
+        stk%high_water = stk%high_water + 1
+        call stk%set(z,stk%high_water)
+    end subroutine rotate_up_stackt
+    
+    module subroutine rotate_down_stackt(stk)
+        class(stack_t(*)), intent(inout) :: stk
+        type(rpn_t) :: z
+        z = stk%peek(stk%high_water)
+        stk%high_water = stk%high_water - 1
+        call stk%push(z)
+    end subroutine rotate_down_stackt
+
+end submodule stack_sm
+
+! Implementation code for rpn_t
+submodule (rpn_stack) rpn_sm
+contains
+
+    module subroutine write_rpns(se, unit, iotype, v_list, iostat, iomsg)
+        class(rpn_t), intent(in)    :: se
+        integer, intent(in)         :: unit
+        character(*), intent(in)    :: iotype
+        integer, intent(in)         :: v_list(:)
+        integer, intent(out)        :: iostat
+        character(*), intent(inout) :: iomsg
+        complex(8) :: z
+        character(len=:), allocatable :: str_re, str_im
+        iostat = 0
+        iomsg = ""
+        z = se%zdata
+        if (complex_mode) then
+            call to_string(z%re,str_re)
+            call to_string(z%im,str_im)
+            if (se%is_cartesian()) then
+                write(6,'(a)') '('//str_re//','//str_im//')'
+            else
+                write(6,'(a)') '('//str_re//','//str_im//') p'
+            end if
+        else
+            call to_string(z%re,str_re)
+            write(6,'(a)') str_re
+        end if
+
+    end subroutine write_rpns
+    
+    ! Convert real to string inserting a leading 0 if necessary
+    module subroutine to_string(x, str)
+        real(8), intent(in) :: x
+        character(len=:), allocatable, intent(out) :: str
+        character(len=32) :: s
+        s = ' '
+        if (f_small == 'i0') then
+            write(s,fmt='('//f_small//')') nint(x)
+        else
+            if (x == 0 .or. (abs(x) < 1.0d7 .and. abs(x) > 1.0d-7)) then
+                write(s(2:),fmt='('//f_small//')') x
+            else
+                write(s(2:),fmt='('//f_large//')') x
+            end if
+            if (s(2:3) == '-.') then
+                s(1:3) = '-0.'
+            else if (s(2:2) == '.') then
+                s(1:2) = '0.'
+            end if
+        end if
+        str = trim(adjustl(s))
+    end subroutine to_string
+
+    module function is_integer_rpns(this) result(r)
+        class(rpn_t), intent(in) :: this
+        logical :: r
+        real(8) :: x
+        x = this%zdata%re
+        r = (abs(nint(x)-x) < eps .and. abs(this%zdata%im) < eps)
+    end function is_integer_rpns
+
+    module function is_cartesian_rpns(this) result(r)
+        class(rpn_t), intent(in) :: this
+        logical :: r
+        r = this%is_cart
+    end function is_cartesian_rpns
+
+    module function is_real_rpns(this) result(r)
+        class(rpn_t), intent(in) :: this
+        logical :: r
+        r = this%zdata%im == 0
+    end function is_real_rpns
+
+    module function is_positive_real_rpns(this) result(r)
+        class(rpn_t), intent(in) :: this
+        logical :: r
+        r = this%zdata%im == 0 .and. this%zdata%re > 0
+    end function is_positive_real_rpns
+
+    module subroutine set_angle_unit_rpns(this, degrees)
+        class(rpn_t), intent(inout) :: this
+        logical, intent(in) :: degrees
+        if (.not. this%is_cart) then
+            this%zdata%im = this%zdata%im*merge(to_deg, to_rad, degrees)
+        end if
+    end subroutine set_angle_unit_rpns
+
+    module function get_value_rpns(this, is_cartesian) result(r)
+        class(rpn_t), intent(in)       :: this
+        logical, optional, intent(out) :: is_cartesian
+        complex(8) :: r
+        r = this%zdata
+        if (present(is_cartesian)) then
+            is_cartesian = this%is_cart
+        end if
+    end function get_value_rpns
+
+    module subroutine set_value_rpns(this, z, is_cartesian)
+        class(rpn_t), intent(inout)      :: this
+        complex(8), optional, intent(in) :: z
+        logical, optional, intent(in)    :: is_cartesian
+        if (present(z)) then
+            this%zdata = z
+        end if
+        if (present(is_cartesian)) then
+            this%is_cart = is_cartesian
+        end if
+    end subroutine set_value_rpns
+
+    module subroutine set_to_rpns(this, z)
+        class(rpn_t), intent(inout) :: this
+        type(rpn_t), intent(in) :: z
+        this%zdata = z%zdata
+        this%is_cart = z%is_cart
+    end subroutine set_to_rpns
+
+    module function add_rpns(a, b) result(r)
+        class(rpn_t), intent(in) :: a
+        type(rpn_t), intent(in)  :: b
+        type(rpn_t) :: r
+        type(rpn_t) :: s
+        logical :: is_cart
+        is_cart = a%is_cartesian() ! The output will be set to this
+        if (a%is_cartesian()) then
+            r = a
+        else
+            r = to_cartesian(a)
+        end if
+        if (b%is_cartesian()) then
+            r%zdata = r%zdata + b%zdata
+        else
+            s = to_cartesian(a)
+            r%zdata = r%zdata + s%zdata
+        end if
+        if (.not. is_cart) then
+            r = to_polar(r)
+        end if
+    end function add_rpns
+    
+    module function subtract_rpns(a, b) result(r)
+        class(rpn_t), intent(in) :: a
+        type(rpn_t), intent(in)  :: b
+        type(rpn_t) :: r
+        type(rpn_t) :: s
+        logical :: is_cart
+        is_cart = a%is_cartesian() ! The output will be set to this
+        if (a%is_cartesian()) then
+            r = a
+        else
+            r = to_cartesian(a)
+        end if
+        if (b%is_cartesian()) then
+            r%zdata = r%zdata - b%zdata
+        else
+            s = to_cartesian(a)
+            r%zdata = r%zdata - s%zdata
+        end if
+        if (.not. is_cart) then
+            r = to_polar(r)
+        end if
+    end function subtract_rpns
+    
+    module function multiply_rpns(a, b) result(r)
+        class(rpn_t), intent(in) :: a
+        type(rpn_t), intent(in)  :: b
+        type(rpn_t) :: r
+        type(rpn_t) :: s
+        logical :: is_cart
+        is_cart = a%is_cartesian() ! The output will be set to this
+        if (a%is_cartesian()) then
+            r = a
+        else
+            r = to_cartesian(a)
+        end if
+        if (b%is_cartesian()) then
+            r%zdata = r%zdata * b%zdata
+        else
+            s = to_cartesian(a)
+            r%zdata = r%zdata * s%zdata
+        end if
+        if (.not. is_cart) then
+            r = to_polar(r)
+        end if
+    end function multiply_rpns
+    
+    module function divide_rpns(a, b) result(r)
+        class(rpn_t), intent(in) :: a
+        type(rpn_t), intent(in)  :: b
+        type(rpn_t) :: r
+        type(rpn_t) :: s
+        logical :: is_cart
+        is_cart = a%is_cartesian() ! The output will be set to this
+        if (a%is_cartesian()) then
+            r = a
+        else
+            r = to_cartesian(a)
+        end if
+        if (b%is_cartesian()) then
+            r%zdata = r%zdata / b%zdata
+        else
+            s = to_cartesian(a)
+            r%zdata = r%zdata / s%zdata
+        end if
+        if (.not. is_cart) then
+            r = to_polar(r)
+        end if
+    end function divide_rpns
+    
+    module function power_rpns(this, x) result(r)
+        class(rpn_t), intent(in) :: this
+        real(8), intent(in) :: x
+        type(rpn_t) :: r
+        type(rpn_t) :: z
+        logical     :: is_cart
+        is_cart = this%is_cartesian()
+        if (.not. is_cart) then
+            z = to_cartesian(this)
+        else
+            z = this
+        end if
+        r%zdata = z%zdata**x
+        if (.not. is_cart) then
+            r = to_polar(r)
+        end if
+    end function power_rpns
+    
+    module function to_cartesian_rpns(stk_z) result(r)
+        type(rpn_t), intent(in) :: stk_z
+        type(rpn_t) :: r
+        real(8) :: s
+        real(8) :: theta
+        if (.not. stk_z%is_cartesian()) then
+            s = stk_z%zdata%re
+            theta = stk_z%zdata%im * merge(to_rad,1.0d0,degrees_mode)
+            r%zdata%re = round(s * cos(theta))
+            r%zdata%im = round(s * sin(theta))
+            r%is_cart = .true.
+        else
+            r = stk_z
+        end if
+    end function to_cartesian_rpns
+
+        
+    module function to_polar_rpns(stk_z) result(r)
+        type(rpn_t), intent(in) :: stk_z
+        type(rpn_t) :: r
+        real(8) :: theta
+        if (stk_z%is_cartesian()) then
+            call r%set_value(to_polar_internal(stk_z%get_value()),is_cartesian = .false.)
+        else
+            r = stk_z
+        end if
+    contains
+        complex(8) function to_polar_internal(z)
+            complex(8), intent(in) :: z
+            real(8) :: r
+            real(8) :: theta
+            r = sqrt(real(z * conjg(z),8))
+            theta = atan2(aimag(z), real(z))
+            to_polar_internal%re = r
+            to_polar_internal%im = theta * merge(1/to_rad,1.0d0,degrees_mode)
+        end function to_polar_internal
+    end function to_polar_rpns
+    
+    module function add_fr(a,b) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t), intent(in) :: b
+        type(rpn_t) :: r
+        r = a + b
+    end function add_fr
+
+    module function subtract_fr(a,b) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t), intent(in) :: b
+        type(rpn_t) :: r
+        r = a - b
+    end function subtract_fr
+    
+    module function multiply_fr(a,b) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t), intent(in) :: b
+        type(rpn_t) :: r
+        r = a * b
+    end function multiply_fr
+    
+    module function divide_fr(a,b) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t), intent(in) :: b
+        type(rpn_t) :: r
+        r = a / b
+    end function divide_fr
+    
+    module function percent_fr(a,b) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t), intent(in) :: b
+        type(rpn_t) :: r
+        r = a * b / rpn_t(cmplx(100.0d0,0.0d0))
+    end function percent_fr
+
+    module function power_fr(a, b) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t), intent(in) :: b
+        type(rpn_t) :: r
+        r = a ** real(b%zdata)
+    end function power_fr
+    
+    module function power_2_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        r = a ** 2.0d0
+    end function power_2_fr
+
+    module function power_3_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        r = a ** 3.0d0
+    end function power_3_fr
+
+    module function sqrt_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        call r%set_value(sqrt(a%zdata))
+    end function sqrt_fr
+
+    module function cbrt_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        call r%set_value(a%zdata ** (1.0d0/3))
+    end function cbrt_fr
+
+    module function reciprocal_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        r = rpn_t(1.0d0)/a
+    end function reciprocal_fr
+
+    module function conj_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        r = a
+        r%zdata%im = -r%zdata%im
+    end function conj_fr
+    
+    module function len_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        type(rpn_t) :: s
+        s = a * conj_fr(a)
+        r = rpn_t(cmplx(sqrt(real(s%zdata%re)),0.0d0,8))
+    end function len_fr
+    
+    module function swap_real_imaginary_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        real(8) :: x
+        r = a
+        x = r%zdata%re
+        r%zdata%re = r%zdata%im
+        r%zdata%im = x
+    end function swap_real_imaginary_fr
+    
+    module function chs_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        r = rpn_t(cmplx(-a%zdata%re,-a%zdata%im))
+    end function chs_fr
+    
+    module function sine_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        r = rpn_t(sin(a%zdata * merge(to_rad,1.0d0,degrees_mode)))
+    end function sine_fr
+
+    module function cosine_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        r = rpn_t(cos(a%zdata%re * merge(to_rad,1.0d0,degrees_mode)))
+    end function cosine_fr
+
+    module function tangent_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        r = rpn_t(tan(a%zdata * merge(to_rad,1.0d0,degrees_mode)))
+    end function tangent_fr
+
+    module function hsine_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        r = rpn_t(sinh(a%zdata))
+    end function hsine_fr
+
+    module function hcosine_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        r = rpn_t(cosh(a%zdata))
+    end function hcosine_fr
+
+    module function htangent_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        r = rpn_t(tanh(a%zdata))
+    end function htangent_fr
+
+    module function asine_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        r = rpn_t(asin(a%zdata) * merge(1/to_rad,1.0d0,degrees_mode))
+    end function asine_fr
+
+    module function acosine_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        r = rpn_t(acos(a%zdata) * merge(1/to_rad,1.0d0,degrees_mode))
+    end function acosine_fr
+
+    module function atangent_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        r = rpn_t(atan(a%zdata) * merge(1/to_rad,1.0d0,degrees_mode))
+    end function atangent_fr
+
+    module function ahsine_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        r = rpn_t(asinh(a%zdata))
+    end function ahsine_fr
+
+    module function ahcosine_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        r = rpn_t(acosh(a%zdata))
+    end function ahcosine_fr
+
+    module function ahtangent_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        r = rpn_t(atanh(a%zdata))
+    end function ahtangent_fr
+ 
+    module function exp_2_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        r = rpn_t(2**a%zdata)
+    end function exp_2_fr
+
+    module function exp_e_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        r = rpn_t(exp(a%zdata))
+    end function exp_e_fr
+
+    module function exp_10_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        r = rpn_t(10**a%zdata)
+    end function exp_10_fr
+
+    module function ln_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        r = rpn_t(log(a%zdata))
+    end function ln_fr
+
+    module function log2_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        r = rpn_t(log(a%zdata)/log(2.0d0))
+    end function log2_fr
+
+    module function lg_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        r = rpn_t(log(a%zdata)/log(10.0d0))
+    end function lg_fr
+
+    module function gamma_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        r = rpn_t(gamma(a%zdata%re))
+    end function gamma_fr
+
+    module function fact_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        if (a%zdata%re == 0) then
+            r = rpn_t(1)
+        else
+            r = rpn_t(a%zdata%re*gamma(a%zdata%re))
+        end if
+    end function fact_fr
+
+    module function ncr_fr(a, b) result(r)
+        type(rpn_t), intent(in) :: a, b
+        type(rpn_t) :: r
+        r = fact_fr(a)/(fact_fr(b)*fact_fr(a-b))
+    end function ncr_fr
+
+    module function npr_fr(a, b) result(r)
+        type(rpn_t), intent(in) :: a, b
+        type(rpn_t) :: r
+        r = fact_fr(a)/fact_fr(b)
+    end function npr_fr
+
+    module function root_fr(a, b) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t), intent(in) :: b
+        type(rpn_t) :: r
+        real(8)     :: bc
+        integer     :: i
+        type(rpn_t) :: base
+        complex(8)  :: z
+        logical     :: a_is_cart
+        real(8)     :: s, delta_theta, theta0, phi
+        real(8), parameter :: two_pi = 8*atan(1.0d0)
+       
+        bc = real(b%get_value())
+        r = power_fr(a, rpn_t(1.0d0/bc))
+        ! If b is an integer >= 2 calculate all roots
+        if (b%is_integer() .and. bc >= 2) then
+            nroots = nint(bc)
+            if (allocated(roots)) then
+                deallocate(roots)
+            end if
+            a_is_cart = a%is_cartesian()
+            base = to_polar_rpns(a)
+            z = base%get_value()
+            s = z%re ** (1.0d0/bc)
+            theta0 = merge(z%im*to_rad,z%im,degrees_mode)/nroots
+            delta_theta = two_pi/nroots
+            allocate(roots(nroots))
+            do i=1, nroots
+                phi = theta0 + (i-1)*delta_theta
+                if (a_is_cart) then
+                    roots(i) = rpn_t(cmplx(round(s*cos(phi)),round(s*sin(phi))),a_is_cart)
+                else
+                    if (degrees_mode) phi = phi*to_deg
+                    roots(i) = rpn_t(cmplx(s,round(phi)),a_is_cart)
+                end if
+            end do
+            r = roots(1)
+            current_root = 1
+        else
+            r = a ** (1.0d0/bc)
+        end if
+        
+    end function root_fr
+
+    module function next_root_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        if (nroots > 0) then
+            if (current_root == nroots) then
+                current_root = 1
+            else
+                current_root = current_root + 1
+            end if
+            r = roots(current_root)
+        else
+            r = a
+        end if
+    end function next_root_fr
+
+    module function previous_root_fr(a) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t) :: r
+        if (nroots > 0) then
+            if (current_root == 1) then
+                current_root = nroots
+            else
+                current_root = current_root - 1
+            end if
+            r = roots(current_root)
+        else
+            r = a
+        end if
+    end function previous_root_fr
+
+    module function atangent2_fr(a, b) result(r)
+        type(rpn_t), intent(in) :: a
+        type(rpn_t), intent(in) :: b
+        type(rpn_t) :: r
+        r%zdata = atan2(real(a%zdata),real(b%zdata)) * merge(to_deg,1.0d0,degrees_mode)
+    end function atangent2_fr
+    
+    module function round(x) result(r)
+        real(8), intent(in) :: x
+        real(8) :: r
+        if (abs(x) < eps) then
+            r = 0
+        else
+            r = x
+        end if
+    end function round
+        
+    module subroutine init(lang)
+        character(5), intent(in), optional :: lang
+        if (present(lang)) decimal = lang
+        call set_places(dec_places)
+    end subroutine init
+
+    module subroutine set_places(n)
+        integer, intent(in) :: n
+        if (n == 0) then
+            f_small = 'i0'
+        else
+            write(f_small,'(2(a,i0),a)')  'f',0,'.',n
+            write(f_large,'(2(a,i0),a)') 'en',10+n,'.',n
+        end if
+        dec_places = n
+    end subroutine set_places
+    
+    module function get_places() result(r)
+        integer :: r
+        r = dec_places
+    end function get_places
+    
+end submodule rpn_sm

test_amap.f90

@@ -0,0 +1,18 @@
+program test_amap
+    use amap
+    implicit none
+    
+    type(amap_t) :: my_amap
+    type(value_t) :: x
+    
+    call my_amap%set('one',1.0d0)
+    call my_amap%set('two',2.d0)
+    call my_amap%set('three',3.0d0)
+    call my_amap%set('four',4.0d0)
+    call my_amap%set('five',5.0d0)
+    call my_amap%print
+    
+    x = my_amap%get('four')
+    write(6,'(f0.6)') my_amap%get_value('four')
+    print *,my_amap%contains('one'),my_amap%contains('ten')
+end program test_amap