program basis_to_flat
implicit none
integer(8),parameter::par = 500
integer(8) i,j,k,l,n,m,nb,rk,binom,non
integer(8) s(par),b(par),c(par,par),d(par),e(par,par),f(par,par,par),nf(par,par,par),g(par,par),h(par,par),sr(par)
write(*,*)'台集合の大きさn=?'
read(*,*)n
write(*,*)'与える基底の個数=?'
read(*,*)nb
!基底族を定める。b(j)が与えた基底を表す。
write(*,*)'与える基底を記述する。'
write(*,*)'例　{1,2,4}　→　1101'
do i = 1, nb
write(*,*)i,'つ目'
read(*,'(b32.32)')b(i)
enddo
!階数を計算する。
rk = 0
do i = 1, n
    if (iand(b(1),2**(i-1)) == 2**(i-1)) rk = rk + 1
enddo
m = 2**n
s(1) = 0  !s(i)が部分集合すべてを表す。
do i= 2,m
  s(i) = s(i-1) + 1
end do
do i= 1,m  !部分集合と基底の共通部分を取る。c(i,j)がs(i)とb(j)の共通部分を表す。
 do j=1,nb
   c(i,j) = iand(s(i),b(j))
 enddo
enddo
do i=1,m
 do j=1,nb
    h(i,j)=0
 enddo
enddo
do i=1,m  !独立集合の階数を数える。
 do j=1,nb
   do l=1,n
     if (iand(c(i,j),2**(l-1)) == 2**(l-1)) h(i,j) = h(i,j)+1
   enddo
 enddo
enddo
do i=1,m
  sr(i)=0
enddo
do i=1,m  !階数を計算する。
 do j=1,nb
    if (sr(i) <= h(i,j)) sr(i) = h(i,j)
 enddo
enddo
non = -1        !-1を64bitで考える。
do i=1,m  !部分集合の補集合を取る。d(i)がs(i)の補集合を表す。
  d(i)= ieor(non,s(i))
end do
do i=1,m  !s(i)に補集合の元を加える。
  do k=1,n
   if (iand(d(i),2**(k-1)) == 2**(k-1)) e(i,k) = s(i)+2**(k-1)
   if (iand(d(i),2**(k-1)) /= 2**(k-1)) e(i,k) = 0
  enddo
enddo
do i = 1,m  !e(i,k)と基底との共通部分を取る。
  do j=1,nb
   do k=1,n
      f(i,j,k) = iand(e(i,k),b(j))
  enddo
 enddo
enddo
do i = 1,m
  do j=1,nb
   do k=1,n
      nf(i,j,k) = 0
  enddo
 enddo
enddo
do i = 1,m  !f(i,j,k)の個数を数える。
  do j=1,nb
   do k=1,n
    do l=1,n
      if (iand(f(i,j,k),2**(l-1)) == 2**(l-1)) nf(i,j,k) = nf(i,j,k) + 1
    enddo
  enddo
 enddo
enddo
do i = 1, m
 do k = 1, n
   g(i,k) = 0
 enddo
enddo
do i=1,m  !階数を計算する。
 do j=1,nb
  do k=1,n
    if (g(i,k) <= nf(i,j,k)) g(i,k) = nf(i,j,k)
  enddo
 enddo
enddo
do i=1,m
 do k=1,n
  if (g(i,k) == sr(i)) s(i) = -1
 enddo
enddo
do i = 0, rk
print'(A6,I3)','rank =',i
 do j = 1, m
  if ((sr(j) == i) .and. (s(j) /= -1)) print'(b8.8)',s(j)
 enddo
enddo
end program basis_to_flat
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function binom(n,k)
implicit none
integer(8) binom,n,k,i
if (n-k > k) k = n-k
binom = n
do i = n-1, k + 1,-1
binom = binom * i
enddo
do i = 1, n-k
binom = binom / i
enddo
end function binom