From 7cb24ceba357b2df9389c867b7878e50b6895841 Mon Sep 17 00:00:00 2001 From: sgeard Date: Fri, 9 Jun 2023 21:07:25 +0100 Subject: [PATCH] Added files --- GNUmakefile | 49 +++ amap.f90 | 195 ++++++++++ hp.f90 | 920 +++++++++++++++++++++++++++++++++++++++++++++++ linked_list.f90 | 97 +++++ rpn_stack.f90 | 467 ++++++++++++++++++++++++ rpn_stack_sm.f90 | 759 ++++++++++++++++++++++++++++++++++++++ test_amap.f90 | 18 + 7 files changed, 2505 insertions(+) create mode 100644 GNUmakefile create mode 100644 amap.f90 create mode 100644 hp.f90 create mode 100644 linked_list.f90 create mode 100644 rpn_stack.f90 create mode 100644 rpn_stack_sm.f90 create mode 100644 test_amap.f90 diff --git a/GNUmakefile b/GNUmakefile new file mode 100644 index 0000000..6daea38 --- /dev/null +++ b/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)" diff --git a/amap.f90 b/amap.f90 new file mode 100644 index 0000000..4732a4e --- /dev/null +++ b/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 diff --git a/hp.f90 b/hp.f90 new file mode 100644 index 0000000..028bde8 --- /dev/null +++ b/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 sw rc cl m+ m- m* m/ 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 diff --git a/linked_list.f90 b/linked_list.f90 new file mode 100644 index 0000000..c3d4a97 --- /dev/null +++ b/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 diff --git a/rpn_stack.f90 b/rpn_stack.f90 new file mode 100644 index 0000000..8380236 --- /dev/null +++ b/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 diff --git a/rpn_stack_sm.f90 b/rpn_stack_sm.f90 new file mode 100644 index 0000000..a830576 --- /dev/null +++ b/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 diff --git a/test_amap.f90 b/test_amap.f90 new file mode 100644 index 0000000..c570c7f --- /dev/null +++ b/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