mirror of
https://github.com/UpsilonNumworks/Upsilon.git
synced 2026-01-19 00:37:25 +01:00
1420 lines
60 KiB
C++
1420 lines
60 KiB
C++
#include <poincare/power.h>
|
|
#include <poincare/addition.h>
|
|
#include <poincare/arithmetic.h>
|
|
#include <poincare/binomial_coefficient.h>
|
|
#include <poincare/constant.h>
|
|
#include <poincare/cosine.h>
|
|
#include <poincare/division.h>
|
|
#include <poincare/float.h>
|
|
#include <poincare/horizontal_layout.h>
|
|
#include <poincare/infinity.h>
|
|
#include <poincare/matrix.h>
|
|
#include <poincare/matrix_identity.h>
|
|
#include <poincare/matrix_inverse.h>
|
|
#include <poincare/naperian_logarithm.h>
|
|
#include <poincare/nth_root.h>
|
|
#include <poincare/opposite.h>
|
|
#include <poincare/parenthesis.h>
|
|
#include <poincare/serialization_helper.h>
|
|
#include <poincare/sine.h>
|
|
#include <poincare/square_root.h>
|
|
#include <poincare/subtraction.h>
|
|
#include <poincare/symbol.h>
|
|
#include <poincare/undefined.h>
|
|
#include <poincare/unreal.h>
|
|
#include <poincare/vertical_offset_layout.h>
|
|
#include <ion.h>
|
|
#include <assert.h>
|
|
#include <cmath>
|
|
#include <math.h>
|
|
#include <utility>
|
|
|
|
namespace Poincare {
|
|
|
|
// Properties
|
|
ExpressionNode::Sign PowerNode::sign(Context * context) const {
|
|
if (Expression::ShouldStopProcessing()) {
|
|
return Sign::Unknown;
|
|
}
|
|
if (childAtIndex(0)->sign(context) == Sign::Positive && childAtIndex(1)->sign(context) != Sign::Unknown) {
|
|
return Sign::Positive;
|
|
}
|
|
if (childAtIndex(0)->sign(context) == Sign::Negative && childAtIndex(1)->type() == ExpressionNode::Type::Rational) {
|
|
RationalNode * r = static_cast<RationalNode *>(childAtIndex(1));
|
|
if (r->isInteger()) {
|
|
assert(!Integer::Division(r->signedNumerator(), Integer(2)).remainder.isOverflow());
|
|
if (Integer::Division(r->signedNumerator(), Integer(2)).remainder.isZero()) {
|
|
return Sign::Positive;
|
|
} else {
|
|
return Sign::Negative;
|
|
}
|
|
}
|
|
}
|
|
return Sign::Unknown;
|
|
}
|
|
|
|
Expression PowerNode::setSign(Sign s, ReductionContext reductionContext) {
|
|
assert(s == ExpressionNode::Sign::Positive);
|
|
return Power(this).setSign(s, reductionContext);
|
|
}
|
|
|
|
int PowerNode::polynomialDegree(Context * context, const char * symbolName) const {
|
|
int deg = ExpressionNode::polynomialDegree(context, symbolName);
|
|
if (deg == 0) {
|
|
return deg;
|
|
}
|
|
int op0Deg = childAtIndex(0)->polynomialDegree(context, symbolName);
|
|
if (op0Deg < 0) {
|
|
return -1;
|
|
}
|
|
if (childAtIndex(1)->type() == ExpressionNode::Type::Rational) {
|
|
RationalNode * r = static_cast<RationalNode *>(childAtIndex(1));
|
|
if (!r->isInteger() || Number(r).sign() == Sign::Negative) {
|
|
return -1;
|
|
}
|
|
Integer numeratorInt = r->signedNumerator();
|
|
if (Integer::NaturalOrder(numeratorInt, Integer(Integer::k_maxExtractableInteger)) > 0) {
|
|
return -1;
|
|
}
|
|
op0Deg *= numeratorInt.extractedInt();
|
|
return op0Deg;
|
|
}
|
|
return -1;
|
|
}
|
|
|
|
Expression PowerNode::getUnit() const {
|
|
return Power(this).getUnit();
|
|
}
|
|
|
|
int PowerNode::getPolynomialCoefficients(Context * context, const char * symbolName, Expression coefficients[], ExpressionNode::SymbolicComputation symbolicComputation) const {
|
|
return Power(this).getPolynomialCoefficients(context, symbolName, coefficients);
|
|
}
|
|
|
|
bool PowerNode::isReal(Context * context) const {
|
|
Expression base(childAtIndex(0));
|
|
Expression index(childAtIndex(1));
|
|
// Both base and index are real and:
|
|
// - either base > 0
|
|
// - or index is an integer
|
|
if (base.isReal(context) &&
|
|
index.isReal(context) &&
|
|
(base.sign(context) == Sign::Positive ||
|
|
(index.type() == ExpressionNode::Type::Rational && static_cast<Rational &>(index).isInteger()))) {
|
|
return true;
|
|
}
|
|
return false;
|
|
}
|
|
|
|
bool PowerNode::childAtIndexNeedsUserParentheses(const Expression & child, int childIndex) const {
|
|
if (childIndex == 0) {
|
|
|
|
if ((child.isNumber() && static_cast<const Number &>(child).sign() == Sign::Negative)
|
|
|| (child.type() == Type::Rational && !static_cast<const Rational &>(child).isInteger()))
|
|
{
|
|
/* ^(-2.3, 4) --> (-2.3)^{4}
|
|
* ^(2/3, 4) --> (2/3)^{4} */
|
|
return true;
|
|
}
|
|
if (child.type() == Type::Conjugate) {
|
|
return childAtIndexNeedsUserParentheses(child.childAtIndex(0), childIndex);
|
|
}
|
|
// ^(2+3,4) --> (2+3)^{4}
|
|
Type types[] = {Type::Power, Type::Subtraction, Type::Opposite, Type::Multiplication, Type::Division, Type::Addition};
|
|
return child.isOfType(types, 6);
|
|
}
|
|
return false;
|
|
}
|
|
|
|
// Private
|
|
|
|
template<typename T>
|
|
Complex<T> PowerNode::computeNotPrincipalRealRootOfRationalPow(const std::complex<T> c, T p, T q) {
|
|
// Assert p and q are in fact integers
|
|
assert(std::round(p) == p);
|
|
assert(std::round(q) == q);
|
|
/* Try to find a real root of c^(p/q) with p, q integers. We ignore cases
|
|
* where the principal root is real as these cases are handled generically
|
|
* later (for instance 1232^(1/8) which has a real principal root is not
|
|
* handled here). */
|
|
if (c.imag() == (T)0.0 && std::pow((T)-1.0, q) < (T)0.0) {
|
|
/* If c real and q odd integer (q odd if (-1)^q = -1), a real root does
|
|
* exist (which is not necessarily the principal root)!
|
|
* For q even integer, a real root does not necessarily exist (example:
|
|
* -2 ^(1/2)). */
|
|
std::complex<T> absc = c;
|
|
absc.real(std::fabs(absc.real()));
|
|
// compute |c|^(p/q) which is a real
|
|
Complex<T> absCPowD = PowerNode::compute(absc, std::complex<T>(p/q), Preferences::ComplexFormat::Real);
|
|
/* As q is odd, c^(p/q) = (sign(c)^(1/q))^p * |c|^(p/q)
|
|
* = sign(c)^p * |c|^(p/q)
|
|
* = -|c|^(p/q) iff c < 0 and p odd */
|
|
return c.real() < (T)0.0 && std::pow((T)-1.0, p) < (T)0.0 ? Complex<T>::Builder(-absCPowD.stdComplex()) : absCPowD;
|
|
}
|
|
return Complex<T>::Undefined();
|
|
}
|
|
|
|
template<typename T>
|
|
Complex<T> PowerNode::compute(const std::complex<T> c, const std::complex<T> d, Preferences::ComplexFormat complexFormat) {
|
|
std::complex<T> result;
|
|
if (c.imag() == (T)0.0 && d.imag() == (T)0.0 && c.real() != (T)0.0 && (c.real() > (T)0.0 || std::round(d.real()) == d.real())) {
|
|
/* pow: (R+, R) -> R+ (2^1.3 ~ 2.46)
|
|
* pow: (R-, N) -> R+ ((-2)^3 = -8)
|
|
* In these cases we rather use std::pow(double, double) because:
|
|
* - pow on complexes is not as precise as pow on double: for instance,
|
|
* pow(complex<double>(2.0,0.0), complex<double>(3.0,0.0) = complex(7.9999999999999982,0.0)
|
|
* and pow(2.0,3.0) = 8.0
|
|
* - Using complex pow, std::pow(2.0, 1000) = (INFINITY, NAN).
|
|
* Openbsd pow of a positive real and another real has a undefined
|
|
* imaginary when the real result is infinity.
|
|
* However, we exclude c == 0 because std:pow(0.0, 0.0) = 1.0 and we would
|
|
* rather have 0^0 = undef. */
|
|
result = std::complex<T>(std::pow(c.real(), d.real()));
|
|
} else {
|
|
result = std::pow(c, d);
|
|
}
|
|
/* Openbsd trigonometric functions are numerical implementation and thus are
|
|
* approximative.
|
|
* The error epsilon is ~1E-7 on float and ~1E-15 on double. In order to
|
|
* avoid weird results as e(i*pi) = -1+6E-17*i, we compute the argument of
|
|
* the result of c^d and if arg ~ 0 [Pi], we discard the residual imaginary
|
|
* part and if arg ~ Pi/2 [Pi], we discard the residual real part.
|
|
* Let's determine when the arg [Pi] (or arg [Pi/2]) is negligeable:
|
|
* With c = r*e^(iθ) and d = x+iy, c^d = r^x*e^(yθ)*e^i(yln(r)+xθ)
|
|
* so arg(c^d) = y*ln(r)+xθ.
|
|
* We consider that arg[π] is negligeable if it is negligeable compared to
|
|
* norm(d) = sqrt(x^2+y^2) and ln(r) = ln(norm(c)).*/
|
|
return Complex<T>::Builder(ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(result, c, d, false));
|
|
}
|
|
|
|
// Layout
|
|
|
|
Layout PowerNode::createLayout(Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const {
|
|
ExpressionNode * indiceOperand = childAtIndex(1);
|
|
HorizontalLayout result = HorizontalLayout::Builder();
|
|
result.addOrMergeChildAtIndex(childAtIndex(0)->createLayout(floatDisplayMode, numberOfSignificantDigits), 0, false);
|
|
result.addChildAtIndex(VerticalOffsetLayout::Builder(
|
|
indiceOperand->createLayout(floatDisplayMode, numberOfSignificantDigits),
|
|
VerticalOffsetLayoutNode::Position::Superscript),
|
|
result.numberOfChildren(),
|
|
result.numberOfChildren(),
|
|
nullptr);
|
|
return std::move(result);
|
|
}
|
|
|
|
// Serialize
|
|
|
|
bool PowerNode::childNeedsSystemParenthesesAtSerialization(const TreeNode * child) const {
|
|
if (static_cast<const ExpressionNode *>(child)->isNumber() && Number(static_cast<const NumberNode *>(child)).sign() == Sign::Negative) {
|
|
return true;
|
|
}
|
|
if (static_cast<const ExpressionNode *>(child)->type() == Type::Rational && !static_cast<const RationalNode *>(child)->isInteger()) {
|
|
return true;
|
|
}
|
|
Type types[] = {Type::Power, Type::Subtraction, Type::Opposite, Type::Multiplication, Type::Division, Type::Addition};
|
|
return static_cast<const ExpressionNode *>(child)->isOfType(types, 6);
|
|
}
|
|
|
|
int PowerNode::serialize(char * buffer, int bufferSize, Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const {
|
|
return SerializationHelper::Infix(this, buffer, bufferSize, floatDisplayMode, numberOfSignificantDigits, "^");
|
|
}
|
|
|
|
// Simplify
|
|
|
|
Expression PowerNode::shallowReduce(ReductionContext reductionContext) {
|
|
return Power(this).shallowReduce(reductionContext);
|
|
}
|
|
|
|
Expression PowerNode::shallowBeautify(ReductionContext reductionContext) {
|
|
return Power(this).shallowBeautify(reductionContext);
|
|
}
|
|
|
|
int PowerNode::simplificationOrderGreaterType(const ExpressionNode * e, bool ascending, bool canBeInterrupted, bool ignoreParentheses) const {
|
|
int baseComparison = SimplificationOrder(childAtIndex(0), e, ascending, canBeInterrupted, ignoreParentheses);
|
|
if (baseComparison != 0) {
|
|
return baseComparison;
|
|
}
|
|
Rational one = Rational::Builder(1);
|
|
return SimplificationOrder(childAtIndex(1), one.node(), ascending, canBeInterrupted, ignoreParentheses);
|
|
}
|
|
|
|
int PowerNode::simplificationOrderSameType(const ExpressionNode * e, bool ascending, bool canBeInterrupted, bool ignoreParentheses) const {
|
|
assert(e->numberOfChildren() > 0);
|
|
int baseComparison = SimplificationOrder(childAtIndex(0), e->childAtIndex(0), ascending, canBeInterrupted, ignoreParentheses);
|
|
if (baseComparison != 0) {
|
|
return baseComparison;
|
|
}
|
|
assert(e->numberOfChildren() > 1);
|
|
return SimplificationOrder(childAtIndex(1), e->childAtIndex(1), ascending, canBeInterrupted, ignoreParentheses);
|
|
}
|
|
|
|
Expression PowerNode::denominator(ReductionContext reductionContext) const {
|
|
return Power(this).denominator(reductionContext);
|
|
}
|
|
|
|
// Evaluation
|
|
template<typename T> MatrixComplex<T> PowerNode::computeOnComplexAndMatrix(const std::complex<T> c, const MatrixComplex<T> n, Preferences::ComplexFormat complexFormat) {
|
|
return MatrixComplex<T>::Undefined();
|
|
}
|
|
|
|
template<typename T> MatrixComplex<T> PowerNode::computeOnMatrixAndComplex(const MatrixComplex<T> m, const std::complex<T> d, Preferences::ComplexFormat complexFormat) {
|
|
if (m.numberOfRows() != m.numberOfColumns()) {
|
|
return MatrixComplex<T>::Undefined();
|
|
}
|
|
T power = Complex<T>::Builder(d).toScalar();
|
|
if (std::isnan(power) || std::isinf(power) || power != (int)power || std::fabs(power) > k_maxApproximatePowerMatrix) {
|
|
return MatrixComplex<T>::Undefined();
|
|
}
|
|
if (power < 0) {
|
|
MatrixComplex<T> inverse = m.inverse();
|
|
if (inverse.isUninitialized()) {
|
|
return MatrixComplex<T>::Undefined();
|
|
}
|
|
Complex<T> minusC = Complex<T>::Builder(-d);
|
|
MatrixComplex<T> result = PowerNode::computeOnMatrixAndComplex(inverse, minusC.stdComplex(), complexFormat);
|
|
return result;
|
|
}
|
|
MatrixComplex<T> result = MatrixComplex<T>::CreateIdentity(m.numberOfRows());
|
|
// TODO: implement a quick exponentiation
|
|
for (int k = 0; k < (int)power; k++) {
|
|
if (Expression::ShouldStopProcessing()) {
|
|
return MatrixComplex<T>::Undefined();
|
|
}
|
|
result = MultiplicationNode::computeOnMatrices<T>(result, m, complexFormat);
|
|
}
|
|
return result;
|
|
}
|
|
|
|
template<typename T> MatrixComplex<T> PowerNode::computeOnMatrices(const MatrixComplex<T> m, const MatrixComplex<T> n, Preferences::ComplexFormat complexFormat) {
|
|
return MatrixComplex<T>::Undefined();
|
|
}
|
|
|
|
template<typename T> Evaluation<T> PowerNode::templatedApproximate(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const {
|
|
/* Special case: c^(p/q) with p, q integers
|
|
* In real mode, c^(p/q) might have a real root which is not the principal
|
|
* root. We return this value in that case to avoid returning "unreal". */
|
|
if (complexFormat == Preferences::ComplexFormat::Real) {
|
|
Evaluation<T> base = childAtIndex(0)->approximate(T(), context, complexFormat, angleUnit);
|
|
if (base.type() != EvaluationNode<T>::Type::Complex) {
|
|
goto defaultApproximation;
|
|
}
|
|
std::complex<T> c = static_cast<Complex<T> &>(base).stdComplex();
|
|
T p = NAN;
|
|
T q = NAN;
|
|
// If the power has been reduced, we look for a rational index
|
|
if (childAtIndex(1)->type() == ExpressionNode::Type::Rational) {
|
|
const RationalNode * r = static_cast<const RationalNode *>(childAtIndex(1));
|
|
p = r->signedNumerator().approximate<T>();
|
|
q = r->denominator().approximate<T>();
|
|
}
|
|
/* If the power has been simplified (reduced + beautified), we look for an
|
|
* index of the for Division(Rational,Rational). */
|
|
if (childAtIndex(1)->type() == ExpressionNode::Type::Division && childAtIndex(1)->childAtIndex(0)->type() == ExpressionNode::Type::Rational && childAtIndex(1)->childAtIndex(1)->type() == ExpressionNode::Type::Rational) {
|
|
const RationalNode * pRat = static_cast<const RationalNode *>(childAtIndex(1)->childAtIndex(0));
|
|
const RationalNode * qRat = static_cast<const RationalNode *>(childAtIndex(1)->childAtIndex(1));
|
|
if (!pRat->denominator().isOne() || !qRat->denominator().isOne()) {
|
|
goto defaultApproximation;
|
|
}
|
|
p = pRat->signedNumerator().approximate<T>();
|
|
q = qRat->signedNumerator().approximate<T>();
|
|
}
|
|
/* We don't handle power that haven't been reduced or simplified as the
|
|
* index can take to many forms and still be equivalent to p/q,
|
|
* with p, q integers. */
|
|
if (std::isnan(p) || std::isnan(q)) {
|
|
goto defaultApproximation;
|
|
}
|
|
Complex<T> result = computeNotPrincipalRealRootOfRationalPow(c, p, q);
|
|
if (!result.isUndefined()) {
|
|
return std::move(result);
|
|
}
|
|
}
|
|
defaultApproximation:
|
|
return ApproximationHelper::MapReduce<T>(this, context, complexFormat, angleUnit, compute<T>, computeOnComplexAndMatrix<T>, computeOnMatrixAndComplex<T>, computeOnMatrices<T>);
|
|
}
|
|
|
|
// Power
|
|
|
|
Expression Power::setSign(ExpressionNode::Sign s, ExpressionNode::ReductionContext reductionContext) {
|
|
assert(s == ExpressionNode::Sign::Positive);
|
|
if (childAtIndex(0).sign(reductionContext.context()) == ExpressionNode::Sign::Negative) {
|
|
Expression result = Power::Builder(childAtIndex(0).setSign(ExpressionNode::Sign::Positive, reductionContext), childAtIndex(1));
|
|
replaceWithInPlace(result);
|
|
return result.shallowReduce(reductionContext);
|
|
}
|
|
return *this;
|
|
}
|
|
|
|
int Power::getPolynomialCoefficients(Context * context, const char * symbolName, Expression coefficients[]) const {
|
|
int deg = polynomialDegree(context, symbolName);
|
|
if (deg <= 0) {
|
|
return Expression::defaultGetPolynomialCoefficients(context, symbolName, coefficients);
|
|
}
|
|
/* Here we only consider the case x^4 as privateGetPolynomialCoefficients is
|
|
* supposed to be called after reducing the expression. */
|
|
if (childAtIndex(0).type() == ExpressionNode::Type::Symbol
|
|
&& strcmp(childAtIndex(0).convert<Symbol>().name(), symbolName) == 0
|
|
&& childAtIndex(1).type() == ExpressionNode::Type::Rational)
|
|
{
|
|
Rational r = childAtIndex(1).convert<Rational>();
|
|
if (!r.isInteger() || r.sign() == ExpressionNode::Sign::Negative) {
|
|
return -1;
|
|
}
|
|
Integer num = r.unsignedIntegerNumerator();
|
|
if (Integer::NaturalOrder(num, Integer(Integer::k_maxExtractableInteger)) > 0) {
|
|
return -1;
|
|
}
|
|
int n = num.extractedInt();
|
|
if (n <= k_maxPolynomialDegree) {
|
|
for (int i = 0; i < n; i++) {
|
|
coefficients[i] = Rational::Builder(0);
|
|
}
|
|
coefficients[n] = Rational::Builder(1);
|
|
return n;
|
|
}
|
|
}
|
|
return -1;
|
|
}
|
|
|
|
Expression Power::shallowReduce(ExpressionNode::ReductionContext reductionContext) {
|
|
|
|
{
|
|
Expression e = Expression::defaultShallowReduce();
|
|
if (e.isUndefined()) {
|
|
return e;
|
|
}
|
|
}
|
|
|
|
Expression base = childAtIndex(0);
|
|
Expression index = childAtIndex(1);
|
|
|
|
// Step 1: There should be no unit in the index!
|
|
if (index.hasUnit()) {
|
|
return replaceWithUndefinedInPlace();
|
|
}
|
|
|
|
// Step 2: Handle matrices
|
|
if (index.deepIsMatrix(reductionContext.context())) {
|
|
return replaceWithUndefinedInPlace();
|
|
}
|
|
ExpressionNode::Type baseType = base.type();
|
|
ExpressionNode::Type indexType = index.type();
|
|
if (base.deepIsMatrix(reductionContext.context())) {
|
|
if (indexType != ExpressionNode::Type::Rational || !static_cast<Rational &>(index).isInteger()) {
|
|
return replaceWithUndefinedInPlace();
|
|
}
|
|
if (baseType != ExpressionNode::Type::Matrix) {
|
|
return *this;
|
|
}
|
|
Matrix matrixBase = static_cast<Matrix &>(base);
|
|
if (matrixBase.numberOfRows() != matrixBase.numberOfColumns()) {
|
|
return replaceWithUndefinedInPlace();
|
|
}
|
|
Integer exponent = static_cast<Rational &>(index).signedIntegerNumerator();
|
|
if (exponent.isNegative()) {
|
|
index.setSign(ExpressionNode::Sign::Positive, reductionContext);
|
|
Expression reducedPositiveExponentMatrix = shallowReduce(reductionContext);
|
|
Expression dummyExpression = Undefined::Builder();
|
|
MatrixInverse inv = MatrixInverse::Builder(dummyExpression);
|
|
reducedPositiveExponentMatrix.replaceWithInPlace(inv);
|
|
inv.replaceChildInPlace(dummyExpression, reducedPositiveExponentMatrix);
|
|
return inv.shallowReduce(reductionContext);
|
|
}
|
|
if (Integer::NaturalOrder(exponent, Integer(k_maxExactPowerMatrix)) > 0) {
|
|
return *this;
|
|
}
|
|
int exp = exponent.extractedInt(); // Ok, because 0 < exponent < k_maxExactPowerMatrix
|
|
if (exp == 0) {
|
|
Matrix id = Matrix::CreateIdentity(matrixBase.numberOfRows());
|
|
replaceWithInPlace(id);
|
|
return std::move(id);
|
|
}
|
|
if (exp == 1) {
|
|
replaceWithInPlace(matrixBase);
|
|
return std::move(matrixBase);
|
|
}
|
|
Expression result = matrixBase.clone();
|
|
// TODO: implement a quick exponentiation
|
|
for (int k = 1; k < exp; k++) {
|
|
result = Multiplication::Builder(result, matrixBase.clone());
|
|
result = result.shallowReduce(reductionContext);
|
|
}
|
|
assert(!result.isUninitialized());
|
|
replaceWithInPlace(result);
|
|
return result;
|
|
}
|
|
|
|
Expression power = *this;
|
|
/* Step 3: if both children are true unresolved complexes, the result is not
|
|
* simplified. TODO? */
|
|
if (!base.isReal(reductionContext.context())
|
|
&& baseType != ExpressionNode::Type::ComplexCartesian
|
|
&& !index.isReal(reductionContext.context())
|
|
&& indexType != ExpressionNode::Type::ComplexCartesian)
|
|
{
|
|
return *this;
|
|
}
|
|
|
|
/* Step 4: Handle simple cases as x^0, x^1, 0^x and 1^x first, for 2 reasons:
|
|
* - we can assert after this step that there is no division by 0:
|
|
* for instance, 0^(-2)->undefined
|
|
* - we save computational time by early escaping for these cases. */
|
|
if (indexType == ExpressionNode::Type::Rational) {
|
|
const Rational rationalIndex = static_cast<Rational &>(index);
|
|
// x^0
|
|
if (rationalIndex.isZero()) {
|
|
// 0^0 = undef or (±inf)^0 = undef
|
|
if (base.isRationalZero() || baseType == ExpressionNode::Type::Infinity) {
|
|
return replaceWithUndefinedInPlace();
|
|
}
|
|
// x^0
|
|
if (reductionContext.target() == ExpressionNode::ReductionTarget::User
|
|
|| base.isNumber()
|
|
|| baseType == ExpressionNode::Type::Constant)
|
|
{
|
|
/* Warning: if the ReductionTarget is User, in all other cases but 0^0,
|
|
* we replace x^0 by one. This is almost always true except when x = 0.
|
|
* However, not substituting x^0 by one would prevent from simplifying
|
|
* many expressions like x/x->1. */
|
|
Expression result = Rational::Builder(1);
|
|
replaceWithInPlace(result);
|
|
return result;
|
|
}
|
|
}
|
|
// x^1
|
|
if (rationalIndex.isOne()) {
|
|
Expression result = childAtIndex(0);
|
|
replaceWithInPlace(result);
|
|
return result;
|
|
}
|
|
}
|
|
if (baseType == ExpressionNode::Type::Rational) {
|
|
Rational rationalBase = static_cast<Rational &>(base);
|
|
// 0^x
|
|
if (rationalBase.isZero()) {
|
|
// 0^x with x > 0 = 0
|
|
ExpressionNode::Sign indexSign = index.sign(reductionContext.context());
|
|
if (indexSign == ExpressionNode::Sign::Positive) {
|
|
Expression result = Rational::Builder(0);
|
|
replaceWithInPlace(result);
|
|
return result;
|
|
} else if (indexSign == ExpressionNode::Sign::Negative) {
|
|
// 0^x with x < 0 = undef
|
|
return replaceWithUndefinedInPlace();
|
|
}
|
|
}
|
|
// 1^x = 1 if x != ±inf
|
|
if (rationalBase.isOne() && !index.recursivelyMatches(Expression::IsInfinity, reductionContext.context())) {
|
|
Expression result = Rational::Builder(1);
|
|
replaceWithInPlace(result);
|
|
return result;
|
|
}
|
|
}
|
|
|
|
/* We do not apply some rules to a^b if the parent node is a logarithm of same
|
|
* base a. In this case there is a simplication of form
|
|
* ln(e^(3^(1/2))->3^(1/2). */
|
|
bool letPowerAtRoot = parentIsALogarithmOfSameBase();
|
|
|
|
/* Step 5: we now bubble up ComplexCartesian, we handle different cases:
|
|
* At least, one child is a ComplexCartesian and the other is either a
|
|
* ComplexCartesian or real. */
|
|
|
|
ComplexCartesian complexBase;
|
|
ComplexCartesian complexIndex;
|
|
ComplexCartesian result;
|
|
/* First, (x+iy)^q with q special values
|
|
* For q = -1, 1/2, -1/2, n with n integer < 10, we avoid introducing arctangent
|
|
* by using the formula (r*e^(i*th))^(a+ib) = r^a*e(-th*b)*e^(b*ln(r)+th*a).
|
|
* Instead, we rather use the cartesian form of the base and the index. */
|
|
if (!letPowerAtRoot && baseType == ExpressionNode::Type::ComplexCartesian) {
|
|
complexBase = static_cast<ComplexCartesian &>(base);
|
|
if (indexType == ExpressionNode::Type::Rational) {
|
|
Rational r = static_cast<Rational &>(index);
|
|
if (r.isMinusOne()) {
|
|
// (x+iy)^(-1)
|
|
result = complexBase.inverse(reductionContext);
|
|
} else if (r.isHalf()) {
|
|
// (x+iy)^(1/2)
|
|
result = complexBase.squareRoot(reductionContext);
|
|
} else if (r.isMinusHalf()) {
|
|
// (x+iy)^(-1/2)
|
|
result = complexBase.squareRoot(reductionContext).inverse(reductionContext);
|
|
} else if (r.isInteger() && r.unsignedIntegerNumerator().isLowerThan(Integer(10))) {
|
|
if (r.sign() == ExpressionNode::Sign::Positive) {
|
|
// (x+iy)^n, n integer positive n < 10
|
|
result = complexBase.powerInteger(r.unsignedIntegerNumerator().extractedInt(), reductionContext);
|
|
} else {
|
|
// (x+iy)^(-n), n integer positive n < 10
|
|
assert(r.sign() == ExpressionNode::Sign::Negative);
|
|
result = complexBase.powerInteger(r.unsignedIntegerNumerator().extractedInt(), reductionContext).inverse(reductionContext);
|
|
}
|
|
}
|
|
if (!result.isUninitialized()) {
|
|
replaceWithInPlace(result);
|
|
return result.shallowReduce();
|
|
}
|
|
}
|
|
}
|
|
// All other cases where one child at least is a ComplexCartesian
|
|
if (!letPowerAtRoot
|
|
&& ((base.isReal(reductionContext.context())
|
|
&& indexType == ExpressionNode::Type::ComplexCartesian)
|
|
|| (baseType == ExpressionNode::Type::ComplexCartesian
|
|
&& index.isReal(reductionContext.context()))
|
|
|| (baseType == ExpressionNode::Type::ComplexCartesian
|
|
&& indexType == ExpressionNode::Type::ComplexCartesian)))
|
|
{
|
|
complexBase = baseType == ExpressionNode::Type::ComplexCartesian ? static_cast<ComplexCartesian &>(base) : ComplexCartesian::Builder(base, Rational::Builder(0));
|
|
complexIndex = indexType == ExpressionNode::Type::ComplexCartesian ? static_cast<ComplexCartesian &>(index) : ComplexCartesian::Builder(index, Rational::Builder(0));
|
|
result = complexBase.power(complexIndex, reductionContext);
|
|
replaceWithInPlace(result);
|
|
return result.shallowReduce();
|
|
}
|
|
|
|
/* Step 6: We look for square root and sum of square roots (two terms maximum
|
|
* so far) at the denominator and move them to the numerator. */
|
|
if (reductionContext.target() == ExpressionNode::ReductionTarget::User) {
|
|
Expression r = removeSquareRootsFromDenominator(reductionContext);
|
|
if (!r.isUninitialized()) {
|
|
return r;
|
|
}
|
|
}
|
|
|
|
#if 0
|
|
/* This was implemented when complex were not handled with the
|
|
* ComplexCartesian node. It was used for instance to turn
|
|
* sqrt(i) -> i^(1/2)
|
|
* -> e(i*π/4) HERE!
|
|
* -> cos(π/4)+i*sin(π/4)
|
|
* -> sqrt(2)/2+i*sqrt(2)/2
|
|
*
|
|
* However, now, ComplexCartesian class handles power of ComplexCartesian
|
|
* when the target of simplification is the User. When the target of
|
|
* simplification is the system, i^(p/q) is a shorter expression than
|
|
* e^(i*Pi*p/2q) so we rather keep it in the first form.
|
|
* We thereby don't need the rule anymore!
|
|
* */
|
|
/* Step 4: we simplify i^(p/q) = e^(i*Pi*p/2q) */
|
|
/*if (indexType == ExpressionNode::Type::Rational) {
|
|
const Rational rationalIndex = static_cast<Rational &>(index);
|
|
// i^(p/q)
|
|
if (baseType == ExpressionNode::Type::Constant && static_cast<Constant &>(base).isIComplex()) {
|
|
Number r = Number::Multiplication(rationalIndex, Rational::Builder(1, 2));
|
|
Expression result = CreateComplexExponent(r, reductionContext);
|
|
replaceWithInPlace(result);
|
|
return result.shallowReduce(reductionContext);
|
|
}
|
|
}*/
|
|
#endif
|
|
|
|
// Step 7: (±inf)^x = 0 or ±inf
|
|
if (baseType == ExpressionNode::Type::Infinity) {
|
|
Expression result;
|
|
ExpressionNode::Sign indexSign = index.sign(reductionContext.context());
|
|
if (indexSign == ExpressionNode::Sign::Negative) {
|
|
// --> 0 if x < 0
|
|
result = Rational::Builder(0);
|
|
} else if (indexSign == ExpressionNode::Sign::Positive) {
|
|
// --> (±inf) if x > 0
|
|
result = Infinity::Builder(false);
|
|
if (base.sign(reductionContext.context()) == ExpressionNode::Sign::Negative) {
|
|
// (-inf)^x --> (-1)^x*inf
|
|
Power p = Power::Builder(Rational::Builder(-1), childAtIndex(1));
|
|
result = Multiplication::Builder(p, result);
|
|
p.shallowReduce(reductionContext);
|
|
}
|
|
}
|
|
if (!result.isUninitialized()) {
|
|
replaceWithInPlace(result);
|
|
return result.shallowReduce(reductionContext);
|
|
}
|
|
}
|
|
|
|
// Step 8: p^q with p, q rationals --> a*b^c*exp(i*pi*d) with a, b, c, d rationals
|
|
if (!letPowerAtRoot && baseType == ExpressionNode::Type::Rational) {
|
|
Rational rationalBase = static_cast<Rational &>(base);
|
|
// p^q with p, q rationals
|
|
// TODO: maybe reduce Number^Rational?
|
|
if (indexType == ExpressionNode::Type::Rational) {
|
|
Rational rationalIndex = static_cast<Rational &>(index);
|
|
if (RationalExponentShouldNotBeReduced(rationalBase, rationalIndex)) {
|
|
return *this;
|
|
}
|
|
return simplifyRationalRationalPower(reductionContext);
|
|
}
|
|
}
|
|
/* Step 9: (a)^(1/2) --> i*(-a)^(1/2)
|
|
* WARNING: this rule true only if:
|
|
* - a real: (-1*i)^(1/2) != i*i^(1/2)
|
|
* - a is negative: (-(-2))^(1/2) != -2^(1/2)
|
|
* We apply this rule only when a is a negative numeral */
|
|
if (!letPowerAtRoot
|
|
&& base.isNumber()
|
|
&& indexType == ExpressionNode::Type::Rational
|
|
&& index.convert<Rational>().isHalf())
|
|
{
|
|
Expression m0 = base.makePositiveAnyNegativeNumeralFactor(reductionContext);
|
|
if (!m0.isUninitialized()) {
|
|
replaceChildAtIndexInPlace(0, m0);
|
|
/* m0 doesn't need to be shallowReduce as
|
|
* makePositiveAnyNegativeNumeralFactor returns a reduced expression */
|
|
Multiplication m1 = Multiplication::Builder();
|
|
replaceWithInPlace(m1);
|
|
// Multiply m1 by i complex
|
|
Constant i = Constant::Builder(UCodePointMathematicalBoldSmallI);
|
|
m1.addChildAtIndexInPlace(i, 0, 0);
|
|
i.shallowReduce(reductionContext);
|
|
m1.addChildAtIndexInPlace(*this, 1, 1);
|
|
shallowReduce(reductionContext);
|
|
return m1.shallowReduce(reductionContext);
|
|
}
|
|
}
|
|
// Step 10: e^(r*i*Pi) with r rational --> cos(pi*r) + i*sin(pi*r)
|
|
if (!letPowerAtRoot && isNthRootOfUnity()) {
|
|
Expression i = index.childAtIndex(index.numberOfChildren()-2);
|
|
static_cast<Multiplication &>(index).removeChildAtIndexInPlace(index.numberOfChildren()-2);
|
|
if (reductionContext.angleUnit() == Preferences::AngleUnit::Degree) {
|
|
index.replaceChildAtIndexInPlace(index.numberOfChildren()-1, Rational::Builder(180));
|
|
}
|
|
Expression cos = Cosine::Builder(index);
|
|
index = index.shallowReduce(reductionContext);
|
|
Expression sin = Sine::Builder(index.clone());
|
|
Expression complexPart = Multiplication::Builder(sin, i);
|
|
sin.shallowReduce(reductionContext);
|
|
Expression a = Addition::Builder(cos, complexPart);
|
|
cos.shallowReduce(reductionContext);
|
|
complexPart.shallowReduce(reductionContext);
|
|
replaceWithInPlace(a);
|
|
return a.shallowReduce(reductionContext);
|
|
}
|
|
// Step 11: x^log(y,x)->y if y > 0
|
|
if (indexType == ExpressionNode::Type::Logarithm) {
|
|
if (index.numberOfChildren() == 2 && base.isIdenticalTo(index.childAtIndex(1))) {
|
|
// y > 0
|
|
if (index.childAtIndex(0).sign(reductionContext.context()) == ExpressionNode::Sign::Positive) {
|
|
Expression result = index.childAtIndex(0);
|
|
replaceWithInPlace(result);
|
|
return result;
|
|
}
|
|
}
|
|
// 10^log(y)
|
|
if (index.numberOfChildren() == 1
|
|
&& baseType == ExpressionNode::Type::Rational
|
|
&& static_cast<Rational &>(base).isTen())
|
|
{
|
|
Expression result = index.childAtIndex(0);
|
|
replaceWithInPlace(result);
|
|
return result;
|
|
}
|
|
}
|
|
/* Step 12: (a^b)^c -> a^(b*c)
|
|
* This rule is not generally true: ((-2)^2)^(1/2) != (-2)^(2*1/2) = -2
|
|
* This rule is true if:
|
|
* - a > 0
|
|
* - c is an integer
|
|
*
|
|
* Warning 1: in real mode only c integer is not enough:
|
|
* ex: ((-2)^(1/2))^2 = unreal != -2
|
|
* We escape that case by returning 'unreal' if the a^b is complex.
|
|
*
|
|
* Warning 2: If we did not apply this rule on expressions of the form
|
|
* (a^b)^(-1), we would end up in infinite loop when factorizing an addition
|
|
* on the same denominator.
|
|
* For ex:
|
|
* 1+[tan(2)^1/2]^(-1) --> (tan(2)^1/2+tan(2)^1/2*[tan(2)^1/2]^(-1))/tan(2)^1/2
|
|
* --> tan(2)+tan(2)*[tan(2)^1/2]^(-1)/tan(2)
|
|
* --> tan(2)^(3/2)+tan(2)^(3/2)*[tan(2)^1/2]^(-1)/tan(2)^3/2
|
|
* --> ...
|
|
* Indeed, we have to apply the rule (a^b)^c -> a^(b*c) as soon as c is -1.
|
|
*/
|
|
if (baseType == ExpressionNode::Type::Power) {
|
|
Power powerBase = static_cast<Power &>(base);
|
|
|
|
bool cInteger = indexType == ExpressionNode::Type::Rational && static_cast<Rational &>(index).isInteger();
|
|
bool applyRule = powerBase.childAtIndex(0).sign(reductionContext.context()) == ExpressionNode::Sign::Positive // a > 0
|
|
|| cInteger; // c integer
|
|
bool cMinusOne = cInteger && static_cast<Rational &>(index).isMinusOne();
|
|
/* If the complexFormat is real, we check that the inner power is defined
|
|
* before applying the rule (a^b)^c -> a^(b*c). Otherwise, we return
|
|
* 'unreal' or we do nothing.
|
|
* We escape this additional check if c = -1 for two reasons:
|
|
* - (a^b)^(-1) has to be reduced to avoid infinite loop discussed above;
|
|
* - if a^b is unreal, a^(-b) also. */
|
|
if (!cMinusOne && reductionContext.complexFormat() == Preferences::ComplexFormat::Real) {
|
|
Expression approximation = powerBase.approximate<float>(reductionContext.context(), reductionContext.complexFormat(), reductionContext.angleUnit());
|
|
if (approximation.type() == ExpressionNode::Type::Unreal) {
|
|
// The inner power is unreal, return "unreal"
|
|
replaceWithInPlace(approximation);
|
|
return approximation;
|
|
} else if (approximation.type() == ExpressionNode::Type::Undefined && !powerBase.isReal(reductionContext.context())) { //TODO Quite dirty way to check this. We want to forbid x, but not a unit.
|
|
/* The inner power is undefined, it can be 'x^(1/2)' for instance. We
|
|
* don't want to simplify this as it could be unreal with x = -2 but
|
|
* also real with x = 2. */
|
|
applyRule = false;
|
|
}
|
|
}
|
|
if (applyRule) {
|
|
return simplifyPowerPower(reductionContext);
|
|
}
|
|
}
|
|
|
|
// Step 13: (a*b*c*...)^r ?
|
|
if (!letPowerAtRoot && baseType == ExpressionNode::Type::Multiplication) {
|
|
Multiplication multiplicationBase = static_cast<Multiplication &>(base);
|
|
// Case 1: (a*b*c*...)^n = a^n*b^n*c^n*... if n integer
|
|
if (indexType == ExpressionNode::Type::Rational && static_cast<Rational &>(index).isInteger()) {
|
|
return simplifyPowerMultiplication(reductionContext);
|
|
}
|
|
// Case 2: (a*b*...)^r -> |a|^r*(sign(a)*b*...)^r if a not -1
|
|
for (int i = 0; i < multiplicationBase.numberOfChildren(); i++) {
|
|
// a is signed and a != -1
|
|
Expression baseChildI = multiplicationBase.childAtIndex(i);
|
|
ExpressionNode::Sign baseChildISign = baseChildI.sign(reductionContext.context());
|
|
if (baseChildISign != ExpressionNode::Sign::Unknown
|
|
&& (baseChildI.type() != ExpressionNode::Type::Rational
|
|
|| !static_cast<Rational &>(baseChildI).isMinusOne()))
|
|
{
|
|
Expression r = index;
|
|
Expression rCopy = r.clone();
|
|
// (sign(a)*b*...)^r
|
|
if (baseChildISign == ExpressionNode::Sign::Negative) {
|
|
multiplicationBase.replaceChildAtIndexInPlace(i, Rational::Builder(-1));
|
|
baseChildI = baseChildI.setSign(ExpressionNode::Sign::Positive, reductionContext);
|
|
} else {
|
|
multiplicationBase.removeChildAtIndexInPlace(i);
|
|
}
|
|
multiplicationBase.shallowReduce(reductionContext);
|
|
|
|
// |a|^r
|
|
Power p = Power::Builder(baseChildI, rCopy);
|
|
|
|
// |a|^r*(sign(a)*b*...)^r
|
|
Power thisRef = *this;
|
|
Multiplication root = Multiplication::Builder(p);
|
|
replaceWithInPlace(root);
|
|
root.addChildAtIndexInPlace(thisRef, 1, 1);
|
|
p.shallowReduce(reductionContext);
|
|
thisRef.shallowReduce(reductionContext);
|
|
return root.shallowReduce(reductionContext);
|
|
}
|
|
}
|
|
}
|
|
/* Step 14: a^(p/q+c+...) -> Rational::Builder(a^p)*a^(1/q+c+...) with a
|
|
* rational and a != 0 and p, q integers */
|
|
if (!letPowerAtRoot
|
|
&& baseType == ExpressionNode::Type::Rational
|
|
&& !static_cast<Rational &>(base).isZero()
|
|
&& indexType == ExpressionNode::Type::Addition)
|
|
{
|
|
Addition additionIndex = static_cast<Addition &>(index);
|
|
// Check is b is rational
|
|
Expression additionIndexChild0 = additionIndex.childAtIndex(0);
|
|
if (additionIndexChild0.type() == ExpressionNode::Type::Rational) {
|
|
const Rational rationalIndex = static_cast<Rational &>(additionIndexChild0);
|
|
if (rationalIndex.unsignedIntegerNumerator().isOne() && !rationalIndex.isInteger()) {
|
|
/* Do not reduce a^(1/q+c+...) to avoid potential infinite loop:
|
|
* a^(1/q+c+...) --> a^(1/q)*a^(c+...) --> a^(1/q+c+...)*/
|
|
/* TODO: do something more sensible here:
|
|
* - add rule (-rational)^x --> (-1)^x*rational^x so we only consider
|
|
* positive rational or (-1)
|
|
* - change simplifyRationalRationalPower to be able to detect when no
|
|
* rational was extracted (ie 2^(1/2) --> 2^(1/2)) to avoid applying
|
|
* this rule in that case
|
|
* Once this is done, we can reduce 4^(1/2+a) --> 2*4^a which is not
|
|
* done so far to avoir the infinite loop:
|
|
* 2^(1/2+s) --> 2^(1/2)*2^s --> 2^(1/2+s)... */
|
|
return *this;
|
|
}
|
|
const Rational rationalBase = static_cast<Rational &>(base);
|
|
if (RationalExponentShouldNotBeReduced(rationalBase, rationalIndex)) {
|
|
return *this;
|
|
}
|
|
Power p1 = Power::Builder(childAtIndex(0).clone(), additionIndex.childAtIndex(0));
|
|
Power thisRef = *this;
|
|
additionIndex.removeChildAtIndexInPlace(0); // p2 = a^(c+...)
|
|
// if addition had only 2 children
|
|
additionIndex.squashUnaryHierarchyInPlace();
|
|
Multiplication m = Multiplication::Builder(p1);
|
|
replaceWithInPlace(m);
|
|
m.addChildAtIndexInPlace(thisRef, 1, 1);
|
|
p1.simplifyRationalRationalPower(reductionContext);
|
|
return m.shallowReduce(reductionContext);
|
|
}
|
|
}
|
|
|
|
/* Step 15: (a0+a1+...am)^n with n integer
|
|
* -> a^n+?a^(n-1)*b+?a^(n-2)*b^2+...+b^n (Multinome)
|
|
* We don't apply this rule when the target is the SystemForApproximation.
|
|
* Indeed, developing the multinome is likely to increase the numbers of
|
|
* operations and lead to precision loss. */
|
|
if (!letPowerAtRoot
|
|
&& reductionContext.target() != ExpressionNode::ReductionTarget::SystemForApproximation
|
|
&& indexType == ExpressionNode::Type::Rational
|
|
&& !static_cast<Rational &>(index).signedIntegerNumerator().isZero()
|
|
&& static_cast<Rational &>(index).isInteger()
|
|
&& baseType == ExpressionNode::Type::Addition)
|
|
{
|
|
// Exponent n
|
|
Rational rationalIndex = static_cast<Rational &>(index);
|
|
Integer n = rationalIndex.unsignedIntegerNumerator();
|
|
/* If n is above 25, the resulting sum would have more than
|
|
* k_maxNumberOfTermsInExpandedMultinome terms so we do not expand it. */
|
|
if (Integer(k_maxNumberOfTermsInExpandedMultinome).isLowerThan(n) || n.isOne()) {
|
|
return *this;
|
|
}
|
|
int clippedN = n.extractedInt(); // Authorized because n < k_maxNumberOfTermsInExpandedMultinome
|
|
assert(clippedN > 0);
|
|
|
|
// Number of terms in addition m
|
|
int m = base.numberOfChildren();
|
|
/* The multinome (a0+a2+...+a(m-1))^n has BinomialCoefficient(n+m-1,n) terms;
|
|
* we expand the multinome only when the number of terms in the resulting
|
|
* sum has less than k_maxNumberOfTermsInExpandedMultinome terms. */
|
|
if (k_maxNumberOfTermsInExpandedMultinome < BinomialCoefficientNode::compute(static_cast<double>(clippedN), static_cast<double>(clippedN+m-1))) {
|
|
return *this;
|
|
}
|
|
|
|
Expression result = base;
|
|
Expression a = result.clone();
|
|
for (int i = 2; i <= clippedN; i++) {
|
|
// result = result * (a0+a1+...+a(m-1) in its expanded form
|
|
if (result.type() == ExpressionNode::Type::Addition) {
|
|
// We need a 'double' distribution and newA will hold the new expanded form
|
|
Expression newA = Addition::Builder();
|
|
for (int j = 0; j < a.numberOfChildren(); j++) {
|
|
Expression m = Multiplication::Builder(result.clone(), a.childAtIndex(j).clone()).distributeOnOperandAtIndex(0, reductionContext);
|
|
if (newA.type() == ExpressionNode::Type::Addition) {
|
|
static_cast<Addition &>(newA).addChildAtIndexInPlace(m, newA.numberOfChildren(), newA.numberOfChildren());
|
|
} else {
|
|
newA = Addition::Builder(newA, m);
|
|
}
|
|
newA = newA.shallowReduce(reductionContext);
|
|
}
|
|
result.replaceWithInPlace(newA);
|
|
result = newA;
|
|
} else {
|
|
// Just distribute result on a
|
|
Multiplication m = Multiplication::Builder(a.clone(), result.clone());
|
|
Expression distributedM = m.distributeOnOperandAtIndex(0, reductionContext);
|
|
result.replaceWithInPlace(distributedM);
|
|
result = distributedM;
|
|
result = result.shallowReduce(reductionContext);
|
|
}
|
|
}
|
|
if (rationalIndex.sign() == ExpressionNode::Sign::Negative) {
|
|
rationalIndex.replaceWithInPlace(Rational::Builder(-1));
|
|
return shallowReduce(reductionContext);
|
|
} else {
|
|
replaceWithInPlace(result);
|
|
return result;
|
|
}
|
|
}
|
|
#if 0
|
|
/* We could use the Newton formula instead which is quicker but not immediate
|
|
* to implement in the general case (Newton multinome). */
|
|
// (a+b)^n with n integer -> a^n+?a^(n-1)*b+?a^(n-2)*b^2+...+b^n (Newton)
|
|
if (!letPowerAtRoot && childAtIndex(1)->type() == ExpressionNode::Type::Rational && static_cast<const Rational *>(childAtIndex(1))->isInteger() && childAtIndex(0)->type() == ExpressionNode::Type::Addition && childAtIndex(0)->numberOfChildren() == 2) {
|
|
Rational * nr = static_cast<Rational *>(childAtIndex(1));
|
|
Integer n = nr->numerator();
|
|
n.setNegative(false);
|
|
if (Integer(k_maxExpandedBinome).isLowerThan(n) || n.isOne()) {
|
|
return this;
|
|
}
|
|
int clippedN = n.extractedInt(); // Authorized because n < k_maxExpandedBinome < k_maxNValue
|
|
Expression * x0 = childAtIndex(0)->childAtIndex(0);
|
|
Expression * x1 = childAtIndex(0)->childAtIndex(1);
|
|
Addition * a = new Addition::Builder();
|
|
for (int i = 0; i <= clippedN; i++) {
|
|
Rational * r = new Rational::Builder(static_cast<int>(BinomialCoefficient::compute(static_cast<double>(i), static_cast<double>(clippedN))));
|
|
Power * p0 = new Power::Builder(x0->clone(), new Rational::Builder(i), false);
|
|
Power * p1 = new Power::Builder(x1->clone(), new Rational::Builder(clippedN-i), false);
|
|
const Expression * operands[3] = {r, p0, p1};
|
|
Multiplication * m = new Multiplication::Builder(operands, 3, false);
|
|
p0->shallowReduce(reductionContext);
|
|
p1->shallowReduce(reductionContext);
|
|
a->addOperand(m);
|
|
m->shallowReduce(reductionContext);
|
|
}
|
|
if (nr->sign(&context) == Sign::Negative) {
|
|
nr->replaceWith(new Rational::Builder(-1), true);
|
|
childAtIndex(0)->replaceWith(a, true)->shallowReduce(reductionContext);
|
|
return shallowReduce(reductionContext);
|
|
} else {
|
|
return replaceWith(a, true)->shallowReduce(reductionContext);
|
|
}
|
|
}
|
|
#endif
|
|
return *this;
|
|
}
|
|
|
|
Expression Power::shallowBeautify(ExpressionNode::ReductionContext reductionContext) {
|
|
// Step 1: X^-y -> 1/(X->shallowBeautify)^y
|
|
Expression p = denominator(reductionContext);
|
|
// If the denominator is initialized, the index of the power is of form -y
|
|
if (!p.isUninitialized()) {
|
|
Division d = Division::Builder(Rational::Builder(1), p);
|
|
replaceWithInPlace(d);
|
|
p.shallowReduce(reductionContext);
|
|
return d.shallowBeautify(reductionContext);
|
|
}
|
|
// Step 2: Turn a^(1/n) into root(a, n)
|
|
if (childAtIndex(1).type() == ExpressionNode::Type::Rational && childAtIndex(1).convert<Rational>().signedIntegerNumerator().isOne()) {
|
|
Integer index = childAtIndex(1).convert<Rational>().integerDenominator();
|
|
// Special case: a^(1/2) --> sqrt(a)
|
|
if (index.isEqualTo(Integer(2))) {
|
|
Expression result = SquareRoot::Builder(childAtIndex(0));
|
|
replaceWithInPlace(result);
|
|
return result;
|
|
}
|
|
Expression result = NthRoot::Builder(childAtIndex(0), Rational::Builder(index));
|
|
replaceWithInPlace(result);
|
|
return result;
|
|
}
|
|
|
|
// Step 4: Force Float(1) in front of an orphan Power of Unit
|
|
if (parent().isUninitialized() && childAtIndex(0).type() == ExpressionNode::Type::Unit) {
|
|
Multiplication m = Multiplication::Builder(Float<double>::Builder(1.0));
|
|
replaceWithInPlace(m);
|
|
m.addChildAtIndexInPlace(*this, 1, 1);
|
|
return std::move(m);
|
|
}
|
|
|
|
return *this;
|
|
}
|
|
|
|
Expression Power::getUnit() const {
|
|
if (childAtIndex(0).type() == ExpressionNode::Type::Unit) {
|
|
return clone();
|
|
}
|
|
Expression baseUnit = childAtIndex(0).getUnit();
|
|
if (baseUnit.isUninitialized()) {
|
|
return baseUnit;
|
|
}
|
|
return Power::Builder(baseUnit, childAtIndex(1).clone());
|
|
}
|
|
|
|
// Private
|
|
|
|
// Simplification
|
|
Expression Power::denominator(ExpressionNode::ReductionContext reductionContext) const {
|
|
if (childAtIndex(0).type() == ExpressionNode::Type::Unit) {
|
|
return Expression();
|
|
}
|
|
Expression pow = clone();
|
|
// If the power is of form x^(-y), denominator should be x^y
|
|
Expression positiveIndex = pow.childAtIndex(1).makePositiveAnyNegativeNumeralFactor(reductionContext);
|
|
if (positiveIndex.isUninitialized()) {
|
|
return Expression();
|
|
}
|
|
// we cannot shallowReduce pow as it is not attached to its parent yet
|
|
// if y was -1, pow is now x^1, denominator is then only x
|
|
if (positiveIndex.isRationalOne()) {
|
|
return pow.childAtIndex(0);
|
|
}
|
|
return pow;
|
|
}
|
|
|
|
Expression Power::simplifyPowerPower(ExpressionNode::ReductionContext reductionContext) {
|
|
// this is p^e = (a^b)^e, we want a^(b*e)
|
|
Expression p = childAtIndex(0);
|
|
Multiplication m = Multiplication::Builder(p.childAtIndex(1), childAtIndex(1));
|
|
replaceChildAtIndexInPlace(0, p.childAtIndex(0));
|
|
replaceChildAtIndexInPlace(1, m);
|
|
m.shallowReduce(reductionContext);
|
|
return shallowReduce(reductionContext);
|
|
}
|
|
|
|
Expression Power::simplifyPowerMultiplication(ExpressionNode::ReductionContext reductionContext) {
|
|
// this is m^r= (a*b*c*...)^r, we want a^r * b^r *c^r * ...
|
|
Expression m = childAtIndex(0);
|
|
Expression r = childAtIndex(1);
|
|
for (int index = 0; index < m.numberOfChildren(); index++) {
|
|
Power p = Power::Builder(m.childAtIndex(index).clone(), r.clone()); // We copy r and factor to avoid inheritance issues
|
|
m.replaceChildAtIndexInPlace(index, p);
|
|
p.shallowReduce(reductionContext);
|
|
}
|
|
replaceWithInPlace(m);
|
|
return m.shallowReduce(reductionContext);
|
|
}
|
|
|
|
Expression Power::simplifyRationalRationalPower(ExpressionNode::ReductionContext reductionContext) {
|
|
// this is a^b with a, b rationals
|
|
Rational a = childAtIndex(0).convert<Rational>();
|
|
Rational b = childAtIndex(1).convert<Rational>();
|
|
if (b.isInteger()) {
|
|
Rational r = Rational::IntegerPower(a, b.signedIntegerNumerator());
|
|
if (r.numeratorOrDenominatorIsInfinity()) {
|
|
return Power::Builder(a, b);
|
|
}
|
|
replaceWithInPlace(r);
|
|
return std::move(r);
|
|
}
|
|
Expression n;
|
|
Expression d;
|
|
if (b.sign() == ExpressionNode::Sign::Negative) {
|
|
b.setSign(ExpressionNode::Sign::Positive);
|
|
n = CreateSimplifiedIntegerRationalPower(a.integerDenominator(), b, false, reductionContext);
|
|
d = CreateSimplifiedIntegerRationalPower(a.signedIntegerNumerator(), b, true, reductionContext);
|
|
} else {
|
|
n = CreateSimplifiedIntegerRationalPower(a.signedIntegerNumerator(), b, false, reductionContext);
|
|
d = CreateSimplifiedIntegerRationalPower(a.integerDenominator(), b, true, reductionContext);
|
|
}
|
|
Multiplication m = Multiplication::Builder(n, d);
|
|
replaceWithInPlace(m);
|
|
return m.shallowReduce(reductionContext);
|
|
}
|
|
|
|
Expression Power::CreateSimplifiedIntegerRationalPower(Integer i, Rational r, bool isDenominator, ExpressionNode::ReductionContext reductionContext) {
|
|
assert(!i.isZero());
|
|
assert(r.sign() == ExpressionNode::Sign::Positive);
|
|
if (i.isOne()) {
|
|
return Rational::Builder(1);
|
|
}
|
|
Integer factors[Arithmetic::k_maxNumberOfPrimeFactors];
|
|
Integer coefficients[Arithmetic::k_maxNumberOfPrimeFactors];
|
|
int numberOfPrimeFactors = Arithmetic::PrimeFactorization(i, factors, coefficients, Arithmetic::k_maxNumberOfPrimeFactors);
|
|
if (numberOfPrimeFactors < 0) {
|
|
/* We could not break i in prime factors (it might take either too many
|
|
* factors or too much time). */
|
|
Expression rClone = r.clone().setSign(isDenominator ? ExpressionNode::Sign::Negative : ExpressionNode::Sign::Positive, reductionContext);
|
|
return Power::Builder(Rational::Builder(i), rClone);
|
|
}
|
|
|
|
Integer r1(1);
|
|
Integer r2(1);
|
|
for (int index = 0; index < numberOfPrimeFactors; index++) {
|
|
Integer n = Integer::Multiplication(coefficients[index], r.signedIntegerNumerator());
|
|
IntegerDivision div = Integer::Division(n, r.integerDenominator());
|
|
r1 = Integer::Multiplication(r1, Integer::Power(factors[index], div.quotient));
|
|
r2 = Integer::Multiplication(r2, Integer::Power(factors[index], div.remainder));
|
|
}
|
|
if (r2.isOverflow() || r1.isOverflow()) {
|
|
// we overflow Integer at one point: we abort
|
|
return Power::Builder(Rational::Builder(i), r.clone());
|
|
}
|
|
Rational p1 = Rational::Builder(r2);
|
|
Integer oneExponent = isDenominator ? Integer(-1) : Integer(1);
|
|
Integer rDenominator = r.integerDenominator();
|
|
Rational p2 = Rational::Builder(oneExponent, rDenominator);
|
|
Power p = Power::Builder(p1, p2);
|
|
if (r1.isEqualTo(Integer(1)) && !i.isNegative()) {
|
|
return std::move(p);
|
|
}
|
|
Integer one(1);
|
|
Rational r3 = isDenominator ? Rational::Builder(one, r1) : Rational::Builder(r1);
|
|
Multiplication m = Multiplication::Builder();
|
|
m.addChildAtIndexInPlace(r3, 0, 0);
|
|
if (!r2.isOne()) {
|
|
m.addChildAtIndexInPlace(p, 1, 1);
|
|
}
|
|
if (i.isNegative()) {
|
|
if (reductionContext.complexFormat() == Preferences::ComplexFormat::Real) {
|
|
/* On real numbers (-1)^(p/q) =
|
|
* - 1 if p is even
|
|
* - -1 if p and q are odd
|
|
* - has no real solution otherwise */
|
|
if (!r.unsignedIntegerNumerator().isEven()) {
|
|
if (r.integerDenominator().isEven()) {
|
|
return Unreal::Builder();
|
|
} else {
|
|
m.addChildAtIndexInPlace(Rational::Builder(-1), 0, m.numberOfChildren());
|
|
}
|
|
}
|
|
} else {
|
|
/* On complex numbers, we pick the first root (-1)^(p/q) = e^(i*pi*p/q) */
|
|
r.setSign(isDenominator ? ExpressionNode::Sign::Negative : ExpressionNode::Sign::Positive);
|
|
Expression exp = CreateComplexExponent(r, reductionContext);
|
|
m.addChildAtIndexInPlace(exp, m.numberOfChildren(), m.numberOfChildren());
|
|
exp.shallowReduce(reductionContext);
|
|
}
|
|
}
|
|
return m.shallowReduce(reductionContext);
|
|
}
|
|
|
|
Expression Power::removeSquareRootsFromDenominator(ExpressionNode::ReductionContext reductionContext) {
|
|
assert(reductionContext.target() == ExpressionNode::ReductionTarget::User);
|
|
Expression result;
|
|
if (childAtIndex(0).type() == ExpressionNode::Type::Rational
|
|
&& childAtIndex(1).type() == ExpressionNode::Type::Rational
|
|
&& (childAtIndex(1).convert<Rational>().isHalf()
|
|
|| childAtIndex(1).convert<Rational>().isMinusHalf()))
|
|
{
|
|
/* We are considering a term of the form sqrt(p/q) (or 1/sqrt(p/q)), with p
|
|
* and q integers. We'll turn those into sqrt(p*q)/q (or sqrt(p*q)/p). */
|
|
Rational castedChild0 = childAtIndex(0).convert<Rational>();
|
|
Rational castedChild1 = childAtIndex(1).convert<Rational>();
|
|
Integer p = castedChild0.signedIntegerNumerator();
|
|
assert(!p.isZero()); // We eliminated case of form 0^(-1/2) at first step of shallowReduce
|
|
Integer q = castedChild0.integerDenominator();
|
|
// We do nothing for terms of the form sqrt(p)
|
|
if (!q.isOne() || castedChild1.isMinusHalf()) {
|
|
Integer pq = Integer::Multiplication(p, q);
|
|
if (pq.isOverflow()) {
|
|
return result;
|
|
}
|
|
Power sqrt = Power::Builder(Rational::Builder(pq), Rational::Builder(1, 2));
|
|
Integer one(1);
|
|
if (castedChild1.isHalf()) {
|
|
result = Multiplication::Builder(Rational::Builder(one, q), sqrt);
|
|
} else {
|
|
result = Multiplication::Builder(Rational::Builder(one, p), sqrt); // We use here the assertion that p != 0
|
|
}
|
|
sqrt.shallowReduce(reductionContext);
|
|
}
|
|
} else if (childAtIndex(1).type() == ExpressionNode::Type::Rational
|
|
&& childAtIndex(1).convert<Rational>().isMinusOne()
|
|
&& childAtIndex(0).type() == ExpressionNode::Type::Addition
|
|
&& childAtIndex(0).numberOfChildren() == 2
|
|
&& TermIsARationalSquareRootOrRational(childAtIndex(0).childAtIndex(0))
|
|
&& TermIsARationalSquareRootOrRational(childAtIndex(0).childAtIndex(1)))
|
|
{
|
|
/* We're considering a term of the form
|
|
*
|
|
* 1/(n1/d1*sqrt(p1/q1) + n2/d2*sqrt(p2/q2))
|
|
*
|
|
* and we want to turn it into
|
|
*
|
|
* n1*q2*d1*d2^2*sqrt(p1*q1) - n2*q1*d2*d1^2*sqrt(p2*q2)
|
|
* -------------------------------------------------------
|
|
* n1^2*d2^2*p1*q2 - n2^2*d1^2*p2*q1
|
|
*/
|
|
const Rational f1 = RationalFactorInExpression(childAtIndex(0).childAtIndex(0));
|
|
const Rational f2 = RationalFactorInExpression(childAtIndex(0).childAtIndex(1));
|
|
const Rational r1 = RadicandInExpression(childAtIndex(0).childAtIndex(0));
|
|
const Rational r2 = RadicandInExpression(childAtIndex(0).childAtIndex(1));
|
|
Integer n1 = f1.signedIntegerNumerator();
|
|
Integer d1 = f1.integerDenominator();
|
|
Integer p1 = r1.signedIntegerNumerator();
|
|
Integer q1 = r1.integerDenominator();
|
|
Integer n2 = f2.signedIntegerNumerator();
|
|
Integer d2 = f2.integerDenominator();
|
|
Integer p2 = r2.signedIntegerNumerator();
|
|
Integer q2 = r2.integerDenominator();
|
|
|
|
// Compute the denominator = n1^2*d2^2*p1*q2 - n2^2*d1^2*p2*q1
|
|
Integer denominator = Integer::Subtraction(
|
|
Integer::Multiplication(
|
|
Integer::Multiplication(
|
|
Integer::Power(n1, Integer(2)),
|
|
Integer::Power(d2, Integer(2))),
|
|
Integer::Multiplication(p1, q2)),
|
|
Integer::Multiplication(
|
|
Integer::Multiplication(
|
|
Integer::Power(n2, Integer(2)),
|
|
Integer::Power(d1, Integer(2))),
|
|
Integer::Multiplication(p2, q1)));
|
|
|
|
// Compute the numerator
|
|
Integer pq1 = Integer::Multiplication(p1, q1);
|
|
Integer pq2 = Integer::Multiplication(p2, q2);
|
|
Power sqrt1 = Power::Builder(Rational::Builder(pq1), Rational::Builder(1, 2));
|
|
Power sqrt2 = Power::Builder(Rational::Builder(pq2), Rational::Builder(1, 2));
|
|
Integer factor1 = Integer::Multiplication(
|
|
Integer::Multiplication(n1, d1),
|
|
Integer::Multiplication(Integer::Power(d2, Integer(2)), q2));
|
|
Multiplication m1 = Multiplication::Builder(Rational::Builder(factor1), sqrt1);
|
|
Integer factor2 = Integer::Multiplication(
|
|
Integer::Multiplication(n2, d2),
|
|
Integer::Multiplication(Integer::Power(d1, Integer(2)), q1));
|
|
Multiplication m2 = Multiplication::Builder(Rational::Builder(factor2), sqrt2);
|
|
Expression numerator;
|
|
if (denominator.isNegative()) {
|
|
numerator = Subtraction::Builder(m2, m1);
|
|
denominator.setNegative(false);
|
|
} else {
|
|
numerator = Subtraction::Builder(m1, m2);
|
|
}
|
|
if (denominator.isOverflow() || factor1.isOverflow() || factor2.isOverflow() || pq1.isOverflow() || pq2.isOverflow()) {
|
|
return result; // Escape
|
|
}
|
|
numerator = numerator.deepReduce(reductionContext);
|
|
Integer one(1);
|
|
result = Multiplication::Builder(numerator, Rational::Builder(one, denominator));
|
|
}
|
|
|
|
if (!result.isUninitialized()) {
|
|
replaceWithInPlace(result);
|
|
result = result.shallowReduce(reductionContext);
|
|
}
|
|
return result;
|
|
}
|
|
bool Power::parentIsALogarithmOfSameBase() const {
|
|
Expression p = parent();
|
|
if (p.isUninitialized()) {
|
|
return false;
|
|
}
|
|
if (p.type() == ExpressionNode::Type::Logarithm && p.childAtIndex(0) == *this) {
|
|
// parent = log(10^x)
|
|
if (p.numberOfChildren() == 1) {
|
|
if (childAtIndex(0).type() == ExpressionNode::Type::Rational && childAtIndex(0).convert<Rational>().isTen()) {
|
|
return true;
|
|
}
|
|
return false;
|
|
}
|
|
// parent = log(x^y,x)
|
|
/* Because p is a logarthim, we simplified p.childAtIndex(1) before
|
|
* p.childAtIndex(0), so p.childAtIndex(1) can be compared to childAtIndex(0). */
|
|
if (childAtIndex(0).isIdenticalTo(p.childAtIndex(1))) {
|
|
return true;
|
|
}
|
|
}
|
|
// parent = ln(e^x)
|
|
if (p.type() == ExpressionNode::Type::NaperianLogarithm
|
|
&& p.childAtIndex(0) == *this
|
|
&& childAtIndex(0).type() == ExpressionNode::Type::Constant
|
|
&& childAtIndex(0).convert<Constant>().isExponential())
|
|
{
|
|
return true;
|
|
}
|
|
return false;
|
|
}
|
|
|
|
bool Power::isNthRootOfUnity() const {
|
|
// We check we are equal to e^(i*pi) or e^(i*pi*rational)
|
|
if (childAtIndex(0).type() != ExpressionNode::Type::Constant || !childAtIndex(0).convert<Constant>().isExponential()) {
|
|
return false;
|
|
}
|
|
if (childAtIndex(1).type() != ExpressionNode::Type::Multiplication) {
|
|
return false;
|
|
}
|
|
if (childAtIndex(1).numberOfChildren() < 2 || childAtIndex(1).numberOfChildren() > 3) {
|
|
return false;
|
|
}
|
|
const Expression i = childAtIndex(1).childAtIndex(childAtIndex(1).numberOfChildren()-2);
|
|
if (i.type() != ExpressionNode::Type::Constant || !static_cast<const Constant &>(i).isIComplex()) {
|
|
return false;
|
|
}
|
|
const Expression pi = childAtIndex(1).childAtIndex(childAtIndex(1).numberOfChildren()-1);
|
|
if (pi.type() != ExpressionNode::Type::Constant || !static_cast<const Constant &>(pi).isPi()) {
|
|
return false;
|
|
}
|
|
if (numberOfChildren() == 2) {
|
|
return true;
|
|
}
|
|
if (childAtIndex(1).childAtIndex(0).type() == ExpressionNode::Type::Rational) {
|
|
return true;
|
|
}
|
|
return false;
|
|
}
|
|
|
|
Expression Power::equivalentExpressionUsingStandardExpression() const {
|
|
if (childAtIndex(1).type() == ExpressionNode::Type::Rational) {
|
|
if (childAtIndex(1).convert<Rational>().isMinusOne()) {
|
|
return Division::Builder(Rational::Builder(1), childAtIndex(0).clone());
|
|
}
|
|
if (childAtIndex(1).convert<Rational>().isHalf()) {
|
|
return SquareRoot::Builder(childAtIndex(0).clone());
|
|
}
|
|
if (childAtIndex(1).convert<Rational>().isMinusHalf()) {
|
|
return Division::Builder(Rational::Builder(1), SquareRoot::Builder(childAtIndex(0).clone()));
|
|
}
|
|
}
|
|
return Expression();
|
|
}
|
|
|
|
Expression Power::CreateComplexExponent(const Expression & r, ExpressionNode::ReductionContext reductionContext) {
|
|
// Returns e^(i*pi*r)
|
|
const Constant exp = Constant::Builder(UCodePointScriptSmallE);
|
|
Constant iComplex = Constant::Builder(UCodePointMathematicalBoldSmallI);
|
|
const Constant pi = Constant::Builder(UCodePointGreekSmallLetterPi);
|
|
Multiplication mExp = Multiplication::Builder(iComplex, pi, r.clone());
|
|
iComplex.shallowReduce(reductionContext);
|
|
Power p = Power::Builder(exp, mExp);
|
|
mExp.shallowReduce(reductionContext);
|
|
return std::move(p);
|
|
#if 0
|
|
const Constant iComplex = Constant::Builder(UCodePointMathematicalBoldSmallI);
|
|
const Constant pi = Constant::Builder(UCodePointGreekSmallLetterPi);
|
|
Expression op = Multiplication::Builder(pi, r).shallowReduce(context, complexFormat, angleUnit, false);
|
|
Cosine cos = Cosine(op).shallowReduce(context, complexFormat, angleUnit, false);;
|
|
Sine sin = Sine(op).shallowReduce(context, complexFormat, angleUnit, false);
|
|
Expression m = Multiplication::Builder(iComplex, sin);
|
|
Expression a = Addition::Builder(cos, m);
|
|
const Expression * multExpOperands[3] = {pi, r->clone()};
|
|
#endif
|
|
}
|
|
|
|
bool Power::TermIsARationalSquareRootOrRational(const Expression & e) {
|
|
if (e.type() == ExpressionNode::Type::Rational) {
|
|
return true;
|
|
}
|
|
if (e.type() == ExpressionNode::Type::Power
|
|
&& e.childAtIndex(0).type() == ExpressionNode::Type::Rational
|
|
&& e.childAtIndex(1).type() == ExpressionNode::Type::Rational
|
|
&& e.childAtIndex(1).convert<Rational>().isHalf())
|
|
{
|
|
return true;
|
|
}
|
|
if (e.type() == ExpressionNode::Type::Multiplication
|
|
&& e.numberOfChildren() == 2
|
|
&& e.childAtIndex(0).type() == ExpressionNode::Type::Rational
|
|
&& e.childAtIndex(1).type() == ExpressionNode::Type::Power
|
|
&& e.childAtIndex(1).childAtIndex(0).type() == ExpressionNode::Type::Rational
|
|
&& e.childAtIndex(1).childAtIndex(1).type() == ExpressionNode::Type::Rational
|
|
&& e.childAtIndex(1).childAtIndex(1).convert<Rational>().isHalf())
|
|
{
|
|
return true;
|
|
}
|
|
return false;
|
|
}
|
|
|
|
const Rational Power::RadicandInExpression(const Expression & e) {
|
|
if (e.type() == ExpressionNode::Type::Rational) {
|
|
return Rational::Builder(1);
|
|
} else if (e.type() == ExpressionNode::Type::Power) {
|
|
assert(e.type() == ExpressionNode::Type::Power);
|
|
assert(e.childAtIndex(0).type() == ExpressionNode::Type::Rational);
|
|
return e.childAtIndex(0).convert<Rational>();
|
|
} else {
|
|
assert(e.type() == ExpressionNode::Type::Multiplication);
|
|
assert(e.childAtIndex(1).type() == ExpressionNode::Type::Power);
|
|
assert(e.childAtIndex(1).childAtIndex(0).type() == ExpressionNode::Type::Rational);
|
|
return e.childAtIndex(1).childAtIndex(0).convert<Rational>();
|
|
}
|
|
}
|
|
|
|
const Rational Power::RationalFactorInExpression(const Expression & e) {
|
|
if (e.type() == ExpressionNode::Type::Rational) {
|
|
return static_cast<const Rational &>(e);
|
|
} else if (e.type() == ExpressionNode::Type::Power) {
|
|
return Rational::Builder(1);
|
|
} else {
|
|
assert(e.type() == ExpressionNode::Type::Multiplication);
|
|
assert(e.childAtIndex(0).type() == ExpressionNode::Type::Rational);
|
|
return e.childAtIndex(0).convert<Rational>();
|
|
}
|
|
}
|
|
|
|
bool Power::RationalExponentShouldNotBeReduced(const Rational & b, const Rational & r) {
|
|
if (r.isMinusOne()) {
|
|
return false;
|
|
}
|
|
/* We check that the simplification does not involve too complex power of
|
|
* integers (ie 3^999, 120232323232^50) that would take too much time to
|
|
* compute:
|
|
* - we cap the exponent at k_maxExactPowerMatrix
|
|
* - we cap the resulting power at DBL_MAX
|
|
* The complexity of computing a power of rational is mainly due to computing
|
|
* the GCD of the resulting numerator and denominator. Euclide algorithm's
|
|
* complexity is apportionned to the number of decimal digits in the smallest
|
|
* integer. */
|
|
Integer maxIntegerExponent = r.unsignedIntegerNumerator();
|
|
if (Integer::NaturalOrder(maxIntegerExponent, Integer(k_maxExactPowerMatrix)) > 0) {
|
|
return true;
|
|
}
|
|
|
|
double index = maxIntegerExponent.approximate<double>();
|
|
double powerNumerator = std::pow(b.unsignedIntegerNumerator().approximate<double>(), index);
|
|
double powerDenominator = std::pow(b.integerDenominator().approximate<double>(), index);
|
|
if (std::isnan(powerNumerator) || std::isnan(powerDenominator) || std::isinf(powerNumerator) || std::isinf(powerDenominator)) {
|
|
return true;
|
|
}
|
|
return false;
|
|
}
|
|
|
|
|
|
template Complex<float> PowerNode::compute<float>(std::complex<float>, std::complex<float>, Preferences::ComplexFormat);
|
|
template Complex<double> PowerNode::compute<double>(std::complex<double>, std::complex<double>, Preferences::ComplexFormat);
|
|
|
|
}
|