#OK to post homework
#George Spahn, 2/14/2021, Assignment 7

with(plots):
# plots the Mandelbrot set by picking all those complex c 
# (of the form x+I*y, with x,y in [-4,4] with resolution 1/R, and deciding
# whether to accept it if iterating z goes to z^2+c (starting at z=0) after
# K iterations stay bounded
DrawMandelbrot:=proc(R,K) local i,j,c,m,S:
S:=[]:

m:=proc(c1,K1) local d,k:
d:=0:
for k from 1 to K1 do
  d:= evalf(d*d+c1):
  if evalf(abs(d)) > evalf(4) then
   RETURN(false):
  fi:
od:
true: 
end:

for i from 1 to 8*R do
for j from 1 to 8*R do
 c:= -4+i/R + I*(-4+j/R):
 if m(c,K) then
  S:=[op(S),c]:
 fi: 
od:
od:
complexplot(S,style=point,axes=none):
end:



#RandSPn(n)  Inputs a positive integer n and outputs a (uniformly) at random
# set partition of {1,...,n}
RandSPn:=proc(n) local Die,Die1:

Die1:=[seq(NuSPnk(n,k),k=1..n)]:
Die:=LD(Die1):
RandSPnk(n,Die):
end:


#NuAdjacentMates(SP,n) Inputs a set partition, SP, of {1, ...n}, and outputs 
#the number of i between 1 and n such that i and i+1 are part-mates

NuAdjacentMates:= proc(SP,n) local i,j,c:
c:=0:
for i from 1 to n-1 do
 for j from 1 to nops(SP) do
  if i in SP[j] and (i+1) in SP[j] then
   c:=c+1:
  fi:
 od:
od:
c:
end:


#EstimateAvAdjacentMates(n,K) uses the above procedure K times, 
# to estimate the average of the number of adjacent mates in a
# random set partition of {1,...n}.


EstimateAvAdjacentMates:=proc(n,K) local T,i:
T:=0:
for i from 1 to K do
 T:=T+NuAdjacentMates(RandSPn(n),n):
od:
evalf(T/K):
end:

#EstimateAvAdjacentMates(100, 1000) gave me ~3.4

#



#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:

#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:

#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: