program gen_base
implicit none
integer(8),parameter::par = 800
integer(8) n,t,tmax,nb,nbmin,nbmax,rkmin,rkmax,i,j,k,l,h,d,d1,d2,u
integer(8) s(par),ns(par)
integer(8) factorial,binom
logical(1) z,t1,t2,loopless,c
real x
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
!!!!!!!!!!!!!!!!!!!!!
z = .false.
do t = 1, tmax
 call random_number(x)   !基底の個数を定める。
 nb = int(x*binom(n,n/2)) + 1
 if (nb > par) 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     !s(i)\s(j)の元を取る。
   c = .false.
    do l = 1, n
     if (iand(d2,2**(l-1)) /= 2**(l-1)) cycle   !s(j)\s(j)の元を取る。
     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
    if (z .eqv. .true.) print'(A6,I32)','階数 = ',ns(1)
print'(A6,I32)','試行回数 = ',t
if (z .eqv. .true.) print'(A6)','成功'
if (z .eqv. .false.) print'(A6)','失敗'
end program gen_base
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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