use functionUtils.frink

/** This program contains routines to predict a sequence.  These are generally
    achieved by finding the differences or quotients of terms, but there are
    other techniques that, say, try to force a polynomial fit onto the data.
*/
fit[list, numAdditionalTerms=1, debug=false] :=
{
   level = 0
   do
   {
      [result, perfectFit, nextRowPoly, symbolic, polyFunc] = polynomialFit[list, numAdditionalTerms-level, debug, true]
      if perfectFit
         return polyFunc
      else
      {
         println["nextRowPoly is $nextRowPoly"]
         [result, perfectFit, nextRowQuot, symbolic, quotFunc] = quotientFit[nextRowPoly, numAdditionalTerms-level-1, debug, true]
         if perfectFit
            return polyFunc[quotFunc[x]]
      }
      
      [result, perfectFit, nextRowQuot, symbolic, quotFunc] = quotientFit[list, numAdditionalTerms-level, debug, true]
      if perfectFit
         return quotFunc
      else
      {
         println["nextRowQuot is $nextRowQuot"]
         [result, perfectFit, nextRowPoly, symbolic, polyFunc] = polynomialFit[nextRowQuot, numAdditionalTerms-level-1, debug, true]
         if perfectFit
            return quotFunc[polyFunc[x]]
      }

      level = level + 1
   } while level < numAdditionalTerms

   return "match not found"
}

/** This is one version that forces a polynomial fit onto the data.  See
    http://extremelearning.com.au/a-simple-formula-for-sequences-and-series/

    although it's not clear from the article that it's fitting a polynomial
    onto the data and that it won't work for geometric terms, or a Fibonacci
    sequence, or whatever.
*/
polynomialFit[list, numAdditionalTerms=10, debug=false, returnExtra=false] :=
{
   list2 = deepCopy[list]
   first = new array
   first@0 = list2@0
   
   origlen = length[list]
   var nextRow
   var symbolic
   level = 0

   LEVEL:
   while (len = length[list2]) >= 2
   {
      diffs = new array
      allZero = true
      for i=0 to len-2
      {
         diff = list2@(i+1) - list2@i
         diffs@i = diff
         if (diff != 0)
            allZero = false
      }

      if debug
         println[diffs]
      
      first.push[diffs@0]
      list2 = diffs
      if returnExtra and level==0
         nextRow = diffs
      
      if allZero == true or len <= 2
      {
         if debug
            println["first is " + first]
         result = new array
         for n = 0 to origlen + numAdditionalTerms - 1
         {
            Tn = 0
            symbolic = 0
            for i = rangeOf[first]
            {
               Tn = Tn + first@i * binomial[n, i]
               symbolic = symbolic + first@i * binomialSymbolic[n,i]
            }
            result.push[Tn]
            if debug
               println["Symbolic is $symbolic"]
         }
         if returnExtra == false
            return result
         else
            return [result, allZero, nextRow, symbolic, toFunction[noEval[n], symbolic]]
       }
       level = level + 1
   }
}

/** This tries the quotient of terms to try and predict a series.
*/
quotientFit[list, numAdditionalTerms=10, debug=false, returnExtra=false] :=
{
   list2 = deepCopy[list]
   first = new array
   first@0 = list2@0
   
   difftri = new array
   difftri.push[list2]
   
   origlen = length[list]

   level:
   while (len = length[list2]) >= 2
   {
      diffs = new array
      allone = true
      for i=0 to len-2
      {
         diff = list2@(i+1) / list2@i
         diffs@i = diff
         if (diff != 1)
            allone = false
      }

      difftri.push[diffs]
      
      if debug
         println[diffs]
      
      first.push[diffs@0]
      list2 = diffs
      
      if allone == true or len <= 2
      {
         if debug
            println["first is " + first]
         if debug
            println["difftri is\n" + join["\n", difftri]]
         result = deepCopy[list]
         for n = 1 to numAdditionalTerms
         {
            // copy last entry in difftri
            lastRow = difftri@(length[difftri] - 1)
            lastRow.push[lastRow@(length[lastRow]-1)]

//            if debug
//               println["New sorta difftri is $difftri"]
            for i=length[difftri]-2 to 0 step -1
            {
               row     = difftri@i
               nextRow = difftri@(i+1)
               row.push[row@(length[row]-1) * nextRow@(length[nextRow]-1)]
            }

            if debug
               println["New difftri is " + difftri]
         }
        
         if returnExtra == false
            return difftri@0
         else
         {
            symbolic = list@0 * (list@1/list@0) ^ noEval[n]
            return [difftri@0, allone, difftri@1, symbolic, toFunction[noEval[n], symbolic]]
         }
      }
   }
}

binomialSymbolic[n,k] :=
{
   if (k<=0 or k>=n)
      return 1
   
   if (n - k) > k
      k = (n-k)

   product = 1 / (n-k)!
   
   for i = 0 to n-k-1
      product = product * (noEval[n] - i)
   return product 
}

toFunction[symbol, symbolic] :=
{
   return constructExpression["AnonymousFunction", [[makeSymbol[symbol]], symbolic]]
}