# The Twelve Days of Christmas

```//
// Written by Grant Jenks
// http://www.grantjenks.com/
//
// DISCLAIMER
// THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE
// LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR
// OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND,
// EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE
// ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.
// SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY
// SERVICING, REPAIR OR CORRECTION.
//
// To view a copy of this license, visit
// or send a letter to Creative Commons, 171 Second Street, Suite 300,
// San Francisco, California, 94105, USA.
//
// On the way to church last weekend, my wife began singing the Twelve Days Of
// Christmas. You probably already know it so I'll simply list here the gifts
// given in the song:
//
//    1  Partridge in a Pear Tree
//    2  Turtle Doves
//    3  French Hens
//    4  Colly Birds
//    5  Gold Rings
//    6  Geese-a-Laying
//    7  Swans-a-Swimming
//    8  Maids-a-Milking
//    10 Lords-a-Leaping
//    11 Pipers Piping
//    12 Drummers Drumming
//
// As she sang the song, I wondered how many total gifts were given. In my mind,
// I could represent this as:
//
//            n
//           ---
//       y = \    i
//           /
//           ---
//           i=0
//
// I then tried to calculate the answer for (n = 12). I remembered, too, that there
// for a while but couldn't remember it. So I went about deriving it. Here's what
// I first wrote out:
//
//      y(0) = 0
//               > 1
//      y(1) = 1     > 1
//               > 2
//      y(2) = 3     > 1
//               > 3
//      y(3) = 6     > 1
//               > 4
//      y(4) = 10
//
// This is a mapping of y(n) values with differeneces taken after each column of
// '>'. These differences may be thought of as the derivatives of the polynomial
// which I am trying to derive. That the second order derivative is constant
// indicates the polynomial is of the form y(n) = A n^2 + B n + C. Furthermore,
// since y(0) = 0, I know C = 0. So the rest of the derivation is solving a
// system of equations. I chose to solve over y(1) and y(2) and validate on y(3)
// and y(4):
//
//      y(1) = 1 = A + B
//      y(2) = 3 = 4 A + 2 B
//
//   for which: A = 0.5 and B = 0.5, so:
//
//      y(n) = 0.5 n^2 + 0.5 n
//
//   which is equivalent to:
//
//             (n + 1)(n)
//      y(n) = ----------
//                  2
//
//   and y(12) = 78.
//
// My wife didn't think of the problem this way. She thought that each day
// brought a new gift and a repeat of all the gifts before it. I could
// represent this as:
//
//            n    i
//           ---  ---
//       y = \    \    j
//           /    /
//           ---  ---
//           i=0  j=0
//
// For which, the program below will derive the following solution:
//
//       f(x) = (1/6) x^3 + (1/2) x^2 + (1/3) x
//
// Which may also be written as:
//
//           (2n + 4)(n + 1)(n)
//       y = ------------------
//                  12
//
// Based on this input:
//
//    Equation.exe 0 1 4 10 20 35
//
// Usage:
//
//    Equation.exe ((-)?[1-9][0-9]*)+
//
// Output:
//
//    A polynomial equation fitting the given values.
//

using System;
using System.Linq;
using System.Text;
using System.Collections.Generic;

// For me, this requires a special reference on the build line. I have to thus
// compile this file with:
// csc.exe /r:c:\Windows\Microsoft.NET\Framework\v4.0.30319\System.Numerics.dll Equation.cs

using System.Numerics;

class Equation
{

static int Main (String[] args)
{
if (args.Length == 0)
{
Console.WriteLine("Error: At least one argument required.");
Console.WriteLine("   Usage:");
Console.WriteLine();
Console.WriteLine("      Equation.exe ((-)?[1-9][0-9]*)+");
Console.WriteLine();
Console.WriteLine("   Output:");
Console.WriteLine();
Console.WriteLine("      A polynomial equation fitting the given values.");
return 2;
}

// Assume the input represents a mapping for function, f(x), starting at x = 0.

List<Fraction> input = new List<Fraction>();

try
{
foreach (String arg in args)
{
}
}
catch (Exception e)
{
Console.WriteLine("ERROR: {0}", e.Message);
return 1;
}

int order = DetermineOrder(input);

if (order == -1)
{
Console.WriteLine("ERROR: Not enough values in input.");
return 1;
}

List<Fraction> coefficients = DetermineCoefficients(input, order);

if (coefficients.Count == 0)
{
Console.WriteLine("ERROR: Failed to calculate coefficients.");
return 1;
}

String solution = PrettyPrintSolution(coefficients);

Console.WriteLine(solution);

return 0;
}

//------------------------------------------------------------------------------
//
// Determine the order of the polynomial for the given data.
//
//------------------------------------------------------------------------------

static int DetermineOrder(List<Fraction> input)
{
int order = 0;
bool allZeros = input.Aggregate(true,
(accumulator, next) => accumulator && (next == 0));

if (allZeros)
{
return order;
}

List<Fraction> diffFrom = new List<Fraction>(input);
List<Fraction> diffTo = new List<Fraction>();

while (!allZeros)
{
order++;

for (int i = 0; i < (diffFrom.Count - 1); i++)
{
}

if (diffTo.Count == 0)
{
return -1;
}

allZeros = diffTo.Aggregate(true,
(accumulator, next) => accumulator && (next == 0));
diffFrom.Clear();
diffTo.Clear();
}

return (order - 1);
}

//------------------------------------------------------------------------------
//
// Determine coefficients.
//
//------------------------------------------------------------------------------

static List<Fraction> DetermineCoefficients(List<Fraction> input, int order)
{
List<Fraction> coefficients = new List<Fraction>();

// Since we assume the input sequence starts at (x = 0), we know one of
// the coefficients immediately.

if (order == 0)
{
return coefficients;
}
else if (order == 1)
{
// y(x) = A x + B

}
else
{
coefficients = SolveHigherOrderSystem(input, order);
}

return coefficients;
}

//------------------------------------------------------------------------------
//
// Solve linear system of equations.
//
//------------------------------------------------------------------------------

static List<Fraction> SolveHigherOrderSystem(List<Fraction> input, int order)
{
// Build a matrix and put it in reduced row echelon form.

int rows = order + 1;
int cols = order + 2;
Fraction[,] matrix = new Fraction[rows, cols];

// Setup matrix mapping to system of linear equations.

for (int i = 0; i < rows; i++)
{
for (int j = 0; j < cols; j++)
{
if (j == (cols - 1))
{
matrix[i, j] = input[i + 1];
}
else
{
matrix[i, j] = BigInteger.Pow(i + 1, rows - j - 1);
}
}
}

TransformMatrixToReducedRowEchelonForm(matrix);

// Collect coefficients from transformed matrix.

List<Fraction> coefficients = new List<Fraction>();

for (int i = (rows - 1); i >= 0; i--)
{
}

return coefficients;
}

//------------------------------------------------------------------------------
//
// Transform matrix to row-echelon form using Gauss-Jordan elimination.
// Some implementations I have seen are more general than this and catch the
// case where the matrix cannot be transformed to reduced-row-echelon form. I
// do not catch that case because it cannot happen here. The algorithm used to
// determine order guarantees enough values in the input and enough regularity
// about them for this to work.
//
//------------------------------------------------------------------------------

static void TransformMatrixToReducedRowEchelonForm(Fraction[,] matrix)
{
int rows = matrix.GetLength(0);
int cols = matrix.GetLength(1);
Fraction temp;

// Transform the matrix to row-echelon form.

for (int i = 0; i < rows; i++)
{
temp = matrix[i, i];

// Scale this row down so that it starts with a 1.

for (int j = 0; j < cols; j++)
{
matrix[i, j] = matrix[i, j] * (1 / temp);
}

// Perform row-subtraction so that the remaining values in this column
// below this row are zero.

for (int j = (i + 1); j < rows; j++)
{
temp = matrix[j, i];

for (int k = 0; k < cols; k++)
{
matrix[j, k] = matrix[j, k] - (temp * matrix[i, k]);
}
}
}

// Backsubstitute to get reduced row-echelon form.

for (int i = (rows - 1); i >= 0; i--)
{
for (int j = (i - 1); j >= 0; j--)
{
temp = matrix[j, i];

for (int k = 0; k < cols; k++)
{
matrix[j, k] = matrix[j, k] - (matrix[i, k] * temp);
}
}
}
}

//------------------------------------------------------------------------------
//
// Pretty-printer for solution.
// This could be a lot prettier if it printed factored polynomials.
//
//------------------------------------------------------------------------------

static String PrettyPrintSolution(List<Fraction> coefficients)
{
coefficients.Reverse();
int exponent = coefficients.Count;
StringBuilder equation = new StringBuilder("f(x) = ");

foreach (Fraction c in coefficients)
{
exponent--;

// Skip zero coefficients unless f(x) = 0.

if (c == (new Fraction(0, 1)) && coefficients.Count > 1)
{
continue;
}

// The first term needs no plus sign.

string op = (exponent == (coefficients.Count - 1) ? "" : " + ");

if (exponent > 1)
{
equation.Append(String.Format("{2}{0} x^{1}", c, exponent, op));
}
else if (exponent == 1)
{
equation.Append(String.Format("{1}{0} x", c, op));
}
else
{
equation.Append(String.Format("{1}{0}", c, op));
}
}

// Now some fix-ups.

equation.Replace("+ -", "- ");
equation.Replace(" 1 x", " x");
equation.Replace("-1 ", "- ");
equation.Replace("+ (-", "- (");
equation.Replace("= (-", "= - (");

return equation.ToString();
}

//------------------------------------------------------------------------------
//
// Fraction struct for doing arbitrary rational arithmetic.
// This object is by no means complete and not well tested. I just didn't like
// the idea of using floating-point for the Gauss-Jordan Elimination. Doing so
// would have required an 'epsilon' that would've given me headaches.
//
//------------------------------------------------------------------------------

struct Fraction
{
private BigInteger numerator;
private BigInteger denominator;

{
numerator = aNumerator;
this.Canonicalize();
}

public static implicit operator Fraction(BigInteger aNumerator)
{
return new Fraction(aNumerator, 1);
}

public static implicit operator Fraction(int aNumerator)
{
return new Fraction(aNumerator, 1);
}

//---------------------------------------------------------------------------
//
// Handles only integers.
//
//---------------------------------------------------------------------------

public static Fraction Parse(String value)
{
return new Fraction(BigInteger.Parse(value), 1);
}

//---------------------------------------------------------------------------
//
// Handles only a subset of operators: -, *, /, ==, !=
//
//---------------------------------------------------------------------------

public static Fraction operator -(Fraction f1, Fraction f2)
{
return new Fraction((f1.numerator * f2.denominator - f2.numerator * f1.denominator),
(f1.denominator * f2.denominator));
}

public static Fraction operator *(Fraction f1, Fraction f2)
{
return new Fraction((f1.numerator * f2.numerator), (f1.denominator * f2.denominator));
}

public static Fraction operator /(Fraction f1, Fraction f2)
{
return new Fraction((f1.numerator * f2.denominator), (f1.denominator * f2.numerator));
}

public static bool operator ==(Fraction f1, Fraction f2)
{
return (f1.numerator == f2.numerator && f1.denominator == f2.denominator);
}

public static bool operator !=(Fraction f1, Fraction f2)
{
return !(f1 == f2);
}

public override String ToString()
{
if (numerator == 0)
{
return "0";
}
else if (denominator == 1 || denominator == -1)
{
return (numerator * denominator).ToString();
}
else
{
if (denominator < 0)
{
return String.Format("({0}/{1})", (-1 * numerator), (-1 * denominator));
}
else
{
return String.Format("({0}/{1})", numerator, denominator);
}
}
}

//---------------------------------------------------------------------------
//
// For equality to work, canonicalization is critical.
//
//---------------------------------------------------------------------------

private void Canonicalize()
{
if (denominator == 0)
{
throw new Exception("Divide by zero error.");
}

if (numerator == 0)
{
denominator = 1;
return;
}

if (numerator > 0 && denominator > 0)
{
// Do nothing.
}
else if (numerator > 0 && denominator < 0)
{
// Do nothing. Denominator will carry sign bit.
}
else if (numerator < 0 && denominator > 0
|| numerator < 0 && denominator < 0)
{
// Flip the sign bits.

numerator *= -1;
denominator *= -1;
}

BigInteger gcd = GCD(numerator, denominator);

while (gcd > 1)
{
numerator /= gcd;
denominator /= gcd;
gcd = GCD(numerator, denominator);
}
}

private BigInteger GCD(BigInteger a, BigInteger b)
{
while (b != 0)
{
BigInteger temp = b;
b = a % b;
a = temp;
}
return BigInteger.Abs(a);
}

//---------------------------------------------------------------------------
//
// Implemented only to appease the compiler.
//
//---------------------------------------------------------------------------

public override bool Equals(Object o)
{
if (o is Fraction)
{
return (this == ((Fraction)o));
}
else
{
return false;
}
}

public override int GetHashCode()
{
int hash = 17;

unchecked
{
hash = hash * 31 + numerator.GetHashCode();
hash = hash * 31 + denominator.GetHashCode();
}
return hash;
}
}
}```