#######################################################################
## Research: Save this file as ScoringPaths. To use it, stay in the   #
## same directory, get into Maple (by typing: maple <Enter> )         #
## and then type:  read "ScoringPaths.txt"; <Enter>                   #
## Then follow the instructions given there                           #
##                                                                    #
## Written by Bryan Ek, Rutgers University,                           #
##  bryan.t.ek@math.rutgers.edu.                                      # 
#######################################################################

print(`This is ScoringPaths`):
print(`by Bryan Ek, advisee to Dr. Zeilberger, Rutgers University.`):
print(``):
print(`If you have comments or notice any errors, please email`):
print(`bryan [dot] t [dot] ek [at] math [dot] rutgers [dot] edu`):
print(`bryan.t.ek@math.rutgers.edu`):
print(``):
print(`The latest version (4/11/2018) of this package is available from`):
print(`http://www.math.rutgers.edu/~bte14/Code/ScoringPaths/ScoringPaths.txt`):
print(``):
print(`This package contains functions for enumerating walks.`):
print(`Given a set of possible steps on a square lattice, how many different walks are there of length n starting from the origin?`):
print(`What if we restrict to having cyclic walks? I.e. we return the x-axis.`):
print(`What if we disallow walks going below the line y=b<0? Or above the line y=a>0?`):
print(`This package uses generating function relations to answer these questions.`):
print(``):
print(`For a detailed list of functions and their purpose type Help().`):
print(`For functions borrowed from other packages, type ezraSCHUTZENBERGER(), ezraEKHAD(), ezraW1D(), and ezraAsyRec.`):
print(``):

###########################################################################################################################################
#Help function
Help:=proc():
	if args=`BoundedScoringPathsEqsVars` then
		print(`BoundedScoringPathsEqsVars(possibleScores,upperLimit,lowerLimit,t,f): inputs possible scores and the upperLimit and lowerLimit for the y-value.`):
		print(`possibleScores is a list or set of scores in the format [x-step,y-step].`):
		print(`Assumes that upperLimit>=0,lowerLimit<=0.`):
		print(`f is a placeholder to keep track of the different g.f.s.`):
		print(`t is the variable in the g.f.`):
		print(`Outputs a set of equations (and the used variables) relating the g.f.s for the number of paths of total x-step n, that start at (0,0) and stay within the limits.`):
		print(`Includes the g.f.s for all pairs of limits (up,low) such that up-low=upperLimit-lowerLimit.`):
		print(`All equations are linear.`):
		print(`Try BoundedScoringPathsEqsVars({[1,0],[1,-2],[2,-3],[0,1]},3,-2,t,f).`):
	elif args=`AllBoundedScoringPaths` then
		print(`AllBoundedScoringPaths(possibleScores,upperLimit,lowerLimit,t,f): inputs possible scores and the upperLimit and lowerLimit for the y-value.`):
		print(`possibleScores is a list or set of scores in the format [x-step,y-step].`):
		print(`Assumes that upperLimit>=0,lowerLimit<=0.`):
		print(`f is a placeholder to keep track of the different g.f.s.`):
		print(`t is the variable in the g.f.`):
		print(`Outputs ALL the g.f.s for the number of paths of total x-step n, that start at (0,0) and stay within its limits.`):
		print(`Includes the g.f.s for all pairs of limits (up,low) such that up-low=upperLimit-lowerLimit.`):
		print(`All equations are linear so the generating functions are all rational!`):
		print(`Uses BoundedScoringPathsEqsVars.`):
		print(`Try AllBoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,1]},3,-2,t,f).`):
	elif args=`BoundedScoringPaths` then
		print(`BoundedScoringPaths(possibleScores,upperLimit,lowerLimit,t): inputs possible scores and the upperLimit and lowerLimit for the y-value.`):
		print(`possibleScores is a list or set of scores in the format [x-step,y-step].`):
		print(`Assumes that upperLimit>=0,lowerLimit<=0.`):
		print(`t is the variable in the g.f.`):
		print(`Outputs only the g.f. for the number of paths of total x-step n, that start at (0,0) and stay within the limits given.`):
		print(`Uses AllBoundedScoringPaths.`):
		print(`Try BoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,1]},3,-2,t).`):
	elif args=`BoundedScoringPathsCoefficients` then
		print(`BoundedScoringPathsCoefficients(possibleScores,upperLimit,lowerLimit,N): inputs possible scores and the upperLimit and lowerLimit for the y-value.`):
		print(`possibleScores is a list or set of scores in the format [x-step,y-step].`):
		print(`Assumes that upperLimit>=0,lowerLimit<=0.`):
		print(`Outputs the first N+1 coefficients of the Taylor series of the g.f. for the number of paths that start at (0,0) and stay within the limits.`):
		print(`Uses BoundedScoringPaths.`):
		print(`Try BoundedScoringPathsCoefficients({[1,0],[1,-2],[2,-3],[0,1]},3,-2,9).`):
	elif args=`BoundedScoringPathsNumber` then
		print(`BoundedScoringPathsNumber(possibleScores,upperLimit,lowerLimit,n): inputs possible scores and the upperLimit and lowerLimit for the y-value.`):
		print(`Assumes that upperLimit>=0,lowerLimit<=0.`):
		print(`possibleScores is a list or set of scores in the format [x-step,y-step].`):
		print(`Outputs the number of paths of x-step length n that start at (0,0) and stay within the bounds.`):
		print(`This function operates empirically.`):
		print(`If there are BOTH [0,-] and [0,+] steps that sum to 0 within the bounds, then you will obtain infinite recursion.`):
		print(`Try seq(BoundedScoringPathsNumber({[1,0],[1,-2],[2,-3],[0,1]},3,-2,n1),n1=0..9).`):

	elif args=`EqualBoundedScoringPathsEqsVars` then
		print(`EqualBoundedScoringPathsEqsVars(possibleScores,upperLimit,lowerLimit,t,f,e): inputs possible scores and the upperLimit and lowerLimit for the y-value.`):
		print(`Assumes that upperLimit>=0,lowerLimit<=0.`):
		print(`possibleScores is a list or set of scores in the format [x-step,y-step].`):
		print(`f,e are placeholder variables for the g.f.s.`):
		print(`f[upperLimit,lowerLimit] is the g.f. for the number of paths of total x-step n, that start at (0,0), end at (n,0), and stay within the limits.`):
		print(`e[upperLimit,lowerLimit,c] is the g.f. for the number of "irreducible" walks that start at (c,0), end at (n,0), stay within the given limits, and touch 0 only at the end.`):
		print(`t is the variable in the g.f.`):
		print(`Outputs a set of equations (and the used variables) relating the g.f.s held in f and e.`):
		print(`Try EqualBoundedScoringPathsEqsVars({[1,0],[1,-2],[2,-3],[0,1]},3,-2,t,f,e).`):
	elif args=`AllEqualBoundedScoringPaths` then
		print(`AllEqualBoundedScoringPaths(possibleScores,upperLimit,lowerLimit,t,f,e): inputs possible scores and the upperLimit and lowerLimit for the y-value.`):
		print(`Assumes that upperLimit>=0,lowerLimit<=0.`):
		print(`possibleScores is a list or set of scores in the format [x-step,y-step].`):
		print(`f,e are placeholder variables for the g.f.s.`):
		print(`f[upperLimit,lowerLimit] is the g.f. for the number of paths of total x-step n, that start at (0,0), end at (n,0), and stay within the limits.`):
		print(`e[upperLimit,lowerLimit,c] is the g.f. for the number of "irreducible" walks that start at (c,0), end at (n,0), stay within the given limits, and touch 0 only at the end.`):
		print(`t is the variable in the g.f.s.`):
		print(`Outputs the g.f.s f[upperLimit,lowerLimit] and e[upperLimit,lowerLimit,c] for all c=lowerLimit..-1,1..upperLimit.`):
		print(`Uses EqualBoundedScoringPathsEqsVars.`):
		print(`Try AllEqualBoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,1]},3,-2,t,f,e).`):
	elif args=`EqualBoundedScoringPaths` then
		print(`EqualBoundedScoringPaths(possibleScores,upperLimit,lowerLimit,t): inputs possible scores and the upperLimit and lowerLimit for the y-value.`):
		print(`Assumes that upperLimit>=0,lowerLimit<=0.`):
		print(`possibleScores is a list or set of scores in the format [x-step,y-step].`):
		print(`t is the variable in the g.f.`):
		print(`Outputs the g.f. for the number of paths of total x-step n, that start at (0,0), end at (n,0), and stay within the limits.`):
		print(`Uses EqualBoundedScoringPathsEqsVars.`):
		print(`Try EqualBoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,1]},3,-2,t).`):
	elif args=`EqualBoundedScoringPathsCoefficients` then
		print(`EqualBoundedScoringPathsCoefficients(possibleScores,upperLimit,lowerLimit,N): inputs possible scores and the upperLimit and lowerLimit for the y-value.`):
		print(`possibleScores is a list or set of scores in the format [x-step,y-step].`):
		print(`Assumes that upperLimit>=0,lowerLimit<=0.`):
		print(`Outputs the first N+1 coefficients of the Taylor series of the g.f. for the number of paths that start at (0,0), end at (n,0), and stay within the limits.`):
		print(`Uses EqualBoundedScoringPaths.`):
		print(`Try EqualBoundedScoringPathsCoefficients({[1,0],[1,-2],[2,-3],[0,1]},3,-2,9).`):
	elif args=`SpecificEqualBoundedScoringPathsNumber` then
		print(`SpecificEqualBoundedScoringPathsNumber(possibleScores,upperLimit,lowerLimit,n,startingLevel): inputs possible scores (x-step,y-step) and the upperLimit and lowerLimit for the y-value.`):
		print(`Assumes that upperLimit>=0,lowerLimit<=0, and startingLevel is somewhere between these 2 numbers.`):
		print(`possibleScores is a list or set of scores in the format [x,y].`):
		print(`Outputs the number of paths that start at (0,startingLevel), end at (n,0), and stay within the bounds.`):
		print(`Does this empirically.`):
		print(`If there are BOTH [0,-] and [0,+] steps that sum to 0 within the bounds, then you will obtain infinite recursion.`):
		print(`Try seq(SpecificEqualBoundedScoringPathsNumber({[1,0],[1,-2],[2,-3],[0,1]},3,-2,n1,0),n1=0..9).`):
	elif args=`SpecificIrreducibleBoundedScoringPathsNumber` then
		print(`SpecificIrreducibleBoundedScoringPathsNumber(possibleScores,upperLimit,lowerLimit,n,startingLevel): inputs possible scores (x-step,y-step) and the upperLimit and lowerLimit for the y-value.`):
		print(`Assumes that upperLimit>=0,lowerLimit<=0, and startingLevel is somewhere between these 2 numbers.`):
		print(`possibleScores is a list or set of scores in the format [x,y].`):
		print(`Outputs the number of paths that start at (0,startingLevel), end at (n,0) while never touching the x-axis previously, and stay within the bounds.`):
		print(`Does this empirically.`):
		print(`If there are BOTH [0,-] and [0,+] steps that sum to 0 within the bounds, then you will obtain infinite recursion.`):
		print(`Try seq(SpecificIrreducibleBoundedScoringPathsNumber({[1,0],[1,-2],[2,-3],[0,1]},3,-2,n1,0),n1=0..9).`):
	elif args=`SpecificIrreducibleBoundedScoringPathsNumberHelp` then
		print(`SpecificIrreducibleBoundedScoringPathsNumberHelp(possibleScores,upperLimit,lowerLimit,n,startingLevel): inputs possible scores (x-step,y-step) and the upperLimit and lowerLimit for the y-value.`):
		print(`Assumes that upperLimit>=0,lowerLimit<=0, and startingLevel is somewhere between these 2 numbers.`):
		print(`possibleScores is a list or set of scores in the format [x,y].`):
		print(`Outputs the number of paths that start at (0,startingLevel), end at (n,0) while never touching the x-axis previously, and stay within the bounds.`):
		print(`Does this empirically.`):
		print(`If there are BOTH [0,-] and [0,+] steps that sum to 0 within the bounds, then you will obtain infinite recursion.`):
		print(`This is a helper function for SpecificIrreducibleBoundedScoringPathsNumber after ensuring that startingLevel<>0.`):
		print(`Try seq(SpecificIrreducibleBoundedScoringPathsNumberHelp({[1,0],[1,-2],[2,-3],[0,1]},3,-2,n1,0),n1=0..9).`):


	elif args=`SemiBoundedScoringPathsEqsVars` then
		print(`SemiBoundedScoringPathsEqsVars(possibleScores,lowerLimit,t,k,f,g): inputs possible scores and lowerLimit for the y-value.`):
		print(`possibleScores is a list or set of scores in the format [x-step,y-step].`):
		print(`Assumes lowerLimit<=0.`):
		print(`f,g,k are placeholders to keep track of the different g.f.s. g is for irreducible (only hits minimum at the respective endpoint) walks.`):
		print(`f[a,b] is the g.f. for the number of walks from (0,a) to (n,b) that stay above the lowerLimit.`):
		print(`g[a,b] is the g.f. for the number of irreducible walks from (0,a) to (n,b) that stay above the lowerLimit.`):
		print(`f,g are the same g.f. as for EqualSemiBoundedScoringPathsEqsVars(possibleScores,lowerLimit,t,f,g).`):
		print(`k[y] is the g.f. for walks above lowerLimit that begin at (0,y) and end anywhere.`):
		print(`t is the variable in the g.f.s.`):
		print(`Outputs a set of equations (and the used variables) relating the g.f.s held in f,g,k.`):
		print(`Try SemiBoundedScoringPathsEqsVars({[1,0],[1,-2],[2,-3],[0,1]},-2,t,k,f,g).`):
	elif args=`AllSemiBoundedScoringPaths` then
		print(`AllSemiBoundedScoringPaths(possibleScores,lowerLimit,t,k,f,g): inputs possible scores and lowerLimit for the y-value.`):
		print(`possibleScores is a list or set of scores in the format [x-step,y-step].`):
		print(`Assumes lowerLimit<=0.`):
		print(`f,g,k are placeholders to keep track of the different g.f.s. g is for irreducible (only hits minimum at the respective endpoint) walks.`):
		print(`f[a,b] is the g.f. for the number of walks from (0,a) to (n,b) that stay above the lowerLimit.`):
		print(`g[a,b] is the g.f. for the number of irreducible walks from (0,a) to (n,b) that stay above the lowerLimit.`):
		print(`f,g are the same g.f. as for EqualSemiBoundedScoringPathsEqsVars(possibleScores,lowerLimit,t,f,g).`):
		print(`k[y] is the g.f. for walks above lowerLimit that begin at (0,y) and end anywhere.`):
		print(`t is the variable in the g.f.s.`):
		print(`Uses SemiBoundedScoringPathsEqsVars to find the relating equations.`):
		print(`Then uses an "optimal" Groebner basis for each variable to find a polynomial of which the variable is a root.`):
		print(`May take a while since it is computing a Groebner basis for each variable.`):
		print(`Try AllSemiBoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,1]},-2,t,k,f,g).`):
	elif args=`SpecificSemiBoundedScoringPaths` then
		print(`SpecificSemiBoundedScoringPaths(possibleScores,lowerLimit,t,k,var): The same as AllSemiBoundedScoringPaths except it only outputs the given g.f. held in var.`):
		print(`inputs possible scores and lowerLimit for the y-value.`):
		print(`possibleScores is a list or set of scores in the format [x-step,y-step].`):
		print(`Assumes lowerLimit<=0.`):
		print(`k[y] is the g.f. for walks above lowerLimit that begin at (0,y) and end anywhere.`):
		print(`t is the variable in the g.f.s.`):
		print(`var should be k[y]. It should be one of the variables returned in SemiBoundedScoringPathsEqsVars.`):
		print(`If var is f[a,b] or g[a,b] it is better (and necessary) to use SpecificEqualSemiBoundedScoringPaths(possibleScores,lowerLimit,t,f,g,var).`):
		print(`Uses SemiBoundedScoringPathsEqsVars to find the relating equations.`):
		print(`Then uses an "optimal" Groebner basis to find a polynomial of which the variable specified in var is a root.`):
		print(`Try SpecificSemiBoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,2]},-2,t,k,k[-1]).`):
	elif args=`SemiBoundedScoringPaths` then
		print(`SemiBoundedScoringPaths(possibleScores,lowerLimit,t,k): This is shorthand for computing SpecificSemiBoundedScoringPaths with var=k[0].`):
		print(`inputs possible scores and lowerLimit for the y-value.`):
		print(`possibleScores is a list or set of scores in the format [x-step,y-step].`):
		print(`Assumes lowerLimit<=0.`):
		print(`k is the g.f. for walks above lowerLimit that begin at the origin, ends anywhere, and stays above y=lowerLimit.`):
		print(`t is the variable in the g.f.`):
		print(`Uses SemiBoundedScoringPathsEqsVars to find the relating equations.`):
		print(`Then uses an "optimal" Groebner basis to find a polynomial with a root equal to the g.f. for the number of walks that start at (0,0), have length n, and do not go below y=lowerLimit.`):
		print(`Try SemiBoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,2]},-2,t,k).`):
	elif args=`SemiBoundedScoringPathsSeries` then
		print(`SemiBoundedScoringPathsSeries(possibleScores,lowerLimit,t,k,f,g,N): inputs possible scores and the lowerLimit for the y-value.`):
		print(`possibleScores is a list or set of scores in the format [x-step,y-step].`):
		print(`Assumes that lowerLimit<=0.`):
		print(`Outputs the Taylor series (truncated at t^(N+1)) of the g.f.s found in SemiBoundedScoringPathsEqsVars.`):
		print(`f,g,k are placeholders to keep track of the different g.f.s. g is for irreducible (only hits minimum at the respective endpoint) walks.`):
		print(`f[a,b] is the g.f. for the number of walks from (0,a) to (n,b) that stay above the lowerLimit.`):
		print(`g[a,b] is the g.f. for the number of irreducible walks from (0,a) to (n,b) that stay above the lowerLimit.`):
		print(`f,g are the same g.f. as for EqualSemiBoundedScoringPathsEqsVars(possibleScores,lowerLimit,t,f,g).`):
		print(`k[y] is the g.f. for walks above lowerLimit that begin at (0,y) and end anywhere.`):
		print(`t is the variable in the g.f.s.`):
		print(`Uses SemiBoundedScoringPathsEqsVars. Iterates to a fixed-point solution starting from {f[a,a]=1,f[a,b]=0,g[0,0]=0,g[a,b]=0,k[y]=1}.`):
		print(`Try SemiBoundedScoringPathsSeries({[1,0],[1,-2],[2,-3],[0,1]},-2,t,k,f,g,9).`):
	elif args=`SemiBoundedScoringPathsCoefficients` then
		print(`SemiBoundedScoringPathsCoefficients(possibleScores,lowerLimit,N): inputs possible scores and the lowerLimit for the y-value.`):
		print(`possibleScores is a list or set of scores in the format [x-step,y-step].`):
		print(`Assumes that lowerLimit<=0.`):
		print(`Outputs the first N+1 coefficients of the Taylor series of the g.f. for the number of paths that start at (0,0), stay above the lowerLimit, and end anywhere.`):
		print(`Uses SemiBoundedScoringPathsSeries or SemiBoundedScoringPaths.`):
		print(`Try SemiBoundedScoringPathsCoefficients({[1,0],[1,-2],[2,-3],[0,1]},-2,9).`):
	elif args=`SpecificSemiBoundedScoringPathsNumber` then
		print(`SpecificSemiBoundedScoringPathsNumber(possibleScores,lowerLimit,n,startingLevel): inputs possible scores (x-step,y-step) in the format {[x1,y1],[x2,y2],...}.`):
		print(`Outputs the number of paths that start at (0,startingLevel), are of length n, and don't go below y=lowerLimit<=0.`):
		print(`There CANNOT be (0,+) steps. Otherwise you get infinity as the answer.`):
		print(`The answer is computed empirically.`):
		print(`Try seq(SpecificSemiBoundedScoringPathsNumber({[1,0],[1,-2],[2,-3],[0,1]},-2,n1,0),n1=0..9).`):
	elif args=`SpecificSemiBoundedScoringPathsNumberHelp` then
		print(`A helper function to compute the number of walks after parsing in SpecificBoundedScoringPathsNumber.`):
	elif args=`SpecificSemiBoundedScoringPathsNumberHelp2` then
		print(`Computes the number of ways to end at the current x value. All of these scores correspond to steps straight down.`):

	elif args=`EqualSemiBoundedScoringPathsEqsVars` then
		print(`EqualSemiBoundedScoringPathsEqsVars(possibleScores,lowerLimit,t,f,g): inputs possible scores and lowerLimit for the y-value.`):
		print(`possibleScores is a list or set of scores in the format [x-step,y-step].`):
		print(`Assumes lowerLimit<=0.`):
		print(`f,g are placeholders to keep track of the different g.f.s. g is for irreducible (only hits minimum at the respective endpoint) walks.`):
		print(`f[a,b] is the g.f. for the number of walks from (0,a) to (n,b) that stay above the lowerLimit.`):
		print(`g[a,b] is the g.f. for the number of irreducible walks from (0,a) to (n,b) that stay above the lowerLimit.`):
		print(`t is the variable in the g.f.s.`):
		print(`Outputs a set of equations (and the used variables) relating the g.f.s held in f and g.`):
		print(`Try EqualSemiBoundedScoringPathsEqsVars({[1,0],[1,-2],[2,-3],[0,1]},-2,t,f,g).`):
	elif args=`AllEqualSemiBoundedScoringPaths` then
		print(`AllEqualSemiBoundedScoringPaths(possibleScores,lowerLimit,t,f,g): inputs possible scores and lowerLimit for the y-value.`):
		print(`possibleScores is a list or set of scores in the format [x-step,y-step].`):
		print(`Assumes lowerLimit<=0.`):
		print(`f,g are placeholders to keep track of the different g.f.s. g is for irreducible (only hits minimum at the respective endpoint) walks.`):
		print(`f[a,b] is the g.f. for the number of walks from (0,a) to (n,b) that stay above the lowerLimit.`):
		print(`g[a,b] is the g.f. for the number of irreducible walks from (0,a) to (n,b) that stay above the lowerLimit.`):
		print(`t is the variable in the g.f.s.`):
		print(`Uses EqualSemiBoundedScoringPathsEqsVars to find the relating equations.`):
		print(`Then uses an "optimal" Groebner basis for each variable to find a polynomial of which the variable is a root.`):
		print(`May take a while since it is computing a Groebner basis and solving it for each variable.`):
		print(`Try AllEqualSemiBoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,1]},-2,t,f,g).`):
	elif args=`SpecificEqualSemiBoundedScoringPaths` then
		print(`SpecificEqualSemiBoundedScoringPaths(possibleScores,lowerLimit,t,f,g,var): The same as AllEqualSemiBoundedScoringPaths except it only outputs the given g.f. held in var.`):
		print(`inputs possible scores and lowerLimit for the y-value.`):
		print(`possibleScores is a list or set of scores in the format [x-step,y-step].`):
		print(`Assumes lowerLimit<=0.`):
		print(`f,g are placeholders to keep track of the different g.f.s. g is for irreducible (only hits minimum at the respective endpoint) walks.`):
		print(`f[a,b] is the g.f. for the number of walks from (0,a) to (n,b) that stay above the lowerLimit.`):
		print(`g[a,b] is the g.f. for the number of irreducible walks from (0,a) to (n,b) that stay above the lowerLimit.`):
		print(`t is the variable in the g.f.s.`):
		print(`var is either f[a,b] or g[a,b]. It should be one of the variables returned in EqualSemiBoundedScoringPathsEqsVars.`):
		print(`Uses EqualSemiBoundedScoringPathsEqsVars to find the relating equations.`):
		print(`Then uses an "optimal" Groebner basis to find a polynomial of which the variable specified in var is a root.`):
		print(`Try SpecificEqualSemiBoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,1]},-2,t,f,g,f[-1,-1]).`):
	elif args=`EqualSemiBoundedScoringPaths` then
		print(`EqualSemiBoundedScoringPaths(possibleScores,lowerLimit,t,f): This is shorthand for computing SpecificEqualSemiBoundedScoringPaths with var=f[0,0].`):
		print(`inputs possible scores and lowerLimit for the y-value.`):
		print(`possibleScores is a list or set of scores in the format [x-step,y-step].`):
		print(`Assumes lowerLimit<=0.`):
		print(`t is the variable in the g.f.`):
		print(`Uses EqualSemiBoundedScoringPathsEqsVars to find the relating equations.`):
		print(`Then uses an "optimal" Groebner basis to find a polynomial with a root equal to the g.f. for the number of walks that start at (0,0), end at (n,0) and do not go below y=lowerLimit.`):
		print(`Try EqualSemiBoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,1]},-2,t,f).`):
	elif args=`EqualSemiBoundedScoringPathsSeries` then
		print(`EqualSemiBoundedScoringPathsSeries(possibleScores,lowerLimit,t,f,g,N): inputs possible scores and the lowerLimit for the y-value.`):
		print(`possibleScores is a list or set of scores in the format [x-step,y-step].`):
		print(`Assumes that lowerLimit<=0.`):
		print(`Outputs the Taylor series (truncated at t^(N+1)) of the g.f.s found in EqualSemiBoundedScoringPathsEqsVars.`):
		print(`f[a,b] is the g.f. for the number of walks from (0,a) to (n,b) that stay above the lowerLimit.`):
		print(`g[a,b] is the g.f. for the number of irreducible walks from (0,a) to (n,b) that stay above the lowerLimit.`):
		print(`Uses EqualSemiBoundedScoringPathsEqsVars. Iterates to a fixed-point solution starting from {f[a,a]=1,f[a,b]=0,g[0,0]=0,g[a,b]=0}.`):
		print(`Try EqualSemiBoundedScoringPathsSeries({[1,0],[1,-2],[2,-3],[0,1]},-2,t,f,g,9).`):
	elif args=`EqualSemiBoundedScoringPathsCoefficients` then
		print(`EqualSemiBoundedScoringPathsCoefficients(possibleScores,lowerLimit,N): inputs possible scores and the lowerLimit for the y-value.`):
		print(`possibleScores is a list or set of scores in the format [x-step,y-step].`):
		print(`Assumes that lowerLimit<=0.`):
		print(`Outputs the first N+1 coefficients of the Taylor series of the g.f. for the number of paths that start at (0,0), end at (n,0), and stay above the lowerLimit.`):
		print(`Uses EqualSemiBoundedScoringPaths.`):
		print(`Try EqualSemiBoundedScoringPathsCoefficients({[1,0],[1,-2],[2,-3],[0,1]},-2,9).`):
	elif args=`SpecificEqualSemiBoundedScoringPathsNumber` then
		print(`SpecificEqualSemiBoundedScoringPathsNumber(possibleScores,lowerLimit,n,startingLevel): inputs possible scores (x-step,y-step) in the format {[x1,y1],[x2,y2],...}.`):
		print(`Outputs the number of paths that start at (0,startingLevel), end at (n,0), and don't go below y=lowerLimit<=0.`):
		print(`There should NOT be BOTH (0,+) and (0,-) steps. Otherwise you get infinity as the answer.`):
		print(`The answer is computed empirically.`):
		print(`Try seq(SpecificEqualSemiBoundedScoringPathsNumber({[1,0],[1,-2],[2,-3],[0,1]},-2,n1,0),n1=0..9).`):
	elif args=`SpecificEqualSemiBoundedScoringPathsNumberN` then
		print(`SpecificEqualSemiBoundedScoringPathsNumberN(possibleScores,lowerLimit,n,startingLevel): A helper function that works if there are no (0,+) steps.`):
		print(`Helps SpecificEqualSemiBoundedScoringPathsNumber.`):
	elif args=`SpecificEqualSemiBoundedScoringPathsNumberP` then
		print(`SpecificEqualSemiBoundedScoringPathsNumberP(possibleScores,lowerLimit,n,startingLevel,maxDescentRate): A helper function that works if there are no (0,-) steps.`):
		print(`Helps SpecificEqualSemiBoundedScoringPathsNumber.`):
		print(`Inputting maxDescentRate just saves some computation time. It is determined by possibleScores.`):
	elif args=`SpecificIrreducibleSemiBoundedScoringPathsNumber` then
		print(`SpecificIrreducibleSemiBoundedScoringPathsNumber(possibleScores,n,startingLevel): inputs possible scores (x-step,y-step) in the format {[x1,y1],[x2,y2],...}.`):
		print(`Outputs the number of paths that start at (0,startingLevel), end at (n,0), and don't go below y=min(startingLevel,0).`):
		print(`There should NOT be BOTH (0,+) and (0,-) steps. Otherwise you get infinity as the answer.`):
		print(`The answer is computed empirically.`):
		print(`Try seq(SpecificIrreducibleSemiBoundedScoringPathsNumber({[1,0],[1,-2],[2,-3],[0,1]},-2,n1,0),n1=0..9).`):
	elif args=`SpecificIrreducibleSemiBoundedScoringPathsNumberN` then
		print(`SpecificIrreducibleSemiBoundedScoringPathsNumberN(possibleScores,lowerLimit,n,startingLevel): A helper function that works if there are no (0,+) steps.`):
		print(`Helps SpecificIrreducibleSemiBoundedScoringPathsNumber.`):
	elif args=`SpecificIrreducibleSemiBoundedScoringPathsNumberP` then
		print(`SpecificIrreducibleSemiBoundedScoringPathsNumberP(possibleScores,lowerLimit,n,startingLevel,maxDescentRate): A helper function that works if there are no (0,-) steps.`):
		print(`Helps SpecificIrreducibleSemiBoundedScoringPathsNumber.`):
		print(`Inputting maxDescentRate just saves some computation time. It is determined by possibleScores.`):


	elif args=`UnBoundedScoringPaths` then
		print(`UnBoundedScoringPaths(possibleScores,t): inputs possible scores and a variable t.`):
		print(`possibleScores is a list or set of scores in the format [x,y]. All scores should have x-step>0.`):
		print(`t is the variable in the generating function.`):
		print(`Outputs the generating function for the number of paths of total x-step, n, that start at (0,0) and go anywhere.`):
		print(`Try UnBoundedScoringPaths({[1,1],[1,-1],[1,2],[2,-3]},t).`):
	elif args=`UnBoundedScoringPathsCoefficients` then
		print(`UnBoundedScoringPathsCoefficients(possibleScores,N): Outputs the first N+1 coefficients of the Taylor series of the g.f. for the number of paths of x-step n.`):
		print(`Uses UnBoundedScoringPaths.`):
		print(`Try UnBoundedScoringPathsCoefficients({[1,1],[1,-1],[1,2],[2,-3]},9).`):
	elif args=`UnBoundedScoringPathsNumber` then
		print(`UnBoundedScoringPathsNumber(possibleScores,n): inputs possible scores (x-step,y-step)`):
		print(`Outputs the number of paths of total x-step, n, that start at (0,0) and go anywhere.`):
		print(`The answer is computed empirically.`):
		print(`Try seq(UnBoundedScoringPathsNumber({[1,1],[1,-1],[1,2],[2,-3]},n1),n1=0..9).`):
	elif args=`UnBoundedScoringPathsNumber1` then
		print(`UnBoundedScoringPathsNumber1(possibleScores,n): A helper function after ensuring no infinite recursion.`):
		print(`Helps UnBoundedScoringPathsNumber.`):

	elif args=`EqualUnBoundedScoringPathsEqsVars` then
		print(`EqualUnBoundedScoringPathsEqsVars(possibleScores,t,GF,h,f,g): inputs possible scores (x-step,y-step), and placeholder variables.`):
		print(`possibleScores is a list or set of scores in the format [x,y].`):
		print(`h,f,g are placeholder generating functions. f,g are the same g.f. as for EqualSemiBoundedScoringPathsEqsVars(possibleScores,0,t,f,g).`):
		print(`GF is the g.f. for the number of walks that begin at (0,0) and end at (n,0).`):
		print(`h[i] is the g.f. for walks that start at (i,0) and return to the x-axis only at the end.`):
		print(`f[a,b] is the g.f. for the number of walks from (0,a) to (n,b) that stay above the x-axis.`):
		print(`g[a,b] is the g.f. for the number of irreducible (only hits min(a,b) at the respective endpoint) walks from (0,a) to (n,b) that stay above the x-axis.`):
		print(`t is the variable in the generating functions.`):
		print(`Outputs many generating functions and equations describing their relations.`):
		print(`The g.f. for the number of walks that begin at (0,0) and end at (n,0) is found by taking 1/(1-h[0]). (Since gf=1+h[0]*gf.)`):
		print(`Alternately, a more "optimal" solution for gf is found after adding the equation in the function EqualUnBoundedScoringPaths.`):
		print(`Try EqualUnBoundedScoringPathsEqsVars({[1,0],[1,-2],[2,-3],[0,1]},t,GF,h,f,g).`):
	elif args=`AllEqualUnBoundedScoringPaths` then
		print(`AllEqualUnBoundedScoringPaths(possibleScores,t,GF,h,f,g): inputs possible scores (x-step,y-step), and placeholder variables.`):
		print(`possibleScores is a list or set of scores in the format [x,y].`):
		print(`h,f,g are placeholder generating functions. f,g are the same g.f. as for EqualSemiBoundedScoringPathsEqsVars(possibleScores,0,t,f,g).`):
		print(`GF is the g.f. for the number of walks that begin at (0,0) and end at (n,0).`):
		print(`h[i] is the g.f. for walks that start at (i,0) and return to the x-axis only at the end.`):
		print(`f[a,b] is the g.f. for the number of walks from (0,a) to (n,b) that stay above the x-axis.`):
		print(`g[a,b] is the g.f. for the number of irreducible (only hits min(a,b) at the respective endpoint) walks from (0,a) to (n,b) that stay above the x-axis.`):
		print(`t is the variable in the generating functions.`):
		print(`Uses EqualUnBoundedScoringPathsEqsVars to find the relating equations.`):
		print(`Then uses an "optimal" Groebner basis for each variable to find a polynomial of which the variable is a root.`):
		print(`May take a while since it is computing a Groebner basis for each variable.`):
		print(`Try AllEqualUnBoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,1]},t,GF,h,f,g).`):
	elif args=`SpecificUnBoundedScoringPaths` then
		print(`SpecificUnBoundedScoringPaths(possibleScores,lowerLimit,t,GF,h,var): The same as AllEqualUnBoundedScoringPaths except it only outputs the given g.f. held in var.`):
		print(`inputs possible scores (x-step,y-step), and placeholder variables.`):
		print(`possibleScores is a list or set of scores in the format [x,y].`):
		print(`GF is the g.f. for the number of walks that begin at (0,0) and end at (n,0).`):
		print(`h[i] is the g.f. for walks that start at (i,0) and return to the x-axis only at the end.`):
		print(`t is the variable in the generating functions.`):
		print(`var is either GF or h[i]. It should be one of the variables returned in EqualUnBoundedScoringPathsEqsVars.`):
		print(`If var is f[a,b] or g[a,b] it is better to use SpecificEqualSemiBoundedScoringPaths(possibleScores,0,t,f,g,var).`):
		print(`Uses EqualUnBoundedScoringPathsEqsVars to find the relating equations.`):
		print(`Then uses an "optimal" Groebner basis to find a polynomial of which var is a root.`):
		print(`Try SpecificUnBoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,1]},t,GF,h,h[1]).`):
	elif args=`EqualUnBoundedScoringPaths` then
		print(`EqualUnBoundedScoringPaths(possibleScores,t,GF): inputs possible scores and a variable t.`):
		print(`possibleScores is a list or set of scores in the format [x-step,y-step].`):
		print(`GF is the g.f. for the number of walks that begin at (0,0) and end at (n,0).`);
		print(`t is the variable in the g.f.`):
		print(`Uses EqualUnBoundedScoringPathsEqsVars to find the relating equations.`):
		print(`Then uses an "optimal" Groebner basis to find a polynomial with a root equal to GF.`):
		print(`Try EqualUnBoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,1]},t,GF).`):
	elif args=`EqualUnBoundedScoringPathsSeries` then
		print(`EqualUnBoundedScoringPathsSeries(possibleScores,t,GF,h,f,g,N): inputs possible scores which is a list or set of scores in the format [x-step,y-step].`):
		print(`Outputs the Taylor series (truncated at t^(N+1)) of the g.f.s found in UnBoundedScoringPathsEqsVars.`):
		print(`h,f,g are placeholder generating functions. f,g are the same g.f. as for EqualSemiBoundedScoringPathsEqsVars(possibleScores,0,t,f,g).`):
		print(`GF is the g.f. for the number of walks that begin at (0,0) and end at (n,0).`):
		print(`h[i] is the g.f. for walks that start at (i,0) and return to the x-axis only at the end.`):
		print(`f[a,b] is the g.f. for the number of walks from (0,a) to (n,b) that stay above the x-axis.`):
		print(`g[a,b] is the g.f. for the number of irreducible (only hits min(a,b) at the respective endpoint) walks from (0,a) to (n,b) that stay above the x-axis.`):
		print(`Uses UnBoundedScoringPathsEqsVars. Iterates to a fixed-point solution starting from {GF=1,f[a,a]=1,f[a,b]=0,h[i]=0,g[0,0]=0,g[a,b]=0}.`):
		print(`Try EqualUnBoundedScoringPathsSeries({[1,0],[1,-2],[2,-3],[0,1]},t,GF,h,f,g,9).`):
	elif args=`EqualUnBoundedScoringPathsCoefficients` then
		print(`EqualUnBoundedScoringPathsCoefficients(possibleScores,N): inputs possible scores which is a list or set of scores in the format [x-step,y-step].`):
		print(`Outputs the first N+1 coefficients of the Taylor series of the g.f. for the number of paths that start at (0,0) and end at (n,0).`):
		print(`Uses EqualUnBoundedScoringPathsSeries since iterating in the minimal polynomial does not work due to some convergence issue.`):
		print(`Try EqualUnBoundedScoringPathsCoefficients({[1,0],[1,-2],[2,-3],[0,1]},9).`):
	elif args=`SpecificUnBoundedScoringPathsNumber` then
		print(`SpecificUnBoundedScoringPathsNumber(possibleScores,n,startingLevel): inputs possible scores (x-step,y-step) in the format {[x1,y1],[x2,y2],...}.`):
		print(`Outputs the number of paths that start at (0,startingLevel) and end at (n,0).`):
		print(`There should NOT be BOTH (0,+) and (0,-) steps. Otherwise you get infinity as the answer.`):
		print(`The answer is computed empirically.`):
		print(`Try seq(SpecificUnBoundedScoringPathsNumber({[1,0],[1,-2],[2,-3],[0,1]},n1,0),n1=0..9).`):
	elif args=`SpecificUnBoundedScoringPathsNumberN` then
		print(`SpecificUnBoundedScoringPathsNumberN(possibleScores,n,startingLevel,maxRate): A helper function that works if there are no (0,+) steps.`):
		print(`Helps SpecificUnBoundedScoringPathsNumber.`):
		print(`Inputting maxRate just saves some computation time. It is determined by possibleScores.`):
	elif args=`SpecificUnBoundedScoringPathsNumberP` then
		print(`SpecificUnBoundedScoringPathsNumberP(possibleScores,n,startingLevel,maxRate): A helper function that works if there are no (0,-) steps.`):
		print(`Helps SpecificUnBoundedScoringPathsNumber.`):
		print(`Inputting maxRate just saves some computation time. It is determined by possibleScores.`):
	elif args=`SpecificIrreducibleUnBoundedScoringPathsNumber` then
		print(`SpecificIrreducibleUnBoundedScoringPathsNumber(possibleScores,n,startingLevel): inputs possible scores (x-step,y-step) in the format {[x1,y1],[x2,y2],...}.`):
		print(`Outputs the number of paths that start at (0,startingLevel) and end at (n,0).`):
		print(`There should NOT be BOTH (0,+) and (0,-) steps. Otherwise you get infinity as the answer.`):
		print(`The answer is computed empirically.`):
		print(`Try seq(SpecificIrreducibleUnBoundedScoringPathsNumber({[1,0],[1,-2],[2,-3],[0,1]},n1,0),n1=0..9).`):
	elif args=`SpecificIrreducibleUnBoundedScoringPathsNumberN` then
		print(`SpecificIrreducibleUnBoundedScoringPathsNumberN(possibleScores,n,startingLevel,maxRate): A helper function that works if there are no (0,+) steps.`):
		print(`Helps SpecificIrreducibleUnBoundedScoringPathsNumber.`):
		print(`Inputting maxRate just saves some computation time. It is determined by possibleScores.`):
	elif args=`SpecificIrreducibleUnBoundedScoringPathsNumberP` then
		print(`SpecificIrreducibleUnBoundedScoringPathsNumberP(possibleScores,n,startingLevel,maxRate): A helper function that works if there are no (0,-) steps.`):
		print(`Helps SpecificIrreducibleUnBoundedScoringPathsNumber.`):
		print(`Inputting maxRate just saves some computation time. It is determined by possibleScores.`):


	elif args=`FindProperRoot` then
		print(`FindProperRoot(possibleScores,P,case,var,t,lowerLimit:=0):`):
		print(`Finding the proper root (minimal polynomial). This is only a helper function for parsing output from`):
		print(`SemiBoundedScoringPaths, EqualSemiBoundedScoringPaths, and EqualUnBoundedScoringPaths.`):
		print(`(And other functions that output minimal polynomial(s).)`):
		print(`S is the step set. P is the polynomial originally found (not necessarily minimal).`):
		print(`var is the g.f. we are interested in.`):
		print(`t is the variable in the g.f.`):
		print(`lowerLimit is applicable for the semi-bounded cases.`):
		print(`case=1: semi bounded free, 2: semi bounded equal, 3: semi bounded irreducible, 4: unbounded equal, 5: unbounded irreducible`):

	elif args=`RatioOfWalks` then
		print(`RatioOfWalks(S): Inputs a step set S.`):
		print(`Outputs the asymptotic ratio of nonnegative excursions to bridges after finding recurrences for excursions and bridges respectively.`):

	elif args=`AsymptoticBase` then
		print(`AsymptoticBase(P,t,f): inputs a minimal polynomial in t and f: the generating function with formal variable t.`):
		print(`Outputs the b such that the coefficients of f follow b^n asymptotically.`):

	elif args=`RandomStepSet` then
		print(`RandomStepSet(numSteps,xBound,yBound): produces a set of steps each with x-value between 0 and xBound and y-value between -yBound and yBound.`):
		print(`There cannot be a [0,0] step, thus we need yBound>=1. And there must be some step with x-value>0 so we should also have xBound>=1.`):
		print(`Try RandomStepSet(8,4,3).`):
	elif args=`RandomZeroStepSet` then
		print(`Producing random step sets for the Bounded, EqualBounded, EqualSemiBounded and EqualUnBounded cases.`):
		print(`The bounded cases could be expanded to include some pairs of (+) and (-) zero steps but I do not currently know the restrictions.`):
		print(`RandomZeroStepSet(numSteps,xBound,yBound): produces a set of steps each with x-value between 0 and xBound and y-value between -yBound and yBound.`):
		print(`There cannot be a [0,0] step, thus we need yBound>=1. And there must be some step with x-value>0 so we should also have xBound>=1.`):
		print(`There also cannot be both (+) and (-) zero steps.`):
		print(`Try RandomZeroStepSet(8,4,3).`):
	elif args=`RandomSemiBoundedStepSet` then
		print(`Producing random step sets for the SemiBounded case.`):
		print(`RandomSemiBoundedStepSet(numSteps,xBound,yBound): produces a set of steps each with x-value between 0 and xBound and y-value between -yBound and yBound.`):
		print(`There cannot be a [0,0] step, thus we need yBound>=1. And there must be some step with x-value>0 so we should also have xBound>=1.`):
		print(`There also cannot be (+) zero steps.`):
		print(`Try RandomSemiBoundedStepSet(8,4,3).`):
	elif args=`RandomUnBoundedStepSet` then
		print(`Producing random step sets for the UnBounded case.`):
		print(`RandomUnBoundedStepSet(numSteps,xBound,yBound): produces a set of steps each with x-value between 1 and xBound and y-value between -yBound and yBound.`):
		print(`There cannot be zero steps.`):
		print(`Try RandomUnBoundedStepSet(8,4,3).`):


	elif args=`PaperBounded1` then
		print(`PaperBounded1(Steps,upperLimit,lowerLimit,K): inputs a set of steps, Steps, and a positive integer K, outputs an article about the`):
		print(`generating function of walks from the origin back to [n,0] for some n that never go below y=lowerLimit or above y=upperLimit.`):
		print(`It also outputs the first K terms, for the sake of the OEIS.`):
		print(`Try:`):
		print(`PaperBounded1({[1,2],[1,-2],[1,3],[1,-3]},8,-2,20);`):
	elif args=`PaperBounded2` then
		print(`PaperBounded2(Steps,K): inputs a set of steps, Steps, and a positive integer K, outputs an article about the`):
		print(`generating function of walks from the origin back to [n,0] for some n that never go further than one step from the x-axis.`):
		print(`It also outputs the first K terms, for the sake of the OEIS.`):
		print(`Try:`):
		print(`PaperBounded2({[1,2],[1,-2],[1,3],[1,-3]},20);`):
	elif args=`PaperSemiBounded` then
		print(`PaperSemiBounded(Steps,K,recurrence:=false,asymptotic:=false): inputs a set of steps, Steps, and a positive integer K, outputs an article about the`):
		print(`generating function of walks from the origin back to [n,0] for some n that never go below the x-axis.`):
		print(`It also outputs the first K terms, for the sake of the OEIS.`):
		print(`If you want a linear recurrence (with polynomial coefficients) for the terms of the g.f., set recurrence=true.`):
		print(`If you want the asymptotic behavior, set asympototic=true. This will not always work.`):
		print(`Try:`):
		print(`PaperSemiBounded({[1,1],[1,-1],[1,0]},40,true);`):
	elif args=`PaperEqualSemiBounded` then
		print(`PaperEqualSemiBounded(Steps,K,recurrence:=false,asymptotics:=false): inputs a set of steps, Steps, and a positive integer K, outputs an article about the`):
		print(`generating function of walks from the origin back to [n,0] for some n that never go below the x-axis.`):
		print(`It also outputs the first K terms, for the sake of the OEIS.`):
		print(`If you want a linear recurrence (with polynomial coefficients) for the terms of the g.f., set recurrence=true.`):
		print(`If you want the asymptotic behavior, set asympototics=true. This will not always work.`):
		print(`Try:`):
		print(`PaperEqualSemiBounded({[1,1],[1,-1],[1,0]},40,true);`):
	elif args=`BookEqualSemiBounded` then
		print(`BookEqualSemiBounded(Steps,K,recurrence:=false,asymptotics:=false): inputs a set of steps, Steps, and a positive integer K.`):
		print(`Outputs a book about the generating functions of walks from the origin back to [n,0] for some n that never go below the x-axis.`):
		print(`Iterates over all subsets of Steps.`):
		print(`If you want linear recurrences (with polynomial coefficients) for the terms of the g.f., set recurrence=true.`):
		print(`If you want the asymptotics, set asympototics=true. This will not always work.`):
		print(`Try:`):
		print(`BookEqualSemiBounded({[1,1],[1,-1],[1,0],[2,2],[2,-2],[3,3],[3,-1]},40,true);`):
	elif args=`PaperEqualSemiBoundedPair` then
		print(`PaperEqualSemiBoundedPair(a,b,K,recurrence:=false,asymptotics:=false): inputs positive integers a and b that are relatively prime, and a positive integer K, outputs an article about the`):
		print(`generating function of walks from the origin back to [(a+b)*n,0] for some n that never go below the x-axis.`):
		print(`It also outputs the first K terms, for the sake of the OEIS.`):
		print(`If you want a linear recurrence (with polynomial coefficients) for the terms of the g.f., set recurrence=true.`):
		print(`If you want the asymptotic behavior, set asympototics=true. This will not always work.`):
		print(`Try:`):
		print(`PaperEqualSemiBoundedPair(2,3,40,true);`):
	elif args=`PaperUnBounded` then
		print(`PaperUnBounded(Steps,K,recurrence:=false,asymptotics:=false): inputs a set of steps, Steps, and a positive integer K, outputs an article about the`):
		print(`generating function of walks from the origin back to [n,0] for some n.`):
		print(`It also outputs the first K terms, for the sake of the OEIS.`):
		print(`If you want a linear recurrence (with polynomial coefficients) for the terms of the g.f., set recurrence=true.`):
		print(`If you want the asymptotic behavior, set asympototics=true. This will not always work.`):
		print(`Try:`):
		print(`PaperUnBounded({[1,1],[1,-1],[1,0]},40,true);`):
	elif args=`PaperUnBoundedPair` then
		print(`PaperUnBoundedPair(a,b,K,recurrence:=false,asymptotics:=false): inputs positive integers a and b that are relatively prime, and a positive integer K, outputs an article about the`):
		print(`generating function of walks from the origin back to [(a+b)*n,0] for some n.`):
		print(`It also outputs the first K terms, for the sake of the OEIS.`):
		print(`If you want a linear recurrence (with polynomial coefficients) for the terms of the g.f., set recurrence=true.`):
		print(`If you want the asymptotic behavior, set asympototics=true. This will not always work.`):
		print(`Try:`):
		print(`PaperUnBoundedPair(2,3,40,true);`):


	else
		print(`Primary functions:`):
		print(`BoundedScoringPaths, EqualBoundedScoringPaths, SemiBoundedScoringPaths, EqualSemiBoundedScoringPaths, UnBoundedScoringPaths, EqualUnBoundedScoringPaths`):
		print(``):
		print(`Secondary functions:`):
		print(`BoundedScoringPathsEqsVars, AllBoundedScoringPaths, EqualBoundedScoringPathsEqsVars, AllEqualBoundedScoringPaths, SemiBoundedScoringPathsEqsVars, AllSemiBoundedScoringPaths, SpecificSemiBoundedScoringPaths, EqualSemiBoundedScoringPathsEqsVars, AllEqualSemiBoundedScoringPaths, SpecificEqualSemiBoundedScoringPaths, EqualUnBoundedScoringPathsEqsVars, AllEqualUnBoundedScoringPaths, SpecificUnBoundedScoringPaths`):
		print(``):
		print(`Enumerating functions:`):
		print(`BoundedScoringPathsCoefficients, BoundedScoringPathsNumber, EqualBoundedScoringPathsCoefficients, SpecificEqualBoundedScoringPathsNumber, SpecificIrreducibleBoundedScoringPathsNumber, SemiBoundedScoringPathsCoefficients, SpecificSemiBoundedScoringPathsNumber, EqualSemiBoundedScoringPathsCoefficients, SpecificEqualSemiBoundedScoringPathsNumber, SpecificIrreducibleSemiBoundedScoringPathsNumber, UnBoundedScoringPathsCoefficients, UnBoundedScoringPathsNumber, EqualUnBoundedScoringPathsCoefficients, SpecificUnBoundedScoringPathsNumber, SpecificIrreducibleUnBoundedScoringPathsNumber`):
		print(``):
		print(`Helper functions:`):
		print(`SpecificIrreducibleBoundedScoringPathsNumberHelp, SemiBoundedScoringPathsSeries, SpecificSemiBoundedScoringPathsNumberHelp, SpecificSemiBoundedScoringPathsNumberHelp2, EqualSemiBoundedScoringPathsSeries, SpecificEqualSemiBoundedScoringPathsNumberN, SpecificEqualSemiBoundedScoringPathsNumberP, SpecificIrreducibleSemiBoundedScoringPathsNumberN, SpecificIrreducibleSemiBoundedScoringPathsNumberP, UnBoundedScoringPathsNumber1, EqualUnBoundedScoringPathsSeries, SpecificUnBoundedScoringPathsNumberN, SpecificUnBoundedScoringPathsNumberP, SpecificIrreducibleUnBoundedScoringPathsNumberN, SpecificIrreducibleUnBoundedScoringPathsNumberP, FindProperRoot, RatioOfWalks, AsymptoticBase, RandomStepSet, RandomZeroStepSet, RandomSemiBoundedStepSet, RandomUnBoundedStepSet`):
		print(``):
		print(`Paper-producing functions:`):
		print(`PaperBounded1, PaperBounded2, PaperSemiBounded, PaperEqualSemiBounded, BookEqualSemiBounded, PaperEqualSemiBoundedPair, PaperUnBounded, PaperUnBoundedPair`):
		print(``):
		print(`If all the y-steps have a common factor, it may be necessary to divide by that factor.`):
		print(`In the unbounded and semi-bounded (equal) cases, you can have steps with x-step=0, as long as they are all in the same direction. E.g. {[0,1],[0,-2]} is NOT allowed.`):
		print(`You cannot have steps with x-step=0 in the free (walk anywhere) unbounded case.`):
		print(`You cannot have steps with x-step=0 and y-step>0 in the free (walk anywhere) semibounded case.`):
		print(`In the bounded case, you can have BOTH [0,-] and [0,+] ONLY IF you cannot sum to 0 within the bounds.`):
		print(`Generating function(s) will be abbreviated as g.f.(s).`):
		print(``):
		print(`For functions borrowed from other packages, type ezraSCHUTZENBERGER(), ezraEKHAD(), ezraW1D(), and ezraAsyRec.`):
	fi:
end:
###########################################################################################################################################


###########################################################################################################################################
###########################################################################################################################################
#Bounded
###########################################################################################################################################
###########################################################################################################################################
#BoundedScoringPathsEqsVars(possibleScores,upperLimit,lowerLimit,t,f): inputs possible scores and the upperLimit and lowerLimit for the y-value.
#possibleScores is a list or set of scores in the format [x-step,y-step].
#Assumes that upperLimit>=0,lowerLimit<=0.
#f is a placeholder to keep track of the different g.f.s.
#t is the variable in the g.f.
#Outputs a set of equations (and the used variables) relating the g.f.s for the number of paths of total x-step n, that start at (0,0) and stay within the limits.
#Includes the g.f.s for all pairs of limits (up,low) such that up-low=upperLimit-lowerLimit.
#All equations are linear.
#Try BoundedScoringPathsEqsVars({[1,0],[1,-2],[2,-3],[0,1]},3,-2,t,f).
BoundedScoringPathsEqsVars:=proc(possibleScores,upperLimit,lowerLimit,t,f) local eqs,vars,i,expr,score:
	eqs:={}:
	vars:={}:
	for i from lowerLimit to upperLimit do
		expr:=1:
		for score in possibleScores do
			if i+score[2]>=lowerLimit and i+score[2]<=upperLimit then
				expr:=expr+t^score[1]*f[upperLimit-i-score[2],lowerLimit-i-score[2]]:
			fi:
		od:
		eqs:=eqs union {f[upperLimit-i,lowerLimit-i]=expr}:
		vars:=vars union {f[upperLimit-i,lowerLimit-i]}:
	od:
	eqs,vars:
end:
###########################################################################################################################################
#AllBoundedScoringPaths(possibleScores,upperLimit,lowerLimit,t,f): inputs possible scores and the upperLimit and lowerLimit for the y-value.
#possibleScores is a list or set of scores in the format [x-step,y-step].
#Assumes that upperLimit>=0,lowerLimit<=0.
#f is a placeholder to keep track of the different g.f.s.
#t is the variable in the g.f.
#Outputs ALL the g.f.s for the number of paths of total x-step n, that start at (0,0) and stay within its limits.
#Includes the g.f.s for all pairs of limits (up,low) such that up-low=upperLimit-lowerLimit.
#All equations are linear so the generating functions are all rational!
#Uses BoundedScoringPathsEqsVars.
#Try AllBoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,1]},3,-2,t,f).
AllBoundedScoringPaths:=proc(possibleScores,upperLimit,lowerLimit,t,f):
	solve(BoundedScoringPathsEqsVars(possibleScores,upperLimit,lowerLimit,t,f)):
end:
###########################################################################################################################################
#BoundedScoringPaths(possibleScores,upperLimit,lowerLimit,t): inputs possible scores and the upperLimit and lowerLimit for the y-value.
#possibleScores is a list or set of scores in the format [x-step,y-step].
#Assumes that upperLimit>=0,lowerLimit<=0.
#t is the variable in the g.f.
#Outputs only the g.f. for the number of paths of total x-step n, that start at (0,0) and stay within the limits given.
#Uses AllBoundedScoringPaths.
#Try BoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,1]},3,-2,t).
BoundedScoringPaths:=proc(possibleScores,upperLimit,lowerLimit,t) local f:
	subs(AllBoundedScoringPaths(possibleScores,upperLimit,lowerLimit,t,f),f[upperLimit,lowerLimit]):
end:
######################################################################################################
#BoundedScoringPathsCoefficients(possibleScores,upperLimit,lowerLimit,N): inputs possible scores and the upperLimit and lowerLimit for the y-value.
#possibleScores is a list or set of scores in the format [x-step,y-step].
#Assumes that upperLimit>=0,lowerLimit<=0.
#Outputs the first N+1 coefficients of the Taylor series of the g.f. for the number of paths that start at (0,0) and stay within the limits.
#Uses BoundedScoringPaths.
#Try BoundedScoringPathsCoefficients({[1,0],[1,-2],[2,-3],[0,1]},3,-2,9).
BoundedScoringPathsCoefficients:=proc(possibleScores,upperLimit,lowerLimit,N) local t,f:
	f:=taylor(BoundedScoringPaths(possibleScores,upperLimit,lowerLimit,t),t,N+1):
	seq(coeff(f,t,i),i=0..N):
end:
######################################################################################################
#BoundedScoringPathsNumber(possibleScores,upperLimit,lowerLimit,n): inputs possible scores and the upperLimit and lowerLimit for the y-value.
#Assumes that upperLimit>=0,lowerLimit<=0.
#possibleScores is a list or set of scores in the format [x-step,y-step].
#Outputs the number of paths of x-step length n that start at (0,0) and stay within the bounds.
#This function operates empirically.
#If there are BOTH [0,-] and [0,+] steps that sum to 0 within the bounds, then you will obtain infinite recursion.
#Try seq(BoundedScoringPathsNumber({[1,0],[1,-2],[2,-3],[0,1]},3,-2,n1),n1=0..9).
BoundedScoringPathsNumber:=proc(possibleScores,upperLimit,lowerLimit,n) local count,score: option remember:
	if n<0 then
		return(0):
	elif n=0 then	#Could end the path now.
		count:=1:
	else
		count:=0:
	end:
	for score in possibleScores do
		if score[1]<=n and score[2]>=lowerLimit and score[2]<=upperLimit then
			count:=count+BoundedScoringPathsNumber(possibleScores,upperLimit-score[2],lowerLimit-score[2],n-score[1]):
		fi:
	od:
	count:
end:

###########################################################################################################################################
#EqualBoundedScoringPathsEqsVars(possibleScores,upperLimit,lowerLimit,t,f,e): inputs possible scores and the upperLimit and lowerLimit for the y-value.
#Assumes that upperLimit>=0,lowerLimit<=0.
#possibleScores is a list or set of scores in the format [x-step,y-step].
#f,e are placeholder variables for the g.f.s.
#f[upperLimit,lowerLimit] is the g.f. for the number of paths of total x-step n, that start at (0,0), end at (n,0), and stay within the limits.
#e[upperLimit,lowerLimit,c] is the g.f. for the number of "irreducible" walks that start at (c,0), end at (n,0), stay within the given limits, and touch 0 only at the end.
#t is the variable in the g.f.
#Outputs a set of equations (and the used variables) relating the g.f.s held in f and e.
#Try EqualBoundedScoringPathsEqsVars({[1,0],[1,-2],[2,-3],[0,1]},3,-2,t,f,e).
EqualBoundedScoringPathsEqsVars:=proc(possibleScores,upperLimit,lowerLimit,t,f,e) local expr,expr2,score,eqs,vars,i:
	expr:=0:
	expr2:=0:
	for score in possibleScores do
		if score[2]=0 then
			expr:=expr+t^score[1]:
		elif score[2]>=lowerLimit and score[2]<=upperLimit then
			expr2:=expr2+t^score[1]*e[upperLimit,lowerLimit,score[2]]:
		fi:
	od:
	eqs:={e[upperLimit,lowerLimit,0]=expr2,f[upperLimit,lowerLimit]=1+(expr+e[upperLimit,lowerLimit,0])*f[upperLimit,lowerLimit]}:	#Could equivalently put in f[upperLimit,lowerLimit]=1/(1-expr-e[u,l,0])
	vars:={f[upperLimit,lowerLimit],seq(e[upperLimit,lowerLimit,i],i in seq(lowerLimit..upperLimit))}:
	for i in {seq(lowerLimit..upperLimit)} minus {0} do	#These e[lowerLimit,upperLimit,c] are "irreducible" (hits 0 only at endpoint) paths from c back to 0.
		expr:=0:
		for score in possibleScores do
			if score[2]=-i then
				expr:=expr+t^score[1]:
			elif i+score[2]>=lowerLimit and i+score[2]<=upperLimit then
				expr:=expr+t^score[1]*e[upperLimit,lowerLimit,i+score[2]]:
			fi:
		od:
		eqs:=eqs union {e[upperLimit,lowerLimit,i]=expr}:
	od:
	eqs,vars:
end:
###########################################################################################################################################
#AllEqualBoundedScoringPaths(possibleScores,upperLimit,lowerLimit,t,f,e): inputs possible scores and the upperLimit and lowerLimit for the y-value.`):
#Assumes that upperLimit>=0,lowerLimit<=0.
#possibleScores is a list or set of scores in the format [x-step,y-step].
#f,e are placeholder variables for the g.f.s.
#f[upperLimit,lowerLimit] is the g.f. for the number of paths of total x-step n, that start at (0,0), end at (n,0), and stay within the limits.
#e[upperLimit,lowerLimit,c] is the g.f. for the number of "irreducible" walks that start at (c,0), end at (n,0), stay within the given limits, and touch 0 only at the end.
#t is the variable in the g.f.s.
#Outputs the g.f.s f[upperLimit,lowerLimit] and e[upperLimit,lowerLimit,c] for all c=lowerLimit..-1,1..upperLimit.
#Uses EqualBoundedScoringPathsEqsVars.
#Try AllEqualBoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,1]},3,-2,t,f,e).
AllEqualBoundedScoringPaths:=proc(possibleScores,upperLimit,lowerLimit,t,f,e):
	solve(EqualBoundedScoringPathsEqsVars(possibleScores,upperLimit,lowerLimit,t,f,e)):
end:
###########################################################################################################################################
#EqualBoundedScoringPaths(possibleScores,upperLimit,lowerLimit,t): inputs possible scores and the upperLimit and lowerLimit for the y-value.
#Assumes that upperLimit>=0,lowerLimit<=0.
#possibleScores is a list or set of scores in the format [x-step,y-step].
#t is the variable in the g.f.
#Outputs the g.f. for the number of paths of total x-step n, that start at (0,0), end at (n,0), and stay within the limits.
#Uses EqualBoundedScoringPathsEqsVars.
#Try EqualBoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,1]},3,-2,t).
EqualBoundedScoringPaths:=proc(possibleScores,upperLimit,lowerLimit,t) local f,e:
	subs(AllEqualBoundedScoringPaths(possibleScores,upperLimit,lowerLimit,t,f,e),f[upperLimit,lowerLimit]):
end:
######################################################################################################
#EqualBoundedScoringPathsCoefficients(possibleScores,upperLimit,lowerLimit,N): inputs possible scores and the upperLimit and lowerLimit for the y-value.
#possibleScores is a list or set of scores in the format [x-step,y-step].
#Assumes that upperLimit>=0,lowerLimit<=0.
#Outputs the first N+1 coefficients of the Taylor series of the g.f. for the number of paths that start at (0,0), end at (n,0), and stay within the limits.
#Uses EqualBoundedScoringPaths.
#Try EqualBoundedScoringPathsCoefficients({[1,0],[1,-2],[2,-3],[0,1]},3,-2,9).
EqualBoundedScoringPathsCoefficients:=proc(possibleScores,upperLimit,lowerLimit,N) local t,f:
	f:=taylor(EqualBoundedScoringPaths(possibleScores,upperLimit,lowerLimit,t),t,N+1):
	seq(coeff(f,t,i),i=0..N):
end:
######################################################################################################
#SpecificEqualBoundedScoringPathsNumber(possibleScores,upperLimit,lowerLimit,n,startingLevel): inputs possible scores (x-step,y-step) and the upperLimit and lowerLimit for the y-value.
#Assumes that upperLimit>=0,lowerLimit<=0, and startingLevel is somewhere between these 2 numbers.
#possibleScores is a list or set of scores in the format [x,y].
#Outputs the number of paths that start at (0,startingLevel), end at (n,0), and stay within the bounds.
#Does this empirically.
#If there are BOTH [0,-] and [0,+] steps that sum to 0 within the bounds, then you will obtain infinite recursion.
#Try seq(SpecificEqualBoundedScoringPathsNumber({[1,0],[1,-2],[2,-3],[0,1]},3,-2,n1,0),n1=0..9).
SpecificEqualBoundedScoringPathsNumber:=proc(possibleScores,upperLimit,lowerLimit,n,startingLevel) local count,score: option remember:
	if n<0 then
		return(0):
	elif n=0 and startingLevel=0 then
		return(1):
	end:
	count:=0:
	for score in possibleScores do
		if score[2]>=lowerLimit and score[2]<=upperLimit then
			count:=count+SpecificEqualBoundedScoringPathsNumber(possibleScores,upperLimit-score[2],lowerLimit-score[2],n-score[1],startingLevel+score[2]):
		fi:
	od:
	count:
end:
######################################################################################################
#SpecificIrreducibleBoundedScoringPathsNumber(possibleScores,upperLimit,lowerLimit,n,startingLevel): inputs possible scores (x-step,y-step) and the upperLimit and lowerLimit for the y-value.
#Assumes that upperLimit>=0,lowerLimit<=0, and startingLevel is somewhere between these 2 numbers.
#possibleScores is a list or set of scores in the format [x,y].
#Outputs the number of paths that start at (0,startingLevel), end at (n,0) while never touching the x-axis previously, and stay within the bounds.
#Does this empirically.
#If there are BOTH [0,-] and [0,+] steps that sum to 0 within the bounds, then you MAY obtain infinite recursion.
#Try seq(SpecificIrreducibleBoundedScoringPathsNumber({[1,0],[1,-2],[2,-3],[0,1]},3,-2,n1,0),n1=0..9).
SpecificIrreducibleBoundedScoringPathsNumber:=proc(possibleScores,upperLimit,lowerLimit,n,startingLevel) local count,score: option remember:
	if startingLevel=0 then	#Must take at least 1 step.
		count:=0:
		for score in possibleScores do
			if score[2]<>0 and score[2]>=lowerLimit and score[2]<=upperLimit then
				count:=count+SpecificIrreducibleBoundedScoringPathsNumberHelp(possibleScores,upperLimit,lowerLimit,n-score[1],score[2]):
			fi:
		od:
		count:
	else
		SpecificIrreducibleBoundedScoringPathsNumberHelp(possibleScores,upperLimit,lowerLimit,n,startingLevel)
	fi:
end:
######################################################################################################
#SpecificIrreducibleBoundedScoringPathsNumberHelp(possibleScores,upperLimit,lowerLimit,n,startingLevel): inputs possible scores (x-step,y-step) and the upperLimit and lowerLimit for the y-value.
#Assumes that upperLimit>=0,lowerLimit<=0, and startingLevel is somewhere between these 2 numbers.
#possibleScores is a list or set of scores in the format [x,y].
#Outputs the number of paths that start at (0,startingLevel), end at (n,0) while never touching the x-axis previously, and stay within the bounds.
#Does this empirically.
#If there are BOTH [0,-] and [0,+] steps that sum to 0 within the bounds, then you MAY obtain infinite recursion.
#This is a helper function for SpecificIrreducibleBoundedScoringPathsNumber after ensuring that startingLevel<>0.
#Try seq(SpecificIrreducibleBoundedScoringPathsNumberHelp({[1,0],[1,-2],[2,-3],[0,1]},3,-2,n1,0),n1=0..9).
SpecificIrreducibleBoundedScoringPathsNumberHelp:=proc(possibleScores,upperLimit,lowerLimit,n,startingLevel) local count,score: option remember:
	if n<0 then
		return(0):
	elif startingLevel=0 then
		if n=0 then
			return(1):
		else
			return(0):
		fi:
	end:
	count:=0:
	for score in possibleScores do
		if score[2]+startingLevel>=lowerLimit and score[2]+startingLevel<=upperLimit then
			count:=count+SpecificIrreducibleBoundedScoringPathsNumberHelp(possibleScores,upperLimit,lowerLimit,n-score[1],startingLevel+score[2]):
		fi:
	od:
	count:
end:


###########################################################################################################################################
###########################################################################################################################################
#Semi-bounded
###########################################################################################################################################
###########################################################################################################################################
#SemiBoundedScoringPathsEqsVars(possibleScores,lowerLimit,t,k,f,g): inputs possible scores and lowerLimit for the y-value.
#possibleScores is a list or set of scores in the format [x-step,y-step].
#Assumes lowerLimit<=0.
#f,g,k are placeholders to keep track of the different g.f.s. g is for irreducible (only hits minimum at the respective endpoint) walks.
#f[a,b] is the g.f. for the number of walks from (0,a) to (n,b) that stay above the lowerLimit.
#g[a,b] is the g.f. for the number of irreducible walks from (0,a) to (n,b) that stay above the lowerLimit.
#f,g are the same g.f. as for EqualSemiBoundedScoringPathsEqsVars(possibleScores,lowerLimit,t,f,g).
#k[y] is the g.f. for walks above lowerLimit that begin at (0,y) and end anywhere.
#t is the variable in the g.f.s.
#Outputs a set of equations (and the used variables) relating the g.f.s held in f,g,k.
#Try SemiBoundedScoringPathsEqsVars({[1,0],[1,-2],[2,-3],[0,1]},-2,t,k,f,g).
SemiBoundedScoringPathsEqsVars:=proc(possibleScores,lowerLimit,t,k,f,g) local eqs,vars,stepsUp,stepsDown,stepsRight,yvalues,score:
	eqs,vars:=EqualSemiBoundedScoringPathsEqsVars(possibleScores,lowerLimit,t,f,g):
	stepsUp:={}:
	stepsDown:={}:
	stepsRight:={}:
	yvalues:={-lowerLimit}:
	for score in possibleScores do
		if score[2]<0 then
			stepsDown:=stepsDown union {score}:
		elif score[2]=0 then
			stepsRight:=stepsRight union {score}:
		else
			stepsUp:=stepsUp union {score}:
			yvalues:=yvalues union {score[2]-1}:
		fi:
	od:
	eqs:=eqs union {k[lowerLimit]=1+(g[lowerLimit,lowerLimit]+add(t^stepR[1],stepR in stepsRight))*k[lowerLimit]+add(t^stepU[1]*k[stepU[2]-1+lowerLimit],stepU in stepsUp),seq(k[y+lowerLimit]=(1+add(g[y-i+lowerLimit,lowerLimit],i=0..(y-1)))*k[lowerLimit],y in yvalues minus {0})}:
	vars:=vars union {k[lowerLimit],seq(k[y+lowerLimit],y in yvalues)}:
	eqs,vars:
end:
###########################################################################################################################################
#AllSemiBoundedScoringPaths(possibleScores,lowerLimit,t,k,f,g): inputs possible scores and lowerLimit for the y-value.
#possibleScores is a list or set of scores in the format [x-step,y-step].
#Assumes lowerLimit<=0.
#f,g,k are placeholders to keep track of the different g.f.s. g is for irreducible (only hits minimum at the respective endpoint) walks.
#f[a,b] is the g.f. for the number of walks from (0,a) to (n,b) that stay above the lowerLimit.
#g[a,b] is the g.f. for the number of irreducible walks from (0,a) to (n,b) that stay above the lowerLimit.
#f,g are the same g.f. as for EqualSemiBoundedScoringPathsEqsVars(possibleScores,lowerLimit,t,f,g).
#k[y] is the g.f. for walks above lowerLimit that begin at (0,y) and end anywhere.
#t is the variable in the g.f.s.
#Uses SemiBoundedScoringPathsEqsVars to find the relating equations.
#Then uses an "optimal" Groebner basis for each variable to find a polynomial of which the variable is a root.
#May take a while since it is computing a Groebner basis for each variable.
#Try AllSemiBoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,1]},-2,t,k,f,g).
AllSemiBoundedScoringPaths:=proc(possibleScores,lowerLimit,t,k,f,g) local eqs,vars,polynomials,var:
	eqs,vars:=SemiBoundedScoringPathsEqsVars(possibleScores,lowerLimit,t,k,f,g):
	eqs:={seq(lhs(eq)-rhs(eq),eq in eqs)}:
	polynomials:={}:
	for var in vars do
		if op(0,var)=k then
			polynomials:=polynomials union {FindProperRoot(possibleScores,Groebner[Basis](eqs,plex(op(vars minus {var}),var))[1],1,var,t,lowerLimit)}:
		elif op(0,var)=f then
			polynomials:=polynomials union {FindProperRoot(possibleScores,Groebner[Basis](eqs,plex(op(vars minus {var}),var))[1],2,var,t,lowerLimit)}:
		else	#op(0,var)=g
			polynomials:=polynomials union {FindProperRoot(possibleScores,Groebner[Basis](eqs,plex(op(vars minus {var}),var))[1],3,var,t,lowerLimit)}:
		fi:
	od:
	polynomials:
#	{seq(Groebner[Basis](eqs,plex(op(vars minus {var}),var))[1],var in vars)}:	#Does not reduce to minimal polynomials
end:
###########################################################################################################################################
#SpecificSemiBoundedScoringPaths(possibleScores,lowerLimit,t,k,var): The same as AllSemiBoundedScoringPaths except it only outputs the given g.f. held in var.
#inputs possible scores and lowerLimit for the y-value.
#possibleScores is a list or set of scores in the format [x-step,y-step].
#Assumes lowerLimit<=0.
#k[y] is the g.f. for walks above lowerLimit that begin at (0,y) and end anywhere.
#t is the variable in the g.f.s.
#var should be k[y]. It should be one of the variables returned in SemiBoundedScoringPathsEqsVars.
#If var is f[a,b] or g[a,b] it is better (and necessary) to use SpecificEqualSemiBoundedScoringPaths(possibleScores,lowerLimit,t,f,g,var).
#Uses SemiBoundedScoringPathsEqsVars to find the relating equations.
#Then uses an "optimal" Groebner basis to find a polynomial of which the variable specified in var is a root.
#Try SpecificSemiBoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,2]},-2,t,k,k[-1]).
SpecificSemiBoundedScoringPaths:=proc(possibleScores,lowerLimit,t,k,var) local eqs,vars,f,g:
	eqs,vars:=SemiBoundedScoringPathsEqsVars(possibleScores,lowerLimit,t,k,f,g):
	FindProperRoot(possibleScores,Groebner[Basis]({seq(lhs(eq)-rhs(eq),eq in eqs)},plex(op(vars minus {var}),var))[1],1,var,t,lowerLimit):
end:
###########################################################################################################################################
#SemiBoundedScoringPaths(possibleScores,lowerLimit,t,k): This is shorthand for computing SpecificSemiBoundedScoringPaths with var=k[0].
#inputs possible scores and lowerLimit for the y-value.
#possibleScores is a list or set of scores in the format [x-step,y-step].
#Assumes lowerLimit<=0.
#k is the g.f. for walks above lowerLimit that begin at the origin, end anywhere, and stay above y=lowerLimit.
#t is the variable in the g.f.
#Uses SemiBoundedScoringPathsEqsVars to find the relating equations.
#Then uses an "optimal" Groebner basis to find a polynomial with a root equal to the g.f. for the number of walks that start at (0,0), end at (n,0) and do not go below y=lowerLimit.
#Try SemiBoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,2]},-2,t,k).
SemiBoundedScoringPaths:=proc(possibleScores,lowerLimit,t,k) local g:
	subs(k[0]=k,SpecificSemiBoundedScoringPaths(possibleScores,lowerLimit,t,k,k[0])):
end:
######################################################################################################
#SemiBoundedScoringPathsSeries(possibleScores,lowerLimit,t,k,f,g,N): inputs possible scores and the lowerLimit for the y-value.
#possibleScores is a list or set of scores in the format [x-step,y-step].
#Assumes that lowerLimit<=0.
#Outputs the Taylor series (truncated at t^(N+1)) of the g.f.s found in SemiBoundedScoringPathsEqsVars.
#f,g,k are placeholders to keep track of the different g.f.s. g is for irreducible (only hits minimum at the respective endpoint) walks.
#f[a,b] is the g.f. for the number of walks from (0,a) to (n,b) that stay above the lowerLimit.
#g[a,b] is the g.f. for the number of irreducible walks from (0,a) to (n,b) that stay above the lowerLimit.
#f,g are the same g.f. as for EqualSemiBoundedScoringPathsEqsVars(possibleScores,lowerLimit,t,f,g).
#k[y] is the g.f. for walks above lowerLimit that begin at (0,y) and end anywhere.
#t is the variable in the g.f.s.
#Uses SemiBoundedScoringPathsEqsVars. Iterates to a fixed-point solution starting from {f[a,a]=1,f[a,b]=0,g[0,0]=0,g[a,b]=0,k[y]=1}.
#Try SemiBoundedScoringPathsSeries({[1,0],[1,-2],[2,-3],[0,1]},-2,t,k,f,g,9).
SemiBoundedScoringPathsSeries:=proc(possibleScores,lowerLimit,t,k,f,g,N) local eqs,vars,vect,var,vect2:
	eqs,vars:=SemiBoundedScoringPathsEqsVars(possibleScores,lowerLimit,t,k,f,g):
	vect:={}:
	for var in vars do
		if (op(0,var)=f and op(var)[1]=op(var)[2]) or op(0,var)=k then
			vect:=vect union {var=1}:
		else
			vect:=vect union {var=0}:
		fi:
	od:
	vect2:={}:
	while vect<>vect2 do
		vect2:=vect:
		vect:={seq(lhs(eq)=convert(taylor(subs(vect2,rhs(eq)),t,N+1),polynom),eq in eqs)}:
	od:
	vect:
end:
######################################################################################################
#SemiBoundedScoringPathsCoefficients(possibleScores,lowerLimit,N): inputs possible scores and the lowerLimit for the y-value.
#possibleScores is a list or set of scores in the format [x-step,y-step].
#Assumes that lowerLimit<=0.
#Outputs the first N+1 coefficients of the Taylor series of the g.f. for the number of paths that start at (0,0), stay above the lowerLimit, and end anywhere.
#Uses SemiBoundedScoringPaths.
#Try SemiBoundedScoringPathsCoefficients({[1,0],[1,-2],[2,-3],[0,1]},-2,9).
SemiBoundedScoringPathsCoefficients:=proc(possibleScores,lowerLimit,N) local p,t,k,k0,i,coe:#	local k0,t,k,f,g:
	p:=SemiBoundedScoringPaths(possibleScores,lowerLimit,t,k):
	k0:=1:
	coe:=coeff(coeff(p,k,1),t,0):
	if coe<>0 then
		p:=p-coe*k:
		for i from 1 to N do
#			k0:=k0+coeff(subs(k=k0,p),t,i)*t^i/(-coe)^2:
			k0:=convert(taylor(subs(k=k0,p),t,i+1),polynom)/(-coe):	#This is faster.
		od:
	else
		k0:=taylor(RootOf(subs(k=_Z,p)),t,N+1):	#This is fastest, but can sometimes give the wrong root. Though all of the "errors" have been in the unbounded case.
#The above is faster (if SemiBoundedScoringPaths is fast) but is not guaranteed to work. (Though it has worked in every case I have checked.)
#		k0:=subs(SemiBoundedScoringPathsSeries(possibleScores,lowerLimit,t,k,f,g,N),k[0]):
	fi:
	[seq(coeff(k0,t,i),i=0..N)]:
end:
######################################################################################################
#SpecificSemiBoundedScoringPathsNumber(possibleScores,lowerLimit,n,startingLevel): inputs possible scores (x-step,y-step) in the format {[x1,y1],[x2,y2],...}.
#Outputs the number of paths that start at (0,startingLevel), are of length n, and don't go below y=lowerLimit<=0.
#There CANNOT be (0,+) steps. Otherwise you get infinity as the answer.
#The answer is computed empirically.
#Try seq(SpecificSemiBoundedScoringPathsNumber({[1,0],[1,-2],[2,-3],[0,1]},-2,n1,0),n1=0..9).
SpecificSemiBoundedScoringPathsNumber:=proc(possibleScores,lowerLimit,n,startingLevel) local score,zeroScores:
	zeroScores:={}:
	for score in possibleScores do
		if score[1]=0 then
			if score[2]>=0 then
				return(infinity):
			else
				zeroScores:=zeroScores union {score[2]}:
			fi:
		fi:
	od:
	SpecificSemiBoundedScoringPathsNumberHelp(possibleScores,lowerLimit,n,startingLevel,zeroScores):
end:
###################################################
#A helper function to compute the number of walks after parsing in SpecificBoundedScoringPathsNumber.
SpecificSemiBoundedScoringPathsNumberHelp:=proc(possibleScores,lowerLimit,n,startingLevel,zeroScores) option remember:
	if n<0 or startingLevel<lowerLimit then
		return(0):
	elif n=0 then
		return(SpecificSemiBoundedScoringPathsNumberHelp2(zeroScores,startingLevel-lowerLimit)):
	fi:
	add(SpecificSemiBoundedScoringPathsNumberHelp(possibleScores,lowerLimit,n-score[1],startingLevel+score[2],zeroScores),score in possibleScores):
end:
###################################################
#Computes the number of ways to end at the current x value. All of these scores correspond to steps straight down.
SpecificSemiBoundedScoringPathsNumberHelp2:=proc(zeroScores,startingLevel) option remember:
	if startingLevel<0 then
		return(0):
	fi:
	1+add(SpecificSemiBoundedScoringPathsNumberHelp2(zeroScores,startingLevel+score),score in zeroScores):
end:

###########################################################################################################################################
#EqualSemiBoundedScoringPathsEqsVars(possibleScores,lowerLimit,t,f,g): inputs possible scores and lowerLimit for the y-value.
#possibleScores is a list or set of scores in the format [x-step,y-step].
#Assumes lowerLimit<=0.
#f,g are placeholders to keep track of the different g.f.s. g is for irreducible (only hits minimum at the respective endpoint) walks.
#f[a,b] is the g.f. for the number of walks from (0,a) to (n,b) that stay above the lowerLimit.
#g[a,b] is the g.f. for the number of irreducible walks from (0,a) to (n,b) that stay above the lowerLimit.
#t is the variable in the g.f.s.
#Outputs a set of equations (and the used variables) relating the g.f.s held in f and g.
#Try EqualSemiBoundedScoringPathsEqsVars({[1,0],[1,-2],[2,-3],[0,1]},-2,t,f,g).
EqualSemiBoundedScoringPathsEqsVars:=proc(possibleScores,lowerLimit,t,f,g) local stepsUp,stepsDown,stepsRight,score,eqs,walks,irrWalks,doneWalks,doneIrrWalks,w,iw:
	stepsUp:={}:
	stepsDown:={}:
	stepsRight:={}:
	for score in possibleScores do
		if score[2]<0 then
			stepsDown:=stepsDown union {score}:
		elif score[2]=0 then
			stepsRight:=stepsRight union {score}:
		else
			stepsUp:=stepsUp union {score}:
		fi:
	od:
	eqs:={f[0,0]=1+(g[0,0]+add(t^stepR[1],stepR in stepsRight))*f[0,0],g[0,0]=add(add(t^(stepU[1]+stepD[1])*f[stepU[2]-1,-stepD[2]-1],stepD in stepsDown),stepU in stepsUp)}:
	walks:={[-lowerLimit,-lowerLimit],seq(seq([stepU[2]-1,-stepD[2]-1],stepU in stepsUp),stepD in stepsDown)} minus {[0,0]}:	#Removing [0,0] since this was already done.
	irrWalks:={}:
	doneWalks:={[0,0]}:
	doneIrrWalks:={[0,0]}:
	while nops(walks)>0 or nops(irrWalks)>0 do	#This will eventually terminate since the coordinates are bounded by [max(stepsUp[2])-1,-max(stepsDown[2])-1].
		if nops(walks)>0 then
			w:=walks[1]:
			if w[1]>w[2] then	#Moving down
				eqs:=eqs union {f[op(w)]=add(g[w[1]-i,0]*f[0,w[2]-i],i=0..w[2])}:
				walks:=walks union {seq([0,w[2]-i],i=0..(w[2]-1))}:
				irrWalks:=irrWalks union {seq([w[1]-i,0],i=0..w[2])}:
			elif w[1]=w[2] then	#Staying flat
				eqs:=eqs union {f[op(w)]=add(g[w[1]-i,0]*f[0,w[2]-i],i=0..(w[2]-1))+f[0,0]}:
				walks:=walks union {seq([0,w[2]-i],i=0..(w[2]-1))}:
				irrWalks:=irrWalks union {seq([w[1]-i,0],i=0..(w[2]-1))}:
			else			#Moving up
				eqs:=eqs union {f[op(w)]=add(g[w[1]-i,0]*f[0,w[2]-i],i=0..(w[1]-1))+f[0,0]*g[0,w[2]-w[1]]}:
				walks:=walks union {seq([0,w[2]-i],i=0..(w[1]-1))}:
				irrWalks:=irrWalks union {seq([w[1]-i,0],i=0..(w[1]-1)),[0,w[2]-w[1]]}:
			fi:
			doneWalks:=doneWalks union {w}:
			walks:=walks minus doneWalks:
			irrWalks:=irrWalks minus doneIrrWalks:
		else	#nops(irrWalks)>0
			iw:=irrWalks[1]:
			if iw[1]>0 then
				eqs:=eqs union {g[op(iw)]=add(t^stepd[1]*f[iw[1]-1,-stepd[2]-1],stepd in stepsDown)}:
				walks:=walks union {seq([iw[1]-1,-stepd[2]-1],stepd in stepsDown)}:
			else	#iw[2]>0
				eqs:=eqs union {g[op(iw)]=add(t^stepu[1]*f[stepu[2]-1,iw[2]-1],stepu in stepsUp)}:
				walks:=walks union {seq([stepu[2]-1,iw[2]-1],stepu in stepsUp)}:
			fi:
			doneIrrWalks:=doneIrrWalks union {iw}:
			irrWalks:=irrWalks minus {iw}:
			walks:=walks minus doneWalks:
		fi:
	od:
#Currently f[a,b] is the g.f. for the number of walks from (0,a) to (n,b) that stay above y=0.
#And g[a,b] is the g.f. for the number of irreducible walks from (0,a) to (n,b).
#The following adjusts this so f[a,b] is the g.f. for the number of walks from (0,a) to (n,b) that stay above the lowerLimit. Similar for g[a,b].
	subs({seq(f[op(w)]=f[w[1]+lowerLimit,w[2]+lowerLimit],w in doneWalks),seq(g[op(iw)]=g[iw[1]+lowerLimit,iw[2]+lowerLimit],iw in doneIrrWalks)},eqs),{seq(f[w[1]+lowerLimit,w[2]+lowerLimit],w in doneWalks),seq(g[iw[1]+lowerLimit,iw[2]+lowerLimit],iw in doneIrrWalks)}:
end:
###########################################################################################################################################
#AllEqualSemiBoundedScoringPaths(possibleScores,lowerLimit,t,f,g): inputs possible scores and lowerLimit for the y-value.
#possibleScores is a list or set of scores in the format [x-step,y-step].
#Assumes lowerLimit<=0.
#f,g are placeholders to keep track of the different g.f.s. g is for irreducible (only hits minimum at the respective endpoint) walks.
#f[a,b] is the g.f. for the number of walks from (0,a) to (n,b) that stay above the lowerLimit.
#g[a,b] is the g.f. for the number of irreducible walks from (0,a) to (n,b) that stay above the lowerLimit.
#t is the variable in the g.f.s.
#Uses EqualSemiBoundedScoringPathsEqsVars to find the relating equations.
#Then uses an "optimal" Groebner basis for each variable to find a polynomial of which the variable is a root.
#May take a while since it is computing a Groebner basis for each variable.
#Try AllEqualSemiBoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,1]},-2,t,f,g).
AllEqualSemiBoundedScoringPaths:=proc(possibleScores,lowerLimit,t,f,g) local eqs,vars,polynomials,var:
	eqs,vars:=EqualSemiBoundedScoringPathsEqsVars(possibleScores,lowerLimit,t,f,g):
	eqs:={seq(lhs(eq)-rhs(eq),eq in eqs)}:
	polynomials:={}:
	for var in vars do
		if op(0,var)=f then
			polynomials:=polynomials union {FindProperRoot(possibleScores,Groebner[Basis](eqs,plex(op(vars minus {var}),var))[1],2,var,t,lowerLimit)}:
		else	#op(0,var)=g
			polynomials:=polynomials union {FindProperRoot(possibleScores,Groebner[Basis](eqs,plex(op(vars minus {var}),var))[1],3,var,t,lowerLimit)}:
		fi:
	od:
	polynomials:
#	{seq(Groebner[Basis](eqs,plex(op(vars minus {var}),var))[1],var in vars)}:	#Does not reduce to minimal polynomials
end:
###########################################################################################################################################
#SpecificEqualSemiBoundedScoringPaths(possibleScores,lowerLimit,t,f,g,var): The same as AllEqualSemiBoundedScoringPaths except it only outputs the given g.f. held in var.
#inputs possible scores and lowerLimit for the y-value.
#possibleScores is a list or set of scores in the format [x-step,y-step].
#Assumes lowerLimit<=0.
#f,g are placeholders to keep track of the different g.f.s. g is for irreducible (only hits minimum at the respective endpoint) walks.
#f[a,b] is the g.f. for the number of walks from (0,a) to (n,b) that stay above the lowerLimit.
#g[a,b] is the g.f. for the number of irreducible walks from (0,a) to (n,b) that stay above the lowerLimit.
#t is the variable in the g.f.s.
#var is either f[a,b] or g[a,b]. It should be one of the variables returned in EqualSemiBoundedScoringPathsEqsVars.
#Uses EqualSemiBoundedScoringPathsEqsVars to find the relating equations.
#Then uses an "optimal" Groebner basis to find a polynomial of which the variable specified in var is a root.
#Try SpecificEqualSemiBoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,1]},-2,t,f,g,f[-1,-1]).
SpecificEqualSemiBoundedScoringPaths:=proc(possibleScores,lowerLimit,t,f,g,var) local eqs,vars:
	eqs,vars:=EqualSemiBoundedScoringPathsEqsVars(possibleScores,lowerLimit,t,f,g):
	if op(0,var)=f then
		FindProperRoot(possibleScores,Groebner[Basis]({seq(lhs(eq)-rhs(eq),eq in eqs)},plex(op(vars minus {var}),var))[1],2,var,t,lowerLimit):
	else	#op(0,var)=g
		FindProperRoot(possibleScores,Groebner[Basis]({seq(lhs(eq)-rhs(eq),eq in eqs)},plex(op(vars minus {var}),var))[1],3,var,t,lowerLimit):
	fi:
end:
###########################################################################################################################################
#EqualSemiBoundedScoringPaths(possibleScores,lowerLimit,t,f): This is shorthand for computing SpecificEqualSemiBoundedScoringPaths with var=f[0,0].
#inputs possible scores and lowerLimit for the y-value.
#possibleScores is a list or set of scores in the format [x-step,y-step].
#Assumes lowerLimit<=0.
#f is the g.f. for walks that stay above lowerLimit, begin at the origin and end at (n,0).
#t is the variable in the g.f.
#Uses EqualSemiBoundedScoringPathsEqsVars to find the relating equations.
#Then uses an "optimal" Groebner basis to find a polynomial with a root equal to the g.f. for the number of walks that start at (0,0), end at (n,0) and do not go below y=lowerLimit.
#Try EqualSemiBoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,1]},-2,t,f).
EqualSemiBoundedScoringPaths:=proc(possibleScores,lowerLimit,t,f) local g:
	subs(f[0,0]=f,SpecificEqualSemiBoundedScoringPaths(possibleScores,lowerLimit,t,f,g,f[0,0])):
end:
######################################################################################################
#EqualSemiBoundedScoringPathsSeries(possibleScores,lowerLimit,t,f,g,N): inputs possible scores and the lowerLimit for the y-value.
#possibleScores is a list or set of scores in the format [x-step,y-step].
#Assumes that lowerLimit<=0.
#Outputs the Taylor series (truncated at t^(N+1)) of the g.f.s found in EqualSemiBoundedScoringPathsEqsVars.
#f[a,b] is the g.f. for the number of walks from (0,a) to (n,b) that stay above the lowerLimit.
#g[a,b] is the g.f. for the number of irreducible walks from (0,a) to (n,b) that stay above the lowerLimit.
#Uses EqualSemiBoundedScoringPathsEqsVars. Iterates to a fixed-point solution starting from {f[a,a]=1,f[a,b]=0,g[0,0]=0,g[a,b]=0}.
#Try EqualSemiBoundedScoringPathsSeries({[1,0],[1,-2],[2,-3],[0,1]},-2,t,f,g,9).
EqualSemiBoundedScoringPathsSeries:=proc(possibleScores,lowerLimit,t,f,g,N) local eqs,vars,vect,var,vect2:
	eqs,vars:=EqualSemiBoundedScoringPathsEqsVars(possibleScores,lowerLimit,t,f,g):
	vect:={}:
	for var in vars do
		if op(0,var)=f and op(var)[1]=op(var)[2] then
			vect:=vect union {var=1}:
		else
			vect:=vect union {var=0}:
		fi:
	od:
	vect2:={}:
	while vect<>vect2 do
		vect2:=vect:
		vect:={seq(lhs(eq)=convert(taylor(subs(vect2,rhs(eq)),t,N+1),polynom),eq in eqs)}:
	od:
	vect:
end:
######################################################################################################
#EqualSemiBoundedScoringPathsCoefficients(possibleScores,lowerLimit,N): inputs possible scores and the lowerLimit for the y-value.
#possibleScores is a list or set of scores in the format [x-step,y-step].
#Assumes that lowerLimit<=0.
#Outputs the first N+1 coefficients of the Taylor series of the g.f. for the number of paths that start at (0,0), end at (n,0), and stay above the lowerLimit.
#Uses EqualSemiBoundedScoringPaths.
#Try EqualSemiBoundedScoringPathsCoefficients({[1,0],[1,-2],[2,-3],[0,1]},-2,9).
EqualSemiBoundedScoringPathsCoefficients:=proc(possibleScores,lowerLimit,N) local p,t,f,f0,i,coe:	#local f0,t,f,g:
	p:=EqualSemiBoundedScoringPaths(possibleScores,lowerLimit,t,f):
	f0:=1:
	coe:=coeff(coeff(p,f,1),t,0):
	if coe<>0 then
		p:=p-coe*f:
		for i from 1 to N do
#			f0:=f0+coeff(subs(f=f0,p),t,i)*t^i/(-coe)^2:
			f0:=convert(taylor(subs(f=f0,p),t,i+1),polynom)/(-coe):	#This is faster.
		od:
	else
		f0:=taylor(RootOf(subs(f=_Z,p)),t,N+1):	#This is fastest, but can sometimes give the wrong root. Though all of the "errors" have been in the unbounded case.
#The above is faster (if EqualSemiBoundedScoringPaths is fast) but is not guaranteed to work. (Though it has worked in every case I have checked.)
#		f0:=subs(EqualSemiBoundedScoringPathsSeries(possibleScores,lowerLimit,t,f,g,N),f[0,0]):
	fi:
	[seq(coeff(f0,t,i),i=0..N)]:
end:
######################################################################################################
#SpecificEqualSemiBoundedScoringPathsNumber(possibleScores,lowerLimit,n,startingLevel): inputs possible scores (x-step,y-step) in the format {[x1,y1],[x2,y2],...}.
#Outputs the number of paths that start at (0,startingLevel), end at (n,0), and don't go below y=lowerLimit<=0.
#There should NOT be BOTH (0,+) and (0,-) steps. Otherwise you get infinity as the answer.
#The answer is computed empirically.
#Try seq(SpecificEqualSemiBoundedScoringPathsNumber({[1,0],[1,-2],[2,-3],[0,1]},-2,n1,0),n1=0..9).
SpecificEqualSemiBoundedScoringPathsNumber:=proc(possibleScores,lowerLimit,n,startingLevel) local zeroPstep,zeroNstep,score,maxDescentRate:
	zeroPstep:=false:
	zeroNstep:=false:
	for score in possibleScores do
		if score[1]=0 then
			if score[2]>0 then
				zeroPstep:=true:
			else	#score[2]<0
				zeroNstep:=true:
			fi:
		fi:
	od:
	if zeroPstep then
		if zeroNstep then
			return(infinity):
		else
			maxDescentRate:=0:
			for score in possibleScores do
				if score[2]<0 then
					maxDescentRate:=min(maxDescentRate,score[2]/score[1]):
				fi:
			od:
			return(SpecificEqualSemiBoundedScoringPathsNumberP(possibleScores,lowerLimit,n,startingLevel,maxDescentRate)):
		fi:
	fi:
	SpecificEqualSemiBoundedScoringPathsNumberN(possibleScores,lowerLimit,n,startingLevel):
end:
#################################################
#A helper function that works if there are no (0,+) steps.
SpecificEqualSemiBoundedScoringPathsNumberN:=proc(possibleScores,lowerLimit,n,startingLevel) option remember:
	if n<0 or startingLevel<lowerLimit then
		return(0):
	elif n=0 and startingLevel=0 then
		return(1):
	fi:
	add(SpecificEqualSemiBoundedScoringPathsNumberN(possibleScores,lowerLimit,n-score[1],startingLevel+score[2]),score in possibleScores):
end:
#################################################
#A helper function that works if there are no (0,-) steps.
#Inputting maxDescentRate just saves some computation time. It is determined by possibleScores.
SpecificEqualSemiBoundedScoringPathsNumberP:=proc(possibleScores,lowerLimit,n,startingLevel,maxDescentRate) option remember:
	if n<0 or startingLevel<lowerLimit or startingLevel+n*maxDescentRate>0 then
		return(0):
	elif n=0 and startingLevel=0 then
		return(1):
	fi:
	add(SpecificEqualSemiBoundedScoringPathsNumberP(possibleScores,lowerLimit,n-score[1],startingLevel+score[2],maxDescentRate),score in possibleScores):
end:
######################################################################################################
#SpecificIrreducibleSemiBoundedScoringPathsNumber(possibleScores,n,startingLevel): inputs possible scores (x-step,y-step) in the format {[x1,y1],[x2,y2],...}.
#Outputs the number of paths that start at (0,startingLevel), end at (n,0), and don't go below y=min(startingLevel,0).
#There should NOT be BOTH (0,+) and (0,-) steps. Otherwise you get infinity as the answer.
#The answer is computed empirically.
#Try seq(SpecificIrreducibleSemiBoundedScoringPathsNumber({[1,0],[1,-2],[2,-3],[0,1]},-2,n1,0),n1=0..9).
SpecificIrreducibleSemiBoundedScoringPathsNumber:=proc(possibleScores,n,startingLevel) local zeroPstep,zeroNstep,score,maxDescentRate,count:
	zeroPstep:=false:
	zeroNstep:=false:
	for score in possibleScores do
		if score[1]=0 then
			if score[2]>0 then
				zeroPstep:=true:
			else	#score[2]<0
				zeroNstep:=true:
			fi:
		fi:
	od:
	if zeroPstep then
		if zeroNstep then
			return(infinity):
		else
			maxDescentRate:=0:
			for score in possibleScores do
				if score[2]<0 then
					maxDescentRate:=min(maxDescentRate,score[2]/score[1]):
				fi:
			od:
			if startingLevel<=0 then
				count:=0:
				for score in possibleScores do
					if score[2]>0 then	#There must be a first step. And it must be positive.
						count:=count+SpecificIrreducibleSemiBoundedScoringPathsNumberP(possibleScores,min(startingLevel+1,0),n-score[1],startingLevel+score[2],maxDescentRate):
					fi:
				od:
				return(count):
			else
				return(SpecificIrreducibleSemiBoundedScoringPathsNumberP(possibleScores,0,n,startingLevel,maxDescentRate)):
			fi:
		fi:
	fi:
	if startingLevel<=0 then
		count:=0:
		for score in possibleScores do
			if score[2]>0 then	#There must be a first step. And it must be positive.
				count:=count+SpecificIrreducibleSemiBoundedScoringPathsNumberN(possibleScores,min(startingLevel+1,0),n-score[1],startingLevel+score[2]):
			fi:
		od:
		count
	else
		SpecificIrreducibleSemiBoundedScoringPathsNumberN(possibleScores,0,n,startingLevel):
	fi:
end:
#################################################
#A helper function that works if there are no (0,+) steps.
SpecificIrreducibleSemiBoundedScoringPathsNumberN:=proc(possibleScores,lowerLimit,n,startingLevel) option remember:
	if n<0 or startingLevel<lowerLimit then
		return(0):
	elif n=0 and startingLevel=0 then
		return(1):
	fi:
	add(SpecificIrreducibleSemiBoundedScoringPathsNumberN(possibleScores,lowerLimit,n-score[1],startingLevel+score[2]),score in possibleScores):
end:
#################################################
#A helper function that works if there are no (0,-) steps.
#Inputting maxDescentRate just saves some computation time. It is determined by possibleScores.
SpecificIrreducibleSemiBoundedScoringPathsNumberP:=proc(possibleScores,lowerLimit,n,startingLevel,maxDescentRate) option remember:
	if n<0 or startingLevel<lowerLimit or startingLevel+n*maxDescentRate>0 then
		return(0):
	elif n=0 and startingLevel=0 then
		return(1):
	fi:
	add(SpecificIrreducibleSemiBoundedScoringPathsNumberP(possibleScores,lowerLimit,n-score[1],startingLevel+score[2],maxDescentRate),score in possibleScores):
end:


###########################################################################################################################################
###########################################################################################################################################
#Unbounded
###########################################################################################################################################
###########################################################################################################################################
#UnBoundedScoringPaths(possibleScores,t): inputs possible scores and a variable t.
#possibleScores is a list or set of scores in the format [x,y]. All scores should have x-step>0.
#t is the variable in the generating function.
#Outputs the generating function for the number of paths of total x-step, n, that start at (0,0) and go anywhere.
#Try UnBoundedScoringPaths({[1,1],[1,-1],[1,2],[2,-3]},t).
UnBoundedScoringPaths:=proc(possibleScores,t):
	1/(1-add(t^score[1],score in possibleScores)):
end:
######################################################################################################
#UnBoundedScoringPathsCoefficients(possibleScores,N): Outputs the first N+1 coefficients of the Taylor series of the g.f. for the number of paths of x-step n.
#Uses UnBoundedScoringPaths.
#Try UnBoundedScoringPathsCoefficients({[1,1],[1,-1],[1,2],[2,-3]},9).
UnBoundedScoringPathsCoefficients:=proc(possibleScores,N) local t,f:
	f:=taylor(UnBoundedScoringPaths(possibleScores,t),t,N+1):
	seq(coeff(f,t,i),i=0..N):
end:
######################################################################################################
#UnBoundedScoringPathsNumber(possibleScores,n): inputs possible scores (x-step,y-step)
#Outputs the number of paths of total x-step, n, that start at (0,0) and go anywhere.
#The answer is computed empirically.
#Try seq(UnBoundedScoringPathsNumber({[1,1],[1,-1],[1,2],[2,-3]},n1),n1=0..9).
UnBoundedScoringPathsNumber:=proc(possibleScores,n):
	if member(0,{seq(score[1],score in possibleScores)}) then
		return(infinity):
	fi:
	UnBoundedScoringPathsNumber1(possibleScores,n):
end:
#################################################
#A helper function after ensuring no infinite recursion.
UnBoundedScoringPathsNumber1:=proc(possibleScores,n) option remember:
	if n<0 then
		return(0):
	elif n=0 then
		return(1):
	fi:
	add(UnBoundedScoringPathsNumber1(possibleScores,n-score[1]),score in possibleScores):
end:

###########################################################################################################################################
#EqualUnBoundedScoringPathsEqsVars(possibleScores,t,GF,h,f,g): inputs possible scores (x-step,y-step), and placeholder variables.
#possibleScores is a list or set of scores in the format [x,y].
#h,f,g are placeholder generating functions. f,g are the same g.f. as for EqualSemiBoundedScoringPathsEqsVars(possibleScores,0,t,f,g).
#GF is the g.f. for the number of walks that begin at (0,0) and end at (n,0).
#h[i] is the g.f. for walks that start at (i,0) and return to the x-axis only at the end.
#f[a,b] is the g.f. for the number of walks from (0,a) to (n,b) that stay above the x-axis.
#g[a,b] is the g.f. for the number of irreducible (only hits min(a,b) at the respective endpoint) walks from (0,a) to (n,b) that stay above the x-axis.
#t is the variable in the generating functions.
#Outputs many generating functions and equations describing their relations.
#Alternately, a more "optimal" solution for gf is found after adding the equation in the function EqualUnBoundedScoringPaths.
#Try EqualUnBoundedScoringPathsEqsVars({[1,0],[1,-2],[2,-3],[0,1]},t,GF,h,f,g).
EqualUnBoundedScoringPathsEqsVars:=proc(possibleScores,t,GF,h,f,g) local stepsUp,stepsDown,stepsRight,score,maxJump,eqs,irrEqs,hEqs,walks,w:
	stepsUp:={}:
	stepsDown:={}:
	stepsRight:={}:
	for score in possibleScores do
		if score[2]<0 then
			stepsDown:=stepsDown union {score}:
		elif score[2]=0 then
			stepsRight:=stepsRight union {score}:
		else
			stepsUp:=stepsUp union {score}:
		fi:
	od:
	maxJump:=max(seq(abs(score[2]),score in possibleScores))-1:

#The h[i] equations written out on multiple lines for ease of reading code.
#	h[0]=2*g[0,0]+add(add(g[ls[2]-i,0]*t^ls[1]*h[i],i=1..(ls[2]-1)),ls in stepsUp)+add(add(g[0,i]*t^ls[1]*h[i+ls[2]],i=1..(-ls[2]-1)),ls in stepsDown)
#	h[0] is partially why irreducible walks do not include steps to the right. If g[0,0] did, then I would have to write h[0]=2*g[0,0]-stepsRight, which is less elegant.
#	seq(h[j]=g[j,0]+add(add(f[j-1,i-1]*t^ls[1]*h[i+ls[2]],i=1..(-ls[2]-1)),ls in stepsDown),j=1..maxJump)
#	seq(h[-j]=g[0,j]+add(add(f[ls[2]-i-1,j-1]*t^ls[1]*h[i],i=1..(ls[2]-1)),ls in stepsUp),j=1..maxJump)
	hEqs:={h[0]=2*g[0,0]+add(add(g[ls[2]-i,0]*t^ls[1]*h[i],i=1..(ls[2]-1)),ls in stepsUp)+add(add(g[0,i]*t^ls[1]*h[i+ls[2]],i=1..(-ls[2]-1)),ls in stepsDown),seq(h[j]=g[j,0]+add(add(f[j-1,i-1]*t^ls[1]*h[i+ls[2]],i=1..(-ls[2]-1)),ls in stepsDown),j=1..maxJump),seq(h[-j]=g[0,j]+add(add(f[ls[2]-i-1,j-1]*t^ls[1]*h[i],i=1..(ls[2]-1)),ls in stepsUp),j=1..maxJump)}:

#We cannot simply call EqualSemiBoundedScoringPathsEqsVars(possibleScores,0,t,f,g) for the remaining equations since there may be some walks missing.
#The irreducible walk equations:
#	g[0,0]=add(add(t^(stepU[1]+stepD[1])*f[stepU[2]-1,-stepD[2]-1],stepD in stepsDown),stepU in stepsUp)
#	seq(g[0,j]=add(t^stepU[1]*f[stepU[2]-1,j-1],stepU in stepsUp),j=1..maxJump)
#	seq(g[j,0]=add(t^stepD[1]*f[j-1,-stepD[2]-1],stepD in stepsDown),j=1..maxJump)
	irrEqs:={g[0,0]=add(add(t^(stepU[1]+stepD[1])*f[stepU[2]-1,-stepD[2]-1],stepD in stepsDown),stepU in stepsUp),seq(g[0,j]=add(t^stepU[1]*f[stepU[2]-1,j-1],stepU in stepsUp),j=1..maxJump),seq(g[j,0]=add(t^stepD[1]*f[j-1,-stepD[2]-1],stepD in stepsDown),j=1..maxJump)}:

	eqs:={f[0,0]=1+(g[0,0]+add(t^stepR[1],stepR in stepsRight))*f[0,0]}:
	walks:={seq(seq([j,i],i=0..(max(seq(-stepD[2],stepD in stepsDown))-1)),j=0..(maxJump-1)),seq(seq([i,j],i=0..(max(seq(stepU[2],stepU in stepsUp))-1)),j=0..(maxJump-1)),
		seq(seq([stepU[2]-1,-stepD[2]-1],stepD in stepsDown),stepU in stepsUp)} minus {[0,0]}: #[0,0] was already done.
	for w in walks do
		if w[1]>w[2] then	#Moving down
			eqs:=eqs union {f[op(w)]=add(g[w[1]-i,0]*f[0,w[2]-i],i=0..w[2])}:
		elif w[1]=w[2] then	#Staying flat
			eqs:=eqs union {f[op(w)]=add(g[w[1]-i,0]*f[0,w[2]-i],i=0..(w[2]-1))+f[0,0]}:
		else			#Moving up
			eqs:=eqs union {f[op(w)]=add(g[w[1]-i,0]*f[0,w[2]-i],i=0..(w[1]-1))+f[0,0]*g[0,w[2]-w[1]]}:
		fi:
	od:

	eqs union irrEqs union hEqs union {GF=1+(h[0]+add(t^sr[1],sr in stepsRight))*GF}, {GF,f[0,0],seq(f[op(w)],w in walks),seq(g[0,j],j=0..maxJump),seq(g[j,0],j=1..maxJump),seq(h[j],j=-maxJump..maxJump)}:
end:
###########################################################################################################################################
#AllEqualUnBoundedScoringPaths(possibleScores,t,GF,h,f,g): inputs possible scores (x-step,y-step), and placeholder variables.
#possibleScores is a list or set of scores in the format [x,y].
#h,f,g are placeholder generating functions. f,g are the same g.f. as for EqualSemiBoundedScoringPathsEqsVars(possibleScores,0,t,f,g).
#GF is the g.f. for the number of walks that begin at (0,0) and end at (n,0).
#h[i] is the g.f. for walks that start at (i,0) and return to the x-axis only at the end.
#f[a,b] is the g.f. for the number of walks from (0,a) to (n,b) that stay above the x-axis.
#g[a,b] is the g.f. for the number of irreducible (only hits min(a,b) at the respective endpoint) walks from (0,a) to (n,b) that stay above the x-axis.
#t is the variable in the generating functions.
#Uses EqualUnBoundedScoringPathsEqsVars to find the relating equations.
#Then uses an "optimal" Groebner basis for each variable to find a polynomial of which the variable is a root.
#May take a while since it is computing a Groebner basis for each variable.
#Try AllEqualUnBoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,1]},t,GF,h,f,g).
AllEqualUnBoundedScoringPaths:=proc(possibleScores,t,GF,h,f,g) local eqs,vars,polynomials,var:
	eqs,vars:=EqualUnBoundedScoringPathsEqsVars(possibleScores,t,GF,h,f,g):
	eqs:={seq(lhs(eq)-rhs(eq),eq in eqs)}:
	polynomials:={}:
	for var in vars do
		if var=GF then
			polynomials:=polynomials union {FindProperRoot(possibleScores,Groebner[Basis](eqs,plex(op(vars minus {var}),var))[1],4,var,t)}:
		elif op(0,var)=h then
			polynomials:=polynomials union {FindProperRoot(possibleScores,Groebner[Basis](eqs,plex(op(vars minus {var}),var))[1],5,var,t)}:
		elif op(0,var)=f then
			polynomials:=polynomials union {FindProperRoot(possibleScores,Groebner[Basis](eqs,plex(op(vars minus {var}),var))[1],2,var,t,0)}:
		else	#op(0,var)=g
			polynomials:=polynomials union {FindProperRoot(possibleScores,Groebner[Basis](eqs,plex(op(vars minus {var}),var))[1],3,var,t,0)}:
		fi:
	od:
	polynomials:
#	{seq(Groebner[Basis](eqs,plex(op(vars minus {var}),var))[1],var in vars)}:	#Does not reduce to minimal polynomials.
end:
###########################################################################################################################################
#SpecificUnBoundedScoringPaths(possibleScores,lowerLimit,t,GF,h,var): The same as AllEqualUnBoundedScoringPaths except it only outputs the given g.f. held in var.
#inputs possible scores (x-step,y-step), and placeholder variables.
#possibleScores is a list or set of scores in the format [x,y].
#GF is the g.f. for the number of walks that begin at (0,0) and end at (n,0).
#h[i] is the g.f. for walks that start at (i,0) and return to the x-axis only at the end.
#t is the variable in the g.f.s.
#var is either GF or h[i]. It should be one of the variables returned in EqualUnBoundedScoringPathsEqsVars.
#If var is f[a,b] or g[a,b] it is better (and necessary) to use SpecificEqualSemiBoundedScoringPaths(possibleScores,0,t,f,g,var).
#Uses EqualUnBoundedScoringPathsEqsVars to find the relating equations.
#Then uses an "optimal" Groebner basis to find a polynomial of which var is a root.
#Try SpecificUnBoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,1]},t,GF,h,h[1]).
SpecificUnBoundedScoringPaths:=proc(possibleScores,t,GF,h,var) local eqs,vars,f,g:
	eqs,vars:=EqualUnBoundedScoringPathsEqsVars(possibleScores,t,GF,h,f,g):
	if var=GF then
		FindProperRoot(possibleScores,Groebner[Basis]({seq(lhs(eq)-rhs(eq),eq in eqs)},plex(op(vars minus {var}),var))[1],4,var,t):
	else	#op(0,var)=h
		FindProperRoot(possibleScores,Groebner[Basis]({seq(lhs(eq)-rhs(eq),eq in eqs)},plex(op(vars minus {var}),var))[1],5,var,t):
	fi:
end:
###########################################################################################################################################
#EqualUnBoundedScoringPaths(possibleScores,t,GF): inputs possible scores and a variable t.
#possibleScores is a list or set of scores in the format [x-step,y-step].
#GF is the g.f. for the number of walks that begin at (0,0) and end at (n,0).
#t is the variable in the g.f.
#Uses EqualUnBoundedScoringPathsEqsVars to find the relating equations.
#Then uses an "optimal" Groebner basis to find a polynomial with a root equal to GF.
#Try EqualUnBoundedScoringPaths({[1,0],[1,-2],[2,-3],[0,1]},t,GF).
EqualUnBoundedScoringPaths:=proc(possibleScores,t,GF) local h:
	SpecificUnBoundedScoringPaths(possibleScores,t,GF,h,GF)
end:
######################################################################################################
#EqualUnBoundedScoringPathsSeries(possibleScores,t,GF,h,f,g,N): inputs possible scores which is a list or set of scores in the format [x-step,y-step].
#Outputs the Taylor series (truncated at t^(N+1)) of the g.f.s found in UnBoundedScoringPathsEqsVars.
#h,f,g are placeholder generating functions. f,g are the same g.f. as for EqualSemiBoundedScoringPathsEqsVars(possibleScores,0,t,f,g).
#GF is the g.f. for the number of walks that begin at (0,0) and end at (n,0).
#h[i] is the g.f. for walks that start at (i,0) and return to the x-axis only at the end.
#f[a,b] is the g.f. for the number of walks from (0,a) to (n,b) that stay above the x-axis.
#g[a,b] is the g.f. for the number of irreducible (only hits min(a,b) at the respective endpoint) walks from (0,a) to (n,b) that stay above the x-axis.
#Uses UnBoundedScoringPathsEqsVars. Iterates to a fixed-point solution starting from {GF=1,f[a,a]=1,f[a,b]=0,h[i]=0,g[0,0]=0,g[a,b]=0}.
#Try EqualUnBoundedScoringPathsSeries({[1,0],[1,-2],[2,-3],[0,1]},t,GF,h,f,g,9).
EqualUnBoundedScoringPathsSeries:=proc(possibleScores,t,GF,h,f,g,N) local eqs,vars,vect,var,vect2:
	eqs,vars:=EqualUnBoundedScoringPathsEqsVars(possibleScores,t,GF,h,f,g):
	vect:={}:
	for var in vars do
		if var=GF or (op(0,var)=f and op(var)[1]=op(var)[2]) then
			vect:=vect union {var=1}:
		else
			vect:=vect union {var=0}:
		fi:
	od:
	vect2:={}:
	while vect<>vect2 do
		vect2:=vect:
		vect:={seq(lhs(eq)=convert(taylor(subs(vect2,rhs(eq)),t,N+1),polynom),eq in eqs)}:
	od:
	vect:
end:
######################################################################################################
#EqualUnBoundedScoringPathsCoefficients(possibleScores,N): inputs possible scores which is a list or set of scores in the format [x-step,y-step].
#Outputs the first N+1 coefficients of the Taylor series of the g.f. for the number of paths that start at (0,0) and end at (n,0).
#Uses EqualUnBoundedScoringPathsSeries.	# or EqualUnBoundedScoringPaths.
#Try EqualUnBoundedScoringPathsCoefficients({[1,0],[1,-2],[2,-3],[0,1]},9).
EqualUnBoundedScoringPathsCoefficients:=proc(possibleScores,N) local GF0,t,GF,h,f,g:	#local p,t,GF,GF0,i,coe:
#	p:=EqualUnBoundedScoringPaths(possibleScores,t,GF):
#	GF0:=1:
#	coe:=coeff(coeff(p,GF,1),t,0):
#	if coe<>0 then
#		p:=p-coe*GF:
#		for i from 1 to N do
##			GF0:=GF0+coeff(subs(GF=GF0,p),t,i)*t^i/(-coe)^2:
#			GF0:=convert(taylor(subs(GF=GF0,p),t,i+1),polynom)/(-coe):	#This is faster.
#		od:
#	else
##		GF0:=taylor(RootOf(subs(GF=_Z,p)),t,N+1):	#This is fastest, but will often give the wrong root.
#The above is faster (if EqualUnBoundedScoringPaths is fast) but is not guaranteed to work.
		GF0:=subs(EqualUnBoundedScoringPathsSeries(possibleScores,t,GF,h,f,g,N),GF):
#	fi:
	[seq(coeff(GF0,t,i),i=0..N)]:
end:
######################################################################################################
#SpecificUnBoundedScoringPathsNumber(possibleScores,n,startingLevel): inputs possible scores (x-step,y-step) in the format {[x1,y1],[x2,y2],...}.
#Outputs the number of paths that start at (0,startingLevel) and end at (n,0).
#There should NOT be BOTH (0,+) and (0,-) steps. Otherwise you get infinity as the answer.
#The answer is computed empirically.
#Try seq(SpecificUnBoundedScoringPathsNumber({[1,0],[1,-2],[2,-3],[0,1]},n1,0),n1=0..9).
SpecificUnBoundedScoringPathsNumber:=proc(possibleScores,n,startingLevel) local zeroPstep,zeroNstep,score,maxRate:
	zeroPstep:=false:
	zeroNstep:=false:
	for score in possibleScores do
		if score[1]=0 then
			if score[2]>0 then
				zeroPstep:=true:
			else	#score[2]<0
				zeroNstep:=true:
			fi:
		fi:
	od:
	if zeroPstep then
		if zeroNstep then
			return(infinity):
		else
			maxRate:=0:
			for score in possibleScores do
				if score[2]<0 then	#Then score[1]<>0 since zeroNstep=false
					maxRate:=min(maxRate,score[2]/score[1]):
				fi:
			od:
			return(SpecificUnBoundedScoringPathsNumberP(possibleScores,n,startingLevel,maxRate)):
		fi:
	fi:
	maxRate:=0:
	for score in possibleScores do
		if score[2]>0 then	#Then score[1]<>0 since zeroPstep=false
			maxRate:=max(maxRate,score[2]/score[1]):
		fi:
	od:
	SpecificUnBoundedScoringPathsNumberN(possibleScores,n,startingLevel,maxRate):
end:
#################################################
#A helper function that works if there are no (0,+) steps.
#Inputting maxRate just saves some computation time. It is determined by possibleScores. maxRate>0.
SpecificUnBoundedScoringPathsNumberN:=proc(possibleScores,n,startingLevel,maxRate) option remember:
	if n<0 or startingLevel+n*maxRate<0 then
		return(0):
	elif n=0 and startingLevel=0 then
		return(1):
	fi:
	add(SpecificUnBoundedScoringPathsNumberN(possibleScores,n-score[1],startingLevel+score[2],maxRate),score in possibleScores):
end:
#################################################
#A helper function that works if there are no (0,-) steps.
#Inputting maxRate just saves some computation time. It is determined by possibleScores. maxRate<0.
SpecificUnBoundedScoringPathsNumberP:=proc(possibleScores,n,startingLevel,maxRate) option remember:
	if n<0 or startingLevel+n*maxRate>0 then
		return(0):
	elif n=0 and startingLevel=0 then
		return(1):
	fi:
	add(SpecificUnBoundedScoringPathsNumberP(possibleScores,n-score[1],startingLevel+score[2],maxRate),score in possibleScores):
end:
######################################################################################################
#SpecificIrreducibleUnBoundedScoringPathsNumber(possibleScores,n,startingLevel): inputs possible scores (x-step,y-step) in the format {[x1,y1],[x2,y2],...}.
#Outputs the number of paths that start at (0,startingLevel), end at (n,0), and never touch the x-axis beforehand.
#There should NOT be BOTH (0,+) and (0,-) steps. Otherwise you get infinity as the answer.
#The answer is computed empirically.
#Try seq(SpecificIrreducibleUnBoundedScoringPathsNumber({[1,0],[1,-2],[2,-3],[0,1]},n1,0),n1=0..9).
SpecificIrreducibleUnBoundedScoringPathsNumber:=proc(possibleScores,n,startingLevel) local zeroPstep,zeroNstep,score,maxRate,count:
	zeroPstep:=false:
	zeroNstep:=false:
	for score in possibleScores do
		if score[1]=0 then
			if score[2]>0 then
				zeroPstep:=true:
			else	#score[2]<0
				zeroNstep:=true:
			fi:
		fi:
	od:
	if zeroPstep then
		if zeroNstep then
			return(infinity):
		else
			maxRate:=0:
			for score in possibleScores do
				if score[2]<0 then	#Then score[1]<>0 since zeroNstep=false
					maxRate:=min(maxRate,score[2]/score[1]):
				fi:
			od:
			if startingLevel=0 then
				count:=0:
				for score in possibleScores do
					if score[2]<>0 then	#Must take at least 1 step that changes altitude.
						count:=count+SpecificIrreducibleUnBoundedScoringPathsNumberP(possibleScores,n-score[1],score[2],maxRate):
					fi:
				od:
				return(count):
			else
				return(SpecificIrreducibleUnBoundedScoringPathsNumberP(possibleScores,n,startingLevel,maxRate)):
			fi:
		fi:
	fi:
	maxRate:=0:
	for score in possibleScores do
		if score[2]>0 then	#Then score[1]<>0 since zeroPstep=false
			maxRate:=max(maxRate,score[2]/score[1]):
		fi:
	od:
	if startingLevel=0 then
		count:=0:
		for score in possibleScores do
			if score[2]<>0 then	#Must take at least 1 step that changes altitude.
				count:=count+SpecificIrreducibleUnBoundedScoringPathsNumberN(possibleScores,n-score[1],score[2],maxRate):
			fi:
		od:
		count:
	else
		SpecificIrreducibleUnBoundedScoringPathsNumberN(possibleScores,n,startingLevel,maxRate):
	fi:
end:
#################################################
#A helper function that works if there are no (0,+) steps.
#Inputting maxRate just saves some computation time. It is determined by possibleScores.
SpecificIrreducibleUnBoundedScoringPathsNumberN:=proc(possibleScores,n,startingLevel,maxRate) option remember:
	if n<0 or startingLevel+n*maxRate<0 then
		return(0):
	elif startingLevel=0 then
		if n=0 then
			return(1):
		else
			return(0):
		fi:
	fi:
	add(SpecificIrreducibleUnBoundedScoringPathsNumberN(possibleScores,n-score[1],startingLevel+score[2],maxRate),score in possibleScores):
end:
#################################################
#A helper function that works if there are no (0,-) steps.
#Inputting maxRate just saves some computation time. It is determined by possibleScores.
SpecificIrreducibleUnBoundedScoringPathsNumberP:=proc(possibleScores,n,startingLevel,maxRate) option remember:
	if n<0 or startingLevel+n*maxRate>0 then
		return(0):
	elif startingLevel=0 then
		if n=0 then
			return(1):
		else
			return(0):
		fi:
	fi:
	add(SpecificIrreducibleUnBoundedScoringPathsNumberP(possibleScores,n-score[1],startingLevel+score[2],maxRate),score in possibleScores):
end:


###########################################################################################################################################
#Finding the proper root (minimal polynomial). This is only a helper function for parsing output from
#SemiBoundedScoringPaths, EqualSemiBoundedScoringPaths, and EqualUnBoundedScoringPaths.
#FindProperRoot(possibleScores,P,case,var,t,lowerLimit:=0):
#S is the step set. P is the polynomial originally found (not necessarily minimal).
#var is the g.f. we are interested in.
#t is the variable in the g.f.
#case=1: semi bounded free, 2: semi bounded equal, 3: semi bounded irreducible, 4: unbounded equal, 5: unbounded irreducible
FindProperRoot:=proc(possibleScores,P,case,var,t,lowerLimit:=0) local PF,currentPower,val,pf,LCM,l:
	PF:=factors(P)[2]:	#This is better than factor, as factor may be of type `+` or `^` or `*`.
	PF:={seq(pf[1],pf in PF)}:	#We only care about the unique factors, not multiplicity.
	currentPower:=0:
	val:=0:
	while nops(PF)>1 do
		if case=1 then		#var=k[y]
			val:=val+SpecificSemiBoundedScoringPathsNumber(possibleScores,lowerLimit,currentPower,op(1,var))*t^currentPower:
		elif case=2 then	#var=f[a,b]
			val:=val+SpecificEqualSemiBoundedScoringPathsNumber(possibleScores,lowerLimit-op(2,var),currentPower,op(1,var)-op(2,var))*t^currentPower:
		elif case=3 then	#var=g[a,b]
			val:=val+SpecificIrreducibleSemiBoundedScoringPathsNumber(possibleScores,lowerLimit-op(2,var),currentPower,op(1,var)-op(2,var))*t^currentPower:
		elif case=4 then	#var=GF
			val:=val+SpecificUnBoundedScoringPathsNumber(possibleScores,currentPower,0)*t^currentPower:
		elif case=5 then	#var=h[j]
			val:=val+SpecificIrreducibleUnBoundedScoringPathsNumber(possibleScores,currentPower,op(1,var))*t^currentPower:
		else
			RETURN(`ERROR! NOT A VALID CASE!`):
		fi:
		for pf in PF do
			if convert(taylor(subs(var=val,pf),t=0,currentPower+1),`polynom`)<>0 then
				PF:=PF minus {pf}:
			fi:
		od:
		currentPower:=currentPower+1:
	od:
#The following parses the polynomial to have only integer coefficients and combines the coefficients of var.
	PF:=PF[1]:	#This will always be at least a 2-term polynomial unless the step set yields trivial results.
	LCM:=1:
	for l in convert(PF,set) do	#Hence this will break into all of the summands.
		LCM:=lcm(LCM,denom(lcoeff(l))):
	od:
	add(factor(LCM*coeff(PF,var,i))*var^i,i=0..degree(PF,var)):#This is better than collect as this tries to factor all of the coefficients of K.
end:
###########################################################################################################################################


###########################################################################################################################################
#RatioOfWalks(S): Inputs a step set S.
#Outputs the asymptotic ratio of nonnegative excursions to bridges after finding recurrences for excursions and bridges respectively.
RatioOfWalks:=proc(S) local t,f,P1,commonValue,R1,E1,P2,R2,E2:
	P1:=EqualSemiBoundedScoringPaths(S,0,t,f):
	commonValue:=igcd(seq(degree(p,t), p in convert(expand(P1),set))):
	if degree(P1,f)<=8 then
		R1:=algtorec(subs(t=t^(1/commonValue),P1),f,t,n,N):
	else
		R1:=Findrec([seq(SpecificEqualSemiBoundedScoringPathsNumber(Steps,0,n1*commonValue,0),n1=1..(2+16)*(1+16)+1)],n,N,40):
		if R1=FAIL then
			print(`The degree of the minimal polynomial for excursions was large (>8) so we tried to guess at the recurrence and came up short.`):
			print(`algtorec should be able to convert the minimal polynomial into a recurrence but it would likely take a lot of time and memory.`):
			return(FAIL):
		fi:
	fi:
	E1:=AsyC(R1,n,N,0,[seq(SpecificEqualSemiBoundedScoringPathsNumber(S,0,n1*commonValue,0),n1=1..degree(R1,N))],1000):
	if E1=FAIL then
		print(`Error computing the asymptotic behavior of excursions.`):
		return(FAIL):
	fi:
	if degree(P1,f)<=8 then
		P2:=EqualUnBoundedScoringPaths(S,t,f):
		R2:=algtorec(subs(t=t^(1/commonValue),P2),f,t,n,N):
	else
		R2:=Findrec([seq(SpecificUnBoundedScoringPathsNumber(Steps,n1*commonValue,0),n1=1..(2+16)*(1+16)+1)],n,N,40):
		if R2=FAIL then
			print(`The degree of the minimal polynomial for excursions was large (>8) so we tried to guess at the recurrence for bridges and came up short.`):
			print(`algtorec should be able to convert the minimal polynomial into a recurrence but it would likely take a lot of time and memory.`):
			return(FAIL):
		fi:
	fi:
	E2:=AsyC(R2,n,N,0,[seq(SpecificUnBoundedScoringPathsNumber(S,n1*commonValue,0),n1=1..degree(R2,N))],1000):
	if E1=FAIL then
		print(`Error computing the asymptotic behavior of bridges.`):
		return(FAIL):
	fi:
	simplify(E1[1]*E1[2]/E2[1]/E2[2]):
end:
###########################################################################################################################################


###########################################################################################################################################
#AsymptoticBase(P,t,f): inputs a minimal polynomial in t and f: the generating function with formal variable t.
#Outputs the b such that the coefficients of f follow b^n asymptotically.
AsymptoticBase:=proc(P,t,f) local radius,l:
	radius:=infinity:
	for l in {evalf(solve(discrim(P,f),t))} do
		if Im(l)=0 and l>0 and l<radius then
			radius:=l:
		fi:
	od:
	1/radius:
end:
###########################################################################################################################################


###########################################################################################################################################
#Producing random step sets.
#RandomStepSet(numSteps,xBound,yBound): produces a set of steps each with x-value between 0 and xBound and y-value between -yBound and yBound.
#There cannot be a [0,0] step, thus we need yBound>=1. And there must be some step with x-value>0 so we should also have xBound>=1.
#Try RandomStepSet(8,4,3).
RandomStepSet:=proc(numSteps,xBound,yBound) local steps,rax,ray,xValue,yValue:
	steps:={}:
	rax:=rand(0..xBound):
	ray:=rand(-yBound..yBound):
	while nops(steps)<numSteps do
		xValue:=rax():
		yValue:=ray():
		if xValue=0 then
			while yValue=0 do
				yValue:=ray():
			od:
		fi:
		steps:=steps union {[xValue,yValue]}:
	od:
	steps:
end:
############################################################################
#Producing random step sets for the Bounded, EqualBounded, EqualSemiBounded and EqualUnBounded cases.
#The bounded cases could be expanded to include some pairs of (+) and (-) zero steps but I do not currently know the restrictions.
#RandomZeroStepSet(numSteps,xBound,yBound): produces a set of steps each with x-value between 0 and xBound and y-value between -yBound and yBound.
#There cannot be a [0,0] step, thus we need yBound>=1. And there must be some step with x-value>0 so we should also have xBound>=1.
#There also cannot be both (+) and (-) zero steps.
#Try RandomZeroStepSet(8,4,3).
RandomZeroStepSet:=proc(numSteps,xBound,yBound) local steps,rax,ray,xValue,yValue,SET:
	steps:={}:
	rax:=rand(0..xBound):
	ray:=rand(-yBound..yBound):
	SET:=0:
	while nops(steps)<numSteps do
		xValue:=rax():
		if xValue=0 then
			if SET=0 then
				yValue:=(-1)^rand(0..1)()*rand(1..yBound)():
				SET:=sign(yValue):
			elif SET=1 then
				yValue:=rand(1..yBound)():
			else	#SET=-1
				yValue:=rand(-yBound..-1)():
			fi:
		else
			yValue:=ray():
		fi:
		steps:=steps union {[xValue,yValue]}:
	od:
end:
############################################################################
#Producing random step sets for the SemiBounded case.
#RandomSemiBoundedStepSet(numSteps,xBound,yBound): produces a set of steps each with x-value between 0 and xBound and y-value between -yBound and yBound.
#There cannot be a [0,0] step, thus we need yBound>=1. And there must be some step with x-value>0 so we should also have xBound>=1.
#There also cannot be (+) zero steps.
#Try RandomSemiBoundedStepSet(8,4,3).
RandomSemiBoundedStepSet:=proc(numSteps,xBound,yBound) local steps,rax,ray,xValue,yValue:
	steps:={}:
	rax:=rand(0..xBound):
	ray:=rand(-yBound..yBound):
	while nops(steps)<numSteps do
		xValue:=rax():
		if xValue=0 then
			yValue:=rand(-yBound..-1)():
		else
			yValue:=ray():
		fi:
		steps:=steps union {[xValue,yValue]}:
	od:
	steps:
end:
############################################################################
#Producing random step sets for the UnBounded case. 
#RandomUnBoundedStepSet(numSteps,xBound,yBound): produces a set of steps each with x-value between 1 and xBound and y-value between -yBound and yBound.
#There cannot be zero steps.
#Try RandomUnBoundedStepSet(8,4,3).
RandomUnBoundedStepSet:=proc(numSteps,xBound,yBound) local steps,rax,ray:
	steps:={}:
	rax:=rand(1..xBound):
	ray:=rand(-yBound..yBound):
	while nops(steps)<numSteps do
		steps:=steps union {[rax(),ray()]}:
	od:
	steps:
end:
###########################################################################################################################################


###########################################################################################################################################
###########################################################################################################################################
#Producing Papers.
###########################################################################################################################################
###########################################################################################################################################
#PaperBounded1(Steps,upperLimit,lowerLimit,K): inputs a set of steps, Steps, and a positive integer K, outputs an article about the
#generating function of walks from the origin back to [n,0] for some n that never go below y=lowerLimit or above y=upperLimit.
#It also outputs the first K terms, for the sake of the OEIS.
#Try:
#PaperBounded1({[1,2],[1,-2],[1,3],[1,-3]},8,-2,20);
PaperBounded1:=proc(Steps,upperLimit,lowerLimit,K) local c,a,b,f,x,gu,mu,t0,i,n,positiveStep,zeroSteps,negativeStep,step:
	t0:=time():
	print(cat(`On the Generating Functions of Enumerating Sequences for Bounded Paths with set of steps `, Steps,`.`)):
	print(``):
	print(`By Bryan Ek's Maple Package ScoringPaths`):
	print(``):
	print(`Theorem:`):
	print(cat(`Let A(n) be the number of ways of walking from (0,0) to (n,0) using the steps `, Steps, ` and never going below y=lowerLimit or above y=upperLimit.`)):
	print(``):

	positiveStep:=false:
#	zeroStep:=false:
	zeroSteps:={}:
	negativeStep:=false:
	for step in (convert(Steps,set) minus {[0,0]}) do	#while not (positiveStep and negativeStep and zeroStep) do	#This doesn’t actually save time since we would need all of the zero steps anyway.
		if step[2]>0 and step[2]<=upperLimit-lowerLimit then	#Ensures this step could feasibly be used.
			positiveStep:=true:
		elif step[2]<0 and step[2]>=lowerLimit-upperLimit then	#Ensures this step could feasibly be used.
			negativeStep:=true:
		else
			zeroSteps:=zeroSteps union {step[1]}:
		fi:
	od:
	if not (negativeStep and positiveStep) then
		if nops(zeroSteps)=0 then
			print(`The problem is trivial. There is no walk back to the x-axis.`):
			print(`----------------------------------------------------------------------------------------------------------------`):
			RETURN(FAIL):
		elif nops(zeroSteps)=1 then
			print(cat(`The problem is trivial. Walks returning to the x-axis consist of only repeated use of the one step `, zeroSteps[1])):
			print(`----------------------------------------------------------------------------------------------------------------`):
			RETURN(FAIL):
		else
			print(`The only relevant steps are:`):
			print(zeroSteps):
			gu:=EqualBoundedScoringPaths(zeroSteps,upperLimit,lowerLimit,x,f):
			print(`The problem is equivalent to summing to n in different combinations using parts of size:`):
			print(seq(st,st in zeroSteps)):
			print(``):
		fi:
	else
		gu:=EqualBoundedScoringPaths(Steps,upperLimit,lowerLimit,x,f):
	fi:

	print(`Let f(x) be its (ordinary) generating function:`):
	print(f(x)=Sum(A(n)*x^n,n=0..infinity)):
	print(``):
	print(`Then f(x) is given by`):
	print(``):
	print(gu):
	print(``):
	print(`and in Maple input format`):
	print(``):
	lprint(gu):
	print(``):
	print(cat(`For the sake of the OEIS, here are the first `, K+1,` terms (starting at n=0).`)):
	mu:=taylor(gu,x=0,K+1):
	lprint([seq( coeff(mu,x,i), i=0..K)]):
#The following might be faster.
#	lprint([seq(SpecificEqualBoundedScoringPathsNumber(Steps,upperLimit,lowerLimit,n1,0),n1=0..K)]):
	print(``):
	print(`----------------------------`):
	print(cat(`This ends this paper that took `, time()-t0, ` seconds to generate.`)):
	print(`----------------------------------------------------------------------------------------------------------------`):
end:

###########################################################################################################################################
#PaperBounded2(Steps,K): inputs a set of steps, Steps, and a positive integer K, outputs an article about the
#generating function of walks from the origin back to [n,0] for some n that never go further than one step from the x-axis.
#It also outputs the first K terms, for the sake of the OEIS.
#Try:
#PaperBounded2({[1,2],[1,-2],[1,3],[1,-3]},20);
PaperBounded2:=proc(Steps,K) local c,a,b,f,x,gu,mu,t0,i,n,lowerLimit,upperLimit,zeroSteps,step,upstep,downstep:
	t0:=time():
	print(cat(`On the Generating Functions of Enumerating Sequences for Bounded Paths with set of steps `, Steps,`.`)):
	print(``):
	print(`By Bryan Ek's Maple Package ScoringPaths`):
	print(``):
	print(`Theorem:`):
	print(cat(`Let A(n) be the number of ways of walking from (0,0) to (n,0) using the steps `, Steps, ` and never going further than one step from the x-axis.`)):
	print(``):

	lowerLimit:=1:
	zeroSteps:={}:
	upperLimit:=-1:
	upstep:=false:
	downstep:=false:
	for step in (convert(Steps,set) minus {[0,0]}) while not (upstep and downstep) do	#while not (positiveStep and negativeStep and zeroStep) do	#This doesn’t actually save time since we would need all of the zero steps anyway.
		if step[2]>0 then
			lowerLimit:=min(lowerLimit,-step[2]):
			if step[1]=0 then
				upstep:=true:
			fi:
		elif step[2]<0 then
			upperLimit:=max(upperLimit,-step[2]):
			if step[1]=0 then
				downstep:=true:
			fi:
		else
			zeroSteps:=zeroSteps union {step}:
		fi:
	od:
	if upstep and downstep then
		print(`There are an infinite number of ways to have walks of length n since you included step(s) directly up and down.`):
		print(`----------------------------------------------------------------------------------------------------------------`):
		RETURN(FAIL):
	elif lowerLimit>0 or upperLimit<0 then
		if nops(zeroSteps)=0 then
			print(`The problem is trivial. There is no walk back to the x-axis.`):
			print(`----------------------------------------------------------------------------------------------------------------`):
			RETURN(FAIL):
		elif nops(zeroSteps)=1 then
			print(cat(`The problem is trivial. Walks returning to the x-axis consist of only repeated use of the one step `, zeroSteps[1])):
			print(`----------------------------------------------------------------------------------------------------------------`):
			RETURN(FAIL):
		else
			print(`The only relevant steps are:`):
			print(zeroSteps):
			gu:=EqualBoundedScoringPaths(zeroSteps,upperLimit,lowerLimit,x,f):
			print(`The problem is equivalent to summing to n in different combinations using parts of size:`):
			print(seq(st[1],st in zeroSteps)):
			print(``):
		fi:
	else
		gu:=EqualBoundedScoringPaths(Steps,upperLimit,lowerLimit,x,f):
	fi:

	print(`Let f(x) be its (ordinary) generating function:`):
	print(f(x)=Sum(A(n)*x^n,n=0..infinity)):
	print(``):
	print(`Then f(x) is given by`):
	print(``):
	print(gu):
	print(``):
	print(`and in Maple input format`):
	print(``):
	lprint(gu):
	print(``):
	print(cat(`For the sake of the OEIS, here are the first `, K+1,` terms (starting at n=0).`)):
	mu:=taylor(gu,x=0,K+1):
	lprint([seq( coeff(mu,x,i), i=0..K)]):
#The following might be faster.
#	lprint([seq(SpecificEqualBoundedScoringPathsNumber(Steps,upperLimit,lowerLimit,n1,0),n1=0..K)]):
	print(``):
	print(`----------------------------`):
	print(cat(`This ends this paper that took `, time()-t0, ` seconds to generate.`)):
	print(`----------------------------------------------------------------------------------------------------------------`):
end:

###########################################################################################################################################
#PaperSemiBounded(Steps,K,recurrence:=false,asymptotics:=false): inputs a set of steps, Steps, and a positive integer K, outputs an article about the
#generating function of walks from the origin back to [n,0] for some n that never go below the x-axis.
#It also outputs the first K terms, for the sake of the OEIS.
#If you would a recurrence or asymptotic behavior for the sequence, set the appropriate boolean to true.
#Try:
#PaperSemiBounded({[1,1],[1,-1],[1,0]},40,true);
PaperSemiBounded:=proc(Steps,K,recurrence:=false,asymptotics:=false) local c,a,b,f,x,gu,mu,t0,i,n,R,EXP,step,commonValue,positiveStep,zeroSteps:
	t0:=time():
	print(cat(`On the Generating Functions of Enumerating Sequences for Generalized Paths with set of steps `, Steps,`.`)):
	print(``):
	print(`By Bryan Ek's Maple Package ScoringPaths`):
	print(``):
	print(`Theorem:`):
	print(cat(`Let A(n) be the number of ways of walking from (0,0) using the steps `, Steps, ` and never going below the x-axis.`)):
	print(``):

	positiveStep:=false:
	zeroSteps:={}:
	for step in (convert(Steps,set) minus {[0,0]}) do
		if step[2]>0 then
			if step[1]=0 then
				print(`There are an infinite number of ways to have meanders of length n since you included step(s) directly up.`):
				print(`----------------------------------------------------------------------------------------------------------------`):
				RETURN(FAIL):
			else
				positiveStep:=true:
			fi:
		elif step[2]=0 then
			zeroSteps:=zeroSteps union {step}:
		fi:
	od:
	if not positiveStep then
		if nops(zeroSteps)=0 then
			print(`The problem is trivial. There is no non-stationary walk.`):
			print(`----------------------------------------------------------------------------------------------------------------`):
			RETURN(FAIL):
		elif nops(zeroSteps)=1 then
			print(cat(`The problem is trivial. Walks consist of only repeated use of the one step `, zeroSteps[1])):
			print(`----------------------------------------------------------------------------------------------------------------`):
			RETURN(FAIL):
		else
			print(`The only relevant steps are:`):
			print(zeroSteps):
			gu:=SemiBoundedScoringPaths(zeroSteps,0,x,f):
			print(`The problem is equivalent to summing to n in different combinations using parts of size:`):
			print(seq(st[1],st in zeroSteps)):
			print(``):
		fi:
	else
		gu:=SemiBoundedScoringPaths(Steps,0,x,f):
	fi:

	print(`Let f(x) be its (ordinary) generating function:`):
	print(f(x)=Sum(A(n)*x^n,n=0..infinity)):
	print(``):
	print(`Then f(x) satisfies`):
	print(``):
	print(gu=0):
	print(``):
	print(`and in Maple input format`):
	print(``):
	lprint(gu=0):
	print(``):
	print(cat(`For the sake of the OEIS, here are the first `, K+1,` terms (starting at n=0).`)):
#	mu:=taylor(RootOf(subs(f=_Z,gu)),x=0,K+1):
#	lprint([seq( coeff(mu,x,i), i=0..K)]):
#The following is faster for small-ish values of K.
	commonValue:=igcd(seq(degree(p,x),p in convert(expand(gu),set))):
	if commonValue>1 then
		print(cat(`The sequence is only non-zero every `,commonValue,` terms so the sequence has been trimmed for ease of reading.`)):
	fi:
	lprint([seq(SpecificSemiBoundedScoringPathsNumber(Steps,0,n1*commonValue,0),n1=0..K)]):
	if recurrence then
		print(``):
		print(`The coefficients satisfy the following recurrence:`):
		if degree(gu,f)<8 then
			R:=algtorec(subs(x=x^(1/commonValue),gu),f,x,n,N):	#This is likely not the lowest order recurrence but it provides an upper bound of search space.
#			R:=Findrec([seq(SpecificSemiBoundedScoringPathsNumber(Steps,0,n1*commonValue,0),n1=1..(2+degree(R,n))*(1+degree(R,N))+1)],n,N,degree(R,N)+degree(R,n)):	#This is likely a lower order recurrence.
		else
			R:=Findrec([seq(SpecificSemiBoundedScoringPathsNumber(Steps,0,n1*commonValue,0),n1=1..(2+16)*(1+16)+1)],n,N,40):
		fi:
		if R<>FAIL then
			print(R):
			print(`and in Maple input format`):
			lprint(R):
			if asymptotics then
				print(``):
				print(`The coefficients are asymptotic to`):
				if degree(R,N)<15 then	#Arbitrary cutoff to make it not run "forever".
					EXP:=AsyC(R,n,N,0,[seq(SpecificSemiBoundedScoringPathsNumber(Steps,0,n1*commonValue,0),n1=1..degree(R,N))],10000):
					if EXP<>FAIL then
						print(factor(simplify(evalf(EXP[1]*EXP[2])))):
					else
						print(`There was an error computing asymptotics from the recurrence.`):
						print(`The following is purely the exponential growth term.`):
						print(AsymptoticBase(gu,x,f)^n):
					fi:
				else
					print(`The following is purely the exponential growth term.`):
					print(AsymptoticBase(gu,x,f)^n):
				fi:
			fi:
		else
			print(`We were unable to find a recurrence.`):
			print(`This is because the conversion would have taken too long; the degree of the minimal polynomial was too high. (Arbitrary cutoff of 8).`):
			print(`We then tried to guess at a recurrence but the search space was not large enough. (Arbitrary cutoff of 40 for degree+order).`):
			if asymptotics then
				print(``):
				print(`The following is purely the exponential growth term.`):
				print(`The coefficients are asymptotic to`):
				print(AsymptoticBase(gu,x,f)^n):
			fi:
		fi:
	elif asymptotics then
		print(`The coefficients are asymptotic to`):
		print(`The following is purely the exponential growth term.`):
		print(AsymptoticBase(gu,x,f)^n):
	fi:
	print(``):
	print(`----------------------------`):
	print(cat(`This ends this paper that took `, time()-t0, ` seconds to generate.`)):
	print(`----------------------------------------------------------------------------------------------------------------`):
end:

###########################################################################################################################################
#PaperEqualSemiBounded(Steps,K,recurrence:=false,asymptotics:=false): inputs a set of steps, Steps, and a positive integer K, outputs an article about the
#generating function of walks from the origin back to [n,0] for some n that never go below the x-axis.
#It also outputs the first K terms, for the sake of the OEIS.
#If you would a recurrence or asymptotic behavior for the sequence, set the appropriate boolean to true.
#Try:
#PaperEqualSemiBounded({[1,1],[1,-1],[1,0]},40,true);
PaperEqualSemiBounded:=proc(Steps,K,recurrence:=false,asymptotics:=false) local c,a,b,f,x,gu,mu,t0,i,n,R,EXP,positiveStep,zeroSteps,negativeStep,step,commonValue,upstep,downstep:
	t0:=time():
	print(cat(`On the Generating Functions of Enumerating Sequences for Generalized Dyck Paths with set of steps `, Steps,`.`)):
	print(``):
	print(`By Bryan Ek's Maple Package ScoringPaths`):
	print(``):
	print(`Theorem:`):
	print(cat(`Let A(n) be the number of ways of walking from (0,0) to (n,0) using the steps `, Steps, ` and never going below the x-axis.`)):
	print(``):

	positiveStep:=false:
#	zeroStep:=false:
	zeroSteps:={}:
	negativeStep:=false:
	upstep:=false:
	downstep:=false:
	for step in (convert(Steps,set) minus {[0,0]}) while not (upstep and downstep) do	#while not (positiveStep and negativeStep and zeroStep) do	#This doesn’t actually save time since we would need all of the zero steps anyway.
		if step[2]>0 then
			positiveStep:=true:
			if step[1]=0 then
				upstep:=true:
			fi:
		elif step[2]<0 then
			negativeStep:=true:
			if step[1]=0 then
				downstep:=true:
			fi:
		else
			zeroSteps:=zeroSteps union {step}:
		fi:
	od:
	if upstep and downstep then
		print(`There are an infinite number of ways to have excursions of length n since you included step(s) directly up and down.`):
		print(`----------------------------------------------------------------------------------------------------------------`):
		RETURN(FAIL):
	elif not (negativeStep and positiveStep) then
		if nops(zeroSteps)=0 then
			print(`The problem is trivial. There is no walk back to the x-axis.`):
			print(`----------------------------------------------------------------------------------------------------------------`):
			RETURN(FAIL):
		elif nops(zeroSteps)=1 then
			print(cat(`The problem is trivial. Walks returning to the x-axis consist of only repeated use of the one step `, zeroSteps[1])):
			print(`----------------------------------------------------------------------------------------------------------------`):
			RETURN(FAIL):
		else
			print(`The only relevant steps are:`):
			print(zeroSteps):
			gu:=EqualSemiBoundedScoringPaths(zeroSteps,0,x,f):
			print(`The problem is equivalent to summing to n in different combinations using parts of size:`):
			print(seq(st[1],st in zeroSteps)):
			print(``):
		fi:
	else
		gu:=EqualSemiBoundedScoringPaths(Steps,0,x,f):
	fi:

	print(`Let f(x) be its (ordinary) generating function:`):
	print(f(x)=Sum(A(n)*x^n,n=0..infinity)):
	print(``):
	print(`Then f(x) satisfies`):
	print(``):
	print(gu=0):
	print(``):
	print(`and in Maple input format`):
	print(``):
	lprint(gu=0):
	print(``):
	print(cat(`For the sake of the OEIS, here are the first `, K+1,` terms (starting at n=0).`)):
#	mu:=taylor(RootOf(subs(f=_Z,gu)),x=0,K+1):
#	lprint([seq( coeff(mu,x,i), i=0..K)]):
#The following is faster for small-ish values of K.
	commonValue:=igcd(seq(degree(p,x),p in convert(expand(gu),set))):
	if commonValue>1 then
		print(cat(`The sequence is only non-zero every `,commonValue,` terms so has been trimmed for ease of reading.`)):
	fi:
	lprint([seq(SpecificEqualSemiBoundedScoringPathsNumber(Steps,0,n1*commonValue,0),n1=0..K)]):
	if recurrence then
		print(``):
		print(`The coefficients satisfy the following recurrence:`):
		if degree(gu,f)<8 then
			R:=algtorec(subs(x=x^(1/commonValue),gu),f,x,n,N):	#This is likely not the lowest order recurrence but it provides an upper bound of search space.
#			R:=Findrec([seq(SpecificEqualSemiBoundedScoringPathsNumber(Steps,0,n1*commonValue,0),n1=1..(2+degree(R,n))*(1+degree(R,N))+1)],n,N,degree(R,N)+degree(R,n)):	#This is likely a lower order recurrence.
		else
			R:=Findrec([seq(SpecificEqualSemiBoundedScoringPathsNumber(Steps,0,n1*commonValue,0),n1=1..(2+16)*(1+16)+1)],n,N,40):
		fi:
		if R<>FAIL then
			print(R):
			print(`and in Maple input format`):
			lprint(R):
			if asymptotics then
				print(``):
				print(`The coefficients are asymptotic to`):
				if degree(R,N)<15 then	#Arbitrary cutoff to make it not run "forever".
					EXP:=AsyC(R,n,N,0,[seq(SpecificEqualSemiBoundedScoringPathsNumber(Steps,0,n1*commonValue,0),n1=1..degree(R,N))],10000):
					if EXP<>FAIL then
						print(factor(simplify(evalf(EXP[1]*EXP[2])))):
					else
						print(`There was an error computing asymptotics from the recurrence.`):
						print(`The following is purely the exponential growth term.`):
						print(AsymptoticBase(gu,x,f)^n):
					fi:
				else
					print(`The following is purely the exponential growth term.`):
					print(AsymptoticBase(gu,x,f)^n):
				fi:
			fi:
		else
			print(`We were unable to find a recurrence.`):
			print(`This is because the conversion would have taken too long; the degree of the minimal polynomial was too high. (Arbitrary cutoff of 8).`):
			print(`We then tried to guess at a recurrence but the search space was not large enough. (Arbitrary cutoff of 40 for degree+order).`):
			if asymptotics then
				print(``):
				print(`The following is purely the exponential growth term.`):
				print(`The coefficients are asymptotic to`):
				print(AsymptoticBase(gu,x,f)^n):
			fi:
		fi:
	elif asymptotics then
		print(`The coefficients are asymptotic to`):
		print(`The following is purely the exponential growth term.`):
		print(AsymptoticBase(gu,x,f)^n):
	fi:
	print(``):
	print(`----------------------------`):
	print(cat(`This ends this paper that took `, time()-t0, ` seconds to generate.`)):
	print(`----------------------------------------------------------------------------------------------------------------`):
end:
###############
with(combinat):
###########################################################################################################################################
#BookEqualSemiBounded(Steps,K,recurrence:=false,asymptotics:=false): inputs a set of steps, Steps, and a positive integer K.
#Outputs a book about the generating functions of walks from the origin back to (n,0) for some n that never go below the x-axis.
#Iterates over all subsets of Steps.
#If you would a recurrence or asymptotic behavior for the sequence, set the appropriate boolean to true.
#Try:
#BookEqualSemiBounded({[1,1],[1,-1],[1,0],[2,2],[2,-2],[3,3],[3,-1]},40,true);
BookEqualSemiBounded:=proc(Steps,K,recurrence:=false,asymptotics:=false) local i,t0,nonTrivialCount,result,steps,sectionCount:
	t0:=time():
	print(`Generating Functions of Excursions Drawn From`):
	print(Steps):
	print(``):
	print(`By Bryan Ek’s Maple Package ScoringPaths`):
	print(``):
	print(`------------------------------------------------------------------------------------------------------------------------------------------------------------------`):
	print(``):
	print(`Chapter 1: Introduction`):
	print(`This book includes results for all excursions built from subsets of size >=2 of the step set:`):
	print(Steps):
	print(``):
	if recurrence then
		print(`Results include finding the minimal polynomial for the generating function (g.f.),`):
		print(`a linear recurrence with polynomial coefficients for the sequence of walks, `):
		if asymptotics then
			print(`the asymptotic behavior of the sequence (only a first order approximation), `):
			print(cat(`and the first `,K+1,` nontrivial terms of each sequence (starting at n=0).`)):
			print(`Nontrivial means that if the sequence is nonzero only every k terms, then we output term n=0,k,2*k,...`):
			print(`The minimal polynomial has NOT been adjusted to accommodate the skipping, but the recurrence and asymptotic has.`):
			print(`For a better asymptotic approximation, you can use AsyC with a higher m (approximation order).`):
		else
			print(cat(`and the first `,K+1,` nontrivial terms of each sequence (starting at n=0).`)):
			print(`Nontrivial means that if the sequence is nonzero only every k terms, then we output term n=0,k,2*k,...`):
			print(`The minimal polynomial has NOT been adjusted to accommodate the skipping, but the recurrence has.`):
		fi:
	else
		print(`Results include finding the minimal polynomial for the generating function (g.f.)`):
		print(cat(`and the first `,K+1,` nontrivial terms of each sequence (starting at n=0).`)):
		print(`Nontrivial means that if the sequence is nonzero only every k terms, then we output term n=0,k,2*k,...`):
		print(`The minimal polynomial has NOT been adjusted to accommodate the skipping.`):
	fi:
	print(`-------------------------------------------------------------------------------------------------------------------------------`):
	print(``):
	nonTrivialCount:=0:
	for i from 2 to nops(Steps) do
		print(cat(`Chapter `,i,`: Walks built from `,i,` steps.`)):
		print(`----------------------------------------------------------------------------------------------------------------`):
		sectionCount:=0:
		for steps in choose(Steps,i) do
			sectionCount:=sectionCount+1:
			print(cat(`Section `,i,`.`,sectionCount)):
			result:=PaperEqualSemiBounded(steps,K,recurrence,asymptotics):
			if result<>FAIL then
				nonTrivialCount:=nonTrivialCount+1:
			fi:
			print(``):
		od:
		print(`-------------------------------------------------------------------------------------------------------------------------------`):
		print(``):
	od:
	print(cat(`This ends this book with `,nops(Steps),` Chapters and `,2^(nops(Steps))-nops(Steps)-1,` sections, `,nonTrivialCount,` of which are nontrivial, that took `,time()-t0, ` seconds to generate.`)):
end:

###########################################################################################################################################
#PaperEqualSemiBoundedPair(a,b,K,recurrence:=false,asymptotics:=false): inputs positive integers a and b that are relatively prime, and a positive integer K, outputs an article about the
#generating function of walks from the origin back to [(a+b)*n,0] for some n that never go below the x-axis.
#It also outputs the first K terms, for the sake of the OEIS.
#If you would a recurrence or asymptotic behavior for the sequence, set the appropriate boolean to true.
#Try:
#PaperEqualSemiBoundedPair(2,3,40,true);
PaperEqualSemiBoundedPair:=proc(a,b,K,recurrence:=false,asymptotics:=false) local c,f,x,gu,mu,t0,i,n,R,EXP:
	if gcd(a,b)<>1 then
		print(`a and b should be coprime.`):
		print(`----------------------------------------------------------------------------------------------------------------`):
		RETURN():
	elif a<=0 or b<=0 then
		print(`a and b should both be positive.`):
		print(`----------------------------------------------------------------------------------------------------------------`):
		RETURN():
	fi:
	t0:=time():
	print(cat(`On the Generating Functions of Enumerating Sequences for Generalized Dyck Paths with steps [1,`,-b,`] and [1,`,a,`].`)):
	print(``):
	print(`By Bryan Ek's Maple Package ScoringPaths`):
	print(``):
	print(`Theorem:`):
	print(cat(`Let A(n) be the number of ways of walking from (0,0)  to `, ((a+b)*n,0),  `using the steps `, {[1,a],[1,-b]}, ` and never going below the x-axis.`)):
	print(``):
	print(`Let f(x) be its (ordinary) generating function`):
	print(f(x)=Sum(A(n)*x^n,n=0..infinity)):
	print(``):
	gu:=EqualSemiBoundedScoringPaths({[1,a],[1,-b]},0,x,f):
	gu:=subs(x=x^(1/(a+b)),gu):
	print(`Then f(x) satisfies`):
	print(``):
	print(gu=0):
	print(``):
	print(`and in Maple input format`):
	print(``):
	lprint(gu=0):
	print(``):
	print(cat(`For the sake of the OEIS, here are the first `, K+1,` terms (starting at n=0).`)):
#	mu:=taylor(RootOf(subs(f=_Z,gu)),x=0,K+1):
#	lprint([seq( coeff(mu,x,i), i=0..K)]):
#The following is faster for small-ish values of K.
	lprint([seq(SpecificEqualSemiBoundedScoringPathsNumber({[1,a],[1,-b]},0,n1*(a+b),0),n1=0..K)]):
	if recurrence then
		print(``):
		print(`The coefficients satisfy the following recurrence:`):
		if degree(gu,f)<8 then
			R:=algtorec(gu,f,x,n,N):	#This is likely not the lowest order recurrence but it provides an upper bound of search space.
#			R:=Findrec([seq(SpecificEqualSemiBoundedScoringPathsNumber({[1,a],[1,-b]},0,n1*(a+b),0),n1=1..(2+degree(R,n))*(1+degree(R,N))+1)],n,N,degree(R,N)+degree(R,n)):	#This is likely a lower order recurrence.
		else
			R:=Findrec([seq(SpecificEqualSemiBoundedScoringPathsNumber({[1,a],[1,-b]},0,n1*(a+b),0),n1=1..(2+20)*(1+20)+1)],n,N,40):
		fi:
		if R<>FAIL then
			print(R):
			print(`and in Maple input format`):
			lprint(R):
			if asymptotics then
				print(``):
				print(`The coefficients are asymptotic to`):
				if degree(R,N)<15 then	#Arbitrary cutoff to make it not run "forever".
					EXP:=AsyC(R,n,N,0,[seq(SpecificEqualSemiBoundedScoringPathsNumber({[1,a],[1,-b]},0,n1*(a+b),0),n1=1..degree(R,N))],10000):
					if EXP<>FAIL then
						print(factor(simplify(evalf(EXP[1]*EXP[2])))):
					else
						print(`There was an error computing asymptotics from the recurrence.`):
						print(`The following is purely the exponential growth term.`):
						print(AsymptoticBase(gu,x,f)^n):
					fi:
				else
					print(`The following is purely the exponential growth term.`):
					print(AsymptoticBase(gu,x,f)^n):
				fi:
			fi:
		else
			print(`We were unable to find a recurrence.`):
			print(`This is because the conversion would have taken too long; the degree of the minimal polynomial was too high. (Arbitrary cutoff of 8).`):
			print(`We then tried to guess at a recurrence but the search space was not large enough. (Arbitrary cutoff of 40 for degree+order).`):
			if asymptotics then
				print(``):
				print(`The following is purely the exponential growth term.`):
				print(`The coefficients are asymptotic to`):
				print(AsymptoticBase(gu,x,f)^n):
			fi:
		fi:
	elif asymptotics then
		print(`The coefficients are asymptotic to`):
		print(`The following is purely the exponential growth term.`):
		print(AsymptoticBase(gu,x,f)^n):
	fi:
	print(``):
	print(`----------------------------`):
	print(cat(`This ends this paper that took `, time()-t0, ` seconds to generate.`)):
	print(`----------------------------------------------------------------------------------------------------------------`):
end:

###########################################################################################################################################
#PaperUnBounded(Steps,K,recurrence:=false,asymptotics:=false): inputs a set of steps, Steps, and a positive integer K, outputs an article about the
#generating function of walks from the origin back to [n,0] for some n.
#It also outputs the first K terms, for the sake of the OEIS.
#If you would a recurrence or asymptotic behavior for the sequence, set the appropriate boolean to true.
#Try:
#PaperUnBounded({[1,1],[1,-1],[1,0]},40,true);
PaperUnBounded:=proc(Steps,K,recurrence:=false,asymptotics:=false) local c,a,b,f,x,gu,mu,t0,i,n,R,EXP,positiveStep,zeroSteps,negativeStep,step,commonValue,upstep,downstep:
	t0:=time():
	print(cat(`On the Generating Functions of Enumerating Sequences for Unbounded Paths with set of steps `, Steps,`.`)):
	print(``):
	print(`By Bryan Ek's Maple Package ScoringPaths`):
	print(``):
	print(`Theorem:`):
	print(cat(`Let A(n) be the number of ways of walking from (0,0) to (n,0) using the steps `, Steps)):
	print(``):

	positiveStep:=false:
#	zeroStep:=false:
	zeroSteps:={}:
	negativeStep:=false:
	upstep:=false:
	downstep:=false:
	for step in (convert(Steps,set) minus {[0,0]}) while not (upstep and downstep) do	#while not (positiveStep and negativeStep and zeroStep) do	#This doesn’t actually save time since we would need all of the zero steps anyway.
		if step[2]>0 then
			positiveStep:=true:
			if step[1]=0 then
				upstep:=true:
			fi:
		elif step[2]<0 then
			negativeStep:=true:
			if step[1]=0 then
				downstep:=true:
			fi:
		else
			zeroSteps:=zeroSteps union {step}:
		fi:
	od:
	if upstep and downstep then
		print(`There are an infinite number of ways to have bridges of length n since you included step(s) directly up and down.`):
		print(`----------------------------------------------------------------------------------------------------------------`):
		RETURN():
	elif not (negativeStep and positiveStep) then
		if nops(zeroSteps)=0 then
			print(`The problem is trivial. There is no walk back to the x-axis.`):
			print(`----------------------------------------------------------------------------------------------------------------`):
			RETURN():
		elif nops(zeroSteps)=1 then
			print(cat(`The problem is trivial. Walks returning to the x-axis consist of only repeated use of the one step `, zeroSteps[1])):
			print(`----------------------------------------------------------------------------------------------------------------`):
			RETURN(FAIL):
		else
			print(`The only relevant steps are:`):
			print(zeroSteps):
			gu:=EqualUnBoundedScoringPaths(zeroSteps,x,f):
			print(`The problem is equivalent to summing to n in different combinations using parts of size:`):
			print(seq(st[1],st in zeroSteps)):
			print(``):
		fi:
	else
		gu:=EqualUnBoundedScoringPaths(Steps,x,f):
	fi:

	print(`Let f(x) be its (ordinary) generating function:`):
	print(f(x)=Sum(A(n)*x^n,n=0..infinity)):
	print(``):
	print(`Then f(x) satisfies`):
	print(``):
	print(gu=0):
	print(``):
	print(`and in Maple input format`):
	print(``):
	lprint(gu=0):
	print(``):
	print(cat(`For the sake of the OEIS, here are the first `, K+1,` terms (starting at n=0).`)):
#	mu:=taylor(RootOf(subs(f=_Z,gu)),x=0,K+1):
#	lprint([seq( coeff(mu,x,i), i=0..K)]):
#The following is faster for small-ish values of K. And the above is typically not accurate.
	commonValue:=igcd(seq(degree(p,x),p in convert(expand(gu),set))):
	if commonValue>1 then
		print(cat(`The sequence is only non-zero every `,commonValue,` terms so has been trimmed for ease of reading.`)):
	fi:
	lprint([seq(SpecificUnBoundedScoringPathsNumber(Steps,n1*commonValue,0),n1=0..K)]):
	if recurrence then
		print(``):
		print(`The coefficients satisfy the following recurrence:`):
		if degree(gu,f)<8 then
			R:=algtorec(subs(x=x^(1/commonValue),gu),f,x,n,N):	#This is likely not the lowest order recurrence but it provides an upper bound of search space.
#			R:=Findrec([seq(SpecificUnBoundedScoringPathsNumber(Steps,n1*commonValue,0),n1=1..(2+degree(R,n))*(1+degree(R,N))+1)],n,N,degree(R,N)+degree(R,n)):	#This is likely a lower order recurrence.
		else
			R:=Findrec([seq(SpecificUnBoundedScoringPathsNumber(Steps,n1*commonValue,0),n1=1..(2+16)*(1+16)+1)],n,N,40):
		fi:
		if R<>FAIL then
			print(R):
			print(`and in Maple input format`):
			lprint(R):
			if asymptotics then
				print(``):
				print(`The coefficients are asymptotic to`):
				if degree(R,N)<15 then	#Arbitrary cutoff to make it not run "forever".
					EXP:=AsyC(R,n,N,0,[seq(SpecificUnBoundedScoringPathsNumber(Steps,n1*commonValue,0),n1=1..degree(R,N))],10000):
					if EXP<>FAIL then
						print(factor(simplify(evalf(EXP[1]*EXP[2])))):
					else
						print(`There was an error computing asymptotics from the recurrence.`):
						print(`The following is purely the exponential growth term.`):
						print(AsymptoticBase(gu,x,f)^n):
					fi:
				else
					print(`The following is purely the exponential growth term.`):
					print(AsymptoticBase(gu,x,f)^n):
				fi:
			fi:
		else
			print(`We were unable to find a recurrence.`):
			print(`This is because the conversion would have taken too long; the degree of the minimal polynomial was too high. (Arbitrary cutoff of 8).`):
			print(`We then tried to guess at a recurrence but the search space was not large enough. (Arbitrary cutoff of 40 for degree+order).`):
			if asymptotics then
				print(``):
				print(`The following is purely the exponential growth term.`):
				print(`The coefficients are asymptotic to`):
				print(AsymptoticBase(gu,x,f)^n):
			fi:
		fi:
	elif asymptotics then
		print(`The coefficients are asymptotic to`):
		print(`The following is purely the exponential growth term.`):
		print(AsymptoticBase(gu,x,f)^n):
	fi:
	print(``):
	print(`----------------------------`):
	print(cat(`This ends this paper that took `, time()-t0, ` seconds to generate.`)):
	print(`----------------------------------------------------------------------------------------------------------------`):
end:

###########################################################################################################################################
#PaperUnBoundedPair(a,b,K,recurrence:=false,asymptotics:=false): inputs positive integers a and b that are relatively prime, and a positive integer K, outputs an article about the
#generating function of walks from the origin back to [(a+b)*n,0] for some n.
#It also outputs the first K terms, for the sake of the OEIS.
#If you would a recurrence or asymptotic behavior for the sequence, set the appropriate boolean to true.
#Try:
#PaperUnBoundedPair(2,3,40,true);
PaperUnBoundedPair:=proc(a,b,K,recurrence:=false,asymptotics:=false) local c,f,x,gu,mu,t0,i,n,R,EXP:
	if gcd(a,b)<>1 then
		print(`a and b should be coprime.`):
		print(`----------------------------------------------------------------------------------------------------------------`):
		RETURN():
	elif a<=0 or b<=0 then
		print(`a and b should both be positive.`):
		print(`----------------------------------------------------------------------------------------------------------------`):
		RETURN():
	fi:
	t0:=time():
	print(cat(`On the Generating Functions of Enumerating Sequences for Unbounded Paths with steps [1,`,-b,`] and [1,`,a,`].`)):
	print(``):
	print(`By Bryan Ek's Maple Package ScoringPaths`):
	print(``):
	print(`Theorem:`):
	print(cat(`Let A(n) be the number of ways of walking from (0,0)  to ` , ((a+b)*n,0), ` using the steps `, {[1,a],[1,-b]})):
	print(cat(`Then we have the simple solution: A(n)=`,(a+b),`*n choose `,a,`*n.`)):
	print(``):
	print(`We also include the following for the completeness of this paper.`):
	print(``):
	print(`Let f(x) be its (ordinary) generating function`):
	print(f(x)=Sum(A(n)*x^n,n=0..infinity)):
	print(``):
	gu:=EqualUnBoundedScoringPaths({[1,a],[1,-b]},x,f):
	gu:=subs(x=x^(1/(a+b)),gu):
	print(`Then f(x) satisfies`):
	print(``):
	print(gu=0):
	print(``):
	print(`and in Maple input format`):
	print(``):
	lprint(gu=0):
	print(``):
	print(`For the sake of the OEIS, here are the first`, K+1, `terms (starting at n=0)`):
#	mu:=taylor(RootOf(subs(f=_Z,gu)),x=0,K+1):
#	lprint([seq( coeff(mu,x,i), i=0..K)]):
#The following is faster for small-ish values of K. And the above is typically not accurate.
	lprint([seq(SpecificUnBoundedScoringPathsNumber({[1,a],[1,-b]},n1*(a+b),0),n1=0..K)]):
	if recurrence then
		print(``):
		print(`The coefficients satisfy the following recurrence:`):
		if degree(gu,f)<8 then
			R:=algtorec(gu,f,x,n,N):	#This is likely not the lowest order recurrence but it provides an upper bound of search space.
#			R:=Findrec([seq(SpecificUnBoundedScoringPathsNumber({[1,a],[1,-b]},n1*(a+b),0),n1=1..(2+degree(R,n))*(1+degree(R,N))+1)],n,N,degree(R,N)+degree(R,n)):	#This is likely a lower order recurrence.
		else
			R:=Findrec([seq(SpecificUnBoundedScoringPathsNumber({[1,a],[1,-b]},n1*(a+b),0),n1=1..(2+16)*(1+16)+1)],n,N,40):
		fi:
		if R<>FAIL then
			print(R):
			print(`and in Maple input format`):
			lprint(R):
			if asymptotics then
				print(``):
				print(`The coefficients are asymptotic to`):
				if degree(R,N)<15 then	#Arbitrary cutoff to make it not run "forever".
					EXP:=AsyC(R,n,N,0,[seq(SpecificUnBoundedScoringPathsNumber({[1,a],[1,-b]},n1*(a+b),0),n1=1..degree(R,N))],10000):
					if EXP<>FAIL then
						print(factor(simplify(evalf(EXP[1]*EXP[2])))):
					else
						print(`There was an error computing asymptotics from the recurrence.`):
						print(`The following is purely the exponential growth term.`):
						print(AsymptoticBase(gu,x,f)^n):
					fi:
				else
					print(`The following is purely the exponential growth term.`):
					print(AsymptoticBase(gu,x,f)^n):
				fi:
			fi:
		else
			print(`We were unable to find a recurrence.`):
			print(`This is because the conversion would have taken too long; the degree of the minimal polynomial was too high. (Arbitrary cutoff of 8).`):
			print(`We then tried to guess at a recurrence but the search space was not large enough. (Arbitrary cutoff of 40 for degree+order).`):
			if asymptotics then
				print(``):
				print(`The following is purely the exponential growth term.`):
				print(`The coefficients are asymptotic to`):
				print(AsymptoticBase(gu,x,f)^n):
			fi:
		fi:
	elif asymptotics then
		print(`The coefficients are asymptotic to`):
		print(`The following is purely the exponential growth term.`):
		print(AsymptoticBase(gu,x,f)^n):
	fi:
	print(``):
	print(`----------------------------`):
	print(cat(`This ends this paper that took `, time()-t0, ` seconds to generate.`)):
	print(`----------------------------------------------------------------------------------------------------------------`):
end:


###########################################################################################################################################
###########################################################################################################################################
#Algebraic to Recurrence Formula from SCHUTZENBERGER package
###########################################################################################################################################
###########################################################################################################################################
ezraSCHUTZENBERGER:=proc():
	if nops([args])=1 and op(1,[args])=`algtorec` then
		print(`algtorec(P,F,x,n,N): Inputs a polynomial P in the two variables F and x,`):
		print(`where P is the ordinary generating function of a sequence, say a(n), that satisfies the algebraic equation.`):
		print(`Outputs a recurrence equation that a(n) satisfies.`):
		print(`N is the shift operator. E.g. for the Catalan numbers when they are defined by P=1+x*P^2, the input`):
		print(`algtorec(F-1-x*F^2,F,x,n,N);`):
		print(`has output`):
		print(`4*n+2+(-n-2)*N`):
		print(`which means (4*n+2)*a(n)-(n+2)*a(n+1)=0.`):
#		print(`To rigorously prove Mathar's conjecture in OEIS  A(x) = (1 - sqrt(1 - 4*(x-x^2)/(1-x+x^2) ))/2, https://oeis.org/A090826 `):
#		print(`The first step is: `):
#		print(`ope:=algtorec(4*F^2-8*F^2*x-4*F^2*x^2-4*F+8*F^2*x^3+4*F*x+4*F^2*x^4+4*F*x^2+4*x,F,x,n,N);`):
		print(`For a recurrence of the number of nonnegative bridges with step set {[1,-2],[1,-1],[1,0],[1,1],[1,2]} type:`):
		print(`algtorec(EqualSemiBoundedScoringPaths({[1,-2],[1,-1],[1,0],[1,1],[1,2]},0,t,F),F,t,n,N);`):
	elif args=`algtodiff` then
		print(`algtodiff(F,P,x,ORDER): Using Comtet's method, given a polynomial F(P,x), finds the linear diff. eq. of order ORDER satisfied by the formal power series f(x) that satisfies the alg.eq. F(f(x),x)=0.`):
		print(`ORDER=degree(F,P) should always work, but you may be able to get away with lower order.`):
		print(`WARNING: it does not always give the MINIMAL operator`):
		print(``):
		print(`For example, for the Catalan number recurrence try: `):
		print(`algtodiff(P-1-x*P^2,P,x,2);`):
	elif args=`simp` then
		print(`A helper function for algtodiff that parses the result.`):
	elif args=`simp1` then
		print(`A helper function for algtodiff that parses the result.`):
	elif args=`difftorec` then
		print(`difftorec(ope1,D,x,n,N): Converts the differential eq. (linear)`):
		print(`with polynomial coefficients satisfied by the g.f f(x) given by the operator`):
		print(`ope(D,x), to the recurrence, in n and N satisfied`):
		print(`by the coefficients of f(x). For example`):
		print(`difftorec(x*D-1,D,x,n,N) will yield n-1`):
	elif args=`pashet` then
		print(`A helper function for difftorec that parses the result.`):
	elif args=`MyGcd` then
		print(`MyGcd(L): inputs a list of polynomials and finds their gcd. Try`):
		print(`MyGcd([n,n^2,n^3]);`):
	elif args=`SeqFromRec` then
		print(`SeqFromRec(ope,Ini,n,N,K): Given the first L-1 terms of the sequence Ini=[f(1), ..., f(L-1)]`):
		print(`satisfied by the recurrence ope(n,N)f(n)=0,`):
		print(`extends it to the first K values.`):
		print(`For example, try:`):
		print(`SeqFromRec(N-n-1,[1],n,N,10);`):
	elif args=`Yafe` then
		print(`Yafeproc(ope,N): Helper function for SeqFromRec.`):
	else
		print(`These are functions taken from the SCHUTZENBERGER package by Doron Zeilberger and Mingjia Yang.`):
		print(`SCHUTZENBERGER: A PACKAGE FOR HANDLING AND CONJECTURING AND PROVING ALGEBRAIC EQUATIONS SATISFIED BY INTEGER SEQUENCES, and RELATED STUFF.`):
		print(`The most current version of the SCHUTZENBERGER package is available at`):
		print(`http://sites.math.rutgers.edu/~zeilberg/tokhniot/SCHUTZENBERGER.txt`):
		print(``):
		print(`Type ezraSCHUTZENBERGER(function) for help:`):
		print(`algtorec, algtodiff, simp, simp1 difftorec, pashet, MyGcd, SeqFromRec, Yafe`):
	fi:
end:

###########################################################################################################################################
###########################################################################################################################################
algtorec:=proc(P,F,x,n,N) local ope,d,i,mu:
	for d from 1 to degree(P,F) do
		ope:=algtodiff(P,F,x,d):
		if ope<>FAIL then
			ope:=difftorec(ope,D,x,n,N):
			mu:=MyGcd([seq(coeff(ope,N,i),i=0..degree(ope,N))]):
			ope:=normal(ope/mu):
			ope:=add(factor(coeff(ope,N,i))*N^i,i=0..degree(ope,N)):
			RETURN(ope):   
		fi:
	od:
	FAIL:
end:
###########################################################################################################################################
algtodiff:=proc(P,F,x,ORDER) local i,gu,a,lu,pa,deg,eq1,lu1,y,var,eq,ope,pip:
	deg:=degree(P,F):
	pa:=array(deg..3*deg+1):
	eq1:=subs(F^deg=y,P):
	pa[deg]:=solve(eq1=0,y):
	pa[deg]:=subs(F=F(x),pa[deg]):
	for i from 1 to 2*deg+1 do
		lu:=pa[deg+i-1]*F(x):
		lu:=expand(lu):
		lu:=subs(F(x)^deg=pa[deg],lu):
		pa[deg+i]:=lu:
	od:

	gu:=a[0]*F(x):
	lu:=subs(F=F(x),P):
	lu1:=diff(lu,x):
	lu1:=subs(diff(F(x),x)=y,lu1):
	lu1:=solve(lu1=0,y):
	lu1:=simp(pa,deg,lu1,F,x):
	lu1:=normal(lu1):

	gu:=gu+a[1]*lu1:
	gu:=normal(gu):
	gu:=simp(pa,deg,gu,F,x):

	var:={a[0],a[1]}:
	ope:=a[0]+a[1]*D:
	lu:=lu1:
	for i from 2 to ORDER do
		lu:=diff(lu,x):
		lu:=subs(diff(F(x),x)=lu1,lu):
		lu:=normal(lu):
		lu:=simp(pa,deg,lu,F,x):
		gu:=gu+a[i]*lu:
		gu:=normal(gu):
		gu:=simp(pa,deg,gu,F,x):
		var:=var union {a[i]}:
		ope:=ope+a[i]*D^i:
	od:

	gu:=normal(gu):
	gu:=numer(gu):
	gu:=expand(gu):

	eq:={}: 
	for i from 0 to deg-1 do
		eq:=eq union {coeff( gu,F(x),i)=0}:
	od:
	var:=solve(eq,var):
	if var=NULL then
		RETURN(FAIL):
	fi:
	ope:=subs(var,ope):
	if ope=0 then
		RETURN(FAIL):
	fi:

	ope:=normal(ope):
	ope:=numer(ope):
	ope:=expand(ope):
	pip:=coeff(ope,D,degree(ope,D)):
	pip:=coeff(pip,x,degree(pip,x)):
	normal(ope/pip):
end:
##################################################################
simp:=proc(pa,deg,bitui,F,x) local mone,mekh,gu:
	gu:=normal(bitui):
	mone:=numer(gu):
	mekh:=denom(gu):
	mone:=simp1(pa,deg,mone,F,x):
	mekh:=simp1(pa,deg,mekh,F,x):
	expand(mone)/expand(mekh):
end:
##################################################################
simp1:=proc(pa,deg,bitui,F,x) local i,gu:
	gu:=expand(bitui):
	for i from deg to 3*deg+1 do
		gu:=subs(F(x)^i=pa[i],gu):
	od:
	gu:
end:
###########################################################################################################################################
difftorec:=proc(ope1,D,x,n,N) local i,j,gu,ord,ope,r,mu:
	ope:=expand(ope1):
	gu:=0:
	ord:=degree(ope,D):
	for i from 0 to  ord do
		mu:=coeff(ope,D,i):
		for j from 0 to degree(mu,x) do
			gu:=gu+coeff(mu,x,j)*product(n+r-j,r=1..i)*N^(i-j):
		od:
	od:
	pashet(gu,n,N):
end:
##################################################################
pashet:=proc(p,n,N) local i,gu1,gu,p1,ra:
	p1:=normal(p):
	gu1:=denom(p1):
	ra:=degree(gu1,N):
	p1:=subs(n=n+ra,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):
end:
###########################################################################################################################################
MyGcd:=proc(L) local gu:
	if L=[] then
		RETURN(FAIL):
	elif nops(L)=1 then
		RETURN(L[1]):
	else
		gu:=MyGcd([op(1..nops(L)-1,L)]):
		RETURN(gcd(gu,L[nops(L)])):
	fi:
end:
###########################################################################################################################################
SeqFromRec:=proc(ope,Ini,n,N,K) local ope1,gu,L,n1,j1:
	ope1:=Yafe(ope,N)[2]:
	L:=degree(ope1,N):
	if nops(Ini)<>L then
		ERROR(`Ini should be of length`, L):
	fi:

	ope1:=expand(subs(n=n-L,ope1)/N^L):
	gu:=Ini:
	for n1 from nops(Ini)+1 to K do
		gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),j1=1..L)]:
	od:
	gu:
end:
##################################################################
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:


###########################################################################################################################################
###########################################################################################################################################
#Finding Recurrence Formula from EKHAD package
###########################################################################################################################################
###########################################################################################################################################
ezraEKHAD:=proc():
	if args=`AZd` then
		print(`AZd(A,x,n,N) finds a recurrence of order ORDER<=8`):
		print(`phrased in terms of n, and the shift in n, N for the integral of A(x,n) with respect to x,`):
		print(`whenever A(x,n) is hypergeometric in (x,n),i.e. A(x,n+1)/A(x,n) and A'(x,n)/A(x,n) are rational functions of x and n.`):
		print(`It follows the algorithm of Gert Almkvist and Doron Zeilberger, "The method of differentiating under the integral sign", J. Symbolic Computation 10(1990), 571-591`):
		print(`A should not be a product of a function of x and a function of n.`):
		print(`Also produces a proof certificate.`):
		print(``):
		print(`AZd(A,x,n,N) returns the expression seq. ope(n,N),cert(x,n)`):
		print(`satisfying ope(n,N)A(x,n)=diff(cert(x,n)*A(x,n),x).`):
		print(`If no recurrence is found, it returns 0.`):
		print(``):
		print(`For example for the recurrence equation and proof certificate for the Legendre polynomials. Try: `):
		print(`AZd(1/sqrt(1-2*x*t+t^2)/t^(n+1),t,n,N);`):
		print(`For a recurrence of the number of unbounded bridges with step set {[1,-2],[1,-1],[1,0],[1,1],[1,2]} type:`):
		print(`AZd((t^(-2)+t^(-1)+1+t+t^2)^n/t,t,n,N);`):
	elif args=`duisd` then
		print(`Actually does the work of AZd.`):
	elif args=`pashetd` then
		print(`A helper function for duisd that parses the result.`):
	elif args=`goremd` then
		print(`A helper function for duisd that parses the result.`):
	else
		print(`These are functions taken from the EKHAD package by Doron Zeilberger.`):
		print(`A Maple package for proving Hypergeometric (Binomial Coeff.) and other kinds of identities.`):
		print(`The most current version of the EKHAD package is available at`):
		print(`http://sites.math.rutgers.edu/~zeilberg/tokhniot/EKHAD.txt`):
		print(``):
		print(`Type ezraEKHAD(function) for help:`):
		print(`AZd, duisd, goremd, pashetd`):
	fi:
end:

###########################################################################################################################################
###########################################################################################################################################
AZd:=proc(A,y,n,N) local ORDER,gu,KAMA:
	KAMA:=8:
	for ORDER from 0 to KAMA do
		gu:=duisd(A,ORDER,y,n,N):
		if gu<>0 then
			RETURN(gu):
		fi:
	od:
	0:
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:=SolveTools[Linear](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:
################################################################
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:
################################################################
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:


###########################################################################################################################################
###########################################################################################################################################
#Finding Algebraic Equation from W1D package
###########################################################################################################################################
###########################################################################################################################################
ezraW1D:=proc():
	if args=`Empir` then
		print(`Empir(gu,x,P): empirically finds an algebraic equation`):
		print(`F(P(x),x)=0`):
		print(`for the formal power series P(x):=sum_{i=0} gu[i]*x^i, where gu is a list.`):
		print(`For example`):
		print(`Empir([seq(1,i=1..20)],x,P)`):
		print(`should yield`):
		print(`P-xP-1=0.`):
		print(`If there is not enough data, it returns 0. You should then try it again with a longer sequence. For example, try:`):
		print(`Empir([seq(binomial(2*i,i)/(i+1),i=0..20)],x,P);`):
	elif args=`Empir1` then
		print(`Empir1(gu,x,P): Helps Empir by checking algebraic equations of varying degree in x and P.`):
	elif args=`empir` then
		print(`empir(gu,degx,degP,x,P): Helps Empire by checking an algebraic equation for specific degx and degP.`):
	elif args=`EmpirF` then
		print(`EmpirF(gu,x,P): a possibly faster version of Empir using the gfun package`):
		print(`For example, try:`):
		print(`EmpirF([seq(binomial(2*i,i)/(i+1),i=0..20)],x,P);`):
	else
		print(`These are functions taken from the W1D package by Doron Zeilberger.`):
		print(`A Maple package for the enumeration and probability of 1 dimensional walks restricted to the positive half line.`):
		print(`The most current version of the W1D package is available at`):
		print(`http://sites.math.rutgers.edu/~zeilberg/tokhniot/W1D.txt`):
		print(``):
		print(`Type ezraW1D(function) for help. Main functions:`):
		print(`Empir, EmpirF.`):
		print(`Helper functions:`):
		print(`Empir1, empir`):
		print(``):
		print(`Since this is empirically done, you may want to increase the degree bound in the gfun package.`):
#		print(`The bound has already been increased to gfun[maxdegeqn]:=14:`):
	fi:
end:
#gfun[maxdegeqn]:=14:

###########################################################################################################################################
###########################################################################################################################################
Empir:=proc(gu,x,P) local i,lu,che:
	for i from 10 by 10 to 10*trunc(nops(gu)/10) do
		lu:=Empir1([op(1..i,gu)],x,P):
		if lu<>FAIL then
			che:=expand(taylor(subs(P=add(gu[i+1]*x^i,i=0..nops(gu)-1),lu),x=0,nops(gu)-1)):
			if {seq(coeff(che,x,i),i=0..nops(gu)-3)}={0} then
				RETURN(lu):
			fi:
		fi:
	od:

	lu:=Empir1(gu,x,P):
	che:=expand(taylor(subs(P=add(gu[i+1]*x^i,i=0..nops(gu)-1),lu),x=0,nops(gu)-1)):
	if {seq(coeff(che,x,i),i=0..nops(gu)-3)}<>{0} then
		RETURN(FAIL):
	else
		RETURN(lu):
	fi:
end:
###########################################################################################################################################
Empir1:=proc(gu,x,P) local degx,degP,L,lu,che,i:
	for L from 1 to (nops(gu)-3)/3 do
		for degP from 1 to L do
			for degx from 0 to min(trunc((nops(gu)-3)/(1+degP))-1,L-degP) do
				lu:=empir(gu,degx,degP,x,P):
				if lu<>FAIL then
					che:=expand(taylor(subs(P=add(gu[i+1]*x^i,i=0..nops(gu)-1),lu),x=0,nops(gu)-1)):
					if {seq(coeff(che,x,i),i=0..nops(gu)-3)}<>{0} then
						RETURN(FAIL):
					else
						RETURN(lu):
					fi:
				fi:
			od:
		od:
	od:
	FAIL:
end:
###########################################################################################################################################
empir:=proc(gu,degx,degP,x,P) local i1,i2,F,a,cand,lu,eq,var,mu,flo,pip,var2,mu1,halev,vu,i,vu1:
	if (1+degx)*(1+degP) > nops(gu)-3 then
		RETURN(`sequence too small`):
	fi:

	F:=0:
	var:={}:
	for i1 from 0 to degx do
		for i2 from 0 to degP do
			F:=F+a[i1,i2]*x^i1*P^i2:
			var:=var union {a[i1,i2]}:
		od:
	od:

	cand:=0: 
	for i1 from 0 to nops(gu)-1 do
		cand:=cand+op(i1+1,gu)*x^i1:
	od:

	lu:=subs(P=cand,F):
	lu:=taylor(lu,x=0,nops(gu)-1):
	eq:={}:
	for i1 from 0 to nops(gu)-2 do
		eq:=eq union {coeff(lu,x,i1)=0}
	od:
	mu:=solve(eq,var):
	F:=subs(mu,F):
	if F=0 then
		RETURN(FAIL):
	fi:

	var2:={}:
	for mu1 in mu do
		if op(1,mu1)=op(2,mu1) then
			var2:=var2 union {op(1,mu1)}:
		fi:
	od:
	if nops(var2)<>1 then
		F:=coeff(F,var2[1],1):
	fi:

	flo:=degree(F,P):
	pip:=coeff(F,P,flo):
	flo:=degree(pip,x):
	pip:=coeff(pip,x,flo):
	F:=F/pip:
	F:=numer(normal(F)):
	vu:=add(gu[i+1]*x^i,i=0..nops(gu)-1):
	halev:=add(factor(coeff(F,P,i))*P^i,i=0..degree(F,P)):
	vu1:=taylor(subs(P=vu,halev),x=0,nops(gu)+1):
	if {seq(expand(coeff(vu1,x,i)),i=1..nops(gu)-1)}<>{0} then
		RETURN(FAIL):
	fi:
	halev:
end:
###########################################################################################################################################
EmpirF:=proc(gu,t,P) local eq1,eq2,vu,lu,i:
	if nops(gu)<20 then
		print(`list is too short `):
		RETURN(FAIL):
	fi:

	eq1:=gfun[listtoalgeq]([op(1..nops(gu)-10,gu)],P(t),[ogf]):
	if eq1=FAIL then
		RETURN(FAIL):
	fi:
	eq1:=subs(P(t)=P,eq1[1]):

	eq2:=gfun[listtoalgeq](gu,P(t),[ogf]):
	if eq2=FAIL then
		RETURN(FAIL):
	fi:
	eq2:=subs(P(t)=P,eq2[1]):

	if eq1<>eq2 then
		RETURN(FAIL):
	fi:

	vu:=add(gu[i+1]*t^i,i=0..nops(gu)-1):
	lu:=taylor(subs(P=vu,eq2),t=0,nops(gu)+4):
	if {seq(expand(coeff(lu,t,i)),i=1..nops(gu)-2)}<>{0} then
		RETURN(FAIL):
	fi:
	eq2:
end:


###########################################################################################################################################
###########################################################################################################################################
#Finding recurrences by guessing and asymptotic from recurrences. From package AsyRec by Doron Zeilberger.
###########################################################################################################################################
###########################################################################################################################################
ezraAsyRec:=proc():
	if nargs=1 and args[1]=Asy then
		print(`Asy(ope,n,N,M): the asymptotic expansion of solutions `):
		print(`to ope(n,N)f(n)=0,where ope(n,N) is  a recurrence operator`):
		print(`up to the M's term`):
		print(`For example, try:`):
		print(`Asy((n+2)^3*N^2-(34*n^3+153*n^2+231*n+117)*N+(n+1)^3,n,N,3);`):
		print(`Asy((n+2)^2*N^2-(7*n^2+21*n+16)*N-8*(n+1)^2,n,N,3);`):
		print(`Asy(N^2-N-(n+1),n,N,3);`):
		print(`Asy(N^3-N^2-(n+2)*N-(n+2)*(n+1),n,N,3);`):
	elif nargs=1 and args[1]=AsyC then
		print(`AsyC(ope,n,N,M,Ini,K): the asymptotic expansion of solutions `):
		print(`to ope(n,N)f(n)=0, with the given initial conditions`):
		print(`Starting at n=1 (NOT at n=0)`):
		print(`where ope(n,N) is  a recurrence operator`):
		print(`up to the M's term`):
		print(`and complete with an empirically derived constant in front`):
		print(`using K terms `):
		print(`For example, try:`):
		print(`AsyC((n+2)^3*N^2-(34*n^3+153*n^2+231*n+117)*N+(n+1)^3,n,N,5,[5,73],1000);`):
		print(`AsyC((n+2)^2*N^2-(7*n^2+21*n+16)*N-8*(n+1)^2,n,N,5,[2,10],1000);`):
		print(`AsyC(N^2-N-(n+1),n,N,5,[1,2],1000);`):

	elif nargs=1 and args[1]=Asy1special then
		print(`Asy1special(ope,n,N,K,x): the asymptotic expansion of solutions `):
		print(`to ope(n,N)f(n)=0,given as a list`):
		print(`[pu,lu,expansion,r] where it is`):
		print(`exp(pu)*lu^n*expansion(x), where x=1/n^(1/r), and r is`):
		print(`a positive integer. It also returns [K,k] (see Asy1)`):
		print(`where ope(n,N) is  a recurrence operator`):
		print(`up to the K's term`):
		print(`For example, try: `):
		print(`Asy1special((n+2)^2*N^2-(7*n^2+21*n+16)*N-8*(n+1)^2,n,N,3,x);`):
		print(`Asy1special(N^2-N-(n+1),n,N,3,x);`):
	elif nargs=1 and args[1]=Atom then
		print(`Atom(s,r,i,x,K): Expanding (n+i)^(s/r)-n^(s/r) in terms of`):
		print(`x=n^(-1/r) up to K terms`):
		print(`For example try:`):
		print(`Atom(2,3,2,x,3);`):
	elif nargs=1 and args[1]=CODV then
		print(`CODV(ope,n,N,k,K): Given a linear recurrence operator ope(n,N)`):
		print(`annihilating a[n], say, outputs the operator`):
		print(`annihilating b[n]:=a[n*K]/n!^k . `):
		print(`For example, try: CODV(N-(n+1),n,N,1,1):`):
	elif nargs=1 and args[1]=CODV1 then
		print(`CODV1(ope,n,N,k): Given a linear recurrence operator ope(n,N)`):
		print(`annihilating a[n], say, outputs the operator`):
		print(`annihilating b[n]:=a[n]/n!^k`):
		print(`For example, try: CODV1(N-(n+1),n,N,1):`):
	elif nargs=1 and args[1]=Finda then
		print(`Finda(ope,N,x,r): finds the first power x^a in the`):
		print(`asymptotic solution of ope(N,n)f(n)=0, where x=1/n^(1/r)`):
		print(`For example, try:`):
		print(`Finda((1+x)-(1+3*x)*N,N,x,1);`):
	elif nargs=1 and args[1]=FindExpP then
		print(`FindExpP(ope,n,N): finds the exponential part of the asymptotics`):
		print(`for the solutions of ope(n,N)a(n)=0 if it is normalized`):
		print(`such that the leading asymp. is 1^n`):
		print(`for example, try:`):
		print(`FindExpP((n+1)*N^2-2*(n+5)*N+n+3,n,N);`):
	elif nargs=1 and args[1]=FindExpP1 then
		print(`FindExpP1(ope,n,N,r,x): finds the exponential part of the asymptotics`):
		print(`in terms of x=n^(1/r)`):
		print(`for the solutions of ope(n,N)a(n)=0 if it is of type r.`):
		print(`for example, try:`):
		print(`FindExpP1((n+1)*N^2-2*(n+5)*N+n+3,n,N,2,x);`):
	elif nargs=1 and args[1]=FindKk then
		print(`FindKk(ope,n,N): Given a linear recurrence operator with polynomial`):
		print(`coefficients, ope(n,N), finds the integers K and k such`):
		print(`that if a(n) is a solution of ope(n,N)a(n)=0, then`):
		print(`b(n):=a(K*n)/n!^k is annihilated by a standard operator`):
		print(`the output is the pair [k,K] and the transformed operator`):
		print(`For exanple, try: FindKk(N^2-n,n,N);`):
	elif nargs=1 and args[1]=Ksect then
		print(`Ksect(ope,n,N,k,r): Given a linear recurrence operator ope(n,N)`):
		print(`annihilating a sequence a[n], and pos. integer k, and`):
		print(`non-neg. integer r, outputs the one, of the same`):
		print(`order, that annihilates  a[k*n+r]`):
		print(`For example try: Ksect(N^2-N-(n+1),n,N,2,0);`):
	elif nargs=1 and args[1]=NakedStirling then 
		print(`NakedStirling(n,K): the asymptotic expansion of`):
		print(`n!/((n/e)^n*sqrt(2*Pi*n). For example, try:`):
		print(`NakedStirling(n,5);`):
	elif nargs=1 and args[1]=NewOpe then
		print(`NewOpe(ope,n,N,K,x): the exponential part + the transformed operator`):
		print(`(up to degree K asymp. in the coefficients), in terms of x=1/n^(1/r)`):
		print(`For example, try:`):
		print(`NewOpe((n+1)*N^2-2*(n+5)*N+n+3,n,N,4,x);`):
	elif nargs=1 and args[1]=Nor then
		print(`Nor(ope,n,N): the Normalizer of the linear recurrence`):
		print(`operator with polynomial coefficients ope(n,N)`):
		print(`followed by its exponential growth constant`):
		print(`For example, try:`):
		print(`Nor(N^2-3*N+2,n,N);`):
	elif nargs=1 and args[1]=OneStepG then
		print(`OneStepG(ope,N,x,r,Cu): Given the asymptotic expansion of`):
		print(`a solution to ope(n,N) expressed in terms of x=1/n^(1/r)`):
		print(`finds the next one`):
		print(`finds one more term`):
		print(`OneStepG((1+x)-(1+3*x)*N,N,x,1,x^2);`):
	elif nargs=1 and args[1]=SeqFromRec then
		ezraSCHUTZENBERGER(SeqFromRec):
	elif nargs=1 and args[1]=Findrec then
		print(`Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with poly cofficients`):
		print(`of minimum ORDER such that DEGREE+ORDER<=MaxC.`):
		print(`e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);`):
	elif nargs=1 and args[1]=findrec then
		print(`findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating`):
		print(`the sequence f of degree DEGREE and order ORDER.`):
		print(`For example, try: findrec([seq(i,i=1..10)],0,2,n,N);`):

	else
		print(`These are functions taken from the AsyRec package by Doron Zeilberger.`):
		print(`AsyRec: A Maple package for finding the asymptotics of solutions of (homog.) linear recurrence equations with polynomial coefficents using (a variant) of Birkhoff-Trjitzinsky method.`):
		print(`The most current version of the AsyRec package is available at`):
		print(`http://sites.math.rutgers.edu/~zeilberg/tokhniot/AsyRec.txt`):
		print(``):
		print(`Type ezraAsyRec(function) for help. Main functions:`):
		print(`Asy, AsyC`):
		print(``):
		print(`Helper functions:`):
		print(`Asy1, Asy1special, Atom, CODV, CODV1, Finda, FindExpP, FindExpP1, FindKk, Ksect, NakedStirling, NewOpe, Nor, OneStepG, SeqFromRec, Findrec, findrec.`):
	fi:
end:
with(SolveTools):	#Needed for Linear()

###########################################################################################################################################
###########################################################################################################################################
#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,evalf(Ini),n,N,K):
	mu:=[seq(evalf(L[i]/subs(n=i,gu)),i=K-10..K)]:
	if abs(mu[nops(mu)])<10^(-7) then
		print(`Something is fishy`):
		RETURN(FAIL):
	fi:

	er:=mu[nops(mu)]-mu[nops(mu)-1]:
	if er=0 then
		D1:=Digits:
	else
		D1:=-trunc(log(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):
	identify(C),gu:
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:
#########################################################################################
#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 from another package not included here)
#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:
#########################################################################################
#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 example, 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 polynomial coefficients 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:
#########################################################################################
#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:
#########################################################################################
#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:
#########################################################################################
rf:=proc(a,n) local i:
	mul(a+i,i=0..n-1):
end:
#########################################################################################
#Ksect(ope,n,N,k,r): Given a linear recurrence operator ope(n,N)
#annihilating a sequence a[n], and pos. integer k, and
#non-neg. integer r, outputs the one, of the same
#order, that annihilates  a[k*n+r]
#For example try: Ksect(N^2-N-(n+1),n,N,2,0);
Ksect:=proc(ope,n,N,K,r) local ORDER,eq,var,lu,mu,c,Ope,var1,i,vu,v:
	ORDER:=degree(ope,N):
	for i from 0 to ORDER do
		lu[i]:=MonoN(ope,n,N,r+K*i):
		lu[i]:=subs(n=n*K,lu[i]):
	od:
	Ope:=add(c[i]*N^i,i=0..ORDER):
	mu:=expand(add(c[i]*lu[i],i=0..ORDER)):
	eq:={seq(coeff(mu,N,i),i=0..ORDER)}:
	var:={seq(c[i],i=1..ORDER)}:
	var1:=solve(eq,var):
	if var1=NULL then
		RETURN(FAIL):
	fi:

	vu:={}:
	for i from 1 to nops(var1) do
		if op(1,var1[i])=op(2,var1[i]) then
			vu:=vu union {op(1,var1[i])}:
		fi:
	od:
	Ope:=subs(var1,Ope):
	for v in vu do
		Ope:=subs(v=0,Ope):
	od:
	Ope:=subs(c[0]=1,numer(normal(Ope))):
end:
#########################################################################################
#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 k<ORDER then
		RETURN(N^k):
	fi:
	if k=ORDER then
		coe0:=coeff(ope,N,ORDER):
		RETURN(add(normal(-coeff(ope,N,i)/coe0)*N^i,i=0..ORDER-1)):
	fi:
	lu1:=expand(MonoN(ope,n,N,k-1)):
	lu2:=expand(MonoN(ope,n,N,ORDER)):
	lu1:=expand(subs(n=n+1,lu1)*N):
	lu1:=subs(N^ORDER=lu2,lu1):
	add(normal(coeff(lu1,N,i))*N^i,i=0..ORDER-1):
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:
#########################################################################################
#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:
#########################################################################################
#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:
#########################################################################################
#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:
#########################################################################################
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:
#########################################################################################
#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:
#########################################################################################
#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:
#########################################################################################
#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:
#########################################################################################
#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:
#########################################################################################
#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:=pashetA(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:
#########################################################################################
pashetA:=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:
#########################################################################################
#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:
#########################################################################################
#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:
###########################################################################################################################################
#Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with polynomial cofficients
#of minimum ORDER such that DEGREE+ORDER<=MaxC.
#e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);
Findrec:=proc(f,n,N,MaxC) local DEGREE, ORDER,ope,L:
	for L from 0 to MaxC do
		for ORDER from 0 to L do
			DEGREE:=L-ORDER:
			if (2+DEGREE)*(1+ORDER)+1>nops(f) then
				print(cat(`Insufficient data for degree `,DEGREE,` and order `,ORDER)):
				RETURN(FAIL):
			fi:
			ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+1,f)],DEGREE,ORDER,n,N):
			if ope<>0 then
				RETURN(ope):
			fi:
		od:
	od:
	FAIL:
end:
###############################################################################################
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrec:=proc(f,DEGREE,ORDER,n,N) local ope,var,eq,a,i,j,n0,kv,var1,eq1,mu,commonDenom:
	if (1+DEGREE)*(1+ORDER)+2+ORDER>nops(f) then
		ERROR(cat(`Insufficient data for a recurrence of order `,ORDER,` degree `,DEGREE)):
	fi:

	ope:=0:
	var:={}: 
	for i from 0 to ORDER do
		for j from 0 to DEGREE do
			ope:=ope+a[i,j]*n^j*N^i:
			var:=var union {a[i,j]}:
		od:
	od:

	eq:={}:
	for n0 from 1 to (1+DEGREE)*(1+ORDER)+2 do
		eq1:=0:
		for i from 0 to ORDER do
			eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
		od:
		eq:= eq union {eq1}:
	od:

	var1:=solve(eq,var):
	ope:=subs(var1,ope):
	if ope=0 then
		RETURN(0):
	fi:

	kv:={}:
	for i from 1 to nops(var1) do
		mu:=op(i,var1):
		if op(1,mu)=op(2,mu) then
			kv:= kv union {op(1,mu)}:
		fi:
	od:
	if nops(kv)>1 then
		print(`The output is possibly not the minimal operator.`):
		print(`Either DEGREE or ORDER may be too high.`):
	fi:
	ope:=subs({kv[1]=1,seq(kv[i]=0,i=2..nops(kv))},ope):
	commonDenom:=ilcm(seq(abs(denom(d)),d in convert(ope,set))):
	add(factor(coeff(ope,N,i)*commonDenom)*N^i,i=0..degree(ope,N)):
end: