!
!    Copyright 2020, 2021 Guy Munhoven
!
!    This file is part of SolveSAPHE v. 2

!    SolveSAPHE is free software: you can redistribute it and/or modify
!    it under the terms of the GNU Lesser General Public License as published by
!    the Free Software Foundation, either version 3 of the License, or
!    (at your option) any later version.
!
!    SolveSAPHE is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU Lesser General Public License for more details.
!
!    You should have received a copy of the GNU Lesser General Public License
!    along with SolveSAPHE.  If not, see <http://www.gnu.org/licenses/>.
!

program sdt


implicit none

integer, parameter :: ilevel_verbose = 0
                                    ! Set the following to true
logical, parameter :: lflag_printzeroes = .true.
logical, parameter :: lflag_printtexcode = .false.

integer :: i, j, k
integer :: n, nn

integer, allocatable, dimension(:)   :: iarr_m, iarr_f
integer, allocatable, dimension(:,:) :: iarr_prodpi
integer,              dimension(3)   :: iswap3
integer                              :: iswap

integer :: ijk, idx, idx_sorted, idx1, idx2
integer :: ijk_test
integer :: m, m_max, m_min, m_curr
integer, allocatable, dimension(:) :: n_m

integer, allocatable, dimension(:)   :: ijk_sorted

integer :: ixp_corr
integer :: ic, ic_curr
integer :: ib, ib_min, ib_max
integer :: ia, ia_min, ia_max
integer :: i_aplusb, idelta_c
integer, allocatable, dimension(:)   :: iarr_fcomp
integer :: i_fcomp

character(len=64) :: c_fmt


write(*,'(A)', advance="no") 'Please enter number of protons in the acid: '
read*, n

nn = (n+1) * ((n+1)*n) / 2 ! k : n+1
                           ! i, j: n + (n-1) + ... + 1
m_max = 3*n - 1  ! for (i,j,k) = (0,1,0)
m_min = 1        ! for (i,j,k) = (n-1, n, n)



allocate(iarr_m(nn))
allocate(iarr_f(nn))
allocate(iarr_prodpi(3,nn))

allocate(ijk_sorted(nn))

allocate(n_m(m_min:m_max))

ijk = 1
do k = 0, n
  do i = 0, n-1
    do j = i + 1, n
      iarr_m(ijk) = 3*n - i - j - k   ! exponent
      iarr_f(ijk) = (i - j)**2 * (i + j + 1 -2*k)
      iarr_prodpi(1,ijk) = i
      iarr_prodpi(2,ijk) = j
      iarr_prodpi(3,ijk) = k
      ijk_sorted(ijk) = ijk
      ijk = ijk + 1
    enddo
  enddo
enddo


if (ilevel_verbose > 1) THEN
  do ijk = 1, nn
    write(*,'(3I3, " : ", i5, " * H^", I0 )') iarr_prodpi(:,ijk), iarr_f(ijk), iarr_m(ijk)
  enddo
endif


! order iarr_prodpi(:, ijk) for each ijk

do ijk = 1, nn

  if (iarr_prodpi(1,ijk) >  iarr_prodpi(2,ijk)) then
    iswap              = iarr_prodpi(1,ijk)
    iarr_prodpi(1,ijk) = iarr_prodpi(2,ijk)
    iarr_prodpi(2,ijk) = iswap
  endif

  if (iarr_prodpi(2,ijk) >  iarr_prodpi(3,ijk)) then
    iswap              = iarr_prodpi(2,ijk)
    iarr_prodpi(2,ijk) = iarr_prodpi(3,ijk)
    iarr_prodpi(3,ijk) = iswap
  endif

  if (iarr_prodpi(1,ijk) >  iarr_prodpi(2,ijk)) then
    iswap              = iarr_prodpi(1,ijk)
    iarr_prodpi(1,ijk) = iarr_prodpi(2,ijk)
    iarr_prodpi(2,ijk) = iswap
  endif

enddo



ijk = 1

n_m(:) = 0

do m = m_max, m_min, -1

  do idx_sorted = ijk, nn

    idx = ijk_sorted(idx_sorted)

    if (iarr_m(idx) == m) then

      if (idx == ijk) then

        n_m(m) = n_m(m) + 1
        ijk = ijk +1
        cycle

      else

        ijk_sorted(idx_sorted) = ijk_sorted(ijk)
        ijk_sorted(ijk)        = idx

        n_m(m) = n_m(m) + 1
        ijk = ijk +1

      endif

    endif      

  enddo

enddo

if (ilevel_verbose > 1) THEN
  write(*,'()')
  write(*,'()')
  
  m_curr = m_max
  do ijk = 1, nn
    idx = ijk_sorted(ijk)
    if (iarr_m(idx) /= m_curr) then
      write(*,'()')
      m_curr = iarr_m(idx)
    endif
    if ((iarr_f(idx) /= 0) .or. lflag_printzeroes) &
    write(*,'(3I3, " : ", i5, " * H^", I0 )') iarr_prodpi(:,idx), iarr_f(idx), iarr_m(idx)
  enddo
endif

! Within each exponent:
!   sort iarr_prodpi(3, :) ascending
!     then sort iarr_prodpi(2, :) ascending within each of the  iarr_prodpi(3, :)
!       then sort iarr_prodpi(1, :) ascending within each of the  iarr_prodpi(2, :)

ijk = 0

do m = m_max, m_min, -1

  do i = 1, n_m(m) - 1
    do j = ijk + 1, ijk + n_m(m) - i
      idx1 = ijk_sorted(j)
      idx2 = ijk_sorted(j+1)
      if     (iarr_prodpi(3,idx2)  < iarr_prodpi(3,idx1)) then
        ijk_sorted(j+1) = idx1
        ijk_sorted(j)   = idx2
      elseif (iarr_prodpi(3,idx2) == iarr_prodpi(3,idx1)) then
        if     (iarr_prodpi(2,idx2)  < iarr_prodpi(2,idx1)) then
          ijk_sorted(j+1) = idx1
          ijk_sorted(j)   = idx2
        elseif (iarr_prodpi(2,idx2) == iarr_prodpi(2,idx1)) then
          if      (iarr_prodpi(2,idx2) < iarr_prodpi(2,idx1)) then
            ijk_sorted(j+1) = idx1
            ijk_sorted(j)   = idx2
          endif
        endif
      endif
    enddo 
  enddo
  ijk = ijk + n_m(m)

enddo


if (ilevel_verbose > 0) THEN
  write(*,'()')
  write(*,'()')
  
  m_curr = m_max
  do ijk = 1, nn
    idx = ijk_sorted(ijk)
    if (iarr_m(idx) /= m_curr) then
      write(*,'()')
      m_curr = iarr_m(idx)
    endif
    if ((iarr_f(idx) /= 0) .or. lflag_printzeroes) &
    write(*,'(3I3, " : ", i5, " * H^", I0 )') iarr_prodpi(:,idx), iarr_f(idx), iarr_m(idx)
  enddo
endif


! Within each exponent: simplify the triplets

ijk = 0

do m = m_max, m_min, -1

  do i = ijk + 1, ijk + n_m(m) - 1
    idx1 = ijk_sorted(i)
    idx2 = ijk_sorted(i+1)
    if (ALL(iarr_prodpi(:,idx1) == iarr_prodpi(:,idx2))) then
      iarr_f(idx2) = iarr_f(idx2) + iarr_f(idx1)
      iarr_f(idx1) = 0   ! Set iarr_f(idx1) to zero
      iarr_m(idx1) = -1  ! and flag the term idx1 obsolete
    endif
  enddo
  ijk = ijk + n_m(m)

enddo

i = maxval(iarr_m(:))
j = maxval(abs(iarr_f(:)))

i = int(log10(dble(i))) + 1
j = int(log10(dble(j))) + 1 + 1  ! plus 1 for the sign
k = int(log10(dble(n))) + 1

write(c_fmt, '(a, i0, a, i0, ".", i0, a, i0, ".", i0, a)') &
  '(2x, i', j, ', 1x, 3("Pi_", i', k, k, ', 1x), "* H^", i', i, i, ')'


write(*,'()')
write(*,'()')
if (lflag_printtexcode) then
  write(*,'("$n = ", i0, "$")') n
else
  write(*,'("n = ", i0)') n
endif
write(*,'()')

m_curr = m_max
do ijk = 1, nn
  idx = ijk_sorted(ijk)
  if (iarr_m(idx) == -1) cycle      ! Skip the obsolete terms
  if (iarr_m(idx) /= m_curr) then
    write(*,'()')
    m_curr = iarr_m(idx)
  endif
  
  if ((iarr_f(idx) /= 0) .or. lflag_printzeroes) then
    if (lflag_printtexcode) then
      write(*,'("$", i0, 1x, 3("\Pi_{",i0,"} "), "H^{", I0, "}$")') iarr_f(idx), iarr_prodpi(:,idx), iarr_m(idx)
    else
      write(*,c_fmt) iarr_f(idx), iarr_prodpi(:,idx), iarr_m(idx)
      !write(*,'(3I3, " : ", i5, " * H^", I0 )') iarr_prodpi(:,idx), iarr_f(idx), iarr_m(idx)
    endif
  endif
enddo

!~ write(*,'()')
!~ do m = m_max, m_min, -1
!~   write(*,'(i3, 2x, i3)') m, n_m(m)
!~ enddo


allocate(iarr_fcomp(nn))

m_curr  = -1
ic_curr = -1
do ijk = 1, nn

  idx = ijk_sorted(ijk)
  if (iarr_m(idx) == -1) cycle      ! Skip the obsolete terms

                                    ! New exponent - reset ic
  if (iarr_m(idx) /= m_curr) then
    if (ilevel_verbose > 1) write(*,'()')
    m_curr = iarr_m(idx)
    ic_curr = - 1
  endif

  if (iarr_prodpi(3,idx) /= ic_curr) then

    ic_curr = iarr_prodpi(3,idx)
    i_aplusb = 3*n - m_curr - ic_curr
    ib_min = iarr_prodpi(2,idx)
    ia_max = iarr_prodpi(1,idx)

    idelta_c = ic_curr - ib_min

    if (ib_min == ia_max) then
                                    ! (b,b,c) chain
      if (ib_min < idelta_c) then
        ia_min = 0
        ib_max = i_aplusb
      else
        ib_max = ic_curr
        ia_min = i_aplusb - ib_max
      endif

    else
                                    ! (b-1,b,c) chain
      if ((ib_min -1) < idelta_c) then
        ia_min = 0
        ib_max = i_aplusb
      else
        ib_max = ic_curr
        ia_min = i_aplusb - ib_max
      endif

    endif
  
!~     write(*,'(i3)') i_aplusb
!~     write(*,'(i3, 2x, 3i3)') m_curr, ia_min, ib_max, ic_curr
!~     write(*,'(i3, 2x, 3i3)') m_curr, ia_max, ib_min, ic_curr
    iarr_fcomp(idx) = iarr_f(idx)
    i_fcomp = 0
    ib      = ib_min

  endif

  if (ilevel_verbose > 1) THEN
    write(*,'(i3, 2x, 3i3, 2x, i12, 2x, i12)') m_curr, i_aplusb - ib, ib, ic_curr, iarr_f(idx), iarr_fcomp(idx)
  endif

  if (ib < ib_max) then
                                    ! Search for the next non-obsolete term
    ijk_test = ijk+1

    do
      idx1 = ijk_sorted(ijk_test)
      if (iarr_m(idx1) /= -1) exit
      ijk_test = ijk_test+1
    enddo

    if (ia_max == ib_min) then      ! (b,b,c) chain
      ixp_corr = 2*(ib-ib_min)+1
    else                            ! (b-1,b,c) chain
      ixp_corr = 2*(ib-ib_min)+2
    endif
    i_fcomp = iarr_fcomp(idx) * 2**(ixp_corr)
    if (iarr_fcomp(idx) > (huge(1)-1)/(2**(ixp_corr))) then
      write (*,'("Overflow error: ", i0, " * ", i0, " /= ", i0)') &
        iarr_fcomp(idx), 2**(ixp_corr), i_fcomp
    endif


    iarr_fcomp(idx1) = iarr_f(idx1) + i_fcomp
    iarr_fcomp(idx)  = 0
    if (ilevel_verbose > 1) THEN
      write(*,'(3x, 2x, 3(3x), 2x, i12, 2x, i12, i3)') iarr_f(idx1), i_fcomp, ixp_corr
    endif
    if (i_fcomp >=0) then
      if ((huge(1) - i_fcomp) < iarr_fcomp(idx1)) then
        write(*,'()')
        write(*,'("Overflow error: ", i0, " + ", i0, " /= ", i0)') &
          iarr_f(idx1), i_fcomp, iarr_fcomp(idx1)
      endif
    endif

    ib = ib + 1

  else

    if (iarr_fcomp(idx) < 0) then

      write(*,'("iarr_fcomp @ b_max < 0!: ", i12, i12)') iarr_f(idx), iarr_fcomp(idx)

    endif

  endif

enddo


i = maxval(iarr_m(:))
j = maxval(iarr_fcomp(:))

i = int(log10(dble(i))) + 1
j = int(log10(dble(j))) + 1
k = int(log10(dble(n))) + 1

write(c_fmt, '(a, i0, a, i0, ".", i0, a, i0, ".", i0, a)') &
  '(2x i', j, ', 1x, 3("Pi_", i', k, k, ', 1x), "* H^", i', i, i, ')'

write(*,'()')
write(*,'()')
write(*,'("Minoring polynomial (assuming K_{i+1} < K_{i}/2, i = 1, ..., ", i0, ")")') n
write(*,'()')

m_curr = m_max
do ijk = 1, nn
  idx = ijk_sorted(ijk)
  if (iarr_m(idx) == -1) cycle      ! Skip the obsolete terms
  if (iarr_m(idx) /= m_curr) then
    write(*,'()')
    m_curr = iarr_m(idx)
  endif
  
  if ((iarr_fcomp(idx) /= 0) .or. lflag_printzeroes) then
    if (lflag_printtexcode) then
      write(*,'("$", i0, 1x, 3("\Pi_{",i0,"} "), "H^{", I0, "}$")') iarr_fcomp(idx), iarr_prodpi(:,idx), iarr_m(idx)
    else
      write(*,c_fmt) iarr_fcomp(idx), iarr_prodpi(:,idx), iarr_m(idx)
      !write(*,'(3I3, " : ", i5, " * H^", I0 )') iarr_prodpi(:,idx), iarr_fcomp(idx), iarr_m(idx)
    endif
  endif
enddo




end

