/* laurent.gp -- Laurent recurrences */
/* 05 Mar 2002 Michael Somos. based on the file with URL
   "http://www.math.harvard.edu/~propp/reach/laurent.html"
   03 Jan 2003 Michael Somos fixed trivial typo in f11.
   03 Jan 2003 Michael Somos modified pm() for new behavior of Str().
 */

undef='undef;
pm(n)=if(n<0,Str("m"(-n)),Str(n))

/* define indexed indeterminate variables */
x1(n)=eval(Str("x_"pm(n)))
x2(n,k)=eval(Str("x_"pm(n)"_"pm(k)))
y2(n,k)=eval(Str("y_"pm(n)"_"pm(k)))
x3(n,i,j)=eval(Str("x_"pm(n)"_"pm(i)"_"pm(j)))

/* ====== The Diamond or Antifrieze recurrence */
/* if x2(n,k)=1 then f1(n,k)=fibonacci(2*n+1) unless n<-1 */

f1(n,k)=
{
  if(n<-1, undef, if(n<1, x2(n,k),
    ( 1 +f1(n-1,k-1)*f1(n-1,k+1) )/f1(n-2,k) ))
} /* end f1() */

/* ====== The Number Wall recurrence */
/* if x2(n,k)=1 then f2(n,k)=2^(n*(n+1)/2) unless n<-1 */

f2(n,k)=
{
  if(n<-1, undef, if(n<1, x2(n,k),
    ( f2(n-1,k)^2 +f2(n-1,k-1)*f2(n-1,k+1) )/f2(n-2,k) ))
} /* end f2() */

/* ====== Robbins-Rumsey condensation */
/* if x3(n,i,j)=1 then f3(n,i,j)=2^(n*(n+1)/2) unless n<-1 */

f3(n,i,j)=
{
  if(n<-1, undef, if(n<1, x3(n,i,j),
    ( f3(n-1,i-1,j)*f3(n-1,i+1,j) +f3(n-1,i,j-1)*f3(n-1,i,j+1) )/f3(n-2,i,j) ))
} /* end f3() */

/* ====== Kuo condensation */
/* "http://www.math.wisc.edu/~propp/somos/kuo3" */
/* "http://www.cs.berkeley.edu/~ekuo/condensation.ps" */

f4(n,i,j)=
{
  if(n<-1, undef, if(n<1, 1,
    ( y2(i+2*n-1,j+2*n-1)*y2(i-2*n+1,j-2*n+1)*f4(n-1,i-2,j+2)*f4(n-1,i+2,j-2)
     +y2(i-2*n+1,j+2*n-1)*y2(i+2*n-1,j-2*n+1)*f4(n-1,i+2,j+2)*f4(n-1,i-2,j-2))
     /f4(n-2,i,j) ))
} /* end f4() */

/* ====== Tilted Diamond Patterns */

f5(n,k)=
{
  if(3*n+k<-3, undef, if(3*n+k<3, x2(n,k),
    ( 1 +f5(n-1,k-1)*f5(n-1,k+1) )/f5(n-2,k) ))
} /* end f5() */

/* ====== Somos-4 recurrence */
/* if x1(n)=1 then f6(n)=Somos_4(n) unless n<0 */

f6(n)=
{
  if(n<0, undef, if(n<4, x1(n),
    ( f6(n-2)^2 +f6(n-1)*f6(n-3) )/f6(n-4) ))
} /* end f6() */

/* ====== S-Arrays */
/* if x2(n,k)=1 then f7(n,k)=Somos_4(n) unless n<0 */

f7(n,k)=
{
  if(n<0, undef, if(n<4, x2(n,k),
    ( f7(n-2,k)^2 +f7(n-1,k+1)*f7(n-3,k-1) )/f7(n-4,k) ))
} /* end f7() */

/* ====== 3D 'Somos-3' face-var */

f8(n,i,j)=
{
  if(n<0, undef, if(n<3, x3(n,i,j),
    ( f8(n-1,i-2,j+2)*f8(n-2,i+2,j-2) +f8(n-1,i+2,j+2)*f8(n-2,i-2,j-2) )
    /f8(n-3,i,j) ))
} /* end f8() */

/* ====== 3D 'Somos-3' edge-var */

f9(n,i,j)=
{
  if(6*n+j<0, undef, if(6*n+j<16, 1,
    ( y2(i+2*n-1,j+2*n-1)*y2(i-2*n+1,j-2*n+1)*f9(n-1,i-2,j+2)*f9(n-1,i+2,j-2)
    + y2(i-2*n+1,j+2*n-1)*y2(i+2*n-1,j-2*n+1)*f9(n-1,i+2,j+2)*f9(n-1,i-2,j-2))
    /f9(n-2,i,j) ))
} /* end f9() */

/* ====== 3D 'Somos-4' face-var */

f10(n,i,j)=
{
  if(n<0, undef, if(n<4, x3(n,i,j),  
    ( f10(n-1,i-2,j+2)*f10(n-3,i+2,j-2) + f10(n-2,i+2,j+2)*f10(n-2,i-2,j-2))
    /f10(n-4,i,j) ))
} /* end f10() */

/* ======= 3D 'Somos-4' edge-var */

f11(n,i,j)=
{
  if(8*n-i+j<0, undef, if(8*n-i+j<16, 1,
    ( y2(-(i+2*n-1),j+2*n-1)*y2(-(i-2*n+1),j-2*n+1)*f11(n-1,i-2,j+2)*f11(n-1,i+2,j-2)
    + y2(-(i-2*n+1),j+2*n-1)*y2(-(i+2*n-1),j-2*n+1)*f11(n-1,i+2,j+2)*f11(n-1,i-2,j-2))
    /f11(n-2,i,j) ))
} /* end f11() */

/* ======= Propp Cube recurrence */
/* if x3(i,j,k)=1 then f12(i,j,k)=3^((i+j+k-1)^2\4) unless i+j+k<0 */

f12(i,j,k)=
{
  if(i+j+k<0, undef, if(i+j+k<3, x3(i,j,k),
    ( f12(i-1,j,k)*f12(i,j-1,k-1)
    + f12(i,j-1,k)*f12(i-1,j,k-1)
    + f12(i,j,k-1)*f12(i-1,j-1,k) )
    /f12(i-1,j-1,k-1) ))
} /* end f12() */