//
// 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.
//
// Copyright
// This work is licensed under the
// Creative Commons Attribution-NonCommercial-NoDerivs 3.0 Unported License.
// To view a copy of this license, visit
// http://creativecommons.org/licenses/by-nc-nd/3.0/
// 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
// 9 Ladies Dancing
// 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
// was a way to reduce that summation to a single equation. I thought about this
// 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)
{
input.Add(Fraction.Parse(arg));
}
}
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++)
{
diffTo.Add(diffFrom[i + 1] - diffFrom[i]);
}
if (diffTo.Count == 0)
{
return -1;
}
allZeros = diffTo.Aggregate(true,
(accumulator, next) => accumulator && (next == 0));
diffFrom.Clear();
diffFrom.AddRange(diffTo);
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.
coefficients.Add(input[0]);
if (order == 0)
{
return coefficients;
}
else if (order == 1)
{
// y(x) = A x + B
coefficients.Add(input[1] - input[0]);
}
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--)
{
coefficients.Add(matrix[i, rows]);
}
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;
public Fraction(BigInteger aNumerator, BigInteger aDenominator)
{
numerator = aNumerator;
denominator = aDenominator;
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;
}
}
}