(* ::Package:: *)

(* ::Text:: *)
(*This package was written by Luis A. Medina and Doron Zeilberger.*)
(*It is one of the packages that accompany the article :*)
(*'An Experimental Mathematics Perspective on the Old, and still Open, Question of When To Stop?'*)
(*	by Luis A. Medina and Doron Zeilberger.*)
(**)
(* Please report bug to : lmedina at math dot rutgers dot edu.*)
(* *)
(*The most current version of this  package and paper*)
(*are  available from http://www.math.rutgers.edu/~lmedina.*)
(**)
(* *)
(* Created: June 20, 2009.*)


BeginPackage["GuessF`"];


Mensaje::usage=" Created: June 20, 2009.

This program was written by Luis A. Medina and Doron Zeilberger,
Rutgers University. [medina, zeilberg] at math dot rutgers dot edu 

This is one of the packages that accompany the article:
'An Experimental Mathematics Perspective on the Old, and 
still Open, Question of When To Stop?'
	by Luis A. Medina and Doron Zeilberger.

Please report bug to: lmedina at math dot rutgers dot edu.

The most current version of this  package and paper
are  available from http://www.math.rutgers.edu/~lmedina.

Functions included in this package are:

GuessRationalFunction[list,m,n,x]
GF[m,n,\[Alpha]]
BUILDER[m,n,\[Alpha]].

If you want to know more information about these functions, 
simpliy input ?nameoffunction.  For example, try ?GF."


GuessRationalFunction::usage="GuessRationalFunction[list,degN,degD,x] tries to match 
a given list of values to a rational function in x whose degree 
in the numerator is at most degN and whose degree in the 
denominator is at most degN."


GF::usage="GF[m,n,\[Alpha]] takes as input a positive integer m and two variables n and \[Alpha], 
It then tries to guess the closed forms of \!\(\*StyleBox[\"f\",\nFontSlant->\"Italic\"]\)(n+\[Alpha],n-\[Alpha]-m+1,2n+1).
For instance, try GF[3,n,\[Alpha]].

GF has an option called Bound.  What Bound does is to 
change the bound on 'n'. For instance, GF[2,n,\[Alpha],Bound->20] 
assumes that n>=20."


BUILDER::usage="BUILDER[m,n,\[Alpha]] takes as input a positive integer m and
two variables n and \[Alpha].  It then find the closed form of \!\(\*StyleBox[\"f\",\nFontSlant->\"Italic\"]\)(n+\[Alpha],n-\[Alpha]-m+1,2n+1).
For instance, try BUILDER[3,n,\[Alpha]].

BUILDER has an option called Bound.  What Bound does is 
to change the bound on 'n'. For instance, BUILDER[2,n,\[Alpha],Bound->20] 
assumes that n>=20."


Bound::usage = "Bound-> 'natural number' is an option for GF[m,n,\[Alpha]] and 
BUILDER[m,n,\[Alpha]]. What it does is to change the bounds on 'n'.  
For instance, GF[2,n,\[Alpha],Bound->20] assumes that n>=20."


Begin["Private`"];


(* ::Text:: *)
(*Error messages.*)


GuessRationalFunction::numelemlist="The number of elements in the list should exceeds the sum of the degrees.";


GF::bound="Try a bigger bound in the last argument."
BUILDER::bound="Try a bigger bound in the last argument."


(* ::Text:: *)
(*____________________________________________________________________________*)
(*AUXILIARY FUNCTION*)
(*____________________________________________________________________________*)


(* ::Text:: *)
(*It takes a rational function in x and returns the same rational function,*)
(*but with its numerator and denominator expanded.*)


NormalRational[R_,x_]:=If[Exponent[Denominator[R],x]==0,Collect[R,x],Collect[Numerator[R],x]/Collect[Denominator[R],x]]


(* ::Text:: *)
(*____________________________________________________________________________*)
(*AUXILIARY FUNCTION*)
(*____________________________________________________________________________*)


(* ::Text:: *)
(*Recursive definition of f. *)


f[h_,t_,n_]/;(h+t==n):=Max[1/2,h/n]


f[0,0,n_]:=1/2 (f[1,0,n]+f[0,1,n])


f[h_,t_,n_]/;h+t<n:=f[h,t,n]=Block[{$RecursionLimit=Infinity},Max[h/(h+t),1/2 (f[h+1,t,n]+f[h,t+1,n])]]


(* ::Text:: *)
(*____________________________________________________________________________*)
(*AUXILIARY FUNCTION*)
(*____________________________________________________________________________*)


(* ::Text:: *)
(*Guess bound for GF and BUILDER *)


A[1]=1;
R[1]=0;
R[2]=1;
R[3]=6;
R[n_Integer]/;n>=4:=Block[{$RecursionLimit=Infinity},2R[n-1]+n-3];


A[n_Integer]/;n>=2:=Block[{$RecursionLimit=Infinity},2A[n-1]+R[n]];


(* ::Text:: *)
(*____________________________________________________________________________*)
(*AUXILIARY FUNCTION*)
(*____________________________________________________________________________*)


(* ::Text:: *)
(*Version of f to work with. *)


F[m_,n_,\[Alpha]_]:=f[n+\[Alpha],n-m+1-\[Alpha],2n+1]


(* ::Text:: *)
(*____________________________________________________________________________*)
(*AUXILIARY FUNCTION*)
(*____________________________________________________________________________*)


(* ::Text:: *)
(*Options for GF and BUILDER.*)


Options[GF]={Bound->-1};
Options[BUILDER]={Bound->-1};


(* ::Text:: *)
(*____________________________________________________________________________*)
(*ONE OF THE MAIN FUNCTION IN THIS PACKAGE!*)
(*____________________________________________________________________________*)


(* ::Text:: *)
(*Guess rational function in one variable.*)
(*	Input:  listValues: a list of rational numbers.*)
(*		degN: a positive integer.*)
(*		degD: a positive integer.*)
(*		x: variable.*)
(*	Output: A rational function in x that matches the values in listValues. *)
(*	    	The  degree in the numerator is at most degN and degree in the *)
(*	    	denominator at most degD.*)


GuessRationalFunction[listValues_List,degN_Integer,degD_Integer,x_Symbol]/;
If[TrueQ[Length[listValues]>= degN+degD+1],True,Message[GuessRationalFunction::numelemlist];False]&&degN>=0&&degD>=0:=
Module[{a,b,G,vars},
G[y_]:=Total[Table[a[i] y^i,{i,0,degN}]]/Total[Table[b[i] y^i,{i,0,degD}]];
vars=Join[Table[a[i],{i,0,degN}],Table[b[i],{i,0,degD}]];
vars=Flatten[Solve[Table[Numerator[listValues[[i]]]Denominator[G[i]]==Numerator[G[i]]Denominator[listValues[[i]]],{i,1,Length[listValues]}],vars]//Quiet];
If[(Denominator[G[x]]/.vars)===0 || vars==={},$Failed,Simplify[G[x]/.vars]]
]


(* ::Text:: *)
(*____________________________________________________________________________*)
(*ONE OF THE MAIN FUNCTION IN THIS PACKAGE!*)
(*____________________________________________________________________________*)


(* ::Text:: *)
(*GF: guess value of f(n+\[Alpha],n-\[Alpha]-m+1,2n+1).*)
(*	Input:  m: an integer.*)
(*		n: a variable.*)
(*		\[Alpha]: a variable.*)
(*	Output: A guess rational function in 'n' for f(n+\[Alpha],n-\[Alpha]-m+1,2n+1).*)


GF[1,n_,\[Alpha]_,OptionsPattern[]]:=Module[{bound},
 bound = OptionValue[Bound];
 If[bound<=0 || !Element[bound,Integers]||Head[bound]==Symbol,
 bound=1
 ];
 Piecewise[{{1/2,\[Alpha]<=-1},{NormalRational[Simplify[GuessRationalFunction[Table[F[1,i,0],{i,1,20}],3,3,n]],n],\[Alpha]==0},{(n+\[Alpha])/(2n),\[Alpha]>=1}}]
]


GF[m_Integer,n_,\[Alpha]_,OptionsPattern[]]/;m>0:=Module[{branch,c=-m,func,bound},
 bound = OptionValue[Bound];
 If[bound<=0 ||!Element[bound,Integers]||Head[bound]==Symbol,
 bound = A[m];
 ];
 branch[-m]=1/2;
 If[0< bound <m,
  Message[GF::bound];$Failed,
  While[(c=c+1)<=0,
   branch[c]=NormalRational[Simplify[GuessRationalFunction[Table[F[m,i,c],{i,bound,bound+2m+5}],m,m,n]/.n->n-bound+1],n]
  ];
  branch[1]=(n+\[Alpha])/(2n+1-m);
  func=Piecewise[Join[{{branch[-m],\[Alpha]<=-m}},Table[{branch[i],\[Alpha]==i},{i,-m+1,0}],{{branch[1],\[Alpha]>=1}}]];
  If[Head[func/.\[Alpha]->0]===Max || Exponent[Numerator[func/.\[Alpha]->0],n]=!=m,
  Message[GF::bound];$Failed,
  Print["The following function was guessed assuming n\[GreaterEqual]",bound,"."];
  func
  ]
 ]
]


(* ::Text:: *)
(*____________________________________________________________________________*)
(*ONE OF THE MAIN FUNCTION IN THIS PACKAGE!*)
(*____________________________________________________________________________*)


(* ::Text:: *)
(*The next function is the builder of the function f(n+\[Alpha],n-\[Alpha]-m+1,2n+1) for n big enough.*)
(*	Input:  m: an integer.*)
(*		n: a variable.*)
(*		\[Alpha]: a variable.*)
(*	Output: A rational function in 'n' for f(n+\[Alpha],n-\[Alpha]-m+1,2n+1).*)


BUILDER[m_Integer,n_,\[Alpha]_,OptionsPattern[]]/;m>0:=Module[{bound,B},
 B = A[m];
 bound = OptionValue[Bound];
 If[bound <=0 || !Element[n,Integers]||Head[bound]==Symbol,
  bound= B
 ];
 If[0< bound <m,
  Message[BUILDER::bound];$Failed,
  If[Head[TheBuilder[m,n,\[Alpha],bound]/.\[Alpha]->0]===Max,
    Message[BUILDER::bound];$Failed,
    If[bound==B,
     Print["The following formula was found assuming n\[GreaterEqual]",B,"."];
     Print["We do not claim this bound on n is minimal!"]; 
     TheBuilder[m,n,\[Alpha],bound],
     Print["The following formula was found assuming n\[GreaterEqual]",bound,"."];
     TheBuilder[m,n,\[Alpha],bound]
    ]
  ]
 ]
]


(* ::Text:: *)
(*____________________________________________________________________________*)
(*AUXILIARY FUNCTION.*)
(*____________________________________________________________________________*)


(* ::Text:: *)
(*TheBuilder is the skeleton of BUILDER.*)


TheBuilder[1,n_,\[Alpha]_,bound_]:=
 Piecewise[{{1/2,\[Alpha]<=-1},{(3+4n)/(4+8n),\[Alpha]==0},{(n+\[Alpha])/(2n),\[Alpha]>=1}}]


TheBuilder[m_Integer,n_,\[Alpha]_,bound_]/;m>0:=TheBuilder[m,n,\[Alpha],bound]=Block[{$RecursionLimit=Infinity},
Module[{c=-m,branch},
 branch[-m]=NormalRational[Together[Assuming[Element[n,Integers] && n>=bound && \[Alpha]<=-m,Refine[Max[(n+\[Alpha])/(2n+1-m),1/2 (TheBuilder[m-1,n,\[Alpha]+1,bound]+TheBuilder[m-1,n,\[Alpha],bound])]]]],n];
 While[(c=c+1)<=0,
 branch[c]=NormalRational[Together[Assuming[Element[n,Integers] && n>= bound && Element[\[Alpha],Integers]&& \[Alpha]==c,Refine[Max[(n+\[Alpha])/(2n+1-m),1/2 (TheBuilder[m-1,n,\[Alpha]+1,bound]+TheBuilder[m-1,n,\[Alpha],bound])]]]],n]
 ];
 branch[1]=NormalRational[Together[Assuming[Element[n,Integers] && n>= bound && \[Alpha]>=1,Refine[Max[(n+\[Alpha])/(2n+1-m),1/2 (TheBuilder[m-1,n,\[Alpha]+1,bound]+TheBuilder[m-1,n,\[Alpha],bound])]]]],n];
 Piecewise[Join[{{branch[-m],\[Alpha]<=-m}},Table[{branch[i],\[Alpha]==i},{i,-m+1,0}],{{branch[1],\[Alpha]>=1}}]]
 ]
]//Quiet


?Mensaje


End[];
EndPackage[];