# Maple package contiguous2f1.mpl for computing contiguous relations of 2F1 series
# (c) Raimundas Vidunas, University of Amsterdam, 2000-2001
#                        Antwerp University, 2002
# Supported by NWO, project number 613-06-565

f21version:= ` Version 3.25, September 10, 2002`:
# Previous version: 3.20, April 29, 2002
# The last changes: correction of the output of contig2f1subs & f21cleargcd
#

#
# Dr. Raimundas Vidunas
# UIA, Antwerp University
# Universiteitsplein 1,
# 2610 Wilrijk
# Belgium 
# email: vidunas@uia.ua.ac.be
# URL:   http://win-www.uia.ac.be/~vidunas
#

contig2f1help:= proc( )
if args=NULL then  
  printf(`  This is a Maple 7.0 package for computing contiguous\n`);
  printf(`      relations of Gauss hypergeometric 2F1 series\n`);
  printf(`             %s\n\n`, f21version);
  printf(`(c) Raimundas Vidunas, University of Amsterdam, 2000-2001\n`);
  printf(`                       Antwerp University, 2002\n`);
  printf(`Supported by the Dutch NWO, project number 613-06-565\n\n`); 
  printf(`By definition, two 2F1 hypergeometric series are contiguous when\n`);
  printf(`their corresponding upper and lower parameters differ by\n`);
  printf(`integers, see G.E.Andrews, R.Askey, R.Roy, "Special Functions",\n`);
  printf(`Cambridge Univ. Press, Cambridge, 1999; chapter 2.5.\n\n`);
  printf(`For any three contiguous 2F1 hypergeometric functions`);
  print('Hypergeom([a+k[1],b+l[1]],[c+m[1]],z)',
    'Hypergeom([a+k[2],b+l[2]],[c+m[2]],z)',
    'Hypergeom([a+k[3],b+l[3]],[c+m[3]],z)');
  printf(`there is a linear (so-called contiguous) relation`);
  print('A*Hypergeom([a+k[1],b+l[1]],[c+m[1]],z)'+
    'B*Hypergeom([a+k[2],b+l[2]],[c+m[2]],z)'+
    'C*Hypergeom([a+k[3],b+l[3]],[c+m[3]],z)'=0);
  printf(`where the coefficients A, B, C are polynomials in\n`);
  printf(`the parameters a, b, c, z. This package computes such\n`);
  printf(`contiguous relations (or their coefficients) for\n`);
  printf(`given integer shifts k[i], l[i], m[i], i=1,2,3).\n\n`);
  printf(`A relevant article: R.Vidunas, "Contiguous relations of\n`);
  printf(`hypergeometric series", http://arXiv.org/abs/math.CA/0109222\n\n`);
  printf(`Available commands: Contiguous2f1, contiguous2f1, reference2f1,\n`);
  printf(`           Contig2f1express, contig2f1express, contig2f1subs,\n`);
  printf(`           contig2f1P, contig2f1Q.\n`);
  printf(`For more information invoke contig2f1help( command_name ).\n\n`);
  printf(`The current contact address:\n`);
  printf(`        Dr. Raimundas Vidunas\n`);
  printf(`        UIA, Antwerp University\n`);
  printf(`        Universiteitsplein 1\n`);
  printf(`        2610 Wilrijk, Belgium\n`);
  printf(`E-mail: vidunas@uia.ua.ac.be\n`);
  printf(`URL:    http://win-www.uia.ac.be/o/vidunas\n\n`);
elif args=`Contiguous2f1` or args=`contiguous2f1` then
  printf(`Given three triples of integers k[i], l[i], m[i] (i=1,2,3),\n`);
  print('Contiguous2f1([k[1],l[1],m[1]], [k[2],l[2],m[2]], [k[3],l[3],m[3]])');
  printf(`returns a contiguous relation\n`);
  print('A*Hypergeom([a+k[1],b+l[1]],[c+m[1]],z)'+
    'B*Hypergeom([a+k[2],b+l[2]],[c+m[2]],z)'+
    'C*Hypergeom([a+k[3],b+l[3]],[c+m[3]],z)'=0,``);
  printf(`whereas function 'contiguous2f1' with the same arguments computes\n`);
  printf(`the coefficient list [A,B,C] of this contiguous relation.\n`);
  printf(`The reference function Hypergeom([a,b],[c],z) is set by\n`);
  printf(`the commando 'reference2f1([a,b],[c],z)'.\n\n`);
  printf(`Both functions can be called with one or two arguments.\n`);
  printf(`The default last arguments are\n`);
  print('[k[2],l[2],m[2]]=[1,0,0] and [k[3],l[3],m[3]]=[0,0,0]');
  printf(`The arguments of function 'Contiguous2f1' can be hypergeometric\n`);
  printf(`functions contiguous to the reference function, say\n`);
  print('Hypergeom([a+2,b],[c-1],z)' or 'hypergeom([b,a+2],[c-1],z)');
  printf(`instead of [2,0,-1].\n\n`); 
  printf(`In general, the resulting contiguous relation is unique.\n`);
  printf(`However, in some special cases (when z=0 or z=1, and/or a, b,\n`);
  printf(`c-a or c-b are integers), the contiguous relation is not unique.\n`);
  printf(`The returned result may not be the most wanted, for example,\n`);
  printf(`a zero coefficient may appear in an undesired place. Note also\n`);
  printf(`that if c is an integer, some hypergeometric series may be\n`);
  printf(`undefined, and contiguous relations would not make sence. See\n`);
  printf(`contig2f1help(reference2f1) for more info on exceptional cases.\n\n`)
elif args=`reference2f1` then
  printf(`The commando 'reference2f1([a,b], [c], z)' sets\n`);
  print('Hypergeom([a,b],[c],z)' and 'Hypergeom([a+1,b],z)');
  printf(`as reference functions for commands Contiguous2f1, contiguous2f1,\n`);
  printf(`Contig2f1express, contig2f1express, contig2f1P and contig2f1Q.\n\n`);
  printf(`Some reference functions may not express all functions\n`);
  printf(`contiguous to them, contiguous relations with them may not be\n`);
  printf(`unique, and/or some contigous functions may be undefined. In the\n`);
  printf(`following three situations 'reference2f1' produces a warning.\n\n`); 
  printf(`If z=0 or z=1, the contiguous relations are not unique.\n`);
  printf(`'contiguous2f1' finds just one contiguous relation, perhaps with\n`);
  printf(`a zero coefficient. The commando 'contig2f1express' may fail\n`);
  printf(`to express a given contiguous 2F1 in the reference functions.\n\n`);
  printf(`If a, b, c-a or c-b are integer(s), the contiguous relations may\n`);
  printf(`not unique (for example, for terminating series), and contiguous\n`);
  printf(`expression of a contiguous 2F1 in the reference functions may\n`);
  printf(`not exist (for example, when the reference functions are\n`);
  printf(`terminating, and the third one does not terminate).\n\n`);
  printf(`When c is an integer, some contiguous hypergeometric functions\n`);
  printf(`(with non-positive lower parameter) are not well-defined.\n`);
  printf(`Then the returned contiguous relation may be meaningless.\n`);
  printf(`If for those functions an upper parameter is a non-positive as\n`);
  printf(`well, and can be interpreted as terminating series, the returned\n`);
  printf(`contiguous relation relation is (accordingly) good.\n\n`);
  printf(`Contiguous relations may not be unique in other cases also,\n`);
  printf(`for example, when z=-1 and c-a+b is an integer (recall Kummer's\n`);
  printf(`formula). However, most contiguous relations are unique in these\n`);
  printf(`other cases, including the one implied by 'contig2f1express'.\n\n`)
elif args=`Contig2f1express` or args=`contig2f1express` or
     args=`contig2f1P` or args=`contig2f1Q` then
  printf(`The commandos 'contig2f1P(k,l,m)' and 'contig2f1Q(k,l,m)' return\n`);
  printf(`the rational coefficients P and Q in a contiguous relation\n`);
  print('Hypergeom([a+k,b+l],[c+m],z)'=
    'P*Hypergeom([a,b],[c],z)'+'Q*Hypergeom([a+1,b],[c],z)');
  printf(`\nThe reference function Hypergeom([a,b],[c],z) is set by\n`);
  printf(`the commando 'reference2f1([a,b],[c],z)'. The returned relation\n`);
  printf(`is the specialization of the corresponding contiguous relation\n`);
  printf(`general a, b, c, z. Check 'contig2f1help(reference2f1)' for\n`);
  printf(`the information on exceptional situations.\n\n`);
  printf(`The commando 'Contig2f1express(k,l,m)' returns the right-hand\n`);
  printf(`of the above expression, and 'contig2f1express(k,l,m)' returns\n`);
  printf(`the coefficient list [contig2f1P(k,l,m), contig2f1Q(k,l,m)].\n`);
  printf(`Function 'Contig2f1express' can be called with one argument,\n`);
  printf(`which must be then a hypergeometric function contiguous to the\n`);
  printf(`reference functions; for example,`);
  print('Contig2f1express(hypergeom([a+k,b+l],[c+m],z))'); print(``)
elif args=`contig2f1subs` then
  printf(`If L is a list of three polynomials, 'contig2f1subs(L)' returns\n`);
  printf(`the same list, but each polynomial is divided by the common GCD.\n`);
  printf(`\nIf the commando 'contig2f1subs' has more arguments, the last\n`);
  printf(`argument is expected to be a list of three polynomials, and\n`);
  printf(`other arguments are interpreted as substitutions. Simplification\n`);
  printf(`of GCD's is performed after each (consequent) substitution.\n`);
  printf(`The function is useful in manipulating output of contiguous2f1.\n\n`)
fi
end:

Contiguous2f1:= proc( )
local F, L, i;
  if nargs=1 then  F:= [args, [1,0,0], [0,0,0]]
    elif nargs=2 then  F:= [args, [0,0,0]]
    elif nargs=3 then  F:= [args]
    else  error("wrong number of arguments, not one, two or three")
  fi;
  L:= F;
  for i to 3 do
    if type(F[i], [integer$3]) then
      F[i]:= subs( f21subs,
        Hypergeom([f21a+F[i][1],f21b+F[i][2]],[f21c+F[i][3]],f21z) )
    else  L[i]:= f21hyper2shift(F[i])  fi
  od;
  L:= contiguous2f1( op(L) );
  L[1]*F[1]+L[2]*F[2]+L[3]*F[3]=0
end:

contiguous2f1:= proc( )
global f21a, f21b, f21c;
local T;
  if nargs=0 or nargs>3
    then ERROR(`wrong number of arguments, not one, two or three`)
  elif not
    convert(map(type,[args],[integer,integer,integer]),`and`)
    then ERROR(`wrong type of arguments, not list(s) of three integers`)
  elif nargs=1 then  f21simplify( f21contig1(op(args[1])) )
  elif nargs=2 then  f21simplify( f21contig2(args[1],args[2]) )
  else  T:= args[3];
    f21simplify( subs(f21a=f21a+T[1], f21b=f21b+T[2], f21c=f21c+T[3],
      f21contig2( args[1]-T, args[2]-T ) ) )
  fi
end:

reference2f1:= proc( a, b, c, z )
global f21a, f21b, f21c, f21z, f21subs;
  f21subs:= [f21a=normal(a), f21b=normal(b), f21c=normal(c), f21z=normal(z)];
  if subs(f21subs,z)=0 then
    WARNING("z=0, hypergeometric functions are trivial, constant 1")
  elif subs(f21subs,z)=1 then
    WARNING("z=1, 2F1's satisfy two term recurrences, recall Gauss formula")
  fi;
  if type(subs(f21subs,a),integer) or type(subs(f21subs,b),integer) then
    WARNING("Terminating 2F1's do not express nonterminating series");
    if type(subs(f21subs,a),negint) or type(subs(f21subs,b),negint) then
      WARNING("The reference functions are terminating")
    fi
  elif type(normal(c-a),integer) or type(normal(c-b),integer) then
    WARNING("Terminating (up to a power of 1-z) 2F1's do not express other series");
    if type(normal(c-a),nonposint) or type(normal(c-b),nonposint) then
      WARNING("The reference functions are terminating up to a power of 1-z")
    fi
  fi;  
  if type(subs(f21subs,c),integer) then
    WARNING("2F1's with lower parameter a negative integer may not be defined");
    if type(subs(f21subs,c),nonposint) then
      WARNING("The reference functions are not properly defined")
    fi
  fi;
  printf(`The current reference functions are the following:`);
  print( Hypergeom([a,b],[c],z), Hypergeom([a+1,b],[c],z) )
end:

Contig2f1express:= proc( )
local L;
  if nargs=1 then  L:= op( f21hyper2shift(args) )
  elif type([args], [integer$3]) then  L:= args
  else  error(`wrong number and/or type or arguments`) 
  fi;
  contig2f1P(L)*subs(f21subs,Hypergeom([f21a,f21b],[f21c],f21z))
  +contig2f1Q(L)*subs(f21subs,Hypergeom([f21a+1,f21b],[f21c],f21z))
end:

contig2f1express:= proc( alpha, beta, gamma )
  [ contig2f1P(alpha,beta,gamma), contig2f1Q(alpha,beta,gamma) ]
end:

contig2f1P:= proc( alpha, beta, gamma )
global f21a, f21b, f21c, f21z;
  try f21simplifyPQ( subs(f21a=f21a+1, f21contsol(alpha-1,beta,gamma))
    *(f21c-f21a-1)/(1-f21z)*simplify( pochhammer(f21c,gamma)
    /pochhammer(f21a+1,alpha-1)/pochhammer(f21b,beta) ) )
  catch "numeric exception: division by zero":
    error("Contiguity coefficients for %1 are not defined with current reference functions",
      subs(f21subs,Hypergeom([f21a+alpha,f21b+beta],[f21c+gamma],f21z)) )
  end try   
end:

contig2f1Q:= proc( alpha, beta, gamma )
global f21a, f21b, f21c, f21z;
  try f21simplifyPQ( f21contsol(alpha,beta,gamma)*simplify(
    pochhammer(f21c,gamma)/pochhammer(f21a+1,alpha-1)/pochhammer(f21b,beta) ))  
  catch "numeric exception: division by zero":
    error("Contiguity coefficients for %1 are not defined with current reference functions",
      subs(f21subs,Hypergeom([f21a+alpha,f21b+beta],[f21c+gamma],f21z)) )
  end try   
end:

contig2f1subs:= proc( )
# Input:  (zero or more) substitutions and a list [F,G,H] of three polynomials
# Output: proportional list of three coprime polynomials, after substitutions
local L, LL, T;
  if args=NULL then error("No arguments") fi;
  LL:= [ seq(args[i], i=1..nargs-1) ];
  L:= f21cleargcd( args[nargs] );
  for T in LL do
    L:= f21cleargcd( subs(T, L) )
  od:
  L
end:

print(`This is a Maple package for computing contiguous relations of`);
print(`Gauss hypergeometric 2F1 series, written by R. Vidunas`);
print(f21version);
print(``);
reference2f1(a,b,c,z);
printf(` You can change them with commando 'reference2f1(a,b,c,z)'\n`);
printf(` The help function is envoked by 'contig2f1help()'\n\n`); 

#
#                  AUXILIARY ROUTINES
#
# Reserved names: f21a, f21b, f21c, f21z, f21subs,
#    f21contsol, f21ra, f21rb, f21rc, f21da, f21db, f21dc,
#    f21dap, f21dcn, f21contig1, f21contig2, f21simplify,
#    f21simplifyPQ, f21cleargcd, f21hyper2shift 
#

f21contsol:= proc( alpha, beta, gamma )
# Computes Q(alpha, beta, gamma)
# Contiguous relations are for 2F1([a,b],[c],z)/Gamma(c)  
global f21a, f21b, f21c, f21z;
local ui, u0, u1, S;
  ui:= f21da( alpha+1 );                   #  alpha+1, 0, 0  #
  u0:= f21da( alpha );                     #   alpha, 0, 0   #
  u1:= normal( ui+(f21b-f21a-alpha)*u0 );  #   alpha, 1, 0   #
  S:= f21a=f21a+alpha;
  ui:= u0*subs(S, f21b=f21b+1, f21rb[3]/f21rb[1]*f21db(beta) )
      +u1*subs(S, f21db(beta+1) );         # alpha, beta+1, 0 #
  u0:= u0*subs(S, f21b=f21b+1, f21rb[3]/f21rb[1]*f21db(beta-1) )
      +u1*subs(S, f21db(beta) );           #  alpha, beta, 0  #
  u1:= normal( ui+(f21c-f21b-beta-1)*u0 ); #  alpha, beta, -1 #
  S:= {S, f21b=f21b+beta};
  normal( u0*subs(S, f21c=f21c-1, f21rc[3]/f21rc[1]*f21dc(gamma+1) )
      +u1*subs(S, f21dc(gamma) ) )         # alpha,beta,gamma #
#print(ui,u0,u1);
end:

### Basic contiguous relations for Gamma(a,b)/Gamma(c)*2F1([a,b],[c],z)
f21ra:= [ 1-f21z, 2*f21a-f21c-f21a*f21z+f21b*f21z,
  (f21c-f21a)*(f21a-1) ]:                     # -FF(a+1),FF,FF(a-1)
f21rb:= subs({f21a=f21b,f21b=f21a}, f21ra):   # -FF(b+1),FF,FF(b-1)
f21rc:= [ 1-f21z, f21c-1-(2*f21c-f21a-f21b-1)*f21z,
  (f21c-f21a)*(f21c-f21b)*f21z ]:             # -FF(c-1),FF,FF(c+1)

f21da:= proc(n::integer)    # a->a+n
global f21a, f21b, f21c, f21z;
local L,j;
if n<0 then
  -subs(f21a=f21a+n, f21dap(-n))/(1-f21z)^n/
    product('(j-f21a)*(f21c-f21a+j-1)','j'=1..-n)
else  f21dap(n)  fi
end:

f21db:= proc(n::integer)    # b->b+n
  subs({f21a=f21b,f21b=f21a}, f21da(n))
end:

f21dc:= proc(n::integer)    # c->c+n
global f21a, f21b, f21c, f21z;
local L;
if n>0 then
  -subs(f21c=f21c+n, f21dcn(n))*(f21z-1)^n/f21z^n/
     product('(j+f21c-f21a)*(j+f21c-f21b)','j'=0..n-1)
else  f21dcn(-n)  fi
end:

f21dap:= proc(n::nonnegint) # a->a+n, n>0
option remember;
global f21a, f21b, f21c, f21z;
local L;
  L:= subs(f21a=f21a+n-1, f21ra);
  normal( (L[2]*f21dap(n-1)+L[3]*f21dap(n-2))/L[1] )
end:
f21dap(0):= 0:
f21dap(1):= 1:

f21dcn:= proc(n::nonnegint) # c->c-n, n>0
option remember;
global f21a, f21b, f21c, f21z;
local L;
  L:= subs(f21c=f21c-n+1, f21rc);
  normal( (L[2]*f21dcn(n-1)+L[3]*f21dcn(n-2))/L[1] )
end:
f21dcn(0):= 0:
f21dcn(1):= 1:

f21contig2:= proc( K, L )
# Input: two triples of integers
# Output: coefficients (to L and reference function) of
#    a contiguous expression to the function given by K
if L=[0,0,0] then  [0,0]   # a very special case
else
  normal( [ f21contsol(op(K))/f21contsol(op(L))*simplify(
      pochhammer(f21a+K[1],L[1]-K[1])*pochhammer(f21b+K[2],L[2]-K[2])
      /pochhammer(f21c+K[3],L[3]-K[3]) ),
    subs(f21a=f21a+K[1], f21b=f21b+K[2], f21c=f21c+K[3],
      f21contsol(op(L-K)) )/f21contsol(op(L))*simplify(
      pochhammer(f21c,K[3])/pochhammer(f21c-f21a,K[3]-K[1])
      /pochhammer(f21c-f21b,K[3]-K[2]) )
      *f21z^(-K[3])*(f21z-1)^(K[3]-K[1]-K[2]) ] )
fi
end:

f21contig1:= proc( alpha, beta, gamma )
# Input:  three integers
# Output: coefficients (to the reference functions) of
#   a contiguous expression of the function given by the arguments
  normal( [ f21contsol(alpha,beta,gamma)*simplify(
    pochhammer(f21c,gamma)/pochhammer(f21a+1,alpha-1)/pochhammer(f21b,beta) ),
  subs(f21a=f21a+1, f21contsol(alpha-1,beta,gamma))*(f21c-f21a-1)/(1-f21z)
    *simplify(pochhammer(f21c,gamma)/pochhammer(f21a+1,alpha-1)
    /pochhammer(f21b,beta) )  ] ) 
end:

f21simplify:= proc( M )
# Input:  list [F,G] of two rational functions in f21a, f21b, f21c, f21z
# Output: list of three coprime polynomials, proportional to [-1,F,G],
#         with substituted f21subs 
local L, T;
  if M=[0,0] then  return [0,-1,1]  fi;  # a very special case
  T:= lcm( denom(M[1]), denom(M[2]) );
  L:= [-T, normal(M[1]*T), normal(M[2]*T)];
  for T in f21subs do
    L:= f21cleargcd( subs(T, L) )
  od
end:

f21simplifyPQ:= proc( F )
# Input:  a rational function
# Output: simplified function after substitution f21subs
global f21a;
local FF, T;
  FF:= normal(F);  
  for T in f21subs do
    FF:= normal( subs(T,FF) ) 
  od 
end:

f21cleargcd:= proc( L )
# Input:  a list of three polynomials
# Output: a proportional list of coprime polynomials
local g;
  g:= gcd(gcd(L[1],L[2]),L[3]);
  if g=0 then  normal(L)
  else  normal( map(`/`, L, g ) )
  fi
end:

f21hyper2shift:= proc( F )
local L;
  if (op(0,F)=`Hypergeom` or op(0,F)=`hypergeom`) and nops(F)=3 then
    if op(3,F)<>subs(f21subs,f21z) then
      error("The argument of the hypergeometric function is\
        different from the argument of the reference function", F)
    fi;
  L:= normal( subs(f21subs,
         [op(1,F)[1]-f21a, op(1,F)[2]-f21b, op(2,F)[1]-f21c]) );
  if not type(L, list(integer)) then
    L:= normal( subs(f21subs,
          [op(1,F)[2]-f21a, op(1,F)[1]-f21b, op(2,F)[1]-f21c]) );
    if not type(L, list(integer)) then
      error("A hypergeometric function is not contiguous to\
            a reference function", F)
    fi
  fi        
  else
    error("An argument does not refer to a hypergeometric function", F)
  fi;
  L
end: