# File     : GUESS
# Format   : Maple text-file
# Authors  : Jean-Francois Beraud / Bruno Gauthier
# About    : Guess package
# Date     : Wed Jan 27 08:32:03 MET 1999
###############################################################################
##
##    Authors:  Jean-Francois Beraud
##              <Jean-Francois.Beraud@sic.sp2mi.univ-poitiers.fr>
##              Bruno Gauthier
##              <Bruno.Gauthier@univ-mlv.fr>
##              http://www-igm.univ-mlv.fr/~gauthier/
##    Address:  Institut Gaspard Monge,
##              Universite of Marne la Vallee, 
##              5, boulevard Descartes, Champs-sur-Marne,
##              F-77454 MARNE LA VALLEE CEDEX 2, FRANCE
##
## All features of this package are due to C. Krattenthaler.
## 
## Modifications:
##
##   1.00   this version                                          Feb 99
##   0.90   preliminary version                                   Jan 99



########################################
# Initialisation function of the package
########################################
`GUESS/init` := proc()
  interface(verboseproc=0);
  with(linalg,linsolve):
  interface(verboseproc=1);
end:

######################
# Product of two lists
######################
`GUESS/ListProd` := proc(l1,l2)
local i;
  [seq(op(i,l1)*op(i,l2), i=1..nops(l1))]
end:

####################
# Opposite of a list
####################
`GUESS/ListMinus` := proc(l)
  map(x->-x,l)
end:

#################################################
# Rational Interpolation
# Gives the rational interpolant to f (a function 
# of variable x), where m and k are respectively 
# the degrees of the numerator and denominator,
# and xlist is a list of m+k+1 abscissas of the
# interpolation points
#################################################
`GUESS/RationalInterpolation` := proc(f, x, m, k, xlist)
local i, xx, fx, mat, tempvec, mymax;

  xx  := xlist;

  fx  := map(f,xx);
  mat := [seq(1,i=1..m+k+2)];
  tempvec := [seq(1,i=1..m+k+1)];
  mat := subsop(1=tempvec,mat);
  mat := subsop(m+2=`GUESS/ListMinus`(fx),mat);

  mymax := max(m,k);
  for i from 1 to mymax
  do
     tempvec := `GUESS/ListProd`(tempvec,xx);
     if (i <= m)
        then mat := subsop(i+1=tempvec,mat);
     fi:
     if (i <= k)
        then mat := subsop(i+m+2=`GUESS/ListMinus`(`GUESS/ListProd`(tempvec,fx)),mat);
     fi:
  od:

  xx := linsolve( transpose(array(mat)), vector([seq(0,i=1..m+k+1)]) );
  if ( xx = NULL )
     then RETURN(NULL);
  fi:
  xx := convert(xx,list);
  xx := subsop(1=subs(vector(0)=0,op(1,xx)),xx); #XXX???
  normal(
          (op(1,xx)+sum('op(i+1,xx)*x^i','i'=1..m))
          / sum('op(i,xx)*x^(i-m-2)','i'=m+2..m+k+2)
        );
end:

#######################
# Intermediate function
#######################
`GUESS/GuessCons` := proc( l, t )
local i,k,func,res,ri,x;
 
  for i from 1 to nops(l)
  do
     func(i) := op(i,l);
  od:

  res := [];

  for i from 0 to nops(l)-2
  do
     ri := `GUESS/RationalInterpolation`(func,x,nops(l)-i-2,i,[seq(k,k=1..nops(l)-1)]);
     if ( ri <> NULL )
        then
             if ( ( simplify(subs(x=nops(l),denom(ri))) <> 0 ) 
                  and
                  ( simplify(subs(x=nops(l),ri)-func(nops(l))) = 0 ) )
                then res := [subs(x=t,ri),op(res)];
             fi:
     fi:
  od:
convert(convert(res,set),list);
end:

################################################
# Main function of the package
# it tries to find a closed form  for a sequence 
# within the  hierarchy  of  expressions  of the 
# form <rational function>, <product of rational
# functions>,  <product  of product  of rational 
# functions>, etc. It may give several answers
################################################
GUESS[Guess] := proc(L)
local i, k, l, lsize, g, mymin, res, v, flag, unev;

  l     := L;
  lsize := nops(l);
  g     := array(1..lsize-1);
  for i from 1 to lsize-1 do g[i] := 0 od:
  res   := [];
  mymin := lsize-1;
  flag  := false;
  unev  := false;


  if ( member('one', [args]) )
     then flag := true;
  fi:
  if ( member('nogamma', [args]) )
     then unev := true;
  fi:
  if ( nargs>=2 and type(args[2],integer) )
     then mymin := min(args[2]+1,lsize-1);
  fi:
  
  for k from 1 to mymin do
        g[k] := l;
        l := [seq(op(i+1,l)/op(i,l),i=1..lsize-k)];
  od:

  for k from 1 to mymin 
  do
        l := `GUESS/GuessCons`(g[k],_i[k-1]);
        if ( l <> [] )
           then for i from 1 to k-1 
                do
                   if ( unev )
                      then l := [seq(op(1,g[k-i])*
                                 Product(op(v,l),_i[k-i]=1.._i[k-i-1]-1),
                                 v=1..nops(l))];
                      else l := [seq(op(1,g[k-i])*
                                 product(op(v,l),_i[k-i]=1.._i[k-i-1]-1),
                                 v=1..nops(l))];
                   fi:
                od:
                if ( flag )
                   then RETURN(op(l));
                fi:
                res := [ op(l), op(res) ];
        fi:
  od:
convert(res,set);
end:

##############################
# Guess function restricted to
# polynomial interpolation
##############################
GUESS[GuessPoly] :=proc(l)
local i,func,x,ri,lsize;

  lsize := nops(l);
  
  for i from 1 to lsize
  do
     func(i) := op(i,l);
  od:
  ri := `GUESS/RationalInterpolation`(func,x,lsize-2,0,[seq(i,i=1..lsize-1)]);
  if ( ri = NULL )
     then RETURN(NULL);
  fi:
  if ( simplify( subs(x=lsize,ri)-func(lsize) ) = 0 )
     then subs(x=_i[0],ri);
     else NULL;
  fi:
end:

##############################
# Guess function restricted to
# rational interpolation
##############################
GUESS[GuessRatio] := proc(l)
  GUESS[Guess](l,0);
end:

#############
# ONLINE HELP
#############

`help/text/GUESS` := TEXT(
` `,
`HELP FOR: A package to guess closed form for a sequence of numbers.`,
` `,
`CALLING SEQUENCE:`,
`  GUESS[<function>](args)`,
`  <function>(args)`,
` `,
`SYNOPSIS:`,
`- This  package  provides  functions  to find a closed form for a sequence`,
`  of  numbers   within  the  hierarchy  of   expressions   of   the   form`,
`  <rational function>, <product of rational function>, <product of product`,
`  of rational function>, etc.`,
` `,
`- Whenever there is a conflict between a function name in GUESS and another`,
`  name used in the same session, use the long form GUESS[<function>].`,
` `,
`- The available functions are:`,
` `,
`     Guess             GuessPoly         GuessRatio`,
` `,
`- For help with a particular function, type either  ?GUESS[<function>]  or`,
`  ?GUESS,<function>  where <function>  is   one   from   the  above  list.`,
` `,
`EXAMPLES:`,
`> with(GUESS):`,
` `,
`> Guess( [1,2,3] );`,
` `,
`                                    {_i[0]}`,
` `,
`> Guess( [1,2,7,42,429,7436,218348,10850216] );`,
` `,
` _i[0] - 1`,
`  --------'         _i[1]   _i[1] 2                                        1/2`,
` '  |  |       (3/4)      (3     )  GAMMA(_i[1] + 4/3) GAMMA(_i[1] + 2/3) 3`,
`{   |  |   3/4 ---------------------------------------------------------------`,
`    |  |                 _i[1] 2`,
`    |  |               (2     )  GAMMA(_i[1] + 3/2) GAMMA(_i[1] + 1/2)`,
` _i[1] = 1`,
` `,
`    }`,
` `
):


`help/GUESS/text/Guess` := TEXT(
` `,
`FUNCTION: Guess - products of rational interpolation`,
`                - for a sequence of numbers`,
` `,
`CALLING SEQUENCE:`,
`  Guess(l,<options>)`,
`  GUESS[Guess](l,<options>)`,
` `,
`PARAMETERS:`,
`  l               - a list of numbers`,
`  level           - an integer (optional, as to be in second position)`,
`  'one'           - the string 'one' (optional)`,
`  'nogamma'       - the string 'nogamma' (optional)`,
` `,
`SYNOPSIS:`,
`- Guess(l) tries to find a closed form for a sequence within the hierarchy`,
`  of expressions  of  the  form  <rational function>, <product of rational`,
`  function>, <product of product of rational function>, etc.`,
` `,
`- Guess(l,level)  does the same thing  as  Guess(l)  but  it searches only`,
`  within the first 'level' levels`,
` `,
`- Guess(l,'one')  does the same thing as Guess(l) but it returns the first`,
`  solution it finds.`,
` `,
`- Guess(l,'nogamma')   does  the  same  thing as  Guess(l)  but it returns`,
`  expressions without GAMMA functions.`,
` `,
`- Whenever there is a conflict between the function name Guess and another`,
`  name used in the same session, use the long form GUESS['Guess'].`,
` `,
`EXAMPLES:`,
`> with(GUESS):`,
`> l := [seq(product(product(GAMMA(i)*i^2,i=1..j),j=1..k),k=1..10)]:`,
` `,
`> Guess( l, 1 );`,
`                                       {}`,
` `,
`> Guess( l );`,
` `,
`                  _i[0] - 1    /_i[1] - 1                   \\`,
`                   --------'   | --------'                 2|`,
`                  '  |  |      |'  |  |    GAMMA(_i[2] + 3) |`,
`                 {   |  |    4 |   |  |    -----------------|}`,
`                     |  |      |   |  |     GAMMA(2 + _i[2])|`,
`                     |  |      |   |  |                     |`,
`                  _i[1] = 1    \\_i[2] = 1                   /`,
` `,
`> Guess( l , 'nogamma' );`,
` `,
`         _i[0] - 1    /_i[1] - 1     /_i[2] - 1                      \\\\`,
`          --------'   | --------'    | --------'                    2||`,
`         '  |  |      |'  |  |       |'  |  |    9 + 6 _i[3] + _i[3] ||`,
`        {   |  |    4 |   |  |    18 |   |  |    --------------------||}`,
`            |  |      |   |  |       |   |  |          2 + _i[3]     ||`,
`            |  |      |   |  |       |   |  |                        ||`,
`         _i[1] = 1    \\_i[2] = 1     \\_i[3] = 1                      //`,
` `,
`SEE ALSO: GUESS, GuessPoly, GuessRatio,`
):


`help/GUESS/text/GuessPoly` := TEXT(
` `,
`FUNCTION: GuessPoly - polynomial interpolation for a sequence of numbers`,
` `,
`CALLING SEQUENCE:`,
`  GuessPoly(l)`,
`  GUESS[GuessPoly](l)`,
` `,
`PARAMETERS:`,
`  l               - a list of numbers`,
` `,
`SYNOPSIS:`,
`- The return value of GuessPoly(l) is a polynomial in the variable _i[0].`,
` `,
`- Whenever there is a conflict between the function  name  GuessPoly and`,
`  another  name  used   in   the   same  session,   use  the  long  form`,
`  GUESS['GuessPoly'].`,
` `,
`EXAMPLES:`,
`> with(GUESS):`,
`> GuessPoly( [1,4,9,16] );`,
` `,
`                                          2`,
`                                     _i[0]`,
` `,
`> GuessPoly( [1,3,6,10] );`,
` `,
`                             1/2 _i[0] (1 + _i[0])`,
` `,
`SEE ALSO: GUESS, Guess, GuessRatio,`
):


`help/GUESS/text/GuessRatio` := TEXT(
` `,
`FUNCTION: GuessRatio - rational interpolation for a sequence of numbers`,
` `,
`CALLING SEQUENCE:`,
`  GuessRatio(l)`,
`  GUESS[GuessRatio](l)`,
` `,
`PARAMETERS:`,
`  l               - a list of numbers`,
` `,
`SYNOPSIS:`,
`- The  return  value  of  GuessRatio(l)   is  a rational function in the`,
`  variable _i[0].`,
` `,
`- Whenever there is a conflict between the function name  GuessRatio and`,
`  another  name  used   in   the   same  session,   use  the  long  form`,
`  GUESS['GuessRatio'].`,
` `,
`EXAMPLES:`,
`> with(GUESS):`,
`> l := [seq((i^2+1)/i,i=1..10)];`,
` `,
`                                                                  101`,
`          l := [2, 5/2, 10/3, 17/4, 26/5, 37/6, 50/7, 65/8, 82/9, ---]`,
`                                                                   10`,
` `,
`> GuessRatio( l );`,
` `,
`                                            2`,
`                                   1 + _i[0]`,
`                                  {----------}`,
`                                      _i[0]`,
` `,
`SEE ALSO: GUESS, Guess, GuessPoly,`
):


save `GUESS.m`;
quit

