[poincare] Step I: add a parameter to approximation routines to indicate

if we're within a reduction routine
This commit is contained in:
Émilie Feral
2020-11-02 17:27:46 +01:00
parent 58867b66ca
commit 75dc415e27
14 changed files with 37 additions and 35 deletions

View File

@@ -260,8 +260,8 @@ public:
/* Approximation Helper */
// These methods reset the sApproximationEncounteredComplex flag. They should not be use to implement node approximation
template<typename U> static U Epsilon();
template<typename U> Expression approximate(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const;
template<typename U> U approximateToScalar(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const;
template<typename U> Expression approximate(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, bool withinReduce = false) const;
template<typename U> U approximateToScalar(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, bool withinReduce = false) const;
template<typename U> static U ApproximateToScalar(const char * text, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, Preferences::UnitFormat unitFormat, ExpressionNode::SymbolicComputation symbolicComputation = ExpressionNode::SymbolicComputation::ReplaceAllDefinedSymbolsWithDefinition);
template<typename U> U approximateWithValueForSymbol(const char * symbol, U x, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const;
/* Expression roots/extrema solver */
@@ -433,7 +433,7 @@ private:
Expression defaultUnaryFunctionDifferential() { return *this; }
/* Approximation */
template<typename U> Evaluation<U> approximateToEvaluation(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const;
template<typename U> Evaluation<U> approximateToEvaluation(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, bool withinReduce = false) const;
/* Properties */
int defaultGetPolynomialCoefficients(Context * context, const char * symbol, Expression expression[]) const;

View File

@@ -239,8 +239,8 @@ public:
typedef float SinglePrecision;
typedef double DoublePrecision;
constexpr static int k_maxNumberOfSteps = 10000;
virtual Evaluation<float> approximate(SinglePrecision p, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const = 0;
virtual Evaluation<double> approximate(DoublePrecision p, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const = 0;
virtual Evaluation<float> approximate(SinglePrecision p, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, bool withinReduce = false) const = 0;
virtual Evaluation<double> approximate(DoublePrecision p, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, bool withinReduce = false) const = 0;
/* Simplification */
/*!*/ virtual void deepReduceChildren(ReductionContext reductionContext);

View File

@@ -52,7 +52,7 @@ Expression AbsoluteValue::shallowReduce(ExpressionNode::ReductionContext reducti
}
// |x| = ±x if x is real
if (c.isReal(reductionContext.context())) {
double app = c.node()->approximate(double(), reductionContext.context(), reductionContext.complexFormat(), reductionContext.angleUnit()).toScalar();
double app = c.node()->approximate<double>(reductionContext.context(), reductionContext.complexFormat(), reductionContext.angleUnit(), true).toScalar();
if (!std::isnan(app)) {
if ((c.isNumber() && app >= 0) || app >= Expression::Epsilon<double>()) {
/* abs(a) = a with a >= 0

View File

@@ -51,7 +51,7 @@ Expression ComplexArgument::shallowReduce(ExpressionNode::ReductionContext reduc
}
bool real = c.isReal(reductionContext.context());
if (real) {
float app = c.node()->approximate(float(), reductionContext.context(), reductionContext.complexFormat(), reductionContext.angleUnit()).toScalar();
float app = c.node()->approximate(float(), reductionContext.context(), reductionContext.complexFormat(), reductionContext.angleUnit(), true).toScalar();
if (!std::isnan(app) && app >= Expression::Epsilon<float>()) {
// arg(x) = 0 if x > 0
Expression result = Rational::Builder(0);

View File

@@ -375,7 +375,7 @@ Expression Expression::defaultHandleUnitsInChildren() {
}
Expression Expression::shallowReduceUsingApproximation(ExpressionNode::ReductionContext reductionContext) {
double approx = node()->approximate(double(), reductionContext.context(), reductionContext.complexFormat(), reductionContext.angleUnit()).toScalar();
double approx = node()->approximate(double(), reductionContext.context(), reductionContext.complexFormat(), reductionContext.angleUnit(), true).toScalar();
/* If approx is capped by the largest integer such as all smaller integers can
* be exactly represented in IEEE754, approx is the exact result (no
* precision were loss). */
@@ -445,11 +445,11 @@ Expression Expression::makePositiveAnyNegativeNumeralFactor(ExpressionNode::Redu
}
template<typename U>
Evaluation<U> Expression::approximateToEvaluation(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const {
Evaluation<U> Expression::approximateToEvaluation(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, bool withinReduce) const {
sApproximationEncounteredComplex = false;
// Reset interrupting flag because some evaluation methods use it
sSimplificationHasBeenInterrupted = false;
Evaluation<U> e = node()->approximate(U(), context, complexFormat, angleUnit);
Evaluation<U> e = node()->approximate(U(), context, complexFormat, angleUnit, withinReduce);
if (complexFormat == Preferences::ComplexFormat::Real && sApproximationEncounteredComplex) {
e = Complex<U>::Undefined();
}
@@ -854,14 +854,13 @@ Expression Expression::setSign(ExpressionNode::Sign s, ExpressionNode::Reduction
/* Evaluation */
template<typename U>
Expression Expression::approximate(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const {
return isUninitialized() ? Undefined::Builder() : approximateToEvaluation<U>(context, complexFormat, angleUnit).complexToExpression(complexFormat);
Expression Expression::approximate(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, bool withinReduce) const {
return isUninitialized() ? Undefined::Builder() : approximateToEvaluation<U>(context, complexFormat, angleUnit, withinReduce).complexToExpression(complexFormat);
}
template<typename U>
U Expression::approximateToScalar(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const {
return approximateToEvaluation<U>(context, complexFormat, angleUnit).toScalar();
U Expression::approximateToScalar(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, bool withinReduce) const {
return approximateToEvaluation<U>(context, complexFormat, angleUnit, withinReduce).toScalar();
}
template<typename U>
@@ -1155,17 +1154,17 @@ void Expression::bracketRoot(const char * symbol, double start, double step, dou
template float Expression::Epsilon<float>();
template double Expression::Epsilon<double>();
template Expression Expression::approximate<float>(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const;
template Expression Expression::approximate<double>(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const;
template Expression Expression::approximate<float>(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, bool withinReduce) const;
template Expression Expression::approximate<double>(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, bool withinReduce) const;
template float Expression::approximateToScalar(Context * context, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit) const;
template double Expression::approximateToScalar(Context * context, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit) const;
template float Expression::approximateToScalar(Context * context, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit, bool withinReduce) const;
template double Expression::approximateToScalar(Context * context, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit, bool withinReduce) const;
template float Expression::ApproximateToScalar<float>(const char * text, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, Preferences::UnitFormat unitFormat, ExpressionNode::SymbolicComputation symbolicComputation);
template double Expression::ApproximateToScalar<double>(const char * text, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, Preferences::UnitFormat unitFormat, ExpressionNode::SymbolicComputation symbolicComputation);
template float Expression::ApproximateToScalar<float>(const char * text, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, Preferences::UnitFormat unitFormat, ExpressionNode::SymbolicComputation symbolicComputation, bool withinReduce);
template double Expression::ApproximateToScalar<double>(const char * text, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, Preferences::UnitFormat unitFormat, ExpressionNode::SymbolicComputation symbolicComputation, bool withinReduce);
template Evaluation<float> Expression::approximateToEvaluation(Context * context, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit) const;
template Evaluation<double> Expression::approximateToEvaluation(Context * context, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit) const;
template Evaluation<float> Expression::approximateToEvaluation(Context * context, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit, bool withinReduce) const;
template Evaluation<double> Expression::approximateToEvaluation(Context * context, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit, bool withinReduce) const;
template float Expression::approximateWithValueForSymbol(const char * symbol, float x, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const;
template double Expression::approximateWithValueForSymbol(const char * symbol, double x, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const;

View File

@@ -45,7 +45,8 @@ Expression HyperbolicTrigonometricFunction::shallowReduce(ExpressionNode::Reduct
&& e.approximateToScalar<double>(
reductionContext.context(),
reductionContext.complexFormat(),
reductionContext.angleUnit()) >= 1.0))
reductionContext.angleUnit(),
true) >= 1.0))
{
result = e;
}
@@ -58,7 +59,8 @@ Expression HyperbolicTrigonometricFunction::shallowReduce(ExpressionNode::Reduct
&& e.approximateToScalar<double>(
reductionContext.context(),
reductionContext.complexFormat(),
reductionContext.angleUnit()) >= 0.0))
reductionContext.angleUnit(),
true) >= 0.0))
{
result = e;
}
@@ -79,7 +81,8 @@ Expression HyperbolicTrigonometricFunction::shallowReduce(ExpressionNode::Reduct
&& std::fabs(e.approximateToScalar<double>(
reductionContext.context(),
reductionContext.complexFormat(),
reductionContext.angleUnit())) < 1.0))
reductionContext.angleUnit()),
true) < 1.0))
{
result = e;
}

View File

@@ -300,7 +300,7 @@ Expression Logarithm::simpleShallowReduce(Context * context, Preferences::Comple
}
bool isNegative = true;
Expression result;
Evaluation<float> baseApproximation = b.node()->approximate(1.0f, context, complexFormat, angleUnit);
Evaluation<float> baseApproximation = b.node()->approximate(1.0f, context, complexFormat, angleUnit, true);
std::complex<float> logDenominator = std::log10(static_cast<Complex<float>&>(baseApproximation).stdComplex());
if (logDenominator.imag() != 0.0f || logDenominator.real() == 0.0f) {
result = Undefined::Builder();

View File

@@ -224,7 +224,7 @@ Matrix Matrix::rowCanonize(ExpressionNode::ReductionContext reductionContext, Ex
float bestPivot = 0.0;
while (iPivot_temp < m) {
// Using float to find the biggest pivot is sufficient.
float pivot = AbsoluteValue::Builder(matrixChild(iPivot_temp, k).clone()).approximateToScalar<float>(reductionContext.context(), reductionContext.complexFormat(), reductionContext.angleUnit());
float pivot = AbsoluteValue::Builder(matrixChild(iPivot_temp, k).clone()).approximateToScalar<float>(reductionContext.context(), reductionContext.complexFormat(), reductionContext.angleUnit(), true);
// Handle very low pivots
if (pivot == 0.0f && matrixChild(iPivot_temp, k).nullStatus(reductionContext.context()) != ExpressionNode::NullStatus::Null) {
pivot = FLT_MIN;

View File

@@ -89,7 +89,7 @@ Expression NAryExpression::checkChildrenAreRationalIntegersAndUpdate(ExpressionN
return replaceWithUndefinedInPlace();
}
// If c was complex but with a null imaginary part, real part is checked.
float app = c.approximateToScalar<float>(reductionContext.context(), reductionContext.complexFormat(), reductionContext.angleUnit());
float app = c.approximateToScalar<float>(reductionContext.context(), reductionContext.complexFormat(), reductionContext.angleUnit(), true);
if (std::isfinite(app) && app != std::round(app)) {
return replaceWithUndefinedInPlace();
}

View File

@@ -779,7 +779,7 @@ Expression Power::shallowReduce(ExpressionNode::ReductionContext reductionContex
* - (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());
Expression approximation = powerBase.approximate<float>(reductionContext.context(), reductionContext.complexFormat(), reductionContext.angleUnit(), true);
if (approximation.type() == ExpressionNode::Type::Unreal) {
// The inner power is unreal, return "unreal"
replaceWithInPlace(approximation);

View File

@@ -67,7 +67,7 @@ Expression SignFunction::shallowReduce(ExpressionNode::ReductionContext reductio
if (s == ExpressionNode::Sign::Negative) {
resultSign = Rational::Builder(-1);
} else {
Evaluation<float> childApproximated = child.node()->approximate(1.0f, reductionContext.context(), reductionContext.complexFormat(), reductionContext.angleUnit());
Evaluation<float> childApproximated = child.node()->approximate(1.0f, reductionContext.context(), reductionContext.complexFormat(), reductionContext.angleUnit(), true);
assert(childApproximated.type() == EvaluationNode<float>::Type::Complex);
Complex<float> c = static_cast<Complex<float>&>(childApproximated);
if (std::isnan(c.imag()) || std::isnan(c.real()) || c.imag() != 0) {

View File

@@ -286,7 +286,7 @@ Expression Trigonometry::shallowReduceInverseFunction(Expression & e, Expression
// Step 1. Look for an expression of type "acos(cos(x))", return x
if (AreInverseFunctions(e.childAtIndex(0), e)) {
float x = e.childAtIndex(0).childAtIndex(0).node()->approximate(float(), reductionContext.context(), reductionContext.complexFormat(), angleUnit).toScalar();
float x = e.childAtIndex(0).childAtIndex(0).nodex()->approximate(float(), reductionContext.context(), reductionContext.complexFormat(), angleUnit, true).toScalar();
if (!(std::isinf(x) || std::isnan(x))) {
Expression result = e.childAtIndex(0).childAtIndex(0);
// We translate the result within [-π,π] for acos(cos), [-π/2,π/2] for asin(sin) and atan(tan)
@@ -314,7 +314,7 @@ Expression Trigonometry::shallowReduceInverseFunction(Expression & e, Expression
// Step 2. Special case for atan(sin(x)/cos(x))
if (e.type() == ExpressionNode::Type::ArcTangent && ExpressionIsEquivalentToTangent(e.childAtIndex(0))) {
float trigoOp = e.childAtIndex(0).childAtIndex(1).childAtIndex(0).node()->approximate(float(), reductionContext.context(), reductionContext.complexFormat(), angleUnit).toScalar();
float trigoOp = e.childAtIndex(0).childAtIndex(1).childAtIndex(0).node()->approximate(float(), reductionContext.context(), reductionContext.complexFormat(), angleUnit, true).toScalar();
if (trigoOp >= -pi/2.0f && trigoOp <= pi/2.0f) {
Expression result = e.childAtIndex(0).childAtIndex(1).childAtIndex(0);
e.replaceWithInPlace(result);

View File

@@ -59,7 +59,7 @@ Expression TrigonometryCheatTable::simplify(const Expression e, ExpressionNode::
}
// Approximate e to quickly compare it to cheat table entries
float eValue = e.node()->approximate(float(), reductionContext.context(), reductionContext.complexFormat(), reductionContext.angleUnit()).toScalar();
float eValue = e.node()->approximation(float(), reductionContext.context(), reductionContext.complexFormat(), reductionContext.angleUnit(), true).toScalar();
if (std::isnan(eValue) || std::isinf(eValue)) {
return Expression();
}

View File

@@ -33,7 +33,7 @@ const Expression VariableContext::expressionForSymbolAbstract(const SymbolAbstra
Symbol unknownSymbol = Symbol::Builder(UCodePointUnknown);
if (m_name != nullptr && strcmp(m_name, unknownSymbol.name()) == 0) {
assert(std::isnan(unknownSymbolValue));
unknownSymbolValue = m_value.approximateToScalar<float>(this, Preferences::sharedPreferences()->complexFormat(),Preferences::sharedPreferences()->angleUnit());
unknownSymbolValue = m_value.approximateToScalar<float>(this, Preferences::sharedPreferences()->complexFormat(), Preferences::sharedPreferences()->angleUnit(), true);
}
return ContextWithParent::expressionForSymbolAbstract(symbol, clone, unknownSymbolValue);
}