#OK to post
#Yuxuan Yang, Feb 14th, Assignment 7 

with(combinat):
with(plots):

#1
NiceValentine:=proc() local plot1, plot2:
plot1 := implicitplot3d((x^2 + 9/4*y^2 + z^2 - 1)^3 - x^2*z^3 - 9/80*y^2*z^3 = 0, x = -2 .. 2, y = -2 .. 2, z = -2 .. 2, style = surface, grid = [100, 100, 100]):
plot2 := textplot3d([-2, 0, -1.5, "I love Experimental Mathematics"], 'align' = 'left', axes = none):
plots[display]([plot1, plot2]):
end:

#2
DrawMandelbrot:=proc(R,K) local i,j,c,k,z,ab,S:
S:={}:
for i from 1 to 8*R+1 do
 for j from 1 to 8*R+1 do
  c:=-4+(i-1)/R+I*(-4+(j-1)/R):
  z:=0:
  ab:=0:
  for k from 1 by 1 while (k<K+1 and ab<10) do
    z:=z^2+c:
    ab:=evalf(abs(z)):
  od:
  if ab<10 then 
   S:=S union {[Re(c),Im(c)]}:
  fi:
od:
od:
pointplot(S, color = red):
end:

#3
RandSPn:=proc(n) local i,k,die1:
die1:=[seq(NuSPnk(n,i),i=1..n)]:
k:=LD(die1):
RandSPnk(n,k):
end:

#4
NuAdjacentMates:=proc(SP,n) local j,i,Nu:
Nu:=0:
for i from 1 to n-1 do 
 for j from 1 to nops(SP) do
  if {i,i+1} subset SP[j] then 
   Nu:=Nu+1:
  fi:
 od:
od:
Nu:
end:

EstimateAvAdjacentMates:=proc(n,K) local i, j,Nu:
Nu:=0:
for i from 1 to K do 
 Nu:=Nu+NuAdjacentMates(RandSPn(n),n):
od:
Nu/K:
end:
#EstimateAvAdjacentMates(100,1000):
#3.464 3.381 3.302
#
#The exact value is (n-1)*NuSPn(n-1)/NuSPn(n), and NuSPn(n)=B_n is just Bell number.
#Using an effecient algorithm of Bell number, it is not hard to give this value
#
#Here is the proof.
#Consider {i,i+1}, the probability of i and i+1 being in the same part is B_{n-1}/B{n}
#This is because if we glue i and i+1 together, we just need to partite n-1 elements.
#By linearity of expectation, the AvAdjacentMates is (n-1) times this probability. 




















#from C7.txt
#LD(L): Inputs a list of positive integers L (of n:=nops(L) members)
#outputs an integer i from 1 to n with the prob. of i being
#proportional to L[i]
#For example LD([1,2,3]) should output 1 with prob. 1/6
#output 2 with prob. 1/3
#output 3 with prob. 3/6=1/2
LD:=proc(L) local n,i,su,r:
n:=nops(L):
r:=rand(1..convert(L,`+`))():
su:=0:

for i from 1 to n do
 su:=su+L[i]:
  if r<=su then
    RETURN(i):
  fi:
od:
end:

#RandSPnk(n,k): a random set partition of {1,..,n} with k parts
RandSPnk:=proc(n,k) local c,SP1,k1:
if n=1 then
 if k=1 then
   RETURN({{1}}):
  else
   RETURN(FAIL):
 fi:
fi:
c:=LD([NuSPnk(n-1,k-1), k*NuSPnk(n-1,k)]):

if c=1 then
 SP1:=RandSPnk(n-1,k-1):
 RETURN(SP1 union {{n}}):
 else
 SP1:=RandSPnk(n-1,k):
  k1:=rand(1..k)():
  RETURN( { op(1..k1-1,SP1), SP1[k1] union {n}, op(k1+1..k,SP1)}):
  
fi:
end:

#NuSPnk(n,k): The number of set-partitions of {1,...,n}
#into exactly k parts (aka Strirling numbers of the second kind)
NuSPnk:=proc(n,k) option remember:

if n=1 then
  if k=1 then
   RETURN(1):
  else
   RETURN(0):
 fi:
fi:
NuSPnk(n-1,k-1)+ k*NuSPnk(n-1,k):
end:

NuSPn:=proc(n) local k:
option remember:
add(NuSPnk(n,k),k=1..n):
end: