######################################################################
##DRC.txt: Save this file as  DRC.txt                                #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read DRC.txt<Enter>                                #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: 
 
print(`Created:  July 14, 2016`):
print(` This is DRC.txt  `):
print(`It is one of the packages that accompany the article `):
print(`On Dave Robbins' Constant`):
print(`by Shalosh B. Ekhad and Doron Zeilberger`):

print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are:  RaSamp `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: AlphaF, Dave, DaveE, DaveF, MomF, SampA  `):
 print(` `):



elif nops([args])=1 and op(1,[args])=AlphaF then
print(`AlphaF(L,n): inputs a list L of positive numbers, of size k, say, and outputs a list of length  N`):
print(`regarding the r.v. `):
print(`"distance between two random points in the L[1]xL[2]x...xL[k] k-Dim box`):
print(`The first entry is the  mean (Robbins' constant), the second the standard-deviation, and the remaining entries`):
print(`are the skewness, kurtosis, and higher alpha coefficients. Try:`):
print(` AlphaF([1,1,1],4); `);
print(``):

elif nops([args])=1 and op(1,[args])=Dave then
print(`Dave(L,n): inputs a list L of positive numbers, of size k, say, and outputs the n-th`):
print(`straight moment of the r.v. "distance between two random points in the L[1]xL[2]x...xL[k] k-Dim box`):
print(`To get the value of the 3D Robbins constanat, type`):
print(` evalf(Dave([1,1,1],1)); `);
print(``):

elif nops([args])=1 and op(1,[args])=DaveE then
print(`DaveE(n): Dave([1,1,1],2*n) done via summation. Try`):
print(`DaveE(4);`):


elif nops([args])=1 and op(1,[args])=DaveF then
print(`DaveF(L,n): a floating point version of Dave(F,n) (q.v.). Try`):
print(` DaveF([1,1,1],1); `);
print(``):
	

elif nops([args])=1 and op(1,[args])=MomF then
print(`MomF(L,n): inputs a list L of positive numbers, of size k, say, and outputs the n-th`):
print(`moment about the mean of the r.v. "distance between two random points in the L[1]xL[2]x...xL[k] k-Dim box`):
print(`To get the value of the 3D Robbins constanat, type`):
print(` MomF([1,1,1],1); `);
print(``):

elif nops([args])=1 and op(1,[args])=RaSamp then
print(`RaSamp(L,K): inputs a list of positive numbers L, and a positive integer K, outputs a list of distances between K random pairs`):
print(`Try:`):
print(`RaSamp([1,1,1],10);`):


elif nops([args])=1 and op(1,[args])=SampA then
print(`SampA(L,K,N): inputs a list of positive numbers L, and  positive integers K and N, outputs a list of length N`):
print(`consisting of (i) the mean (ii) standard deviation, the 3rd-through N-th alpha coefficients. `):
print(`Try:`):
print(`SampA([1,1,1],1000,4);`):


else
print(`There is no ezra for`,args):
fi:
 
end:


Dave:=proc(L,n) local k,u,i,gu,lu:
option remember:
k:=nops(L):
gu:=mul((L[i]-u[i])*2*L[i],i=1..nops(L)):
lu:=(add(u[i]^2,i=1..k))^(n/2)*gu:
for i from 1 to k do
 lu:=int(lu,u[i]=0..L[i]):
od:

lu:

end:

DaveF:=proc(L,n) local k,u,i,gu,lu:
option remember:
k:=nops(L):
gu:=mul((L[i]-u[i])*2/L[i]^2,i=1..nops(L)):
lu:=(add(u[i]^2,i=1..k))^(n/2)*gu:
for i from 1 to k do
 lu:=Int(lu,u[i]=0..L[i]):
od:

evalf(lu):

end:



MomF:=proc(L,n) local mu,k,u,i,gu,lu:
option remember:
k:=nops(L):
mu:=DaveF(L,1):
gu:=mul((L[i]-u[i])*2/L[i]^2,i=1..nops(L)):
lu:=((add(u[i]^2,i=1..k))^(1/2)-mu)^n*gu:
for i from 1 to k do
 lu:=Int(lu,u[i]=0..L[i]):
od:

evalf(lu):

end:




AlphaF:=proc(L,N) local mu,sig,i:

mu:=DaveF(L,1):
sig:=sqrt(DaveF(L,2)-mu^2):

[mu,sig, seq(MomF(L,i)/sig^i,i=3..N)]:

end:




#Dis(P,Q): The Euclidean distance between points P and Q
Dis:=proc(P,Q) local i:
sqrt(add((P[i]-Q[i])^2,i=1..nops(P))):
end:


#RaSamp(L,K): inputs a list of positive numbers L, and a positive integer K, outputs a list of distances between K random pairs
#Try:
#RaSamp([1,1,1],10);
RaSamp:=proc(L,K) local k,ra,gu,P1,P2,i,i1:

k:=nops(L):



for i from 1 to nops(k) do
 ra:=rand(0..10^20*L[i]):
od:

gu:=[]:

for i from 1 to K do
 P1:=[seq(ra[i1]()/10^20,i1=1..k)]:
 P2:=[seq(ra[i1]()/10^20,i1=1..k)]:
 gu:=[op(gu),evalf(Dis(P1,P2))]:
od:

gu:

end:






#SampA(L,K,N): inputs a list of positive numbers L, and a positive integer K, outputs the sample mean, s.d.
#sample skewmess,etc. up to the N-th alpha coeff. 
#Try:
#SampA([1,1,1],10000,4);
SampA:=proc(L,K,N) local gu,mu,sig,lu,i,j:

gu:=RaSamp(L,K):

mu:=convert(gu,`+`)/nops(gu):
sig:=sqrt(add((gu[i]-mu)^2,i=1..nops(gu))/nops(gu)):

lu:=[mu,sig]:

for j from 3 to N do
 lu:=[op(lu),add((gu[i]-mu)^j,i=1..nops(gu))/nops(gu)/sig^j]:
od:

lu:

end:

#DaveE(n): Dave([1,1,1],2*n) done via summation. Try
#DaveE(4);
DaveE:=proc(n) local a,b:

add(add(n!/(a+1)!/(b+1)!/(n-a-b+1)!/(1+2*a)/(1+2*b)/(1+2*(n-a-b)),b=0..n-a),a=0..n):

end: