#(n,dn-1) core partitions
#Anthony Zaleski, 2016-17

print(`(n,dn-1)-CORE PARTITIONS`):
print(`-----`):
print(`Execute Help(); to begin.`):
print(`-----`):
print(`IMPORTANT: do not assign to the variable names`):
print(`N1, N2, NRatio,`):
print(`which are used globally in this package.`):

################
#HELP PROCEDURES
################



Help:=proc(): 
if nargs=0 then
	print(`qrec.txt contains the following primary procedures.`):
	print(`Type Help(<procedure name>); for help on a procedure.`):
	print(`Type Help1(); to list secondary procedures.`):
	print(`Ndn, NRatioLim, Gdn`):
	print(`mdnK, Mdnk, Mcdnk, Msdnk`):
else 
	LookupProc(args[1]):
fi:
end:


Help1:=proc(): 
if nargs=0 then
	print(`qrec.txt contains the following secondary procedures.`):
	print(`Type Help(<procedure name>); for help on a procedure.`):
	print(`Type Help(); to list primary procedures.`):
	print(`GdnMoms1,GdnMoms2`):
	print(`GuessRec, Rec2Seq, qGuessRec, qRec2Seq`):
	print(`GuessRecPol, GuessMoms, GuessqRecMoms`):
	print(`RecPol2Seq `):
	print(`AnalyzePolSeq `):
	print(`GuessRec1, qGuessRec1`):
	print(`Moms, Straight2Central`):
	print(`GuessPoly, GuessRecPol1 `):
	
else 
	LookupProc(args[1]):
fi:
end:


LookupProc:=proc(a) local EG:
EG:=`For example, try`:
if a=Help then
	print(`Help(): The primary help procedure.`):
elif a=Help1 then
	print(`Help1(): The secondary help procedure.`):
elif a=Ndn then
	print(`Ndn(d,n): The number of (n,dn-1)-core partitions with distinct parts.`):
	print(`d may be numeric or symbolic; n must be numeric.`):
	print(EG):
	print(`Ndn(d,6);`):
elif a=Gdn then
	print(`Gdn(q,d,n): The g.f. (according to size) of (s,dn-1)-core partitions with distinct parts.`):
	print(EG):
	print(`Gdn(q,3,7);`):
elif a=GuessRec1 then
	print(`GuessRec1(L,r): inputs a sequence of numbers, L, and a pos.`):
	print(`integer r, and tries to guess a linear recurrence equation`):
	print(`of order r such that L(n)=c1*L(n-1)+...+cr*L(n-r)=0`):
	print(`for all n within the range.`):
	print(`The  output consists of [initial conditions, recursion].`):
	print(`For example, the Fibonacci sequence is [[0,1],[1,1]].`):
	print(EG):
	print(`GuessRec1([1,1,2,3,5,8,13,21,34],2);`):
elif a=GuessRec then
	print(`GuessRec(L): inputs a sequence of numbers, L,`):
	print(`and tries to guess a linear recurrence equation`):
	print(`of order r such that L(n)=c1*L(n-1)+...+cr*L(n-r)`):
	print(`for all n within the range.`):
	print(`The  output consists of [initial conditions, recursion].`):
	print(`For example, the Fibonacci sequence is [[0,1],[1,1]].`):
	print(EG):
	print(`GuessRec([0,1,1,2,3,5,8,13,21,34]);`):
elif a=Rec2Seq then 
	print(`Rec2Seq(R,N): Computes the first N terms of the recursion`):
	print(`specified by R=[initial conditions, recursion]. For example,`):
	print(`R=[[0,1],[1,1]] gives the Fibonacci sequence.`):
	print(EG):
	print(`Rec2Seq([[0,1],[1,1]],15);`):
elif a=qGuessRec1 then
	print(`qGuessRec1(L,q,n,k,d): inputs list L of polys in q.`):
	print(`Outputs [ini,rec] s.t. ini=L[1..k] and`):
	print(`rec is a list of k+1 polys in q^n of degree <=d`):
	print(`satisfying rec[1](n)*L[n-k]+...+rec[k+1](n)*L[n]=0`):
	print(`for all n.  If not possible, returns FAIL.`):
	print(EG):
	print(`qGuessRec1([seq(q^i,i=1..50)],q,n,1,0);`):
elif a=qGuessRec then
	print(`qGuessRec(L,q,n,dmax:=4): inputs list L of polys in q.`):
	print(`Outputs [ini,rec] s.t. ini=L[1..k] and`):
	print(`rec is a list of k+1 polys in q^n of degree <=dmax`):
	print(`satisfying rec[1](n)*L[n-k]+...+rec[k+1](n)*L[n]=0`):
	print(`for all n.  If not possible, returns FAIL.`):
	print(EG):
	print(`qGuessRec([seq(q^i,i=1..50)],q,n);`): 
elif a=qRec2Seq then
	print(`qRec2Seq(Rec,q,n,N): inputs Rec=[ini,rec] specifying`):
	print(`a q-recurrence and outputs the first N terms.`):
	print(EG):
	print(`qRec2Seq([[q],[q,-1]],q,n,5);`):
elif a=GuessRecPol1 then
	print(`GuessRecPol1(L,R,n,d): Inputs list L, and recursion`):
	print(`R=[ini,rec] of order k specifying sequence F(n). outputs [A[1],...,A[k]], where A[i] is`):
	print(`a poly. in n of degree <=d s.t. L[n]=F(n)A[1](n)+...+F(n+k-1)A[k](n).`):
	print(EG):
	print(`GuessRecPol1([1,1,2,3,5,8,13],[[1,1],[1,1]],n,1);`):
elif a=GuessRecPol then
	print(`GuessRecPol(L,R,n,n0:=1,mind:=0): Inputs a list L=[a(n0),a(n0+1),...] specifying sequence a(n), and a recursion`):
	print(`R=[ini,rec] of order k.  outputs [A[1],...,A[k]], where A[i] is`):
	print(`a poly. in n of degree >=mind, s.t. a(n)=F(n)A[1](n)+...+F(n+k-1)A[k](n).   `):
	print(EG):
	print(`GuessRecPol([1,1,2,3,5,8,13],[[1,1],[1,1]],n);`):
elif a=Moms then
	print(`Moms(L,t,K): inputs a list L of poly.s in t`):
	print(`(or a single poly). Outputs list M s.t. M[n][k] is the`):
	print(`kth moment of L[n] (or s.t. M[k] is the kth mom. of L.`):
	print(`DOES NOT divide by p(1).`):
	print(EG):
	print(`Moms([1+t^2,1+t^3],t,3);`): 
elif a=GuessMoms then
	print(`GuessMoms(L,q,K,n,rec:=0): inputs list L of polys in q.`):
	print(`Outputs pair R,M, where R specifies an order d recursion satisfied`):
	print(`by the sequence {L[n](1)}, and M[k]=[p1,p2,...pd] is`):
	print(`a list of d polys in n s.t. the kth pre-moment of L[n] is`):
	print(`L[n](1)*p1+...+L[n+d-1](1)*pd.`):
	print(`If desired, there is an optional input to specify the recursion manually.`):
	print(EG):
	print(`GuessMoms([seq((1+q)^n,n=0..20)],q,3,n);`):
elif a=GuessqRecMoms then
	print(`GuessqRecMoms(Rec,q,n,K,N:=50): q-recursion Rec, and symbols q,n specifying a poly. seq. {f[n](q)}.`):
	print(`Using the first N f_n(q), guesses pair R,M, where R specifies an order d recursion satisfied`):
	print(`by the sequence {f_n(1)}, and M[k]=[p1,p2,...pd] is`):
	print(`a list of d polys in n s.t. the kth moment of f_n[n] is`):
	print(`fL[n](1)*p1+...+f[n+d-1](1)*pd.`):
	print(EG):
	print(`GuessqRecMoms([[1,1+q,q^2+q+1,2*q^3+q^2+q+1],[-q^n/q,-q^n/q,0, -1, 1]],q,n,3,50);`):
elif a=RecPol2Seq then
	print(`RecPol2Seq(R,A,N,n): given order d recursion R=[ini,rec] specifying seq. F(n) and list of`):
	print(`d polys A in n, computes the first N terms of the squence L s.t.`):
	print(`L[n]=F(n)A[1](n)+...+F(n+k-1)A[k](n)`):
	print(EG):
	print(`RecPol2Seq([[0,1],[1,1]],[1,0],5,n);`):
elif a=AnalyzePolSeq then
	print(`AnalyzePolSeq(L,q,K,n,a,r,verbose:=1,rec:=0): Given the sequence L=[p1,p2,...] of`):
	print(`polynomials in q, tries to guess closed forms for up to the Kth moment of pn in`):
	print(`terms of n. n,a,r are symbols. If verbose=1, prints explanation; if verbose=0,`):
	print(`prints only the data:`): 
	print(`R,M,Mas,r1,Ml`):
	print(`where R is the recursion satisfied by an=pn(1), M is the list`):
	print(`[mu, central moments...], Mas is an asymptotic version of M in terms of r,`):
	print(`the limit of a(n+1)/a(n), r1 is the numeric value of r, and Ml is the list`):
	print(`[limit of mu/sigma,limit of standardized central moments...].`):
	print(`optional input rec allows manually forcing the recursion for an.`):
	print(EG):
	print(`AnalyzePolSeq([seq((1+q)^n,n=0..20)],q,4,n,a,r)`):
elif a=GuessPoly then
	print(`GuessPoly(X,Y,x,MaxDeg:=10,MinDeg:=0): looks for polynomial in x of degree`):
	print(`between MinDeg and MaxDeg interpolating the points with x-values X and y-values Y.`):
	print(EG):
	print(`GuessPoly([seq(i,i=1..10)],[seq(i^2,i=1..10)],x)`):
elif a=GdnMoms1 then 
	print(`GdnMoms1(d,n,K)`):
	print(`Inputs symbolic or numeric d, numeric n; finds polynomial expressions for`):
	print(`first K straight moments of Gdn (BEFORE dividing by Ndn).`):
	print(EG):
	print(`GdnMoms1(d,3,2);`):
elif a=GdnMoms2 then 
	print(`GdnMoms2(d,n,K)`):
	print(`Inputs symbolic or numeric d, symbolic n; finds the`):
	print(`first K straight moments of Gdn (BEFORE dividing by Ndn).`):
	print(`outputs a list of pairs [ai,bi] such that the ith moment is`):
	print(`ai*Ndn(d,n)+bi*Ndn(d,n+1).`):
	print(EG):
	print(`GdnMoms2(d,n,1);`):
elif a=Straight2Central then 
	print(`Straight2Central(L): converts list of straight moments (starting with first)`):
	print(`to central ones.`):
	print(EG):
	print(`Straight2Central(GdnStraightMoms1(d,3,2));`):
elif a=NRatioLim then
	print(`NRatioLim(d): The limit of Ndn(d,n+1)/Ndn(d,n) as n tends to infinity.`):
	print(`d may be numeric or symbolic.`):
	print(`For example, to get the Golden Ratio, try`):
	print(`NRatioLim(1);`):
elif a=mdnK then 
	print(`mdnK(d,n,K): the first K pre-moments of X_{d,n}, the random variable`):
	print(`size of an (n,dn-1)-core partition with distinct parts.`):
	print(`d and n may independently be numeric or symbolic.`):
	print(`Output may be include the symbols N1:=Ndn(d,n) and N2:=Ndn(d,n+1).`):
	print(EG):
	print(`mdnK(3,n,4);`):
elif a=MdnK then 
	print(`MdnK(d,n,K): the first K straight moments of X_{d,n}, the random variable`):
	print(`size of an (n,dn-1)-core partition with distinct parts.`):
	print(`d and n may independently be numeric or symbolic.`):
	print(`Output may include the symbol NRatio:=Ndn(d,n)/Ndn(d,n+1).`):
	print(EG):
	print(`MdnK(d,3,4);`):
elif a=McdnK then 
	print(`McdnK(d,n,K): the first K central moments of X_{d,n}, the random variable`):
	print(`size of an (n,dn-1)-core partition with distinct parts.`):
	print(`d and n may independently be numeric or symbolic.`):
	print(`Output may include the symbol NRatio:=Ndn(d,n)/Ndn(d,n+1).`):
	print(EG):
	print(`McdnK(d,3,2);`):
elif a=MsdnK then 
	print(`McdnK(d,n,K): the first K standardized moments of X_{d,n}, the random variable`):
	print(`size of an (n,dn-1)-core partition with distinct parts.`):
	print(`d and n may independently be numeric or symbolic.`):
	print(`Output may include the symbol NRatio:=Ndn(d,n)/Ndn(d,n+1).`):
	print(EG):
	print(`MsdnK(d,3,3);`):
else
	print(`No such procedure exists.`):
fi:
end:



#################
#################

Ndn:=proc(d,s) option remember:
if s=1 then return 1: fi:
if s=2 then return d: fi:
expand(Ndn(d,s-1)+d*Ndn(d,s-2)):
end:


Gdn1:=proc(q,t,d,s,a) local i,j: option remember:
if a=s-1 then return add(q^(s*i^2/2-s*i/2+a*i)*t^i,i=0..d-1): fi:
if a>=s then return 1: fi:
expand( Gdn1(q,t,d,s,a+1)+add(q^(s*i^2/2-s*i/2+a*i)*t^i,i=1..d)*Gdn1(q,t,d,s,a+2) ):
end:


Gdn:=proc(q,d,s) local g,t: #option remember:
g:=Gdn1(q,t,d,s,1):
expand(add(q^(k*(1-k)/2)*coeff(g,t,k),k=0..degree(g,t))):
end:

#################
#################
 
GuessRec1:=proc(L,r) local c,i,var,eq,n:

if nops(L)<=2*r+3 then
 RETURN(FAIL):
fi:


var:={seq(c[i],i=1..r)}:


eq:={seq( L[n]=add(c[i]*L[n-i],i=1..r) , n=r+1..nops(L))}:

var:=solve(eq,var):

if var=NULL then
  RETURN(FAIL):
fi:

[L[1..r], [seq(subs(var, c[i]),i=1..r)]   ] :

end:


GuessRec:=proc(L) local r,ans:

for r from 1 to (nops(L)-3)/2 do
 ans:=GuessRec1(L,r):
 if ans<>FAIL then
  RETURN(ans):
 fi:
od:

FAIL:
end:


Rec2Seq:=proc(R,N) local ini,rec,Se,k,n,nt:

ini:=R[1]:
rec:=R[2]:
k:=nops(ini):

if nops(rec)<>k then
 RETURN(FAIL):
fi:

Se:=ini:

for n from k+1 to N do

#nt:=c1*Se[n]+c2*Se[n-1]+...+ ck*se[n-k+1]
nt:=add(rec[i]*Se[n-i],i=1..k) :
nt:=normal(nt):
nt:=expand(nt):
Se:=[op(Se),nt]:

od:


Se:

end:



qGuessRec1:=proc(L,q,n,k,d) local a,c,i,j,p,rec,var,var1,N,eq:
N:=k*(d+1)+4:

if nops(L)<N then return FAIL: fi:

var:={seq(seq(c[i][j],j=0..d),i=0..k)}:
rec:=[seq(add(q^(n*j)*c[i][j],j=0..d),i=0..k)]:
eq:={seq(add(subs(n=j,rec[-i])*L[j+1-i],i=1..k+1),j=k+1..nops(L))}:

N:=min(nops(eq),N):

var1:=solve(subs(q=3,eq[1..N]),var):
if var1=NULL then return FAIL: fi:

var1:=solve(eq[1..N],var):
var:={seq(a=1,a in var)}:

rec:=subs(var,subs(var1,rec)):
if var1=NULL or rec[-1]=0  then return FAIL: fi:

for i from N+1 to nops(eq) do
  if expand(subs(var1,eq[i]))<>0 then return FAIL: fi:
od:

[L[1..k],rec]:
end:


qGuessRec:=proc(L,q,n,dmax:=4) local k,Rec,d:
for k from 1 to nops(L) do
  for d from 0 to dmax do
    Rec:=qGuessRec1(L,q,n,k,d):
    if Rec<>FAIL then return Rec: fi:
  od:
od:
FAIL:
end:


qRec2Seq:=proc(Rec,q,n,N) local L,rec,i,k:
L:=Rec[1]:
rec:=Rec[2]:
k:=nops(rec)-1:
for i from nops(Rec[1])+1 to N do
  L:=[op(L),
    expand(-add(subs(n=i,rec[-j-1])*L[-j],j=1..k)/subs(n=i,rec[-1])) ]:
od:
L:
end:


qGuessRecOpt:=proc(L,q,n,R:=4) local k,Rec,d,r:

for r from 1 to R do
  for k from 1 to r do
    Rec:=qGuessRec1(L,q,n,k,r-k):
    if Rec<>FAIL then return Rec: fi:
  od:
od:

FAIL:
end:




GuessRecPol1:=proc(L,R,n,d,n0:=1) local K,k,A,a,f,i,eq,var,F:
K:=nops(R[1]):
if nops(L)<K*d+4 then
  return FAIL:
fi:

for k from 1 to K do
  A[k]:=add(a[k][i]*n^i,i=0..d):
od:
var:={seq(seq(a[k][i],i=0..d),k=1..K)}:
F:=Rec2Seq(R,nops(L)+K-1+(n0-1)):
eq:={seq(L[i]=add(subs(n=i+(n0-1),A[k])*F[i+k-1+(n0-1)],k=1..K),i=1..K*d+4)}:
var:=solve(eq,var):
for k from 1 to K do
  A[k]:=subs(var,A[k]):
od:

if var=NULL then return FAIL: fi:

for i from K*d+5 to nops(L) do
  if  L[i]<>expand(normal(add(subs(n=i+(n0-1),A[k])*F[i+k-1+(n0-1)],k=1..K))) then return FAIL: fi:
od:

[seq(A[k],k=1..K)]:
end:


GuessRecPol:=proc(L,R,n,n0:=1,mind:=0) local d,A:
for d from mind to (nops(L)-4)/nops(R[1]) do
  A:=GuessRecPol1(L,R,n,d,n0):
  if A<>FAIL then return A: fi:
od:
FAIL:
end:


Moms:=proc(L,t,K) local tD,N,M1,f,k,M,n:
if not type(L,`list`) then
	return op(Moms([L],t,K)):
fi:

N:=nops(L):
M:=[]:
for n from 1 to N do
	M1:=[]:
	f:=L[n]:
	for k from 1 to K do
		f:=t*diff(f,t):
		M1:=[op(M1),subs(t=1,f)]:
	od:
	M:=[op(M),M1]:
od:
M:
end:


Straight2Central:=proc(L) local k,j:
[seq(normal(add(L[k-j]*L[1]^j*(-1)^j*binomial(k,j),j=0..k-1)+(-L[1])^k),k=1..nops(L))]:
end:

Straight2Standard:=proc(L) local M:
M:=Straight2Central(L):
[seq()]:
end:


GuessMoms:=proc(L,t,K,n,rec:=0) local N,M0,M,R,k,nn:
N:=nops(L):
M0:=[seq(subs(t=1,L[nn]),nn=1..N)]:
if rec<>0 then 
	R:=rec:
else
	R:=GuessRec(M0):
fi:
M:=Moms(L,t,K):
R,[seq(GuessRecPol([seq(M[nn][k],nn=1..N)],R,n),k=1..K)]:
end:



GuessqRecMoms:=proc(Rec,q,n,K,N:=50) 
GuessMoms(qRec2Seq(Rec,q,n,N),q,K,n):
end:


RecPol2Seq:=proc(R,A,N,n) local L,n1,i:
L:=Rec2Seq(R,N+nops(A)):
seq(expand(normal(subs(n=n1,add(L[n1+i-1]*A[i],i=1..nops(A))))),n1=1..N):
end:


RecRoot:=proc(R) local rec,r:
rec:=R[2]:
max(solve(r^nops(rec)=add(rec[k]*r^(nops(rec)-k),k=1..nops(rec)),r)):
end:

NRatioLim:=proc(d)
1/2+(1/2)*sqrt(4*d+1):
end:


Straight2Central:=proc(L) local k,j:
[seq(expand(add(L[k-j]*L[1]^j*(-1)^j*binomial(k,j),j=0..k-1)+(-L[1])^k),k=1..nops(L))]:
end:

#################
#GUESS POLYNOMIAL
#################

GuessPoly1:=proc(X,Y,x,Deg) local eq,var,a,P,i:
P:=add(x^i*a[i],i=0..Deg):
var:={seq(a[i],i=0..Deg)}:
eq:={seq(subs(x=X[i],P)=Y[i], i=1..Deg+1)}:

var:=solve(eq,var):
P:=subs(var,P):

if var=NULL or CheckAnsatz(X,Y,P,x)=false then 
	return FAIL: 
fi:
P:
end:


GuessPoly:=proc(X,Y,x,MaxDeg:=10,MinDeg:=0) local d,P:
if nops(X)<>nops(Y) then
	print(`ERROR in GuessPoly: X,Y must be lists of the same length.`):
	return(FAIL):
fi:

for d from MinDeg to MaxDeg do
	if nops(X)<d+4 then
		#print(`ERROR in GuessPoly: the lists X and Y are too short`):
		#print(`to make a reliable guess.  Try more data points.`):
		return FAIL:
	fi:

	P:=GuessPoly1(X,Y,x,d):
	if P<>FAIL then
		return P:
	fi:
od:
FAIL:
end:



AnalyzePolSeq:=proc(L,q,K,n,a,r,verbose:=1,rec:=0) local RM,R,M1,M,Mas,Ml,r1,i,j,p:
RM:=GuessMoms(L,q,K,n,rec):
R:=RM[1]: M:=RM[2]:
if R=FAIL or M[1]=FAIL then return FAIL: fi:
for i from 2 to nops(M) do
	if M[i]=FAIL then 
		M:=M[1..i-1]:
		break:
	fi:
od:

M1:=M:

M:=[seq(expand(add( M[j][i]*a[n+i-1],i=1..nops(R[2]) )/a[n]),j=1..nops(M))]:
M:=[M[1],op(2..nops(M),Straight2Central(M))]:

Mas:=[seq(expand(add( M1[j][i]*r^(i-1),i=1..nops(R[2]) )),j=1..nops(M))]:
Mas:=[Mas[1],op(2..nops(Mas),Straight2Central(Mas))]:

r1:=RecRoot(R):

Ml:=[seq(limit(subs(r=r1,Mas[i]/sqrt(Mas[2])^i),n=infinity),i=1..nops(Mas))]:

if verbose=0 then return R,M,Mas,r1,Ml: fi:

print(`The sequence`,a[n]=p[n](1),`satisfies the recursion`):
print(seq(a[i]=R[1][i],i=1..nops(R[1]))):
print(a[n]=add(R[2][i]*a[n-i],i=1..nops(R[2]))):
print(`The limit of the ratio`,a[n+1]/a[n],`is `):
print(r=r1):
print(``):
print(`The expected value of the random variable given by`,p[n](q),`is `):
print( M[1] ):
print(`which is asymptotic to`):
print(Mas[1]):

for j from 2 to nops(M) do
	print(``):
	if j=2 then print(`The variance is`):
	else print(`Central moment`,j,`is `): fi:
	print( M[j] ):
	print(`which is asymptotic to`):
	print(Mas[j]):
	if j=2 then 
		if Ml[1]=infinity then print(`and there is concentration about the mean.`):
		else print(`and there is no concentration about the mean.`): fi:
	fi:
	
	if j>=3 then
		print(`and the limit of standardized central moment`,j,`is `):
		print(Ml[j]):
	fi:
od:
end:


############
############

#numeric s, numeric/symbolic d
GdnMoms1:=proc(d,s,K) local L,q,d1,M,k,dmax:
if not type(d,symbol) then
	return Moms(Gdn(q,d,s),q,K):
fi:
dmax:=2*K+floor(s/2)+5:
L:=[seq(Gdn(q,d1,s),d1=1..dmax)]:
M:=Moms(L,q,K):
[seq(GuessPoly([seq(d1,d1=1..dmax)],[seq(M[d1][k],d1=1..dmax)],d,2*k+floor(s/2),2*k+floor(s/2)),k=1..K)]:
end:




#symbolic s, numeric/symbolic d
GdnMoms2:=proc(d,s,K) local M,k,ss,L:
M:=[seq(GdnMoms1(d,ss,K) , ss=2..2+4*K+4)]:
L:=[]:
for k from 1 to K do
L:=[op(L),GuessRecPol([seq(M[ss][k],ss=1..5+4*k)],[[1,d],[1,d]],s,2,2*k)]:
od:
#L:
[seq(L[i][1]*N1+L[i][2]*N2,i=1..nops(L))]:
end:

#numeric/symbolic s, numeric/symbolic d
mdnK:=proc(d,n,K)
if type(n,`symbol`) then GdnMoms2(d,n,K):
else GdnMoms1(d,n,K):
fi:
end:


MdnK:=proc(d,n,K) local L,i:
L:=mdnK(d,n,K):
if type(n,`symbol`) then
	subs({N1=1,N2=NRatio},L):
else
	[seq(normal(L[i]/Ndn(d,n)),i=1..nops(L))]:
fi:
end:

McdnK:=proc(d,n,K) 
if K<2 then print(`ERROR: need K>=2.`): return FAIL: fi:
Straight2Central(MdnK(d,n,K)):
end:

MsdnK:=proc(d,n,K) local L,i:
if K<2 then print(`ERROR: need K>=2.`): return FAIL: fi:
L:=McdnK(d,n,K):
[seq(normal(L[i]/L[2]^(i/2)),i=1..nops(L))]:
end: