[poincare] In Power::approximation and SquareRoot::approximation, the

real (or imaginary) part negligence should depend on the argument value
relatively to some other values (norms of the parameters) and not of its
value absolutely!
This commit is contained in:
Émilie Feral
2020-02-25 14:39:29 +01:00
committed by LeaNumworks
parent ebdac63681
commit 515405a5df
5 changed files with 30 additions and 11 deletions

View File

@@ -10,7 +10,7 @@ namespace Poincare {
namespace ApproximationHelper {
template <typename T> int PositiveIntegerApproximationIfPossible(const ExpressionNode * expression, bool * isUndefined, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit);
template <typename T> std::complex<T> TruncateRealOrImaginaryPartAccordingToArgument(std::complex<T> c);
template <typename T> std::complex<T> TruncateRealOrImaginaryPartAccordingToArgument(std::complex<T> c, T norm1, T norm2 = 1.0);
template <typename T> using ComplexCompute = Complex<T>(*)(const std::complex<T>, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit);
template<typename T> Evaluation<T> Map(const ExpressionNode * expression, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, ComplexCompute<T> compute);

View File

@@ -15,6 +15,10 @@ template <typename T> T absMod(T a, T b) {
return result > b/2 ? b-result : result;
}
template <typename T> bool isNegligeable(T x, T precision, T norm1, T norm2) {
return x <= precision && x/norm1 <= precision && x/norm2 <= precision;
}
static inline int absInt(int x) { return x < 0 ? -x : x; }
template <typename T> int ApproximationHelper::PositiveIntegerApproximationIfPossible(const ExpressionNode * expression, bool * isUndefined, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) {
@@ -27,13 +31,15 @@ template <typename T> int ApproximationHelper::PositiveIntegerApproximationIfPos
return absInt((int)scalar);
}
template <typename T> std::complex<T> ApproximationHelper::TruncateRealOrImaginaryPartAccordingToArgument(std::complex<T> c) {
template <typename T> std::complex<T> ApproximationHelper::TruncateRealOrImaginaryPartAccordingToArgument(std::complex<T> c, T norm1, T norm2) {
T arg = std::arg(c);
T precision = 10*Expression::Epsilon<T>();
if (absMod<T>(arg, (T)M_PI) <= precision) {
T precision = 10.0*Expression::Epsilon<T>();
T argModPi = absMod<T>(arg, (T)M_PI);
T argModHalfPi = absMod<T>(arg-(T)M_PI/2.0, (T)M_PI);
if (isNegligeable(argModPi, precision, norm1, norm2)) {
c.imag(0);
}
if (absMod<T>(arg-(T)M_PI/2.0, (T)M_PI) <= precision) {
if (isNegligeable(argModHalfPi, precision, norm1, norm2)) {
c.real(0);
}
return c;
@@ -106,8 +112,8 @@ template<typename T> MatrixComplex<T> ApproximationHelper::ElementWiseOnComplexM
template int Poincare::ApproximationHelper::PositiveIntegerApproximationIfPossible<float>(Poincare::ExpressionNode const*, bool*, Poincare::Context*, Poincare::Preferences::ComplexFormat, Poincare::Preferences::AngleUnit);
template int Poincare::ApproximationHelper::PositiveIntegerApproximationIfPossible<double>(Poincare::ExpressionNode const*, bool*, Poincare::Context*, Poincare::Preferences::ComplexFormat, Poincare::Preferences::AngleUnit);
template std::complex<float> Poincare::ApproximationHelper::TruncateRealOrImaginaryPartAccordingToArgument<float>(std::complex<float>);
template std::complex<double> Poincare::ApproximationHelper::TruncateRealOrImaginaryPartAccordingToArgument<double>(std::complex<double>);
template std::complex<float> Poincare::ApproximationHelper::TruncateRealOrImaginaryPartAccordingToArgument<float>(std::complex<float>,float,float);
template std::complex<double> Poincare::ApproximationHelper::TruncateRealOrImaginaryPartAccordingToArgument<double>(std::complex<double>,double,double);
template Poincare::Evaluation<float> Poincare::ApproximationHelper::Map(const Poincare::ExpressionNode * expression, Poincare::Context * context, Poincare::Preferences::ComplexFormat, Poincare::Preferences::AngleUnit angleUnit, Poincare::ApproximationHelper::ComplexCompute<float> compute);
template Poincare::Evaluation<double> Poincare::ApproximationHelper::Map(const Poincare::ExpressionNode * expression, Poincare::Context * context, Poincare::Preferences::ComplexFormat, Poincare::Preferences::AngleUnit angleUnit, Poincare::ApproximationHelper::ComplexCompute<double> compute);
template Poincare::Evaluation<float> Poincare::ApproximationHelper::MapReduce(const Poincare::ExpressionNode * expression, Poincare::Context * context, Poincare::Preferences::ComplexFormat, Poincare::Preferences::AngleUnit angleUnit, Poincare::ApproximationHelper::ComplexAndComplexReduction<float> computeOnComplexes, Poincare::ApproximationHelper::ComplexAndMatrixReduction<float> computeOnComplexAndMatrix, Poincare::ApproximationHelper::MatrixAndComplexReduction<float> computeOnMatrixAndComplex, Poincare::ApproximationHelper::MatrixAndMatrixReduction<float> computeOnMatrices);

View File

@@ -151,8 +151,13 @@ Complex<T> PowerNode::compute(const std::complex<T> c, const std::complex<T> d,
* 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. */
return Complex<T>::Builder(ApproximationHelper::TruncateRealOrImaginaryPartAccordingToArgument(result));
* 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::TruncateRealOrImaginaryPartAccordingToArgument(result, std::log(std::abs(c)), std::abs(d)));
}
// Layout

View File

@@ -35,8 +35,12 @@ Complex<T> SquareRootNode::computeOnComplex(const std::complex<T> c, Preferences
* The error epsilon is ~1E-7 on float and ~1E-15 on double. In order to avoid
* weird results as sqrt(-1) = 6E-16+i, we compute the argument of the result
* of sqrt(c) and if arg ~ 0 [Pi], we discard the residual imaginary part and
* if arg ~ Pi/2 [Pi], we discard the residual real part.*/
return Complex<T>::Builder(ApproximationHelper::TruncateRealOrImaginaryPartAccordingToArgument(result));
* 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θ), so arg(sqrt(c)) = θ/2.
* We consider that arg[Pi] is negligeable if it is negligeable compared to
* θ. */
return Complex<T>::Builder(ApproximationHelper::TruncateRealOrImaginaryPartAccordingToArgument(result, std::arg(c)));
}
Expression SquareRootNode::shallowReduce(ReductionContext reductionContext) {

View File

@@ -149,6 +149,10 @@ QUIZ_CASE(poincare_approximation_power) {
assert_expression_approximates_to_scalar<float>("2^3", 8.0f);
assert_expression_approximates_to_scalar<double>("(3+𝐢)^(4+𝐢)", NAN);
assert_expression_approximates_to_scalar<float>("[[1,2][3,4]]^2", NAN);
assert_expression_approximates_to<float>("(-10)^0.00000001", "unreal", Radian, Real);
assert_expression_approximates_to<float>("(-10)^0.00000001", "1+3.141593ᴇ-8×𝐢", Radian, Cartesian);
}
QUIZ_CASE(poincare_approximation_subtraction) {