###################################################################### ##BallsInBoxes: Save this file as BallsInBoxes # ## To use it, stay in the # ##same directory, get into Maple (by typing: maple ) # ##and then type: read BallsInBoxes # ##Then follow the instructions given there # ## # ##Written by Doron Zeilberger, Rutgers University , # #zeilberg at math dot rutgers dot edu # ###################################################################### #Created: June 13, 2011 Digits:=30: print(`Created: June 13, 2011 [Herb's Wilf 80th Birthday]`): print(` This is BallsInBoxes `): print(`that investigate the "maximum number of balls in a box"`): print(`when r balls are placed, uniformly at random, in n boxes.`): print(`It accompanies the article `): print(`Balls in Boxes: Variations on a Theme of Warren Ewens`): print(`and Herbert Wilf `): print(`by Shalosh B. Ekhad and Doron Zeilberger`): print(`More accurately: that article accompanies this package!`): print(`Both the article and the package are available from`): print(`http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/bib.html`): print(``): print(`Please report bugs to zeilberg at math dot rutgers dot edu`): print(``): print(`The most current version of this package is in the following url:`): print(`http://www.math.rutgers.edu/~zeilberg/tokhniot/BallsInBoxes .`): print(`------------------------------------`): print(`For a list of the GENERAL procedures type ezraG();, for help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(``): print(`-------------------------------------------`): print(`------------------------------------`): print(`For a list of the Guessing Formulas `): print(` procedures type ezraF();, for help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(``): print(`-------------------------------------------`): print(`------------------------------------`): print(`For a list of the Distribution procedures type ezraDistr();,`): print(` for help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(``): print(`-------------------------------------------`): print(`------------------------------------`): print(`------------------------------------`): print(`For a list of the Sequence procedures type ezraG();, for help with`): print(`a specific procedure, type ezraSeq(procedure_name); .`): print(``): print(`-------------------------------------------`): print(`------------------------------------`): print(`For a list of the MAIN procedures type ezra();, for help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(``): print(`-------------------------------------------`): print(`------------------------------------`): print(`For a list of the Story procedures type ezraStory();, for help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(``): print(`-------------------------------------------`): print(`------------------------------------`): print(`For a list of the supporting procedures type ezra1();, for help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(``): print(`-------------------------------------------`): print(`------------------------------------`): print(`For a list of the checking procedures type ezraCheck(); .`): print(` For help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(``): print(`-------------------------------------------`): print(`------------------------------------`): print(`For a list of the PreComputed procedures type ezraPC();, `): print(` for help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(``): print(`-------------------------------------------`): with(combinat): ezraPC:=proc() if args=NULL then print(` The Pre-Computed procedures are `): print(` AsyabmPC `): else ezra(args): fi: end: ezraGarbage:=proc() if args=NULL then print(` The Garbage procedures, that are abandoned are: `): print(` Smallestm`): print(` EmpForAve, EmpForMom, `): else ezra(args): fi: end: ezraDistr:=proc() if args=NULL then print(` The Probability-Distributions procedures are: `): print(`Pdistr, PdistrEmpir, PdistrEmpirS, PdistrPA `): print(`PdistrC, PdistrEmpirC, PdistrPAC `): print(`PdistrCv, PdistrEmpirCv, PdistrPACv `): else ezra(args): fi: end: ezraG:=proc() if args=NULL then print(` The checking procedures are: `): print(` ArnmG, AsyabmG, AsyabmGPA,`): print(`AsyabmGV,CheckRecabmG, RecabmG, RecabmGV,PrnmG, PrnmGPA`): else ezra(args): fi: end: ezraSeq:=proc() if args=NULL then print(` The Integer Sequences procedures are: `): print(` RecabmSeq, RecabmSeqV `): else ezra(args): fi: end: ezraCheck:=proc() if args=NULL then print(` The checking procedures are: `): print(` CheckRecabm `): else ezra(args): fi: end: ezraF:=proc() if args=NULL then print(` The guessing formulas procedures are: AveFormula, `): print(` NuskhaEmpir1, NuskhaPA,NuskhaPA1 `): else ezra(args): fi: end: ezraStory:=proc() if args=NULL then print(` The story procedures are: `): print(` AsyabmV, AsyDBv, AveFormulaStory, PrnmExactStory`): print(` RecabmV , SeferRecs `): else ezra(args): fi: end: ezra1:=proc() if args=NULL then print(` The supporting procedures are: `): print(` Arnm, , CFAveab, CFAveabEmpir `): print(`CFAveabPA, CFMomab, CFMomabPA`): print(` CFAveabPANor, CFMomabPA, Cnr, Cnr1, DB227, HowGoodPA1, `): print(` HowGoodPA1exact,Pdistr1, PrnmExact, PrnmGexact, RecAm `): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(`The main procedures are: `): print(` AsyDB, Asyabm, AsyabmEmpir, AsyabmPA, `): print(` LargestmPA, Prnm, PrnmEmpir, PrnmExact,`): print(` PrnmPA, PrnmReliable, Recabm, , SmallestmPA `): elif nops([args])=1 and op(1,[args])=Arnm then print(`Arnm(r,n,m): The exact coefficient of x^r in`): print(`(1+x+x^2/2+...+x^m/m!)^n`): print(`For example, try:`): print(`Arnm(100,100,7);`): elif nops([args])=1 and op(1,[args])=ArnmG then print(`ArnmG(r,n,m1,m2): The exact coefficient of x^r in`): print(`(x^m1/m1!+ ... + x^m2/m2!)^n`): print(`For example, try:`): print(`ArnmG(100,100,1,7);`): elif nops([args])=1 and op(1,[args])=Asyabm then print(`Asyabm(a,b,m,M,K,n): An asymptotic formula, to order M`): print(`(using K terms of the sequence) for Prnm(a*n,b*n,m);`): print(`For example, try:`): print(`Asyabm(1,1,4,3,1000,n); `): elif nops([args])=1 and op(1,[args])=AsyabmEmpir then print(`AsyabmEmpir(a,b,m,K,n): The estimated asymptotic formula`): print(`(using K terms of the sequence) for Prnm(a*n,b*n,m);`): print(`For example, try:`): print(`AsyabmEmpir(1,1,4,100,n);`): elif nops([args])=1 and op(1,[args])=AsyabmGPA then print(`AsyabmGPA(a,b,m1,m2,n): The Poisson approximation to`): print(`AsyabmG(a,b,m1,m2,M,K,n) .`): print(`For example, try:`): print(`AsyabmGPA(1,1,1,4,n); `): elif nops([args])=1 and op(1,[args])=AsyabmPA then print(`AsyabmPA(a,b,m,n): The Poisson approximation to`): print(`Asyabm(a,b,m,M,K,n) .`): print(`For example, try:`): print(`AsyabmPA(1,1,4,n); `): elif nops([args])=1 and op(1,[args])=AsyabmG then print(`AsyabmG(a,b,m1,m2,M,K,n): An asymptotic formula, to order M`): print(`(using K terms of the sequence) for PrnmG(a*n,b*n,m1,m2);`): print(`For example, try:`): print(`AsyabmG(2,1,1,4,3,1000,n); `): elif nops([args])=1 and op(1,[args])=AsyabmGV then print(`AsyabmGV(a,b,m1,m2,M,K,n): `): print(`Verbose version of AsyabmG(a,b,m1,m2,M,K,n) .`): print(`For example, try:`): print(`AsyabmGV(2,1,1,4,3,1000,n); `): elif nops([args])=1 and op(1,[args])=AsyabmPC then print(`AsyabmPC(a,b,m,n): The prec-computed asymptotic formula, to order 4`): print(`(using 5000 terms of the sequence) for Prnm(a*n,b*n,m);`): print(`for 1<=a,b<=2 and 2<=m<=7`): print(`For example, try:`): print(`AsyabmPC(1,1,4,n);`): elif nops([args])=1 and op(1,[args])=AsyabmV then print(`AsyabmV(a,b,m,M,K,n): `): print(`Verbose version of Asyabm(a,b,m,M,K,n) .`): print(`For example, try:`): print(`AsyabmV(1,1,4,3,1000,n); `): elif nops([args])=1 and op(1,[args])=AsyDB then print(`AsyDB(Maxa,Maxb,Maxm,M,K,n): creates a data base for the`): print(`asymptotic formulas, in n, of order M, for the probability of no box`): print(`getting more than m boxes when you place a*n balls in `): print(`b*n boxes for a and b relatively prime a<=Maxa, b<=Maxb`): print(`m<=Maxm. K is a large positive integer for estimation purposes.`): print(`The output is a list of pairs [[a,b,m],formula]`): print(`For example, try:`): print(`AsyDB(2,2,4,4,5000,n);`): elif nops([args])=1 and op(1,[args])=AsyDBv then print(`AsyDBv(Maxa,Maxb,Maxm,M,K,n): A verbose version of`): print(`AsyDB(Maxa,Maxb,Maxm,M,K,n).`): print(`For example, try:`): print(`AsyDBv(2,2,4,4,5000,n);`): elif nops([args])=1 and op(1,[args])=AveFormula then print(`AveFormula(a,b,n,d,L,k): A conjectured expression of degree d in`): print(`log(n), for the`): print(`average of the r.v. "maximum number of balls".`): print(`when a*n balls are placed in b*n boxes. Using`): print(`the AsyabmEmpir(a,b,m,K,n) and least-squares`): print(`For example, try:`): print(`AveFormula(2,1,n,4,100,1000,10);`): elif nops([args])=1 and op(1,[args])=AveFormulaStory then print(`AveFormulaStory(A): All the average formulas for the`): print(`random variable "max. number of balls" `): print(`Rn balls n boxes 1<=R<=A and n balls in Rn boxes`): print(`for example try:`): print(`AveFormulaStory(5,n,4,100,1000,10);`): elif nops([args])=1 and op(1,[args])=CheckRecabm then print(`CheckRecabm(a,b,m,K): Checks the recurrence of RecCm(a,b,m,n,N)`): print(`to K terms`): print(`For example, try:`): print(`CheckRecabm(2,1,3,300);`): elif nops([args])=1 and op(1,[args])=CheckRecabmG then print(`CheckRecabmG(a,b,m1,m2,K): Checks the recurrence of RecCm(a,b,m,n,N)`): print(`to K terms`): print(`For example, try:`): print(`CheckRecabmG(2,1,1,3,300);`): elif nops([args])=1 and op(1,[args])=CFAveab then print(`CFAveab(a,b,n,eps,d,L,Sta): `): print(` guesses a polynomial in log(n) of degree d`): print(`that best fits the quantity: "Average Number of the maximum number`): print(`of balls that lend in any box" if a*n balls are placed (uniformly`): print(`at random) in b*n boxes. using L data points`): print(`Sta is the starting place for the estimate.`): print(`For example, try:`): print(`CFAveab(1,1,n,10^(-10),1,4, 1000);`): elif nops([args])=1 and op(1,[args])=CFAveabEmpir then print(`CFAveabEmpir(a,b,n,eps,d,L,Sta,K): Like`): print(`CFAveab(a,b,n,eps,d,L,Sta): but using PdistrEmpir instead of`): print(`Pdistr instead`): print(`For example, try:`): print(`CFAveabEmpir(1,1,n,10^(-10),1,5,1000,50);`): elif nops([args])=1 and op(1,[args])=CFAveabPA then print(`CFAveabPA(a,b,n,eps,d,L,Sta): The Poisson approximation to`): print(`CFAveab(a,b,n,eps,d,L,Sta) `): print(`For example, try:`): print(`CFAveabPA(1,1,n,10^(-10),1,10,1000);`): elif nops([args])=1 and op(1,[args])=CFAveabPANor then print(`CFAveabPANor(a,b,n,eps,d,L,Sta): bn -> n in `): print(`CFAveabPA(a,b,n,eps,d,L,Sta) `): print(`For example, try:`): print(`CFAveabPANor(1,1,n,10^(-10),1,10,1000);`): elif nops([args])=1 and op(1,[args])=CFMomab then print(`CFMomab(a,b,n,eps,d,L,Sta,K): `): print(` guesses a polynomial in log(n) of degree d`): print(`that best fits the quantity: "Average Number of the maximum number`): print(`of balls that lend in any box" if a*n balls are placed (uniformly`): print(`at random) in b*n boxes. using L data points`): print(`as well as the standad-deviation and the higher moments`): print(`up to the K-th . `): print(`Sta is the starting place for the estimate.`): print(`For example, try:`): print(`CFMomab(1,1,n,10^(-10),1,4, 1000,4);`): elif nops([args])=1 and op(1,[args])=CFMomabPA then print(`CFMomabPA(a,b,n,eps,d,L,Sta,K): `): print(`The Poisson approximation to`): print(`CFMomab(a,b,n,eps,d,L,Sta,K): `): print(`For example, try:`): print(`CFMomabPA(1,1,n,10^(-10),1,4, 1000,4);`): elif nops([args])=1 and op(1,[args])=CFMomabPANor then print(`CFMomabPANor(a,b,n,eps,d,L,Sta,K): `): print(`normalized version b*n->n of`): print(`CFMomabPA(a,b,n,eps,d,L,Sta,K): `): print(`For example, try:`): print(`CFMomabPANor(1,1,n,10^(-10),1,4, 1000,4);`): elif nops([args])=1 and op(1,[args])=Cnr then print(`Cnr(h,x,n,r): the coeff. of x^r in the poly. h(x)^n,`): print(`by the naive way.`): print(`For example, try:`): print(`Cnr(1+x+x^2,x,10,5);`): print(``): elif nops([args])=1 and op(1,[args])=Cnr1 then print(`Cnr1(f,x,n,r): an appx. for the coeff. of x^r in the poly. f(x)^n,`): print(`by numerical integration`): print(`For example, try:`): print(`Cnr1(1+x+x^2,x,10,5);`): print(``): elif nops([args])=1 and op(1,[args])=DB227 then print(`DB227(n): The data base of all asymptotic formulas for`): print(`Asyabm(a,b,m,4,5000,n) for 1<=a,b<=2, 2<=m<=7`): print(`For example, try:`): print(`DB227(n);`): elif nops([args])=1 and op(1,[args])=EmpForAve then print(`EmpForAve(n,R,eps,d,L,Sta,d1,M1,M2): An empirically-derived formula`): print(`as a polynomial of degree d in log(n) and of degree d1`): print(`in R,`): print(`using Least Squares and the Poisson Approximation`): print(`for the quantity:`): print(`The average of the random variable "Max. Number of Balls that lended`): print(`in an individual box" upon placing, uniformaly, n*R balls in n boxes`): print(`the data set is from CFAveabPANor(a,b,n,eps,d,L,Sta) with a between`): print(`1 and M1 and b between 1 and M2`): print(`For example, try: EmpForAve(n,R,10^(-10),1,10,2000,1,3,3);`): elif nops([args])=1 and op(1,[args])=EmpForMom then print(`EmpForMom(n,R,eps,L,Sta,M1,M2,d,d1): An empirically-derived formula`): print(`as a polynomial of degree d in log(n) and of degree d1`): print(`in R,`): print(`using Least Squares and the Poisson Approximation`): print(`for the quantity:`): print(`The average of the random variable "Max. Number of Balls that lended`): print(`in an individual box" upon placing, uniformaly, n*R balls in n boxes`): print(`followed by the standard-deviation and the (normalized) higher`): print(`moments (alpha coefficients), up to the K-th.`): print(`The data set is from CFAveabPANor(a,b,n,eps,d,L,Sta) with a between`): print(`1 and M1 and b between 1 and M2`): print(`For example, try: `): print(`For example, try: EmpForMom(n,R,10^(-10),10,1000,4,3,3,1,1);`): elif nops([args])=1 and op(1,[args])=HowGoodPA1 then print(`HowGoodPA1(R0,N0,Incr,M0,m,eps): explores the ratio of`): print(`PrnmPA(r,n,m) to Prnm(r,n,m) for r in [R0, R0+Incr*M0]`): print(`and n in [N0,N0+Incr*M0] for a specfic m provided`): print(`Prnm(r,n,m) is larger than eps. `): print(`Returns the one the ratio furthest away from 1 together with the`): print(`(n,r) that yield it`): print(`For example, try:`): print(`HowGoodPA1(1000,1000,1000,3,5,1/100);`): elif nops([args])=1 and op(1,[args])=HowGoodPA1exact then print(`HowGoodPA1exact(R0,N0,Incr,M0,m,eps): Like`): print(`HowGoodPA1(R0,N0,Incr,M0,m,eps) but using EXACT`): print(`rational arithmetic, rather than floating points.`): print(`In other words using`): print(`PrnmExact(r,n,m) instead of Prnm(r,n,m)`): print(`For example, try:`): print(`HowGoodPA1exact(1000,1000,1000,3,5,1/100);`): elif nops([args])=1 and op(1,[args])=Largestm then print(`Largestm(r,n,conf): The largest m such that the prob.`): print(`that all boxes would have >=m balls in them with prob.>=conf`): print(`For example, try:`): print(`Largestm(1000,50,.99);`): elif nops([args])=1 and op(1,[args])=LargestmPA then print(`LargestmPA(r,n,conf): The largest m such that the prob.`): print(`that all boxes would have >=m balls in them with prob.>=conf`): print(`using the Poisson Approximation`): print(`For example, try:`): print(`LargestmPA(1000,50,.99);`): elif nops([args])=1 and op(1,[args])=NuskhaEmpir1 then print(`NuskhaEmpir1(R,n,K,d,K1,CU): `): print(`A guessed-formula, as a polynomial in log(n)`): print(`of degree d, for the average, standard-deviation and first`): print(`centralized-normalized K moments of the prob. distr.`): print(` for the random variable`): print(`"maximum number of balls amongst all boxes" when n*R balls are`): print(`thrown, uniformly at random, in n boxes, for R an integer.`): print(`Using the Empirical Asympototic approach as`): print(`given by PdistrEmpir `): print(`For example, try:`): print(`NuskhaEmpir1(1,n,5,3,100,20);`): elif nops([args])=1 and op(1,[args])=NuskhaPA then print(`NuskhaPA(R,n,K,d,d1): A guessed-formula, as a polynomial in log(n)`): print(`of degree d, and in R of degree d1`): print(`for the average, standard-deviation and first`): print(`centralized-normalized K moments of the prob. distr. `): print(`for the random variable`): print(`"maximum number of balls amongst all boxes" when n*R balls are`): print(`thrown, uniformly at random, in n boxes, for R an integer.`): print(`Using the Poisson Approximation`): print(`For example, try:`): print(`NuskhaPA(R,n,5,3,1);`): elif nops([args])=1 and op(1,[args])=NuskhaPA1 then print(`NuskhaPA1(R,n,K,d): A guessed-formula, as a polynomial in log(n)`): print(`of degree d, for the average, standard-deviation and first`): print(`centralized-normalized K moments of the prob. distr.`): print(` for the random variable`): print(`"maximum number of balls amongst all boxes" when n*R balls are`): print(`thrown, uniformly at random, in n boxes, for R an integer.`): print(`Using the Poisson Approximation`): print(`For example, try:`): print(`NuskhaPA1(1,n,5,3);`): elif nops([args])=1 and op(1,[args])=Pdistr then print(`Pdistr(r,n,eps,K): The list such that`): print(`its i-th entry is the probability that`): print(`in (uniformly) placing r balls in n boxes`): print(`the max. number of balls in any box is i`): print(`The length of the list is such that the prob. of `): print(`having a box with than its length is<=eps .`): print(`It also returns the average, the variance, and the`): print(`alpha coefficients up to order K`): print(`For example, try:`): print(`Pdistr(14400,9000,10^(-10),10);`): elif nops([args])=1 and op(1,[args])=PdistrC then print(`PdistrC(r,n,eps,K): The central part where`): print(`most (1-eps) of`): print(`the probability distribution Pdistr(r,n,eps,K) is located`): print(`followed by the average, variance, and normalized moments`): print(`followed by the average, variance, skewness, kurtosis and`): print(`higher alpha coefficients up to K of them`): print(`and centralized-normalized`): print(`distribution. `): print(`For example, try:`): print(`PdistrC(14400,9000,10^(-3),4);`): elif nops([args])=1 and op(1,[args])=PdistrCv then print(`PdistrCv(r,n,eps,K): Verbose version of`): print(` PdistrC(r,n,eps,K): `): print(`For example, try:`): print(`PdistrCv(300,300,10^(-3),4);`): elif nops([args])=1 and op(1,[args])=PdistrEmpir then print(`PdistrEmpir(r,n,eps,K,K1,CU): The Emprical approximation to`): print(`Pdistr(r,n,eps,K)`): print(`using PrnmEmpir(r,n,m,K1,CU) instead of Prnm(r,n,m)`): print(`For example, try:`): print(`PdistrEmpir(14400,9000,10^(-10),10,100,30);`): elif nops([args])=1 and op(1,[args])=PdistrEmpirC then print(`PdistrEmpirC(r,n,eps,K,K1): The central part where`): print(`most (1-eps) of`): print(`the probability distribution PdistrEmpir(r,n,eps,K) is located`): print(`i.e. using the Empirical Approximation`): print(`followed by the average, variance, skewness, kurtosis etc.`): print(`up to K alpha coefficinets`): print(`followed by the centralized-normalized`): print(`distribution. `): print(`For example, try:`): print(`PdistrEmpirC(14400,9000,10^(-3),8,100);`): elif nops([args])=1 and op(1,[args])=PdistrEmpirCv then print(`PdistrEmpirCv(r,n,eps,K,K1): verbose version of`): print(`PdistrEmpirC(r,n,eps,K,K1)`): print(`For example, try:`): print(`PdistrEmpirCv(300,300,10^(-3),4,100);`): elif nops([args])=1 and op(1,[args])=PdistrEmpirS then print(`PdistrEmpirS(a,b,n,M1,K1,M2): The first M1 members of the`): print(`as (asymptotic) formulas in n of the prob. distribution`): print(`Max number of Balls, using AsyabmEmpir(a,b,m,K1,n)`): print(`followed by the average, s.d. and the normalized moments`): print(`through the M2`): print(`For example, try:`): print(`PdistrEmpirS(2,1,n,10,100,10);`): elif nops([args])=1 and op(1,[args])=PdistrPA then print(`PdistrPA(r,n,eps,K): The Poisson approximation to`): print(`Pdistr(r,n,eps,K)`): print(`For example, try:`): print(`PdistrPA(14400,9000,10^(-10),10);`): elif nops([args])=1 and op(1,[args])=PdistrPAC then print(`PdistrPAC(r,n,eps,K): The central part where`): print(`most (1-eps) of`): print(`the probability distribution PdistrPA(r,n,eps,K) is located`): print(`i.e. using the Poisson Approximation`): print(`followed by the average, variance, and normalized moments`): print(`followed by the average, variance, skewness, kurtosis and`): print(`higher alpha coefficients up to K of them`): print(`and centralized-normalized`): print(`distribution. `): print(`For example, try:`): print(`PdistrPAC(14400,9000,10^(-3),6);`): elif nops([args])=1 and op(1,[args])=PdistrPACv then print(`PdistrPACv(r,n,eps,K): Verbose version of`): print(`PdistrPAC(r,n,eps,K)`): print(`For example, try:`): print(`PdistrPACv(300,300,10^(-3),4);`): elif nops([args])=1 and op(1,[args])=Pdistr1 then print(`Pdistr1(r,n,L): The list of length L such that`): print(`its i-th entry is the probability that`): print(`in (uniformly) placing r balls in n boxes`): print(`the max. number of balls in any box is i`): print(`For example, try:`): print(`Pdistr1(14400,9000,7);`): elif nops([args])=1 and op(1,[args])=Prnm then print(`Prnm(r,n,m): The probability that if you put r balls in n boxes`): print(`each box received at most m balls.`): print(`WARNING: It uses a fixed floating-point arithmetic with `, Digits): print(`digits. It is much safer to use PrnmReliable(r,n,m,k) .`): print(`For example, try:`): print(`Prnm(14400,9000,7);`): elif nops([args])=1 and op(1,[args])=PrnmEmpir then print(`PrnmEmpir(r,n,m,K,CU): Like PrnmEmpir(r,n,m)`): print(`but using extrapolation of the empirically`): print(`derived Asyabm(a,b,m,K): CU is the cut-off `): print(`of whether to use r/n=a/b or use an approximation.`): print(`For example, try:`): print(`PrnmEmpir(14400,9000,7,100,30);`): elif nops([args])=1 and op(1,[args])=PrnmExact then print(`PrnmExact(r,n,m): The Exact (as a rational number)`): print(`probability that if you put r balls in n boxes`): print(`each box received at most m balls.`): print(`For example, try:`): print(`PrnmExaxt(14400,9000,7);`): elif nops([args])=1 and op(1,[args])=PrnmExactStory then print(`PrnmExactStory(r,n,m1,m2,K): The story of`): print(`the exact evaluation of Prnm(r,n,m) with timings`): print(`and the number of digits. For example, try:`): print(`PrnmExactStory(14400,9000,6,12,30);`): elif nops([args])=1 and op(1,[args])=PrnmReliable then print(`PrnmReliable(r,n,m,eps): computes Prnm(r,n,m)`): print(`with floating-point reliable to k digits`): print(`For example, try:`): print(`PrnmReliable(14400,9000,8,20);`): elif nops([args])=1 and op(1,[args])=PrnmG then print(`PrnmG(r,n,m1,m2): The probability that if you put r balls in n boxes`): print(`each box received between m1 and m2 balls (inclusive).`): print(`For example, try:`): print(`PrnmG(14400,9000,0,7);`): elif nops([args])=1 and op(1,[args])=PrnmGexaxt then print(`PrnmGexaxt(r,n,m1,m2): The EXACT probability (as a rational number)`): print(`that if you put r balls in n boxes`): print(`each box received between m1 and m2 balls (inclusive).`): print(`For example, try:`): print(`PrnmGexaxt(14400,9000,0,7);`): elif nops([args])=1 and op(1,[args])=PrnmGPA then print(`PrnmGPA(r,n,m1,m2): The Poisson approximation for PrnmG(r,n,m1,m2)`): print(`For example, try:`): print(`PrnmGPA(14400,9000,0,7);`): elif nops([args])=1 and op(1,[args])=PrnmPA then print(`PrnmPA(r,n,m): The Poisson approximation for Prnm(r,n,m)`): print(`For example, try:`): print(`PrnmPA(14400,9000,7);`): elif nops([args])=1 and op(1,[args])=RecAm then print(`RecAm(a,m,n,N): The linear recurrence (with poly. coeffs.) `): print(`satistfied by`): print(`the coeff. of x^(a*n) in (1+x+x^2/2!+...+x^m/m!)^n`): print(`where a in an integer (a*n balls in n boxes).`): print(`For example, try:`): print(`RecAm(2,3,n,N);`): elif nops([args])=1 and op(1,[args])=Recabm then print(`Recabm(a,b,m,n,N): The linear recurrence (with poly. coeffs.) `): print(`satistfied by`): print(`the coeffs. of x^(a*n) in (1+x+x^2/2!+...+x^m/m!)^(b*n)`): print(`where a and b are integers (a*n balls in b*n boxes).`): print(`For example, try:`): print(`Recabm(2,1,3,n,N);`): elif nops([args])=1 and op(1,[args])=RecabmSeq then print(`RecabmSeq(a,b,m,n,N,L): The linear recurrence (with poly. coeffs.) `): print(`satistfied by`): print(`the sequence that counts the number of ways of placing`): print(`a*n balls in b*n boxes such that no box receives more than m`): print(`m balls.`): print(`also the initial values and the first L terms`): print(`For example, try:`): print(`RecabmSeq(2,1,3,n,N,30);`): elif nops([args])=1 and op(1,[args])=RecabmSeqV then print(`RecabmSeqV(a,b,m,n,N,L): verbose form of`): print(`RecabmSeq(a,b,m,n,N,L)`): print(`For example, try:`): print(`RecabmSeqV(2,1,3,n,N,30);`): elif nops([args])=1 and op(1,[args])=RecabmG then print(`RecabmG(a,b,m1,m2,n,N): The linear recurrence (with poly. coeffs.) `): print(`satistfied by`): print(`the coeffs. of x^(a*n) in (x^m1/m1!+ ...+ x^m2/m2!)^(b*n)`): print(`where a iand b are integers (a*n balls in b*n boxes).`): print(`For example, try:`): print(`RecabmG(2,1,1,3,n,N);`): elif nops([args])=1 and op(1,[args])=RecabmGV then print(`RecabmGV(a,b,m1,m2,n,N): verbose version of RecabmG(a,b,m1,m2,n,N)`): print(`For example, try:`): print(`RecabmGV(2,1,1,3,n,N);`): elif nops([args])=1 and op(1,[args])=RecabmV then print(`RecabmV(a,b,m,n,N): verbose version of RecabmV(a,b,m,n,N)`): print(`For example, try:`): print(`RecabmV(2,1,3,n,N);`): elif nops([args])=1 and op(1,[args])=SeferRecs then print(`SeferRecs(Maxa,Maxb,Maxm,L,n,N): Inputs positive integers`): print(`Maxa, Maxb and Maxm, as well as symbols n and N (shift in n)`): print(`and outputs a web-book with the recurrences, and the first`): print(`L terms for the sequences enumerating the number of ways`): print(`of placing a*n balls in b*n boxes such that no box has more`): print(`than m balls, for a and b relatively prime `): print(`a<=Maxa, b<=Maxb, m<=Maxm. `): print(`For example, try:`): print(`SeferRecs(2,2,4,30,n,N);`): elif nops([args])=1 and op(1,[args])=Smallestm then print(`Smallestm(r,n,conf): The smallest m such that the prob.`): print(`that all boxes would have <=m balls in them with prob.>=conf`): print(`For example, try:`): print(`Smallestm(14400,9000,.99);`): elif nops([args])=1 and op(1,[args])=SmallestmPA then print(`SmallestmPA(r,n,conf): The smallest m such that the prob.`): print(`that all boxes would have <=m balls in them with prob.>=conf`): print(`USING THE (excellent!) Poisson Approximation!`): print(`For example, try:`): print(`SmallestmPA(14400,9000,.99);`): else print(`There is no ezra for`,args): fi: end: #Cnr(h,x,n,r): the coeff. of x^r in the poly. h(x)^n, #by the naive way #For example, try: #Cnr(1+x+x^2,x,10,5); Cnr:=proc(h,x,n,r): coeff(expand(h^n),x,r): end: #Cnr1(h,x,n,r): the coeff. of x^r in the poly. h(x)^n, #by numerical integration #For example, #Cnr1(1+x+x^2,x,10,5); Cnr1:=proc(h,x,n,r) local t: evalf(Int(subs(x=exp(I*t),h)^n/exp(r*I*t)/2/Pi,t=0..2*Pi)): end: #CnrF(f,x,n,r): the coeff. of x^r in the poly. h(x)^n, #by the fast way (that goes back to Euler) #For example, try: #CnrF(1+x+x^2,x,10,5); CnrF:=proc(f,x,n,r) local gu,a,L,m,hs,s,kat: kat:=ldegree(f,x): if r0 then RETURN(CnrF(normal(f/x^(kat)),x,n,r-n*kat)): fi: m:=degree(f,x): for L from 0 to m do a[L]:=coeff(f,x,L): od: gu:=[a[0]^n]: for s from 1 to m do hs:=add(((n+1)*L-s)*gu[nops(gu)-L+1]*a[L],L=1..s)/s/a[0]: gu:=[op(gu),hs]: od: if r<=m then RETURN(gu[r+1]): fi: for s from m+1 to r do hs:=add(((n+1)*L-s)*gu[nops(gu)-L+1]*a[L],L=1..m)/s/a[0]: gu:=[op(2..nops(gu),gu),hs]: od: gu[nops(gu)]: end: #PrnmExaxt(r,n,m): The probability that if you put r balls in n boxes #each box received at most m balls: #For example, try: #PrnmExact(14400,9000,7); PrnmExact:=proc(r,n,m) local f,i,x: f:=add(x^i/i!,i=0..m): r!/n^r*CnrF(f,x,n,r): end: #Arnm(r,n,m): The exact coefficient of x^r in #(1+x+x^2/2+...+x^m/m!)^n #For example, try: #Arnm(100,100,7); Arnm:=proc(r,n,m) local f,i,x: f:=add(x^i/i!,i=0..m): CnrF(f,x,n,r): end: #ArnmG(r,n,m1,m2): The exact coefficient of x^r in #(x^m1/m1!+ ...+x^m2/m2!)^n #For example, try: #ArnmG(100,100,0,7); ArnmG:=proc(r,n,m1,m2) local f,i,x: f:=add(x^i/i!,i=m1..m2): CnrF(f,x,n,r): end: #Prnm(r,n,m): The probability (in floating point) #that if you put r balls in n boxes #each box received at most m balls: #For example, try: #Prnm(14400,9000,7); Prnm:=proc(r,n,m) local f,i,x: if r>n*m then RETURN(0): fi: f:=evalf(add(x^i/i!,i=0..m)): r!/n^r*CnrF(f,x,n,r): end: #Pdistr1(r,n,L): The list of length L such that #its i-th entry is the probability that #in (uniformly) placing r balls in n boxes #the max. number of balls in any box is i #For example, try: #Pdistr1(14400,9000,9); Pdistr1:=proc(r,n,L) local gu,i,mu1,mu2: mu1:=Prnm(r,n,1): gu:=[mu1]: for i from 2 to L do mu2:=Prnm(r,n,i): gu:=[op(gu),mu2-mu1]: mu1:=mu2: od: gu: end: #Pdistr(r,n,eps,K): The list such that #its i-th entry is the probability that #in (uniformly) placing r balls in n boxes #the max. number of balls in any box is i #the list is such that the probability of getting #more than its length is <=eps. It also returns #the average, standard deviation, and normalized higher moments #to order K. #For example, try: #Pdistr(14400,9000,1/10^5,4); Pdistr:=proc(r,n,eps,K) local gu,i,mu1,mu2,av,m2,hakol,j: mu1:=Prnm(r,n,1): gu:=[mu1]: for i from 2 while convert(gu,`+`)<1-eps do mu2:=Prnm(r,n,i): gu:=[op(gu),mu2-mu1]: mu1:=mu2: od: hakol:=convert(gu,`+`): av:=add(i*gu[i],i=1..nops(gu))/hakol: m2:=add(gu[i]*(i-av)^2,i=1..nops(gu))/hakol: gu, [av,sqrt(m2),seq(add(gu[i]*(i-av)^j,i=1..nops(gu))/hakol/m2^(j/2),j=3..K)]: end: #PdistrPA(r,n,eps,K): Pdistr(r,n,eps,K) using the Poisson Approximation #For example, try: #PdistrPA(14400,9000,1/10^5,4); PdistrPA:=proc(r,n,eps,K) local gu,i,mu1,mu2,av,m2,hakol,j: mu1:=PrnmPA(r,n,1): gu:=[mu1]: for i from 2 while convert(gu,`+`)<1-eps do mu2:=PrnmPA(r,n,i): gu:=[op(gu),mu2-mu1]: mu1:=mu2: od: hakol:=convert(gu,`+`): av:=add(i*gu[i],i=1..nops(gu))/hakol: m2:=add(gu[i]*(i-av)^2,i=1..nops(gu))/hakol: gu, [av,sqrt(m2),seq(add(gu[i]*(i-av)^j,i=1..nops(gu))/hakol/m2^(j/2),j=3..K)]: end: #PdistrEmpir(r,n,eps,K,K1,CU): Pdistr(r,n,eps,K) using the #Empirical approximation PrnmEmpir(r,n,i,K1,CU): #For example, try: #PdistrEmpir(14400,9000,1/10^5,10,100,20); PdistrEmpir:=proc(r,n,eps,K,K1,CU) local gu,i,mu1,mu2,av,m2,hakol,j: mu1:=PrnmEmpir(r,n,1,K1,CU): gu:=[mu1]: for i from 2 while convert(gu,`+`)<1-eps do mu2:=PrnmEmpir(r,n,i,K1,CU): gu:=[op(gu),mu2-mu1]: mu1:=mu2: od: hakol:=convert(gu,`+`): av:=add(i*gu[i],i=1..nops(gu))/hakol: m2:=add(gu[i]*(i-av)^2,i=1..nops(gu))/hakol: gu, [av,sqrt(m2),seq(add(gu[i]*(i-av)^j,i=1..nops(gu))/hakol/m2^(j/2),j=3..K)]: end: ##################From EKHAD############ pashetd:=proc(p,N) local gu,p1,mekh,i: p1:=normal(p): mekh:=denom(p1): p1:=numer(p1): p1:=expand(p1): gu:=0: for i from 0 to degree(p1,N) do gu:=gu+factor(coeff(p1,N,i))*N^i: od: RETURN(gu,mekh): end: goremd:=proc(N,ORDER) local i, gu: gu:=0: for i from 0 to ORDER do gu:=gu+(bpc[i])*N^i: od: RETURN(gu): end: duisd:= proc(A,ORDER,y,n,N) local gorem, conj, yakhas,lu1a,LU1,P,Q,R,S,j1,g,Q1,Q2,l1,l2,mumu, mekb1,fu,meka1,k1,gugu,eq,ia1,va1,va2,degg,shad,KAMA,i1: KAMA:=10: gorem:=goremd(N,ORDER): conj:=gorem: yakhas:=0: for i1 from 0 to degree(conj,N) do yakhas:=yakhas+coeff(conj,N,i1)*simplify(subs(n=n+i1,A)/A): od: lu1a:=yakhas: LU1:=numer(yakhas): yakhas:=1/denom(yakhas): yakhas:=normal(diff(yakhas,y)/yakhas+diff(A,y)/A): P:=LU1: Q:=numer(yakhas): R:=denom(yakhas): j1:=0: while j1 <= KAMA do g:=gcd(R,Q-j1*diff(R,y)): if g <> 1 then Q2:=(Q-j1*diff(R,y))/g: R:=normal(R/g): Q:=normal(Q2+j1*diff(R,y)): P:=P*g^j1: j1:=-1: fi: j1:=j1+1: od: P:=expand(P): R:=expand(R): Q:=expand(Q): Q1:=Q+diff(R,y): Q1:=expand(Q1): l1:=degree(R,y): l2:=degree(Q1,y): meka1:=coeff(R,y,l1): mekb1:=coeff(Q1,y,l2): if l1-1 <>l2 then k1:=degree(P,y)-max(l1-1,l2): else mumu:= -mekb1/meka1: if type(mumu,integer) and mumu > 0 then k1:=max(mumu, degree(P,y)-l2): else k1:=degree(P,y)-l2: fi: fi: fu:=0: if k1 < 0 then RETURN(0): fi: for ia1 from 0 to k1 do fu:=fu+apc[ia1]*y^ia1: od: gugu:=Q1*fu+R*diff(fu,y)-P: gugu:=expand(gugu): degg:=degree(gugu,y): for ia1 from 0 to degg do eq[ia1+1]:=coeff(gugu,y,ia1)=0: od: for ia1 from 0 to k1 do va1[ia1+1]:=apc[ia1]: od: for ia1 from 0 to ORDER do va2[ia1+1]:=bpc[ia1]: od: eq:=convert(eq,set): va1:=convert(va1,set): va2:=convert(va2,set): va1:=va1 union va2 : va1:=solve(eq,va1): fu:=subs(va1,fu): gorem:=subs(va1,gorem): if fu=0 and gorem=0 then RETURN(0): fi: for ia1 from 0 to k1 do gorem:=subs(apc[ia1]=1,gorem): od: for ia1 from 0 to ORDER do gorem:=subs(bpc[ia1]=1,gorem): od: fu:=normal(fu): shad:=pashetd(gorem,N): S:=lu1a*R*fu*shad[2]/P: S:=subs(va1,S): S:=normal(S): S:=factor(S): for ia1 from 0 to k1 do S:=subs(apc[ia1]=1,S): od: for ia1 from 0 to ORDER do S:=subs(bpc[ia1]=1,S): od: RETURN(shad[1],S): end: AZd:= proc(A,y,n,N) local ORDER,gu,KAMA: KAMA:=20: for ORDER from 0 to KAMA do gu:=duisd(A,ORDER,y,n,N): if gu<>0 then RETURN(gu): fi: od: 0: end: ##################End From EKHAD############ #RecAm(a,m,n,N): The linear recurrence (with poly. coeffs.) #satistfied by #the coeffs. of x^(a*n) in (1+x+x^2/2!+...+x^m/m!)^n #where a in an integer (a*n balls in n boxes). #For example, try: #RecAm(2,3,n,N): RecAm:=proc(a,m,n,N) local i,x: AZd(add(x^i/i!,i=0..m)^n/x^(a*n+1),x,n,N)[1]: end: #Recabm(a,b,m,n,N): The linear recurrence (with poly. coeffs.) #satistfied by #the coeffs. of x^(a*n) in (1+x+x^2/2!+...+x^m/m!)^(b*n) #where a and b are pos. integers (a*n balls in b*n boxes). #For example, try: #Recabm(2,1,3,n,N): Recabm:=proc(a,b,m,n,N) local i,x: option remember: AZd(add(x^i/i!,i=0..m)^(n*b)/x^(n*a+1),x,n,N)[1]: end: #RecabmG(a,b,m1,m2,n,N): The linear recurrence (with poly. coeffs.) #satistfied by #the coeffs. of x^(a*n) in (x^m1/m1!+...+x^m2/m2!)^(b*n) #where a and b are pos. integers (b*n balls in a*n boxes). #For example, try: #RecabmG(2,1,0,3,n,N): RecabmG:=proc(a,b,m1,m2,n,N) local i,x: option remember: AZd(add(x^i/i!,i=m1..m2)^(n*b)/x^(n*a+1),x,n,N)[1]: end: #PrnmPA(r,n,m): Poisson approximation for the #The probability (in floating point) #that if you put r balls in n boxes #each box received at most m balls: #For example, try: #PrnmPA(14400,9000,7); PrnmPA:=proc(r,n,m) local a,i: a:=evalf(r/n): (exp(-a)*add(a^i/i!,i=0..m))^n: end: #CheckRecabm(a,b,m,K): Checks the recurrence of Recabm(a,b,m,n,N) #to K terms #For example, try: #CheckRecabm(2,1,3,300); CheckRecabm:=proc(a,b,m,K) local ope,n,N,gu,i,n1,ord1: ope:=Recabm(a,b,m,n,N): ord1:=degree(ope,N): gu:=[seq(Arnm(a*i,b*i,m),i=1..K)]: evalb( {seq(add(subs(n=n1,coeff(ope,N,i))*gu[n1+i],i=0..ord1),n1=1..K-ord1)}={0} ): end: #CheckRecabmG(a,b,m1,m2,K): Checks the recurrence of Recabm(a,b,m,n,N) #to K terms #For example, try: #CheckRecabmG(2,1,1,3,300); CheckRecabmG:=proc(a,b,m1,m2,K) local ope,n,N,gu,i,n1,ord1: ope:=RecabmG(a,b,m1,m2,n,N): ord1:=degree(ope,N): gu:=[seq(ArnmG(a*i,b*i,m1,m2),i=1..K)]: print(`The sequence is`,gu): print(`the recurrence is`, ope): evalb( {seq(add(subs(n=n1,coeff(ope,N,i))*gu[n1+i],i=0..ord1),n1=1..K-ord1)}={0} ): end: #begin AsyRec with(SolveTools): with(numtheory): #MonoN(ope,n,N,k): Given a linear recurrence operator with #poly coeffs. ope(n,N), and a pos. integer k, outputs the expression of #N^k as a linear combination of 1, N, ..., N^(ORDER-1) #For example, try #MonoN(N-n-1,n,N,3); MonoN:=proc(ope,n,N,k) local ORDER,coe0,i,lu1,lu2: ORDER:=degree(ope,N): if kL 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: #HafelOper(ope,n,N,L): applies the operator ope(n,N) to #the sequence L, for example, try: #HafelOper(N-(n+1),n,N,[1,2,6,24,120]); HafelOper:=proc(ope,n,N,L) local i,n1: [seq(add(subs(n=n1,coeff(ope,N,i))*L[n1+i],i=0..degree(ope,N)), n1=1..nops(L)-degree(ope,N))]: end: #TestKsect(ope,n,N,K,r,M,Ini): tests procedure Ksect #with initial conditions Ini up to M terms TestKsect:=proc(ope,n,N,K,r,M,Ini) local gu,Ope,n1: gu:=SeqFromRec(ope,n,N,Ini,M*K+r): Ope:=Ksect(ope,n,N,K,r): gu:=[seq(gu[K*n1+r],n1=1..M-1)]: evalb({op(HafelOper(Ope,n,N,gu))}={0}): end: rf:=proc(a,n) local i:mul(a+i,i=0..n-1):end: #CODV1(ope,n,N,k): Given a linear recurrence operator ope(n,N) #annihilating a[n], say, outputs the operator #annihilating b[n]:=a[n]/n!^k #For example, try: CODV1(N-(n+1),n,N,1): CODV1:=proc(ope,n,N,k) local i: add(coeff(ope,N,i)*(rf(n+1,i))^k*N^i,i=0..degree(ope,N)): end: #CODV(ope,n,N,k,K): Given a linear recurrence operator ope(n,N) #annihilating a[n], say, outputs the operator #annihilating b[n]:=a[n*K]/n!^k #For example, try: CODV(N-(n+1),n,N,1,1): CODV:=proc(ope,n,N,k,K) local Ope: Ope:=Ksect(ope,n,N,K,0): CODV1(Ope,n,N,k): end: #FindKk(ope,n,N): Given a linear recurrence operator with polynomial #coefficients, ope(n,N), finds the integers K and k such #that if a(n) is a solution of ope(n,N)a(n)=0, then #b(n):=a(K*n)/n!^k is annihilated by a standard operator #the output is the pair [k,K] and the transformed operator #For exanple, try: FindKk(N^2-n,n,N); FindKk:=proc(ope,n,N) local ope1,Ld,deg,sp,yakhas,K,k,pu: pu:=Aope1(ope,n,N): if {solve(pu,N)}<>{} and {solve(pu,N)}<>{0} then RETURN([1,0],ope): fi: ope1:=expand(ope): Ld:=ldegree(ope,N): ope1:=expand(subs(n=n-Ld,ope1)/N^Ld ): deg:=degree(ope,N): sp:=degree(coeff(ope,N,0),n)-degree(coeff(ope,N,deg),n): yakhas:=sp/deg: k:=numer(yakhas): K:=denom(yakhas): [K,k],CODV(ope,n,N,k,K): end: #Aope1(ope,n,N): the asymptotics of the difference operator with poly. coeffs. ope(n,N) Aope1:=proc(ope,n,N) local gu: gu:=expand(numer(normal(ope))): gu:=coeff(gu,n,degree(gu,n)): gu:=expand(gu/coeff(gu,N,degree(gu,N))): factor(gu): end: C3:=proc(a,b,c) option remember: if [a,b,c]=[0,0,0] then RETURN(1): fi: if a<0 or b<0 or c<0 then RETURN(0): fi: C3(a-1,b,c)+C3(a,b-1,c)+C3(a,b,c-1)-4*C3(a-1,b-1,c-1): end: #Aluf1(pit): is it positive dominant? Aluf1:=proc(pit) local aluf,shi,i,makom: shi:=evalf(abs(pit[1])): aluf:={evalf(evalc(pit[1]))}: makom:={1}: for i from 2 to nops(pit) do if evalf(abs(pit[i]))>evalf(shi) then shi:=abs(evalf(evalc(pit[i]))): aluf:={pit[i]}: makom:={i}: elif evalf(abs(pit[i]))=shi then aluf:=aluf union {evalf(evalc(pit[i]))}: makom:=makom union {i}: fi: od: aluf,shi,makom: end: #OneStepAS1(ope1,n,N,alpha,f,S1): Given the partial asymptotic #expansion of solutions of ope1(n,N) with exponent alpha #extends it by one term OneStepAS1:=proc(ope1,n,N,alpha,f,S1) local x,f1,L,F,A,mu,i,A1: f1:=subs(n=1/x,f): L:=degree(f1,x): F:=f1+A*x^(L+S1): mu:=add( subs(n=1/x,coeff(ope1,N,i))*(1+i*x)^alpha*subs(x=x/(1+i*x),F) ,i=ldegree(ope1,N)..degree(ope1,N)): mu:=normal(mu): mu:=taylor(mu,x=0,L+11): mu:=simplify(mu): for i from L+1 to L+9 while coeff(mu,x,i)=0 do od: if i=L+10 then RETURN(FAIL): fi: mu:=coeff(mu,x,i): A1:=subs(Linear({mu},{A}),A): if A1=0 then RETURN(FAIL): fi: subs({x=1/n,A=A1},F): end: #OneStepA(ope1,n,N,alpha,f): Given the partial asymptotic #expansion of solutions of ope1(n,N) with exponent alpha #extends it by one term OneStepA:=proc(ope1,n,N,alpha,f) local S1,pu: for S1 from 1 to 5 while OneStepAS1(ope1,n,N,alpha,f,S1)=FAIL do od: pu:=OneStepAS1(ope1,n,N,alpha,f,S1): if pu=FAIL then RETURN(f): else RETURN(pu): 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: #Nor(ope,n,N): the Normalizer of the linear recurrence #operator with polynomial coefficients ope(n,N) #followed by its exponential growth constant #For example, try: #Nor(N^2-3*N+2,n,N); Nor:=proc(ope,n,N) local gu,pit,lu,makom,ope1: gu:=Aope1(ope,n,N): if degree(gu,N)=0 then print(`Not of exponential growth, case not implemented`): RETURN(FAIL): fi: if not type(expand(normal(ope/gu)),`+`) then pit:=[solve(gu,N)]: makom:=Aluf1(pit)[3]: lu:=pit[makom[1]]: ope1:=Yafe(subs(N=lu*N,ope),N)[2]: RETURN(ope1,lu): fi: pit:=[solve(gu,N)]: makom:=Aluf1(pit)[3]: if nops(makom)<>1 then lu:=evalf(evalc(pit[makom[1]])): if coeff(lu,I,1)>(1/10)^(Digits-2) then print(`Dominant roots are complex`): RETURN(FAIL): elif pit[makom[1]]+pit[makom[2]]=0 then print(`Dominant real roots are negatives of each other`): RETURN(FAIL): fi: fi: lu:=evalf(evalc(pit[makom[1]])): if coeff(lu,I,1)>(1/10)^(Digits-2) then print(`Dominant root is complex`): RETURN(FAIL): fi: lu:=pit[makom[1]]: ope1:=Yafe(subs(N=lu*N,ope),N)[2]: ope1,lu: end: #Atom(s,r,i,x,K): Expanding (n+i)^(s/r)-n^(s/r) in terms of #x=n^(-1/r) up to K terms #For example try: #Atom(2,3,2,x,3); Atom:=proc(s,r,i,x,K) local gu,y,i1: gu:=(1+i*y)^(s/r)-1: gu:=taylor(gu,y=0,K+1): gu:=add(coeff(gu,y,i1)*y^i1,i1=0..K): expand(subs(y=x^r,gu)/x^s): end: #FindExpP1(ope,n,N,r,x): finds the exponential part of the asymptotics #in terms of x=n^(1/r) #for the solutions of ope(n,N)a(n)=0 if it is of type r #For example, try: #FindExpP1((n+1)*N^2-2*(n+5)*N+n+3,n,N,2,x); FindExpP1old:=proc(ope,n,N,r,x) local c,eq,var,s,sof,gu,ope1,deg,i,i1,v: sof:=add(c[s]*x^s,s=1..r-1): var:={seq(c[i],i=1..r-1)}: deg:=degree(ope,n): ope1:=expand(ope/n^deg): ope1:=expand(subs(n=1/x^r,ope1)): gu:=0: for i from 0 to degree(ope1,N) do gu:=gu+coeff(ope1,N,i)*exp(add(c[s]*Atom(s,r,i,x,2),s=1..r-1) ): od: gu:=taylor(gu,x=0,2*r-1): eq:={seq(coeff(gu,x,i1),i1=r..2*r-2)}: var:=[solve(eq,var)][1]: for v in var do if op(1,v)=op(2,v) then RETURN(FAIL): fi: od: subs(var,sof): end: #Findr(ope,n,N): Given an operator ope(n,N) #such that the leading terms in n is a multiple of #(N-1), finds the largest r such that it is # a multiple of (N-1)^r. For example, try: #Findr((N-1)^4,n,N); Findr:=proc(ope,n,N) local ope1,r: ope1:=coeff(expand(ope),n,degree(ope,n)): if expand(subs(N=1,ope1))<>0 then RETURN(FAIL): fi: for r from 1 while expand(subs(N=1,diff(ope1,N$r)))=0 do od: r: end: FindExpP:=proc(ope,n,N) local r,x,lu: r:=Findr(ope,n,N): if r=FAIL then RETURN(FAIL): fi: lu:=FindExpP1(ope,n,N,r,x) : if lu=FAIL then RETURN(FAIL): fi: subs(x=n^(1/r),lu): end: #NewOpe(ope,n,N,K): the exponential part + the transformed operator #(up to degree K asymp. in the coefficients) #For example, try: #NewOpe((n+1)*N^2-2*(n+5)*N+n+3,n,N,4); NewOpe:=proc(ope,n,N,K,x) local r,lu,lu1,ope1,deg,s,i,gu,mu: r:=Findr(ope,n,N): if r=FAIL then RETURN(FAIL): fi: lu:=FindExpP1(ope,n,N,r,x) : lu1:=FindExpP(ope,n,N) : if lu=FAIL then RETURN(FAIL): fi: deg:=degree(ope,n): ope1:=expand(ope/n^deg): ope1:=subs(n=1/x^r,ope1): ope1:=expand(ope1): gu:=0: for i from 0 to degree(ope1,N) do mu:=0: for s from 1 to r-1 do mu:=mu+coeff(lu,x,s)*Atom(s,r,i,x,K+3): od: gu:=gu+coeff(ope1,N,i)*exp(mu)*N^i: od: gu:=taylor(gu,x=0,K+3): lu1,add(coeff(gu,x,i)*x^i,i=0..K): end: #Bdok(ope,n,N,K): the exponential part + the transformed operator #(up to degree K asymp. in the coefficients) #For example, try: #Bdok((n+1)*N^2-2*(n+5)*N+n+3,n,N,4); Bdok:=proc(ope,n,N,K) local r,x,lu,lu1,ope1,deg,s,i,gu,mu: r:=Findr(ope,n,N): if r=FAIL then RETURN(FAIL): fi: lu:=FindExpP1(ope,n,N,r,x) : lu1:=FindExpP(ope,n,N) : if lu=FAIL then RETURN(FAIL): fi: deg:=degree(ope,n): ope1:=expand(ope/n^deg): ope1:=subs(n=1/x^r,ope1): ope1:=expand(ope1): gu:=0: for i from 0 to degree(ope1,N) do mu:=0: for s from 1 to r-1 do mu:=mu+coeff(lu,x,s)*Atom(s,r,i,x,K+3): od: gu:=gu+coeff(ope1,N,i)*exp(mu)*N^i: od: gu:=taylor(gu,x=0,K+3): lu1,add(coeff(gu,x,i)*x^i,i=0..K): end: OpePerN:=proc(N,n,r) local i,j,ope: ope:=1-add(N^(-i)*mul(n-j,j=1..i-1),i=1..r): ope:=expand(N^r*subs(n=n+r,ope)): ope:=add(factor(coeff(ope,N,i))*N^i,i=0..r): ope:=expand(FindKk(ope,n,N)[2]): ope:=add(factor(coeff(ope,N,i))*N^i,i=0..r): ope:=numer(Nor(ope,n,N)[1]): ope:=add(factor(coeff(ope,N,i))*N^i,i=0..r): end: OpePer:=proc(N,n,r) local i,j,ope: ope:=1-add(N^(-i)*mul(n-j,j=1..i-1),i=1..r): ope:=expand(N^r*subs(n=n+r,ope)): end: OpePerG:=proc(N,n,S) local i,j,ope,r,gu,t: ope:=1-add(N^(-i)*mul(n-j,j=1..i-1),i in S): r:=max(op(S)): ope:=expand(N^r*subs(n=n+r,ope)): gu:=exp(add(t^i/i,i in S)): gu:=taylor(gu,t=0,degree(ope,N)+2): ope,[seq(coeff(gu,t,i)*i!,i=1..degree(ope,N))]: end: #Finda(ope,N,x,r): finds the first power x^a in the #asymptotic solution of ope(N,n)f(n)=0, where x=1/n^(1/r) #For example, try: #Finda((1+x)-(1+3*x)*N,N,x,1); Finda:=proc(ope,N,x,r) local gu,i,a,a1: gu:=add(coeff(ope,N,i)*(1+i*x^r)^(a/r),i=0..degree(ope,N)): gu:=taylor(gu,x=0,5*r+5): gu:=expand(add(coeff(gu,x,i)*x^i,i=0..5*r+4)): for i from 0 to 5*r+4 while expand(simplify(coeff(gu,x,i)))=0 do od: if i=5*r+5 then RETURN(FAIL): fi: a1:=[solve(coeff(gu,x,i),a)][1]: if a1=NULL then RETURN(FAIL): elif not type(a1,integer) then RETURN(FAIL): else RETURN(-a1): fi: end: #OneStepG(ope,N,x,r,Cu): Given a partial #asympt. expansion, in terms of x=1/n^(1/r), #for solutions of the recurrence equation #ope(n,N)a(n)=0, where ope is normalized of type #r (i.e. its leading coeff. is a multiple of (N-1)^r) #finds the next term #OneStepG((1+x)-(1+3*x)*N,N,x,1,x^2); OneStepG:=proc(ope,N,x,r,Cu) local gu,i,a,c,c1,Cu1,K,ka: K:=degree(Cu,x): Cu1:=Cu+c*x^(K+1): gu:=add(coeff(ope,N,i)* add(coeff(Cu1,x,a)*x^a*(1+i*x^r)^(-a/r),a=ldegree(Cu1,x)..degree(Cu1,x)), i=0..degree(ope,N)): gu:=expand(gu): gu:=taylor(gu,x=0,5*r+8+K): gu:= simplify(expand(add(coeff(gu,x,i)*x^i,i=0..5*r+7+K))): gu:=pashet(gu,x,5*r+7+K,c): ka:=expand(simplify(coeff(gu,x,ldegree(gu,x)))): ka:=simplify(ka): c1:=[solve(simplify(coeff(gu,x,ldegree(gu,x))),c)][1]: if c1=NULL then RETURN(FAIL): else RETURN(subs(c=c1,Cu1)): fi: end: #Asy1(ope,n,N,K): the asymptotic expansion of solutions #to ope(n,N)f(n)=0,where ope(n,N) is a recurrence operator #up to the K's term Asy1:=proc(ope,n,N,K) local gu,lu,alpha,mu,ope1,ku,i,f,x,ka,vu,ope1A,r,Kk,pu,a: gu:=FindKk(ope,n,N): Kk:=gu[1]: ope1:=gu[2]: vu:=Nor(ope1,n,N): if vu=FAIL then RETURN(FAIL): fi: ope1:=vu[1]: lu:=vu[2]: ope1A:=Aope1(ope1,n,N): for r from 1 while subs(N=1,diff(ope1A,N$r))=0 do od: if r=1 then ###r=1 case mu:=add( subs(n=1/x,coeff(ope1,N,i))*(1+i*x)^alpha,i=ldegree(ope1,N)..degree(ope1,N)): mu:=normal(mu): ku:=factor(coeff(taylor(mu,x=0,2),x,1)): ka:=simplify([solve(ku,alpha)]): if coeff(ka[1],I,1)<>0 then RETURN(FAIL): fi: alpha:=max(op(ka)): if normal(subs(N=1,ope1))=0 then RETURN(lu^n*n^alpha,Kk): fi: f:=1: for i from 1 to K do f:=OneStepA(ope1,n,N,alpha,f): od: if f=FAIL then RETURN(FAIL): fi: RETURN(lu^n*n^alpha*f,Kk): #end r=1 case else ope1:=numer(ope1): if expand(diff(ope1,n))=0 then ka:=n^(r-1)*lu^n: RETURN(ka,Kk): fi: pu:=FindExpP(ope1,n,N): if pu=FAIL then RETURN(FAIL): fi: ope1:=NewOpe(ope1,n,N,K+5,x)[2]: a:=Finda(ope1,N,x,r): if a=FAIL then RETURN(FAIL): fi: ka:=x^a: for i from a to K do ka:=OneStepG(ope1,N,x,r,ka): od: ka:=exp(pu)*subs(x=1/n^(1/r),ka)*lu^n: RETURN(ka,Kk): fi: end: pashet:=proc(gu,x,K,c) local gu1,i,pu: Digits:=300: gu1:=0: for i from 0 to K do pu:=evalf(coeff(gu,x,i)): if not (abs(coeff(pu,c,1))<10^(-20) and abs(coeff(pu,c,0))<10^(-20) ) then gu1:=gu1+coeff(gu,x,i)*x^i: fi: od: gu1: end: #NakedStirling(n,K): the asymptotic expansion of #n!/((n/e)^n*sqrt(2*Pi*n). For example, try: #NakedStirling(n,5); NakedStirling:=proc(n,K) local gu,i: gu:=n!/(n/exp(1))^n/sqrt(2*Pi*n): gu:=expand(asympt(gu,n,K+1)): (n/exp(1))^n*sqrt(2*Pi*n), add(coeff(op(i,gu),n,-(i-1))/n^(i-1),i=1..K+1): end: #AsyF(ope,n,N,M): the asymptotic expansion of solutions #to ope(n,N)f(n)=0,where ope(n,N) is a recurrence operator #up to the M's term, in terms of the factorial function AsyF:=proc(ope,n,N,M) local K,k,gu: gu:=Asy1(ope,n,N,M): if gu=FAIL then RETURN(FAIL): fi: K:=gu[2][1]: k:=gu[2][2]: gu:=gu[1]: (n/K)!^k*subs(n=n/K,gu): end: #AsyFC(ope,n,N,M,Ini,K): the asymptotic expansion of solutions #to ope(n,N)f(n)=0, with the given initial conditions #where ope(n,N) is a recurrence operator #up to the M's term, in terms of the factorial function #and complete with an empirically derived constant in front #using K terms AsyFC:=proc(ope,n,N,M,Ini,K) local gu,L,i,mu,er,C,D1: Digits:=100: gu:=AsyF(ope,n,N,M): L:=SeqFromRec(ope,n,N,evalf(Ini),K): mu:=[seq(evalf(L[i]/subs(n=i,gu)),i=K-10..K)]: er:=mu[nops(mu)]-mu[nops(mu)-1]: D1:=-trunc(log(abs(er))/log(10))-3: if D1<2 then print(`can't determine the constant`): RETURN(gu): fi: C:=evalf(mu[nops(mu)],D1): C,gu: end: ###Asy #Asy1special(ope,n,N,K,x): the asymptotic expansion of solutions #to ope(n,N)f(n)=0,given as a list #[pu,lu,expansion,r] where it is #exp(pu)*lu^n*expansion(x), where x=1/n^(1/r), and r is #a positive integer. It also returns [K,k] (see Asy1) #where ope(n,N) is a recurrence operator #up to the K's term Asy1special:=proc(ope,n,N,K,x) local gu,lu,alpha,mu,ope1,ku,i,f,ka,vu,ope1A,r,Kk,pu,a: gu:=FindKk(ope,n,N): Kk:=gu[1]: ope1:=gu[2]: vu:=Nor(ope1,n,N): if vu=FAIL then RETURN(FAIL): fi: ope1:=vu[1]: lu:=vu[2]: ope1A:=Aope1(ope1,n,N): for r from 1 while subs(N=1,diff(ope1A,N$r))=0 do od: if degree(ope1,n)=0 then RETURN([0,lu,x^(-r+1),1,1],Kk): fi: if r=1 then ###r=1 case mu:=add( subs(n=1/x,coeff(ope1,N,i))*(1+i*x)^alpha,i=ldegree(ope1,N)..degree(ope1,N)): mu:=normal(mu): ku:=factor(coeff(taylor(mu,x=0,2),x,1)): ka:=simplify([solve(ku,alpha)]): if coeff(ka[1],I,1)<>0 then RETURN(FAIL): fi: alpha:=max(op(ka)): if normal(subs(N=1,ope1))=0 then RETURN([0,lu,x^(-alpha),1,1],Kk): fi: f:=1: for i from 1 to K do f:=OneStepA(ope1,n,N,alpha,f): od: if f=FAIL then RETURN(FAIL): fi: f:=subs(n=1/x,f): RETURN([0,lu,x^(-alpha),f,1],Kk): #end r=1 case else ope1:=numer(ope1): if expand(diff(ope1,n))=0 then ka:=n^(r-1)*lu^n: RETURN(ka,Kk): fi: pu:=FindExpP(ope1,n,N): if pu=FAIL then RETURN(FAIL): fi: ope1:=NewOpe(ope1,n,N,K+5,x)[2]: a:=Finda(ope1,N,x,r): if a=FAIL then RETURN(FAIL): fi: ka:=x^a: for i from a to K do ka:=OneStepG(ope1,N,x,r,ka): od: RETURN([pu,lu,1,ka,r],Kk): fi: end: #IsMispar(a): is a (generalized) numeric type IsMispar:=proc(a) : if type(a,numeric) then true: elif type(a,`^`) and type(op(1,a),numeric) and type(op(2,a),numeric) then true: elif a=Pi then true: elif type(a,`^`) and op(1,a)=Pi and type(op(2,a),numeric) then true: elif type(a,function) and type(op(1,a),numeric) then true: else false: fi: end: #GetRidOfConst(P): gets rid of the constant in front GetRidOfConst:=proc(P) local i,P1: if not type(P,`*`) then RETURN(P): fi: P1:=1: for i from 1 to nops(P) do if not IsMispar(op(i,P)) then P1:=P1*op(i,P): fi: od: P1: end: #Asy(ope,n,N,M): the asymptotic expansion of solutions #to ope(n,N)f(n)=0,where ope(n,N) is a recurrence operator #up to the M's term Asy:=proc(ope,n,N,M) local K,k,gu,vu,vu1,vu2,pu,lu,ka1,ka2,r,x,i,vu2a: gu:=Asy1special(ope,n,N,M,x): if gu=FAIL then RETURN(FAIL): fi: K:=gu[2][1]: k:=gu[2][2]: gu:=gu[1]: pu:=subs(n=n/K,gu[1]): lu:=gu[2]: ka1:=gu[3]: ka2:=gu[4]: r:=gu[5]: ka2:=subs(x=K^(1/r)*x,ka2): vu:=NakedStirling(n,M+1): vu1:=vu[1]: vu2:=vu[2]: vu1:=subs(n=n/K,vu1)^k: vu2:=subs(n=n/K,vu2)^k: vu2:=expand(subs(n=1/x^r,vu2)): vu2:=expand(ka2*vu2): vu2:=add(simplify(coeff(vu2,x,i))*x^i,i=0..M*r): vu2:=subs(x=1/n^(1/r),vu2): ka1:=subs(x=1/n^(1/r),ka1): ka1:=subs(n=n/K,ka1): vu2a:=op(1,vu2): vu2:=expand(normal(vu2/vu2a)): gu:=simplify(vu2a*vu1*lu^(n/K))*exp(pu)*ka1*vu2: GetRidOfConst(gu): end: #AsyC(ope,n,N,M,Ini,K): the asymptotic expansion of solutions #to ope(n,N)f(n)=0, with the given initial conditions #where ope(n,N) is a recurrence operator #up to the M's term, #and complete with an empirically derived constant in front #using K terms AsyC:=proc(ope,n,N,M,Ini,K) local gu,L,i,mu,er,C,D1: Digits:=100: gu:=Asy(ope,n,N,M): if gu=FAIL then RETURN(FAIL): fi: L:=SeqFromRec(ope,n,N,evalf(Ini),K): mu:=[seq(evalf(L[i]/subs(n=i,gu)),i=K-10..K)]: er:=mu[nops(mu)]-mu[nops(mu)-1]: if er=0 then D1:=Digits: else D1:=-trunc(log(evalf(abs(er)))/log(10.))-3: fi: if D1<2 then print(`can't determine the constant`): RETURN(gu): fi: C:=evalf(mu[nops(mu)],D1): C,gu: end: #FindExpP1g(ope,n,N,r,x,ds): finds the exponential part of the asymptotics #in terms of x=n^(1/r) as a poly. of degree ds in x #for the solutions of ope(n,N)a(n)=0 if it is of type r #For example, try: #FindExpP1g((n+1)*N^2-2*(n+5)*N+n+3,n,N,2,x,1); FindExpP1g:=proc(ope,n,N,r,x,ds) local c,eq,var,s,sof,gu,ope1, deg,i,i1,v,ka,i11: sof:=add(c[s]*x^s,s=1..ds): var:={seq(c[i],i=1..ds)}: deg:=degree(ope,n): ope1:=expand(ope/n^deg): ope1:=expand(subs(n=1/x^r,ope1)): gu:=0: for i from 0 to degree(ope1,N) do gu:=gu+coeff(ope1,N,i)*exp(add(c[s]*Atom(s,r,i,x,2),s=1..ds) ): od: gu:=taylor(gu,x=0,3*r+10): for i1 from 0 to 3*r+8 while coeff(gu,x,i1)=0 do od: eq:={seq(coeff(gu,x,i11),i11=i1..i1+ds-1)}: eq:={seq(coeff(gu,x,i11),i11=i1..i1+ds)}: lprint(eq,var): ka:=[solve(eq,var)]: if ka=[] then RETURN(FAIL): fi: var:=ka[1]: for v in var do if op(1,v)=op(2,v) then RETURN(FAIL): fi: od: subs(var,sof): end: #FindExpP1(ope,n,N,r,x): finds the exponential part of the asymptotics #in terms of x=n^(1/r) #for the solutions of ope(n,N)a(n)=0 if it is of type r #For example, try: #FindExpP1((n+1)*N^2-2*(n+5)*N+n+3,n,N,2,x); FindExpP1:=proc(ope,n,N,r,x) local s,gu,nosaf: for s from 1 to r-1 do for nosaf from 0 to 1 do gu:=FindExpP1gExtra(ope,n,N,r,x,s,nosaf): if gu<>FAIL then RETURN(gu): fi: od: od: FAIL: end: #FindExpP1gExtra(ope,n,N,r,x,ds,nosaf) #: finds the exponential part of the asymptotics #in terms of x=n^(1/r) as a poly. of degree ds in x #for the solutions of ope(n,N)a(n)=0 if it is of type r #For example, try: #FindExpP1gExtra((n+1)*N^2-2*(n+5)*N+n+3,n,N,2,x,1); FindExpP1gExtra:=proc(ope,n,N,r,x,ds,nosaf) local c,eq,var,s,sof,gu,ope1, deg,i,i1,v,ka,i11,ku,varf: sof:=add(c[s]*x^s,s=1..ds): var:={seq(c[i],i=1..ds)}: deg:=degree(ope,n): ope1:=expand(ope/n^deg): ope1:=expand(subs(n=1/x^r,ope1)): gu:=0: for i from 0 to degree(ope1,N) do gu:=gu+coeff(ope1,N,i)*exp(add(c[s]*Atom(s,r,i,x,2),s=1..ds) ): od: gu:=taylor(gu,x=0,3*r+10): for i1 from 0 to 3*r+8 while coeff(gu,x,i1)=0 do od: eq:={seq(coeff(gu,x,i11),i11=i1..i1+ds-1+nosaf)}: ka:=[solve(eq,var)]: if ka=[] then RETURN(FAIL): fi: var:=ka[1]: for v in var do if op(1,v)=op(2,v) then RETURN(FAIL): fi: od: #varf:=evalf(var): #ku:=[seq(subs(varf,c[s]),s=1..ds)]: #print(`ku is`, ku): #add(ku[s]*x^s,s=1..ds): subs(var,sof): end: ###end from AsyRec #Asyabm(a,b,m,M,K,n): An asymptotic formula, to order M #(using K terms of the sequence) for Prnm(a*n,b*n,m); #For example, try: #Asyabm(1,1,4,3,1000,n); Asyabm:=proc(a,b,m,M,K,n) local N,ope,lu,Ini,i,ku,C,gu,ru: option remember: if not {type(a,integer), type(b,integer), type(m,integer)}={true} then print(`Bad input`): RETURN(FAIL): fi: if m<2 then print(`the third argument`, m, `should be at least 2 `): fi: ope:=Recabm(a,b,m,n,N): Ini:=[seq(Arnm(a*i,b*i,m),i=1..degree(ope,N))]: lu:=AsyC(ope,n,N,M,Ini,K): if lu=FAIL then RETURN(FAIL): fi: C:=lu[1]: lu:=evalf(lu[2]): gu:=lu: ku:=1: for i from 1 to nops(lu) do if type(op(i,lu),`^`) and type(op(1,op(i,lu)),float) then ku:=ku*op(i,lu): gu:=gu/op(i,lu): fi: od: ku:=simplify(ku): ru:=asympt(exp(n)*n!/n^n,n,M): ru:=ru-op(nops(ru),ru): ru:=subs(n=n*a,ru): if type(gu,`^`) then RETURN(ru*gu): fi: gu:=expand(ru*gu): if not type(gu,`+`) then print(`Something is wrong`): RETURN(FAIL): fi: gu:=add(C*coeff(gu,n,-i)/n^i,i=0..M): ku:=simplify(ku*exp(-n*a)*((a/b)^a)^n): evalf(ku*gu): end: #AsyabmPA(a,b,m,n): A fast version of Asyabm(a,b,m,M,K) #using the Poisson approxination and not completely #justified assumtions of "independence" #For example, try: #AsyabmPA(1,1,4,n); AsyabmPA:=proc(a,b,m,n) local c,r,i: c:=a/b: c:=evalf(exp(-c)*add(c^i/i!,i=0..m)): c:=evalf(log(c)): exp(c*n): end: #AsyabmGPA(a,b,m1,m2,n): A fast version of AsyabmG(a,b,m1,m2,M,K) #using the Poisson approxination and not completely #justified assumtions of "independence" #For example, try: #AsyabmGPA(2,1,1,4,n); AsyabmGPA:=proc(a,b,m1,m2,n) local c,r,i: c:=a/b: c:=evalf(exp(-c)*add(c^i/i!,i=m1..m2)): c:=evalf(log(c)): exp(c*n): end: #RecabmV(a,b,m,n,N): Verbose version of Recabm(a,b,m,n,N) #For example, try: #RecabmV(2,1,3,n,N): RecabmV:=proc(a,b,m,n,N) local ope,Ini,P,Q,i: print(`Theorem: Let P(n) be the probability that if you place (uniformly)`): print(a*n,`balls in `, b*n , `boxes then no box has more than `, m, `balls. `): print(`Then we have: `): print(P(n)=(a*n)!/(b*n)^(a*n)*Q(n)): print(``): ope:=Recabm(a,b,m,n,N): if ope=FAIL then RETURN(FAIL): fi: Ini:=[seq(Arnm(a*i,b*i,m),i=1..degree(ope,N))]: print(`where Q(n) is the sequence satisfying the recurrece`): print(add(coeff(ope,N,i)*Q(n+i),i=0..degree(ope,N))=0): print(``): print(`subject to the initial conditions`): print(seq(Q(i)=Ini[i],i=1..nops(Ini))): print(``): print(`In operator notation, using `, N, `as the shift operator in `, n): print(`the recurrence operator annihilating Q(n) is`): lprint(ope): end: #AsyabmV(a,b,m,M,K,n): Verbose version of #Asyabm(a,b,m,M,K,n) and with a theorem number #For example do #AsyabmV(1,1,4,5,5000,n); AsyabmV:=proc(a,b,m,M,K,n) local lu,mu,ku: lu:=Asyabm(a,b,m,M,K,n): print(`Theorem: Let P(n) be the probability that if you `): print(`place (uniformly)`, a*n,`balls in `, b*n , `boxes then `): print(`no box has more than `, m, `balls. `): print(`Then we have the following asymptotic formula to order `, M): print(P(n)=lu): print(``): print(`Comment 1.: compare this EXACT result with the naive`): print(` Poisson approximation`): mu:=AsyabmPA(a,b,m,n): print(mu): print(`Comment 2: Let's see how good the asymptotic formula is`): print(`For n=100, it gives you`): print(evalf(subs(n=100,lu))): print(`whereas the exact value, using the Ewens-Wilf approach is`): ku:=evalf(PrnmExact(100*a,100*b,m)): print(ku): print(`The ratio is:`, evalf(subs(n=100,lu)/ku)): print(`For n=300, it gives you`): print(evalf(subs(n=300,lu))): print(`whereas the exact value, using the Ewens-Wilf approach is`): ku:=Prnm(300*a,300*b,m): print(ku): print(`The ratio is:`, evalf(subs(n=300,lu)/ku)): print(``): print(`In Maple input format, the asymptotic formula is`): lprint(lu): end: #AsyabmV1(a,b,m,M,K,n,co): Verbose version of #Asyabm(a,b,m,M,K,n) and with a theorem number #For example do #Asyabm(1,1,4,5,5000,n); AsyabmV1:=proc(a,b,m,M,K,n,co) local lu,mu,ku: lu:=Asyabm(a,b,m,M,K,n): printf("Theorem Number %d :\n" , co): print(`Let P(n) be the probability that if you `): print(`place (uniformly)`, a*n,`balls in `, b*n , `boxes then `): print(`no box has more than `, m, `balls. `): print(`Then we have the following asymptotic formula to order `, M): print(P(n)=lu): print(``): print(`Comment 1.: compare this EXACT result with the naive`): print(` Poisson approximation`): mu:=AsyabmPA(a,b,m,n): print(mu): print(`Comment 2: Let's see how good the asymptotic formula is`): print(`For n=100, it gives you`): print(evalf(subs(n=100,lu))): print(`whereas the exact value, using the Ewens-Wilf approach is`): ku:=evalf(PrnmExact(100*a,100*b,m)): print(ku): print(`The ratio is:`, evalf(subs(n=100,lu)/ku)): print(`For n=300, it gives you`): print(evalf(subs(n=300,lu))): print(`whereas the exact value, using the Ewens-Wilf approach is`): ku:=Prnm(300*a,300*b,m): print(ku): print(`The ratio is:`, evalf(subs(n=300,lu)/ku)): print(``): print(`In Maple input format, the asymptotic formula is`): lprint(lu): print(`------------------------------`): print(``): end: #AsyDB(Maxa,Maxb,Maxm,M,K,n): creates a data base for the #asymptotic formulas, in n, of order M, for the probability of no box #getting more than m boxes when you place a*n balls in #b*n boxes for a and b relatively prime a<=Maxa, b<=Maxb #m<=Maxm. K is a large positive integer for estimation purposes. #The output is a list of pairs [[a,b,m],formula] #For example, try: #AsyDB(2,2,4,4,5000,n); AsyDB:=proc(Maxa,Maxb,Maxm,M,K,n) local a,b,yakhas,kv,m,lu,gu,a1,b1: kv:={}: gu:={}: for a from 1 to Maxa do for b from 1 to Maxb do yakhas:=a/b: a1:=numer(yakhas): b1:=denom(yakhas): if not member([a1,b1],kv) then kv:=kv union {[a1,b1]}: for m from max(2,trunc(a/b)+1) to Maxm do lu:=Asyabm(a1,b1,m,M,K,n): if lu<>FAIL then gu:=[op(gu),[[a1,b1,m],lu]]: fi: od: fi: od: od: gu: end: #CFAveab(a,b,n,eps,d,L,Sta): guesses a polynomial in log(n) of degree d #that best fits the quantity: "Average Number of the maximum number #of balls that lend in any box" if a*n balls are placed (uniformly #at random) in b*n boxes. using L data points. #eps should be small (we only cover the 1-eps part of the distr. ). #Sta is the starting point #For example, try: #CFAveab(1,1,n,10^(-10),1,5,1000); CFAveab:=proc(a,b,n,eps,d,L,Sta) local gu,k,v,u: gu:=[]: for k from 0 to L-1 do v:=Pdistr(2^k*Sta*a,2^k*Sta*b,eps,2)[2][1]: gu:=[op(gu),[evalf(log(2^k*Sta)),v]]: od: gu:=CurveFitting[LeastSquares](gu,u): subs(u=log(n),gu): end: #CFAveabEmpir(a,b,n,eps,d,L,Sta,K): Like #CFAveab(a,b,n,eps,d,L,Sta): but using PdistrEmpir instead of #Pdistr instead #For example, try: #CFAveabEmpir(1,1,n,10^(-10),1,5,1000,50); CFAveabEmpir:=proc(a,b,n,eps,d,L,Sta,K) local gu,k,v,u: gu:=[]: for k from 0 to L-1 do v:=PdistrEmpir(2^k*Sta*a,2^k*Sta*b,eps,2,K,2*b)[2][1]: gu:=[op(gu),[evalf(log(2^k*Sta)),v]]: od: gu:=CurveFitting[LeastSquares](gu,u): subs(u=log(n),gu): end: #CFAveabPA(a,b,n,eps,d,L,Sta): guesses a polynomial in log(n) of degree d #that best fits the quantity: "Average Number of the maximum number #of balls that lend in any box" if a*n balls are placed (uniformly #at random) in b*n boxes. using L data points. #eps should be small (we only cover the 1-eps part of the distr. ). #For example, try: #CFAveabPA(1,1,n,10^(-10),1,5,1000); CFAveabPA:=proc(a,b,n,eps,d,L,Sta) local gu,k,v,u,a1: gu:=[]: for k from 0 to L-1 do v:=PdistrPA(2^k*Sta*a,2^k*Sta*b,eps,2)[2][1]: gu:=[op(gu),[evalf(log(2^k*Sta)),v]]: od: gu:=CurveFitting[LeastSquares](gu,u,curve=add(a1[i]*u^i,i=0..d)): subs(u=log(n),gu): end: #CFMomab(a,b,n,eps,d,L,Sta,K): guesses a polynomial in log(n) of degree d #that best fits the quantity: "Average Number of the maximum number #of balls that lend in any box" and the standard deviation and the #higher (normalized) moments, up to the Kth #if a*n balls are placed (uniformly #at random) in b*n boxes. using L data points. #eps should be small (we only cover the 1-eps part of the distr. ). #For example, try: #CFMomab(1,1,n,10^(-10),1,5,1000,4); CFMomab:=proc(a,b,n,eps,d,L,Sta,K) local gu,k,v,u,a1,j,i: for j from 1 to K do gu[j]:=[]: od: for k from 0 to L-1 do v:=Pdistr(2^k*Sta*a,2^k*Sta*b,eps,K)[2]: for j from 1 to K do gu[j]:=[op(gu[j]),[evalf(log(2^k*Sta)),v[j]]]: od: od: for j from 1 to K do gu[j]:=CurveFitting[LeastSquares](gu[j],u,curve=add(a1[i]*u^i,i=0..d)): gu[j]:=subs(u=log(n),gu[j]): od: [seq(gu[j],j=1..K)]: end: #CFMomabPA(a,b,n,eps,d,L,Sta,K): guesses a polynomial in log(n) of degree d #that best fits the quantity: "Average Number of the maximum number #of balls that lend in any box" and the standard deviation and the #higher (normalized) moments, up to the Kth #if a*n balls are placed (uniformly #at random) in b*n boxes. using L data points. #eps should be small (we only cover the 1-eps part of the distr. ). #For example, try: #CFMomabPA(1,1,n,10^(-10),1,5,1000,4); CFMomabPA:=proc(a,b,n,eps,d,L,Sta,K) local gu,k,v,u,a1,j,i: for j from 1 to K do gu[j]:=[]: od: for k from 0 to L-1 do v:=PdistrPA(2^k*Sta*a,2^k*Sta*b,eps,K)[2]: for j from 1 to K do gu[j]:=[op(gu[j]),[evalf(log(2^k*Sta)),v[j]]]: od: od: for j from 1 to K do gu[j]:=CurveFitting[LeastSquares](gu[j],u,curve=add(a1[i]*u^i,i=0..d)): gu[j]:=subs(u=log(n),gu[j]): od: [seq(gu[j],j=1..K)]: end: #CFAveabPANor(a,b,n,eps,d,L,Sta): guesses a polynomial in log(n) of degree d #that best fits the quantity: "Average Number of the maximum number #of balls that lend in any box" if (a/b)*n balls are placed (uniformly #at random) in n boxes. using L data points. #eps should be small (we only cover the 1-eps part of the distr. ). #For example, try: #CFAveabPANor(1,1,n,10^(-10),1,5,1000); CFAveabPANor:=proc(a,b,n,eps,d,L,Sta) local gu: gu:=CFAveabPA(a,b,n,eps,d,L,Sta): evalf(expand(subs(n=n/b,gu))): end: #EmpForAve(n,R,eps,d,L,Sta,d1,M1,M2): An empirically-derived formula #as a polynomial of degree d in log(n) and of degree d1 #in R, #using Least Squares and the Poisson Approximation #for the quantity: #The average of the random variable "Max. Number of Balls that lended #in an individual box" upon placing, uniformaly, n*R balls in n boxes #the data set is from CFAveabPANor(a,b,n,eps,d,L,Sta) with a between #1 and M1 and b between 1 and M2 #For example, try: EmpForAve(n,R,10^(-10),1,10,100,1,3,3); EmpForAve:=proc(n,R,eps,d,L,Sta,d1,M1,M2) local gu,a,b,j,lu,mu,i: gu:=[]: for a from 1 to M1 do for b from 1 to M2 do if gcd(a,b)=1 then gu:=[op(gu), [a/b,CFAveabPANor(a,b,n,eps,d,L,Sta)]]: fi: od: od: mu:=0: for j from 0 to d do lu:=[seq([gu[i][1],coeff(gu[i][2],log(n),j)],i=1..nops(gu))]: lu:=CurveFitting[LeastSquares](lu,R,curve=add(a1[i]*R^i,i=0..d1)): mu:=mu+lu*log(n)^j: od: mu: end: #CFMomabPANor(a,b,n,eps,d,L,Sta,K): like CFMombPA(a,b,n,eps,d,L,Sta) #but normalized bn->n #For example, try: #CFMomabPANor(1,1,n,10^(-10),1,5,1000,6); CFMomabPANor:=proc(a,b,n,eps,d,L,Sta,K) local gu: gu:=CFMomabPA(a,b,n,eps,d,L,Sta,K): evalf(expand(subs(n=n/b,gu))): end: #EmpForMom(n,R,eps,L,Sta,K,M1,M2,d,d1): An empirically-derived formula #as a polynomial of degree d in log(n) and of degree d1 #in R, #using Least Squares and the Poisson Approximation #for the quantity: #The average of the random variable "Max. Number of Balls that lended #in an individual box" upon placing, uniformaly, n*R balls in n boxes #the standard deviation, and the higher normalized moments to order K. #The data set is from CFMomabPANor(a,b,n,eps,d,L,Sta,K) with a between #1 and M1 and b between 1 and M2 #For example, try: EmpForMom(n,R,10^(-10),10,1000,4,3,3,1,1); EmpForMom:=proc(n,R,eps,L,Sta,K,M1,M2,d,d1) local gu,a,b,j,lu,mu,i,pu,j1: gu:=[]: for a from 1 to M1 do for b from 1 to M2 do if gcd(a,b)=1 then gu:=[op(gu), [a/b,CFMomabPANor(a,b,n,eps,d,L,Sta,K)]]: fi: od: od: pu:=[]: for j1 from 1 to K do mu:=0: for j from 0 to d do lu:=[seq([gu[i][1],coeff(gu[i][2][j1],log(n),j)],i=1..nops(gu))]: lu:=CurveFitting[LeastSquares](lu,R,curve=add(a1[i]*R^i,i=0..d1)): mu:=mu+lu*log(n)^j: od: pu:=[op(pu),mu]: od: pu: end: #HowGoodPA1(R0,N0,Incr,M0,m,eps): explores the ratio of #PrnmPA(r,n,m) to Prnm(r,n,m) for r in [R0, R0+Incr*M0] #and n in [N0,N0+Incr*M0] for a specfic m provided #Prnm(r,n,m) is larger than eps. #Returns the one the ratio furthest away from 1 together with the #(n,r) that yield it #For example, try: #HowGoodPA1(1000,1000,1000,3,5,1/100); HowGoodPA1:=proc(R0,N0,Incr,M0,m,eps) local n,r,aluf, si,lu,vu: si:=0: for n from N0 by Incr to N0+M0*Incr do for r from R0 by Incr to R0+M0*Incr do lu:=PrnmPA(r,n,m): if lu>eps then vu:=abs(1-lu/Prnm(r,n,m)): if vu>si then si:=vu: aluf:=[r,n]: fi: fi: od: od: si,aluf: end: #HowGoodPA1exact(R0,N0,Incr,M0,m,eps): explores the ratio of #PrnmPA(r,n,m) to Prnm(r,n,m) for r in [R0, R0+Incr*M0] #and n in [N0,N0+Incr*M0] for a specfic m provided #Prnm(r,n,m) is larger than eps. #Returns the one the ratio furthest away from 1 together with the #(n,r) that yield it #For example, try: #HowGoodPA1exact(1000,1000,1000,3,5,1/100); HowGoodPA1exact:=proc(R0,N0,Incr,M0,m,eps) local n,r,aluf, si,lu,vu: si:=0: for n from N0 by Incr to N0+M0*Incr do for r from R0 by Incr to R0+M0*Incr do lu:=PrnmPA(r,n,m): if lu>eps then vu:=abs(1-lu/evalf(PrnmExact(r,n,m))): if vu>si then si:=vu: aluf:=[r,n]: fi: fi: od: od: si,aluf: end: #Smallestm(r,n,conf): The smallest m such that the prob. #that all boxes would have <=m balls in them with prob.>=conf #For example, try: #Smallestm(14400,9000,.99); Smallestm:=proc(r,n,conf) local m: for m from 1 while Prnm(r,n,m)<=conf do od: m: end: #SmallestmExaxt(r,n,conf): The smallest m such that the prob. #that all boxes would have <=m balls in them with prob.>=conf #For example, try: #SmallestmExact(14400,9000,.99); SmallestmExact:=proc(r,n,conf) local m: for m from 1 while PrnmExact(r,n,m)<=conf do od: m: end: #SmallestmPA(r,n,conf): The smallest m such that the prob. #that all boxes would have <=m balls in them with prob.>=conf #For example, try: #SmallestmPA(14400,9000,.99); SmallestmPA:=proc(r,n,conf) local m: for m from 1 while PrnmPA(r,n,m)<=conf do od: m: end: #AsyDBv(Maxa,Maxb,Maxm,M,K,n): tells a story and #creates a data base for the #asymptotic formulas, in n, of order M, for the probability of no box #getting more than m boxes when you place a*n balls in #b*n boxes for a and b relatively prime a<=Maxa, b<=Maxb #m<=Maxm. K is a large positive integer for estimation purposes. #The output is a list of pairs [[a,b,m],formula] #For example, try: #AsyDBv(2,2,4,4,5000,n); AsyDBv:=proc(Maxa,Maxb,Maxm,M,K,n) local a,b,yakhas,kv,m,lu,gu,a1,b1,co: co:=0: print(`Asymptotic Formulas for the Probability that No Box has more than`): print(`a given number of Balls`): print(``): print(`By Shalosh B. Ekhad `): print(``): print(`If you place a*n balls, uniformly at random, in b*n boxes`): print(`What is the probability that NO box received more than m balls?`): print(`In this web-book we will answer this question, for large n, and`): print(`(fixed) integers, a, b, m`): print(`by deriving ASYMPTOTIC FORMULAS to order`,M): print(`(perfectly rigorously, except for the constant in front,`): print(` but we omit the proofs),`): print(`for a and b relatively prime and`): print(`a between 1 and `, Maxa): print(`b between 1 and `, Maxb): print(`m between trunc(a/b)+1 and `, Maxm): kv:={}: gu:={}: for a from 1 to Maxa do for b from 1 to Maxb do yakhas:=a/b: a1:=numer(yakhas): b1:=denom(yakhas): if not member([a1,b1],kv) then kv:=kv union {[a1,b1]}: for m from trunc(a/b)+1 to Maxm do lu:=Asyabm(a1,b1,m,M,K,n): if lu<>FAIL then co:=co+1: gu:=[op(gu),[[a1,b1,m],lu]]: AsyabmV1(a1,b1,m,M,K,n,co): print(``): print(`------------------------------------------------------------`): print(``): fi: od: fi: od: od: print(`To summarize, the data base is`): lprint(gu): end: #DB227(n): The data base of all asymptotic formulas for #Asyabm(a,b,m,4,5000,n) for 1<=a,b<=2, 2<=m<=7 #For example, try: #DB227(n); DB227:=proc(n): [[[1, 1, 2], exp(-.118626412980456974767390675020*n)*(1.30656296487637652785891\ 855200-.102275552582151910126929076373/n-.129093989631668878163464807408e-1/n^2 +.518527231242236406208120583522e-2/n^3+.987187467122550170659871016510e-2/n^4) ], [[1, 1, 3], exp(-.215520483132030233649237984900e-1*n)*(1.089297330557798411\ 89394201670-.609771155454836427003291783474e-1/n-.\ 166803929805670701249851082582e-1/n^2-.856371311058873694851147110011e-2/n^3-.\ 566485093052986914550579046959e-2/n^4)], [[1, 1, 4], exp(-.\ 379213754697383021450919646200e-2*n)*(1.02558297302972311535159114898-.\ 382153991415876152544667757666e-1/n-.845389808366167671364344240166e-2/n^2-.\ 592985020694844715544346650005e-2/n^3-.109269443116002154787894832435e-1/n^4)], [[1, 1, 5], exp(-.599134814791224472197124666000e-3*n)*(1.006290019152518297853\ 42623937-.182512269841929640156941646556e-1/n+.526083375206458002814652873706e-\ 2/n^2+.113775751204163036111803894336e-1/n^3+.204167231429803924337916878644e-1 /n^4)], [[1, 1, 6], exp(-.833755469705394457811151510000e-4*n)*(1.0012838676267\ 2538586418357402-.634722307434075093298452306344e-2/n+.\ 851741762946014859363568075217e-2/n^2+.826120120023067383214683526217e-2/n^3+.\ 441148658943100564583125702368e-2/n^4)], [[1, 1, 7], exp(-.\ 102519144623715440570021380000e-4*n)*(1.00021916209812449681217264332-.\ 165911874459576518131741593032e-2/n+.475229699163670655413925905078e-2/n^2-.\ 931320437186505460202028695655e-3/n^3-.860850885256682260693965291460e-2/n^4)], [[1, 2, 2], exp(-.325864706845311860268546635380e-1*n)*(1.081951615066790453447\ 52282674-.275001911682903497108940606133e-1/n-.576134190347567591648821006914e-\ 2/n^2-.153467679804589823508707684201e-2/n^3+.396656301986066105281091653585e-3 /n^4)], [[1, 2, 3], exp(-.358982603542055603493713404800e-2*n)*(1.0169933799505\ 9206910141682025-.171035242182198533423389019604e-1/n-.\ 322270068160471359074433351899e-2/n^2-.150899118431601473701645196452e-2/n^3-.\ 140145508782939459320402026469e-2/n^4)], [[1, 2, 4], exp(-.\ 345517023903891872412247008000e-3*n)*(1.00279614876761989177820362127-.\ 635829784988417352906204295844e-2/n+.198975304200808384370231322641e-2/n^2+.\ 244427728172200728032239660753e-2/n^3+.274187608122772234401816972099e-2/n^4)], [[1, 2, 5], exp(-.283425600474604208444417080000e-4*n)*(1.000355893118482286965\ 02460731-.147387854965165978959554890496e-2/n+.193153843065887697277545139334e-\ 2/n^2+.610238255223039859107652234558e-3/n^3-.316649563178029465230030471994e-3 /n^4)], [[1, 2, 6], exp(-.200484684420283360708856800000e-5*n)*(1.0000362021330\ 8942387309942257-.237103852239994029500380597019e-3/n+.\ 641212449137374823474971859366e-3/n^2-.436208097130998975528012215847e-3/n^3-.\ 351592966319437884623416599646e-3/n^4)], [[1, 2, 7], exp(-.\ 124394263117810851978508000000e-6*n)*(1.00000305562949794557629185726-.\ 289263914328242805465200397702e-4/n+.128446005768584932712136554353e-3/n^2-.\ 241498919855944162098377290939e-3/n^3+.258852430902102600047110130855e-3/n^4)], [[2, 1, 2], (2*Pi^(1/2)/(1/n)^(1/2)+1/12*Pi^(1/2)*(1/n)^(1/2)+1/576*Pi^(1/2)*(1 /n)^(3/2)-139/207360*Pi^(1/2)*(1/n)^(5/2))*2.^(-1.*n)], [[2, 1, 3], exp(-.24682\ 4033659199894088503318400*n)*(1.52051625216732015792751537789-.1598115734953169\ 39033066851191/n-.188924800749149893344030565708e-1/n^2+.\ 638692935981624798976175590129e-3/n^3+.651615839425598827173916955817e-2/n^4)], [[2, 1, 4], exp(-.670057779617542672546094274100e-1*n)*(1.186304296087581525011\ 09061902-.831939290253426848992184323406e-1/n-.134527781947433090001860979724e-\ 1/n^2-.210066571714976867601923193491e-2/n^3+.183475553118729417542515683240e-2 /n^4)], [[2, 1, 5], exp(-.182739505794619657907309184600e-1*n)*(1.0725981785501\ 6164757528020215-.565946453079145026973998361801e-1/n-.\ 109576103789346899473625939016e-1/n^2-.600247880844895203045943853005e-2/n^3-.\ 666194203675480057923935110833e-2/n^4)], [[2, 1, 6], exp(-.\ 469908623776670954665815100000e-2*n)*(1.02681201553604489176876982220-.\ 365146732216364540883218339048e-1/n-.357921919041620797611649076366e-2/n^2-.\ 100918458452456658515547093071e-2/n^3-.236745347412175579607986198368e-2/n^4)], [[2, 1, 7], exp(-.110939746243074916158884632000e-2*n)*(1.008902722003517003509\ 69278329-.200597214853875174704267578953e-1/n+.474004350086863654051141174786e-\ 2/n^2+.764793038265199202491636565046e-2/n^3+.110500955748151575370781011020e-1 /n^4)]]: end: #AsyabmPC(a,b,m,n): The prec-computed asymptotic formula, to order 4 #(using 5000 terms of the sequence) for Prnm(a*n,b*n,m); #for 1<=a,b<=2 and 2<=m<=7 #For example, try: #AsyabmPC(1,1,4,n); AsyabmPC:=proc(a,b,m,n) local gu,i: gu:=DB227(n): if not (1<=a and a<=2 and 1<=b and b<=2 and 2<=m and m<=7) then print([a,b,m], `has not been pre-computed you must use Asyabm`): RETURN(FAIL): fi: for i from 1 to nops(gu) do if gu[i][1]=[a,b,m] then RETURN(gu[i][2]): fi: od: FAIL: end: ###Begin generalized program, between m1 and m2 (inclusive) #PrnmG(r,n,m1,m2): The probability (in floating point) #that if you put r balls in n boxes #each box received between m1 and m2 balls (inclusive) #For example, try: #PrnmG(14400,9000,2,7); PrnmG:=proc(r,n,m1,m2) local f,i,x: if r>n*m2 then RETURN(0): fi: if r=m2 then print(`The third argument must be less than the fourth argument`): RETURN(FAIL): fi: if m1<0 or m2<=2 then RETURN(FAIL): fi: ope:=RecabmG(a,b,m1,m2,n,N): Ini:=[seq(ArnmG(a*i,b*i,m1,m2),i=1..degree(ope,N))]: lu:=AsyC(ope,n,N,M,Ini,K): if lu=FAIL then RETURN(FAIL): fi: C:=lu[1]: lu:=evalf(lu[2]): gu:=lu: ku:=1: for i from 1 to nops(lu) do if type(op(i,lu),`^`) and type(op(1,op(i,lu)),float) then ku:=ku*op(i,lu): gu:=gu/op(i,lu): fi: od: ku:=simplify(ku): ru:=asympt(exp(n)*n!/n^n,n,M): ru:=ru-op(nops(ru),ru): ru:=subs(n=n*a,ru): if type(gu,`^`) then RETURN(ru*gu): fi: gu:=expand(ru*gu): if not type(gu,`+`) then print(`Something is wrong`): RETURN(FAIL): fi: gu:=add(C*coeff(gu,n,-i)/n^i,i=0..M): ku:=simplify(ku*exp(-n*a)*((a/b)^a)^n): evalf(ku*gu): end: #AsyabmGV(a,b,m1,m2,M,K,n): Verbose version of #AsyabmG(a,b,m1,m2,M,K,n) and with a theorem number #For example do #AsyabmGV(1,1,1,4,5,5000,n); AsyabmGV:=proc(a,b,m1,m2,M,K,n) local lu,mu,ku: lu:=AsyabmG(a,b,m1,m2,M,K,n): print(`Theorem: Let P(n) be the probability that if you `): print(`place (uniformly)`, a*n,`balls in `, b*n , `boxes then `): print(`no box has more than `, m, `balls. `): print(`Then we have the following asymptotic formula to order `, M): print(P(n)=lu): print(``): print(`Comment 1.: compare this EXACT result with the naive`): print(` Poisson approximation`): mu:=AsyabmPA(a,b,m,n): print(mu): print(`Comment 2: Let's see how good the asymptotic formula is`): print(`For n=100, it gives you`): print(evalf(subs(n=100,lu))): print(`whereas the exact value, using the Ewens-Wilf approach is`): ku:=evalf(PrnmExact(100*a,100*b,m)): print(ku): print(`The ratio is:`, evalf(subs(n=100,lu)/ku)): print(`For n=200, it gives you`): print(evalf(subs(n=200,lu))): print(`whereas the exact value, using the Ewens-Wilf approach is`): ku:=evalf(PrnmExact(1000*a,1000*b,m)): print(ku): print(`The ratio is:`, evalf(subs(n=1000,lu)/ku)): print(``): print(`In Maple input format, the asymptotic formula is`): lprint(lu): end: #RecabmGV(a,b,m1,m2,n,N): Verbose version of RecabmG(a,b,m1,m2,n,N) #For example, try: #RecabmGV(2,1,1,3,n,N): RecabmGV:=proc(a,b,m1,m2,n,N) local ope,Ini,P,Q,i: print(`Theorem: Let P(n) be the probability that if you place (uniformly)`): print(a*n,`balls in `, b*n , `boxes then each box contains between `): print(m1, `balls and `, m2, `balls (inclusive). `) : print(`Then we have: `): print(P(n)=(a*n)!/(b*n)^(a*n)*Q(n)): print(``): ope:=RecabmG(a,b,m1,m2,n,N): if ope=FAIL then RETURN(FAIL): fi: Ini:=[seq(ArnmG(a*i,b*i,m1,m2),i=1..degree(ope,N))]: print(`where Q(n) is the sequence satisfying the recurrece`): print(add(coeff(ope,N,i)*Q(n+i),i=0..degree(ope,N))=0): print(``): print(`subject to the initial conditions`): print(seq(Q(i)=Ini[i],i=1..nops(Ini))): print(``): print(`In operator notation, using `, N, `as the shift operator in `, n): print(`the recurrence operator annihilating Q(n) is`): lprint(ope): end: #AsyabmGV(a,b,m1,m2,M,K,n): Verbose version of #AsyabmGV(a,b,m1,m2,M,K,n): #For example do #AsyabmGV(1,1,1,4,5,5000,n); AsyabmGV:=proc(a,b,m1,m2,M,K,n) local lu,mu,ku: lu:=AsyabmG(a,b,m1,m2,M,K,n): print(`Theorem: Let P(n) be the probability that if you place (uniformly)`): print(a*n,`balls in `, b*n , `boxes each box has at least`, m1 , `balls `): print(`and at most `, m2, `balls `): print(`Then we have the following asymptotic formula to order `, M): print(P(n)=lu): print(``): print(`Comment 1.: compare this EXACT result with the naive`): print(` Poisson approximation`): mu:=AsyabmGPA(a,b,m1,m2,n): print(mu): print(`Comment 2: Let's see how good the asymptotic formula is`): print(`For n=100, it gives you`): print(evalf(subs(n=100,lu))): print(`whereas the exact value, using the Ewens-Wilf approach is`): ku:=evalf(PrnmG(100*a,100*b,m1,m2)): print(ku): print(`The ratio is:`, evalf(subs(n=100,lu)/ku)): print(`For n=1000, it gives you`): print(evalf(subs(n=1000,lu))): print(`whereas the exact value, using the Ewens-Wilf approach is`): ku:=evalf(PrnmG(1000*a,1000*b,m1,m2)): print(ku): print(`The ratio is:`, evalf(subs(n=1000,lu)/ku)): print(``): print(`In Maple input format, the asymptotic formula is`): lprint(lu): end: ###End generalized program, between m1 and m2 (inclusive) #LargestmPA(r,n,conf): The largest m such that the prob. #that all boxes would have >=m balls in them with prob.>=conf #For example, try: #LargestmPA(1000,100,.99); LargestmPA:=proc(r,n,conf) local m,M1: for m from 1 while PrnmPA(r,n,m)<=0.99999 do od: M1:=m: for m from 1 while PrnmGPA(r,n,m,M1)>conf do od: m: end: #Largestm(r,n,conf): The largest m such that the prob. #that all boxes would have >=m balls in them with prob.>=conf #For example, try: #Largestm(1000,100,.99); Largestm:=proc(r,n,conf) local m,M1: for m from 1 while evalf(PrnmGexact(r,n,0,m))<=0.99999 do od: M1:=m: for m from 1 while PrnmGexact(r,n,m,M1)>conf do od: m: end: #Largestm(r,n,conf): The largest m such that the prob. #that all boxes would have >=m balls in them with prob.>=conf #For example, try: #Largestm(1000,100,.99); Largestm:=proc(r,n,conf) local m,M1: for m from 1 while evalf(PrnmGexact(r,n,0,m))<=0.99999 do od: M1:=m: for m from 1 while PrnmGexact(r,n,m,M1)>conf do od: m: end: ###Start Sequences #RecabmSeq(a,b,m,n,N,L): The linear recurrence (with poly. coeffs.) #satistfied by #(a*n)! times the coeffs. of x^(a*n) in (1+x+x^2/2!+...+x^m/m!)^(b*n) #where a and b are pos. integers (a*n balls in b*n boxes). #also the initial values and the first L terms #For example, try: #RecabmSeq(2,1,3,n,N,30): RecabmSeq:=proc(a,b,m,n,N,L) local i,x,Ini,ope,Ini1: option remember: ope:=AZd((a*n)!*add(x^i/i!,i=0..m)^(n*b)/x^(n*a+1),x,n,N)[1]: Ini:=[seq((a*i)!*Arnm(a*i,b*i,m),i=1..degree(ope,N))]: Ini1:=[seq((a*i)!*Arnm(a*i,b*i,m),i=1..20)]: if SeqFromRec(ope,n,N,Ini,20)<>Ini1 then print(`Something terrible happened`): RETURN(FAIL): fi: Ini1:=[seq((a*i)!*Arnm(a*i,b*i,m),i=1..L)]: ope,Ini,Ini1: end: #RecabmSeqV(a,b,m,n,N,L): Verbose version of RecabmSeq(a,b,m,n,N,L) #For example, try: #RecabmSeqV(2,1,3,n,N,30): RecabmSeqV:=proc(a,b,m,n,N,L) local ope,Ini,P,Q,i,gu,Ini1: print(`Theorem: Let A(n) be the number of ways of placing`): print(a*n,`balls in `,b*n,`boxes such that no box has more than `,m,`balls. `): gu:=RecabmSeq(a,b,m,n,N,L): if gu=FAIL then RETURN(FAIL): fi: ope:=gu[1]: Ini:=gu[2]: Ini1:=gu[3]: print(`then A(n) is the sequence satisfying the recurrece`): print(add(coeff(ope,N,i)*A(n+i),i=0..degree(ope,N))=0): print(``): print(`subject to the initial conditions`): print(seq(A(i)=Ini[i],i=1..nops(Ini))): print(``): print(`The first (starting at n=1)`,L,`terms of the enumerating sequence are`): print(`Ini1`): print(`In operator notation, using `, N, `as the shift operator in `, n): print(`the recurrence operator annihilating Q(n) is`): lprint(ope): end: #RecabmSeqV1(a,b,m,n,N,L): Like RecabmSeqV1(a,b,m,n,N,L,co): #but with a numbered theorem. #For example, try: #RecabmSeqV1(2,1,3,n,N,30,2): RecabmSeqV1:=proc(a,b,m,n,N,L,co) local ope,Ini,P,Q,i,gu,Ini1: printf("Theorem Number %d :\n" , co): print(` Let A(n) be the number of ways of placing`): print(a*n,`balls in `,b*n,`boxes such that no box has more than `,m,`balls. `): gu:=RecabmSeq(a,b,m,n,N,L): if gu=FAIL then RETURN(FAIL): fi: ope:=gu[1]: Ini:=gu[2]: Ini1:=gu[3]: print(`then A(n) is the sequence satisfying the recurrece`): print(add(coeff(ope,N,i)*A(n+i),i=0..degree(ope,N))=0): print(``): print(`subject to the initial conditions`): print(seq(A(i)=Ini[i],i=1..nops(Ini))): print(``): print(`The first (starting at n=1)`,L,`terms of the enumerating sequence are`): print(`Ini1`): print(`In operator notation, using `, N, `as the shift operator in `, n): print(`the recurrence operator annihilating Q(n) is`): lprint(ope): end: #SeferRecs(Maxa,Maxb,Maxm,L,n,N): Inputs positive integers #Maxa, Maxb and Maxm, as well as symbols n and N (shift in n) #and outputs a web-book with the recurrences, and the first #L terms for the sequences enumerating the number of ways #of placing a*n balls in b*n boxes such that no box has more #than m balls, for a and b relatively prime #a<=Maxa, b<=Maxb, m<=Maxm. #For example, try: #SeferRecs(2,2,4,30,n,N); SeferRecs:=proc(Maxa,Maxb,Maxm,L,n,N) local a,b,yakhas,kv,m,lu,gu,a1,b1,co: co:=0: print(`Recurrences for the Number of Ways of Placing Balls in Boxes`): print(` Such that No Box has more than a given number of Balls`): print(``): print(`By Shalosh B. Ekhad `): print(``): print(`If you place a*n balls in b*n boxes`): print(`In how many ways can you do it in such a way that NO box received`): print(` more than m balls?`): print(`In this web-book we will ANSWER this question `): print(`(in the sense of Herb Wilf) by exhibiting linear recurrences, in `, n): print(`with polynomial coefficients`): print(`for (fixed) integers, a, b, m`): print(`for a and b relatively prime and`): print(`a between 1 and `, Maxa): print(`b between 1 and `, Maxb): print(`m between a/b and `, Maxm): print(`We also give, for Sloane's sake, the first `, L, `terms of each sequence `): print(` Everything is perfectly rigorous `): print(` but we omit the proofs.`): kv:={}: gu:={}: for a from 1 to Maxa do for b from 1 to Maxb do yakhas:=a/b: a1:=numer(yakhas): b1:=denom(yakhas): if not member([a1,b1],kv) then kv:=kv union {[a1,b1]}: for m from trunc(a/b)+1 to Maxm do lu:=RecabmSeq(a1,b1,m,n,N,L): if lu<>FAIL then co:=co+1: gu:=[op(gu),[[a1,b1,m],[lu[1],lu[2]]]]: RecabmSeqV1(a1,b1,m,n,N,L,co): print(``): print(`------------------------------------------------------------`): print(``): fi: od: fi: od: od: print(`To summarize, the data base is`): lprint(gu): end: ###End Sequences #AsyabmEmpir(a,b,m,K,n): The estimated asymptotic formula #(using K terms of the sequence) for Prnm(a*n,b*n,m); #For example, try: #AsyabmEempir(1,1,4,100,n); AsyabmEmpir:=proc(a,b,m,K,n) local gu1,gu2,mu,C: option remember: if m<=a/b then RETURN(0): fi: gu1:=PrnmExact(a*(K-1),b*(K-1),m): gu2:=PrnmExact(a*K,b*K,m): mu:=evalf(gu2/gu1): C:=evalf(gu2/mu^K): mu:=log(mu): C*exp(mu*n): end: #PrnmEmpir(r,n,m,K,CU): Like PrnmEmpir(r,n,m) #but using extrapolation of the empirically #derived Asyabm(a,b,m,K): CU is the cut-off #of whether to use r/n=a/b or use an approximation #For example, try: #PrnmEmpir(14400,9000,7,100); PrnmEmpir:=proc(r,n,m,K,CU) local a,b,y,n1,lu,N,cf,ka,i,y1: if r>n*m then RETURN(0): fi: y:=r/n: a:=numer(y): b:=denom(y): if b>CU then cf:=convert(y,confrac,ka): for i from 1 to nops(ka) while denom(ka[i])10^(-10) do od: m1:=m: gu:=[seq(AsyabmEmpir(a,b,m,K,n),m=1..m1)]: gu:=[0,seq(gu[i]-gu[i-1],i=2..nops(gu))]: gu:=add(i*gu[i],i=1..nops(gu)): mu:=[seq([evalf(log(2^i*L)),evalf(subs(n=2^i*L,gu))],i=0..k)]: gu:=CurveFitting[LeastSquares](mu,u,curve=add(a1[i]*u^i,i=0..d)): subs(u=log(n),gu): end: #PdistrEmpirS(a,b,n,M1,K1,M2): The first M1 members of the #as (asymptotic) formulas in n of the prob. distribution #Max number of Balls, using AsyabmEmpir(a,b,m,K1,n) #followed by the average, s.d. and the normalized moments #through the M2 #For example, try: #PdistrEmpirS(2,1,n,10,100,10); PdistrEmpirS:=proc(a,b,n,M1,K1,M2) local m,gu,i,av,m2,j: option remember: gu:=[seq(AsyabmEmpir(a,b,m,K1,n),m=1..M1)]: gu:=[gu[1],seq(gu[i]-gu[i-1],i=2..nops(gu))]: av:=expand(add(i*gu[i],i=1..nops(gu))): m2:=expand(add(gu[i]*(i-av)^2,i=1..nops(gu))): gu, [av,sqrt(m2),seq(add(gu[i]*(i-av)^j,i=1..nops(gu))/m2^(j/2),j=3..M2)]: end: #PdistrC(r,n,eps): The central part where`): #most (1-eps) of`): #the probability distribution Pdistr(r,n,eps,K) is located`): #followed by the average, variance, and centralized-normalized`): #distribution. `): #For example, try:`): #PdistrC(14400,9000,10^(-3));`): PdistrC:=proc(r,n,eps,K) local gu,ave1,sd1,i1,i2,i,mu1,mu2,ku: gu:=Pdistr(r,n,10^(-10),K); ku:=gu[2]: ave1:=gu[2][1]: sd1:=gu[2][2]: gu:=gu[1]: for i from 1 while convert([op(1..i,gu)],`+`)2 then print(`The skewness of M is`, gu[2][3]): fi: if K>3 then print(`The kurtosis of M is`, gu[2][4]): fi: for i from 5 to K do print(`The `, i , `-th alpha coefficient is`, gu[2][i]): od: print(`----------------------------------------`): print(`Ignoring tails, the normalized distribution is`): print(gu[3]): end: #PdistrPACv(r,n,eps,K): Verbose version of PdistrPAC #For example: try: #PdistrPACv(1000,2000,10^(-10),6): PdistrPACv:=proc(r,n,eps,K) local gu,i: gu:=PdistrPAC(r,n,eps,K): print(`Suppose that`, r, ` balls are thrown, uniformly at random, `): print(`in `, n, `boxes, and consider the random variable, let's call it M, `): print(`"Maximum number of balls in the same box" `): print(`Using the (very good!) Poisson approximation we have the following`): print(`estimates.`): print(`The expected value is`, gu[2][1]): print(`The standard deviation is`, gu[2][2]): print(``): print(`Ignoring tails on either end whose probability is less than`,eps): print(`The prob. distribution is as followed `): for i from 1 to nops(gu[1]) do print(`The probability that M equals`, gu[1][i][1], `equals `, gu[1][i][2]): od: print(`----------------------------`): if K>2 then print(`The skewness of M is`, gu[2][3]): fi: if K>3 then print(`The kurtosis of M is`, gu[2][4]): fi: for i from 5 to K do print(`The `, i , `-th alpha coefficient is`, gu[2][i]): od: print(`----------------------------------------`): print(`Ignoring tails, the normalized distribution is`): print(gu[3]): end: #PdistrEmpirCv(r,n,eps,K,K1): Verbose version of PdistrEmpirC #For example: try: #PdistrEmpirv(1000,2000,10^(-6),100,6): PdistrEmpirCv:=proc(r,n,eps,K,K1) local gu,i: gu:=PdistrEmpirC(r,n,eps,K,K1): print(`Suppose that`, r, ` balls are thrown, uniformly at random, `): print(`in `, n, `boxes, and consider the random variable, let's call it M, `): print(`"Maximum number of balls in the same box" `): print(`Using the (excellent!) Empirical asymptotic Approximation`): print(`employing `, K1, `terms for estimating the asymptotics,`): print(` we have the following estimates.`): print(`The expected value is`, gu[2][1]): print(`The standard deviation is`, gu[2][2]): print(``): print(`Ignoring tails on either end whose probability is less than`,eps): print(`The prob. distribution is as followed `): for i from 1 to nops(gu[1]) do print(`The probability that M equals`, gu[1][i][1], `equals `, gu[1][i][2]): od: print(`----------------------------`): if K>2 then print(`The skewness of M is`, gu[2][3]): fi: if K>3 then print(`The kurtosis of M is`, gu[2][4]): fi: for i from 5 to K do print(`The `, i , `-th alpha coefficient is`, gu[2][i]): od: print(`----------------------------------------`): print(`Ignoring tails, the normalized distribution is`): print(gu[3]): end: #AveFormulaStory(A): All the average formulas for the #random variable "max. number of balls" #Rn balls n boxes 1<=R<=A and n balls in Rn boxes #for example try: #AveFormulaStory(5,n,4,100,1000,10); AveFormulaStory:=proc(A,n,d,K,L,k) local gu,R,lu,lu1,mu: for R from 1 to A do print(`-------------------------------------`): gu:=AveFormula(R,1,n,d,K,L,k): print(`When `, R*n, `balls are thrown into`, n, ` boxes `): print(`The expected value of the random variable : max number of balls in`): print(` a box is roughly`): print(gu): print(`and in Maple input form:`): lprint(gu): print(`Let's see how good it is`): print(`----------------------`): lu:=Pdistr(R*100,100,10^(-10),2)[2][1]: print(`for n=100, the true value of the average is`, lu): lu1:=evalf(subs(n=100,gu)): print(`While our appx. gives:`, lu1): mu:=lu1/lu: print(`Their ratio is`, mu): if abs(mu-1)>0.3 then print(`Oops, there must be an overflow error, in the computation of the`): print(`real thing!. The value of`, lu, `for the average does not make sense`): print(`Let's take instead the Poisson approximation`): lu:=PdistrPA(R*100,100,10^(-10),2)[2][1]: print(lu): mu:=lu1/lu: print(`Now their ratio is`, mu): fi: print(`----------------------`): lu:=Pdistr(R*200,200,10^(-10),2)[2][1]: print(`for n=200, the true value of the average is`, lu): lu1:=evalf(subs(n=200,gu)): print(`While our appx. gives:`, lu1): mu:=lu1/lu: print(`Their ratio is`, mu): if abs(mu-1)>0.3 then print(`Oops, there must be an overflow error, in the computation of the`): print(`real thing!. The value of`, lu, `for the average does not make sense`): print(`Let's take instead the Poisson approximation`): lu:=PdistrPA(R*200,200,10^(-10),2)[2][1]: print(lu): mu:=lu1/lu: print(`Now their ratio is`, mu): fi: od: for R from 2 to A do print(`-------------------------------------`): gu:=AveFormula(1,R,n,d,K,L,k): print(`When `, n, `balls are thrown into`, R*n, ` boxes `): print(`The expected value of the random variable : max number of balls in`): print(` a box is roughly`): print(gu): print(`and in Maple input form:`): lprint(gu): print(`Let's see how good it is`): print(`----------------------`): lu:=Pdistr(1000,R*1000,10^(-10),2)[2][1]: print(`for n=1000, the true value of the average is computed as`, lu): lu1:=evalf(subs(n=100,gu)): print(`While our appx. gives:`, lu1): mu:=lu1/lu: print(`Their ratio is`, mu): if abs(mu-1)>0.3 then print(`Oops, there must be an overflow error, in the computation of the`): print(`real thing!. The value of`, lu, `for the average does not make sense`): print(`Let's take instead the Poisson approximation`): lu:=PdistrPA(100,R*100,10^(-10),2)[2][1]: print(lu): mu:=lu1/lu: print(`Now their ratio is`, mu): fi: print(`----------------------`): od: end: #PrnmReliable(r,n,m,eps): computes Prnm(r,n,m) #with floating-point reliable to k digits #For example, try: #PrnmReliable(14400,9000,8,20); PrnmReliable:=proc(r,n,m,k) local lu1,lu2,lu3,i: lu1:=evalf(Prnm(r,n,m),100+k): lu2:=evalf(Prnm(r,n,m),200+k): if abs(lu1-lu2)<10^(-k-4) then RETURN(evalf(lu2,k)): fi: for i from 3 while abs(lu1-lu2)>10^(-k-4) do lu3:=evalf(Prnm(r,n,m),100*i+k): lu1:=lu2: lu2:=lu3: od: RETURN(evalf(lu2,k)): end: #PrnmExactStory(r,n,m1,m2,K): The story of #the exact evaluation of Prnm(r,n,m) with timings #and the number of digits. For example, try: #PrnmExactStory(14400,9000,6,12,40); PrnmExactStory:=proc(r,n,m1,m2,K) local m,lu,t0,t1: print(`This is the story, and timings, of the EXACT`): print(`values of the probability that when `, r, `balls `): print(`are thrown into`, n, `boxes `): print(`each box received at most m balls for m from`, m1, `to `, m2): print(`because of the size of the rational numbers involved, we do not`): print(`dispaly them but rather give the floating-point version`): print(`to `, K, `digits `): print(`Note that these numbers are EXACT, not numerical estimates`): for m from m1 to m2 do t0:=time(): lu:=PrnmExact(r,n,m): t1:=time(): print(`When m=`, m, `the probability is a certain rational number `): print(`that has `, trunc(evalf(log[10](numer(lu))))): print(`(decimal) digits in its numerator `): print(`and `, trunc(evalf(log[10](denom(lu)))), `in its denominator `): print(`Its floating-point representation, to`, K, `digits is`): print(evalf(lu,K)): print(`The exact computation of this rational number took`,t1-t0, `seconds. `): od: end: