###################################################################### ##HISTABRUT: Save this file as HISTABRUT # ## To use it, stay in the # ##same directory, get into Maple (by typing: maple ) # ##and then type: read HISTABRUT2 # ##Then follow the instructions given there # ## # ##Written by Doron Zeilberger, Rutgers University , # #zeilberg at math dot rutgers dot edu # ###################################################################### #Created: Nov. 2009-Feb. 2010 #First posted: Aug. 19, 2010 print(` This is HISTABRUT `): print(`Written by Doron Zeilberger, Rutgers University`): print(`zeilberg at math dot rutgers dot edu`): print(`First posted: Aug. 20, 2010 `): print(`It is accompanied by the article `): print(`HISTABRUT: A Maple package for Symbol-Crunching in Probability Theory`): print(`by Doron Zeilberger`): print(`that is exclusively published in the Personal Journal of`): print(`Shalosh B. Ekhad and Doron Zeilberger, and arxiv.org .`): 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(`For a list of the supporting procedures type EzraSupporting(); `): print(` , for help with`): print(`a specific procedure, type Ezra(procedure_name); .`): print(`-----------------------------`): print(`For a list of the MAIN procedures type Ezra();, for help with`): print(`a specific procedure, type Ezra(procedure_name); .`): print(`-----------------------------`): print(`For a list of the guessing procedures type EzraGuess();,`): print(` for help with`): print(`a specific procedure, type Ezra(procedure_name); .`): print(`-------------------------------------------------`): print(`For a list of procedures with probability generating functions`): print(`type EzraPGF(); `): print(` for help with`): print(`a specific procedure, type Ezra(procedure_name); .`): print(`-------------------------------------------------`): print(`-----------------------------`): print(`For a list of the Most Important procedures for uni-variate `): print(` distributions, type: `): print(` EzraUni();, for help with`): print(`a specific procedure, type Ezra(procedure_name); .`): print(`-----------------------------`): print(``): print(`For a list of the Most Important procedures for bi-variate `): print(` distributions, type: `): print(` EzraBV();, for help with`): print(`a specific procedure, type Ezra(procedure_name); .`): print(`-----------------------------`): print(``): with(combinat): EzraPGF:=proc() if args=NULL then print(` The procedures with prob. generating functions for `): print(` famous (and less famous) random variables are: CheckZ80`): print(` ChungFeller, Faxy, Feller, Fn, FnSlow,GnD,Cor, GJst, GnD, GRpgf, `): print(` qbin, qfac, RSseq, RS50 `): else Ezra(args): fi: end: EzraGuess:=proc() if args=NULL then print(` The guessing procedures are: Findrec `): print(` GuessPOL, GuessPOL1V1, GuessPOL1V, , GuessPOL1Ve, GuessQPOL1Vk`): print(` RatToPol `): else Ezra(args): fi: end: EzraSupporting:=proc() if args=NULL then print(` The supporting procedures are: AlphaA, AlphaA2,`): print(` AveAndMomsSgfQuasiL, FaM `): print(` FaMrgf, GuessPOL1V1, GuessPOL1V, LQuasi , Milim `): else Ezra(args): fi: end: EzraUni:=proc() if args=NULL then print(` The most important procedures for uni-variate distributions are: `): print(` YofD, YofR, Yofseq, plotDist `): else Ezra(args): fi: end: ezra2:=proc() if args=NULL then print(` The most important procedures about bi-variate distributions are: `): print(` AlphaD2,AlphaAseq2, Argf2, CheckANrgf2, TwoWordsStat `): else Ezra(args): fi: end: EzraBV:=proc() if args=NULL then print(` The procedures dealing with bi-variate distributions are: `): print(` AlphaD2,AlphaAseq2, AnalyseMoms2, AnalyseMoms2V, Argf2`): print(` AveAndMoms2, AveAndMomsD2, AveAndMomsSrgf2 `): print(` AveAndMomsS2, AveAndMomsSseq2,`): print(` CheckANrgf2, CorContestP, CorContestW, FaM2 `): print(` FaMrgf2, LD, TwoWordsStat `): else Ezra(args): fi: end: Ezra:=proc() if args=NULL then print(`The main procedures are: Alpha, `): print(`AlphaD,AlphaAseq,`): print(` AlphaAgf, AlphaArgf, AlphaAgfQuasi `): print(` AveAndMoms, AveAndMomsD, `): print(` AveAndMomsS,AveAndMomsSseq, AveAndMomsSgf `): print(` AveAndMomsSgfQuasi, AveAndMomsSrgf, YofD, Yofseq, YofR `): print(` plotDist `): elif nops([args])=1 and op(1,[args])=Alpha then print(`Alpha(f,x,N): Given a probability generating function`): print(`f(x) (a polynomial, rational function and possibly beyond)`): print(`returns a list, of length N, whose `): print(`(i) First entry is the average`): print(`(ii): Second entry is the variance`): print(`for i=3...N, the i-th entry is the so-called alpha-coefficients`): print(`that is the i-th moment about the mean divided by the`): print(`variance to the power i/2 (For example, i=3 is the skewness`): print(`and i=4 is the Kurtosis)`): 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(`Alpha(((1+x)/2)^100,x,4);`): elif nops([args])=1 and op(1,[args])=AlphaA then print(`AlphaA(F,x,N,n,Deg,ORD1): Given an explicit (symbolic)`): print(`expression F(x,n), denoting a sequence of discrete`): print(`prob. distribution, parametrized by the discrete variable`): print(`n, finds the asymptotic alpha coefficients starting with the`): print(`3rd (the asymptic skewness):`): print(`to order ORD1.`): print(`It is done by Guessing (so it is SLOW!) `): print(`For example, try:`): print(`AlphaA(((1+x)/2)^n,x,4,n,10,2);`): print(`AlphaA(GRpgf(x,A,2),x,10,A,20,1);`): print(`AlphaA(GRpgf(x,A,3),x,10,A,20,1);`): elif nops([args])=1 and op(1,[args])=AlphaD then print(`AlphaD(F,x,N,n,ORD1): Given an explicit (symbolic)`): print(`expression F(x,n), denoting a sequence of discrete`): print(`prob. distribution, parametrized by the discrete variable`): print(`n, finds the asymptotic alpha coefficients starting with the`): print(`3rd (the asymptic skewness):`): print(`to order ORD1.`): print(`It is done direcly (by computing the factorial moments via`): print(`Taylor series `): print(` (so it is faster.) `): print(`For example, try:`): print(`AlphaD(((1+x)/2)^n,x,4,n,2);`): print(`AlphaD(GRpgf(x,A,2),x,20,A,4);`): elif nops([args])=1 and op(1,[args])=AlphaAseq then print(`AlphaAseq(F,x,N,n,Deg,ORD1): Given a list F`): print(`whose n-th entry is `): print(`prob. distribution, parametrized by the discrete variable`): print(`n, tries to find asymptotic expressions in n, of order ORD1`): print(`using guessed polynomial expressions for the normalized `): print(`moments of degree <=Deg`): print(`for the first N moments.`): print(`For example, try:`): print(`AlphaAseq([seq(((1+x)/2)^i,i=1..20)],x,20,n,12,2);`): print(`AlphaAseq([seq(qfac(i,q),i=1..30)],q,10,n,20,2);`): print(` AlphaAseq([seq(qfac(2*i,q)/qfac(i,q)/qfac(i+1,q),i=1..30)],q,8,n,16,2);`): print(` AlphaAseq([seq(qbin(2*i,i,q),i=1..30)],q,8,n,16,2);`): elif nops([args])=1 and op(1,[args])=AlphaA2 then print(`AlphaA2(F,x1,x2,N,n,ORD1,Deg): `): print(`Given an expression F(n,x1,x2) such that`): print(`F(n,x1,x2) is a bi-variate`): print(`prob. distribution, parametrized by the discrete variable`): print(`n , depending on the continuous variables x1 and x2`): print(`tries to find asymptotic expressions in n, of order ORD1`): print(`using guessed expressions for the normalized mixed-moments`): print(`moments of degree <=Deg`): print(`for the first N^2 mixed moments.`): print(`For example, try:`): print(`Much slower than AlphaD2(F,x1,x2,N,n,ORD1)`): print(`only use for checking!`): print(`AlphaA2(((1+x1+x2)/3)^n,x1,x2,4,n,2,10);`): elif nops([args])=1 and op(1,[args])=AlphaD2 then print(`AlphaD2(F,x1,x2,N,n,ORD1): `): print(`Given an expression F(n,x1,x2) such that`): print(`F(n,x1,x2) is a bi-variate`): print(`prob. distribution, parametrized by the discrete variable`): print(`n , depending on the continuous variables x1 and x2`): print(`tries to find asymptotic expressions in n, of order ORD1`): print(`for the first N^2 mixed moments.`): print(`For example, try:`): print(`AlphaD2(((1+x1+x2)/3)^n,x1,x2,4,n,2,10);`): elif nops([args])=1 and op(1,[args])=AlphaAseq2 then print(`AlphaAseq2(F,x1,x2,N,n,Deg,ORD1): `): print(`Given a sequence F`): print(`whose n-th item is a bi-variate`): print(`prob. distribution, parametrized by the discrete variable`): print(`n , depending on the continuous variables x1 and x2`): print(`tries to find asymptotic expressions in n, of order ORD1`): print(`using guessed expressions for the normalized mixed-moments`): print(`moments of degree <=Deg`): print(`for the first N^2 mixed moments.`): print(`For example, try:`): print(`AlphaAseq2([seq(((1+x1+x2)/3)^i,i=1..20)],x1,x2,4,n,10,2);`): print(`To get the first 6 normalized moments of the (inv,maj) `): print(`bi-distribution on permutations, type:`): print(`AlphaAseq2([seq(GnD(i,1,q1,q2),i=1..20)],q1,q2,4,n,10,1);`): elif nops([args])=1 and op(1,[args])=AlphaAgf then print(`AlphaAgf(F,s,t,n,N,Deg,Ord1): given a rational function F(s,t)`): print(`such that the coeff. of s^n is the prob. gen. function `): print(`of a discrete r.v. parametrized by n`): print(`finds the asymptotic alpha up `): print(`to the N-th coefficients as n goes to infinity`): print(`(starting with alpha_3 (the skewness)), with order Ord1 .`): print(`it is done by pure guessing of polynomials.`): print(`For example, try:`): print(`AlphaAgf(1/(1-s*(1+t)/2),s,t,n,6,10,1);`): elif nops([args])=1 and op(1,[args])=AlphaAgfQuasi then print(`AlphaAgfQquasi(F,s,t,n,N,Deg,k): given a rational function F(s,t)`): print(`such that the coeff. of s^n is the prob. gen. function `): print(`of a discrete r.v. parametrized by n`): print(`finds the asymptotic alpha up `): print(`to the N-th coefficients as n goes to infinity`): print(`(starting with alpha_3 (the skewness)) .`): print(`it is done by pure guessing of quasi-polynomials of period k`): print(`For example, try:`): print(`AlphaAgfQquasi(1/(1-s*(1+t)/2),s,t,n,6,10,1);`): print(`AlphaAgfQquasi(Feller(z,t,1,1,1),z,t,n,4,20,2);`): elif nops([args])=1 and op(1,[args])=AlphaArgf then print(`AlphaArgf(F,s,t,n,N,ORD1): given a rational function F(s,t)`): print(`such that the coeff. of s^n is the prob. gen. function `): print(`of a discrete r.v. parametrized by n`): print(`finds the asymptotic alpha up `): print(`to the N-th coefficients as n goes to infinity`): print(`up to order ORD1`): print(`Without GUESSING`): print(`(starting with alpha_3 (the skewness)) .`): print(`For example, try:`): print(`AlphaArgf(1/(1-s*(1+t)/2),s,t,n,6,2);`): elif nops([args])=1 and op(1,[args])=AnalyseMoms2 then print(`AnalyseMoms2(L,n,ORD1,r,s): Given a list of lists L of polynomials`): print(`in n, such that L[i][j] equals the mixed moment M[i-1,j-1]`): print(`and we have that M[2*r,2*s+1]=0 and M[2*r+1,2*s]=0 `): print(`and the dist. is believed to be `): print(`a bivariate distribution that is joint-independently`): print(`asympotically normal, finds the asymptotics,`): print(`as expressions, to order ORD1, of M[2*r,2*s] and M[2*r+1,2*s+1]`): print(`expressed in terms of r and s`): print(`The output is`): print(`(i) the variance of the first r.v.`): print(`(ii) the variance of the second r.v.`): print(`(iii) the covariance`) : print(`(iv) the asympotics to order ORD1 of the [2*r,2*s] moment`): print(`(v) the asympotics to order ORD1 of the [2*r+1,2*s+1] moment`): print(`(vi) the asymptotics to order ORD1 `): print(` of the normalized moments [2*r,2*s]`): print(`(vii) the asymptotics to order ORD1 of`): print(` the normalized moments [2*r-1,2*s-1]`): print(`For example, do:`): print(`read oHISTABRUT20:`): print(`AnalyseMoms2(InvMajMoms10[3],n,1,r,s);`): elif nops([args])=1 and op(1,[args])=AnalyseMoms2V then print(`AnalyseMoms2V(L,n,ORD1,r,s): Verbose version of`): print(`AnalyseMoms2(L,n,ORD1,r,s). For example, try:`): print(`read oHISTABRUT20: AnalyseMoms2V(InvMajMoms10[3],n,1,r,s);`): elif nops([args])=1 and op(1,[args])=Argf2 then print(`Argf2(F,s,x1,x2,N,n,ORD1): given a rational function`): print(` F(s,x1,x2)`): print(`such that the coeff. of s^n is the prob. gen. function `): print(`of a bi-variate discrete r.v. parametrized by n`): print(`finds the asymptotic alpha up `): print(`to the [N,N]-th coefficients as n goes to infinity`): print(`up to order ORD1`): print(`Without GUESSING`): print(`Argf2(1/(1-s*(1+x1+x2)/3),s,x1,x2,6,n,2);`): print(`Argf2(GJst({1,2},{[1,1],[2,2]},s,t),s,t[1,1],t[2,2],6,n,1);`): print(`Argf2(GJst({1,2},{[1,1,1],[2,2,2]},s,t),s,t[1,1,1],t[2,2,2],6,n,1);`): elif nops([args])=1 and op(1,[args])=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)^100,x,4);`): elif nops([args])=1 and op(1,[args])=AveAndMoms2 then print(`AveAndMoms2(f,x1,x2,N): Given a specific`): print(`probability generating function`): print(`for a bivariate distribution`): print(`f(x1,x2) (a polynomial, rational function and possibly beyond)`): print(`returns the averages of the first and second r.v. followed`): print(`by a list-of-lists L`): print(`such that L[r+1][s+1] is the (r,s)-mixed moment about the means.`); print(`For example, try:`): print(`AveAndMoms2(((1+x1+x2)/3)^100,x1,x2,4);`): elif nops([args])=1 and op(1,[args])=AveAndMomsD then print(`AveAndMomsD(F,x,N): Given an explicit (symbolic)`): print(`expression F(x,n), denoting a sequence of discrete`): print(`prob. distribution, parametrized by the discrete variable n`): print(`for the first N moments.`): print(`For example, try: `): print(`AveAndMomsD(((1+x)/2)^n,x,4);`): elif nops([args])=1 and op(1,[args])=AveAndMomsD2 then print(`AveAndMomsD2(F,x1,x2,N): Given an explicit (symbolic)`): print(`expression F(x1,x2,n), denoting a sequence of discrete`): print(`prob. distribution, parametrized by the discrete variable`): print(`n, finds the averages of the r.vs corresponding to x1 and x2`): print(`followed by the list-of-lists whose L such that L[r+1][s+1]`): print(`0<=r,s<=N is the`): print(`(r,s)-th mixed moment (about the means)`): print(`For example, try:`): print(`AveAndMomsD2(((1+x1+x2)/3)^n,x1,x2,4);`): elif nops([args])=1 and op(1,[args])=AveAndMomsSrgf2 then print(`AveAndMomsSrgf2(F,s,x1,x2,n,N): Given a rational generating function`): print(`F(s,x1,x2), such that viewed as a formal-power-series `): print(`in s, the coeff. of s^n is the bi-variate prob. generating function`): print(` of a sequence of discrete `): print(`prob. distribution, parametrized by the discrete variable`): print(`n, in x1, x2,finds the (linear)`): print(`expression (in n) for the average followed by the list`): print(`of lists (let's call it L) such that the (i,j)-th moments`): print(`about the mean`): print(`is L[i+1][j+1], 0 <=i,j <=N `): print(`For example, try:`): print(`AveAndMomsSrgf2(1/(1-(1+x1+x2)/3*s),s,x1,x2,n,4);`): elif nops([args])=1 and op(1,[args])=AveAndMomsSgf then print(`AveAndMomsSgf(F,z,x,n,N,Deg): Given a generating function in z`): print(`F(x,z), such that the Maclaurin coeff. of z^n is the p.g.f.`): print(` for a discrete`): print(`prob. distribution, parametrized by the discrete variable`): print(`n, tries to find explicit polynomial expressions, in n`): print(`of degree <=Deg (by pure guessing), `): print(`for the first N moments.`): print(`For example, try:`): print(`AveAndMomsSgf(1/(1-s*(1+t)/2),s,t,n,4,10);`): elif nops([args])=1 and op(1,[args])=AveAndMomsSrgf then print(`AveAndMomsSrgf(F,s,t,n,N): given a Rational function F(s,t) `): print(`such that the coeff. of s^n is the prob. gen. function `): print(`of a discrete r.v. parametrized by n`): print(`finds polynomial expressions in n for the first N moments.`): print(`Without GUESSING!`): print(`For example, try:`): print(`AveAndMomsSrgf(1/(1-s*(1+t)/2),s,t,n,6);`): elif nops([args])=1 and op(1,[args])=AveAndMomsSgfQuasi then print(`AveAndMomsSgfQuasi(F,x,N,n,Deg,z,k): Given an explicit (symbolic)`): print(`generating function F(x,z), such that the Maclaurin coeff. of z^n `): print(`is a prob. distribution, parametrized by the discrete variable`): print(`n, tries to find explicit quasi-polynomial expressions, in n`): print(`(of period k)`): print(`of degree <=Deg (by pure guessing)`): print(`for the first N moments`): print(`For example, try:`); print(`AveAndMomsSgfQuasi(1/(1-z*(1+x)/2),z,x,n,4,10,1);`): print(`AveAndMomsSgfQuasi(Feller(z,t,1,1,1),z,t,n,4,20,2);`): elif nops([args])=1 and op(1,[args])=AveAndMomsSgfQuasiL then print(`AveAndMomsSgfQuasiL(F,x,N,n,Deg,z,k): `): print(`the leading terms of AveAndMomsSgfQuasi(F,x,N,n,Deg,z,k))`): print(`if they exist, or it returs FAIL. `): print(`For example, try:`); print(`AveAndMomsSgfQuasiL(1/(1-z*(1+x)/2),z,x,n,4,10,1);`): print(`AveAndMomsSgfQuasiL(Feller(z,t,1,1,1),z,t,n,4,20,2);`): elif nops([args])=1 and op(1,[args])=AveAndMomsS then print(`AveAndMomsS(F,x,N,n,Deg): Given an explicit (symbolic)`): print(`expression F(x,n), denoting a sequence of discrete`): print(`prob. distribution, parametrized by the discrete variable`): print(`n, tries to find explicit polynomial expressions, in n`): print(`of degree <=Deg`): print(`for the first N moments.`): print(`For example, try:`): print(`AveAndMomsS(((1+x)/2)^n,x,4,n,10);`): elif nops([args])=1 and op(1,[args])=AveAndMomsS2 then print(`AveAndMomsS2(F,x1,x2,N,n,Deg): Given an explicit (symbolic)`): print(`expression F1(x1,x2,n), denoting a sequence of bivariate discrete`): print(`prob. distributions, parametrized by the discrete variable`): print(`n, tries to find explicit polynomial expressions, in n`): print(`of degree <=Deg`): print(`for the first N moments`): print(`For example, try:`): print(`AveAndMomsS2(((1+x1+x2)/3)^n,x1,x2,4,n,10);`): elif nops([args])=1 and op(1,[args])=AveAndMomsSseq then print(`AveAndMomsSseq(F,x,N,n,Deg): Given a list F`): print(`whose n-th entry is `): print(`prob. distribution, parametrized by the discrete variable`): print(`n, tries to find explicit polynomial expressions, in n`): print(`of degree <=Deg`): print(`for the first N moments.`): print(`For example, try:`): print(`AveAndMomsSseq([seq(((1+x)/2)^i,i=1..20)],x,4,n,10);`): elif nops([args])=1 and op(1,[args])=AveAndMomsSseq2 then print(`AveAndMomsSseq2(F,x1,x2,N,n,Deg): Given a list F`): print(`whose n-th entry is a bi-variate`): print(`prob. distribution, parametrized by the discrete variable`): print(`n, tries to find explicit polynomial expressions, in n`): print(`of degree <=Deg`): print(`for the first (N+1)^2 mixed moments moments.`): print(`For example, try:`): print(`AveAndMomsSseq2([seq(((1+x1+x2)/2)^i,i=1..20)],x1,x2,4,n,10);`): elif nops([args])=1 and op(1,[args])=CheckANrgf2 then print(`CheckANrgf2(F,s,x1,x2,N)`): print(`Given an expression F(n,x1,x2)`): print(`whose n-th item is a bi-variate`): print(`prob. distribution, parametrized by the discrete variable`): print(`n , depending on the continuos variables x1 and x2`): print(`checks whether it seems to be joint-asymptotically`): print(`normal with the right correlation.`): print(`For example, try:`): print(`CheckANrgf2(1/(1-(1+x1+x2)/3*s),s,x1,x2,4);`): print(`CheckANrgf2(GJst({1,2},{[1,2,2],[2,1,1]},s,t),s,t[1,2,2],t[2,1,1],8);`): elif nops([args])=1 and op(1,[args])=CheckANrgf2Plus then print(`CheckANrgf2Plus(F,s,x1,x2,N,r1,r2,ORD1)`): print(`Given an expression F(n,x1,x2)`): print(`whose n-th item is a bi-variate`): print(`prob. distribution, parametrized by the discrete variable`): print(`n , depending on the continuos variables x1 and x2`): print(`checks whether it seems to be joint-asymptotically`): print(`normal with the right correlation.`): print(`and also guesses polynomial expressions`); print(`for the mixed-moments`): print(`For example, try:`): print(`CheckANrgf2Plus(1/(1-(1+x1+x2)/3*s),s,x1,x2,8,r1,r2,1);`): elif nops([args])=1 and op(1,[args])=CheckZ80 then print(`CheckZ80(n,t,q1,q2): `): print(`checks Zeilberger's recurrence for Fn(t,q1,q2) from`): print(`his 1980 article for n. For example, try:`): print(`CheckZ80(4,t,q1,q2);`): elif nops([args])=1 and op(1,[args])=ChungFeller then print(`ChungFeller(z,t1,t2,t3,t4): The Chung-Feller formal-power series`): print(`whose coeff. of z^n is the weight-enumerator of all`): print(`sequences of {1,-1} of length n that sum to 0, with the weight`): print(`weight(Sequenec):=t1^(#losing times)*t2^(#returns to the origin)`): print(`*t3^(last visit to the origin)*t4^(#sign-changes)`): print(`Do: ChungFeller(z,t1,t2,t3,t4);`): elif nops([args])=1 and op(1,[args])=CorContestP then print(`CorContestP(r,k):Finds the pair of words with distinct letters`): prin(`that are most correlated`): print(`and least correlated among k-letter words in {1,...,r}.`): print(`of course r>=k`): print(`For example, try:`): print(`CorContestP(3,2);`): elif nops([args])=1 and op(1,[args])=CorContestW then print(`CorContestW(r,k):Finds the pair of words that are most correlated`): print(`and least correlated among k-letter words in {1,...,r}.`): print(`For example, try:`): print(`CorContestW(2,3);`): elif nops([args])=1 and op(1,[args])=FaM then print(`FaM(F,x,N): Given an explicit (symbolic)`): print(`expression F(x,n), denoting a sequence of discrete`): print(`prob. distribution, parametrized by the discrete variable`): print(`n (not needed as input),`): print(` finds the average followed by the 2nd- through the`): print(`N-th factorial moment (about the mean)`): print(`For example, try:`): print(`FaM(((1+x)/2)^n,x,4);`): elif nops([args])=1 and op(1,[args])=FaM2 then print(`FaM2(F,x1,x2,N): Given an explicit (symbolic)`): print(`expression F(x1,x2,n), denoting a sequence of discrete`): print(`prob. distribution, parametrized by the discrete variable`): print(`n, finds the averages of the r.vs corresponding to x1 and x2`): print(`followed by the list-of-lists whose L such that L[r+1][s+1]`): print(`0<=r,s<=N is the`): print(`(r,s)-th mixed factorial moment (about the means)`): print(`For example, try:`): print(`FaM2(((1+x1+x2)/3)^n,x1,x2,4);`): elif nops([args])=1 and op(1,[args])=FaMrgf then print(`FaMrgf(F,s,x,n,N): Given a rational generating function`): print(`F(s,x), such that viewed as a formal-power-series `): print(`in s, the coeff. of s^n is the prob. generating function`): print(` of a sequence of discrete `): print(` prob. distribution, parametrized by the discrete variable `): print(` n, finds the (linear) `): print(` expression (in n) for the average followed by the 2nd- through the `): print(` N-th factorial moment (about the mean) `): print(` For example, try: `): print(` FaMrgf(1/(1-(1+x)/2*s),s,x,n,4); `): elif nops([args])=1 and op(1,[args])=FaMrgf2 then print(`FaMrgf2(F,s,x1,x2,n,N): Given a rational generating function`): print(`F(s,x1,x2), such that viewed as a formal-power-series `): print(`in s, the coeff. of s^n is the bi-variate prob. generating function`): print(` of a sequence of discrete `): print(`prob. distribution, parametrized by the discrete variable`): print(`n, in x1, x2,finds the (linear)`): print(`expression (in n) for the average followed by the list`): print(`of lists (let's call it L) such that the (i,j)-th factorial moment`): print(`is L[i+1][j+1], 0 <=i,j <=N `): print(`For example, try:`): print(`FaMrgf2(1/(1-(1+x1+x2)/3*s),s,x1,x2,n,4);`): elif nops([args])=1 and op(1,[args])=Faxy then print(`Faxy(a,x,y): the normalized bi-variate distributin`): print(`of two normal r.v.s whose correlation is a`): elif nops([args])=1 and op(1,[args])=Feller then print(`Feller(z,t1,t2,t3,t4): The Feller formal-power series`): print(`whose coeff. of z^n is the weight-enumerator of all`): print(`sequences of {1,-1} of length n with the weight`): print(`weight(Sequenec):=t1^(#losing times)*t2^(#returns to the origin)`): print(`*t3^(last visit to the origin)*t4^(#sign-changes)`): print(`Do: Feller(z,t1,t2,t3,t4);`): elif nops([args])=1 and op(1,[args])=Findrec then print(`Findrec(f,n,N,MaxC): `): print(`Given a list f tries to find a linear recurrence equation with`): print(`poly coffs.`): print(`of maximum DEGREE+ORDER<=MaxC`): print(`e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);`): print(`Findrec(AlphaA(((1+x)/2)^n,x,20,n,10),r,R,3);`): elif nops([args])=1 and op(1,[args])=Fn then print(`Fn(n,t,q1,q2):`); print(`The sum of q1^inv(pi)*q2^maj(pi)*t^des(pi) over all n-perms`): print(`For example, try: Fn(3,t,q1,q2); `): elif nops([args])=1 and op(1,[args])=FnSlow then print(`FnSlow(n,t,q1,q2):`); print(`The sum of q1^inv(pi)*q2^maj(pi)*t^des(pi) over all n-perms`): print(`using Zeilberger's 1980 recurrence `): print(`For example, try: FnSlow(3,t,q1,q2); `): elif nops([args])=1 and op(1,[args])=GCor then print(`GCor(a,r,s): the (r,s)-moment of Faxy(a,x,y) divided by the s.d.`): print(`to the power (r+s)/2. For example, try:`): print(`GCor(1/3,2,2);`): elif nops([args])=1 and op(1,[args])=GJst then print(`GJst(alphabet,MISTAKES,s,t): given an alphabet alphabet`): print(`and a set of words in that alphabet, MISTAKES, and symbols`): print(`s and t, uses the Goulden-Jackson method to find the`): print(`ratinal function in s whose coeff. of s^n is the`): print(`weight-enumerator of all n-letter words in the alphabet`): print(`with the weight being `): print(`t[m1]^(#of times m1 showed up)*t[m2]^(#of times m2 showed up)*... .`): print(`For example, try:`): print(`GJst({1,2},{[1,2,1]},s,t);`): elif nops([args])=1 and op(1,[args])=GnD then print(`GnD(n,t,q1,q2): the tri-variate prob. gen. function`): print(`for des,inv,maj . `): print(`For example, try: GnD(4,t,1,1); `): elif nops([args])=1 and op(1,[args])=GRpgf then print(`GRpgf(x,A,k): the prob. generating function, in x, `): print(`for "duration of life" `): print(`for a fair gambler in [0,k*A] currently with A dollars`): print(`Note: A can be symbolic, but k has to be a pos. integer`): print(`For example, try:`): print(`GRpgf(x,10,2);`): print(`GRpgf(x,A,2);`): elif nops([args])=1 and op(1,[args])=GRpgfOld then print(` GRpgfOld(x,n,A): the prob. generating function, `): print(`in x, for "duration of life"`): print(`for a fair gambler in [0,A] currently with n dollars.`): print(`For example, try:`): print(`GRpgfOld(x,10,20);`): elif nops([args])=1 and op(1,[args])=GuessPOL1V then print(`GuessPOL1V(L,x): Given a list L, tries to find`): print(`a polynomial P(x),such that L[i]=P(i) for i=1..nops(L)`): print(`if there is no polynomial of degree <=nops(L)-3 that it`): print(`returns FAIL.`): print(`For example, try:`): print(`GuessPOL1V([seq(i^2,i=1..10)],x);`): elif nops([args])=1 and op(1,[args])=GuessPOL1Ve then print(`GuessPOL1Ve(L,x): Given a list L, tries to find`): print(`a polynomial P(x), and a starting place s`): print(`such that L[i]=P(i) for i=s..nops(L)`): print(`it also returns the starting place.`): print(`If there is no polynomial of degree <=nops(L)-3 that it`): print(`for starting place <=nops(L)/2, it`): print(`returns FAIL.`): print(`For example, try:`): print(`GuessPOL1Ve([0,0,0,seq(i^2,i=4..10)],x);`): elif nops([args])=1 and op(1,[args])=GuessPOL1V1 then print(`GuessPOL1V1(L,x,d): Given a list L, tries to find`): print(`a polynomial P(x), of degree d, such that L[i]=P(i) for i=1..nops(L)`): print(`For example, try:`): print(`GuessPOL1V1([seq(i^2,i=1..10)],x,2);`): elif nops([args])=1 and op(1,[args])=GuessPOL then print(`GuessPOL(L,VarList): Tries to guess a polynomial in`): print(`VarList from the data L that is a list of lists of lists...`): print(`such that P(i1,i2,i3, ...)=L[i1][i2][i3] ... .`): print(`For example, try:`): print(`GuessPOL([seq([seq(i+j,i=1..10)],j=1..10)],[x,y]);`): elif nops([args])=1 and op(1,[args])=GuessQPOL1Vk then print(`GuessQPOL1Vk(L,x,k): Given a list L, a variable x, and `): print(`a positive integer k, tries to find`): print(`polynomials [P_1(x),P_2(x), ..., P_k(x)]`): print(`such that L[k*i+j]=P_j(i) for i=1..nops(L)/k .`): print(`For example, try:`): print(`GuessQPOL1Vk([1,2,1,4,1,6,1,8,1,10],x,2);`): elif nops([args])=1 and op(1,[args])=LD then print(`LD(f,x,k,n,t): Given a sequence of discrete probability g.f.`): print(`distributions that has closed form f^n, where`): print(`f is the probability generating function for one step in `): print(`x[1], ..., x[k]`): print(`and a letter t for the continuous variables t[1], ..., t[k], finds `): print(`the expectation, variance, and the inverse-Fourier-transfrom`): print(`of the limiting distribution`): print(`in t[1], ..., t[k]. For example, try:`): print(`LD((1+x[1]+x[2]),x,2,n,t);`): elif nops([args])=1 and op(1,[args])=LQuasi then print(`LQuasiP(L,n): Given a quasi-polynomial L in n returns`): print(`the leading term, if it is consistent, otherwise returns FAIL`); print(`For example, try: `): print(`LQuasi([3*n+1,3*n+2],n);`): elif nops([args])=1 and op(1,[args])=Milim then print(`Milim(r,k): `): print(`The list of all k-letter words in the alphabet {1,2,...,r}`): print(`for example, try: `): print(`Milim(2,3);`): elif nops([args])=1 and op(1,[args])=RS50 then print(`RS50(t): the first 50 terms in the PGF for Alternating permutations`): print(`in other words, RSseq(50,y); pre-computed `): print(`For example, try: RS50(t);`): elif nops([args])=1 and op(1,[args])=RSseq then print(`RSseq(N,t): the first N terms in the PGF for Alternating permutations`): print(`For example, try: RSseq(10,t);`): elif nops([args])=1 and op(1,[args])=TwoWordsStat then print(`TwoWordsStat(alphab,w1,w2,N,n,ORD1): given two words of the same`): print(`length w1,w2, in the alphabet alphab,`): print(`finds the average, variance (as expressions in n) and the`): print(`asymptotic expressions (in n), to order ORD1, of`): print(`normalized mixed-moments (up to the Nth) of the random variables`): print(`number of occurences of w1 and number of occurrences of w2`): print(`(as factors) in a (uniform-at-random) n-letter word.`): print(`For example, try:`): print(`TwoWordsStat({1,2},[1,2],[2,1],4,n,1);`): elif nops([args])=1 and op(1,[args])=YofD then print(`YofD(F,x,N,n,ORD1,r): Given the inputs`): print(`(i) a discrete prob. generating function `): print(`parametrized by n phrased in terms of the variable x`): print(`(ii) a symbol x for the variable name`): print(`(iii) a pos. integer N for the number of moments to use for guessing`): print(`(iv) a symbol n for (i)`): print(`(v) a pos. intger ORD1, for the order of the asymptotic expansions`): print(`(vi) a symbol r`): print(`outputs:`): print(`(i) a polynomial expression for the average in terms of n`): print(`(ii) a polynomial expression for the variance in terms of n`): print(`(iii) the first ORD1 terms in the asymptotic expansion of the`): print(`even (2r)-th alpha-coeff. (m_{2r}/m_2^r)`): print(`(iv) the first ORD1 terms in the asymptotic expansion of the`); print(`even (2r+1)-th alpha-coeff. (m_{2r+1}/m_2^(r+1/2))`): print(`if it is not asymptotically normal, it returns FAIL.`): print(`For example try:`): print(`YofD(((1+x)/2)^n,x,16,n,2,r);`): elif nops([args])=1 and op(1,[args])=Yofseq then print(`Yofseq(F,x,N,n,Deg,ORD1,r): Given the inputs`): print(`(i) a list discrete probability distributions F`): print(`parametrized by n phrased in terms of the variable x`): print(`(ii) a symbol x for the variable name`): print(`(iii) a pos. integer N for the number of moments to use for guessing`): print(`(iv) a symbol n for (i)`): print(`(v) a pos. integer Deg for the max. degree of `): print(` polynomials in n to be used`): print(`(vi) a pos. intger ORD1, for the order of the asymptotic expansions`): print(`(vii) a symbol r`): print(`outputs:`): print(`(i) a polynomial expression for the average in terms of n`): print(`(ii) a polynomial expression for the variance in terms of n`): print(`(iii) the first ORD1 terms in the asymptotic expansion of`): print(` the even [(2r)-th] alpha-coeff. (m_{2r}/m_2^r), as expressions`): print(`in BOTH n and r`): print(`(iv) the first ORD1 terms in the asymptotic expansion of the`): print(`odd (2r+1)-th alpha-coeff. (m_{2r+1}/m_2^(r+1/2))`): print(`if it is not asymptotically normal, it returns FAIL.`): print(`For example try:`): print(`Yofseq([seq(((1+x)/2)^i,i=1..20)],x,16,n,10,2,r);`): print(`Yofseq(RS50(t),t,16,n,45,1,r);`): print(` Yofseq([seq(qfac(i,q),i=1..26)],q,12,n,20,1,r);`): print(` Yofseq([seq(qbin(2*i,i,q),i=1..40)],q,16,n,32,2,r);`): print(` Yofseq([seq(normal(qbin(2*i,i,q)/(1-q^(i+1))*(1-q)*(i+1)),i=1..30)],q,12,n,24,1,r);`): elif nops([args])=1 and op(1,[args])=YofR then print(`YofR(f,s,t,N,n,ORD1,r): Given the inputs`): print(`(i) a rational generating function f in the variables s and t`): print(`such that the Maclaurin coefficient of s^n is the probability`): print(`generating function of a discete r.v. parametrized by n`): print(`, in the variable t`): print(`(ii) symbos s and t for the variable names to describe f`): print(`(iii) a pos. integer N for the number of moments to use for guessing`): print(`(iv) a symbol n for (i)`): print(`(v) a pos. intger ORD1, for the order of the asymptotic expansions`): print(`(vi) a symbol r`): print(`outputs:`): print(`(i) a polynomial expression for the average in terms of n`): print(`(ii) a polynomial expression for the variance in terms of n`): print(`(iii) the first ORD1 terms in the asymptotic expansion of`): print(` the even [(2r)-th] alpha-coeff. (m_{2r}/m_2^r), as expressions`): print(`in BOTH n and r`): print(`(iv) the first ORD1 terms in the asymptotic expansion of the`): print(`odd (2r+1)-th alpha-coeff. (m_{2r+1}/m_2^(r+1/2))`): print(`if it is not asymptotically normal, it returns FAIL.`): print(`For example try:`): print(`YofR(1/(1-s*(1+t)/2),s,t,12,n,1,r);`): print(`YofR(GJst({1,2},{[2,2,2]},s,t),s,t[2,2,2],30,n,2,r);`): elif nops([args])=1 and op(1,[args])=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);`): print(`plotDist(normal(GRpgf(x,50,2)),x,3,3);`): print(`plotDist(normal(GRpgf(x,30,3)),x,3,3);`): elif nops([args])=1 and op(1,[args])=qbin then print(`qbin(a,b,q): the normalized q-binomial coefficient q-a choose b`): print(`for example, try:`): print(`qbin(10,5,q);`): elif nops([args])=1 and op(1,[args])=qfac then print(`qfac(n1,q): the normalized q-factorial for numeric n1`): print(`for example, try:`): print(`qfac(5,q);`): elif nops([args])=1 and op(1,[args])=RatToPol then print(`RatToPol(f,s,n): given a rational function of the form`): print(`P(s)/(1-s)^k, finds the polynomial expression in n,`): print(`for the coeff. of s^n (for n sufficiently large)`): print(`for example, try:`): print(`RatToPol(s^5/(1-s)^3,s,n);`): else print(`There is no Ezra for`,args): fi: end: #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: #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: #GuessPOL1V1(L,x,d): Given a list L, tries to find #a polynomial P(x) of degree d such that L[i]=P(i) for i=1..nops(L) #For example, try: #GuessPOL1V1([seq(i^2,i=1..10)],x,2); GuessPOL1V1:=proc(L,x,d) local eq,var,a,P,i: if nops(L)-d<3 then print(`Insufficient data`): RETURN(FAIL): fi: P:=add(a[i]*x^i,i=0..d): var:={seq(a[i],i=0..d)}: eq:={seq(expand(subs(x=i,P)-L[i]),i=1..nops(L))}: var:=solve(eq,var): if var=NULL then RETURN(FAIL): fi: subs(var,P): end: #GuessPOL1V(L,x): Given a list L, tries to find #a polynomial P(x) such that L[i]=P(i) for i=1..nops(L) . #For example, try: #GuessPOL1V([seq(i^2,i=1..10)],x); GuessPOL1V:=proc(L,x) local d,gu,i: for d from 0 to nops(L)-3 do gu:=GuessPOL1V1(L,x,d): if gu<>FAIL then if {seq(expand(subs(x=i,gu)-L[i]),i=1..nops(L))}={0} then RETURN(gu): fi: fi: od: FAIL: end: #GuessPOL(L,VarList): Tries to guess a polynomial in #VarList from the data L that is a list of lists of lists... #such that P(i1,i2,i3, ...)=L[i1][i2][i3] ... #For example, try: #GuessPOL([seq([seq(i+j,i=1..10)],j=1..10)],[x,y]); GuessPOL:=proc(L,VarList) local gu,i,VarList1,mu: if nops(VarList)=1 then RETURN(GuessPOL1V(L,VarList[1])): fi: gu:=[]: VarList1:=[op(2..nops(VarList),VarList)]: for i from 1 to nops(L) do mu:=GuessPOL(L[i],VarList1): if mu=FAIL then RETURN(FAIL): else gu:=[op(gu),mu]: fi: od: GuessPOL1V(gu,VarList[1]): end: #GRpgfOld(x,n,A): the prob. generating function, in x, for "duration of life" #for a fair gambler in [0,A] currently with n dollars. #For example, try: #GRpgfOld(x,10,20); GRpgfOld:=proc(x,n,A) local a,b: a:=(1+sqrt(1-x^2))/x: b:=(1-sqrt(1-x^2))/x: simplify((1-b^A)/(a^A-b^A)*a^n+(a^A-1)/(a^A-b^A)*b^n): end: #GRpgf(x,A,k): the prob. generating function, in x, for "duration of life" #for a fair gambler in [0,k*A] currently with A dollars #For example, try: #GRpgf(x,10,2); GRpgf:=proc(x,A,k) local a,b,i: a:=(1+sqrt(1-x^2))/x: b:=(1-sqrt(1-x^2))/x: (1+add(a^(A*i)*b^((k-2-i)*A),i=0..k-2))/ add(a^(A*i)*b^((k-1-i)*A),i=0..k-1): end: ########From DAVID_IAN #GJst(alphabet,MISTAKES,s,t): given an alphabet alphabet #and a set of words in that alphabet, MISTAKES, and symbols #s and t, uses the Goulden-Jackson method to find the #ratinal function in s whose coeff. of s^n is the #weight-enumerator of all n-letter words in the alphabet #with the weight being #t[m1]^(#of times m1 showed up)*t[m2]^(#of times m2 showed up)*... #For example, try: #GJst({1,2},{[1,2,1]},s,t); GJst:=proc(alphabet,MISTAKES,s,t) local v,eq, var,i,lu,C: eq:={}: var:={}: for i from 1 to nops(MISTAKES) do v:=op(i,MISTAKES): eq:= eq union {findeqzDetail(v,MISTAKES,C,t,s)}: var:=var union {C[op(v)]}: od: var:=solve(eq,var): lu:=1-s*nops(alphabet): for i from 1 to nops(MISTAKES) do v:=op(i,MISTAKES): lu:=lu-subs(var,C[op(v)]): od: lu:=normal(1/lu): for i from 1 to nops(MISTAKES) do v:=op(i,MISTAKES): lu:=subs(t[op(v)]=t[op(v)]-1,lu): od: normal(subs(s=s/nops(alphabet),lu)): end: findeqzDetail:=proc(v,MISTAKES,C,t,s) local eq,i,u: eq:=t[op(v)]*s^(nops(v)): for i from 1 to nops(MISTAKES) do u:=op(i,MISTAKES): eq:=eq+t[op(v)]*overlapz(u,v,s)*C[op(u)]: od: C[op(v)]-eq=0: end: #overlapz is a procedure that given two words u and v, and a variable s #computes the weight-enumerator of all v\suffix(u), #for all suffixes of u that are prefixes of v, but with uniform weight s overlapz:=proc(u,v,s) local i,j,lu,gug: lu:=0: for i from 2 to nops(u) do for j from i to nops(u) while (j-i+1<=nops(v) and op(j,u)=op(j-i+1,v)) do : od: if j-i=nops(v) and u<>v then ERROR(v,`is a subword of`,u,`illegal input`): fi: if j=nops(u)+1 and (i>1 or j>2) then gug:=1: gug:=gug*s^(nops(v)-(j-i)): lu:=lu+gug: fi: od: lu: end: ########End stuff From DAVID_IAN ######### Begin stuff From FELLER#### #Feller(z,t1,t2,t3,t4): The Feller formal-power series #whose coeff. of z^n is the weight-enumerator of all #sequences of {1,-1} of length n with the weight #weight(Sequenec):=t1^(#losing times)*t2^(#returns to the origin) #*t3^(last visit to the origin)*t4^(#sign-changes) #Do: Feller(z,t1,t2,t3,t4); Feller:=proc(z,t1,t2,t3,t4) local G,P1,N1,P,N,psi,alpha,beta: G:=normal( 1+P1+N1+( P*(1+P1+t4*N1) +N*(1+t4*P1+N1) +t4*P*N*(1+t4*P1+N1) +t4*P*N*(1+P1+t4*N1))/(1-t4^2*P*N)): psi:=(1-sqrt(1-4*z^2))/2; alpha:=z/(1-z-psi); beta:=z/(1-z-psi): G:=normal(subs({ P=t2*subs(z=z*t3,psi)/(1-t2*subs(z=z*t3,psi)), N=t2*subs(z=z*t1*t3,psi)/(1-t2*subs(z=z*t1*t3,psi)), P1=alpha, N1=subs(z=t1*z,beta)},G)): normal(subs(z=z/2,G)): end: #ChungFeller(z,t1,t2,t3,t4): The Chung-Feller formal-power series #whose coeff. of z^n is the weight-enumerator of all #sequences of {1,-1} of length n that sum to 0, with the weight #weight(Sequenec):=t1^(#losing times)*t2^(#returns to the origin) #*t3^(last visit to the origin)*t4^(#sign-changes) #Do: ChungFeller(z,t1,t2,t3,t4); ChungFeller:=proc(z,t1,t2,t3,t4) local F,P,N,psi,G: F:=normal(1+(P+N+2*t4*P*N)/(1-t4^2*P*N)): psi:=(1-sqrt(1-4*z^2))/2; G:= normal(subs({P=t2*subs(z=z*t3,psi)/(1-t2*subs(z=z*t3,psi)), N=t2*subs(z=z*t1*t3,psi)/(1-t2*subs(z=z*t1*t3,psi)) },F)): normal(subs(z=z/2,G)): end: ######### End stuff From FELLER#### RS50:=proc(t): [t, 1/2*t+1/2*t^2, 1/6*t+1/2*t^2+1/3*t^3, 1/24*t+7/24*t^2+11/24*t^3+5/24*t^4, 1 /120*t+1/8*t^2+43/120*t^3+3/8*t^4+2/15*t^5, 1/720*t+31/720*t^2+37/180*t^3+67/ 180*t^4+211/720*t^5+61/720*t^6, 1/5040*t+1/80*t^2+2/21*t^3+4/15*t^4+589/1680*t^ 5+53/240*t^6+17/315*t^7, 1/40320*t+127/40320*t^2+503/13440*t^3+2057/13440*t^4+ 4033/13440*t^5+4159/13440*t^6+6551/40320*t^7+277/8064*t^8, 1/362880*t+17/24192* t^2+4661/362880*t^3+599/8064*t^4+24737/120960*t^5+827/2688*t^6+94631/362880*t^7 +2839/24192*t^8+62/2835*t^9, 1/3628800*t+73/518400*t^2+7123/1814400*t^3+57133/ 1814400*t^4+17749/151200*t^5+36619/151200*t^6+535453/1814400*t^7+385123/1814400 *t^8+303271/3628800*t^9+50521/3628800*t^10, 1/39916800*t+31/1209600*t^2+21629/ 19958400*t^3+343/28800*t^4+585967/9979200*t^5+16123/100800*t^6+5252539/19958400 *t^7+162853/604800*t^8+6712403/39916800*t^9+7909/134400*t^10+1382/155925*t^11, 1/479001600*t+2047/479001600*t^2+130807/479001600*t^3+279151/68428800*t^4+ 6270401/239500800*t^5+21989951/239500800*t^6+46991731/239500800*t^7+64207501/ 239500800*t^8+113180237/479001600*t^9+62447027/479001600*t^10+19665491/ 479001600*t^11+540553/95800320*t^12, 1/6227020800*t+1/1520640*t^2+43831/ 691891200*t^3+5837/4561920*t^4+199501/18869760*t^5+106979/2280960*t^6+15792995/ 124540416*t^7+56377/253440*t^8+35932321/138378240*t^9+916711/4561920*t^10+ 206102713/2075673600*t^11+25933/912384*t^12+21844/6081075*t^13, 1/87178291200*t +8191/87178291200*t^2+49481/3632428800*t^3+1346231/3632428800*t^4+22739263/ 5811886080*t^5+125990593/5811886080*t^6+158657089/2179457280*t^7+347887039/ 2179457280*t^8+1378378093/5811886080*t^9+1403731603/5811886080*t^10+75627697/ 454053600*t^11+33834127/454053600*t^12+154945121/7925299200*t^13+199360981/ 87178291200*t^14, 1/1307674368000*t+5461/435891456000*t^2+892709/326918592000*t ^3+3625667/36324288000*t^4+583605751/435891456000*t^5+29637871/3228825600*t^6+ 308647777/8172964800*t^7+55328753/544864320*t^8+48884101807/261534873600*t^9+ 6993962249/29059430400*t^10+23680002257/108972864000*t^11+1636094839/ 12108096000*t^12+72214858159/1307674368000*t^13+5829251611/435891456000*t^14+ 929569/638512875*t^15, 1/20922789888000*t+4681/2988969984000*t^2+195071/ 380414361600*t^3+525879863/20922789888000*t^4+2975820103/6974263296000*t^5+ 5013773861/1394852659200*t^6+15007880821/836911595520*t^7+18765011411/ 321889075200*t^8+109105709723/836911595520*t^9+864245314297/4184557977600*t^10+ 1635670551113/6974263296000*t^11+1326060148151/6974263296000*t^12+34648177031/ 321889075200*t^13+847138153657/20922789888000*t^14+190473830831/20922789888000* t^15+3878302429/4184557977600*t^16, 1/355687428096000*t+257/1394852659200*t^2+ 32219497/355687428096000*t^3+307147/51661209600*t^4+45302657173/355687428096000 *t^5+609774727/464950886400*t^6+254526164999/32335220736000*t^7+42686545229/ 1394852659200*t^8+5829056722823/71137485619200*t^9+72807469793/464950886400*t^ 10+7035623087249/32335220736000*t^11+7900707533/35765452800*t^12+57790948302343 /355687428096000*t^13+23548852891/278970531840*t^14+950967285197/32335220736000 *t^15+2869299781/464950886400*t^16+6404582/10854718875*t^17, 1/6402373705728000 *t+131071/6402373705728000*t^2+48362021/3201186852864000*t^3+4246015451/ 3201186852864000*t^4+114539532959/3201186852864000*t^5+1435829923169/ 3201186852864000*t^6+10304312128577/3201186852864000*t^7+47505797204927/ 3201186852864000*t^8+342261931489/7275424665600*t^9+8551100610389/ 80029671321600*t^10+570085188430847/3201186852864000*t^11+705073145844737/ 3201186852864000*t^12+647102994052289/3201186852864000*t^13+435500327242559/ 3201186852864000*t^14+16080699334727/246245142528000*t^15+9689659411763/ 457312407552000*t^16+26684005437391/6402373705728000*t^17+2404879675441/ 6402373705728000*t^18, 1/121645100408832000*t+73/33874993152000*t^2+48383869/ 20274183401472000*t^3+42690407/152437469184000*t^4+21395152243/2252687044608000 *t^5+1694421301/11725959168000*t^6+2803175677/2271617187840*t^7+48643113599/ 7258927104000*t^8+2776942269809/111396612096000*t^9+92187137527/1385795174400*t ^10+18068794435037/137919614976000*t^11+4213145411347/21776781312000*t^12+ 89170125085309/413758844928000*t^13+9168376316221/50812489728000*t^14+ 2268917222591923/20274183401472000*t^15+7609755023657/152437469184000*t^16+ 40968397162423/2703224453529600*t^17+855344118107/304874938368000*t^18+ 443861162/1856156927625*t^19, 1/2432902008176640000*t+524287/ 2432902008176640000*t^2+96796867/270322445352960000*t^3+99173957/ 1766813368320000*t^4+10769264633/4505374089216000*t^5+991980572707/ 22526870446080000*t^6+12910988859473/28963119144960000*t^7+7467438836077/ 2633010831360000*t^8+237367227106469/19308746096640000*t^9+3781891109257/ 99019210752000*t^10+3695401323221/42066984960000*t^11+981962757392051/ 6436248698880000*t^12+836456227682749/4137588449280000*t^13+847356682421827/ 4137588449280000*t^14+2130265242137729/13516122267648000*t^15+2041854814655273/ 22526870446080000*t^16+3402497099271779/90107481784320000*t^17+970533809451421/ 90107481784320000*t^18+4581126864886571/2432902008176640000*t^19+14814847529501 /97316080327065600*t^20, 1/51090942171709440000*t+1271/61928414753587200*t^2+ 74686849/1459741204905984000*t^3+14737529/1376186994524160*t^4+162379927603/ 283838567620608000*t^5+109562931757/8601168715776000*t^6+7633359210539/ 50089158991872000*t^7+833819705143/737243032780800*t^8+1381776899322931/ 243290200817664000*t^9+50116300046863/2457476775936000*t^10+22006589545419103/ 405483668029440000*t^11+53844286226089/491495355187200*t^12+20686980685346663/ 121645100408832000*t^13+30120438239329/147448606556160*t^14+161546928602968291/ 851515702861824000*t^15+1161826134493823/8601168715776000*t^16+ 11744460972928483/162193467211776000*t^17+194697391650479/6880934972620800*t^18 +4574683392106439/601069907902464000*t^19+390287129333417/309642073767936000*t^ 20+18888466084/194896477400625*t^21, 1/1124000727777607680000*t+42799/ 22938790362808320000*t^2+392158387/56200036388880384000*t^3+109557850957/ 56200036388880384000*t^4+465828947053/3568256278659072000*t^5+87630055506181/ 24977793950613504000*t^6+231752678205757/4683336365740032000*t^7+ 1999240533939907/4683336365740032000*t^8+696327422258959/281704443052032000*t^9 +54574000480780771/5352384417988608000*t^10+69494303123450287/ 2230160174161920000*t^11+161506526942189137/2230160174161920000*t^12+ 696553956773441911/5352384417988608000*t^13+977032379836419721/ 5352384417988608000*t^14+938905094564015617/4683336365740032000*t^15+ 114982080052359241/669048052248576000*t^16+2843327334355420441/ 24977793950613504000*t^17+1427267728778345191/24977793950613504000*t^18+ 24114043507292443/1146939518140416000*t^19+15825473721430903/ 2957896652046336000*t^20+946075012113714451/1124000727777607680000*t^21+ 69348874393137901/1124000727777607680000*t^22, 1/25852016738884976640000*t+ 60787/374666909259202560000*t^2+2941450049/3231502092360622080000*t^3+ 1059482269/3122224243826688000*t^4+29440830752587/1034080669555399065600*t^5+ 7679965195277/8325931316871168000*t^6+412273297807207/26929184103005184000*t^7+ 59741185203781/390278030478336000*t^8+877308170100184019/861733891296165888000* t^9+2854869924510949/594709379776512000*t^10+1288260297648644513/ 76940526008586240000*t^11+2360706655107271/53099051765760000*t^12+ 56387333525875431743/615524208068689920000*t^13+264509622044867497/ 1784128139329536000*t^14+255298027005099457/1346459205150259200*t^15+ 2078917544211271/10841056402176000*t^16+20234235200704648199/ 132574444814794752000*t^17+786813512812958027/8325931316871168000*t^18+ 1517077879119447923/34015811498532864000*t^19+145222302210189707/ 9366672731480064000*t^20+96776543179971529031/25852016738884976640000*t^21+ 70044011000242679/124888969753067520000*t^22+113927491862/2900518163668125*t^23 , 1/620448401733239439360000*t+8388607/620448401733239439360000*t^2+70598995501 /620448401733239439360000*t^3+7022734486007/124089680346647887872000*t^4+ 105409913979101/17727097192378269696000*t^5+1515583685340887/ 6531035807718309888000*t^6+37434283027125737/8272645356443192524800*t^7+ 2163947879461016179/41363226782215962624000*t^8+8238765996308894603/ 20681613391107981312000*t^9+44223001430213061749/20681613391107981312000*t^10+ 6583935698207967871/777504262823608320000*t^11+19847478724606621057/ 777504262823608320000*t^12+126296233708093424477/2110368713378365440000*t^13+ 326997447001724756401/2954516198729711616000*t^14+3368041876058881832653/ 20681613391107981312000*t^15+3957198028380550809907/20681613391107981312000*t^ 16+1483303154279295342031/8272645356443192524800*t^17+5508198174909959164661/ 41363226782215962624000*t^18+80712127254778993403/1042770423081074688000*t^19+ 4276906011310885123411/124089680346647887872000*t^20+7045037713402647386891/ 620448401733239439360000*t^21+1618197631575447463157/620448401733239439360000*t ^22+21036995093331899701/56404400157567221760000*t^23+238685140977801337/ 9545360026665222144000*t^24, 1/15511210043330985984000000*t+ 124453664010157182881/757568256084541440000*t^18+241/222814192966041600000*t^2+ 6418344701/470036667979726848000000*t^3+6161110399226077/ 5170403347776995328000000*t^5+102949064663/11363523841268121600000*t^4+ 127219680744641/2272704768253624320000*t^6+3973415326428899671/ 3102242008666197196800000*t^7+58870668456604/3698160658676859375*t^25+ 4315392183382813/252522752028180480000*t^8+1113359781669034831/ 7493338185184051200000*t^9+9791321053566251/10822403658350592000*t^10+ 3490440127784163133081/861733891296165888000000*t^11+26115640727068674913/ 1893920640211353600000*t^12+13496291480661395491627/369314524841213952000000*t^ 13+16124854743753673231/210435626690150400000*t^14+93907330281523951681/ 11363523841268121600000*t^22+5826696865230506310809/45354415331377152000000*t^ 15+354761545447111192491421/3102242008666197196800000*t^19+9489154115649116857/ 151513651216908288000*t^20+165488085718720789/668442578898124800000*t^24+ 65621268586277236439/378784128042270720000*t^16+64873305979853692893401/ 344693556518466355200000*t^17+136482966873659175371897/ 5170403347776995328000000*t^21+9358556843465453416061/5170403347776995328000000 *t^23, 1/403291461126605635584000000*t+1621587003944651587124201/ 8962032469480125235200000*t^18+33554431/403291461126605635584000000*t^2+ 4604586251/2922401892221779968000000*t^3+7713486764992423/ 33607621760550469632000000*t^5+93719014051843/67215243521100939264000000*t^4+ 435923357719107673/33607621760550469632000000*t^6+14030747737004722793/ 40329146112660563558400000*t^7+66164529094650835995511/ 403291461126605635584000000*t^25+215694800764702012823/ 40329146112660563558400000*t^8+474768290594155759121/8962032469480125235200000* t^9+3269888981133024721391/8962032469480125235200000*t^10+ 20611528671393515644403/11202540586850156544000000*t^11+4152009146866639978847/ 589607399307902976000000*t^12+25185872494635188974709/1200272205733945344000000 *t^13+59485890222226029710219/1200272205733945344000000*t^14+ 673895565993586123856083/33607621760550469632000000*t^22+ 1052268024112442332113533/11202540586850156544000000*t^15+ 5960492574083086104779783/40329146112660563558400000*t^19+ 3904017217328290096358393/40329146112660563558400000*t^20+ 4948919732663793374429/3953837854182408192000000*t^24+70123163019831635586421/ 487066982036963328000000*t^16+1605390517828543543968791/ 8962032469480125235200000*t^17+4087072509293123892361/ 403291461126605635584000000*t^26+1684393785753620180117293/ 33607621760550469632000000*t^21+401953255086217764035923/ 67215243521100939264000000*t^23, 1/10888869450418352160768000000*t+ 1618393280906714557781699/8962032469480125235200000*t^18+22369621/ 3629623150139450720256000000*t^2+190633226237/1088886945041835216076800000*t^3+ 394035829031137/9259242729947578368000000*t^5+124994031969623/ 604937191689908453376000000*t^4+145793585679678263/50411432640825704448000000*t ^6+21507113490168216917/236714553269964177408000000*t^7+ 9394061631649897576726147/10888869450418352160768000000*t^25+ 25363062528114194341/15780970217997611827200000*t^8+39420789553662199980581/ 2177773890083670432153600000*t^9+11345479101661070341441/ 80658292225321127116800000*t^10+240903002786870762871493/ 302468595844954226688000000*t^11+6749460090063058334963/ 1976918927091204096000000*t^12+1033184089963131546958649/ 90740578753486268006400000*t^13+650880397527665896209943/ 21604899703211016192000000*t^14+3002515193687830448895257/ 75617148961238556672000000*t^22+133147193591142889161989/ 2076443449736527872000000*t^15+26509259485340797380515831/ 155555277863119316582400000*t^19+47472656354218781722806293/ 362962315013945072025600000*t^20+868030781551303827961291/ 201645730563302817792000000*t^24+657255131770584305009039/ 5930756781273612288000000*t^16+11129959886717539666134907/ 71169081375283347456000000*t^17+394002679933228635193471/ 3629623150139450720256000000*t^26+55046367898063200855482911/ 680554340651147010048000000*t^21+8374643517010684/1298054391195577640625*t^27+ 99769650387295936806223/6599314818435364945920000*t^23, 1/ 304888344611713860501504000000*t+1119357101274141146034821101/ 6775296546926974677811200000*t^18+19173961/43555477801673408643072000000*t^2+ 5719063895999/304888344611713860501504000000*t^3+386537719875472609/ 50814724101952310083584000000*t^5+9001478311796801/ 304888344611713860501504000000*t^4+31577270117506511647/ 50814724101952310083584000000*t^6+3483488226654834964273/ 152444172305856930250752000000*t^7+940204895551171029049135189/ 304888344611713860501504000000*t^25+14165475208038111134147/ 30488834461171386050150400000*t^8+10349843356972300301629/ 1742219112066936345722880000*t^9+3169202096950866609731753/ 60977668922342772100300800000*t^10+11164122204683851621273117/ 33876482734634873389056000000*t^11+53424744083494156400521699/ 33876482734634873389056000000*t^12+148975278010315841498204051/ 25407362050976155041792000000*t^13+23132607756104369425118591/ 1337229581630323949568000000*t^14+10180983927179001118788611141/ 152444172305856930250752000000*t^22+1044762603280942022671866361/ 25407362050976155041792000000*t^15+2170100063139310175719089089/ 12195533784468554420060160000*t^19+9607274500644009707932552507/ 60977668922342772100300800000*t^20+575233125914127663328781021/ 50814724101952310083584000000*t^24+2172423827308083137010253/ 27290399625108652032000000*t^16+4281603631731839242755459167/ 33876482734634873389056000000*t^17+180747257277752951251546027/ 304888344611713860501504000000*t^26+17382472308083322736107656123/ 152444172305856930250752000000*t^21+21866924546405967976005251/ 304888344611713860501504000000*t^27+226263612816737992108653317/ 7259246300278901440512000000*t^23+13181680435827682794403/ 3209350995912777478963200000*t^28, 1/8841761993739701954543616000000*t+ 315825673742495974240905947/2258432182308991559270400000*t^18+617093/ 20325889640780924033433600000*t^2+17157325905751/ 8841761993739701954543616000000*t^3+828948289369526789/ 631554428124264425324544000000*t^5+27595123218961/6775296546926974677811200000* t^4+436519918018762597/3387648273463487338905600000*t^6+24505820499607771894873 /4420880996869850977271808000000*t^7+74471879835779020049018262421/ 8841761993739701954543616000000*t^25+1315371485921688988733/ 10162944820390462016716800000*t^8+16588964544831257154906611/ 8841761993739701954543616000000*t^9+555485474590450762933/ 30112429097453220790272000*t^10+50239449474364719219896251/ 384424434510421824110592000000*t^11+224615130435508175541823/ 322633168901284508467200000*t^12+2119089730306445664948276179/ 736813499478308496211968000000*t^13+15961714643990822954910367/ 1693824136731743669452800000*t^14+332191271048966698473646093/ 3387648273463487338905600000*t^22+1079152221572638293036959377/ 43341970557547558600704000000*t^15+215006765584280110010132103479/ 1263108856248528850649088000000*t^19+139649714686551470427761101/ 813035585631236961337344000*t^20+82187145209228119775238023/ 3387648273463487338905600000*t^24+30312997393092651876366391/ 564608045577247889817600000*t^16+14777870890622530955878507919/ 155118631469117578149888000000*t^17+6385953565422598253359631/ 2903698520111560576204800000*t^26+633186361824532196496489682363/ 4420880996869850977271808000000*t^21+3591768906684227173751093371/ 8841761993739701954543616000000*t^27+241143615074246710460021756113/ 4420880996869850977271808000000*t^23+106877320322149043860079/ 2258432182308991559270400000*t^28+689005380505609448/263505041412702261046875*t ^29, 1/265252859812191058636308480000000*t+9756564626228325879410795638751/ 88417619937397019545436160000000*t^18+536870911/ 265252859812191058636308480000000*t^2+1608507692273/ 8289151869130970582384640000000*t^3+3415461057773009933/ 15603109400717121096253440000000*t^5+346306989847691/ 637627066856228506337280000000*t^4+6848407483120014246691/ 265252859812191058636308480000000*t^6+43062416325720742238437/ 33156607476523882329538560000000*t^7+4970028825781521419359088873911/ 265252859812191058636308480000000*t^25+164825026315943521858381/ 4736658210931983189934080000000*t^8+151523310198215925777731621/ 265252859812191058636308480000000*t^9+1671665026358831567641079771/ 265252859812191058636308480000000*t^10+1650683053798440026374760731/ 33156607476523882329538560000000*t^11+9765177894998426018849920741/ 33156607476523882329538560000000*t^12+119177569219925873013230597051/ 88417619937397019545436160000000*t^13+432185334329824471962177047621/ 88417619937397019545436160000000*t^14+4858430190318136769638122255023/ 37893265687455865519472640000000*t^22+79008149514015475837245767411/ 5526101246087313721589760000000*t^15+2496276314037308596520100997613/ 16578303738261941164769280000000*t^19+149604910997156232627610456577/ 872542302013786377093120000000*t^20+7820604631095698144889349861/ 177308061371785467002880000000*t^24+188857790346637385485237174541/ 5526101246087313721589760000000*t^16+5958839059000172796514575945761/ 88417619937397019545436160000000*t^17+71851742112639699427580715167/ 11532733035312654723317760000000*t^26+43162991434746207194463715876151/ 265252859812191058636308480000000*t^21+2725727029146723852426530573/ 1745084604027572754186240000000*t^27+2761325912505708363515131500217/ 33156607476523882329538560000000*t^23+441543893249023104553682821/ 265252859812191058636308480000000*t^30+9205418322566705050766647817/ 33156607476523882329538560000000*t^28+8272418863876146191617128571/ 265252859812191058636308480000000*t^29, 1/8222838654177922817725562880000000*t+ 91148542851669743738496881/1115919123817058795520000000*t^18+49981/ 382760259469251166863360000000*t^2+109983814337/5856722688160913687838720000000 *t^3+13831411589006124173/391563745437043943701217280000000*t^5+124966567177/ 1783044686968561336320000000*t^4+5746718809259100583/ 1148280778407753500590080000000*t^6+3086473587593157026557/ 10488314609920819920568320000000*t^7+2441887363638461077862008850857/ 69099484488890107711979520000000*t^25+24034014265718935409/ 2658057357425355325440000000*t^8+7306155967521527967172097/ 43507082826338215966801920000000*t^9+340559738388648313705549/ 164040111201107642941440000000*t^10+1120419191482192031123693/ 61373627811448893997056000000*t^11+34356090032028105340504541/ 287070194601938375147520000000*t^12+5875974553618015749322026307/ 9708192035629188686807040000000*t^13+18177165498293843639398883/ 7505103126848062095360000000*t^14+174020411978088844701392309893/ 1148280778407753500590080000000*t^22+563451175674249750476801383/ 72099805820162763202560000000*t^15+3313416937165693279364972571673/ 26697528097980268888719360000000*t^19+15137255113712313964285502897/ 95690064867312791715840000000*t^20+716871255302433409392596699/ 10252506950069227683840000000*t^24+61669426919948996571065801/ 2990314527103524741120000000*t^16+3735877107027662118468199092401/ 83058976304827503209349120000000*t^17+1833001563683766206234063537/ 127586753156417055621120000000*t^26+95598663117114832597407365339/ 565027049692704103464960000000*t^21+49853330588778405917973918373/ 10876770706584553991700480000000*t^27+2769523518291646312572691948543/ 24472734089815246481326080000000*t^23+23572884907472995411324393/ 1148280778407753500590080000000*t^30+2218379002973589534677929/ 2007483878335233392640000000*t^28+30519118018114665817049565733/ 161232130474076917994618880000000*t^29+129848163681107301953/ 122529844256906551386796875*t^31, 1/263130836933693530167218012160000000*t+ 151546816612822564581438158715839/2657887241754480102699171840000000*t^18+ 2147483647/263130836933693530167218012160000000*t^2+154417633243099/ 87710278977897843389072670720000000*t^3+9883529570392771091/ 1790005693426486599776993280000000*t^5+768459907318075429/ 87710278977897843389072670720000000*t^4+11772359060599621040507/ 12530039853985406198438952960000000*t^6+2426753036118849646402591/ 37590119561956218595316858880000000*t^7+44559906153425913314976196515089/ 767145297182779971332997120000000*t^25+5027330964919726035039953/ 2211183503644483446783344640000000*t^8+85443516059034117407100791/ 1790005693426486599776993280000000*t^9+8276347829363495519554327487/ 12530039853985406198438952960000000*t^10+1050526673186015620864623317/ 162727790311498781797908480000000*t^11+53323635826554562490278685741/ 1139094532180491472585359360000000*t^12+892035414522870990022821967553/ 3417283596541474417756078080000000*t^13+3928464994689296237819322272063/ 3417283596541474417756078080000000*t^14+9823488666881348172728552485153/ 59952343798973235399229440000000*t^22+10842503087477828081191947387679/ 2657887241754480102699171840000000*t^15+4256988535796972036762447209499/ 44380306448590576853975040000000*t^19+42228322659573884896851862111043/ 310662145140134037977825280000000*t^20+72674761660804346347427373253841/ 737061167881494482261114880000000*t^24+31475809528116294421148920659169/ 2657887241754480102699171840000000*t^16+75608372891626530743505711350849/ 2657887241754480102699171840000000*t^17+45880580527492279674188689729177/ 1634353024432879069361602560000000*t^26+185122761200894996790923712567949/ 1139094532180491472585359360000000*t^21+137107089494945765686511290582747/ 12530039853985406198438952960000000*t^27+248717764373922846995708350933801/ 1790005693426486599776993280000000*t^23+11293563781580186824572769953659/ 87710278977897843389072670720000000*t^30+2088463430347521052196056349/ 3095656905102276825496682496000000*t^32+3233669862460580265647134572409/ 963849219537338938341457920000000*t^28+68326977245435195514549141477509/ 87710278977897843389072670720000000*t^29+3551888312140990281660844836191/ 263130836933693530167218012160000000*t^31, 1/ 8683317618811886495518194401280000000*t+4322959508983390070519967710827/ 114653959448232475018395648000000*t^18+16843009/ 34052225956125045080463507456000000*t^2+47922787816261/ 299424745476271948121317048320000000*t^3+484447536552982208519/ 578887841254125766367879626752000000*t^5+12054878699530907/ 11350741985375015026821169152000000*t^4+30808894514256970781/ 180170507704365317886050304000000*t^6+3405638458256031817829651/ 248094789108911042729091268608000000*t^7+21059991069863836711350703018683661/ 248094789108911042729091268608000000*t^25+2694360816752660238639199/ 4864603708017863582923358208000000*t^8+22475804120640075549927991/ 1710998545578696846407525990400000*t^9+329589411440303046110425153/ 1621534569339287860974452736000000*t^10+911828846103423219294433504337/ 413491315181518404548485447680000000*t^11+867699362470174967671248689/ 49137411192099632150740992000000*t^12+249381909370449376460588173841/ 2301435891548339913998991360000000*t^13+33103621002670191635243728247/ 63176671532699527050952704000000*t^14+24115244003837448895893884644109/ 147412233576298896452222976000000*t^22+321975576380774741259736694985131/ 157878502160216118100330807296000000*t^15+2206674851399271014403617453331457/ 31575700432043223620066161459200000*t^19+48415144319189912496486293946847/ 442236700728896689356668928000000*t^20+41060949528887980198192236289/ 326988217249301847343104000000*t^24+2233152805737210533917824273503/ 343961878344697425055186944000000*t^16+300228140234578078888201823842829/ 17542055795579568677814534144000000*t^17+232436569861044809048427601470017/ 4864603708017863582923358208000000*t^26+16388491029244466540492774031710759/ 112770358685868655785950576640000000*t^21+34093137157575177292112533920061/ 1540961423036714551112368128000000*t^27+64496523300888074746111081190460767/ 413491315181518404548485447680000000*t^23+2070836236115675839146783299927/ 3783580661791671675607056384000000*t^30+60387365907582249796918225747/ 6810445191225009016092701491200000*t^32+1736640792209901647222/ 4043484860477916195764296875*t^33+83416627053836118557903982523/ 10071643287821663732760576000000*t^28+11319191087879416436997935230933/ 4631102730033006130943037014016000*t^29+39943790748948843096639303846341/ 457016716779572973448326021120000000*t^31, 1/ 295232799039604140847618609643520000000*t+441931315954656120812717079518359/ 18638434282803291720177942528000000*t^18+53353631/ 1833744093413690315823718072320000000*t^2+2084643417491017/ 147616399519802070423809304821760000000*t^3+242285526110028046199/ 1968218660264027605650790730956800000*t^5+61687823937383333/ 493700332842147392721770250240000000*t^4+186343619168650617911/ 6169964452238331052196836147200000*t^6+411857007344430868277159/ 145434876374189231944639709184000000*t^7+4107099859152734594835947232490001/ 36674881868273806316474361446400000*t^25+42455780599326421888434217/ 324431647296268286645734735872000000*t^8+453535113994235708026549/ 129434138863019417719642521600000*t^9+51096610114251611330795502151/ 843522282970297545278910313267200000*t^10+465660103225211142144841332251/ 639032032553255716120386600960000000*t^11+45232576183219534410087384820951/ 7029352358085812877324252610560000000*t^12+11879901384461845976569055964571/ 273870871094252449765879971840000000*t^13+62949855227866442212730387141221/ 273870871094252449765879971840000000*t^14+291372125515971645692119458512931517/ 1917096097659767148361159802880000000*t^22+47768287038668640657240467965237/ 48798809758612254685556794982400000*t^15+25881967832058011785535071623766457/ 536786907344734801541124744806400000*t^19+8920164667081285044070403182285531/ 107357381468946960308224948961280000*t^20+1027367209583738091251134463199104041 /7029352358085812877324252610560000000*t^24+261343468722037603876932080292671/ 76683843906390685934446392115200000*t^16+183178760307561705026522294788969/ 18638434282803291720177942528000000*t^17+3208722853602371413031946896492227/ 44395909630015660277837384908800000*t^26+21219905272067762534415142376770217/ 174281463423615195305559982080000000*t^21+23466932917557376760077758341067523/ 602515916407355389484935938048000000*t^27+162296921474139020792515662228756913/ 1004193194012258982474893230080000000*t^23+498467493937465315543881618112913/ 281174094323432515092970104422400000*t^30+8740222019457905981538006701777737/ 147616399519802070423809304821760000000*t^32+156175308565000080182689281967421/ 26839345367236740077056237240320000000*t^33+948073377927549762406172325504023/ 54774174218850489953175994368000000*t^28+12267253631991067278689181387142409/ 1968218660264027605650790730956800000*t^29+56586697536055983956123465935412407/ 147616399519802070423809304821760000000*t^31+80723299235887898062168247453281/ 295232799039604140847618609643520000000*t^34, 1/ 10333147966386144929666651337523200000000*t+24750186995067430390978686873066311 /1739587199728307227216607969280000000*t^18+5726623061/ 3444382655462048309888883779174400000000*t^2+893419221062909/ 738081997599010352119046524108800000000*t^3+90875581574563737162719/ 5166573983193072464833325668761600000000*t^5+2732619341109885583/ 191354591970113794993826876620800000000*t^4+594852250797750369624133/ 114812755182068276996296125972480000000*t^6+20215374339931590318229091/ 35631544711676361826436728750080000000*t^7+ 99612250964597887877591262831483923401/738081997599010352119046524108800000000* t^25+1476970550804608114801337573/49205466506600690141269768273920000000*t^8+ 133880384390061467768816552803/147616399519802070423809304821760000000*t^9+ 44313405725243179977776251/2531536065575998875406172160000000*t^10+ 172159337934795706200102928556513/738081997599010352119046524108800000000*t^11+ 186044165785654019250266990978791/82009110844334483568782947123200000000*t^12+ 12380915032743901682663821746176933/738081997599010352119046524108800000000*t^ 13+2174546859707153821714695213505787/22366121139363950064213531033600000000*t^ 14+328653588374608994597140572185198093/2485124571040438896023725670400000000*t ^22+11176481340278149362601593574424527/24720449680349629018341271142400000000* t^15+156360350461889087101462612702729117/4944089936069925803668254228480000000 *t^19+1869233974886820094438459033490632307/ 31312569595109530089898943447040000000*t^20+ 12871088334343813248676030975649594041/82009110844334483568782947123200000000*t ^24+779312898048935536622827761264073/453805356450862754926071644160000000*t^16 +36249034756685078537673963534093733/6709836341809185019264059310080000000*t^17 +4854021729900950138700616515088925017/49205466506600690141269768273920000000*t ^26+45133345345340487068233737479935553639/ 469688543926642951348484151705600000000*t^21+ 8990927346733649797547936667255674447/147616399519802070423809304821760000000*t ^27+114973417445647425612033397858791387007/ 738081997599010352119046524108800000000*t^23+ 535583412733397520529458457615869163/114812755182068276996296125972480000000*t^ 30+41910773922042004182045184194830827/156562847975547650449494717235200000000* t^32+17986620272362901948092739929139293/ 449267302886354127376810927718400000000*t^33+ 516327122441389022108059053723137641/16401822168866896713756589424640000000*t^ 28+13896424871341243816873391897839776701/ 1033314796638614492966665133752320000000*t^29+ 945157595704507860401258330526495331/738081997599010352119046524108800000000*t^ 31+4380378738451868149141703258018237/1148127551820682769962961259724800000000* t^34+418781231495293038913922/2405873491984360136479756640625*t^35, 1/ 371993326789901217467999448150835200000000*t+ 27629705873825596905096300715149365663/3381757516271829249709085892280320000000 *t^18+1108378657/11999784735158103789290304779059200000000*t^2+1072103556128897 /10628380765425749070514269947166720000000*t^3+579658021327722156653/ 237240642085396184609693525606400000000*t^5+590258286115765849693/ 371993326789901217467999448150835200000000*t^4+1922287924633441058731/ 2224840471231466611650714402816000000*t^6+513702489344092496293762373/ 4649916584873765218349993101885440000000*t^7+ 1992913539235905311761133311936720537739/ 13285475956782186338142837433958400000000*t^25+2393005334889330694733896679/ 357685891144135786026922546298880000000*t^8+86549953899028507052904825893/ 379585027336633895375509640970240000000*t^9+686660658196099615839706174127/ 139847115334549329875187762462720000000*t^10+962038481129842017174749460423169/ 13285475956782186338142837433958400000000*t^11+ 10293962058970505834344725912886783/13285475956782186338142837433958400000000*t ^12+4168952925704417624917113897611639/664273797839109316907141871697920000000* t^13+131941588316884918662457284689807533/ 3321368989195546584535709358489600000000*t^14+ 57257894067390600245368672473793171669/528399611917473320267044670668800000000* t^22+212768467280721088884672267040068049/ 1056799223834946640534089341337600000000*t^15+ 9559621910451463068033527298203201821/483108216610261321387012270325760000000*t ^19+12531414255768818129104887836387696831/ 307432501479257204519007808389120000000*t^20+ 520507500297284551933201979725769192117/ 3321368989195546584535709358489600000000*t^24+ 9261019672556915894565551706317873/11124202356157333058253572014080000000*t^16+ 1376857231978331099943578539907462231/483108216610261321387012270325760000000*t ^17+326691847072600595358975877030547449009/ 2657095191356437267628567486791680000000*t^26+634185543684887724959185841096293 /5722974258306172576430760740782080000000*t^36+ 313405602938410109050544561954843587/4366938941466721655099542732800000000*t^21 +228270802740353585352539332612068915961/ 2657095191356437267628567486791680000000*t^27+ 579958433943267139126282783423880111/4125924210180803210603365662720000000*t^23 +48275525913251391638685152268253518799/ 4649916584873765218349993101885440000000*t^30+ 5355629634269448548693021118343668197/5812395731092206522937491377356800000000* t^32+1981336606429447286773938236358837043/ 10628380765425749070514269947166720000000*t^33+ 944629191990522157740193930464641033/18581085254240820053346625781760000000*t^ 28+117377107674973495968102567348528620593/ 4649916584873765218349993101885440000000*t^29+ 2884344633844278609987055559161414733/830342247298886646133927339622400000000*t ^31+10052421396794935226488440685900288007/ 371993326789901217467999448150835200000000*t^34+ 929760750445541410257052628280893531/371993326789901217467999448150835200000000 *t^35, 1/13763753091226345046315979581580902400000000*t+ 2934986666051068995932122494523489/651967903657572633450758799360000000*t^18+ 7957/1593699320051844215101854842880000000*t^2+1210439868314759/ 147997345066949946734580425608396800000000*t^3+136048562732730627481/ 411989735728757933618174676172800000000*t^5+2460486631605791/ 14343293880466597935916693585920000000*t^4+15707855902383180233/ 112056983441145296374349168640000000*t^6+6000478678974816816419111273/ 286744856067215521798249574616268800000000*t^7+ 1953382782742357267835572770037043888677/ 12604169497460022936186794488627200000000*t^25+4762135414207939974047/ 3283721126480448245402173440000000*t^8+12784185145216577178675755760463/ 229395884853772417438599659693015040000000*t^9+4798815327869902153533178589/ 3585823470116649483979173396480000000*t^10+1192160640154608176843556801296201/ 54618067822326766056809442784051200000000*t^11+31741100730234143195936962007/ 123649085176436189102730117120000000*t^12+341353078379232337833383112271823/ 150049636874524082573652315340800000000*t^13+4692443309914126234100633209621/ 298818622509720790331597783040000000*t^14+44450831684654391253767193368053/ 529194726994782981697044480000000*t^22+238760232134769903687940146806501533/ 2757162077569380017290861294387200000000*t^15+ 114013956293614478557388567475682129583/ 9625002161696744787633552154951680000000*t^19+640926871612534256672528211258377 /24146959394724912350028103680000000*t^20+131288917273148281102862810167184471/ 896455867529162370994793349120000000*t^24+2609229315321892398651502067/ 6706384789104392625192960000000*t^16+88888594944072679289853497550251111/ 61462338197297220866114637004800000000*t^17+33773310317415540820766527994915471 /239054898007776632265278226432000000*t^26+3352396095117218222217025913201/ 2049041982923799705130956226560000000*t^36+ 51059037190567160752722110569077676531/1002604391843410915378495016140800000000 *t^21+369878003652419696794286773599224264047/ 3343963336060822411641394456166400000000*t^27+ 877892577955929532359270003865065728183/ 7352432206851680046108963451699200000000*t^23+ 18011716023101103095183516922172091/896455867529162370994793349120000000*t^30+ 384632409216331015423676798866903/149409311254860395165798891520000000*t^32+ 131770087780151310285670383581014790597/ 199474682481541232555304051906969600000000*t^33+ 265436939059576549532554057719588799/3585823470116649483979173396480000000*t^28 +803954969003380946927994719205132426899/ 19116323737814368119883304974417920000000*t^29+ 380732641131509063919209991327326516833/ 47790809344535920299708262436044800000000*t^31+ 371442180313369136188310324228273/2868658776093319587183338717184000000*t^34+ 9282755234721512951731436944200350041/ 509768633008383149863554799317811200000000*t^35+56518638202982204522669764/ 801155872830791925447758961328125*t^37, 1/ 523022617466601111760007224100074291200000000*t+ 484714319080329738969658275783857897647/ 203194490080264612183374989937868800000000*t^18+137438953471/ 523022617466601111760007224100074291200000000*t^2+3126970296104579/ 4842802013579639923703770593519206400000000*t^3+2525324090772827042745427/ 58113624162955679084445247122230476800000000*t^5+262344312509332528247/ 14528406040738919771111311780557619200000000*t^4+1286758361547228262517437357/ 58113624162955679084445247122230476800000000*t^6+42046575973038627459030864103/ 10896304530554189828333483835418214400000000*t^7+ 71749913616361705304818050454638288355577/ 478958440903480871575098190567833600000000*t^25+3334003658500492354250855132953 /10896304530554189828333483835418214400000000*t^8+ 4279050129936007489395550461073/322853467571975994913584706234613760000000*t^9+ 1030792636840409943305524050198247/2905681208147783954222262356111523840000000* t^10+1659011112267225689929857839733677/ 259435822156052138769844853224243200000000*t^11+ 21426651348813931111780277964432467/259435822156052138769844853224243200000000* t^12+383146054533985779361793533630791467/ 478958440903480871575098190567833600000000*t^13+ 125241022153380685861132867823078899/20824280039281777025004269155123200000000* t^14+12566455529928188490797773539523928931757/ 203194490080264612183374989937868800000000*t^22+ 3355766024401736817599903497471483291/93130807953454613917380203721523200000000 *t^15+621998364273605663928287465187999302459/ 91437520536119075482518745472040960000000*t^19+ 1511687687548374130896972370384510958789/ 91437520536119075482518745472040960000000*t^20+ 35984549738202240568533822542367053230373/ 279392423860363841752140611164569600000000*t^24+ 49061573048264401064443619656717821743/ 279392423860363841752140611164569600000000*t^16+ 143733760389132537131901710932493238097/ 203194490080264612183374989937868800000000*t^17+ 72248513530557494221964632114224894918407/ 478958440903480871575098190567833600000000*t^26+ 177994855399718957287029358900492736627/ 14528406040738919771111311780557619200000000*t^36+ 101459156222665281715182754612110613943/ 2944847682322675538889492607795200000000*t^21+ 630575142102062916441562230747274293733/ 4804367076963928495737867652300800000000*t^27+ 864517401149233752629303447140253297093/ 9012658834205285217810987456921600000000*t^23+ 100248752149222270329777226202790899524597/ 2905681208147783954222262356111523840000000*t^30+ 66235464217098185062012310249893880509243/ 10896304530554189828333483835418214400000000*t^32+ 36791824521331812274238248505823138317339/ 19371208054318559694815082374076825600000000*t^33+ 411545532053490856818316222148072162311/ 4184448744452453851126529890713600000000*t^28+ 183497812340398257378356323235649727978507/ 2905681208147783954222262356111523840000000*t^29+ 173052185076595633861873300331809199246533/ 10896304530554189828333483835418214400000000*t^31+ 27427751968487931680938760262610026001327/ 58113624162955679084445247122230476800000000*t^34+ 23489580527043108252017828576198947741/ 523022617466601111760007224100074291200000000*t^38+ 1303414394978406351650953944681071142797/ 14528406040738919771111311780557619200000000*t^35+ 559713744390266935691173813049295342691/ 523022617466601111760007224100074291200000000*t^37, 1/ 20397882081197443358640281739902897356800000000*t+ 3217858890869099418538949529849806847877/ 2641528371043439958383874869192294400000000*t^18+91625968981/ 6799294027065814452880093913300965785600000000*t^2+63321157086052319/ 1274867630074840209915017608743931084800000000*t^3+33495148924248532911839/ 6011754223754035767356404874713497600000000*t^5+29149628641268260561/ 15739106544133829752037254428937420800000000*t^4+2574370390733547890457070387/ 755477114118423828097788212588996198400000000*t^6+5356053171761311897813910101/ 7726470485302061878272833992387461120000000*t^7+ 59572095579044982123200929548871968223309/ 437309880824917317525089652257587200000000*t^25+8907265464715420277013658850623 /141651958897204467768335289860436787200000000*t^8+ 745370813989558265548216722554897/242831929538064801888574782617891635200000000 *t^9+1152567760982391937699311921910639/ 12591285235307063801629803543149936640000000*t^10+ 23470203835463171920158664610878103/ 12877450808836769797121389987312435200000000*t^11+ 174237205926756949259704669025612599/ 6745331376057355608015966183830323200000000*t^12+ 66340962049364293305900381734897559827/ 242831929538064801888574782617891635200000000*t^13+ 13913078150888545106670676771651651787/ 6226459731745251330476276477381836800000000*t^14+ 3474042534744636530017954775284267749619/ 80046314274043635102541662702796800000000*t^22+ 43166545847433629166270548834988754241/ 2971719417423869953181859227841331200000000*t^15+ 33485936001154896182675698077466959751697/ 8915158252271609859545577683523993600000000*t^19+ 5869696831160675988765776559155489926177/ 594343883484773990636371845568266240000000*t^20+ 388963751603633388924623721915894522079163/ 3632101510184729942777827945139404800000000*t^24+ 10313334271343508744829412449863345719/ 134522278154989997880660294264422400000000*t^16+ 5824617724524731675099263126446847374943/ 17434087248886703725333574136669143040000000*t^17+ 937946819709236776240236742147098277921037/ 6226459731745251330476276477381836800000000*t^26+ 50487041893471581170866050470846028251/ 814091717800025676829513160117452800000000*t^36+ 1589721988474715573758386987075664470540863/ 71321266018172878876364621468191948800000000*t^21+ 1756502572474888059676146816131370822616847/ 12141596476903240094428739130894581760000000*t^27+ 113902959995098027292273668824693415653113/ 1556614932936312832619069119345459200000000*t^23+ 2016244183798835625672686985305196602934667/ 37773855705921191404889410629449809920000000*t^30+ 56989780693552749093974021660345595282583/ 4569418028942079605430170640659251200000000*t^32+ 2851707639898893654451887855217309498546297/ 618117638824164950261826719390996889600000000*t^33+ 270690986503776719414935873522637311689583/ 2248443792019118536005322061276774400000000*t^28+ 4455592763957317129749726076004059019196529/ 51509803235347079188485559949249740800000000*t^29+ 917969363417114565382635601365363797625721/ 32688913591662569485000451506254643200000000*t^31+ 117127922053323223554038019349733869455293/ 83941901568713758677532023620999577600000000*t^34+32207686319158956594455462/ 1126482925555250126673224649609375*t^39+ 4755982464399223324644655046593444131331/ 6799294027065814452880093913300965785600000000*t^38+ 23805722971420062657051438753290845648531/ 70825979448602233884167644930218393600000000*t^35+ 436067475310383493665479257217655453559/ 52981511899214138593870861662085447680000000*t^37, 1/ 815915283247897734345611269596115894272000000000*t+ 2460917983289039209542385310829357002479/ 4097317802323549641676515660791808000000000*t^18+78536544841/ 116559326178271104906515895656587984896000000000*t^2+434202259286916899/ 116559326178271104906515895656587984896000000000*t^3+12628511464257376135137143 /18131450738842171874346917102135908761600000000*t^5+13737517093183502012929/ 74174116658899794031419206326919626752000000000*t^4+ 1495207794144408248940911539/2924427538522930947475309210021920768000000000*t^6 +33017744894023016274632391712549/ 271971761082632578115203756532038631424000000000*t^7+ 122477937548900758410605842653435623180373703/ 1046045234933202223520014448200148582400000000*t^25+ 3425867853247919204566786989433051/ 271971761082632578115203756532038631424000000000*t^8+ 4282086825567995946238402486752457/ 6181176388241649502618267193909968896000000000*t^9+ 312899515155406032019191538545935737/ 13598588054131628905760187826601931571200000000*t^10+ 395297246985781100866953030301028393/ 781528049088024649756332633712754688000000000*t^11+ 178167893471013405896966889678789520219/ 22664313423552714842933646377669885952000000000*t^12+ 30382665988973228204272348467349447231/ 334940592466296278466999700162609152000000000*t^13+ 7832240196595317676228522837349594881757/ 9713277181522592075542991304715665408000000000*t^14+ 83070812529585422089556307047072397930256061/ 2852850640726915155054584858727677952000000000*t^22+ 5942173673833190960652250357517165664833/ 1046045234933202223520014448200148582400000000*t^15+ 5700960057436235754826191426461138930895493/ 2852850640726915155054584858727677952000000000*t^19+ 3235937174783528182562170742003550396716927/ 570570128145383031010916971745535590400000000*t^20+ 294736636129145777629172986032575724842868763/ 3486817449777340745066714827333828608000000000*t^24+ 835451569627772482088202127458775246289/ 25764660958945867574384592320200704000000000*t^16+ 10860506029980123693533495742761746079219/ 71159539791374300919728874027220992000000000*t^17+ 67270131455030887227391195362001038196945783/ 475475106787819192509097476454612992000000000*t^26+ 21634065899015372936089657694709451772645771/ 90657253694210859371734585510679543808000000000*t^36+ 39446275137893091556806131608082326232770883/ 2852850640726915155054584858727677952000000000*t^21+ 1448204786136595477524252539967607758020021989/ 9713277181522592075542991304715665408000000000*t^27+ 185464128540466910956322551912154305292443621/ 3486817449777340745066714827333828608000000000*t^23+ 48762194791470731034233096332232392815293887/ 647551812101506138369532753647711027200000000*t^30+ 67074767422272622217913131920189456865606333/ 2956214794376441066469606049261289472000000000*t^32+ 2647471687681383524014574276346442953001989059/ 271971761082632578115203756532038631424000000000*t^33+ 1329238632229501606258089046123661473907601691/ 9713277181522592075542991304715665408000000000*t^28+ 2473189660887273211101368675406567090334901987/ 22664313423552714842933646377669885952000000000*t^29+ 39491470314353908251551806608065451496617881/ 883025198320235643231181027701424128000000000*t^31+ 947544637077785422429624881310717201427270717/ 271971761082632578115203756532038631424000000000*t^34+ 594046028724599200715086271256233067377/ 32636611329915909373824450783844635770880000000*t^40+ 372784833053968051810137467220769308961271/ 815915283247897734345611269596115894272000000000*t^39+ 145325094967241319920076111972032394585907/ 26319847846706378527277782890197286912000000000*t^38+ 500337468655166797890800922923834435511709/ 490039209157896537144511273030700236800000000*t^35+ 34900977370555148986005229223454074206912723/ 815915283247897734345611269596115894272000000000*t^37, 1/ 33452526613163807108170062053440751665152000000000*t+ 1210408075732196118095369696752860846173/ 4226445393669503933414199790707671040000000*t^18+61681/ 1876637992635035902123193075949895680000000*t^2+1823649598956213761/ 6690505322632761421634012410688150333030400000000*t^3+ 568313260920863960847895247/6690505322632761421634012410688150333030400000000*t ^5+1191335386128169141/65932548141244261361261516735039668224000000*t^4+ 41119349108489774021843873/549437901177035511343845972791997235200000000*t^6+ 2011042596100306774017819076183/96963845255547266980203078415770294681600000000 *t^7+585176640731091991347306139556807179136237953/ 6126836376037327309188656053743727411200000000*t^25+ 3996910848422583825280377409/1623954387715375895597574303818711040000000*t^8+ 85011091817983853672601852033496007/ 557542110219396785136167700890679194419200000000*t^9+ 774549675076498295639553271769723/137359475294258877835961493197999308800000000 *t^10+54474879221241175758205328991256634089/ 398244364442426275097262643493342281728000000000*t^11+ 12801911707745861488778736975322991/5494379011770355113438459727919972352000000 *t^12+16342338835891703033250590153964210887761/ 557542110219396785136167700890679194419200000000*t^13+ 3330982255363394884135604449507551709/ 11773669310936475243082413702685655040000000*t^14+ 1474926605912993112682093259165466592853/ 78753640875829265840016145168465920000000*t^22+ 1202994481129818134831432100078045201531517/ 557542110219396785136167700890679194419200000000*t^15+ 264101387100054483737613876908970494464132603/ 257327127793567746985923554257236551270400000000*t^19+ 20577012571001432947732745938777699420609/ 6558277335004402655297896226960179200000000*t^20+ 11679939935461737022432228505709777892589/ 183758495376934953626704338726420480000000*t^24+ 15619112147340644516125707380014684277/ 1174012609352639981503944386307686400000000*t^16+ 5793028264929500495606215702857831523633973/ 85775709264522582328641184752412183756800000000*t^17+ 3980092616553264774898309349827242045300401/ 31698340452521279500606498430307532800000000*t^26+ 5310907129363953366099913823771866342691/ 7135557158143318329140856789506457600000000*t^36+ 964271624975833514457089637435430520420578531/ 116966876269803521357237979207834796032000000000*t^21+ 2777560844745270355696147362605129291922839747/ 19225590007565406384005782789333765324800000000*t^27+ 9509825009904349269263273825884588751083552153/ 257327127793567746985923554257236551270400000000*t^23+ 363125690235343224355376892653814749269117/ 3712418251196185887458418735081062400000000*t^30+ 3063065087007898550589750342968942866944107/ 82415685176555326701576895918799585280000000*t^32+ 40601560001544754241966113434729495805865237543/ 2230168440877587140544670803562716777676800000000*t^33+ 570644435424503318056512714510593012582009/ 3924556436978825081027471234228551680000000*t^28+ 71251761292245893812866535844057854631262337141/ 557542110219396785136167700890679194419200000000*t^29+ 181001428902747221721393803776753347949084505653/ 2787710551096983925680838504453395972096000000000*t^31+ 276688899305637637641474535055170462526453/ 36629193411802367422923064852799815680000000*t^34+ 1410211493828985228276049834684/121699582862361447435141825020548828125*t^41+ 28917544888620868911699325190913469967/ 96959629619476854943031642257411276800000000*t^40+ 3535617672707276625603811667517292282925453/ 955786474661823060233430344384021476147200000000*t^39+ 29118117120665745736913010085322941577069/ 988988222118663920418922751025595023360000000*t^38+ 5839780355894817027912090548706260356451228291/ 2230168440877587140544670803562716777676800000000*t^35+ 1130270436471188492269583616772221497448912023/ 6690505322632761421634012410688150333030400000000*t^37, 1/ 1405006117752879898543142606244511569936384000000000*t+ 476496291703383118637955130438611248805246953/ 3602579789109948457802929759601311717785600000000*t^18+2199023255551/ 1405006117752879898543142606244511569936384000000000*t^2+210421116029652571/ 10807739367329845373408789278803935153356800000000*t^3+ 177605454041331995033698883/17562576471910998731789282578056394624204800000000* t^5+241782428443799500621153/ 140500611775287989854314260624451156993638400000000*t^4+ 380937458263622985004495121/35623887366959429476246009286118447513600000000*t^6 +161975790584655900488576855144467/ 46833537258429329951438086874817052331212800000000*t^7+ 267322795078082799503152819583227783877776752833/ 3602579789109948457802929759601311717785600000000*t^25+ 21981798534358671345128075119855597/ 46833537258429329951438086874817052331212800000000*t^8+ 3067460878441153664365287023617324013/ 93667074516858659902876173749634104662425600000000*t^9+ 126277064556250486939121444337164667923/ 93667074516858659902876173749634104662425600000000*t^10+ 2111943480191062741439715681151343357129/ 58541921573036662439297608593521315414016000000000*t^11+ 1271178510362667248789466518070328863209/ 1888449083001182659332180922371655335936000000000*t^12+ 1459712322867755271004568365903637977597/ 158221409656855844430534077279787338956800000000*t^13+ 565589197599980859227818279936635470560879/ 5854192157303666243929760859352131541401600000000*t^14+ 7108029516266576188845992097023024852632562837/ 614076100416468487125499390841132679168000000000*t^22+ 9330120684527486636070369787094053019504879/ 11708384314607332487859521718704263082803200000000*t^15+ 2754360042833456993674273887931500598569704309/ 5403869683664922686704394639401967576678400000000*t^19+ 9053353616479839235048425178300867394017715339/ 5403869683664922686704394639401967576678400000000*t^20+ 246811155411617340203078812083941116140594430949/ 5403869683664922686704394639401967576678400000000*t^24+ 62071474942919925012515151452704737959292689/ 11708384314607332487859521718704263082803200000000*t^16+ 104434622000029783134268571781493949570549783/ 3602579789109948457802929759601311717785600000000*t^17+ 380900222067623401953880210881882299298129002303/ 3602579789109948457802929759601311717785600000000*t^26+ 91759566142520890385477541600288665418642550147/ 46833537258429329951438086874817052331212800000000*t^36+ 100248745753003703437921701202809681649310503/ 21175037945395465073293082442797678592000000000*t^21+ 10364622733519612119397957304745185976310201/ 1405006117752879898543142606244511569936384000000000*t^42+ 1548426589686286001018723743634319424366367367329/ 11708384314607332487859521718704263082803200000000*t^27+ 132986621018462191629945209178180675219120684059/ 5403869683664922686704394639401967576678400000000*t^23+ 690300687433062317037463449039302455409449175269/ 5854192157303666243929760859352131541401600000000*t^30+ 3247952441796948471682367903912434717025497417529/ 58541921573036662439297608593521315414016000000000*t^32+ 2870749452527388398534990817840839739165736419563/ 93667074516858659902876173749634104662425600000000*t^33+ 1700278924053120620416020563430873125898280345439/ 11708384314607332487859521718704263082803200000000*t^28+ 62981911195390245170489689374659762945379688903/ 450322473638743557225366219950163964723200000000*t^29+ 5079413250601495129726008953534982183372389447879/ 58541921573036662439297608593521315414016000000000*t^31+ 1359351354855055784557833932233743664392818147093/ 93667074516858659902876173749634104662425600000000*t^34+ 273363299173912297184785358991042651206919751/ 1405006117752879898543142606244511569936384000000000*t^41+ 347758702335421239405217881400097116006204783/ 140500611775287989854314260624451156993638400000000*t^40+ 2841644716427572128889195782774158241401584273/ 140500611775287989854314260624451156993638400000000*t^39+ 261827871161505491424783528224633956436102641/ 2195322058988874841473660322257049328025600000000*t^38+ 11870673828323921249575070733682372446383434379/ 2036240750366492606584264646731176188313600000000*t^35+ 1186841508565149023024180623405892758505959951/ 2195322058988874841473660322257049328025600000000*t^37, 1/ 60415263063373835637355132068513997507264512000000000*t+ 1451975901313761322372659726514496809324369/ 24507345504149309236754624214974909644800000000*t^18+231927781/ 3185954915539410200778101147946738253824000000000*t^2+1953910415490184721/ 1438458644366043705651312668297952321601536000000000*t^3+ 50720782828337470035571/43132193234364129104986886605635751772160000000*t^5+ 6652307692498452233/41555933680948828705801319321044412006400000000*t^4+ 713173077338632007395005991/477893237330911530116715172192010738073600000000*t^ 6+28597030175021528349441190434121/ 50769128624683895493575741234045376056524800000000*t^7+ 731323529725171655002301058333738048041172323389/ 13278079794148095744473655399673406045552640000000*t^25+ 64981318964001006795920070342959/ 743389480292529046848223601187572259225600000000*t^8+ 1317259514684763154048192368697892381/ 191794485915472494086841689106393642880204800000000*t^9+ 200495628657060973375581148473882959/ 637190983107882040155620229589347650764800000000*t^10+ 1113027039431771214500064640011314864459/ 119871553697170308804276055691496026800128000000000*t^11+ 75544358586652157732938409949384112297/ 398244364442426275097262643493342281728000000000*t^12+ 63625830845545254314545343644396165277429/ 22475916318219432900801760442155505025024000000000*t^13+ 124438257103014135418266970204747627927/ 3871820209856922119001164589518605516800000000*t^14+ 30993653126484603738296723471358752195381323/ 4498726010377058513739922277224415232000000000*t^22+ 21317890806797703484270064608945051490557/ 74421363714429460011523937014676025507840000000*t^15+ 8137697653770267775043315429737493050703981309/ 33195199485370239361184138499183515113881600000000*t^19+ 10598286076768023728608825721127398272209719/ 12253672752074654618377312107487454822400000000*t^20+ 165151706234551399320648435109020268869227701/ 5251574036603423407875990903208909209600000000*t^24+ 7113257498035156895600383928704328887969/ 3462994473412402392150109943420367667200000000*t^16+ 1158856685509868110924714995548911406335009511/ 95897242957736247043420844553196821440102400000000*t^17+ 30110036785079638684809512454503990515434551/ 355178920349989988938472814709781299200000000*t^26+ 9982032445849861994295408475548240508205120027/ 2230168440877587140544670803562716777676800000000*t^36+ 72578180100136523518341274486526550840012590729/ 27662666237808532800986782082652929261568000000000*t^21+ 8487100222248240082285556139304853579552651/ 66905053226327614216340124106881503330304000000000*t^42+ 2753096574824455845543087859307908428088476369971/ 23974310739434061760855211138299205360025600000000*t^27+ 871070481287446125098832029199537093961411222381/ 55325332475617065601973564165305858523136000000000*t^23+ 2644339340049342433258377897252472904750939221/ 19912218222121313754863132174667114086400000000*t^30+ 374706248386372925576429334636942482099297087/ 4916597091881805865398304240658546688000000000*t^32+ 45078364059729821187981962670767450801681940951499/ 958972429577362470434208445531968214401024000000000*t^33+ 76371038771528328099965748046020649516080758759/ 557542110219396785136167700890679194419200000000*t^28+ 6034586567555971291187175208869927146017506357353/ 41955043794009608081496619492023609380044800000000*t^29+ 38686321080455212361429783870450353088047176983287/ 359614661091510926412828167074488080400384000000000*t^31+ 15987333426095627372317549592847188704465988309/ 637190983107882040155620229589347650764800000000*t^34+ 4758801132491037164062107238850403195313965991/ 2876917288732087411302625336595904643203072000000000*t^41+ 13254220507386462191084083889303011505978809/ 955786474661823060233430344384021476147200000000*t^40+ 24160495735040921516457265753545448921356168571/ 287691728873208741130262533659590464320307200000000*t^39+ 516098083439704913515348955653804/109894723324712387033933067993555591796875*t^ 43+62342335858733505235385140783248507518128697/ 159297745776970510038905057397336912691200000000*t^38+ 1191966056401623283994558949138479303875463971/ 103672695089444591398292804922374942097408000000*t^35+ 1847351246290383226062516025428279417679354759/ 1265506138738454873007019942784122277068800000000*t^37, 1/ 2658271574788448768043625811014615890319638528000000000*t+ 2939075571292006792952721827843793165258373769/ 114039964598389050538122085414612436307148800000000*t^18+8796093022207/ 2658271574788448768043625811014615890319638528000000000*t^2+ 11723462702371894571/126584360704211846097315514810219804300935168000000000*t^3 +1691584176882493116392579677/ 12658436070421184609731551481021980430093516800000000*t^5+ 9957626778573454243069/684239787590334303228732512487674617842892800000000*t^4+ 16725340399200996339454199/82177892781741949075427017411543852646400000000*t^6+ 27234985187792659633874581519657/ 303802465690108430633557235544527530322244403200000*t^7+ 214650995566158278926587689510608593807251756767/ 5480633310905405372953478776600655403417600000000*t^25+ 604080295701990639170444262673113851/ 37975308211263553829194654443065941290280550400000000*t^8+ 11875990008559383574042602698328555307/ 8438957380280789739821034320681320286729011200000000*t^9+ 2392869103093609573541164456562381497/ 33355562767908259841189858975024981370470400000000*t^10+ 98472059508621153377124130652127543630457/ 42194786901403948699105171603406601433645056000000000*t^11+ 2201791617432607852809687048756449225008519/ 42194786901403948699105171603406601433645056000000000*t^12+ 462387585084917580429320036815854567520603/ 545622244414706233178084115561292259917824000000000*t^13+ 6597691239985448739501271084186435345840217/ 632921803521059230486577574051099021504675840000000*t^14+ 9631254193131416920152674371500466334231232676601/ 2434314628927150886486836823273457775017984000000000*t^22+ 25542189464317500069222843307947092719831121/ 254623714060196242149772587261936387961651200000000*t^15+ 334025333899106780870598391449032913097241176561/ 2921177554712581063784204187928149330021580800000000*t^19+ 54914814891417307493035972433871050403893304521/ 127007719770112220164530616866441275218329600000000*t^20+ 10126477350954216544315323514486078099866346571711/ 486862925785430177297367364654691555003596800000000*t^24+ 5719936892226154930321583827936301661310226131/ 7384087707745691022343405030596155250887884800000000*t^16+ 4133338696122521718437116937435182328208019159/ 843895738028078973982103432068132028672901120000000*t^17+ 189658880016253284996717868923201991471580701335989/ 2921177554712581063784204187928149330021580800000000*t^26+ 76419195907987230361336760543226573525971022279023/ 8438957380280789739821034320681320286729011200000000*t^36+ 311058430181195665991521363068591369137064994933/ 221301329902468262407894256661223434092544000000000*t^21+ 139772115315571112338159099146841162658974698967/ 126584360704211846097315514810219804300935168000000000*t^42+ 80107221709434938734987028064563754081220726492757/ 843895738028078973982103432068132028672901120000000*t^27+ 638087052436939880239289887952629495474913827489/ 65792287268301375310455049277661020946432000000000*t^23+ 1037548345061750027361142320804578190386402123863937/ 7384087707745691022343405030596155250887884800000000*t^30+ 976829474482150683722700834669155898695179317223/ 10059151359203102836722466211873792458752000000000*t^32+ 2798633992880916161156989240852003532521517413089827/ 42194786901403948699105171603406601433645056000000000*t^33+ 22479080098026858870143918735126313849744658330049/ 183455595223495429126544224362637397537587200000000*t^28+ 1031346540887061751236627854633249508388280399276159/ 7384087707745691022343405030596155250887884800000000*t^29+ 1971183110082119262491430629872731863911365211695517/ 15823045088026480762164439351277475537616896000000000*t^31+ 333418592632192200819123250564337605308549448639993/ 8438957380280789739821034320681320286729011200000000*t^34+ 41428478605551815020000004118671950086038467709/ 4364977955317649865424672924490338079342592000000000*t^41+ 32464895802101629503158150306673379580088833167/ 550366785670486287379632673087912192612761600000000*t^40+ 188138773271976387891905830911946970432908994333/ 666233477390588663670081656895893706847027200000000*t^39+ 1589515884519518540721608102017614123903854761/ 531654314957689753608725162202923178063927705600000000*t^44+ 19975799973393935798635658529657099453127409401/ 241661052253495342549420528274055990029058048000000000*t^43+ 41136508921716806494602931814620605821447398751929/ 37975308211263553829194654443065941290280550400000000*t^38+ 172148717632719533853679199725698440329344436506257/ 8438957380280789739821034320681320286729011200000000*t^35+ 25972809452356707526458295695663418708100423981531/ 7595061642252710765838930888613188258056110080000000*t^37, 1/ 119622220865480194561963161495657715064383733760000000000*t+ 189951857334574917461732280166447812483481/ 17427148781561264284843186507327576473600000000*t^18+50991843607/ 346731074972406361049168584045384681346039808000000000*t^2+ 738578159045522380223/119622220865480194561963161495657715064383733760000000000 *t^3+42290525742870331352955814513/ 2848148115844766537189599083229945596771041280000000000*t^5+ 21358421853984494776261/16511003570114588621388980192637365778382848000000000*t ^4+14936837205073470411899102569/ 550366785670486287379632673087912192612761600000000*t^6+ 74728962126785227595468019689267/ 5357018399708024834839371316419960370102272000000000*t^7+ 2511382676940705078784133425244862659044026707929489/ 93894992830047248478777991754833371322122240000000000*t^25+ 14018378085182415571998047657421227/ 4953301071034376586416694057791209733514854400000000*t^8+ 963415274897701132877405046969016034819/ 3417777739013719844627518899875934716125249536000000000*t^9+ 3519227514899134799853073757624660089/ 220146714268194514951853069235164877045104640000000*t^10+ 11959317936942241707520356765749031755027/ 20865553962232721884172887056629638071582720000000000*t^11+ 25722597844274368632706641738170305576733/ 1834555952234954291265442243626373975375872000000000*t^12+ 176364383169970342539699197639680918767072431/ 712037028961191634297399770807486399192760320000000000*t^13+ 6808761372091481051859743203240549403484857/ 2063875446264323577673622524079670722297856000000000*t^14+ 697304492263583437727716285467667964056190591041/ 317519299425280550411326542166103188045824000000000*t^22+ 170727187294608527814733486015769578975827257179/ 4984259202728341440081798395652404794349322240000000000*t^15+ 88504976154108477024992156135566859187822352353829/ 1708888869506859922313759449937967358062624768000000000*t^19+ 3194465450948584898756366196912890706892781239/ 15240926372413466419743674023972953026199552000000*t^20+ 1401617483385085925729518245184141954270164752761/ 105839766475093516803775514055367729348608000000000*t^24+ 274112770156749546220226807485149404756182019/ 963141874923351002914357177903846337072332800000000*t^16+ 197549245932289931724941273690277900910655707321/ 102241214414940337232447146577485226550755328000000000*t^17+ 625758834899935408085994467991078014067398801569/ 13138729631390919327365236227562890539827200000000*t^26+ 2015273976374293666683747209315694409130368993543/ 122303730148996952751029482908424931691724800000000*t^36+ 8261988873811009155771940430587791768234051035641/ 11332154307074667919852516246272993090600960000000000*t^21+ 35690395034539602412259111657727849852998820709/ 5503667856704862873796326730879121926127616000000000*t^42+ 128054858034455046846930767306207719081132232836199929/ 1708888869506859922313759449937967358062624768000000000*t^27+ 631341118652809758499643273875541783375060343587789/ 109544158301721789891907657047305599875809280000000000*t^23+ 9011463666834562257863493940906655331249450063719/ 64209458328223400194290478526923089138155520000000*t^30+ 238456005476719295359975563675120895339779637795219/ 2063875446264323577673622524079670722297856000000000*t^32+ 494360230850634487457617187654589121395592342440160873/ 5696296231689533074379198166459891193542082560000000000*t^33+ 57415056730721904973803719993470772950442787918119/ 550366785670486287379632673087912192612761600000000*t^28+ 3884525626268074453837710994368811663208539740235517/ 30207631531686917818677566034256998753632256000000000*t^29+ 675362931763335024399298848786932626711925372633022019/ 4984259202728341440081798395652404794349322240000000000*t^31+ 314926038167788774250069017538249355802015854474277/ 5503667856704862873796326730879121926127616000000000*t^34+ 235457221704862051341816918833985447866983714350761/ 5696296231689533074379198166459891193542082560000000000*t^41+ 6098306122452419683902902382013581728316300127/ 30020006491117433857070873077522483233423360000000*t^40+ 195488190391115011457764277611442625125331276072539/ 244126981358122846044822778562566765437517824000000000*t^39+ 18666249674074243741873704899818643573681701181/ 346731074972406361049168584045384681346039808000000000*t^44+ 103537504005512749467288942622106408/ 54397888045732631581796868656810017939453125*t^45+ 6774892170875807002568924237801932712643077986871/ 9201709297344630350920243191973670389567979520000000000*t^43+ 12881572286882789789667612368087121048378145143553/ 4953301071034376586416694057791209733514854400000000*t^38+ 62578092470576957201216357854190953941914715175783037/ 1898765410563177691459732722153297064514027520000000000*t^35+ 932375322645928811545754439617926067323376747334019/ 131452989962066147870289188456766719850971136000000000*t^37, 1/ 5502622159812088949850305428800254892961651752960000000000*t+ 274308995052618045011722014966572788161975089745223/ 61140246220134321665003393653336165477351686144000000000*t^18+162139963543/ 25357705805585663363365462805531128538993786880000000000*t^2+ 276966811841094148141/687827769976511118731288178600031861620206469120000000000 *t^3+422912627561399256271300715621/ 262029626657718521421443115657154994902935797760000000000*t^5+ 7033725044402115547026521/ 62529797270591919884662561690911987420018769920000000000*t^4+ 927662473601425377852318791375771/ 262029626657718521421443115657154994902935797760000000000*t^6+ 192300882812501146583250582510551/ 90563235019948797726305224305468316210692096000000000*t^7+ 18328725541331574500274364478345379388766564219168777/ 1042558196250869448626431494657115364335288320000000000*t^25+ 25673985961752824826403908988454281/ 52127909812543472431321574732855768216764416000000000*t^8+ 8682221883915887843699148602274111549149/ 157217775994631112852865869394292996941761478656000000000*t^9+ 547681376155503085879620216074891458777379/ 157217775994631112852865869394292996941761478656000000000*t^10+ 751306181105821438095884874038354612546627/ 5458950555369135862946731576190729060477829120000000000*t^11+ 20101583122763098670706655593396331338568637/ 5458950555369135862946731576190729060477829120000000000*t^12+ 18532770905074085426262474628623516863341449403/ 262029626657718521421443115657154994902935797760000000000*t^13+ 267137120471088163766128357000033896661711400773/ 262029626657718521421443115657154994902935797760000000000*t^14+ 1150918548279631334810538048103819396743395943528103/ 975296377137910129360210107905043405345914880000000000*t^22+ 59474013408245504561139458871446775846605659361/ 5210816439215993323721880140909332285001564160000000000*t^15+ 127103948110642592174303607177561222991446986071/ 5573517299866389423314870582611067673772032000000000*t^19+ 969376122920728757762627170361163519382739484664283/ 9826110999664444553304116837143312308860092416000000000*t^20+ 10245157638012469682997576302358893305557956285520341/ 1259757820469800583756938056044014398571806720000000000*t^24+ 834892945585539274421321172167740857995699723203/ 8188425833053703794420097364286093590716743680000000000*t^16+ 1565602636972733017500475943724382031435346915917/ 2108284352418424885000117022528833292322471936000000000*t^17+ 92205258487221361786629464829433964533312965783980753/ 2748562517388655819106046667732395051429396480000000000*t^26+ 298007379089116010723721824337298608059191014174605869/ 10917901110738271725893463152381458120955658240000000000*t^36+ 226243133314592890729367618327621783210654058691447/ 617024238597453347146255374388905011545374720000000000*t^21+ 176108991337977482402455793896065133621324004559587/ 6093712247853919102824258503654767323324088320000000000*t^42+ 139156680756186344137414379829249937030696007293431377/ 2456527749916111138326029209285828077215023104000000000*t^27+ 101696324375425090583213095106830330216195570101811/ 30725800499263428872120440391317424355409920000000000*t^23+ 1160601888391406606685942963626743322481073270874033899/ 8734320888590617380714770521905166496764526592000000000*t^30+ 255137638293747375851879449715891915688738570684399629/ 1976516580392273329687609708620781211552317440000000000*t^32+ 27777833101400043731683681328935696704583145062968465553/ 262029626657718521421443115657154994902935797760000000000*t^33+ 7183315759232068258050451209746215943412151392061243/ 84707853445383142700897558940890623352242176000000000*t^28+ 6888080073605127059465685954885743666065267301069576083/ 61140246220134321665003393653336165477351686144000000000*t^29+ 7972728351287283802444118763672863334757117506760408711/ 57318980831375926560940681550002655135017205760000000000*t^31+ 1548646845421523828525651472268378590527678743498648171/ 20156125127516809340111008896704230377148907520000000000*t^34+ 38179246331390792344122068264782448013040728887634911/ 262029626657718521421443115657154994902935797760000000000*t^41+ 11590814361595975510624635634013855325771688710841817/ 19652221999328889106608233674286624617720184832000000000*t^40+ 38694862992731389251540858943794551105634270041266727/ 19652221999328889106608233674286624617720184832000000000*t^39+ 3442116679210442967253053866215327275991883089397/ 7018650714046031823788654883673794506328637440000000000*t^44+ 192832714641173486049725471168421909968816649540811/ 5502622159812088949850305428800254892961651752960000000000*t^45+ 89490081186778333402864344201716076557900889589099/ 20230228528720915256802593488236231224123719680000000000*t^43+ 350923027193976051443948867093355168290742509719/ 289611692621688892097384496252644994366402723840000000000*t^46+ 51108709716404566236739426682367889059248012636216257/ 9248104470272418403109757023193705702456557568000000000*t^38+ 76362922197153903196237769402959569666090711976904917/ 1559700158676895960841923307483065445850808320000000000*t^35+ 189060686920504398440398866186564466180245587783998189/ 14292525090421010259351442672208454267432861696000000000*t^37, 1/ 258623241511168180642964355153611979969197632389120000000000*t+ 2820502131604197962459839615020126141853908027263/ 1567698621029085170897522914188106807111581696000000000*t^18+499069107643/ 1834207386604029649950101809600084964320550584320000000000*t^2+ 332360175968531582213/ 12931162075558409032148217757680598998459881619456000000000*t^3+ 44406444957027412567991100434557/ 258623241511168180642964355153611979969197632389120000000000*t^5+ 1463282169837817679478977/ 152850615550335804162508484133340413693379215360000000000*t^4+ 181926299708197584114211207819/ 402503266755327989894689885802081405380853760000000000*t^6+ 182609947583774437416255995620362929/ 577284021230286117506616864182169598145530429440000000000*t^7+ 15835801014042311435170134846514797529642088152294038061/ 1421006821489935058477826127217648241588997980160000000000*t^25+ 68685508811735493529588928937905167/ 818842583305370379442009736428609359071674368000000000*t^8+ 78230449235873574690597044007414476398741/ 7389235471747662304084695861531770856262789496832000000000*t^9+ 10732036181792684328504705271051310341/ 14448835216857927842373483080074717116235776000000000*t^10+ 298603402200230778156091021012137218551084157/ 9236544339684577880105869826914713570328486871040000000000*t^11+ 625802463892862370970793751811445377173029/ 661690976408380104599603827417058067936706560000000000*t^12+ 57551027725470655454108207391472941745980703/ 2914885787671661658416053594292611777618457395200000000*t^13+ 26888356320259984079309593862110077183488262949/ 87343208885906173807147705219051664967645265920000000000*t^14+ 2065159464049129840679754457999412239969118773369243/ 3359354187919468223351834816117371729524817920000000000*t^22+ 5002715948717281171291341765303925720970964625681/ 1346996049537334274182106016425062395672904335360000000000*t^15+ 63186120126750773394518721152962006683970974847454713/ 6465581037779204516074108878840299499229940809728000000000*t^19+ 6034447010181044796007697129548961751488666560597/ 133688585029448225215021997784262752501497856000000000*t^20+ 17469540695057487924002216168208112262857681797263/ 3619993736982185585508442689781650570608640000000000*t^24+ 113513219735660421314093857954341758754684410233/ 3184387823965329253385593419444591951945400320000000000*t^16+ 11995629727474328536036784462727205707578989289678567/ 43103873585194696773827392525601996661532938731520000000000*t^17+ 229160852488755383858912792898834223575435031881766879/ 10078062563758404670055504448352115188574453760000000000*t^26+ 23272242878154019196775105981132490688732595010218613/ 559892364653244703891972469352895288254136320000000000*t^36+ 80664400990895760263434555058375775751435698581967919/ 450563138521198920980774137898278710747731066880000000000*t^21+ 9102817259684086286616511225263719896417638141770173/ 87343208885906173807147705219051664967645265920000000000*t^42+ 4630797152079998394681648036024568009278874356938365877/ 112640784630299730245193534474569677686932766720000000000*t^27+ 380002547466811847761901276220979523137359685856189/ 206541689170048700360149146397913988603052032000000000*t^23+ 2439291595601018051575213974705848786694335638796286269/ 20380082073378107221667797884445388492450562048000000000*t^30+ 1300012960882916816648979642839817765795812782069705049/ 9553163471895987760156780258333775855836200960000000000*t^32+ 2096071356675580778396530304576693607440254461277166915947/ 17241549434077878709530957010240798664613175492608000000000*t^33+ 48038947660909110238678533938340799994416593859957867/ 727860074049218115059564210158763874730377216000000000*t^28+ 2436702551681238758573675365836978247614082308062833132703/ 25862324151116818064296435515361197996919763238912000000000*t^29+ 182274902477778058121189276139325403054583941830239943163/ 1346996049537334274182106016425062395672904335360000000000*t^31+ 3309735083498451510186445529470595508095924530967639/ 34373557216019745693485913112574445087621120000000000*t^34+ 940837735952335027599025901540662919306616047833379429/ 2173304550514018324730792900450520840077291028480000000000*t^41+ 101298426613138891526381245412574089697378696312841/ 68236881942114198286834144702384113255972864000000000*t^40+ 9512816503618090900872838340304698380651869488673893/ 2220323158578023528871603323777575377482809344000000000*t^39+ 106271990693074442678888540008785865408699856176937/ 35273218973154416345194265569232403160010588160000000000*t^44+ 45361105584983995647044252937847808918/ 58804116977436974739922415018011629392548828125*t^47+ 4442105151777100748173755015910111631801132182954133/ 13611749553219377928577071323874314735220928020480000000000*t^45+ 260726163676739312511106805973678754422527045697653731/ 12931162075558409032148217757680598998459881619456000000000*t^43+ 516310164907686115722311990476634926522213316453/ 22644535637086785801853108760494876102722846720000000000*t^46+ 553255011501812792312014183471155071069364706877751257/ 52405925331543704284288623131430998980587159552000000000*t^38+ 207564555456986332269149852626263812118263523977068347593/ 3078848113228192626701956608971571190109495623680000000000*t^35+ 829551069838664042918261857611952042685868431360968797031/ 36946177358738311520423479307658854281313947484160000000000*t^37, 1/ 12413915592536072670862289047373375038521486354677760000000000*t+ 58220892115692273019486832665268276240564090581869617/ 82759437283573817805748593649155833590143242364518400000000*t^18+ 140737488355327/12413915592536072670862289047373375038521486354677760000000000* t^2+19941610628480639110489/ 12413915592536072670862289047373375038521486354677760000000000*t^3+ 222034700951045383279776972949309/ 12413915592536072670862289047373375038521486354677760000000000*t^5+ 9903500372669036018833355431/ 12413915592536072670862289047373375038521486354677760000000000*t^4+ 701410785921151482756788640303265987/ 12413915592536072670862289047373375038521486354677760000000000*t^6+ 81826222153177205291401521545502113363/ 1773416513219438952980327006767625005503069479239680000000000*t^7+ 465498609371378404721338035836011162321604240979349686277/ 68208327431516882806935654106447115596271903047680000000000*t^25+ 708716941840155498005064020089885879279/ 50669043234841112942295057336217857300087699406848000000000*t^8+ 21180686925545970422868265353304990823/ 10659152596360264180197307328430503413992904458240000000*t^9+ 54986787628685553661217484538827479744849923/ 354683302643887790596065401353525001100613895847936000000000*t^10+ 1013911974946390805740858183995246850315481741/ 136416654863033765613871308212894231192543806095360000000000*t^11+ 24795144188527844867769938251315125721394227303/ 104318618424672879587078059221625000323709969367040000000000*t^12+ 3186970596585060782616894073293948533248386567559/ 591138837739812984326775668922541668501023159746560000000000*t^13+ 53728428709362410497548539925839277286680654082681/ 591138837739812984326775668922541668501023159746560000000000*t^14+ 104456149698661011301322254069721614328515316836223527/ 336001612963137353728747064563778894562915778560000000000*t^22+ 4887016554386033546317845171776166954884899852086073/ 4137971864178690890287429682457791679507162118225920000000000*t^15+ 1012919257678534229083756980416637846336075963269306021/ 248278311850721453417245780947467500770429727093555200000000*t^19+ 2268147838824093800994363148769326055419849824549909941/ 112853778113964297007838991339757954895649875951616000000000*t^20+ 189035000041306730529138906818538168839527405909265694773/ 68208327431516882806935654106447115596271903047680000000000*t^24+ 50364432629365876749452221415301280058194599609799879/ 4137971864178690890287429682457791679507162118225920000000000*t^16+ 210585001440891980494935244295303687013620493878146871/ 2068985932089345445143714841228895839753581059112960000000000*t^17+ 92151222185239492643059411846809274656340422675420179921/ 6200757039228807527903241282404283236024718458880000000000*t^26+ 291379167388346953999433414386353045171089067476932437793/ 4967553258317756170813240915315476205890950922240000000000*t^36+ 75340778751168128228290993186859578991524518412214048759/ 886708256609719476490163503383812502751534739619840000000000*t^21+ 51058548516787511130319946306137756653473571032838978281/ 161219683019948995725484273342511364136642679930880000000000*t^42+ 25454260722038222580722383219003582370378265277518457678009/ 886708256609719476490163503383812502751534739619840000000000*t^27+ 67619646445694541054052427599728364730380461214739906507/ 68208327431516882806935654106447115596271903047680000000000*t^23+ 127732011473673023577055859072934022359698037924645887068169/ 1241391559253607267086228904737337503852148635467776000000000*t^30+ 281238597543242563559311152048791454320004156757840444777927/ 2068985932089345445143714841228895839753581059112960000000000*t^32+ 544232590947723001304965711318810173474077357090295261292343/ 4137971864178690890287429682457791679507162118225920000000000*t^33+ 8743787438707846165524686754581072205928632867480640616411/ 177341651321943895298032700676762500550306947923968000000000*t^28+ 18727182881878939826403087811718862156963717678783031977163/ 248278311850721453417245780947467500770429727093555200000000*t^29+ 23552443775735377235682872769256637949311659946685196641611/ 188089630189940495013064985566263258159416459919360000000000*t^31+ 468674127685805796107612689974813754973518083464689877129417/ 4137971864178690890287429682457791679507162118225920000000000*t^34+ 1976957055025835696388757154751387618314836973003758729213/ 1773416513219438952980327006767625005503069479239680000000000*t^41+ 1172515423763184811386401696950902742403823005339170109357/ 354683302643887790596065401353525001100613895847936000000000*t^40+ 594359715187336588751721029791632872522816946213132488311/ 70936660528777558119213080270705000220122779169587200000000*t^39+ 174265927176669043527330458924657338622683516804389777581/ 12413915592536072670862289047373375038521486354677760000000000*t^44+ 184073286012359886016543759331048085438071236829063951/ 12413915592536072670862289047373375038521486354677760000000000*t^47+ 5516994249383296071214195242422482492286460673697/ 11234312753426310109377637146944230804091842854912000000000*t^48+ 2312298112242744359481701435621987614924429001611580389/ 1128537781139642970078389913397579548956498759516160000000000*t^45+ 922936086137887597959927836575402927801589743990796603219/ 12413915592536072670862289047373375038521486354677760000000000*t^43+ 86877758427920057395248446210461090365846860295344567/ 400448890081808795834267388624947581887789882408960000000000*t^46+ 6508765844290426071056374892223375029576620040698657830629/ 354683302643887790596065401353525001100613895847936000000000*t^38+ 51205566526184803007167392304188124200195553906440854491049/ 591138837739812984326775668922541668501023159746560000000000*t^35+ 367564345164971768363846527830220031479142993447266233007/ 10493588835617981970297792939453402399426446622720000000000*t^37, 1/ 608281864034267560872252163321295376887552831379210240000000000*t+ 1172842808686587755537825705076959641292763406570837/ 4368899404110231709202414446033867534924605554688000000000*t^18+4043309297/ 8737798808220463418404828892067735069849211109376000000000*t^2+ 6647203558464378409649/ 67586873781585284541361351480143930765283647931023360000000000*t^3+ 9025881375227415152166930553447/ 4945381008408679356684976937571507129167096190074880000000000*t^5+ 1707133007209362054748411/ 26213396424661390255214486676203205209547633328128000000000*t^4+ 48352397254655830505663578529/ 6988375479781762264786586690536711599452848128000000000*t^6+ 4010195843783745105712912374779743641061/ 608281864034267560872252163321295376887552831379210240000000000*t^7+ 27417832445168621450129615603878349584460740285520295161/ 6779326661550359548762367243845656519710594826240000000000*t^25+ 950593877766344626663821819428549249/ 416085657534307781828801375812749289040438624256000000000*t^8+ 3527065672579120352198832615165730247334697/ 9655267683083612077337335925734847252183378275860480000000000*t^9+ 4747003700823256387135851797832632528027/ 149790836712350801458368495292589744054557904732160000000*t^10+ 48467259963393194570527963777430841863305675651/ 28965803049250836232012007777204541756550134827581440000000000*t^11+ 218997533280906429161617318755552729039864043/ 3744770917808770036459212382314743601363947618304000000000*t^12+ 125201062377971557862372310595238063668952524716373/ 86897409147752508696036023331613625269650404482744320000000000*t^13+ 31847094963244051701321932035537903960912129/ 1213077718758914815827409258929298218776788992000000000*t^14+ 286095601917575081171265602788020202364526341076951509/ 1872385458904385018229606191157371800681973809152000000000*t^22+ 8273572790324375551110184044897821850547346340991313/ 22528957927195094847120450493381310255094549310341120000000000*t^15+ 958420530246648260867694603163683758292115930317321139/ 577117518059077382231738295371247985661814830530560000000000*t^19+ 49141849043634443839545866536059446377407648238507/ 5637289553690621560261179930366280690225297489920000000*t^20+ 2888391711790806101625638682478787169358510949747173309/ 1872385458904385018229606191157371800681973809152000000000*t^24+ 35488269970138690594094804581725765821169079477723/ 8737798808220463418404828892067735069849211109376000000000*t^16+ 12707906280775267156614051686723142474981628355323/ 349542681355750910442605692446880557128661073920000000000*t^17+ 1951919586389121371511712641622681591421327287144220219/ 208042828767153890914400687906374644520219312128000000000*t^26+ 7414921104501598614874001025922636300731530956216329943/ 96019767123301795806646471341403682086255067136000000000*t^36+ 78001366466496583821563331911455774401280072028208869801/ 1987849228870155427687098572945409728390695527383040000000000*t^21+ 63732414654045405673524478121994416129856379478158007/ 76423896281811633397126783312545787782937706496000000000*t^42+ 5487756525944858824438836049164429985687774581623588230273/ 283978461267165061098156938992201389770099361054720000000000*t^27+ 441728139070568110708818080330470294654271996499368730389/ 851935383801495183294470816976604169310298083164160000000000*t^23+ 44384628313197016760328785665916527389836400151524080657/ 524267928493227805104289733524064104190952666562560000000*t^30+ 188043846112576914325966303788137186331588054491184486217/ 1456299801370077236400804815344622511641535184896000000000*t^32+ 178632461885531864420443694073988896904098836062455477491279/ 1325232819246770285124732381963606485593797018255360000000000*t^33+ 4244330264086481689328715573949902863813157869056693/ 119786671288106008459446368828441673640968192000000000*t^28+ 345588558350129896260819471237255670577486322809308507523233/ 5963547686610466283061295718836229185172086582149120000000000*t^29+ 1978130700519891320442960285228274721088357894853183869933019/ 17890643059831398849183887156508687555516259746447360000000000*t^31+ 1097986987875409700618907736954268665408912958894231644757/ 8737798808220463418404828892067735069849211109376000000000*t^34+ 73549584097341973359336094510892934126853747939268256104831/ 28965803049250836232012007777204541756550134827581440000000000*t^41+ 198164316258338006653356469638833052108654523338278897/ 29958167342470160291673699058517948810911580946432000000*t^40+ 13084309789805937593703036098427656205831883308886519115317/ 877751607553055643394303265975895204743943479623680000000000*t^39+ 51361354587799812708580180233659682860661342623405803/ 970866534246718157600536543563081674427690123264000000000*t^44+ 29216830885952333428308801965536381703145585951101280037/ 202760621344755853624084054440431792295850943793070080000000000*t^47+ 87176517890549500795745183943750553204/ 278845328893007589895761129278958371635634765625*t^49+ 252666049710437124246461897868464734303445397512317/ 26213396424661390255214486676203205209547633328128000000000*t^48+ 659281547047362832363379094714950649668691621395122787359/ 67586873781585284541361351480143930765283647931023360000000000*t^45+ 140488655635987754192360633421132926455212835638317032997491/ 608281864034267560872252163321295376887552831379210240000000000*t^43+ 36474323269104088080553537418081824310131454816481269/ 26213396424661390255214486676203205209547633328128000000000*t^46+ 36587692851249580310102334841009622731487158495963679447/ 1248256972602923345486404127438247867121315872768000000000*t^38+ 7059175596938975850812703158094352262044196382299042715607189/ 67586873781585284541361351480143930765283647931023360000000000*t^35+ 4399419307595415647436458367321630737772786537466213960435783/ 86897409147752508696036023331613625269650404482744320000000000*t^37, 1/ 30414093201713378043612608166064768844377641568960512000000000000*t+ 537528683943931831449168831455135143802452955916263911/ 5372565483432852507262428575528134401055933857792000000000000*t^18+ 562949953421311/ 30414093201713378043612608166064768844377641568960512000000000000*t^2+ 29912416060002198961871/ 5069015533618896340602101361010794807396273594826752000000000000*t^3+ 19684243484178137890189345350503/ 107851394332316943417065986404484995902048374358016000000000000*t^5+ 1148232938797233633514917367/ 220391979722560710460960928739599774234620591079424000000000000*t^4+ 382624773968657203961051315710282661/ 460819593965354212782009214637344982490570326802432000000000000*t^6+ 452832716802293553985492770118287434303/ 490549890350215774896977551065560787812542605950976000000000000*t^7+ 1426973655431793366016723241028431116840144520982617473219/ 614376478703001334106589531473512622098772656128000000000000*t^25+ 5560845379955533302451117031110560567188063/ 15207046600856689021806304083032384422188820784480256000000000000*t^8+ 23826242286031130270191277097545948655322019/ 362072538115635452900150097215056771956876685344768000000000000*t^9+ 2298203200697384566752285798024672649627237469/ 362072538115635452900150097215056771956876685344768000000000000*t^10+ 7222717017820017998945052119013040992447970499/ 19571488546791105562170275525138203889560901910528000000000000*t^11+ 237344400242818013259331173360263542774798462451/ 16840583168169090832565120800700314974738450481152000000000000*t^12+ 819110458154620255040287777745222272373084856727439/ 2172435228693812717400900583290340631741260112068608000000000000*t^13+ 47290091937036561736949565065437072116828748929709/ 6370777796756049024636072091760529711851202674688000000000000*t^14+ 10894157499458992260080649353192062810572352147937813015579/ 149088692165261657076532392970905729629302164553728000000000000*t^22+ 188795386367855856612837228465712355330906893645444487/ 1689671844539632113534033787003598269132091198275584000000000000*t^15+ 295063102083126429237031390882554751303326829499172817937/ 447266076495784971229597178912717188887906493661184000000000000*t^19+ 1834123237045857739332159527552337068053933480147723813/ 497515101775066708820464047733834470398116233216000000000000*t^20+ 17748825050691790306324895211624827185048708833934425276367/ 21298384595037379582361770424415104232757452079104000000000000*t^24+ 2239475934866971908886722412335768199867682954541861497/ 1689671844539632113534033787003598269132091198275584000000000000*t^16+ 2522963942216503464395554046395485537229413704736856477/ 198784922887015542768709857294540972839069552738304000000000000*t^17+ 3520988101371341758553989999918035947421091509781244329789/ 614376478703001334106589531473512622098772656128000000000000*t^26+ 14652099025649966592828298953406788584434918838826827262180277/ 153606531321784737594003071545781660830190108934144000000000000*t^36+ 2631591192661390264452117410482144819593850296027919471589/ 149088692165261657076532392970905729629302164553728000000000000*t^21+ 703120855380548457780687346729109317970516871566486081976779/ 362072538115635452900150097215056771956876685344768000000000000*t^42+ 267596785981212789544738551955530406152156834526747719327247/ 21298384595037379582361770424415104232757452079104000000000000*t^27+ 5612968547546761107047250839903417630969738134507000086577/ 21298384595037379582361770424415104232757452079104000000000000*t^23+ 212012797624246250355534560290001231986419362861958472471627/ 3172099833303439512266646658955441055942599245824000000000000*t^30+ 1276583059352015166124939146787309681196919875910634815701737/ 10908928695019145639746272656407736314339182772224000000000000*t^32+ 26110496228634554835505654863579614196544451738356761398741107/ 198784922887015542768709857294540972839069552738304000000000000*t^33+ 522758122810998002395738447128300412523898825190578894971377/ 21298384595037379582361770424415104232757452079104000000000000*t^28+ 6384357282295158638220044676666924958742767186102098320770619/ 149088692165261657076532392970905729629302164553728000000000000*t^29+ 41785995329156136453163578463580629042783220632088984835868847/ 447266076495784971229597178912717188887906493661184000000000000*t^31+ 905133482281857870886342411438343095249557985904034551110193/ 6854652513345363543748615768777274925485156990976000000000000*t^34+ 1881238825477315218459436456002445888595770124923812054222389/ 362072538115635452900150097215056771956876685344768000000000000*t^41+ 8717489346550264984982966567185946097710210240594200234762223/ 724145076231270905800300194430113543913753370689536000000000000*t^40+ 17643100671237575391968782803974272391740469740650321562361873/ 724145076231270905800300194430113543913753370689536000000000000*t^39+ 2553629700142243229718616233152344741817031814023063912914273/ 15207046600856689021806304083032384422188820784480256000000000000*t^44+ 6053285248188621896314383785111649088103498225146815121/ 30414093201713378043612608166064768844377641568960512000000000000*t^50+ 4783349253347316890397072220444240323555504399724358112881/ 5069015533618896340602101361010794807396273594826752000000000000*t^47+ 190482409667483187018578496941878335879395307173682453871/ 30414093201713378043612608166064768844377641568960512000000000000*t^49+ 44064243331000172401262727781487668244178600097996462861/ 460819593965354212782009214637344982490570326802432000000000000*t^48+ 190353441111129271319450436184841423541366929132580413607991/ 5069015533618896340602101361010794807396273594826752000000000000*t^45+ 9452304342321544181349169798261145378552094707635447320994463/ 15207046600856689021806304083032384422188820784480256000000000000*t^43+ 34294292048983636751324515068682914095984019151681746661321/ 5069015533618896340602101361010794807396273594826752000000000000*t^46+ 94211734585047516506227263572798193764856942579153310188547439/ 2172435228693812717400900583290340631741260112068608000000000000*t^38+ 200664082071337225089906958121509268586114732672238523225665337/ 1689671844539632113534033787003598269132091198275584000000000000*t^35+ 148277070825450184254122756497526767126684584312676841466909329/ 2172435228693812717400900583290340631741260112068608000000000000*t^37]: end: ###start DesMajInv GnD:=proc(n,t,q1,q2) local gu: option remember: gu:= Fn(n,t,q1,q2)/t^((n-1)/2)/q1^(n*(n-1)/4)/q2^(n*(n-1)/4): gu/n!: end: #The sum of q1^inv(pi)*q2^maj(pi)*t^des(pi) over all n-perms Fn:=proc(n,t,q1,q2) local i: option remember: if n=0 then RETURN(1): fi: add(Fni(n,i,t,q1,q2),i=1..n): end: #The sum of q1^inv(pi)*q2^maj(pi)*t^des(pi) over all n-perms that end in i Fni:=proc(n,i,t,q1,q2) option remember: if n<1 then RETURN(FAIL): fi: if n=1 then if i=1 then RETURN(1): else RETURN(0): fi: fi: if i=n then RETURN(Fn(n-1,t,q1,q2)): fi: expand(q1*Fni(n,i+1,t,q1,q2)+ q1^(n-i)*(q2^(n-1)*t-1)*Fni(n-1,i,t,q1,q2)): end: ####MajiInvDes #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,scaling=constrained): end: #AveAndMomsS(F,x,N,n,Deg): Given an explicit (symbolic) #expression F(x,n), denoting a sequence of discrete #prob. distribution, parametrized by the discrete variable #n, tries to find explicit polynomial expressions, in n #of degree <=Deg #for the first N moments #For example, try: #AveAndMomsS(((1+x)/2)^n,x,4,n,10); AveAndMomsS:=proc(F,x,N,n,Deg) local mu,gu,n0,i,j,lu,mu1: mu:=[seq(AveAndMoms(simplify(subs(n=n0,F)),x,N),n0=1..Deg+5)]: gu:=[]: for j from 1 to N do mu1:=[seq(mu[i][j],i=1..nops(mu))]: lu:=GuessPOL1V(mu1,n): if lu=FAIL then print(`couldn't guess the`, j , `-th moment . `): RETURN(lu): else gu:=[op(gu),lu]: fi: od: gu: end: #AlphaA(F,x,N,n,Deg,ORD1): Given an explicit (symbolic) #expression F(x,n), denoting a sequence of discrete #prob. distribution, parametrized by the discrete variable #n, finds the asymptotic alpha coefficients starting with the #3rd (the asymptic skewness): #For example, try: #AlphaA(((1+x)/2)^n,x,4,n,10,2); AlphaA:=proc(F,x,N,n,Deg,ORD1) local mu,d2,di,i,gu,z,lu1,lu2,lu,i1: if nargs<>6 then ERROR(`We need six arguments`): fi: mu:=AveAndMomsS(F,x,N,n,Deg): if nops(mu)<3 then RETURN(FAIL): fi: d2:=degree(mu[2],n): gu:=[]: for i from 3 to N do if mu[i]=0 then gu:=[op(gu),0]: else di:=degree(mu[i],n): #gu:=[op(gu),coeff(mu[i],n,di)*n^(di-d2*i/2)/coeff(mu[2],n,d2)^(i/2)]: lu1:=normal(mu[i]/n^di): lu2:=normal(mu[2]/n^d2): lu:=lu1/lu2^(i/2): lu:=normal(subs(n=1/z,lu)): lu:=taylor(lu,z=0,ORD1+1): lu:=add(coeff(lu,z,i1)*z^i1,i1=0..ORD1): lu:=subs(z=1/n,lu)*n^(di-d2*i/2): gu:=[op(gu),lu]: fi: od: gu: end: #AlphaD(F,x,N,n,ORD1): Given an explicit (symbolic) #expression F(x,n), denoting a sequence of discrete #prob. distribution, parametrized by the discrete variable #n, finds the asymptotic alpha coefficients starting with the #3rd (the asymptotic skewness): #For example, try: #AlphaD(((1+x)/2)^n,x,4,n,2); AlphaD:=proc(F,x,N,n,ORD1) local mu,d2,di,i,gu,z,lu1,lu2,lu,i1: mu:=AveAndMomsD(F,x,N): if nops(mu)<3 then RETURN(FAIL): fi: d2:=degree(mu[2],n): gu:=[]: for i from 3 to N do if mu[i]=0 then gu:=[op(gu),0]: else di:=degree(mu[i],n): lu1:=normal(mu[i]/n^di): lu2:=normal(mu[2]/n^d2): lu:=lu1/lu2^(i/2): lu:=normal(subs(n=1/z,lu)): lu:=taylor(lu,z=0,ORD1+1): lu:=add(coeff(lu,z,i1)*z^i1,i1=0..ORD1): lu:=subs(z=1/n,lu)*n^(di-d2*i/2): gu:=[op(gu),lu]: fi: od: gu: end: #AveAndMomsSgf(F,x,N,n,Deg,z): Given an explicit (symbolic) #generating function F(x,z), such that the Maclaurin coeff. of z^n #is a prob. distribution, parametrized by the discrete variable #n, tries to find explicit polynomial expressions, in n #of degree <=Deg (by pure guessing) #for the first N moments #For example, try: #AveAndMomsSgf(1/(1-z*(1+x)/2),z,x,4,n,10); AveAndMomsSgf:=proc(F,z,x,n,N,Deg) local mu,gu,n0,i,j,lu,mu1,F1: F1:=taylor(F,z=0,Deg+10): mu:=[seq(AveAndMoms( expand(coeff(F1,z,n0)),x,N),n0=1..Deg+5)]: gu:=[]: for j from 1 to N do mu1:=[seq(mu[i][j],i=1..nops(mu))]: lu:=GuessPOL1Ve(mu1,n): if lu=FAIL then print(`couldn't guess the`, j , `-th moment . `): RETURN(lu): else gu:=[op(gu),lu[1]]: fi: od: gu: end: #GuessPOL1VeOld(L,x): Given a list L, tries to find #a polynomial P(x) of degree d such that L[i]=P(i) for i=s..nops(L) #for some starting place s. Also returns the starting place s #For example, try: #GuessPOL1VeOld([0,0,0,seq(i^2,i=4..10)],x); GuessPOL1VeOld:=proc(L,x) local s,mu: for s from 1 to nops(L)/2 do mu:=GuessPOL1V([op(s..nops(L),L)],x): if mu<>FAIL then RETURN(expand(subs(x=x-s+1,mu)),s): fi: od: FAIL: end: #RatToPol(f,s,n): given a rational function of the form #P(s)/(1-s)^k, finds the polynomial expression in n, #for the coeff. of s^n (for n sufficiently large) #for example, try: #RatToPol(s^5/(1-s)^3,s,n); RatToPol:=proc(f,s,n) local f1,gu,k,gu1,t,i: f1:=normal(f): gu:=denom(f1): k:=degree(gu,s): gu1:=normal(gu/(1-s)^k): if not type(gu1,numeric) then print(`Bad input`): RETURN(FAIL): fi: f1:=normal((1-s)^k*f1): f1:=expand(subs(s=1-t,f1)): expand(add(expand(binomial(k-i-1+n,k-i-1))*coeff(f1,t,i),i=0..k-1)): end: #AveAndMomsSrgfSlow(F,s,t,n,N): given a rational function F(s,t) #such that the coeff. of s^n is the prob. gen. function #of a discrete r.v. parametrized by n #finds polynomial expressions in n for the first N moments. #For example, try: #AveAndMomsSrgfSlow(1/(1-s*(1+t)/2),s,t,n,6); AveAndMomsSrgfSlow:=proc(F,s,t,n,N) local gu,F1,G1,i,mu: if normal(subs(t=1,F)-1/(1-s))<>0 then print(normal(subs(t=1,F))): print(`Bad input`): RETURN(FAIL): fi: F1:=t*diff(F,t): G1:=normal(subs(t=1,F1)): mu:=RatToPol(G1,s,n): if mu=FAIL then RETURN(FAIL): fi: F1:=1/t^coeff(mu,n,0)*subs(s=s/t^coeff(mu,n,1),F): gu:=[mu]: F1:=t*diff(F1,t): for i from 2 to N do F1:=t*diff(F1,t): G1:=normal(subs(t=1,F1)): mu:=RatToPol(G1,s,n): if mu=FAIL then RETURN(mu): else gu:=[op(gu),mu]: fi: od: gu: end: #AlphaArgf(F,s,t,n,N,ORD1): given a rational function F(s,t) #such that the coeff. of s^n is the prob. gen. function #of a discrete r.v. parametrized by n #finds the asymptotic alpha up to the N-th coefficient as n goes to infinity #up to order ORD1 #(starting with alpha_3 (the skewness)) . #For example, try: #AlphaArgf(1/(1-s*(1+t)/2),s,t,n,6,2); AlphaArgf:=proc(F,s,t,n,N,ORD1) local gu,d2,di,i,mu,z,i1,lu,lu1,lu2: gu:=AveAndMomsSrgf(F,s,t,n,N): if gu=FAIL then RETURN(gu): fi: d2:=degree(gu[2],n): mu:=[]: for i from 3 to N do di:=degree(gu[i],n): # mu:=[op(mu),coeff(gu[i],n,di)/coeff(gu[2],n,d2)^(i/2)*n^(di-d2*i/2)]: lu1:=normal(gu[i]/n^di): lu2:=normal(gu[2]/n^d2): lu:=lu1/lu2^(i/2): lu:=normal(subs(n=1/z,lu)): lu:=taylor(lu,z=0,ORD1+1): lu:=add(coeff(lu,z,i1)*z^i1,i1=0..ORD1): lu:=subs(z=1/n,lu)*n^(di-d2*i/2): mu:=[op(mu),lu]: od: mu: end: #GuessQPOL1Vk(L,n,k): Given a list L, a variable x, and #a positive integer k, tries to find #polynomials [P_1(x),P_2(x), ..., P_k(x)] # such that L[k*i+j]=P_j(i) for i=1..nops(L)/k . #For example, try: #GuessQPOL1Vk([1,2,1,4,1,6,1,8,1,10],n,2); GuessQPOL1Vk:=proc(L,n,k) local gu,j,L1,mu,n1,i1: gu:=[]: for j from 1 to k do L1:=[seq(L[k*i1+j],i1=1..trunc((nops(L)-j)/k))]: mu:=GuessPOL1V(L1,n1): if mu=FAIL then RETURN(FAIL): else gu:=[op(gu),expand(subs(n1=(n-j)/k,mu))]: fi: od: gu: end: #AveAndMomsSgfQuasi(F,x,N,n,Deg,z,k): Given an explicit (symbolic) #generating function F(x,z), such that the Maclaurin coeff. of z^n #is a prob. distribution, parametrized by the discrete variable #n, tries to find explicit quasi-polynomial expressions, in n #(of period k) #of degree <=Deg (by pure guessing) #for the first N moments #For example, try: #AveAndMomsSgfQuasi(1/(1-z*(1+x)/2),z,x,n,4,10,1); #AveAndMomsSgfQuasi(Feller(z,t,1,1,1),z,t,n,4,10,2); AveAndMomsSgfQuasi:=proc(F,z,x,n,N,Deg,k) local mu,gu,n0,i,j,lu,mu1,F1: F1:=taylor(F,z=0,Deg+10): mu:=[seq(AveAndMoms( expand(coeff(F1,z,n0)),x,N),n0=1..Deg+5)]: gu:=[]: for j from 1 to N do mu1:=[seq(mu[i][j],i=1..nops(mu))]: lu:=GuessQPOL1Vk(mu1,n,k): if lu=FAIL then print(`couldn't guess the`, j , `-th moment . `): RETURN(lu): else gu:=[op(gu),lu]: fi: od: gu: end: #AlphaAgf(F,s,t,n,N,ORD1): given a rational function F(s,t) #such that the coeff. of s^n is the prob. gen. function #of a discrete r.v. parametrized by n #finds the asymptotic alpha up to the N-th coefficients as n goes to infinity #to order ORD1 #(starting with alpha_3 (the skewness)) . #by guessing #For example, try: #AlphaAgf(1/(1-s*(1+t)/2),s,t,n,6,10,2); AlphaAgf:=proc(F,s,t,n,N,Deg,ORD1) local gu,d2,di,i,mu,lu,lu1,lu2,z,i1: gu:=AveAndMomsSgf(F,s,t,n,N,Deg): if gu=FAIL then RETURN(gu): fi: d2:=degree(gu[2],n): mu:=[]: for i from 3 to N do di:=degree(gu[i],n): # mu:=[op(mu),coeff(gu[i],n,di)/coeff(gu[2],n,d2)^(i/2)*n^(di-d2*i/2)]: lu1:=normal(gu[i]/n^di): lu2:=normal(gu[2]/n^d2): lu:=lu1/lu2^(i/2): lu:=normal(subs(n=1/z,lu)): lu:=taylor(lu,z=0,ORD1+1): lu:=add(coeff(lu,z,i1)*z^i1,i1=0..ORD1): lu:=subs(z=1/n,lu)*n^(di-d2*i/2): mu:=[op(mu),lu]: od: mu: end: #LQuasi(L,n): Given a quasi-polynomial L in n returns #the leading term, if it is consistent, otherwise returns FAIL #For example, try: #LQuasi([3*n+1,3*n+2],n); LQuasi:=proc(L,n) local k,deg1,i,coef1: k:=nops(L): deg1:={seq(degree(L[i],n),i=1..k)}: if nops(deg1)<>1 then RETURN(FAIL): fi: deg1:=deg1[1]: coef1:={seq(coeff(L[i],n,deg1),i=1..k)}: if nops(coef1)<>1 then RETURN(FAIL): fi: coef1[1]*n^deg1: end: #AveAndMomsSgfQuasiL(F,x,N,n,Deg,z,k): #the leading coefficients of AveAndMomsSgfQuasi(F,x,N,n,Deg,z,k) #if they are the same, otherwise, returns FAIL #For example, try: #AveAndMomsSgfQuasiL(1/(1-z*(1+x)/2),z,x,n,4,10,1); #AveAndMomsSgfQuasiL(Feller(z,t,1,1,1),z,t,n,4,10,2); AveAndMomsSgfQuasiL:=proc(F,z,x,n,N,Deg,k) local gu,gu1,mu1,i: gu:=AveAndMomsSgfQuasi(F,z,x,n,N,Deg,k): if gu=FAIL then RETURN(FAIL): fi: gu1:=[]: for i from 1 to nops(gu) do mu1:=LQuasi(gu[i],n): if mu1=FAIL then RETURN(FAIL): fi: gu1:=[op(gu1),mu1]: od: gu1: end: #AlphaAgfQuasi(F,s,t,n,N,k): given a rational function F(s,t) #such that the coeff. of s^n is the prob. gen. function #of a discrete r.v. parametrized by n #finds the asymptotic alpha up to the N-th coefficients as n goes to infinity #(starting with alpha_3 (the skewness)) . #by guessing a quasi-polynomial of period k #For example, try: #AlphaAgfQuasi(1/(1-s*(1+t)/2),s,t,n,6,10,1); #AlphaAgfQuasi(Feller(z,t,1,1,1),z,t,n,4,10,2); AlphaAgfQuasi:=proc(F,s,t,n,N,Deg,k) local gu,d2,di,i,mu: gu:=AveAndMomsSgfQuasiL(F,s,t,n,N,Deg,k): if gu=FAIL then RETURN(gu): fi: d2:=degree(gu[2],n): mu:=[]: for i from 3 to N do di:=degree(gu[i],n): mu:=[op(mu),coeff(gu[i],n,di)/coeff(gu[2],n,d2)^(i/2)*n^(di-d2*i/2)]: od: mu: end: ###Findrec #findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating #the sequence f of degree DEGREE and order ORDER #For example, try: findrec([seq(i,i=1..10)],0,2,n,N); findrecVerbose:=proc(f,DEGREE,ORDER,n,N) local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a: if (1+DEGREE)*(1+ORDER)+3+ORDER>nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): fi: ope:=0: var:={}: for i from 0 to ORDER do for j from 0 to DEGREE do ope:=ope+a[i,j]*n^j*N^i: var:=var union {a[i,j]}: od: od: eq:={}: for n0 from 1 to (1+DEGREE)*(1+ORDER)+2 do eq1:=0: for i from 0 to ORDER do eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f): od: eq:= eq union {eq1}: od: var1:=solve(eq,var): kv:={}: for i from 1 to nops(var1) do mu:=op(i,var1): if op(1,mu)=op(2,mu) then kv:= kv union {op(1,mu)}: fi: od: ope:=subs(var1,ope): if ope=0 then RETURN(FAIL): fi: ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}: if nops(ope)>1 then print(`There is some slack, there are `, nops(ope)): print(ope): RETURN(Yafe(ope[1],N)[2]): elif nops(ope)=1 then RETURN(Yafe(ope[1],N)[2]): else RETURN(FAIL): fi: end: #findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating #the sequence f of degree DEGREE and order ORDER #For example, try: findrec([seq(i,i=1..10)],0,2,n,N); findrec:=proc(f,DEGREE,ORDER,n,N) local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a: option remember: if not findrecEx(f,DEGREE,ORDER,ithprime(20)) then RETURN(FAIL): fi: if not findrecEx(f,DEGREE,ORDER,ithprime(40)) then RETURN(FAIL): fi: if not findrecEx(f,DEGREE,ORDER,ithprime(80)) then RETURN(FAIL): fi: if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): fi: ope:=0: var:={}: for i from 0 to ORDER do for j from 0 to DEGREE do ope:=ope+a[i,j]*n^j*N^i: var:=var union {a[i,j]}: od: od: eq:={}: for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do eq1:=0: for i from 0 to ORDER do eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f): od: eq:= eq union {eq1}: od: var1:=solve(eq,var): kv:={}: for i from 1 to nops(var1) do mu:=op(i,var1): if op(1,mu)=op(2,mu) then kv:= kv union {op(1,mu)}: fi: od: ope:=subs(var1,ope): if ope=0 then RETURN(FAIL): fi: ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}: if nops(ope)>1 then RETURN(Yafe(ope[1],N)[2]): elif nops(ope)=1 then RETURN(Yafe(ope[1],N)[2]): else RETURN(FAIL): fi: end: Yafe:=proc(ope,N) local i,ope1,coe1,L: if ope=0 then RETURN(1,0): fi: ope1:=expand(ope): L:=degree(ope1,N): coe1:=coeff(ope1,N,L): ope1:=normal(ope1/coe1): ope1:=normal(ope1): ope1:= convert( [seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`): factor(coe1),ope1: end: #Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with #poly coffs. #of maximum DEGREE+ORDER<=MaxC #e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2); Findrec:=proc(f,n,N,MaxC) local DEGREE, ORDER,ope,L: for L from 0 to MaxC do for ORDER from 0 to L do DEGREE:=L-ORDER: if (2+DEGREE)*(1+ORDER)+4>=nops(f) then print(`Insufficient data for degree`, DEGREE, `and order `,ORDER): RETURN(FAIL): fi: ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N): if ope<>FAIL then RETURN(ope): fi: od: od: FAIL: end: #SeqFromRec(ope,n,N,Ini,K): Given the first L-1 #terms of the sequence Ini=[f(1), ..., f(L-1)] #satisfied by the recurrence ope(n,N)f(n)=0 #extends it to the first K values SeqFromRec:=proc(ope,n,N,Ini,K) local ope1,gu,L,n1,j1: ope1:=Yafe(ope,N)[2]: L:=degree(ope1,N): if nops(Ini)<>L then ERROR(`Ini should be of length`, L): fi: ope1:=expand(subs(n=n-L,ope1)/N^L): gu:=Ini: for n1 from nops(Ini)+1 to K do gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)), j1=1..L)]: od: gu: end: #end Findrec with(linalg): #findrecEx(f,DEGREE,ORDER,m1): Explores whether thre #is a good chance that there is a recurrence of degree DEGREE #and order ORDER, using the prime m1 #For example, try: findrecEx([seq(i,i=1..10)],0,2,n,N,1003); findrecEx:=proc(f,DEGREE,ORDER,m1) local ope,var,eq,i,j,n0,eq1,a,A1, D1,E1,Eq,Var,f1,n,N: option remember: f1:=f mod m1: if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): fi: ope:=0: var:={}: for i from 0 to ORDER do for j from 0 to DEGREE do ope:=ope+a[i,j]*n^j*N^i: var:=var union {a[i,j]}: od: od: eq:={}: for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do eq1:=0: for i from 0 to ORDER do eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f1) mod m1: od: eq:= eq union {eq1}: od: Eq:= convert(eq,list): Var:= convert(var,list): D1:=nops(Var): E1:=nops(Eq): if E10 then RETURN(false): fi: if E1-D1>=1 then for j from 1 to nops(Var) do A1[D1,j]:=coeff(Eq[D1+1],Var[j]): od: if det(A1) mod m1 <>0 then RETURN(false): fi: fi: if E1-D1>=2 then for j from 1 to nops(Var) do A1[D1,j]:=coeff(Eq[D1+2],Var[j]): od: if det(A1) mod m1 <>0 then RETURN(false): fi: fi: true: end: #######end Findrec #AveAndMomsSseq(F,x,N,n,Deg): #Given a sequence F #whose n-th item is the #prob. distribution, parametrized by the discrete variable #n, tries to find explicit polynomial expressions, in n #of degree <=Deg #for the first N moments #For example, try: #AveAndMomsSseq([seq(((1+x)/2)^i,i=1..20)],x,4,n,10); AveAndMomsSseq:=proc(F,x,N,n,Deg) local mu,gu,n0,i,j,lu,mu1: if not type(F,list) then print(`first argument should be a list`): RETURN(FAIL): fi: if nops(F)[0$nops(guo)] then print(`Not asymptotically normal`): RETURN(FAIL): fi: gue:=[seq(limit(L[2*i],n=infinity),i=1..trunc(nops(L)/2))]: gue1:=[seq((2*r1)!/r1!/2^r1,r1=2..nops(gue)+1)]: if gue<>gue1 then print(`Not asymptotically normal`): RETURN(FAIL): fi: Le:=[seq(L[2*i]/gue[i],i=1..trunc(nops(L)/2))]: Lo:=[seq(L[2*i-1]/gue[i],i=1..trunc(nops(L))/2)]: vu:=1: for i from 1 to ORD1 do Le1:=[seq(coeff(Le[i1],n,-i),i1=1..nops(Le))]: pol:=GuessPOL1V(Le1,r): if pol=FAIL then print(`Insufficient data to go beyond order`, i-1): zugi:=vu: break: fi: pol:=expand(subs(r=r-1,pol)): vu:=vu+factor(pol)/n^i: od: zugi:=(2*r)!/r!/2^r*vu: ka1:={seq(eval(subs(r=i1+1,zugi)-L[2*i1]),i1=1..trunc(nops(L)/2))}: if ka1<>{0} then print(ka1): print(zugi, `is a lemon`): RETURN(FAIL): fi: if Lo=[0$(nops(Lo))] then RETURN(gu[1],gu[2],zugi,0): fi: ka:=coeff(normal(Lo[1]*sqrt(n)),n,0)/sqrt(n): Lo:=expand([seq(normal(Lo[i]/ka),i=1..nops(Lo))]): vu:=0: for i from 0 to ORD1 do Lo1:=[seq(coeff(Lo[i1],n,-i),i1=1..nops(Lo))]: pol:=GuessPOL1V(Lo1,r): if pol=FAIL then print(`Insufficient data to go beyond order`, i-1): peret:=vu: break: fi: pol:=expand(subs(r=r-1,pol)): vu:=vu+factor(pol)/n^i: od: peret:=(2*r)!/r!/2^r*ka*vu: ka1:={seq(normal(eval(subs(r=i1,peret)-L[2*i1-3])), i1=2..trunc((nops(L)+3)/2))}: if ka1<>{0} then print(ka1): print(peret, `is a lemon`): RETURN(zugi,FAIL): fi: gu[1],gu[2],zugi,peret: end: #YofR(f,s,t,N,n,ORD1,r): Given the inputs #(i) a rational function f of the variables s and t #such that the Maclaurin coeff. of s^n is the prob. #generating function for a discrete r.v. parametrized by n #(ii) variables s and t for representing f #(iii) a pos. integer N for the number of moments to use for guessing #(iv) a symbol n for (i) #(v) a pos. integer Deg for the max. degree of polynomials in n to be used #(vi) a pos. intger ORD1, for the order of the asymptotic expansions #(vii) a symbol r #outputs: #(i) a polynomial expression for the average in terms of n #(ii) a polynomial expression for the variance in terms of n #(iii) the first ORD1 terms in the asymptotic expansion of the #even (2r)-th alpha-coeff. (m_{2r}/m_2^r) #(iv) the first ORD1 terms in the asymptotic expansion of the #even (2r+1)-th alpha-coeff. (m_{2r+1}/m_2^(r+1/2)) #if it is not asymptotically normal, it returns FAIL. #For example try: #YofR(1/(1-(1+t)*s/2),s,t,12,n,1,r); YofR:=proc(f,s,t,N,n,ORD1,r) local guo,gue,gue1, i,r1,Le,Lo,Le1,pol, zugi,peret,ka,i1,vu,ka1,Lo1,L,gu: if nargs<7 then ERROR(`Needs 7 arguments`): fi: gu:=AveAndMomsSrgf(f,s,t,n,2): L:=AlphaArgf(f,s,t,n,N,ORD1): guo:=[seq(limit(L[2*i-1],n=infinity),i=1..trunc(nops(L)/2))]: if guo<>[0$nops(guo)] then print(`Not asymptotically normal`): RETURN(FAIL): fi: gue:=[seq(limit(L[2*i],n=infinity),i=1..trunc(nops(L)/2))]: gue1:=[seq((2*r1)!/r1!/2^r1,r1=2..nops(gue)+1)]: if gue<>gue1 then print(`Not asymptotically normal`): RETURN(FAIL): fi: Le:=[seq(L[2*i]/gue[i],i=1..trunc(nops(L)/2))]: Lo:=[seq(L[2*i-1]/gue[i],i=1..trunc(nops(L))/2)]: vu:=1: for i from 1 to ORD1 do Le1:=[seq(coeff(Le[i1],n,-i),i1=1..nops(Le))]: pol:=GuessPOL1V(Le1,r): if pol=FAIL then print(`Insufficient data to go beyond order`, i-1): zugi:=vu: break: fi: pol:=expand(subs(r=r-1,pol)): vu:=vu+factor(pol)/n^i: od: zugi:=(2*r)!/r!/2^r*vu: ka1:={seq(eval(subs(r=i1+1,zugi)-L[2*i1]),i1=1..trunc(nops(L)/2))}: if ka1<>{0} then print(ka1): print(zugi, `is a lemon`): RETURN(FAIL): fi: if Lo=[0$(nops(Lo))] then RETURN(gu[1],gu[2],zugi,0): fi: ka:=coeff(normal(Lo[1]*sqrt(n)),n,0)/sqrt(n): Lo:=expand([seq(normal(Lo[i]/ka),i=1..nops(Lo))]): vu:=0: for i from 0 to ORD1 do Lo1:=[seq(coeff(Lo[i1],n,-i),i1=1..nops(Lo))]: pol:=GuessPOL1V(Lo1,r): if pol=FAIL then print(`Insufficient data to go beyond order`, i-1): peret:=vu: break: fi: pol:=expand(subs(r=r-1,pol)): vu:=vu+factor(pol)/n^i: od: peret:=(2*r)!/r!/2^r*ka*vu: ka1:={seq(normal(eval(subs(r=i1,peret)-L[2*i1-3])), i1=2..trunc((nops(L)+3)/2))}: if ka1<>{0} then print(ka1): print(peret, `is a lemon`): RETURN(zugi,FAIL): fi: gu[1],gu[2],zugi,peret: end: ### #FaM(F,x,N): Given an explicit (symbolic) #expression F(x,n), denoting a sequence of discrete #prob. distribution, parametrized by the discrete variable #n, finds the average followed by the 2nd- through the #N-th factorial moment (about the mean) #For example, try: #FaM(((1+x)/2)^n,x,4); FaM:=proc(F,x,N) local F1,ku,y,ave1,i: ku:=subs(x=1,F): if ku=0 then RETURN(FAIL): fi: F1:=F/ku: #ave1:=subs(x=1,normal(simplify(diff(F1,x)))): ave1:=limit(normal(simplify(diff(F1,x))),x=1): F1:=F1/x^ave1: F1:=subs(x=1+y,F1): F1:=series(F1,y=0,N+1): [ave1,seq(simplify(expand(i!*coeff(F1,y,i))),i=2..N)]: end: #AveAndMomsD(F,x,N): Given an explicit (symbolic) #expression F(x,n), denoting a sequence of discrete #prob. distribution, parametrized by the discrete variable #n, #for the first N moments #For example, try: #AveAndMomsD(((1+x)/2)^n,x,4); AveAndMomsD:=proc(F,x,N) local ave1,gu,i,j: gu:=FaM(F,x,N): ave1:=gu[1]: gu:=[0,op(2..nops(gu),gu)]: [ave1, seq(simplify(add(stirling2(i,j)*gu[j],j=1..i)),i=2..N)]: end: #YofD(F,x,N,n,ORD1,r): Given the inputs #(i) a discrete prob. generating function #parametrized by n phrased in terms of the variable x #(ii) a symbol x for the variable name #(iii) a pos. integer N for the number of moments to use for guessing #(iv) a symbol n for (i) #(v) a pos. intger ORD1, for the order of the asymptotic expansions #(vi) a symbol r #outputs: #(i) a polynomial expression for the average in terms of n #(ii) a polynomial expression for the variance in terms of n #(iii) the first ORD1 terms in the asymptotic expansion of the #even (2r)-th alpha-coeff. (m_{2r}/m_2^r) #(iv) the first ORD1 terms in the asymptotic expansion of the #even (2r+1)-th alpha-coeff. (m_{2r+1}/m_2^(r+1/2)) #if it is not asymptotically normal, it returns FAIL. #For example try: #YofD(((1+x)/2)^n,x,16,n,2,r); YofD:=proc(F,x,N,n,ORD1,r) local guo,gue,gue1, i,r1,Le,Lo,Le1,pol, zugi,peret,ka,i1,vu,ka1,Lo1,L,gu: gu:=AveAndMomsD(F,x,2): L:=AlphaD(F,x,N,n,ORD1): guo:=[seq(limit(L[2*i-1],n=infinity),i=1..trunc(nops(L)/2))]: if guo<>[0$nops(guo)] then print(`Not asymptotically normal`): RETURN(FAIL): fi: gue:=[seq(limit(L[2*i],n=infinity),i=1..trunc(nops(L)/2))]: gue1:=[seq((2*r1)!/r1!/2^r1,r1=2..nops(gue)+1)]: if gue<>gue1 then print(`Not asymptotically normal`): RETURN(FAIL): fi: Le:=[seq(L[2*i]/gue[i],i=1..trunc(nops(L)/2))]: Lo:=[seq(L[2*i-1]/gue[i],i=1..trunc(nops(L))/2)]: vu:=1: for i from 1 to ORD1 do Le1:=[seq(coeff(Le[i1],n,-i),i1=1..nops(Le))]: pol:=GuessPOL1V(Le1,r): if pol=FAIL then print(`Insufficient data to go beyond order`, i-1): zugi:=vu: break: fi: pol:=expand(subs(r=r-1,pol)): vu:=vu+factor(pol)/n^i: od: zugi:=(2*r)!/r!/2^r*vu: ka1:={seq(eval(subs(r=i1+1,zugi)-L[2*i1]),i1=1..trunc(nops(L)/2))}: if ka1<>{0} then print(ka1): print(zugi, `is a lemon`): RETURN(FAIL): fi: if Lo=[0$(nops(Lo))] then RETURN(gu[1],gu[2],zugi,0): fi: ka:=coeff(normal(Lo[1]*sqrt(n)),n,0)/sqrt(n): Lo:=expand([seq(normal(Lo[i]/ka),i=1..nops(Lo))]): vu:=0: for i from 0 to ORD1 do Lo1:=[seq(coeff(Lo[i1],n,-i),i1=1..nops(Lo))]: pol:=GuessPOL1V(Lo1,r): if pol=FAIL then print(`Insufficient data to go beyond order`, i-1): peret:=vu: break: fi: pol:=expand(subs(r=r-1,pol)): vu:=vu+factor(pol)/n^i: od: peret:=(2*r)!/r!/2^r*ka*vu: ka1:={seq(normal(eval(subs(r=i1,peret)-L[2*i1-3])), i1=2..trunc((nops(L)+3)/2))}: if ka1<>{0} then print(ka1): print(peret, `is a lemon`): RETURN(zugi,FAIL): fi: simplify(gu[1]),simplify(gu[2]),zugi,peret: end: #FaMrgf(F,s,x,n,N): Given a rational generating function #F(s,x), such that viewed as a formal-power-series #in s, the coeff. of s^n is the prob. generating function #ofa sequence of discrete #prob. distribution, parametrized by the discrete variable #n, finds the (linear) #expression (in n) for the average followed by the 2nd- through the #N-th factorial moment (about the mean) #For example, try: #FaMrgf(1/(1-(1+x)/2*s),s,x,n,4); FaMrgf:=proc(F,s,x,n,N) local F1,ku,y,mu,i,G1: ku:=normal(simplify(subs(x=1,F))): if normal(ku-1/(1-s))<>0 then print(`Bad input`): RETURN(FAIL): fi: F1:=x*diff(F,x): G1:=normal(subs(x=1,F1)): mu:=RatToPol(G1,s,n): if mu=FAIL then RETURN(FAIL): fi: F1:=1/x^coeff(mu,n,0)*subs(s=s/x^coeff(mu,n,1),F): F1:=subs(x=1+y,F1): F1:=taylor(F1,y=0,N+1): [mu,seq(RatToPol(expand(i!*coeff(F1,y,i)),s,n),i=2..N)]: end: #AveAndMomsSrgf(F,s,t,n,N): given a rational function F(s,t) #such that the coeff. of s^n is the prob. gen. function #of a discrete r.v. parametrized by n #finds polynomial expressions in n for the first N moments. #For example, try: #AveAndMomsSrgf(1/(1-s*(1+t)/2),s,t,n,6); AveAndMomsSrgf:=proc(F,s,t,n,N) local gu,ave1,i,j: gu:=FaMrgf(F,s,t,n,N): ave1:=gu[1]: gu:=[0,op(2..nops(gu),gu)]: [ave1, seq(add(stirling2(i,j)*gu[j],j=1..i),i=2..N)]: end: ####Start Bivariate######################## #FaM2(F,x1,x2,N): Given an explicit (symbolic) #expression F(x1,x2,n), denoting a sequence of discrete #prob. distribution, parametrized by the discrete variable #n, finds the averages of the r.vs corresponding to x1 and x2 #followed by the list-of-lists whose L such that L[r][s] #1<=r,s<=N is the #(r-1,s-1)-th mixed factorial moment (about the means) #For example, try: #FaM2(((1+x1+x2)/3)^n,x1,x2,4); FaM2:=proc(F,x1,x2,N) local F1,ku,y1,y2,ave1,ave2,gu,r,s,F11: ku:=subs({x1=1,x2=1},F): if ku=0 then RETURN(FAIL): fi: F1:=F/ku: ave1:=limit(normal(simplify(diff(subs(x2=1,F1),x1))),x1=1): ave2:=limit(normal(simplify(diff(subs(x1=1,F1),x2))),x2=1): F1:=F1/x1^ave1/x2^ave2: F1:=subs({x1=1+y1,x2=1+y2},F1): gu:=[]: F1:=series(F1,y1=0,N+1): for r from 0 to N do F11:=r!*coeff(F1,y1,r): F11:=series(F11,y2=0,N+1): gu:=[op(gu),[seq(expand(s!*coeff(F11,y2,s)),s=0..N)]]: od: [ave1,ave2,gu]: end: #AveAndMomsD2(F,x1,x2,N): Given an explicit (symbolic) #expression F(x1,x2,n), denoting a sequence of bi-variate discrete #prob. distribution, parametrized by the discrete variable #n, using the variables x1, x2, returns the #averages followed by a list of lists L such that #L[r+1][s+1] is the (r,s)-mixed moment (about the mean) #for 0<=r,s<=N #For example, try: #AveAndMomsD2(((1+x1+x2)/3)^n,x1,x2,4); AveAndMomsD2:=proc(F,x1,x2,N) local ave1,ave2,gu,i1,i2,j1,j2,T: gu:=FaM2(F,x1,x2,N): ave1:=gu[1]: ave2:=gu[2]: gu:=gu[3]: for i1 from 0 to N do for i2 from 0 to N do T[i1,i2]:= add(add(stirling2(i1,j1)*stirling2(i2,j2)*gu[j1+1][j2+1], j1=0..N),j2=0..N): od: od: [ave1,ave2,[seq([seq(T[i1,i2],i2=0..N)],i1=0..N)]]: end: #AveAndMoms2(f,x1,x2,N): Given a specific #probability generating function #for a bivariate distribution #f(x1,x2) (a polynomial, rational function and possibly beyond) #returns the averages of the first and second r.v. followed #by a list-of-lists L #such that L[r+1][s+1] is the (r,s)-mixed moment about the mean #For example, try: #AveAndMoms2(((1+x1+x2)/3)^100,x1,x2,4); AveAndMoms2:=proc(f,x1,x2,N) local mu,F,memu1,memu2,i1,i2,T,F1,F11: mu:=simplify(subs({x1=1,x2=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({x1=1,x2=1},diff(F,x1))): memu2:=simplify(subs({x1=1,x2=1},diff(F,x2))): F:=F/x1^memu1/x2^memu2: F1:=F: for i1 from 0 to N do F11:=F1: for i2 from 0 to N do T[i1,i2]:=subs({x1=1,x2=1},F11): F11:=normal(x2*diff(F11,x2)): od: F1:=normal(x1*diff(F1,x1)): od: [memu1,memu2,[seq([seq(T[i1,i2],i2=0..N)],i1=0..N)]]: end: #AveAndMomsS2(F,x1,x2,N,n,Deg): Given an explicit (symbolic) #expression F1(x1,x2,n), denoting a sequence of bivariate discrete #prob. distributions, parametrized by the discrete variable #n, tries to find explicit polynomial expressions, in n #of degree <=Deg #for the first N moments #For example, try: #AveAndMomsS2(((1+x1+x2)/3)^n,x1,x2,4,n,10); AveAndMomsS2:=proc(F,x1,x2,N,n,Deg) local mu,n0,i,j1,j2,lu,mu1,ave1,ave2,T: mu:=[seq(AveAndMoms2(simplify(subs(n=n0,F)),x1,x2,N),n0=1..Deg+5)]: ave1:=GuessPOL1V([seq(mu[n0][1],n0=1..Deg+5)],n): if ave1=FAIL then RETURN(FAIL): fi: ave2:=GuessPOL1V([seq(mu[n0][2],n0=1..Deg+5)],n): if ave2=FAIL then RETURN(FAIL): fi: mu:=[seq(mu[n0][3],n0=1..Deg+5)]: for j1 from 0 to N do for j2 from 0 to N do mu1:=[seq(mu[i][j1+1][j2+1],i=1..nops(mu))]: lu:=GuessPOL1V(mu1,n): if lu=FAIL then print(`couldn't guess the`, [j1,j2] , `-th moment . `): REETURN(FAIL): else T[j1,j2]:=lu: fi: od: od: [ave1,ave2, [seq([seq(T[j1,j2],j2=0..N)],j1=0..N)] ]: end: #AveAndMomsSseq2(L,x1,x2,N,n,Deg): Given a list L of bi-variate #prob. distributions #F(x1,x2,n), denoting a sequence of bivariate discrete #prob. distributions, parametrized by the discrete variable #n, tries to find explicit polynomial expressions, in n #of degree <=Deg #for the first N moments #For example, try: #AveAndMomsSseq2([seq(((1+x1+x2)/3)^n1,n1=1..20)],,x1,x2,4,n,10); AveAndMomsSseq2:=proc(L,x1,x2,N,n,Deg) local mu,n0,i,j1,j2,lu,mu1,ave1,ave2,T: mu:=[seq(AveAndMoms2(L[i],x1,x2,N),i=1..nops(L))]: if nops(mu)0 then print(`Bad input`): RETURN(ku,FAIL): fi: F1:=x1*diff(F,x1): G1:=normal(subs({x1=1,x2=1},F1)): mu1:=RatToPol(G1,s,n): if mu1=FAIL then RETURN(FAIL): fi: F2:=x2*diff(F,x2): G2:=normal(subs({x1=1,x2=1},F2)): mu2:=RatToPol(G2,s,n): if mu2=FAIL then RETURN(FAIL): fi: F12:=1/x1^coeff(mu1,n,0)/x2^coeff(mu2,n,0) *subs({s=s/x1^coeff(mu1,n,1)/x2^coeff(mu2,n,1)},F): F12:=subs({x1=1+y1,x2=1+y2},F12): F12:=taylor(F12,y1=0,N+1): for i from 0 to N do F12a:=coeff(F12,y1,i): F12a:=taylor(F12a,y2=0,N+1): for j from 0 to N do ka:=normal(coeff(F12a,y2,j)): ka:=RatToPol(expand(i!*j!*ka),s,n): T[i,j]:=ka: od: od: [mu1,mu2,[seq([seq(T[i,j],j=0..N)],i=0..N)]]: end: #AveAndMomsSrgf2(F,s,x1,x2,n,N): Given a rational generating function #F(s,x1,x2), such that viewed as a formal-power-series #in s, the coeff. of s^n is the bi-variate prob. generating function #ofa sequence of discrete #prob. distribution, parametrized by the discrete variable #n, in x1, x2,finds the (linear) #expression (in n) for the average followed by #the list of lists for mixed moments (about the mean) #For example, try: #AveAndMomsSrgf2(1/(1-(1+x1+x2)/3*s),s,x1,x2,n,4); AveAndMomsSrgf2:=proc(F,s,x1,x2,n,N) local ave1,ave2,gu,i1,i2,j1,j2,T: gu:=FaMrgf2(F,s,x1,x2,n,N): ave1:=gu[1]: ave2:=gu[2]: gu:=gu[3]: for i1 from 0 to N do for i2 from 0 to N do T[i1,i2]:= add(add(stirling2(i1,j1)*stirling2(i2,j2)*gu[j1+1][j2+1], j1=0..N),j2=0..N): od: od: [ave1,ave2,[seq([seq(T[i1,i2],i2=0..N)],i1=0..N)]]: end: #AlphaD2(F,x1,x2,N,n,ORD1): #Given an expression F(n,x1,x2) #whose n-th item is a bi-variate #prob. distribution, parametrized by the discrete variable #n , depending on the continuos variables x1 and x2 #tries to find asymptotic expressions in n, of order ORD1 #using guessed expressions for the normalized mixed-moments #moments of degree <=Deg #for the first N^2 mixed moments #For example, try: #AlphaD2((1+x1+x2)/3)^n,x1,x2,4,n,2); AlphaD2:=proc(F,x1,x2,N,n,ORD1) local mu,d1,d2,di,i,z,i1,lu,lu1,lu2a,lu2b,T,j,ave1,ave2,var1,var2: mu:=AveAndMomsD2(F,x1,x2,N): ave1:=mu[1]: ave2:=mu[2]: mu:=mu[3]: var1:=mu[1][3]: var2:=mu[3][1]: if mu=FAIL then RETURN(FAIL): fi: d1:=degree(mu[3][1],n): d2:=degree(mu[1][3],n): for i from 0 to N do for j from 0 to N do if mu[i+1][j+1]=0 then T[i,j]:=0: else di:=degree(mu[i+1][j+1],n): lu1:=normal(mu[i+1][j+1]/n^di): lu2a:=normal(mu[3][1]/n^d1): lu2b:=normal(mu[1][3]/n^d2): lu:=lu1/lu2a^(i/2)/lu2b^(j/2): lu:=normal(subs(n=1/z,lu)): lu:=taylor(lu,z=0,ORD1+1): lu:=add(coeff(lu,z,i1)*z^i1,i1=0..ORD1): lu:=subs(z=1/n,lu)*n^(di-d1*i/2-d2*j/2): T[i,j]:=lu: fi: od: od: [ave1,var1],[ave2,var2],[seq([seq(T[i,j],j=0..N)],i=0..N)]: end: #Argf2(F,s,x1,x2,N,n,ORD1) #Given an expression F(n,x1,x2) #whose n-th item is a bi-variate #prob. distribution, parametrized by the discrete variable #n , depending on the continuos variables x1 and x2 #tries to find asymptotic expressions in n, of order ORD1 #using guessed expressions for the normalized mixed-moments #moments of degree <=Deg #for the first N^2 mixed moments #For example, try: #Argf2(1/(1-(1+x1+x2)/3*s),s,x1,x2,4,n,2); Argf2:=proc(F,s,x1,x2,N,n,ORD1) local mu,d1,d2,di,i,z,i1,lu,lu1,lu2a,lu2b,T,j,ave1,ave2,var1,var2: mu:=AveAndMomsSrgf2(F,s,x1,x2,n,N): ave1:=mu[1]: ave2:=mu[2]: mu:=mu[3]: var1:=mu[1][3]: var2:=mu[3][1]: if mu=FAIL then RETURN(FAIL): fi: d1:=degree(mu[3][1],n): d2:=degree(mu[1][3],n): for i from 0 to N do for j from 0 to N do if mu[i+1][j+1]=0 then T[i,j]:=0: else di:=degree(mu[i+1][j+1],n): lu1:=normal(mu[i+1][j+1]/n^di): lu2a:=normal(mu[3][1]/n^d1): lu2b:=normal(mu[1][3]/n^d2): lu:=lu1/lu2a^(i/2)/lu2b^(j/2): lu:=normal(subs(n=1/z,lu)): lu:=taylor(lu,z=0,ORD1+1): lu:=add(coeff(lu,z,i1)*z^i1,i1=0..ORD1): lu:=subs(z=1/n,lu)*n^(di-d1*i/2-d2*j/2): T[i,j]:=lu: fi: od: od: [ave1,var1],[ave2,var2],[seq([seq(T[i,j],j=0..N)],i=0..N)]: end: #TwoWordsStat(alphab,w1,w2,N,n,ORD1): given two words of the same #length w1,w2, in the alphabet alphab, #finds the average, variance (as expressions in n) and the #asymptotic expressions (in n), to order ORD1, of #normalized mixed-moments (up to the Nth) of the random variables #number of occurences of w1 and number of occurrences of w2 #(as factors) in a (uniform-at-random) n-letter word #For example, try: #TwoWordsStat({1,2},[1,2],[2,1],4,n,1); TwoWordsStat:=proc(alphab,w1,w2,N,n,ORD1) local gu,s,t: if nops(w1)<>nops(w2) then print(w1,w2, `should be of the same size, but they are not`): RETURN(FAIL): fi: if convert(w1,set) minus alphab<>{} then print(w1, `is not a word in the alphabet`, alphab): RETURN(FAIL): fi: if convert(w2,set) minus alphab<>{} then print(w2, `is not a word in the alphabet`, alphab): RETURN(FAIL): fi: gu:=GJst(alphab,{w1,w2},s,t): Argf2(gu,s,t[op(w1)],t[op(w2)],N,n,ORD1): end: #Milim(r,k): The list of all k-letter words in the alphabet {1,2,...,r} #for example, try: #Milim(2,3); Milim:=proc(r,k) local i,lu,w: option remember: if k=0 then RETURN({[]}): fi: lu:=Milim(r,k-1): [seq(seq([i,op(w)], w in lu),i=1..r)]: end: #CorContestW(r,k):Finds the pair of words that are most correlated #and least correlated among k-letter words in {1,...,r}. #For example, try: #CorContestW(2,3); CorContestW:=proc(r,k) local gu,w1,w2,alufG,alufK,siG,siK,alphab,i,halev,n,j: alphab:={seq(i,i=1..r)}: gu:=Milim(r,k): w1:=gu[1]: w2:=gu[2]: alufG:={[w1,w2]}; alufK:={[w1,w2]}: siG:=coeff(TwoWordsStat(alphab,w1,w2,2,n,0)[3][2][2],n,0): siK:=coeff(TwoWordsStat(alphab,w1,w2,2,n,0)[3][2][2],n,0): for i from 1 to nops(gu) do w1:=gu[i]: for j from i+1 to nops(gu) do w2:=gu[j]: halev:=coeff(TwoWordsStat(alphab,w1,w2,2,n,0)[3][2][2],n,0): if evalf(halev)>evalf(siG) then siG:=halev: alufG:={[w1,w2]}: elif halev=siG then alufG:=alufG union {[w1,w2]}: fi: if evalf(halev)evalf(siG) then siG:=halev: alufG:={[w1,w2]}: elif halev=siG then alufG:=alufG union {[w1,w2]}: fi: if evalf(halev)FAIL then mu:=expand(subs(x=x-s+1,mu)): for s from nops(L) to 1 by -1 while subs(x=s,mu)=L[s] do od: RETURN(expand(mu),s+1): fi: od: FAIL: end: #Faxy(a,x,y): the normalized bi-variate distributin #of two normal r.v.s whose correlation is a Faxy:=proc(a,x,y) exp(-x^2/2-y^2/2+x*y*a)/2/sqrt(2)/Pi*sqrt(2-2*a^2): end: #GCor(a,r,s): the (r,s)-moment of Faxy(a,x,y) divided by the s.d. #to the power (r+s)/2. For example, try: #GCor(1/3,2,2); GCor:=proc(a,r,s) local x,y: normal(int(int(x^r*y^s*Faxy(a,x,y), x=-infinity..infinity),y=-infinity..infinity) *(1-a^2)^((r+s)/2)): end: #CheckANrgf2(F,s,x1,x2,N) #Given an expression F(n,x1,x2) #whose n-th item is a bi-variate #prob. distribution, parametrized by the discrete variable #n , depending on the continuos variables x1 and x2 #checks whether it seems to be joint-asymptotically #normal with the right correlation. #For example, try: #CheckANrgf2(1/(1-(1+x1+x2)/3*s),s,x1,x2,4); CheckANrgf2:=proc(F,s,x1,x2,N) local gu,n,a,r1,r2: gu:=Argf2(F,s,x1,x2,N,n,0)[3]: if gu[1][1]<>1 then RETURN(FAIL): fi: a:=coeff(gu[2][2],n,0): if a=1 then print(`Correlation is 1`): RETURN(FAIL): fi: evalb( {seq(seq(limit(gu[1+r1][1+r2],n=infinity)-GCor(a,r1,r2),r1=1..N),r2=1..N)} ={0}): end: #CheckANrgf2Plus(F,s,x1,x2,N,r1,r2,ORD1) #Given an expression F(n,x1,x2) #whose n-th item is a bi-variate #prob. distribution, parametrized by the discrete variable #n , depending on the continuos variables x1 and x2 #checks whether it seems to be joint-asymptotically #normal with the right correlation. #and also guesses polynomial expressions #for the mixed-moments #For example, try: #CheckANrgf2Plus(1/(1-(1+x1+x2)/3*s),s,x1,x2,8,r1,r2,1); CheckANrgf2Plus:=proc(F,s,x1,x2,N,r1,r2,ORD1) local gu,a,i1,n, i2,luee,lu1,mu,luee1: gu:=Argf2(F,s,x1,x2,N,n,ORD1)[3]: if gu[1][1]<>1 then RETURN(FAIL): fi: luee:=[]: for i1 from 0 to (nops(gu)-1)/2 do lu1:=[]: for i2 from 0 to (nops(gu[2*i1+1])-1)/2 do mu:=gu[2*i1+1][2*i2+1]: a:=limit(mu,n=infinity): if a=0 then RETURN(FAIL): fi: mu:=expand(normal(mu/a)): lu1:=[op(lu1),mu]: od: luee:=[op(luee),lu1]: od: luee1:= [seq([seq(coeff(luee[i1][i2],n,-1),i2=1..nops(luee[i1]))],i1=1..nops(luee))]: end: #GuessPOLtight(L,VarList): Tries to guess a polynomial in #VarList from the data L that is a list of lists of lists... #such that P(i1,i2,i3, ...)=L[i1][i2][i3] ... #For example, try: #GuessPOLtight([seq([seq(i+j,i=1..10)],j=1..10)],[x,y]); GuessPOLtight:=proc(L,VarList) local gu,i,VarList1,mu: if nops(VarList)=1 then RETURN(GuessPOL1Vtight(L,VarList[1])): fi: gu:=[]: VarList1:=[op(2..nops(VarList),VarList)]: for i from 1 to nops(L) do mu:=GuessPOLtight(L[i],VarList1): if mu=FAIL then RETURN(FAIL): else gu:=[op(gu),mu]: fi: od: GuessPOL1Vtight(gu,VarList[1]): end: #GuessPOL1Vtight(L,x): Given a list L, tries to find #a polynomial P(x) such that L[i]=P(i) for i=1..nops(L) . #For example, try: #GuessPOL1Vtight([seq(i^2,i=1..10)],x); GuessPOL1Vtight:=proc(L,x) local d,gu,i: for d from 0 to nops(L)-2 do gu:=GuessPOL1V1tight(L,x,d): if gu<>FAIL then if {seq(expand(subs(x=i,gu)-L[i]),i=1..nops(L))}={0} then RETURN(gu): fi: fi: od: FAIL: end: #GuessPOL1V1tight(L,x,d): Given a list L, tries to find #a polynomial P(x) of degree d such that L[i]=P(i) for i=1..nops(L) #For example, try: #GuessPOL1V1tight([seq(i^2,i=1..10)],x,2); GuessPOL1V1tight:=proc(L,x,d) local eq,var,a,P,i: if nops(L)-d<2 then print(`Insufficient data`): RETURN(FAIL): fi: P:=add(a[i]*x^i,i=0..d): var:={seq(a[i],i=0..d)}: eq:={seq(expand(subs(x=i,P)-L[i]),i=1..nops(L))}: var:=solve(eq,var): if var=NULL then RETURN(FAIL): fi: subs(var,P): end: #CheckZ80(n,t,q1,q2): checks Zeilberger's recurrence for Fn(t,q1,q2) from #his 1980 article for n. For example, try: #CheckZ80(4,t,q1,q2); CheckZ80:=proc(n,t,q1,q2) local k,i: expand( t*q2^n*Fn(n,t,q1,q2)- add( (-1)^(n-k)*Fn(k,t,q1,q2)*qbinPure(n,k,q1)*mul(1-t*q2^i,i=k+1..n)* q1^binomial(n-k,2), k=0..n ) ): end: #FnSlow(n,t,q1,q2): Fn(n,t,q1,q2) using #Doron Zeilberger's 1980 recurrence #(bottom of p. 198 of http://www.math.rutgers.edu/~zeilberg/mamarimY/ #Zeilberger_y1980_p192.pdf with a typo corrected #(there should be p^((n-k)*(n-k-1)/2) inside the summand #For example, try: #FnSlow(4,t,q1,q2); FnSlow:=proc(n,t,q1,q2) local k,i: option remember: if n<0 then RETURN(0): fi: if n=0 then RETURN(1): fi: normal( expand( add( (-1)^(n-k)*FnSlow(k,t,q1,q2)*qbinPure(n,k,q1)*mul(1-t*q2^i,i=k+1..n)* q1^binomial(n-k,2), k=0..n-1 )) /(t*q2^n-1) ): end: #AnalyseMoms2(L,n,ORD1,r,s): Given a list of lists L of polynomials #in n, such that L[i][j] equals the mixed moment M[i-1,j-1] #and we have that M[2*r,2*s+1]=0 and M[2*r+1,2*s]=0 #and the dist. is believed to be #a bivariate distribution that is joint-independently #asympotically normal, finds the asymptotics, #of expressions, to order ORD1 of M[2*r,2*s] and M[2*r+1,2*s+1] #The output is #(i) the variance of the first r.v. #(ii) the variance of the second r.v. #(iii) the covariance #(iv) the asympotic to order ORD1 of the [2*r,2*s] moment #(v) the asympotic to order ORD1 of the [2*r+1,2*s+1] moment #(vi) the asymptotics to order ORD1 of the normalized moments [2*r,2*s] #(vii) the asymptotics to order ORD1 of the normalized moments [2*r-1,2*s-1] #For example, do: #read oHISTABRUT20: #AnalyseMoms2(InvMajMoms10,n,1,r,s); AnalyseMoms2:=proc(L,n,ORD1,r,s) local d1,d2,c1,c2,lu,i,gue,guo,r1,s1, gue1,guo1,lu2a,lu2b,Gue,Guo,x: d1:=degree(L[1][3],n): d2:=degree(L[3][1],n): c1:=coeff(L[1][3],n,d1): c2:=coeff(L[3][1],n,d2): if {seq(seq(degree(L[2*r1+1][2*s1+1],n)-(d1*r1+d2*s1),r1=1..nops(L)/2-1), s1=1..nops(L)/2-1)}<>{0} then print(`Something is wrong`): RETURN(FAIL): fi: if {seq(seq(degree(L[2*r1][2*s1],n)-(d1*r1+d2*s1-(d1+d2)/2-1),r1=1..nops(L)/2-1), s1=1..nops(L)/2-1)}<>{0} then print(`Something is wrong`): RETURN(FAIL): fi: if {seq( seq(coeff(L[2*r1+1][2*s1+1],n,d1*r1+d2*s1)/c1^r1/c2^s1/(2*r1)!*r1!*2^r1/(2*s1)!*s1!*2^s1, r1=1..nops(L)/2-1),s1=1..nops(L)/2-1)}<>{1} then print(`Not asymptotically normal`): RETURN(FAIL): fi: gue:=[]: for i from 1 to ORD1 do lu:=[seq( [seq(coeff(L[2*r1+1][2*s1+1],n,d1*r1+d2*s1-i)/c1^r1/c2^s1/(2*r1)!*r1!*2^r1/(2*s1)!*s1!*2^s1, r1=0..nops(L)/2-1)],s1=0..nops(L)/2-1)]: gue1:=GuessPOL(lu,[r,s]): if gue1=FAIL then RETURN(FAIL,[gue,guo]): fi: gue1:=expand(subs({r=r+1,s=s+1},gue1)): gue:=[op(gue),gue1]: od: gue:=1+add(gue[i]/n^i,i=1..ORD1): guo:=[]: for i from 0 to ORD1 do lu:=[seq( [seq(coeff(L[2*r1][2*s1],n,d1*r1+d2*s1-(d1+d2)/2-1-i)/c1^r1/c2^s1/ eval((2*r1)!/r1!/2^r1)/eval((2*s1)!/s1!/2^s1), r1=1..nops(L)/2)],s1=1..nops(L)/2)]: guo1:=GuessPOL(lu,[r,s]): if guo1=FAIL then RETURN(FAIL,[gue,guo]): fi: guo:=[op(guo),guo1]: od: guo:=add(guo[i+1]/n^i,i=0..ORD1): for r1 from 1 to nops(L)/2-1 do for s1 from 1 to nops(L)/2-1 do if degree( expand(L[2*r1+1][2*s1+1]-c1^r1*c2^s1* eval((2*r1)!/2^r1/r1!)*eval((2*s1)!/2^s1/s1!)*subs({r=r1,s=s1},gue)* n^(d1*r1+d2*s1)),n) >d1*r1+d2*s1-ORD1-1 then print(`Something is wrong`): RETURN(FAIL): fi: od: od: for r1 from 1 to nops(L)/2-1 do for s1 from 1 to nops(L)/2-1 do if degree( expand(L[2*r1][2*s1]-c1^r1*c2^s1* eval((2*r1)!/2^r1/r1!)*eval((2*s1)!/2^s1/s1!)*subs({r=r1,s=s1},guo)* n^(d1*r1+d2*s1-(d1+d2)/2-1)),n) >d1*r1+d2*s1-(d1+d2)/2-1-ORD1 then print(`kaki,Something is wrong`): RETURN(FAIL): fi: od: od: lu2a:=expand(L[3][1]/n^d1/c1): lu2b:=expand(L[1][3]/n^d2/c2): Gue:=gue/lu2a^r/lu2b^s: Guo:=guo/lu2a^(r-1/2)/lu2b^(s-1/2)*sqrt(c1*c2): Gue:=normal(subs(n=1/x,Gue)): Guo:=normal(subs(n=1/x,Guo)): Gue:=taylor(Gue,x=0,ORD1+1): Guo:=taylor(Guo,x=0,ORD1+1): Gue:=(2*r)!/r!/2^r*(2*s)!/s!/2^s*add(coeff(Gue,x,i1)/n^i1,i1=0..ORD1): Guo:= (2*r)!/r!/2^r*(2*s)!/s!/2^s*add(coeff(Guo,x,i1)/n^i1,i1=0..ORD1)/n: L[1][3],L[3][1],L[2][2], gue,guo,Gue,Guo: end: #AnalyseMoms2V(L,n,ORD1,r,s): Verbose version of #AnalyseMoms2(L,n,ORD1,r,s). For example, try: #read oHISTABRUT20: #AnalyseMoms2V(InvMajMoms10,n,1,r,s); AnalyseMoms2V:=proc(L,n,ORD1,r,s) local gu,mu,c1,c2,d1,d2: gu:=AnalyseMoms2(L,n,ORD1,r,s): if gu=FAIL then RETURN(FAIL): fi: print(`The variance of the first r.v. is`): print(gu[1]): print(`The variance of the second r.v. is`): print(gu[2]): print(`The covariance is`): print(gu[3]): mu:=gu[3]/gu[2]^2: mu:=normal(subs(n=1/x,mu)): mu:=taylor(mu,x=0,ORD1+1): mu:=add(coeff(mu,x,i)/n^i,i=0..ORD1): print(`The asympotic correlation to order`, ORD1, `is : `): print(mu): d1:=degree(gu[1],n): d2:=degree(gu[2],n): c1:=coeff(gu[1],n,d1): c2:=coeff(gu[1],n,d2): print(`The asympotics of the `, [2*r,2*s], `mixed moment is`): print((2*r)!/r!/2^r* (2*s)!/s!/2^s*c1^r*c2^s*n^(d1*r+d2*s)*gu[4]): print(`The asympotics of the `, [2*r-1,2*s-1], `mixed moment is`): print((2*r)!/r!/2^r* (2*s)!/s!/2^s*c1^r*c2^s*n^(d1*r+d2*s-(d1+d2)/2-1)*gu[5]): print(`The asymptotics of the `, [2*r,2*s], ` mixed alpha coefficient is` ): print(gu[6]): print(`The asymptotics of the `, [2*r+1,2*s+1], `mixed alpha coefficient is` ): print(gu[7]): end: