package org.matheclipse.core.reflection.system;

import static org.matheclipse.core.expression.F.C0;
import static org.matheclipse.core.expression.F.List;
import static org.matheclipse.core.expression.F.evalExpandAll;

import java.util.List;

import org.matheclipse.core.basic.Config;
import org.matheclipse.core.convert.ExprVariables;
import org.matheclipse.core.convert.JASIExpr;
import org.matheclipse.core.eval.exception.JASConversionException;
import org.matheclipse.core.eval.interfaces.AbstractFunctionEvaluator;
import org.matheclipse.core.expression.ASTRange;
import org.matheclipse.core.expression.ExprRingFactory;
import org.matheclipse.core.expression.F;
import org.matheclipse.core.interfaces.IAST;
import org.matheclipse.core.interfaces.IExpr;
import org.matheclipse.core.polynomials.QuarticSolver;

import edu.jas.poly.GenPolynomial;
import edu.jas.poly.Monomial;

/**
 * Determine the roots of a univariate polynomial
 * 
 * See Wikipedia entries for: Quadratic equation , Cubic function and Quartic function
 */
public class Roots extends AbstractFunctionEvaluator {

	public Roots() {
	}

	@Override
	public IExpr evaluate(final IAST ast) {
		if (ast.size() != 2) {
			return null;
		}
		return roots(ast, false);
	}

	protected static IAST roots(final IAST ast, boolean numericSolutions) {
		ExprVariables eVar = new ExprVariables(ast.arg1());
		if (!eVar.isSize(1)) {
			// factor only possible for univariate polynomials
			return null;
		}
		IExpr expr = evalExpandAll(ast.arg1());
		IAST variables = eVar.getVarList();
		IExpr denom = F.C1;
		if (expr.isAST()) {
			expr = Together.together((IAST) expr);

			// split expr into numerator and denominator
			denom = F.eval(F.Denominator(expr));
			if (!denom.isOne()) {
				// search roots for the numerator expression
				expr = F.eval(F.Numerator(expr));
			}
		}
		return rootsOfVariable(expr, denom, variables, numericSolutions);
	}

	protected static IAST rootsOfVariable(final IExpr expr, final IExpr denom, final IAST variables, boolean numericSolutions) {

		if (numericSolutions) {
			IAST result = List();
			IAST resultList = RootIntervals.croots(expr, true);
			if (resultList == null) {
				return null;
			}
			if (resultList.size() > 0) {
				result.addAll(resultList);
			}
			return result;
		}

		IAST result = null;
		ASTRange r = new ASTRange(variables, 1);
		List varList = r.toList();
		try {
			result = F.List();
			IAST factors = Factor.factorComplex(expr, varList, F.List, true);
			for (int i = 1; i < factors.size(); i++) {
				IAST quarticResultList = QuarticSolver.solve(F.evalExpand(factors.get(i)), variables.arg1());
				if (quarticResultList != null) {
					for (int j = 1; j < quarticResultList.size(); j++) {
						result.add(quarticResultList.get(j));
					}
				}
			}
			result = QuarticSolver.createSet(result);
			return result;
		} catch (JASConversionException e) {
			try {
				// try to generate a common expression polynomial
				JASIExpr eJas = new JASIExpr(varList, new ExprRingFactory());
				GenPolynomial ePoly = eJas.expr2IExprJAS(expr);
				result = rootsOfPolynomial(ePoly);
			} catch (JASConversionException e2) {
				if (Config.SHOW_STACKTRACE) {
					e2.printStackTrace();
				}
			}
		}
		if (result != null) {
			if (!denom.isNumber()) {
				// eliminate roots from the result list, which occur in the denominator
				int i = 1;
				while (i < result.size()) {
					IExpr temp = denom.replaceAll(F.Rule(variables.arg1(), result.get(i)));
					if (temp != null && F.eval(temp).isZero()) {
						result.remove(i);
						continue;
					}
					i++;
				}
			}
			return result;
		}
		return null;
	}

	private static IAST rootsOfPolynomial(GenPolynomial ePoly) {
		long varDegree = ePoly.degree(0);
		IAST result = List();
		if (ePoly.isConstant()) {
			return result;
		}
		IExpr a;
		IExpr b;
		IExpr c;
		IExpr d;
		IExpr e;
		if (varDegree <= 4) {
			// solve quartic equation: a*x^2 + b*x + c = 0
			a = C0;
			b = C0;
			c = C0;
			d = C0;
			e = C0;
			for (Monomial monomial : ePoly) {
				IExpr coeff = monomial.coefficient();
				long lExp = monomial.exponent().getVal(0);
				if (lExp == 4) {
					a = coeff;
				} else if (lExp == 3) {
					b = coeff;
				} else if (lExp == 2) {
					c = coeff;
				} else if (lExp == 1) {
					d = coeff;
				} else if (lExp == 0) {
					e = coeff;
				} else {
					throw new ArithmeticException("Roots::Unexpected exponent value: " + lExp);
				}
			}
			result = QuarticSolver.quarticSolve(a, b, c, d, e);
			if (result != null) {
				result = QuarticSolver.createSet(result);
				return result;
			}

		}

		return null;
	}

}
 

Ads help maintain this website.