# -*- maple -*-
#
# $Id: guess.maple,v 1.6 2004/05/28 13:12:25 wiz Exp $
#
# File     : guess.maple
# Format   : Maple text-file
# Authors  : Jean-Francois Beraud / Bruno Gauthier / Thomas Klausner
# About    : Guess package
# Date     : $Date: 2004/05/28 13:12:25 $
###############################################################################
##
##    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
##
##    Updated for Maple 8 by:
##              Thomas Klausner
##              <klausner@dmg.tuwien.ac.at>
##              Vienna University of Technology
##              Institute for Discrete Mathematics and Geometry
##              Wiedner Hauptstraße 8-10/104
##              1040 Wien
##              Austria
##
## All features of this package are due to C. Krattenthaler.
## 
## Modifications:
##
##   1.6    updated for Maple 8                                   May 04
##   1.00   this version                                          Feb 99
##   0.90   preliminary version                                   Jan 99


GUESS := module()
    export Guess, GuessPoly, GuessRatio;
    local GuessCons, ListMinus, ListProd, RationalInterpolation;
    option package;

    ######################
    # Product of two lists
    ######################
    ListProd := proc(l1, l2)
	local i;

	[seq(op(i,l1)*op(i,l2), i=1..nops(l1))]
    end proc;

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

    #################################################
    # 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
    #################################################
    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=ListMinus(fx), mat);

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

	xx := linalg[linsolve](transpose(array(mat)), vector([seq(0, i=1..m+k+1)]));
	if (xx = NULL) then
	    NULL
	else
	    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 if
    end proc;

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

	res := [];

	for i from 0 to nops(l)-2 do
	    ri := 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)]
		end if
	    end if
	end do;
	convert(convert(res, set), list)
    end proc;

    ################################################
    # 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 := 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
	end do;
	res := [];
	mymin := lsize-1;
	flag := false;
	unev := false;

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

	for k from 1 to mymin do
	    l := 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))]
		    end if
                end do;
                if (flag) then
                   return op(l)
		end if;

		res := [ op(l), op(res) ]
	    end if
	end do;

	convert(res,set)
    end proc;

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

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

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

savelib('GUESS');

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

`help/GUESS/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`
):

INTERFACE_HELP(insert, topic="GUESS", library=libname[1],
               text=`help/GUESS/text/GUESS`);

INTERFACE_HELP(insert, topic="GUESS/Guess", library=libname[1],
               text=`help/GUESS/text/Guess`);

INTERFACE_HELP(insert, topic="GUESS/GuessPoly", library=libname[1],
               text=`help/GUESS/text/GuessPoly`);

INTERFACE_HELP(insert, topic="GUESS/GuessRatio", library=libname[1],
               text=`help/GUESS/text/GuessRatio`);
	       
quit;

