read(`C10.txt`): #C10 should be in the same folder
Randomize():  

#################
# Simulations
###################
# SimulateNAND(10,5, 1000) = 696.576
#                             87072
#                             -----
#                              125
# 696.576, so about 68% of all true/false vectors output true.

# SimulateNAND(10,10, 1000) = 661.712
# 82714/125
# about 65% of all T/F vectors output true


#SimulateNAND(10,20, 1000) = 645.7160000
#
#                             161429
#                             ------
#                              250  


# SimulateNAND(10,50, 1000) = 635.1870000
#
#                             635187
#                             ------
#                              1000 

######################################################
# SimulateMonotone(10,5, 1/2,1000) = 515.9360000
# [SimulateMonotone(10,10, 1/2, 1000), SimulateMonotone(10,20, 1/2, 1000), SimulateMonotone(10,50, 1/2, 1000)]
# [65197/125, 25347/50, 529133/1000] = [521.5760000, 506.9400000, 529.1330000]
#
# [SimulateMonotone(10,5, 3/4,1000), SimulateMonotone(10,10, 3/4, 1000), SimulateMonotone(10,20, 3/4, 1000), SimulateMonotone(10,50, 3/4, 1000)]
# [44118/125, 81723/250, 142367/500, 22741/100] = [352.9440000, 326.8920000, 284.7340000, 227.4100000]
#
# [SimulateMonotone(10,5, 1,1000), SimulateMonotone(10,10, 1, 1000), SimulateMonotone(10,20, 1, 1000), SimulateMonotone(10,50,1, 1000)]
# [5382/25, 20163/125, 13903/125, 10651/200] = [215.2800000, 161.3040000, 111.2240000, 53.25500000]
###############################################################################################################









SimulateNAND:=proc(n,k,K) local P,truetot,i,V,v:
 truetot:=0:
 V:=TF(n):
 for i from 1 to K do: 
  P:=RSLP(n,k):
  for v in V do: # doing this will bypass the construction of the set in EvalSP
   if EvalSP1(P,v) then: truetot++: fi:
  od:
 od:
 truetot/K:
end:

SimulateMonotone:=proc(n,k,p,K) local P,truetot,i,V,v:
 truetot:=0:
 V:=TF(n):
 for i from 1 to K do: 
  P:=Monotone(n,k,p):
  for v in V do: # doing this will bypass the construction of the set in EvalSP
   if EvalMSP1(P,v) then: truetot++: fi:
  od:
 od:
 truetot/K:
end:

# Outputs one of AND with prob p, OR with prob 1-p
LDMon:=proc(p):
if rand(1..denom(p))()<=numer(p) then: return(AND) fi:

OR:
end:

# Modified RSLP. This version doesn't append to list, which is slow. Shouldn't matter for small n+k though
Monotone:=proc(n,k,p) local P,i,i1,j1:
 P:=[seq(i,i=1..n),seq(-1,i=n+1..k+n)]:

 for i from n+1 to k+n do
  i1:=rand(1..i-1)():
  j1:=rand(1..i-1)():
  P[i]:=[i1,j1,LDMon(p)]:
 od:
P:

end:

# Evaluates Monotone RSLP for a given vector
EvalMSP1:=proc(P,IN) local n,P1,i,i1,j1,oper:
P1:=IN:
n:=nops(IN):
for i from n+1 to nops(P) do
 i1:=P[i][1]:
 j1:=P[i][2]:
 oper:=P[i][3]:
 P1:=[op(P1), oper(P1[i1],P1[j1]) ]:
od:

P1[-1]:
end: