program basis_to_flat
implicit none
integer(8) i,j,k,l,n,m,nb,rk,binom,rkmin,rkmax,tmax,y,t,non
integer(8),parameter::param = 500
integer(8) s(param),b(param),c(param,param),d(param),e(param,param),f(param,param,param),nf(param,param,param),g(param,param),h(param,param),sr(param)
logical(1) loopless,t1,t2,z
!試行回数
write(*,*)'マトロイドの基底とフラットを無作為に生成します。'
write(*,*)'試行回数の上限?'
write(*,*)'注:31以下に設定してください。'
read(*,*)t
tmax = 2**t - 2
!台集合の大きさ。
write(*,*)'台集合の元の個数?'
read(*,*)n
!ループの存在を許さないか
write(*,*)'ループなしにする？ T or F'
read(*,*)loopless
!階数の最小値を設けるか
write(*,*)'階数の最小値を設ける? T or F'
read(*,*)t1
if(t1 .eqv. .true.) then
write(*,*)'最小値?'
read(*,*)rkmin
endif
!階数の最大値を設けるか
write(*,*)'階数の最大値を設ける? T or F'
read(*,*)t2
if (t2 .eqv. .true.) then
write(*,*)'最大値?'
read(*,*)rkmax
endif
!!!!!!!!!!!!!!!!
m = 2**n
y = 0
if (m > param) y = y + 1
if (m > param) goto 13
!基底族を定める。b(j)が与えた基底を表す。
call genbase_unlim_sub(rk,nb,b,n,loopless,t1,t2,rkmin,rkmax,tmax,z)
if (z .eqv. .false.) goto 14
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
13 if (y == 1) print*,'台集合が大きすぎます。'
14 continue
end program basis_to_flat
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine genbase_unlim_sub(rk,nb,s,n,loopless,t1,t2,rkmin,rkmax,tmax,z)
implicit none
integer(8) n,m,t,rk,tmax,nb,nbmin,nbmax,rkmin,rkmax,i,j,k,l,h,d,d1,d2,u
integer(8),parameter:: param = 500
integer(8) s(param),ns(param)
integer(8) factorial,binom
logical(1) z,t1,t2,loopless,c
real x
!!!!! 与える条件
!台集合の元の個数
!n =
!!!!!!!!!!!!!!!!!!!!!
z = .false.
m     = 2*n        !部分集合の個数
do t = 1, tmax
 call random_number(x)   !基底の個数を定める。
 nb = int(x*binom(n,n/2)) + 1
 if (nb > param) cycle
 do i = 1, nb   !部分集合を無作為に与える。
 call random_number(x)
 s(i) = int(x*(2**n))
 enddo
 do i = 1, nb    !同じ集合がある場合はやり直し。
  do j = 1, nb
    if (i == j) cycle
    if (s(i) == s(j)) goto 10
  enddo
 enddo
 if (loopless .eqv. .true.) then !ループがあるかを判定する。
   u = s(1)
   do i = 2, nb
      u = ior(u,s(i))
   enddo
   if (u /= 2**n-1) goto 10
 endif
 do i = 1, nb    !部分集合の大きさを測る。
      ns(i) = 0
 enddo
 do i = 1, nb
  do j = 1, n
         if (iand(s(i),2**(j-1)) == 2**(j-1)) ns(i) = ns(i) + 1
  enddo
 enddo
 do i = 1, nb   !元の個数が一様でないならばやり直し。
  do j = 1, nb
     if (ns(i) /= ns(j)) goto 10
  enddo
 enddo
 if ((t1 .eqv. .true.) .and. (ns(1) < rkmin)) goto 10     !基底の最小値より小さいならばやり直し。
 if ((t2 .eqv. .true.) .and. (ns(1) > rkmax)) goto 10     !基底の最大値より大きいならばやり直し。
!基底族の公理を満足するか判定する。
 do i = 1, nb
  do j = 1, nb
  if (i == j) cycle
  d = ieor(s(i),s(j))        !s(i)とs(j)の対称差を取る。
  d1 = iand(d,s(i))     !s(i)からs(j)を引く。
  d2 = iand(d,s(j))     !s(j)からs(i)を引く。
   do k = 1, n
   if (iand(d1,2**(k-1)) /= 2**(k-1)) cycle
   c = .false.
    do l = 1, n
     if (iand(d2,2**(l-1)) /= 2**(l-1)) cycle
     do h = 1, nb
            if ((s(i) - 2**(k-1)) + 2**(l-1) == s(h)) c = .true.
     enddo
    enddo
        if (c .eqv. .false.) goto 10
   enddo
  enddo
 enddo
    z = .true.
    exit
10 enddo
    do i = 1, nb
      print'(b32.32)',s(i)
    enddo
    rk = ns(1)
    if (z .eqv. .true.) print'(A6,I32)','階数 = ',rk
print'(A6,I32)','試行回数 = ',t
if (z .eqv. .true.) print'(A6)','成功'
if (z .eqv. .false.) print'(A6)','失敗'
end subroutine genbase_unlim_sub

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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