## Jiang Wenrong ##
## 2017.07.02 ##

###########################################################################
##    
##    Input: Multiplicity
##    Output: The value of d
##
###########################################################################

Valueofd:=proc(mu)
	local G,S,i,j,p,d3,d1,d2,k,l,mm;
     if mu=2 then 
     RETURN(fsolve(1-2*d^2-2*d*(1-d^2)^(1/2)-d,d,0..1));
     elif mu=3 then 
   RETURN(fsolve((-8*d^2-2*d+1)*(-d^2+1)^(1/2)-9*d-d^2+6*d^3,d,0..1));
     fi;

	G := Matrix(mu, mu); 
	S := Matrix(mu+1, mu+1);
	for i to mu do 
		for j to mu do 
			G[i, j] := factorial(i-1+j-1)/(factorial(i-1)*factorial(j-1)); 
		end do; 
	end do;
	for i to mu do 
		for j from 2 to mu+1 do 
			S[i, j] := factorial(i-1+j-1)/(factorial(i-1)*factorial(j-1));
		end do;
	end do;

	for i from 4 to mu+1 do
		for j from 2 to i-1 do 
			for l from 4 to mu+5-i do 
				for k to l-1 do 
					S[i-j+k-1, j+l-k-2] := G[k, l-k]*S[i-j, j]+S[i-j+k-1, j+l-k-2] 
				end do;
			end do;
		end do;
	end do;

#### t_ij=S`(i+1, j+2), c_ij = S(i+1, j+1) ####

	for mm from 4 to mu do 
		p := (-d^2+1)^((1/2)*mm)-add(S[mm-j+1, j+1]*d^j*(-d^2+1)^((mm-j)*(1/2)), j = 1 .. mm)-d*(add(S[i+1, 2], i = 1 .. mm-2)+add(add(S[m-j+1, j+2]*d^j*(-d^2+1)^((m-j)*(1/2)), j = 1 .. m), m = 1 .. mm-2)+1); 
		
		d1 := evalf(sqrt(1/(S[mm, 2]^2+1)));
		d2 := evalf(sqrt(1/(mm-1)));
             d3 := fsolve(p, d = 0 .. min(d1,d2));
		userinfo(3,Valueofd,print(`p(d)=`,p));
		userinfo(3,Valueofd,print(`mu=`,mm));
		userinfo(3,Valueofd,print(`d1=`,d1));
		userinfo(3,Valueofd,print(`d2=`,d2));
		userinfo(3,Valueofd,print(`d3=`,d3));
	end do;
   RETURN(d3);
    
end proc: