module bemsa
  implicit none
  integer, parameter :: wp  = kind(1.0D0)
  
  private
  public :: bemsav, emsav
contains
  function emsav(x,c) result(v)
    implicit none
    real(wp),dimension(1:6)::x
    real(wp),dimension(0:195)::c
    real(wp)::v
    ! ::::::::::::::::::::
    real(wp),dimension(0:195)::p
    
    call bemsav(x,p)
    v = dot_product(p,c)
    
    return
  end function emsav
  
  subroutine bemsav(x,p)
    implicit none
    real(wp),dimension(1:6),intent(in)::x
    real(wp),dimension(0:195),intent(out)::p
    ! ::::::::::::::::::::
    real(wp),dimension(0:32)::m
    
    call evmono(x,m)
    call evpoly(m,p)
    
    return
  end subroutine bemsav
  
  subroutine evmono(x,m)
    implicit none
    real(wp),dimension(1:6),intent(in)::x
    real(wp),dimension(0:32),intent(out)::m
    !::::::::::::::::::::
    
    m(0) = 1.0D0
    m(1) = x(6)
    m(2) = x(5)
    m(3) = x(3)
    m(4) = x(4)
    m(5) = x(2)
    m(6) = x(1)
    m(7) = m(1)*m(2)
    m(8) = m(1)*m(3)
    m(9) = m(2)*m(3)
    m(10) = m(3)*m(4)
    m(11) = m(2)*m(5)
    m(12) = m(1)*m(6)
    m(13) = m(4)*m(5)
    m(14) = m(4)*m(6)
    m(15) = m(5)*m(6)
    m(16) = m(1)*m(9)
    m(17) = m(1)*m(10)
    m(18) = m(2)*m(10)
    m(19) = m(1)*m(11)
    m(20) = m(3)*m(11)
    m(21) = m(2)*m(12)
    m(22) = m(3)*m(12)
    m(23) = m(2)*m(13)
    m(24) = m(3)*m(13)
    m(25) = m(1)*m(14)
    m(26) = m(3)*m(14)
    m(27) = m(1)*m(15)
    m(28) = m(2)*m(15)
    m(29) = m(4)*m(15)
    m(30) = m(2)*m(24)
    m(31) = m(1)*m(26)
    m(32) = m(1)*m(28)

    return
  end subroutine evmono

  subroutine evpoly(m,p)
    implicit none
    real(wp),dimension(0:32),intent(in)::m
    real(wp),dimension(0:195),intent(out)::p
    !::::::::::::::::::::
    
    p(0) = m(0)
    p(1) = m(1) + m(2) + m(3)
    p(2) = m(4) + m(5) + m(6)
    p(3) = m(7) + m(8) + m(9)
    p(4) = m(10) + m(11) + m(12)
    p(5) = p(1)*p(2) - p(4)
    p(6) = m(13) + m(14) + m(15)
    p(7) = p(1)*p(1) - p(3) - p(3)
    p(8) = p(2)*p(2) - p(6) - p(6)
    p(9) = m(16)
    p(10) = m(17) + m(18) + m(19) + m(20) + m(21) + m(22)
    p(11) = p(2)*p(3) - p(10)
    p(12) = m(23) + m(24) + m(25) + m(26) + m(27) + m(28)
    p(13) = m(29)
    p(14) = p(1)*p(6) - p(12)
    p(15) = p(1)*p(3) - p(9) - p(9) - p(9)
    p(16) = p(1)*p(4) - p(10)
    p(17) = p(2)*p(7) - p(16)
    p(18) = p(2)*p(4) - p(12)
    p(19) = p(1)*p(8) - p(18)
    p(20) = p(2)*p(6) - p(13) - p(13) - p(13)
    p(21) = p(1)*p(7) - p(15)
    p(22) = p(2)*p(8) - p(20)
    p(23) = p(9)*p(2)
    p(24) = m(30) + m(31) + m(32)
    p(25) = p(3)*p(6) - p(24)
    p(26) = p(13)*p(1)
    p(27) = p(9)*p(1)
    p(28) = p(3)*p(4) - p(23)
    p(29) = p(1)*p(10) - p(23) - p(28) - p(23)
    p(30) = p(1)*p(11) - p(23)
    p(31) = p(1)*p(12) - p(25) - p(24) - p(24)
    p(32) = p(1)*p(14) - p(25)
    p(33) = p(4)*p(5) - p(25) - p(31)
    p(34) = p(2)*p(11) - p(25)
    p(35) = p(4)*p(6) - p(26)
    p(36) = p(2)*p(12) - p(26) - p(35) - p(26)
    p(37) = p(13)*p(2)
    p(38) = p(2)*p(14) - p(26)
    p(39) = p(3)*p(3) - p(27) - p(27)
    p(40) = p(3)*p(7) - p(27)
    p(41) = p(1)*p(16) - p(28)
    p(42) = p(2)*p(21) - p(41)
    p(43) = p(1)*p(18) - p(33)
    p(44) = p(7)*p(8) - p(43)
    p(45) = p(6)*p(6) - p(37) - p(37)
    p(46) = p(2)*p(18) - p(35)
    p(47) = p(1)*p(22) - p(46)
    p(48) = p(6)*p(8) - p(37)
    p(49) = p(1)*p(21) - p(40)
    p(50) = p(2)*p(22) - p(48)
    p(51) = p(9)*p(6)
    p(52) = p(13)*p(3)
    p(53) = p(9)*p(4)
    p(54) = p(9)*p(5)
    p(55) = p(1)*p(24) - p(51)
    p(56) = p(3)*p(12) - p(51) - p(55) - p(51)
    p(57) = p(3)*p(14) - p(51)
    p(58) = p(13)*p(7)
    p(59) = p(9)*p(8)
    p(60) = p(2)*p(24) - p(52)
    p(61) = p(13)*p(4)
    p(62) = p(4)*p(14) - p(58)
    p(63) = p(6)*p(11) - p(52)
    p(64) = p(13)*p(5)
    p(65) = p(9)*p(3)
    p(66) = p(9)*p(7)
    p(67) = p(3)*p(16) - p(53)
    p(68) = p(4)*p(15) - p(54) - p(67)
    p(69) = p(2)*p(39) - p(68)
    p(70) = p(1)*p(29) - p(54) - p(68)
    p(71) = p(7)*p(11) - p(53)
    p(72) = p(1)*p(31) - p(56) - p(55)
    p(73) = p(1)*p(32) - p(57)
    p(74) = p(3)*p(18) - p(59)
    p(75) = p(2)*p(29) - p(57) - p(55)
    p(76) = p(1)*p(34) - p(59)
    p(77) = p(6)*p(16) - p(58)
    p(78) = p(1)*p(36) - p(63) - p(60)
    p(79) = p(2)*p(32) - p(58)
    p(80) = p(6)*p(12) - p(64) - p(61) - p(61)
    p(81) = p(13)*p(6)
    p(82) = p(1)*p(45) - p(80)
    p(83) = p(2)*p(33) - p(62) - p(60)
    p(84) = p(2)*p(34) - p(63)
    p(85) = p(6)*p(18) - p(61)
    p(86) = p(2)*p(36) - p(64) - p(80)
    p(87) = p(13)*p(8)
    p(88) = p(8)*p(14) - p(61)
    p(89) = p(1)*p(39) - p(65)
    p(90) = p(3)*p(21) - p(66)
    p(91) = p(1)*p(41) - p(67)
    p(92) = p(2)*p(49) - p(91)
    p(93) = p(1)*p(43) - p(74)
    p(94) = p(8)*p(21) - p(93)
    p(95) = p(1)*p(46) - p(83)
    p(96) = p(7)*p(22) - p(95)
    p(97) = p(2)*p(45) - p(81)
    p(98) = p(2)*p(46) - p(85)
    p(99) = p(1)*p(50) - p(98)
    p(100) = p(6)*p(22) - p(87)
    p(101) = p(1)*p(49) - p(90)
    p(102) = p(2)*p(50) - p(100)
    p(103) = p(9)*p(13)
    p(104) = p(9)*p(12)
    p(105) = p(9)*p(14)
    p(106) = p(13)*p(15)
    p(107) = p(9)*p(20)
    p(108) = p(13)*p(10)
    p(109) = p(13)*p(11)
    p(110) = p(9)*p(16)
    p(111) = p(9)*p(10)
    p(112) = p(9)*p(11)
    p(113) = p(9)*p(17)
    p(114) = p(3)*p(24) - p(104)
    p(115) = p(7)*p(24) - p(105)
    p(116) = p(3)*p(31) - p(104) - p(115)
    p(117) = p(6)*p(39) - p(114)
    p(118) = p(3)*p(32) - p(105)
    p(119) = p(13)*p(21)
    p(120) = p(9)*p(18)
    p(121) = p(9)*p(19)
    p(122) = p(4)*p(24) - p(103) - p(103) - p(103)
    p(123) = p(1)*p(60) - p(107) - p(122)
    p(124) = p(13)*p(16)
    p(125) = p(14)*p(16) - p(119)
    p(126) = p(2)*p(56) - p(106) - p(125)
    p(127) = p(4)*p(32) - p(119)
    p(128) = p(1)*p(63) - p(107) - p(126)
    p(129) = p(13)*p(17)
    p(130) = p(6)*p(24) - p(108)
    p(131) = p(13)*p(12)
    p(132) = p(3)*p(45) - p(130)
    p(133) = p(13)*p(14)
    p(134) = p(9)*p(22)
    p(135) = p(8)*p(24) - p(109)
    p(136) = p(13)*p(18)
    p(137) = p(14)*p(18) - p(124)
    p(138) = p(6)*p(34) - p(109)
    p(139) = p(13)*p(19)
    p(140) = p(9)*p(9)
    p(141) = p(9)*p(15)
    p(142) = p(9)*p(21)
    p(143) = p(3)*p(41) - p(110)
    p(144) = p(4)*p(39) - p(112)
    p(145) = p(1)*p(68) - p(111) - p(144)
    p(146) = p(1)*p(69) - p(112)
    p(147) = p(1)*p(70) - p(113) - p(145)
    p(148) = p(11)*p(21) - p(110)
    p(149) = p(1)*p(72) - p(116) - p(115)
    p(150) = p(1)*p(73) - p(118)
    p(151) = p(3)*p(43) - p(120)
    p(152) = p(15)*p(18) - p(121) - p(151)
    p(153) = p(2)*p(69) - p(117)
    p(154) = p(1)*p(75) - p(121) - p(152)
    p(155) = p(7)*p(34) - p(120)
    p(156) = p(6)*p(41) - p(119)
    p(157) = p(1)*p(78) - p(126) - p(123)
    p(158) = p(2)*p(73) - p(119)
    p(159) = p(1)*p(80) - p(132) - p(130) - p(130)
    p(160) = p(13)*p(13)
    p(161) = p(1)*p(82) - p(132)
    p(162) = p(3)*p(46) - p(134)
    p(163) = p(2)*p(75) - p(127) - p(123)
    p(164) = p(1)*p(84) - p(134)
    p(165) = p(6)*p(43) - p(124)
    p(166) = p(1)*p(86) - p(138) - p(135)
    p(167) = p(8)*p(32) - p(124)
    p(168) = p(4)*p(45) - p(133)
    p(169) = p(2)*p(80) - p(131) - p(168)
    p(170) = p(13)*p(20)
    p(171) = p(2)*p(82) - p(133)
    p(172) = p(2)*p(83) - p(137) - p(135)
    p(173) = p(2)*p(84) - p(138)
    p(174) = p(6)*p(46) - p(136)
    p(175) = p(2)*p(86) - p(139) - p(169)
    p(176) = p(13)*p(22)
    p(177) = p(14)*p(22) - p(136)
    p(178) = p(3)*p(39) - p(141)
    p(179) = p(1)*p(89) - p(141) - p(178) - p(178)
    p(180) = p(3)*p(49) - p(142)
    p(181) = p(1)*p(91) - p(143)
    p(182) = p(2)*p(101) - p(181)
    p(183) = p(1)*p(93) - p(151)
    p(184) = p(8)*p(49) - p(183)
    p(185) = p(1)*p(95) - p(162)
    p(186) = p(21)*p(22) - p(185)
    p(187) = p(6)*p(45) - p(170)
    p(188) = p(1)*p(98) - p(172)
    p(189) = p(7)*p(50) - p(188)
    p(190) = p(2)*p(97) - p(170) - p(187) - p(187)
    p(191) = p(2)*p(98) - p(174)
    p(192) = p(1)*p(102) - p(191)
    p(193) = p(6)*p(50) - p(176)
    p(194) = p(1)*p(101) - p(180)
    p(195) = p(2)*p(102) - p(193)

    return
  end subroutine evpoly

end module bemsa