/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.apache.commons.math.analysis.polynomials;

import java.util.ArrayList;

import org.apache.commons.math.fraction.BigFraction;

/**
 * A collection of static methods that operate on or return polynomials.
 * 
 * @version $Revision$ $Date$
 * @since 2.0
 */
public class PolynomialsUtils {

	/** Coefficients for Chebyshev polynomials. */
	private static final ArrayList<BigFraction> CHEBYSHEV_COEFFICIENTS;

	/** Coefficients for Hermite polynomials. */
	private static final ArrayList<BigFraction> HERMITE_COEFFICIENTS;

	/** Coefficients for Laguerre polynomials. */
	private static final ArrayList<BigFraction> LAGUERRE_COEFFICIENTS;

	/** Coefficients for Legendre polynomials. */
	private static final ArrayList<BigFraction> LEGENDRE_COEFFICIENTS;

	static {

		// initialize recurrence for Chebyshev polynomials
		// T0(X) = 1, T1(X) = 0 + 1 * X
		CHEBYSHEV_COEFFICIENTS = new ArrayList<BigFraction>();
		CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);
		CHEBYSHEV_COEFFICIENTS.add(BigFraction.ZERO);
		CHEBYSHEV_COEFFICIENTS.add(BigFraction.ONE);

		// initialize recurrence for Hermite polynomials
		// H0(X) = 1, H1(X) = 0 + 2 * X
		HERMITE_COEFFICIENTS = new ArrayList<BigFraction>();
		HERMITE_COEFFICIENTS.add(BigFraction.ONE);
		HERMITE_COEFFICIENTS.add(BigFraction.ZERO);
		HERMITE_COEFFICIENTS.add(BigFraction.TWO);

		// initialize recurrence for Laguerre polynomials
		// L0(X) = 1, L1(X) = 1 - 1 * X
		LAGUERRE_COEFFICIENTS = new ArrayList<BigFraction>();
		LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
		LAGUERRE_COEFFICIENTS.add(BigFraction.ONE);
		LAGUERRE_COEFFICIENTS.add(BigFraction.MINUS_ONE);

		// initialize recurrence for Legendre polynomials
		// P0(X) = 1, P1(X) = 0 + 1 * X
		LEGENDRE_COEFFICIENTS = new ArrayList<BigFraction>();
		LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);
		LEGENDRE_COEFFICIENTS.add(BigFraction.ZERO);
		LEGENDRE_COEFFICIENTS.add(BigFraction.ONE);

	}

	/**
	 * Private constructor, to prevent instantiation.
	 */
	private PolynomialsUtils() {
	}

	/**
	 * Create a Chebyshev polynomial of the first kind.
	 * <p>
	 * <a href=
	 * "http://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html">Chebyshev
	 * polynomials of the first kind</a> are orthogonal polynomials. They can be
	 * defined by the following recurrence relations:
	 * 
	 * <pre>
	 *  T<sub>0</sub>(X)   = 1
	 *  T<sub>1</sub>(X)   = X
	 *  T<sub>k+1</sub>(X) = 2X T<sub>k</sub>(X) - T<sub>k-1</sub>(X)
	 * </pre>
	 * </p>
	 * 
	 * @param degree
	 *            degree of the polynomial
	 * @return Chebyshev polynomial of specified degree
	 */
	public static PolynomialFunction createChebyshevPolynomial(final int degree) {
		return buildPolynomial(degree, CHEBYSHEV_COEFFICIENTS, new RecurrenceCoefficientsGenerator() {
			private final BigFraction[] coeffs = { BigFraction.ZERO, BigFraction.TWO, BigFraction.ONE };

			/** {@inheritDoc} */
			public BigFraction[] generate(int k) {
				return coeffs;
			}
		});
	}

	/**
	 * Create a Hermite polynomial.
	 * <p>
	 * <a href="http://mathworld.wolfram.com/HermitePolynomial.html">Hermite
	 * polynomials</a> are orthogonal polynomials. They can be defined by the
	 * following recurrence relations:
	 * 
	 * <pre>
	 *  H<sub>0</sub>(X)   = 1
	 *  H<sub>1</sub>(X)   = 2X
	 *  H<sub>k+1</sub>(X) = 2X H<sub>k</sub>(X) - 2k H<sub>k-1</sub>(X)
	 * </pre>
	 * </p>
	 * 
	 * @param degree
	 *            degree of the polynomial
	 * @return Hermite polynomial of specified degree
	 */
	public static PolynomialFunction createHermitePolynomial(final int degree) {
		return buildPolynomial(degree, HERMITE_COEFFICIENTS, new RecurrenceCoefficientsGenerator() {
			/** {@inheritDoc} */
			public BigFraction[] generate(int k) {
				return new BigFraction[] { BigFraction.ZERO, BigFraction.TWO, new BigFraction(2 * k) };
			}
		});
	}

	/**
	 * Create a Laguerre polynomial.
	 * <p>
	 * <a href="http://mathworld.wolfram.com/LaguerrePolynomial.html">Laguerre
	 * polynomials</a> are orthogonal polynomials. They can be defined by the
	 * following recurrence relations:
	 * 
	 * <pre>
	 *        L<sub>0</sub>(X)   = 1
	 *        L<sub>1</sub>(X)   = 1 - X
	 *  (k+1) L<sub>k+1</sub>(X) = (2k + 1 - X) L<sub>k</sub>(X) - k L<sub>k-1</sub>(X)
	 * </pre>
	 * </p>
	 * 
	 * @param degree
	 *            degree of the polynomial
	 * @return Laguerre polynomial of specified degree
	 */
	public static PolynomialFunction createLaguerrePolynomial(final int degree) {
		return buildPolynomial(degree, LAGUERRE_COEFFICIENTS, new RecurrenceCoefficientsGenerator() {
			/** {@inheritDoc} */
			public BigFraction[] generate(int k) {
				final int kP1 = k + 1;
				return new BigFraction[] { new BigFraction(2 * k + 1, kP1), new BigFraction(-1, kP1),
						new BigFraction(k, kP1) };
			}
		});
	}

	/**
	 * Create a Legendre polynomial.
	 * <p>
	 * <a href="http://mathworld.wolfram.com/LegendrePolynomial.html">Legendre
	 * polynomials</a> are orthogonal polynomials. They can be defined by the
	 * following recurrence relations:
	 * 
	 * <pre>
	 *        P<sub>0</sub>(X)   = 1
	 *        P<sub>1</sub>(X)   = X
	 *  (k+1) P<sub>k+1</sub>(X) = (2k+1) X P<sub>k</sub>(X) - k P<sub>k-1</sub>(X)
	 * </pre>
	 * </p>
	 * 
	 * @param degree
	 *            degree of the polynomial
	 * @return Legendre polynomial of specified degree
	 */
	public static PolynomialFunction createLegendrePolynomial(final int degree) {
		return buildPolynomial(degree, LEGENDRE_COEFFICIENTS, new RecurrenceCoefficientsGenerator() {
			/** {@inheritDoc} */
			public BigFraction[] generate(int k) {
				final int kP1 = k + 1;
				return new BigFraction[] { BigFraction.ZERO, new BigFraction(k + kP1, kP1), new BigFraction(k, kP1) };
			}
		});
	}

	/**
	 * Get the coefficients array for a given degree.
	 * 
	 * @param degree
	 *            degree of the polynomial
	 * @param coefficients
	 *            list where the computed coefficients are stored
	 * @param generator
	 *            recurrence coefficients generator
	 * @return coefficients array
	 */
	private static PolynomialFunction buildPolynomial(final int degree, final ArrayList<BigFraction> coefficients,
			final RecurrenceCoefficientsGenerator generator) {

		final int maxDegree = (int) Math.floor(Math.sqrt(2 * coefficients.size())) - 1;
		synchronized (PolynomialsUtils.class) {
			if (degree > maxDegree) {
				computeUpToDegree(degree, maxDegree, generator, coefficients);
			}
		}

		// coefficient for polynomial 0 is l [0]
		// coefficients for polynomial 1 are l [1] ... l [2] (degrees 0 ... 1)
		// coefficients for polynomial 2 are l [3] ... l [5] (degrees 0 ... 2)
		// coefficients for polynomial 3 are l [6] ... l [9] (degrees 0 ... 3)
		// coefficients for polynomial 4 are l[10] ... l[14] (degrees 0 ... 4)
		// coefficients for polynomial 5 are l[15] ... l[20] (degrees 0 ... 5)
		// coefficients for polynomial 6 are l[21] ... l[27] (degrees 0 ... 6)
		// ...
		final int start = degree * (degree + 1) / 2;

		final double[] a = new double[degree + 1];
		for (int i = 0; i <= degree; ++i) {
			a[i] = coefficients.get(start + i).doubleValue();
		}

		// build the polynomial
		return new PolynomialFunction(a);

	}

	/**
	 * Compute polynomial coefficients up to a given degree.
	 * 
	 * @param degree
	 *            maximal degree
	 * @param maxDegree
	 *            current maximal degree
	 * @param generator
	 *            recurrence coefficients generator
	 * @param coefficients
	 *            list where the computed coefficients should be appended
	 */
	private static void computeUpToDegree(final int degree, final int maxDegree,
			final RecurrenceCoefficientsGenerator generator, final ArrayList<BigFraction> coefficients) {

		int startK = (maxDegree - 1) * maxDegree / 2;
		for (int k = maxDegree; k < degree; ++k) {

			// start indices of two previous polynomials Pk(X) and Pk-1(X)
			int startKm1 = startK;
			startK += k;

			// Pk+1(X) = (a[0] + a[1] X) Pk(X) - a[2] Pk-1(X)
			BigFraction[] ai = generator.generate(k);

			BigFraction ck = coefficients.get(startK);
			BigFraction ckm1 = coefficients.get(startKm1);

			// degree 0 coefficient
			coefficients.add(ck.multiply(ai[0]).subtract(ckm1.multiply(ai[2])));

			// degree 1 to degree k-1 coefficients
			for (int i = 1; i < k; ++i) {
				final BigFraction ckPrev = ck;
				ck = coefficients.get(startK + i);
				ckm1 = coefficients.get(startKm1 + i);
				coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])).subtract(ckm1.multiply(ai[2])));
			}

			// degree k coefficient
			final BigFraction ckPrev = ck;
			ck = coefficients.get(startK + k);
			coefficients.add(ck.multiply(ai[0]).add(ckPrev.multiply(ai[1])));

			// degree k+1 coefficient
			coefficients.add(ck.multiply(ai[1]));

		}

	}

	/** Interface for recurrence coefficients generation. */
	private static interface RecurrenceCoefficientsGenerator {
		/**
		 * Generate recurrence coefficients.
		 * 
		 * @param k
		 *            highest degree of the polynomials used in the recurrence
		 * @return an array of three coefficients such that P<sub>k+1</sub>(X) =
		 *         (a[0] + a[1] X) P<sub>k</sub>(X) - a[2] P<sub>k-1</sub>(X)
		 */
		BigFraction[] generate(int k);
	}

}
