a=readline("Enter highest value of n to be computed up to")
a=as.numeric(a)

z=matrix(nrow=1, ncol=a)

for(n in seq(from=3, to=a, by=1))
{
  q=matrix(nrow=1,ncol=n)
  for(j in 1:n)
  {
    m=matrix(nrow=1,ncol=n)
    prob1r=0
    for (k in 0:(j-1))
    {
      prob1r=prob1r+(((dbinom(k,(j-1),(1/n))))*(1/(k+1)))
    }
    m[1,1]=1-((j/n)*prob1r)
    prob2r=0
    for(k in 0:(j-1))
    {
      probir=0
      for (i in 0:k)
      {
        probir=probir+((dbinom(i,k,m[1,1])))*(1/(1+i))
      }
      prob2r=prob2r+((dbinom(k,j-1,(1/n))))*probir
    }
    m[1,2]=1-((j/n)*(((n-1)/n)^(j-1))*prob2r)
    for(h in 3:n)
    {
      prob1=1
      for(y in 2:(h-1))
      {
        prob=0
        proby1=1
        for(k in 1:(y-1))
        {
          proby1=proby1*m[1,k]
        }
        for(k in 0:(j-1))
        {
          prob=prob+((dbinom(k,j-1,(1/n))))*((1-proby1)^k)
        }
        prob1=prob1*prob
      }
      prob2=0
      for(k in 0:(j-1))
      {
        probrej=1
        for (g in 1:(h-1))
        {
          probrej=probrej*(m[1,g])
        }
        probi=0
        for (i in 0:k)
        {
          probi=probi+((dbinom(i,k,probrej)))*(1/(1+i))
        }
        prob2=prob2+((dbinom(k,j-1,(1/n))))*probi
      }
      m[1,h]=1-((j/n)*(((n-1)/n)^(j-1))*prob1*prob2)
    }
    probno=1
    for(p in 1:n)
    {
      probno=probno*(m[1,p])
    }
    q[1,j]=(1-probno)*((((n*(n+1))/2)-(((n-j)*(n-j+1))/2))/j)
  }
  z[1,n]=((which.max(q))/n)*100
}

y=1/(z[1,3:a])
x=3:a
lm.r=lm(y~x)

plot(x,y)
abline(lm.r)