###################################################################### ## Noga12.txt Save this file as Noga12.txt to use it, # # stay in the # ## same directory, get into Maple (by typing: maple ) # ## and then type: read `Noga12.txt`: # ## Then follow the instructions given there # ## # ## Written by Doron Zeilberger, Rutgers University , # ## DoronZeil at gmail dot com # ###################################################################### print(`First Written: Sept. 2021: tested for Maple 2020 `): print(`This version, Oct. 14, 2021, thanks to T. Kyle Petersen.`): print(`Previous version, Oct. 4, 2021, thanks to Stoyan Dimitrov, who told us about the paper by T. Kyle Patersen nand Mathieu Guay-Paquet`): print(): print(`This is Noga12.txt, A Maple package`): print(`accompanying Shalsoh B. Ekhad and Doron Zeilberger's article: `): print(` An Experimental (yet fully rigorous!) Study of a certain "Measure Of Disarray" that 12-year Noga Alon Proved was always Even `): print(): print(`The most current version is available on WWW at:`): print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/Noga12.txt .`): print(`Please report all bugs to: DoronZeil at gmail dot com .`): print(): print(`For general help, and a list of the MAIN functions,`): print(` type "ezra();". For specific help type "ezra(procedure_name);" `): print(``): print(`-------------------------`): print(`For a list of the supporting functions type: ezra1();`): print(): print(`-------------------------`): print(``): print(`-------------------------`): print(`For a list of the STORY functions type: ezraSt();`): print(): print(`-------------------------`): print(``): print(`-------------------------`): print(`For a list of the PRE-Computed functions type: ezraPC();`): print(): print(`-------------------------`): print(``): print(`-------------------------`): print(`For a list of the SIMULATIONS functions type: ezraSi();`): print(): print(`-------------------------`): ezraSi:=proc() if args=NULL then print(`The SIMULATIONS procedures are`): print(` NogaSi`): else ezra(args): fi: end: ezraSt:=proc() if args=NULL then print(`The STORY procedures are`): print(` JointStory, NogaStory, Persi`): else ezra(args): fi: end: ezra1:=proc() if args=NULL then print(`The SUPPORTING procedures are`): print(` AveAndMoms, AveAndMomsJ, D1D2, Fabt, Fc, SeqFcf, GP, inv, LeadMomsEven, LeadMomsOdd, NPg, NPgJ, PlotDist, SEAj1, SEz `): else ezra(args): fi: end: ezraPC:=proc() if args=NULL then print(`The PRE-COMPUTED procedures are`): print(` SEAjPC, SEApc `): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(` Noga12.txt: A Maple package for investigating Spearman's footrule "measure of disarrary" and its relation to the number of inversions `): print(`The MAIN procedures are`): print(` LimNor, na, NP, NPc, NPJ, NPcJ, SEA, SEAj, SeqF, SeqFg,`): elif nargs=1 and args[1]=AveAndMoms then print(`AveAndMoms(f,x,N): Given a probability generating function`): print(`f(x) (a polynomial, rational function and possibly beyond)`): print(`returns a list whose first entry is the average `): print(`(alias expectation, alias mean)`): print(`followed by the variance, the third moment (about the mean) and`): print(`so on, until the N-th moment (about the mean).`): print(`If f(1) is not 1, than it first normalizes it by dividing`): print(`by f(1) (if it is not 0) .`): print(`For example, try:`): print(`AveAndMoms(((1+x)/2)^100,x,4);`): elif nargs=1 and args[1]=AveAndMomsJ then print(`AveAndMomsJ(f,x,y,N): Given a probability generating function`): print(`f(x,y) (a polynomial, rational function and possibly beyond)`): print(`returns a double L such that`): print(`L[1][1]=1, L[1][2] is the average w.r.t to the x-variable, L[1][i] (for 3<=i<=N+1) is the (i-1)-th moment about the mean`): print(`L[2][1] is the average w.r.t to the y-variable, L[1][i] (3<=i<=N+1) is the (i-1)-th moment about the mean of the y-variable`): print(`L[2][2] is the covariance, and more genrally L[i][j] (2<=i,j<=N+1) is the mixed (i-1,j-1) moment about the means`): print(`For example, try:`): print(`AveAndMomsJ(((x+y+1)^100,x,y,4);`): elif nargs=1 and args[1]=D1D2 then print(`D1D2(f,x,y,i,j): Given an expression f in x and y, outputs (x*Dx)^i*(y*Dy)^j f . Try:`): print(`D1D2(x^4*y^4,x,y,2,3);`): elif nargs=1 and args[1]=Fabt then print(`Fabt(a,b,t): The weight-enumerator accordring to the appraoch of Petersen and Guay-Paquet, of truncated Motzkin paths from [0,0] to [a,b].Try:`): print(`Fabt(5,0,t);`): elif nargs=1 and args[1]=Fc then print(`Fc(c,x,y): sqrt(1-c^2)/2/Pi*exp(-x^2/2-y^2/2+c*x*y)`): elif nargs=1 and args[1]=GP then print(`GP(L,n): inputs a list L , a variable name n, guesses a polynomial of degree <=nops(L)-3 in n,`): print(`let's call it P(n) such that P(i)=L[i] for all i. If it fails it returns FAIL.`): print(`Try: `): print(` GP([seq(i^3,i=1..10)],n); `): elif nargs=1 and args[1]=inv then print(`inv(pi): the number of inversions of the permutations pi. Try:`): print(`inv([2,1,4,3]);`): elif nargs=1 and args[1]=JointStory then print(`JointStory(K): An article about the mixed mooments of the two measures of disarray "number of inversions" and Spearman's footrule, and a partial (elementary!) proof`): print(`that they are joint asymptotically normal, with correlation 3/sqrt(10). Up to mixed moments (i,j) with 1<=i,j<=K. Try:`): print(`JointStory(3):`): elif nargs=1 and args[1]=LeadMomsEven then print(`LeadMomsEven(n,r): The three leading terms of the (2*r)-th moment about the mean of Spearman's Footrule on S_n. Try:`): print(`LeadMomsEven(n,r);`): elif nargs=1 and args[1]=LeadMomsOdd then print(`LeadMomsOdd(n,r): The two leading terms of the (2*r+1)-th moment about the mean of Spearman's Footrule on S_n. Try:`): print(`LeadMomsOdd(n,r);`): elif nargs=1 and args[1]=SEAjPC then print(`SEAjPC(n): The pre-comuted SEAj(n,10). A list of length 4, each with four elements, such that calling the list L, L[i][j] is the (i,j) mixed moment. TYPE:`): print(`SEAjPC(n);`): elif nargs=1 and args[1]=LimNor then print(`LimNor(c,N): the list of lists of length N, let's call it L, such that its L[i][j] is the scaled momement of the joint normal distribution with`): print(` pdf sqrt(1-c^2)/2/Pi*exp(-x^2/2-y^2/2-c*x*y). Try: `): print(` LimNor(1/2,4); `): elif nargs=1 and args[1]=na then print(`na(pi): The Noga Alon permutation index (alias Spearman's footrule). Try:`): print(`na([3,1,2,4]);`): elif nargs=1 and args[1]=NP then print(`NP(n,q): The weight-enumerator of S_n acccording to q^(na(pi)) done direcly (by pure brute force). Try:`): print(`NP(7,q);`): elif nargs=1 and args[1]=NPc then print(`NPc(n,q): The weight-enumerator of S_n acccording to q^(na(pi)) done cleverly, via dynamical programming. Try:`): print(`NPc(10,q);`): elif nargs=1 and args[1]=NPcJ then print(`NPcJ(n,p,q): The JOINT weight-enumerator of Sn acccording to p^inv(pi)*q^(na(pi)), done cleverly`. Try): print(`NPcJ(10,p,q);`): elif nargs=1 and args[1]=NogaSi then print(`NogaSi(n1,K,m): picks K permutations of {1,...,n1} and computes their empiricacl moments with respect to na(pi) (Spearman's footrule), up to the m moments. It also returns the theoretical value. Try:`): print(`NogaSi(100,100,4);`): elif nargs=1 and args[1]=NogaStory then print(`NogaStory(K): An article about the first K moments of 12-year-old Noga Alon's "measure of disarray" (aka Spearman footrule). try:`): print(`NogaStory(10):`): elif nargs=1 and args[1]=NPg then print(`NPg(L1,L2,q): The weight-enumerator according to the Spearman footrule statistics for bijections from the list L1 to the list L2, using Dynamical programming. Try:`): print(`It is needed for the clever way for computing the genearting function for Spearman's footrule for the symmetric group. Try:`): print(`NPg([1,2],[3,2],q); `): elif nargs=1 and args[1]=NPgJ then print(`NPgJ(L1,L2,p,q): The Joint weight-enumerator according to the Spearman footrule statistics for bijections from the list L1 to the list L2, using Dynamical programming. Try:`): print(`It is needed for the clever way for computing the genearting function for Spearman's footrule for the symmetric group. Try:`): print(`NPgJ([1,2],[3,2],p,q); `): elif nargs=1 and args[1]=NPJ then print(`NPJ(n,p,q): The weight-enumerator of S_n acccording to p^(inv(pi))*q^(na(pi)) done direcly. Try:`): print(`NPJ(7,p,q);`): elif nargs=1 and args[1]=Persi then print(`Persi(): Outputs an article theorem about the leading asympototic for the even and odd moments`): print(`Try: `): print(`Persi():`): elif nargs=1 and args[1]=plotDist then print(`plotDist(f,x,K): Given a prob. gen. function f(x) that has a `): print(`Taylor series`): print(`for a discrete r.v.`): print(`plots its normalized version (X-mu)/sig between mu-K1*sig and`): print(`mu+K2*sig`): print(`For example, try:`): print(`plotDist((1+x)^40,x,5,5);`): elif nargs=1 and args[1]=SEA then print(`SEA(n,K): The symbolic first few moments, in n, of Spearman's footrule, using the generating functions up to n=K. Try:`): print(`SEA(n,10);`): elif nargs=1 and args[1]=SEAj then print(`SEAj(n,N): The list of length N, each of length N , let's call it Lm such that L[i][j] is the mixed (i,j), moment of (inv,na) on Sn`): print(`SEAj(n,3);`): elif nargs=1 and args[1]=SEAj1 then print(`SEAj1(n,K,i,j): The symbolic mixed (i,j) moment of (inv,na) on S_n, using K data points. RETURNS FAIL if it does. Try:`): print(`SEAj1(n,10,1,1);`): elif nargs=1 and args[1]=SEApc then print(`SEApc(n):The pre-computed lisf, of length 16 whose i-th entry is the explicit expression for the ith moment(about the mean) of Spearman's footrule. `): print(`In other words it is the output of SEA(n,16). Try:`): print(`SEApc(n);`): elif nargs=1 and args[1]=SEz then print(`SEz(n,z,K): The first few terms in the Noga12 polynomials with q=1+z, divided by n!, in terms of z , using K data points. Try:`): print(`SEz(n,z,10);`): elif nargs=1 and args[1]=SeqF then print(`SeqF(q,N): The first N terms of the weight-enumerators for Spearman's footrule using dynamical programming inspired by`): print(`T. Kyle Petersen and Mathieu Guay-Paquet in Elc. J. Combinatorics 21 (2014), https://arxiv.org/abs/1404.4674. Try:`): print(`SeqF(q,25);`): elif nargs=1 and args[1]=SeqFcf then print(`SeqFcf(q,N): The first N terms of the weight-enumerators for Spearman's footrule using the continued fraction given by`): print(`T. Kyle Petersen and Mathieu Guay-Paquet in Elc. J. Combinatorics 21 (2014), https://arxiv.org/abs/1404.4674. Try:`): print(`SeqFcf(q,25);`): elif nargs=1 and args[1]=SeqFg then print(`SeqFg(p,q,N): The first N terms of the bi-variate weight-enumerator of S_n according to (inv, SpearmanFootRule). Try:`): print(`SeqFg(p,q,22);`): else print(`There is no such thing as`, args): fi: end: with(combinat): #start pre-computed functions #SEApc(n):The pre-computed lisf, of length 16 whose i-th entry is the explicit expression for the ith moment(about the mean) of Spearman's footrule. #In other words it is the output of SEA(n,16). Try: #SEApc(n); SEApc:=proc(n): [1/3*(n-1)*(n+1), 1/45*(n+1)*(2*n^2+7), -2/945*(n+2)*(n+1)*(2*n^2+31), 1/4725*( n+1)*(28*n^5+180*n^3+160*n^2+887*n+1265), -4/93555*(n+2)*(n+1)*(44*n^5-10*n^4+ 788*n^3+86*n^2+3587*n+8555), 1/127702575*(n+1)*(168168*n^8-145288*n^7+1800148*n ^6+2180892*n^5+18508182*n^4+32547228*n^3+112117257*n^2+385870348*n+368963105), -2/18243225*(n+2)*(n+1)*(8008*n^8-11648*n^7+171164*n^6-88560*n^5+1645002*n^4+ 2988888*n^3+4890161*n^2+46078520*n+73541545), 1/13956067125*(n+1)*(5717712*n^11 -14041456*n^10+111237120*n^9+63288800*n^8+1347724536*n^7+1996817312*n^6+ 17098013040*n^5+53375545600*n^4+125630091477*n^3+758605059019*n^2+2016696623115 *n+1690532291725), -8/5568470782875*(n+2)*(n+1)*(325909584*n^11-1185159352*n^10 +9805807128*n^9-16336808896*n^8+137996736504*n^7+46247875310*n^6+758421298152*n ^5+6095266855096*n^4+6075969746037*n^3+57909124031467*n^2+310708670730195*n+ 412170672282775), 1/102088631019375*(n+1)*(16730025312*n^14-77618212672*n^13+ 613234099632*n^12-579463689632*n^11+7657390179296*n^10+4566116902864*n^9+ 146567796905816*n^8+264922488779784*n^7+1660374425319302*n^6+11150634041788208* n^5+25328167884268777*n^4+156189484242864248*n^3+877883775625389940*n^2+ 1953357983931868000*n+1484183419415591125), -2/1152673452055125*(n+2)*(n+1)*( 164910249504*n^14-1109464214144*n^13+8370043078672*n^12-25751449518528*n^11+ 147736685187168*n^10-233209438306272*n^9+1401213332279736*n^8+5054820945342736* n^7-239641999825958*n^6+153223555779363936*n^5+437783312650259667*n^4+ 1127558386069441992*n^3+14546767596974857836*n^2+55349216195607176680*n+ 62918965273617289375), 1/7866996310276228125*(n+1)*(630286973604288*n^17-\ 4587803141201280*n^16+41490610779574528*n^15-121537834976646400*n^14+ 693899509395747536*n^13-569631402030782000*n^12+13267658973097781616*n^11-\ 4727800684317350640*n^10+247762156313623294564*n^9+967190710075660559840*n^8+ 2232743539438550658884*n^7+37931516754466359868480*n^6+144423980180862942383187 *n^5+566284059362434848409065*n^4+5218420201406733448141222*n^3+ 23349960033038734660186310*n^2+44614282739339475367534175*n+ 30910149914294298623964625), -4/605153562328940625*(n+2)*(n+1)*(30013665409728* n^17-320631452421280*n^16+2765546787907904*n^15-12975721868574400*n^14+ 65861252373298576*n^13-205623213769012000*n^12+945177411045640128*n^11-\ 336696405594869920*n^10+337102024604167364*n^9+88208812941229726190*n^8-\ 15641836046451497468*n^7+1268098670377959554690*n^6+12368049439985433951807*n^5 +24527455379512376607965*n^4+253127980781587818546686*n^3+ 2144082379078624724170880*n^2+6344155088392481577695275*n+ 6306057216658003432593875), 1/112817914119895359375*(n+1)*(5222377781292672*n^ 20-53467201094186880*n^19+568252355975761600*n^18-2910619369445686080*n^17+ 16240399684224410592*n^16-39556643830936578240*n^15+265622171313324675600*n^14-\ 685499406641656136160*n^13+7794281123984819955192*n^12-933748100246656188600*n^ 11+74394190104692822469900*n^10+1129907569789722516387060*n^9+ 1261371992814884079691902*n^8+30633006976385773149989940*n^7+ 254360895531730006950917625*n^6+857312906603211474145450680*n^5+ 7581361379779479840463886517*n^4+59594008411886192132856738780*n^3+ 216542892109305292050288122775*n^2+363873321281031483900954769500*n+ 232236148189948798919827855625), -2/161577816602514133696875*(n+2)*(n+1)*( 12465815763945608064*n^20-192679213766366046720*n^19+2052316442727432426688*n^ 18-13477912523363713860864*n^17+75285694339635219562080*n^16-\ 315292696252870168650240*n^15+1423251234545239570368816*n^14-\ 4016151458766672742516608*n^13+10715263580950195873506216*n^12+ 66465583070204167091269920*n^11-301656625762452829380092196*n^10+ 3041307788119924663676812608*n^9+7587123517483961581586627622*n^8+ 8523244654101660223249268760*n^7+745976376419625214625455244967*n^6+ 3213970814653533365091564762864*n^5+13020135347555134636994492122143*n^4+ 183112499399233038412705397900280*n^3+1086945387261633335459029516677225*n^2+ 2659931093482712374723833829044000*n+2363800180002987643438613102924375), 1/ 13734114411213701364234375*(n+1)*(423837735974150674176*n^23-\ 5671352562320778068736*n^22+69753445930917236423936*n^21-\ 512167138776107483725056*n^20+3464292937539684538762496*n^19-\ 14457503785273070066576896*n^18+68012402267830184451207936*n^17-\ 239703833035265438442364416*n^16+2073728258794139339012331936*n^15-\ 7676891024843401632452380896*n^14+46201354387869519975007747776*n^13+ 140236740228477031253470637184*n^12-368027233894590183663602299344*n^11+ 13872437454476600462673921405024*n^10+38294765208587475835849132624896*n^9+ 226527976367498963381443521223104*n^8+4734267436385491833576886884338391*n^7+ 23813518998475304409168147114413529*n^6+136942590800833244788154719236284331*n^ 5+1527369118300809152967091477946982309*n^4+ 9426174034386926559344802656500349845*n^3+ 28986560620941201396296959387271570475*n^2+ 43812040367415770828724496403813793625*n+26060087576309929501650068898306484375 )]: end: #SEAjPC(n): The pre-comuted SEAj(n,10). TYPE: #SEAjPC(n); SEAjPC:=proc(n): [[1/30*(n+1)*(n^2+1), -1/630*(n+2)*(n+1)*(2*n^2+17), 1/3150*(n+1)*(14*n^5+55*n^3+55*n^2+231*n+395), -1/155925*(n+2)*(n+1)*(220*n^5-50*n^4+2796*n^3+342*n^2+9707*n+31115)], [-1/630*(n+2)*(n+1)*(n^2+5), 1/113400*(n+1)*(392*n^5+90*n^4+380*n^3+1155*n^2+1703*n+6000), -1/1871100* (n+2)*(n+1)*(1804*n^5-270*n^4+15424*n^3+4203*n^2+33379*n+189900), 1/1702701000*(n+1)*(1289288*n^8-920348*n^7+6729828*n^6+12232052*n^5+57839362*n^4+112856303*n^3+299424692*n^2+1446621703*n+1681765120)], [1/5040*(n+1)*(14*n^5+9*n^4-33*n^3+21*n^2+n+84), -1/831600*(n+2)*(n+1)* (484*n^5+30*n^4+2256*n^3+1557*n^2+833*n+37860), 1/378378000*(n+1)*(224224*n^8-96954*n^7+406084*n^6+1779911*n^5+4930576*n^4+11406164*n^3+25523106*n^2+150493949*n+198087540), -1/1702701000*(n+2)*(n+1)*(560560*n^8-659360*n^7+4995326*n^6-5122*n^5+27772760*n^4+120508115*n^3+ 5964334*n^2+1340518327*n+3106600140)], [-1/83160*(n+2)*(n+1)*(22*n^5+9*n^4+19*n^3+81*n^2-167*n+1260), 1/13621608000*(n+1)*(6446440*n^8-82368*n^7-11985750*n^6+33565392*n^5+66314535*n^4+146271258*n^3+336035945*n^2+2204607348*n+3324827520), -1/20432412000*(n+2)*(n+1)*(4380376 *n^8-3814200*n^7+21293790*n^6+14423784*n^5+92764569*n^4+766906350*n^3-103586305*n^2+7305979236*n+20179983600), 1/3473510040000*(n+1)*(493629136*n^11-925530320*n^10+2847632940*n^9+6071639640*n^8+18325935628*n^7+39174467640*n^6+270446041895*n^5+1331971412135*n^4+ 2202855668961*n^3+17979024108405*n^2+65216298539540*n+64518734534400)]]: end: #end pre-computed functions ##Start From Histabrut #AveAndMoms(f,x,N): Given a probability generating function #f(x) (a polynomial, rational function and possibly beyond) #returns a list whose first entry is the average #(alias expectation, alias mean) #followed by the variance, the third moment (about the mean) and #so on, until the N-th moment (about the mean). #If f(1) is not 1, than it first normalizes it by dividing #by f(1) (if it is not 0) . #For example, try: #AveAndMoms(((1+x)/2)^100,x,4); AveAndMoms:=proc(f,x,N) local mu,F,memu1,gu,i: mu:=simplify(subs(x=1,f)): if mu=0 then print(f, `is neither a prob. generating function nor can it be made so`): RETURN(FAIL): fi: F:=f/mu: memu1:=simplify(subs(x=1,diff(F,x))): gu:=[memu1]: F:=F/x^memu1: F:=x*diff(F,x): for i from 2 to N do F:=x*diff(F,x): gu:=[op(gu),simplify(subs(x=1,F))]: od: gu: end: #D1D2(f,x,y,i,j): Given an expression f in x and y, outputs (x*Dx)^i*(y*Dy)^j f . try: #D1D2(x^4*y^4,x,y,2,3); D1D2:=proc(f,x,y,i,j) local gu,i1,j1,f1: f1:=f: for i1 from 1 to i do f1:=expand(x*diff(f1,x)): od: for j1 from 1 to j do f1:=expand(y*diff(f1,y)): od: f1: end: #AveAndMomsJ(f,x,y,N): Given a probability generating function #f(x,y) (a polynomial, rational function and possibly beyond) #returns a double L such that #L[1][1]=1, L[1][2] is the average w.r.t to the x-variable, L[1][i] (for 3<=i<=N+1) is the (i-1)-th moment about the mean #L[2][1] is the average w.r.t to the y-variable, L[1][i] (3<=i<=N+1) is the (i-1)-th moment about the mean of the y-variable #L[2][2] is the covariance, and more genrally L[i][j] (2<=i,j<=N+1) is the mixed (i-1,j-1) moment about the means #For example, try: #AveAndMomsJ(((x+y+1)^100,x,y,4); AveAndMomsJ:=proc(f,x,y,N) local mu,F,memu1,memu2,gu,i,j,T: mu:=simplify(subs({x=1,y=1},f)): if mu=0 then print(f, `is neither a prob. generating function nor can it be made so`): RETURN(FAIL): fi: F:=f/mu: memu1:=simplify(subs({x=1,y=1},diff(F,x))): memu2:=simplify(subs({x=1,y=1},diff(F,y))): T[1,1]:=1: T[1,2]:=memu1: T[2,1]:=memu2: F:=F/(x^memu1*y^memu2): for i from 3 to N+1 do T[i,1]:=subs({x=1,y=1},D1D2(F,x,y,i-1,0)): od: for i from 3 to N+1 do T[1,i]:=subs({x=1,y=1},D1D2(F,x,y,0,i-1)): od: for i from 2 to N+1 do for j from 2 to N+1 do T[i,j]:=subs({x=1,y=1},D1D2(F,x,y,i-1,j-1)): od: od: [seq([seq(T[i,j],j=1..N+1)],i=1..N+1)]: end: #Alpha(f,x,N): Given a probability generating function #f(x) (a polynomial, rational function and possibly beyond) #returns a list, of length N, whose #(i) First entry is the average #(ii): Second entry is the variance #for i=3...N, the i-th entry is the so-called alpha-coefficients #that is the i-th moment about the mean divided by the #variance to the power i/2 (For example, i=3 is the skewness #and i=4 is the Kurtosis) #If f(1) is not 1, than it first normalizes it by dividing #by f(1) (if it is not 0) . #For example, try: #Alpha(((1+x)/2)^100,x,4); Alpha:=proc(f,x,N) local gu,i: gu:=AveAndMoms(f,x,N): if gu=FAIL then RETURN(gu): fi: if gu[2]=0 then print(`The variance is 0`): RETURN(FAIL): fi: [gu[1],gu[2],seq(gu[i]/gu[2]^(i/2),i=3..N)]: end: ##End From Histabrut #plotDist(f,x,K): Given a prob. gen. function f(x) that has a #Taylor series #for a discrete r.v. #plots its normalized version (X-mu)/sig between mu-K1*sig and #mu+K2*sig #For example, try: #plotDist((1+x)^40,x,5,5); plotDist:=proc(f,x,K1,K2) local mu,f1,lu,gu,sig,i,j1: f1:=f/subs(x=1,f): mu:=subs(x=1,diff(f1,x)): gu:=f1/x^mu: sig:=sqrt(subs(x=1,x*diff(x*diff(gu,x),x))): lu:=taylor(f1,x=0,trunc(mu)+K2*trunc(sig)+10): lu:=[seq([i,coeff(lu,x,i)],i=max(0,trunc(mu-K1*sig))..trunc(mu+K2*sig))]: lu:=evalf([seq([(lu[j1][1]-mu)/sig,lu[j1][2]*sig],j1=1..nops(lu))]): plot(lu): end: #Start GuessPol #GP1(L,n,d): inputs a list L , a variable name n, and a positive integer d, guesses a polynomial of degree <=d in n, #let's call it P(n) such that P(i)=L[i]. The length of L must be at least d+3. If it fails it returns FAIL. #Try: #GP1([seq(i^3,i=1..10)],n,3); GP1:=proc(L,n,d) local a,i,var,eq,P: #Any d+1 data points can fit a polynomial of degree d. We want at least 2 "confirmation" points if nops(L)<=d+2 then print(`The list must be of length at least`, d+3): RETURN(FAIL): fi: #We define a generic polynomial P of n, of degree d, with UNDETERMINED coefficients. Our goal is to determine them P:=add(a[i]*n^i,i=0..d): #This is the set of unknowns with d+1 unknowns var:={seq(a[i],i=0..d)}: #We set-up a system with d+3 equations fitting the polynomial to the data eq:={seq(subs(n=i,P)-L[i],i=1..d+3)}: #Maple tries to solve the system, we rename var the solution var:=solve(eq,var): #If var=NULL then there is no solution, in other words, no polynomial of degree d fits the data if var=NULL then RETURN(FAIL): fi: #Now we plug the solution var into the generic polynomial P, making it an actual polynomial (that fits the data) P:=subs(var,P): #If nops(L) is larger than d+3, we check that the polynomial fits the remaining data for i from d+4 to nops(L) do if expand(subs(n=i,P))-L[i]<>0 then print(P, `does not agree with n=`, i): RETURN(FAIL): fi: od: #If it agrees we return the polynomial P: end: #GP(L,n): inputs a list L , a variable name n, guesses a polynomial of degree <=nops(L)-3 in n, #let's call it P(n) such that P(i)=L[i] for all i. If it fails it returns FAIL. #Try: #GP([seq(i^3,i=1..10)],n); GP:=proc(L,n) local d,P: for d from 0 to nops(L)-3 do #We start with d=0 and keep trying with higher and higher degree P:=GP1(L,n,d): if P<>FAIL then RETURN(P): fi: od: #If there is no such polynomial of degree <=nopd(L)-3 it returns FAIL FAIL: end: #End GuessPol #na(pi): The Noga Alon permutation index (alias Spearman's footrule). Try: #na([3,1,2,4]); na:=proc(pi) local i: add(abs(pi[i]-i),i=1..nops(pi)): end: #NP(n,q): The weight-enumerator of S_n acccording to q^(na(pi)) done direcly. Try: #NP(10,q); NP:=proc(n,q) local gu,pi: gu:=permute(n): add(q^na(pi),pi in gu): end: #NPg(L1,L2,q): The weight-enumerator according to the Spearman footrule statistics for bijections from the list L1 to the list L2, using Dynamical programming. Try: #It is needed for the clever way for computing the genearting function for Spearman's footrule for the symmetric group. Try: #NPg([1,2],[3,2],q); NPg:=proc(L1,L2,q) local i,gu,t1,k,L1a,L2a: option remember: if nops(L1)<>nops(L2) then RETURN(FAIL): fi: if L1=[] then RETURN(1): fi: gu:=0: k:=nops(L1): t1:=L1[1]: for i from 1 to k do L1a:=[op(2..k,L1)]: L2a:=[op(1..i-1,L2),op(i+1..k,L2)]: gu:=expand(gu+q^abs(t1-L2[i])*NPg(L1a,L2a,q)): od: gu: end: #NPc(n,q): The weight-enumerator of Sn acccording to q^(na(pi)), done cleverly NPc:=proc(n,q) local i: option remember: NPg([seq(i,i=1..n)],[seq(i,i=1..n)],q): end: #SE(n,q,K): The first few terms in the Noga12 polynomials, using K data points. Try: #SE(n,q,18); SE:=proc(n,q,K) local gu,P,j,mu,i: gu:=[seq(NPc(i,q),i=1..K)]: P:=1: for j from 1 to K do mu:=[seq(coeff(gu[i],q,j),i=j..K)]: mu:=GP(mu,n): if mu=FAIL then RETURN(P): else mu:=factor(subs(n=n-j+1,mu)): P:=P+mu*q^j: fi: od: P: end: #SEz(n,z,K): The first few terms in the Noga12 polynomials with q=1+z, divided by n!, in terms of z , using K data points. Try: #SEz(n,z,18); SEz:=proc(n,z,K) local gu,P,j,mu,i,q: gu:=[seq(NPc(i,q)/i!,i=1..K)]: gu:=expand(subs(q=1+z,gu)): P:=1: for j from 1 to K do mu:=[seq(coeff(gu[i],z,j),i=j..K)]: mu:=GP(mu,n): if mu=FAIL then RETURN(P): else mu:=factor(subs(n=n-j+1,mu)): P:=P+mu*z^j: fi: od: P: end: #SEAold(n,K): The symbolic first few moments, in n, of Spearman's footrule, using the generating functions up to n=K. Try: #SEAold(n,10); SEAold:=proc(n,K) local gu,lu,j,mu,i,q: gu:=[seq(NPc(i,q)/i!,i=1..K)]: lu:=[]: for j from 1 to K do mu:=[seq(AveAndMoms(gu[i],q,j)[j],i=j..K)]: mu:=GP(mu,n): if mu=FAIL then RETURN(lu): else mu:=factor(subs(n=n-j+1,mu)): lu:=[op(lu),mu]: fi: od: lu: end: #SEA(n,K): The symbolic first K moments, in n, of Spearman's footrule, using the generating functions up to n=K. Try: #SEA(n,10); SEA:=proc(n,K1) local K,gu,lu,j,mu,mu1,i,q: K:=K1+3: gu:=SeqF(q,3*K): gu:[seq(gu[i]/i!,i=1..nops(gu))]: lu:=[seq(AveAndMoms(gu[i],q,K),i=1..3*K)]: mu:=[]: for j from 1 to K1 do mu1:=[seq(lu[i][j],i=j+1..nops(lu))]: mu1:=GP(mu1,n): if mu1=FAIL then RETURN(mu): else mu1:=factor(subs(n=n-j,mu1)): mu:=[op(mu),mu1]: fi: od: mu: end: #####START JOINT DISTIBUTION #inv(pi): the number of inversions of the permutations pi. Try: #inv([2,1,4,3]); inv:=proc(pi) local i,j,co: co:=0: for i from 1 to nops(pi) do for j from i+1 to nops(pi) do if pi[i]>pi[j] then co:=co+1: fi: od: od: co: end: #NPJ(n,p,q): The weight-enumerator of S_n acccording to p^(inv(pi))*q^(na(pi)) done direcly. Try: #NPJ(10,p,q); NPJ:=proc(n,p,q) local gu,pi: gu:=permute(n): add(p^(inv(pi))*q^na(pi),pi in gu): end: #NPgJ(L1,L2,p,q): The Joint weight-enumerator according to the Spearman footrule statistics for bijections from the list L1 to the list L2, using Dynamical programming. Try: #It is needed for the clever way for computing the genearting function for Spearman's footrule for the symmetric group. Try: #NPgJ([1,2],[3,2],p,q); NPgJ:=proc(L1,L2,p,q) local i,gu,t1,k,L1a,L2a,co,i1: option remember: if nops(L1)<>nops(L2) then RETURN(FAIL): fi: if L1=[] then RETURN(1): fi: gu:=0: k:=nops(L1): t1:=L1[1]: for i from 1 to k do L1a:=[op(2..k,L1)]: L2a:=[op(1..i-1,L2),op(i+1..k,L2)]: co:=0: for i1 from 1 to k-1 do if L2[i]>L2a[i1] then co:=co+1: fi: od: gu:=expand(gu+p^co*q^abs(t1-L2[i])*NPgJ(L1a,L2a,p,q)): od: gu: end: #NPcJ(n,p,q): The JOINT weight-enumerator of Sn acccording to p^inv(pi)*q^(na(pi)), done cleverly NPcJ:=proc(n,p,q) local i: option remember: NPgJ([seq(i,i=1..n)],[seq(i,i=1..n)],p,q): end: #SEAj1(n,K,i,j): The symbolic mixed (i,j) moment of (inv,na) on S_n, using K data points. RETURNS FAIL if it does. Try: #SEAj1(n,K,i,j): Try #SEAj1(n,10,1,1); SEAj1:=proc(n,K,i,j) local p,q,i1,gu: option remember: gu:=[seq(NPcJ(i1,p,q)/i1!,i1=1..K)]: factor(subs(n=n-i-j+1,GP([seq(AveAndMomsJ(gu[i1],p,q, max(i+1,j+1))[i+1][j+1],i1=i+j..K)],n))): end: #SEAjOld(n,N): The list of length N, each of length N , let's call it L, such that L[i][j] is the mixed (i,j), moment of (inv,na) on Sn #SEAjOld(n,N); SEAjOld:=proc(n,N) local i,j,K: K:=5*N+3: [seq([seq(SEAj1(n,K,i,j),j=1..N)],i=1..N)]: end: #Fc(c,x,y): sqrt(1-c^2)/2/Pi*exp(-x^2/2-y^2/2+c*x*y) Fc:=proc(c,x,y): sqrt(1-c^2)/2/Pi*exp(-x^2/2-y^2/2+c*x*y):end: #LimNor(c,N): the list of lists of length N, let's call it L, such that its L[i][j] is the scaled momement of the joint normal distribution with #pdf sqrt(1-c^2)/2/Pi*exp(-x^2/2-y^2/2-c*x*y). Try: #LimNor(1/2,4); LimNor:=proc(c,N) local x,y,i,j,f,mu20,mu02: f:=Fc(c,x,y): mu20:=int(int(x^2*f,x=-infinity..infinity),y=-infinity..infinity): mu02:=int(int(y^2*f,x=-infinity..infinity),y=-infinity..infinity): [seq([seq( simplify(int(int(x^i*y^j*f,x=-infinity..infinity),y=-infinity..infinity)/mu20^(i/2)/mu02^(j/2)) ,j=1..N)],i=1..N)]: end: #START STORY PROCEDURES #NogaStory(K): An article about the first K moments of 12-year-old Noga Alon's "measure of disarray" (aka Spearman footrule). try: #NogaStory(A): NogaStory:=proc(K) local n,i,A,pi,gu: print(`Explicit Expressions for the first`, K, `moments of Spearman's Footrule and a partial (elementary!) proof of its Asymptotic Normality`): print(``): print(`By Shalosh B. Ekhad `): print(``): print(`Let's define the following "measure of disarray" on permutations of length n, let's call it A(pi)`): print(``): print(A(pi)= Sum(abs(pi[i]-i),i=1..n)): gu:=SEA(n,K): print(``): print(`The expectation of A is`): print(``): print(gu[1]): print(``): print(`and in Maple notation `): print(``): lprint(gu[1]): print(``): print(``): print(`The variance of A is`): print(``): print(gu[2]): print(``): print(`and in Maple notation `): print(``): lprint(gu[2]): print(``): for i from 3 to nops(gu) do print(`The`, i, `-th moment about the mean of A is`): print(``): print(gu[i]): print(``): print(`and in Maple notation `): print(``): lprint(gu[i]): print(``): od: print(`Let's now look at the scaled moments`): print(``): for i from 3 to nops(gu) do print(`The limit of the`, i, `-th scaled moment about the mean of A is`): print(``): lprint(limit(gu[i]/gu[2]^(i/2),n=infinity)): print(``): od: print(`This proves asymptotic normality up to the 10th moment by fully elementary way, confirming the fact that A is asymptotically normal proved`): print(` (using fancy methods) by Persi Diaconis and Ron Graham in the following paper:`): print(` Spearman's Footrule as a measure of disarray, J. Royal Stat. Soc., Section B, vol. 39, No. 2 (1977), pp. 262-268 . See Theorem 1`): print(``): print(`-------------------------------------`): print(``): print(`This ends this article that took`, time(), `seconds. to produce`): end: #JointStory(K): An article about the mixed mooments of the two measures of disarray "number of inversions" and Spearman's footrule, and a partial (elementary!) proof #that they are joint asymptotically normal, with correlation 3/sqrt(10). Try: #JointStory(): JointStory:=proc(K) local n,i,j,A,pi,gu,m20,m02,rho,lu1,lu2,x,y: print(`Explicit Expressions for Many mixed moments of the Number of Inversions and Spearman's Footrule and a partial (elementary!) proof that they are Jointly asymptotically normal with Correlation 3/sqrt(10)`): print(``): print(`By Shalosh B. Ekhad `): print(``): print(`Let's define two "measures of disarray" on permutations of length n, let's call them inv(pi) and A(pi)`): print(``): print(`inv(pi) is the familiar number of inversions, the number of pairs 1<=ipi[j] , and the less natural one, called Spearman's footrule`): print(``): print(A(pi)= Sum(abs(pi[i]-i),i=1..n)): gu:=SEAj(n,K): print(`The variance of inv is famously`): m20:=gu[1][3]: print(``): print(m20): print(``): print(`and in Maple notation`): print(``): lprint(m20): print(``): print(`The variance of A, less famously (first found by Diaconis and Graham) is`): m02:=gu[3][1]: print(``): print(m02): print(``): print(`and in Maple notation`): print(``): lprint(m02): print(``): print(`The covariance is`): print(``): print(gu[2][2]): print(``): print(`and in Maple notation`): print(``): lprint(gu[2][2]): print(``): rho:=limit(gu[2][2]/sqrt(m02*m20),n=infinity): print(``): print(`hence the asymptotic correlation is`): print(``): print(rho): print(``): for i from 1 to K-1 do print(`--------------------------------`): for j from 1 to nops(gu[i])-1 do print(``): print(`The mixed `, (i,j), `moment (about the means) of The number of inveresions and Spearman's footrule is`): print(``): print(gu[i+1][j+1]): print(``): print(`and in Maple notation `): print(``): lprint(gu[i+1][j+1]): print(``): od: od: print(``): lu1:= [seq([seq( simplify(limit(gu[i+1][j+1]/(m20^(i/2)*m02^(j/2)),n=infinity)), j=1..K-1)],i=1..K-1)]: lu2:=LimNor(3/sqrt(10),K-1): print(`hence the limit of the scaled mixed moments up to the`, K-1, `are `): print(``): print(matrix(lu1)): print(``): print(`and this happens to be EXACTLY the same as those of the joint normal distibution with correlation 3/sqrt(10), whose pdf is`): print(``): print(Fc(3/sqrt(10),x,y)): print(``): print(`that happen to be`): print(``): print(matrix(lu2)): print(``): print(`as you can see, they are indeed the same!`): print(``): print(`This gives a partial (elementary!) proof of the joint asymptotic normality of (inv,Spearman) with correlattion 3/sqrt(10)`): print(``): end: #NogaSi(n1,K,m): picks K permutations of {1,...,n} and computes their empiricacl average, variance, and scaled moments with respect to na(pi) (Spearman's footrule), up to the m moments. Try: #NogaSi(100,100,4); NogaSi:=proc(n1,K,m) local n,f,x,i,lu: lu:=SEApc25(n): f:=add(x^na(randperm(n1)),i=1..K)/K: Alpha(f,x,m),subs(n=n1,[op(1..2,lu)]): end: #Start procedures from :The Generating function for total displacement by Kyle and Mathieu #CF(A): Evalates the finite contrunued ffaction 1-A[1]/(1-A[2])... where A is a list Try # CF:=proc(A) : if A=[] then RETURN(1): fi: normal(1-A[1]/CF([op(2..nops(A),A)])): end: #KM1(q,z,N): The truncated Kyle Petersen- Mathuieu Guay-Paquet comtinued fractin KM1:=proc(q,z,N) local L,i: L:=[]: for i from 1 to N do L:=[op(L),i*q^(i-1)*z,i*q^i*z]: od: subs(q=q^2,normal(1/CF(L))): end: #SeqFcf(q,N): The first N terms of the weight-enumerators for Spearman's footrule using the continued fraction given by #T. Kyle Petersen and Mathieu Guay-Paquet in Elc. J. Combinatorics 21 (2014), https://arxiv.org/abs/1404.4674. Try: #SeqFcf(q,25); SeqFcf:=proc(q,N) local f,i,z: f:=taylor(KM1(q,z,N),z=0,N+1): [seq(expand(coeff(f,z,i)),i=1..N)]: end: #Fabt(a,b,t): The weight-enumerator accordring to the appraoch of Petersen and Guay-Paquet, of truncated Motzkin paths from [0,0] to [a,b].Try: #Fabt(5,0,t); Fabt:=proc(a,b,t) option remember: if a=0 then if b=0 then RETURN(1): else RETURN(0): fi: fi: if b<0 then RETURN(0): fi: expand(t^b*(2*b+1)*Fabt(a-1,b,t)+t^b*b*Fabt(a-1,b-1,t)+t^b*(b+1)*Fabt(a-1,b+1,t)): end: SeqF:=proc(q,N) local t,i: subs(t=q^2,[seq(Fabt(i,0,t),i=1..N)]):end: #End procedures from :The Generating function for total displacement by Kyle and Mathieu ###start with inversions #KM1g(p,q,z,N): The truncated Kyle Petersen- Mathuieu Guay-Paquet comtinued fractin KM1g:=proc(p,q,z,N) local L,i: L:=[]: for i from 1 to N do L:=[op(L),normal((1-p^i)/(1-p))*(p*q)^(i-1)*z,normal((1-p^i)/(1-p))*(p*q)^i*z]: od: subs(q=q^2,normal(1/CF(L))): end: #SeqFgCF(p,q,N): The first N terms of the bi-variate weight-enumerators for the pair (inv,Spearman's footrule) using the continued fraction given by #T. Kyle Petersen and Mathieu Guay-Paquet in Elc. J. Combinatorics 21 (2014), https://arxiv.org/abs/1404.4674. Try: #SeqFgCF(p,q,25); SeqFgCF:=proc(p,q,N) local f,i,z: f:=taylor(KM1g(p,q,z,N),z=0,N+1): [seq(expand(coeff(f,z,i)),i=1..N)]: end: #FabtG(a,b,p,t): The weight-enumerator accordring to the appraoch of Petersen and Guay-Paquet, of truncated Motzkin paths from [0,0] to [a,b].Try: #FabtG(5,0,p,t); FabtG:=proc(a,b,p,t) local i: option remember: if a=0 then if b=0 then RETURN(1): else RETURN(0): fi: fi: if b<0 then RETURN(0): fi: expand( (p*t)^b*(add(p^i,i=0..b-1)+ add(p^i,i=0..b))*FabtG(a-1,b,p,t) +(p*t)^b*add(p^i,i=0..b-1)*FabtG(a-1,b-1,p,t) +(p*t)^b*add(p^i,i=0..b)*FabtG(a-1,b+1,p,t) ): end: #SeqFg(p,q,N): The first N terms of the bi-variate weight-enumerator of S_n according to (inv, SpearmanFootRule). Try: #SeqFg(p,q,22); SeqFg:=proc(p,q,N) local t,i: subs(t=q^2,[seq(FabtG(i,0,p,t),i=1..N)]):end: #SEAj(n,K): The symbolic mixed moments {i,j) for 1<=i,j<=K #SEAj(n,10); SEAj:=proc(n,K1) local K,gu,lu,j1,j2,mu11,mu,mu1,i,p,q: K:=K1+5: gu:=SeqFg(p,q,3*K): gu:[seq(gu[i]/i!,i=1..nops(gu))]: lu:=[seq(AveAndMomsJ(gu[i],p,q,K),i=1..3*K)]: mu:=[]: for j1 from 1 to K1 do mu1:=[]: for j2 from 1 to K1 do mu11:=[seq(lu[i][j1][j2],i=j1+j2+1..nops(lu))]: mu11:=GP(mu11,n): if mu11=FAIL then print(`could not do the`, [j1,j2], `mixed moment `): RETURN(mu): else mu11:=factor(subs(n=n-j1-j2,mu11)): mu1:=[op(mu1),mu11]: fi: od: mu:=[op(mu),mu1]: od: mu: end: #GP1tight(L,n,d): inputs a list L , a variable name n, and a positive integer d, guesses a polynomial of degree <=d in n, #let's call it P(n) such that P(i)=L[i]. The length of L must be at least d+3. If it fails it returns FAIL. #Try: #GP1tight([seq(i^3,i=1..10)],n,3); GP1tight:=proc(L,n,d) local a,i,var,eq,P: #Any d+1 data points can fit a polynomial of degree d. We want at least 2 "confirmation" points if nops(L)<=d+1 then print(`The list must be of length at least`, d+2): RETURN(FAIL): fi: #We define a generic polynomial P of n, of degree d, with UNDETERMINED coefficients. Our goal is to determine them P:=add(a[i]*n^i,i=0..d): #This is the set of unknowns with d+1 unknowns var:={seq(a[i],i=0..d)}: #We set-up a system with d+2 equations fitting the polynomial to the data eq:={seq(subs(n=i,P)-L[i],i=1..d+2)}: #Maple tries to solve the system, we rename var the solution var:=solve(eq,var): #If var=NULL then there is no solution, in other words, no polynomial of degree d fits the data if var=NULL then RETURN(FAIL): fi: #Now we plug the solution var into the generic polynomial P, making it an actual polynomial (that fits the data) P:=subs(var,P): #If nops(L) is larger than d+2, we check that the polynomial fits the remaining data for i from d+3 to nops(L) do if expand(subs(n=i,P))-L[i]<>0 then print(P, `does not agree with n=`, i): RETURN(FAIL): fi: od: #If it agrees we return the polynomial P: end: #LeadMomsEven(n,r): The three leading terms of the (2*r)-th moment about the mean of Spearman's Footrule on S_n. Try: #LeadMomsEven(n,r); LeadMomsEven:=proc(n,r) local gu,i,p0,p1,p2,lu: gu:=SEApc(n): lu:=[seq(coeff(gu[2*i],n,3*i)/((2/45)^i*(2*i)!/(i!*2^i)),i=1..8)]: p0:=GP1tight(lu,r,0): lu:=[seq(coeff(gu[2*i],n,3*i-1)/((2/45)^i*(2*i)!/(i!*2^i)),i=1..8)]: p1:=factor(GP1tight(lu,r,3)): lu:=[seq(coeff(gu[2*i],n,3*i-2)/((2/45)^i*(2*i)!/(i!*2^i)),i=1..8)]: p2:=factor(GP1tight(lu,r,6)): (2/45)^r*(2*r)!/((2^r*r!))*n^(3*r)*(p0+p1/n+p2/n^2): end: #LeadMomsOdd(n,r): The three leading terms of the (2*r+1)-th moment about the mean of Spearman's Footrule on S_n. Try: #LeadMomsOdd(n,r); LeadMomsOdd:=proc(n,r) local gu,i,p0,p1,lu: gu:=SEApc(n): lu:=[seq(coeff(gu[2*i+1],n,3*i+1)/((2/45)^i*(2*i)!/(i!*2^i)),i=1..7)]: p0:=factor(GP1tight(lu,r,2)): lu:=[seq(coeff(gu[2*i+1],n,3*i)/((2/45)^i*(2*i)!/(i!*2^i)),i=1..7)]: p1:=factor(GP1tight(lu,r,5)): (2/45)^r*(2*r)!/((2^r*r!))*n^(3*r+1)*(p0+p1/n): end: #Persi(): Outputs an article theorem about the leading asympototic for the even and odd moments #Try: #Persi(): Persi:=proc() local gu,n,r: print(`Leading Asympotitcs to order 3 of the moments of the Spearman's Footrule `): print(``): print(`By Shalosh B. Ekhad `): print(``): gu:=LeadMomsEven(n,r): print(``): print(`Thereom: The first leading terms in the expression of the (2r)-th moment of Spearman's Footrule are`): print(``): print(gu): print(``): print(`and in Maple notation`): print(``): lprint(gu): print(``): gu:=LeadMomsOdd(n,r): print(`Also the first two terms in the exrpression for the (2r+1)-th moments are`): print(``): print(gu): print(``): print(`and in Maple notation`): print(``): lprint(gu): print(``): print(`Note that this immediately implies asymptotic normality, as first proved by Persin Diaconis, in a joint paper with Ron Graham`): print(``): print(` Spearman's Footrule as a measure of disarray, J. Royal Stat. Soc., Section B, vol. 39, No. 2 (1977), pp. 262-268 . See Theorem 1`): print(``): print(`but gives much more!`): print(``): end: