p101.cpp

/**
 * If we are presented with the first k terms of a sequence {u[n]} it is
 * impossible to deduce with certainty the value of the next term. Instead,
 * based on the principle of simplicity we fit the known terms with a
 * polynomial of the least degree to predict the next term.
 *
 * Let OP(k, n) denote the n-th term of the optimum (lowest-degree) polynomial
 * generating function for the first k terms of a sequence. It should be clear
 * that OP(k, n) will accurately generate the terms of the sequence for n ≤ k;
 * potentially the first incorrect term (FIT) will be OP(k, k+1).
 *
 * For example, the following OPs are obtained for the cubic sequence
 * u[n] = n^3:
 *
 *   OP(1, n) = 1                = 1, <1>, 1, 1, ...
 *   OP(2, n) = 7*n − 6          = 1, 8, <15>, ...
 *   OP(3, n) = 6*n^2 − 11*n + 6 = 1, 8, 27, <58>, ...
 *   OP(4, n) = n^3              = 1, 8, 27, 64, 125, ...
 *
 * The term in bracket is the first incorrect term for each k. Clearly no FITs
 * exist for k ≥ 4. The sum of FITs for u[n] = n^3 is 1 + 15 + 58 = 74.
 *
 * Consider the sequence generated by the tenth degree polynomial:
 *
 *   u[n] = 1 − n + n^2 − n^3 + n^4 − n^5 + n^6 − n^7 + n^8 − n^9 + n^10
 *
 * Find the sum its FITs.
 */

#include <cstdint>
#include <iostream>
#include <numeric>
#include <type_traits>
#include <vector>
#include "euler.h"

BEGIN_PROBLEM(101, solve_problem_101)
  PROBLEM_TITLE("Optimum polynomial")
  PROBLEM_ANSWER("37076114526")
  PROBLEM_DIFFICULTY(1)
  PROBLEM_FUN_LEVEL(2)
  PROBLEM_TIME_COMPLEXITY("n^3")
  PROBLEM_SPACE_COMPLEXITY("n")
  PROBLEM_KEYWORDS("curve fitting")
END_PROBLEM()

// @timecomplexity: O(n^2)
// @spacecomplexity: O(n)
template <class Iter>
typename std::iterator_traits<Iter>::value_type
next_in_sequence(Iter first, Iter last)
{
  using T = typename std::iterator_traits<Iter>::value_type;
  if (first == last)
  {
    return T(0);
  }

  // Diff the sequence until there is only one term.
  std::vector<T> seq(first, last);
  for (size_t i = seq.size() - 1; i >= 1; --i)
  {
    // in-place first difference.
    for (size_t j = 0; j < i; j++)
    {
      seq[j] = seq[j+1] - seq[j];
    }
  }

  // Computes the next term.
  return std::accumulate(seq.begin(), seq.end(), T(0));
}

static void solve_problem_101()
{
#if 0
  const int ORDER = 3;
  typedef int value_t;
  auto u = [](int n) -> value_t { return n*n*n; };
#else
  const int ORDER = 10;
  typedef int64_t value_t;
  auto u = [](int n) -> value_t {
    value_t y = 0;
    for (int i = 0; i <= 10; i++)
    {
      y = 1 - n * y;
    }
    return y;
  };
#endif

  // Display the true polynomial.
  std::vector<value_t> y(ORDER);
  for (int i = 0; i < ORDER; i++)
  {
    y[i] = u(i + 1);
    if (verbose())
    {
      std::cout << "u(" << (i + 1) << ") = " << y[i] << std::endl;
    }
  }

  // Fit 1 to 10 terms, and find the sum of First-Incorrect-Terms.
  value_t total_incorrect = 0;
  for (int terms = 1; terms <= ORDER; terms++)
  {
    value_t fit = next_in_sequence(y.begin(), y.begin() + terms);
    total_incorrect += fit;
  }
  std::cout << total_incorrect << std::endl;
}