Merge remote-tracking branch 'upstream/master' into omega-hotfix

2026-01-18 16:27:34 +01:00 · 2020-03-17 20:59:23 +01:00
parent eadcf0b6f9 5f8f1d96ee
commit bfe2fc1a04
32 changed files with 496 additions and 404 deletions
--- a/apps/probability/distribution/distribution.cpp
+++ b/apps/probability/distribution/distribution.cpp
@@ -25,9 +25,12 @@ double Distribution::rightIntegralFromAbscissa(double x) const {
 }

 double Distribution::finiteIntegralBetweenAbscissas(double a, double b) const {
-  if (b <= a) {
+  if (b < a) {
    return 0.0;
  }
+  if (a == b) {
+    return evaluateAtDiscreteAbscissa(a);
+  }
  if (isContinuous()) {
    return cumulativeDistributiveFunctionAtAbscissa(b) - cumulativeDistributiveFunctionAtAbscissa(a);
  }
@@ -100,6 +103,7 @@ double Distribution::rightIntegralInverseForProbability(double * probability) {
 }

 double Distribution::evaluateAtDiscreteAbscissa(int k) const {
+  assert(isContinuous()); // Discrete distributions override this method
  return 0.0;
 }

--- a/apps/probability/test/distributions.cpp
+++ b/apps/probability/test/distributions.cpp
@@ -19,7 +19,13 @@ void assert_cumulative_distributive_function_direct_and_inverse_is(Probability::
  quiz_assert(!std::isnan(r));
  quiz_assert(!std::isinf(r));
  quiz_assert(std::fabs(r-x) < FLT_EPSILON || std::fabs(r-x)/x < FLT_EPSILON);
+}

+void assert_finite_integral_between_abscissas_is(Probability::Distribution * distribution, double a, double b, double result) {
+  double r = distribution->finiteIntegralBetweenAbscissas(a, b);
+  quiz_assert(!std::isnan(r));
+  quiz_assert(!std::isinf(r));
+  quiz_assert(std::fabs(r-result) < FLT_EPSILON || std::fabs(r-result)/result < FLT_EPSILON);
 }

 //TODO other distributions
@@ -37,6 +43,13 @@ QUIZ_CASE(binomial_distribution) {
  distribution.setParameterAtIndex(0.1, 1);
  assert_cumulative_distributive_function_direct_and_inverse_is(&distribution, 0.0, 0.166771816996665822596668249389040283858776092529296875);
  assert_cumulative_distributive_function_direct_and_inverse_is(&distribution, 1.0, 0.4817852491014791294077213024138472974300384521484375);
+
+  // B(21, 0.2)
+  distribution.setParameterAtIndex(21.0, 0);
+  distribution.setParameterAtIndex(0.2, 1);
+  assert_finite_integral_between_abscissas_is(&distribution, 4.0, 4.0, 0.21563235015849934848);
+  assert_finite_integral_between_abscissas_is(&distribution, 5.0, 4.0, 0.0);
+  assert_finite_integral_between_abscissas_is(&distribution, 4.0, 5.0, 0.398919847793223794688);
 }

 QUIZ_CASE(chi_squared_distribution) {
@@ -58,6 +71,12 @@ QUIZ_CASE(chi_squared_distribution) {
  assert_cumulative_distributive_function_direct_and_inverse_is(&distribution, 1.3, 0.047059684573231390369851823152202996425330638885498046875);
  assert_cumulative_distributive_function_direct_and_inverse_is(&distribution, 2.9874567, 0.250530060451470470983537097708904184401035308837890625);
  assert_cumulative_distributive_function_direct_and_inverse_is(&distribution, 4.987, 0.53051693435084168459781039928202517330646514892578125);
+
+  // Chi Squared distribution with 6 degrees of freedom
+  distribution.setParameterAtIndex(6.0, 0);
+  assert_finite_integral_between_abscissas_is(&distribution, 4.0, 4.0, 0.0);
+  assert_finite_integral_between_abscissas_is(&distribution, 5.0, 4.0, 0.0);
+  assert_finite_integral_between_abscissas_is(&distribution, 4.0, 5.0, 0.1328633002997339414718);
 }

 QUIZ_CASE(student_distribution) {
@@ -79,6 +98,12 @@ QUIZ_CASE(student_distribution) {
  assert_cumulative_distributive_function_direct_and_inverse_is(&distribution, -4.987, 0.00167496657737900025118837898929768925881944596767425537109375);
  assert_cumulative_distributive_function_direct_and_inverse_is(&distribution, 1.3, 0.876837383157582639370275501278229057788848876953125);
  assert_cumulative_distributive_function_direct_and_inverse_is(&distribution, 2.9874567, 0.98612148076325445433809591122553683817386627197265625);
+
+  // Student distribution with 6 degrees of freedom
+  distribution.setParameterAtIndex(6.0, 0);
+  assert_finite_integral_between_abscissas_is(&distribution, 4.0, 4.0, 0.0);
+  assert_finite_integral_between_abscissas_is(&distribution, 5.0, 4.0, 0.0);
+  assert_finite_integral_between_abscissas_is(&distribution, 4.0, 5.0, 0.002333318101494775250);
 }

 QUIZ_CASE(geometric_distribution) {
@@ -92,6 +117,12 @@ QUIZ_CASE(geometric_distribution) {
  distribution.setParameterAtIndex(0.2, 0);
  assert_cumulative_distributive_function_direct_and_inverse_is(&distribution, 7.0, 0.8322278399999998299563230830244719982147216796875);
  assert_cumulative_distributive_function_direct_and_inverse_is(&distribution, 3.0, 0.5904);
+
+  // Geometric distribution with probability of success 0.4
+  distribution.setParameterAtIndex(0.4, 0);
+  assert_finite_integral_between_abscissas_is(&distribution, 1.0, 1.0, 0.24);
+  assert_finite_integral_between_abscissas_is(&distribution, 2.0, 1.0, 0.0);
+  assert_finite_integral_between_abscissas_is(&distribution, 1.0, 2.0, 0.384);
 }

 QUIZ_CASE(fisher_distribution) {
@@ -115,4 +146,10 @@ QUIZ_CASE(fisher_distribution) {
  assert_cumulative_distributive_function_direct_and_inverse_is(&distribution, 1.4, 0.94560850441205857);
  assert_cumulative_distributive_function_direct_and_inverse_is(&distribution, 1.425, 0.95425004959692871775);

+  // Fisher distribution with d1 = 4 and d2 = 2
+  distribution.setParameterAtIndex(4.0, 0);
+  distribution.setParameterAtIndex(2.0, 1);
+  assert_finite_integral_between_abscissas_is(&distribution, 1.0, 1.0, 0.0);
+  assert_finite_integral_between_abscissas_is(&distribution, 2.0, 1.0, 0.0);
+  assert_finite_integral_between_abscissas_is(&distribution, 1.0, 2.0, 0.19555555555555555);
 }
--- a/build/targets.all.mak
+++ b/build/targets.all.mak
@@ -29,32 +29,34 @@ all_official:
 	$(call file_check,$(ANDROID_GRADLE_PROPERTIES))
 	$(call file_check,$(IOS_MOBILE_PROVISION))
 	$(call command_check,$(EMCC))
-	$(Q) rm -rf output/stable_release
-	$(Q) mkdir -p output/stable_release
+	$(Q) rm -rf output/all_official
+	$(Q) mkdir -p output/all_official
 	$(Q) echo "BUILD_FIRMWARE    DEVICE N0110"
 	$(Q) $(MAKE) clean
 	$(Q) $(MAKE) epsilon.official.onboarding.dfu
-	$(Q) cp output/release/device/n0110/epsilon.official.onboarding.dfu output/stable_release/epsilon.device.n0110.dfu
+	$(Q) cp output/release/device/n0110/epsilon.official.onboarding.dfu output/all_official/epsilon.device.n0110.dfu
 	$(Q) echo "BUILD_FIRMWARE    DEVICE N0100"
 	$(Q) $(MAKE) MODEL=n0100 clean
 	$(Q) $(MAKE) MODEL=n0100 epsilon.official.onboarding.dfu
-	$(Q) cp output/release/device/n0100/epsilon.official.onboarding.dfu output/stable_release/epsilon.device.n0100.dfu
+	$(Q) cp output/release/device/n0100/epsilon.official.onboarding.dfu output/all_official/epsilon.device.n0100.dfu
 	$(Q) echo "BUILD_FIRMWARE    SIMULATOR WEB ZIP"
 	$(Q) $(MAKE) DEBUG=0 PLATFORM=simulator TARGET=web clean
 	$(Q) $(call source_emsdk); $(MAKE) DEBUG=0 PLATFORM=simulator TARGET=web output/release/simulator/web/simulator.official.zip
-	$(Q) cp output/release/simulator/web/simulator.official.zip output/stable_release/simulator.web.zip
+	$(Q) cp output/release/simulator/web/simulator.official.zip output/all_official/simulator.web.zip
 	$(Q) echo "BUILD_FIRMWARE    SIMULATOR WEB JS"
 	$(Q) $(call source_emsdk); $(MAKE) DEBUG=0 PLATFORM=simulator TARGET=web epsilon.official.js
-	$(Q) cp output/release/simulator/web/epsilon.official.js output/stable_release/epsilon.js
+	$(Q) cp output/release/simulator/web/epsilon.official.js output/all_official/epsilon.js
+	$(Q) cp output/release/simulator/web/epsilon.official.js.mem output/all_official/epsilon.js.mem
 	$(Q) echo "BUILD_FIRMWARE    SIMULATOR WEB PYTHON JS"
 	$(Q) $(MAKE) DEBUG=0 PLATFORM=simulator TARGET=web clean
 	$(Q) $(call source_emsdk); $(MAKE) DEBUG=0 PLATFORM=simulator TARGET=web EPSILON_GETOPT=1 EPSILON_APPS=code epsilon.official.js
-	$(Q) cp output/release/simulator/web/epsilon.official.js output/stable_release/epsilon.python.js
+	$(Q) cp output/release/simulator/web/epsilon.official.js output/all_official/epsilon.python.js
+	$(Q) cp output/release/simulator/web/epsilon.official.js.mem output/all_official/epsilon.python.js.mem
 	$(Q) echo "BUILD_FIRMWARE    SIMULATOR ANDROID"
 	$(Q) $(MAKE) PLATFORM=simulator TARGET=android clean
 	$(Q) $(MAKE) PLATFORM=simulator TARGET=android epsilon.official.apk
-	$(Q) cp output/release/simulator/android/app/outputs/apk/release/android-release-unsigned.apk output/stable_release/epsilon.apk
+	$(Q) cp output/release/simulator/android/app/outputs/apk/release/android-release-unsigned.apk output/all_official/epsilon.apk
 	$(Q) echo "BUILD_FIRMWARE    SIMULATOR IOS"
 	$(Q) $(MAKE) PLATFORM=simulator TARGET=ios clean
-	$(Q) $(MAKE) PLATFORM=simulator TARGET=ios IOS_PROVISIONNING_PROFILE="~/Downloads/NumWorks_Graphing_Calculator_Distribution.mobileprovision" output/release/simulator/ios/app/epsilon.official.ipa
-	$(Q) cp output/release/simulator/ios/app/epsilon.official.ipa output/stable_release/epsilon.ipa
+	$(Q) $(MAKE) PLATFORM=simulator TARGET=ios IOS_PROVISIONNING_PROFILE=$(IOS_MOBILE_PROVISION) output/release/simulator/ios/app/epsilon.official.ipa
+	$(Q) cp output/release/simulator/ios/app/epsilon.official.ipa output/all_official/epsilon.ipa
--- a/liba/include/math.h
+++ b/liba/include/math.h
@@ -35,13 +35,13 @@ typedef double double_t;
 #define FP_SUBNORMAL 0x08
 #define FP_ZERO 0x10

-#define fpclassify(x) __builtin_fpclassify(FP_NAN, FP_INFINITE, FP_NORMAL, FP_SUBNORMAL, FP_ZERO, x)
-#define signbit(x) __builtin_signbit(x)
 #define finite(x) __builtin_finite(x)
+#define fpclassify(x) __builtin_fpclassify(FP_NAN, FP_INFINITE, FP_NORMAL, FP_SUBNORMAL, FP_ZERO, x)
 #define isfinite(x) __builtin_isfinite(x)
 #define isnormal(x) __builtin_isnormal(x)
 #define isnan(x) __builtin_isnan(x)
 #define isinf(x) __builtin_isinf(x)
+#define signbit(x) __builtin_signbit(x)

 float acosf(float x);
 float acoshf(float x);
@@ -60,19 +60,19 @@ float fabsf(float x);
 float fmaxf(float x, float y);
 float floorf(float x);
 float fmodf(float x, float y);
-float frexpf(float x, int *eptr);
+float frexpf(float x, int *exp);
 float hypotf(float x, float y);
-float ldexpf(float x, int n);
+float ldexpf(float x, int exp);
 float lgammaf(float x);
 float lgammaf_r(float x, int *signgamp);
 float log1pf(float x);
 float log10f(float x);
 float logf(float x);
-float modff(float value, float *iptr);
+float modff(float x, float *iptr);
 float nearbyintf(float x);
 float powf(float x, float y);
 float roundf(float x);
-float scalbnf(float x, int n);
+float scalbnf(float x, int exp);
 float sinf(float x);
 float sinhf(float x);
 float sqrtf(float x);
@@ -99,7 +99,7 @@ double fabs(double x);
 double fmax(double x, double y);
 double floor(double x);
 double fmod(double x, double y);
-double frexp(double x, int *eptr);
+double frexp(double x, int *exp);
 double hypot(double x, double y);
 double lgamma(double x);
 double lgamma_r(double x, int *signgamp);
@@ -108,14 +108,13 @@ double log1p(double x);
 double log10(double x);
 double log2(double x);
 double logb(double x);
-double modf(double value, double *iptr);
+double modf(double x, double *iptr);
 double nearbyint(double x);
 double pow(double x, double y);
-double nearbyint(double x);
 double rint(double x);
 double round(double x);
-double scalb(double x, double fn);
-double scalbn(double x, int n);
+double scalb(double x, double exp);
+double scalbn(double x, int exp);
 double sin(double x);
 double sinh(double x);
 double sqrt(double x);
@@ -124,6 +123,8 @@ double tanh(double x);
 double tgamma(double x);
 double trunc(double x);

+extern int signgam;
+
 /* The C99 standard says that any libc function can be re-declared as a macro.
 * (See N1124 paragraph 7.1.4). This means that C files willing to actually
 * implement said functions should either re-define the prototype or #undef the
@@ -145,18 +146,19 @@ double trunc(double x);
 #define fabsf(x) __builtin_fabsf(x)
 #define floorf(x) __builtin_floorf(x)
 #define fmodf(x, y) __builtin_fmodf(x, y)
-#define frexpf(x, y) __builtin_frexpf(x, y)
-#define ldexpf(x, n) __builtin_ldexpf(x, n)
+#define frexpf(x, exp) __builtin_frexpf(x, exp)
+#define ldexpf(x, exp) __builtin_ldexpf(x, exp)
 #define lgammaf(x) __builtin_lgammaf(x)
 #define lgammaf_r(x, signgamp) __builtin_lgammaf_r(x, signgamp)
 #define log1pf(x) __builtin_log1pf(x)
 #define log10f(x) __builtin_log10f(x)
 #define logf(x) __builtin_logf(x)
-#define nanf(s) __builtin_nanf(s)
+#define modff(x, iptr) __builtin_modff(x, iptr)
+#define nanf(tagp) __builtin_nanf(tagp)
 #define nearbyintf(x) __builtin_nearbyintf(x)
 #define powf(x, y) __builtin_powf(x, y)
 #define roundf(x) __builtin_roundf(x)
-#define scalbnf(x, n) __builtin_scalbnf(x, n)
+#define scalbnf(x, exp) __builtin_scalbnf(x, exp)
 #define sinf(x) __builtin_sinf(x)
 #define sinhf(x) __builtin_sinhf(x)
 #define sqrtf(x) __builtin_sqrtf(x)
@@ -182,7 +184,8 @@ double trunc(double x);
 #define fabs(x) __builtin_fabs(x)
 #define floor(x) __builtin_floor(x)
 #define fmod(x, y) __builtin_fmod(x, y)
-#define ldexp(x, n) __builtin_scalbn(x, n)
+#define frexp(x, exp) __builtin_frexp(x, exp)
+#define ldexp(x, exp) __builtin_scalbn(x, exp)
 #define lgamma(x) __builtin_lgamma(x)
 #define lgamma_r(x, signgamp) __builtin_lgamma_r(x, signgamp)
 #define log(x) __builtin_log(x)
@@ -190,11 +193,14 @@ double trunc(double x);
 #define log10(x) __builtin_log10(x)
 #define log2(x) __builtin_log2(x)
 #define logb(x) __builtin_logb(x)
-#define nan(s) __builtin_nan(s)
+#define modf(x, iptr) __builtin_modf(x, iptr)
+#define nan(tagp) __builtin_nan(tagp)
+#define nearbyint(x) __builtin_nearbyint(x)
 #define pow(x, y) __builtin_pow(x, y)
 #define rint(x) __builtin_rint(x)
 #define round(x) __builtin_round(x)
-#define scalbn(x, n) __builtin_scalbn(x, n)
+#define scalb(x, exp) __builtin_scalb(x, exp)
+#define scalbn(x, exp) __builtin_scalbn(x, exp)
 #define sin(x) __builtin_sin(x)
 #define sinh(x) __builtin_sinh(x)
 #define sqrt(x) __builtin_sqrt(x)
@@ -203,8 +209,6 @@ double trunc(double x);
 #define tgamma(x) __builtin_tgamma(x)
 #define trunc(x) __builtin_trunc(x)

-extern int signgam;
-
 LIBA_END_DECLS

 #endif
--- a/liba/test/math.c
+++ b/liba/test/math.c
@@ -1,210 +1,97 @@
 // This file tests that each math fuction links.
 #include <math.h>

-int test_fpclassifyf(float x) {
-  return fpclassify(x);
-}
-int test_isfinitef(float x) {
-  return isfinite(x);
-}
-int test_isnormalf(float x) {
-  return isnormal(x);
-}
-int test_isnanf(float x) {
-  return isnan(x);
-}
-int test_isinff(float x) {
-  return isinf(x);
-}
+int test_fpclassifyf(float x) { return fpclassify(x); }
+int test_signbitf(float x) { return signbit(x); }
+int test_finitef(float x) { return finite(x); }
+int test_isfinitef(float x) { return isfinite(x); }
+int test_isnormalf(float x) { return isnormal(x); }
+int test_isnanf(float x) { return isnan(x); }
+int test_isinff(float x) { return isinf(x); }

-float test_acosf(float x) {
-  return acosf(x);
-}
-float test_acoshf(float x) {
-  return acoshf(x);
-}
-float test_asinf(float x) {
-  return asinf(x);
-}
-float test_asinhf(float x) {
-  return asinhf(x);
-}
-float test_atanf(float x) {
-  return atanf(x);
-}
-float test_atan2f(float y, float x) {
-  return atan2f(y, x);
-}
-float test_atanhf(float x) {
-  return atanhf(x);
-}
-float test_ceilf(float x) {
-  return ceilf(x);
-}
-float test_copysignf(float x, float y) {
-  return copysignf(x, y);
-}
-float test_cosf(float x) {
-  return cosf(x);
-}
-float test_coshf(float x) {
-  return coshf(x);
-}
-float test_expf(float x) {
-  return expf(x);
-}
-float test_expm1f(float x) {
-  return expm1f(x);
-}
-float test_fabsf(float x) {
-  return fabsf(x);
-}
-float test_floorf(float x) {
-  return floorf(x);
-}
-float test_fmodf(float x, float y) {
-  return fmodf(x, y);
-}
-float test_lgammaf(float x) {
-  return lgammaf(x);
-}
-float test_lgammaf_r(float x, int *signgamp) {
-  return lgammaf_r(x, signgamp);
-}
-float test_log1pf(float x) {
-  return log1pf(x);
-}
-float test_log10f(float x) {
-  return log10f(x);
-}
-float test_logf(float x) {
-  return logf(x);
-}
-float test_nanf(const char *s) {
-  return nanf(s);
-}
-float test_nearbyintf(float x) {
-  return nearbyintf(x);
-}
-float test_powf(float x, float y) {
-  return powf(x, y);
-}
-float test_roundf(float x) {
-  return roundf(x);
-}
-float test_scalbnf(float x, int n) {
-  return scalbnf(x, n);
-}
-float test_sinf(float x) {
-  return sinf(x);
-}
-float test_sinhf(float x) {
-  return sinhf(x);
-}
-float test_sqrtf(float x) {
-  return sqrtf(x);
-}
-float test_tanf(float x) {
-  return tanf(x);
-}
-float test_tanhf(float x) {
-  return tanhf(x);
-}
+float test_acosf(float x) { return acosf(x); }
+float test_acoshf(float x) { return acoshf(x); }
+float test_asinf(float x) { return asinf(x); }
+float test_asinhf(float x) { return asinhf(x); }
+float test_atanf(float x) { return atanf(x); }
+float test_atan2f(float y, float x) { return atan2f(y, x); }
+float test_atanhf(float x) { return atanhf(x); }
+float test_ceilf(float x) { return ceilf(x); }
+float test_copysignf(float x, float y) { return copysignf(x, y); }
+float test_cosf(float x) { return cosf(x); }
+float test_coshf(float x) { return coshf(x); }
+float test_expf(float x) { return expf(x); }
+float test_expm1f(float x) { return expm1f(x); }
+float test_fabsf(float x) { return fabsf(x); }
+float test_floorf(float x) { return floorf(x); }
+float test_fmodf(float x, float y) { return fmodf(x, y); }
+float test_frexpf(float x, int *exp) { return frexpf(x, exp); }
+float test_ldexpf(float x, int exp) { return ldexpf(x, exp); }
+float test_lgammaf(float x) { return lgammaf(x); }
+float test_lgammaf_r(float x, int *signgamp) { return lgammaf_r(x, signgamp); }
+float test_log1pf(float x) { return log1pf(x); }
+float test_log10f(float x) { return log10f(x); }
+float test_logf(float x) { return logf(x); }
+float test_modff(float x, float *iptr) { return modff(x, iptr); }
+float test_nanf(const char *s) { return nanf(s); }
+float test_nearbyintf(float x) { return nearbyintf(x); }
+float test_powf(float x, float y) { return powf(x, y); }
+float test_roundf(float x) { return roundf(x); }
+float test_scalbnf(float x, int exp) { return scalbnf(x, exp); }
+float test_sinf(float x) { return sinf(x); }
+float test_sinhf(float x) { return sinhf(x); }
+float test_sqrtf(float x) { return sqrtf(x); }
+float test_tanf(float x) { return tanf(x); }
+float test_tanhf(float x) { return tanhf(x); }
+float test_truncf(float x) { return truncf(x); }

-int test_fpclassify(double x) {
-  return fpclassify(x);
-}
-int test_isfinite(double x) {
-  return isfinite(x);
-}
-int test_isnormal(double x) {
-  return isnormal(x);
-}
-int test_isnan(double x) {
-  return isnan(x);
-}
-int test_isinf(double x) {
-  return isinf(x);
-}
+int test_fpclassify(double x) { return fpclassify(x); }
+int test_signbit(double x) { return signbit(x); }
+int test_finite(double x) { return finite(x); }
+int test_isfinite(double x) { return isfinite(x); }
+int test_isnormal(double x) { return isnormal(x); }
+int test_isnan(double x) { return isnan(x); }
+int test_isinf(double x) { return isinf(x); }

-double test_acos(double x) {
-  return acos(x);
-}
-double test_acosh(double x) {
-  return acosh(x);
-}
-double test_asin(double x) {
-  return asin(x);
-}
-double test_asinh(double x) {
-  return asinh(x);
-}
-double test_atan(double x) {
-  return atan(x);
-}
-double test_atanh(double x) {
-  return atanh(x);
-}
-double test_ceil(double x) {
-  return ceil(x);
-}
-double test_copysign(double x, double y) {
-  return copysign(x, y);
-}
-double test_cos(double x) {
-  return cos(x);
-}
-double test_cosh(double x) {
-  return cosh(x);
-}
-double test_exp(double x) {
-  return exp(x);
-}
-double test_expm1(double x) {
-  return expm1(x);
-}
-double test_fabs(double x) {
-  return fabs(x);
-}
-double test_floor(double x) {
-  return floor(x);
-}
-double test_lgamma(double x) {
-  return lgamma(x);
-}
-double test_lgamma_r(double x, int *signgamp) {
-  return lgamma_r(x, signgamp);
-}
-double test_log1p(double x) {
-  return log1p(x);
-}
-double test_log10(double x) {
-  return log10(x);
-}
-double test_log(double x) {
-  return log(x);
-}
-double test_pow(double x, double y) {
-  return pow(x, y);
-}
-double test_round(double x) {
-  return round(x);
-}
-double test_scalbn(double x, int n) {
-  return scalbn(x, n);
-}
-double test_sin(double x) {
-  return sin(x);
-}
-double test_sinh(double x) {
-  return sinh(x);
-}
-double test_sqrt(double x) {
-  return sqrt(x);
-}
-double test_tan(double x) {
-  return tan(x);
-}
-double test_tanh(double x) {
-  return tanh(x);
-}
+double test_acos(double x) { return acos(x); }
+double test_acosh(double x) { return acosh(x); }
+double test_asin(double x) { return asin(x); }
+double test_asinh(double x) { return asinh(x); }
+double test_atan(double x) { return atan(x); }
+double test_atan2(double y, double x) { return atan2(y, x); }
+double test_atanh(double x) { return atanh(x); }
+double test_ceil(double x) { return ceil(x); }
+double test_copysign(double x, double y) { return copysign(x, y); }
+double test_cos(double x) { return cos(x); }
+double test_cosh(double x) { return cosh(x); }
+double test_erf(double x) { return erf(x); }
+double test_erfc(double x) { return erfc(x); }
+double test_exp(double x) { return exp(x); }
+double test_expm1(double x) { return expm1(x); }
+double test_fabs(double x) { return fabs(x); }
+double test_floor(double x) { return floor(x); }
+double test_fmod(double x, double y) { return fmod(x, y); }
+double test_frexp(double x, int *exp) { return frexp(x, exp); }
+double test_ldexp(double x, int exp) { return ldexp(x, exp); }
+double test_lgamma(double x) { return lgamma(x); }
+double test_lgamma_r(double x, int *signgamp) { return lgamma_r(x, signgamp); }
+double test_log(double x) { return log(x); }
+double test_log1p(double x) { return log1p(x); }
+double test_log10(double x) { return log10(x); }
+double test_log2(double x) { return log2(x); }
+double test_logb(double x) { return logb(x); }
+double test_modf(double x, double *iptr) { return modf(x, iptr); }
+double test_nan(const char *s) { return nan(s); }
+double test_nearbyint(double x) { return nearbyint(x); }
+double test_pow(double x, double y) { return pow(x, y); }
+double test_rint(double x) { return rint(x); }
+double test_round(double x) { return round(x); }
+double test_scalb(double x, double exp) { return scalb(x, exp); }
+double test_scalbn(double x, int exp) { return scalbn(x, exp); }
+double test_sin(double x) { return sin(x); }
+double test_sinh(double x) { return sinh(x); }
+double test_sqrt(double x) { return sqrt(x); }
+double test_tan(double x) { return tan(x); }
+double test_tanh(double x) { return tanh(x); }
+double test_tgamma(double x) { return tgamma(x); }
+double test_trunc(double x) { return trunc(x); }
--- a/libaxx/include/cmath
+++ b/libaxx/include/cmath
@@ -3,11 +3,13 @@

 #include <math.h>

+#undef finite
 #undef fpclassify
 #undef isfinite
 #undef isnormal
 #undef isnan
 #undef isinf
+#undef signbit

 #undef acosf
 #undef acoshf
@@ -26,12 +28,15 @@
 #undef fmaxf
 #undef floorf
 #undef fmodf
+#undef frexpf
 #undef hypotf
+#undef ldexpf
 #undef lgammaf
 #undef lgammaf_r
 #undef log1pf
 #undef log10f
 #undef logf
+#undef modff
 #undef nanf
 #undef nearbyintf
 #undef powf
@@ -42,12 +47,14 @@
 #undef sqrtf
 #undef tanf
 #undef tanhf
+#undef truncf

 #undef acos
 #undef acosh
 #undef asin
 #undef asinh
 #undef atan
+#undef atan2
 #undef atanh
 #undef ceil
 #undef copysign
@@ -61,88 +68,129 @@
 #undef fmax
 #undef floor
 #undef fmod
+#undef frexp
 #undef hypot
+#undef ldexp
 #undef lgamma
 #undef lgamma_r
+#undef log
 #undef log1p
 #undef log10
-#undef log
+#undef log2
+#undef logb
+#undef modf
+#undef nan
+#undef nearbyint
 #undef pow
+#undef rint
 #undef round
+#undef scalb
 #undef scalbn
 #undef sin
 #undef sinh
 #undef sqrt
 #undef tan
 #undef tanh
+#undef tgamma
+#undef trunc

 namespace std {

+inline constexpr int fpclassify(float x) { return __builtin_fpclassify(FP_NAN, FP_INFINITE, FP_NORMAL, FP_SUBNORMAL, FP_ZERO, x); }
+inline constexpr bool isfinite(float x) { return __builtin_isfinite(x); }
+inline constexpr bool isinf(float x) { return __builtin_isinf(x); }
+inline constexpr bool isnan(float x) { return __builtin_isnan(x); }
+inline constexpr bool isnormal(float x) { return __builtin_isnormal(x); }
+inline constexpr bool signbit(float x) { return __builtin_signbit(x); }

-static inline double acos(double x) { return __builtin_acos(x); }
-static inline float acos(float x) { return __builtin_acosf(x); }
-static inline double acosh(double x) { return __builtin_acosh(x); }
-static inline float acosh(float x) { return __builtin_acoshf(x); }
-static inline double asin(double x) { return __builtin_asin(x); }
-static inline float asin(float x) { return __builtin_asinf(x); }
-static inline double asinh(double x) { return __builtin_asinh(x); }
-static inline float asinh(float x) { return __builtin_asinhf(x); }
-static inline double atan(double x) { return __builtin_atan(x); }
-static inline float atan(float x) { return __builtin_atanf(x); }
-static inline double atanh(double x) { return __builtin_atanh(x); }
-static inline float atanh(float x) { return __builtin_atanhf(x); }
-static inline double ceil(double x) { return __builtin_ceil(x); }
-static inline float ceil(float x) { return __builtin_ceilf(x); }
-static inline double copysign(double x, double y) { return __builtin_copysign(x, y); }
-static inline float copysign(float x, float y) { return __builtin_copysignf(x, y); }
-static inline double cos(double x) { return __builtin_cos(x); }
-static inline float cos(float x) { return __builtin_cosf(x); }
-static inline double cosh(double x) { return __builtin_cosh(x); }
-static inline float cosh(float x) { return __builtin_coshf(x); }
-static inline double erf(double x) { return __builtin_erf(x); }
-static inline double erfc(double x) { return __builtin_erfc(x); }
-static inline double exp(double x) { return __builtin_exp(x); }
-static inline float exp(float x) { return __builtin_expf(x); }
-static inline double fabs(double x) { return __builtin_fabs(x); }
-static inline float fabs(float x) { return __builtin_fabsf(x); }
-static inline double fmax(double x, double y) { return __builtin_fmax(x, y); }
-static inline float fmax(float x, float y) { return __builtin_fmaxf(x, y); }
-static inline double floor(double x) { return __builtin_floor(x); }
-static inline float floor(float x) { return __builtin_floorf(x); }
-static inline double fmod(double x, double y) { return __builtin_fmod(x, y); }
-static inline float fmod(float x, float y) { return __builtin_fmodf(x, y); }
-static inline int fpclassify(double x) { return __builtin_fpclassify(FP_NAN, FP_INFINITE, FP_NORMAL, FP_SUBNORMAL, FP_ZERO, x); }
-static inline int fpclassify(float x) { return __builtin_fpclassify(FP_NAN, FP_INFINITE, FP_NORMAL, FP_SUBNORMAL, FP_ZERO, x); }
-static inline double hypot(double x, double y) { return __builtin_hypot(x, y); }
-static inline float hypotf(float x, float y) { return __builtin_hypotf(x, y); }
-static inline bool isfinite(double x) { return __builtin_isfinite(x); }
-static inline bool isfinite(float x) { return __builtin_isfinite(x); }
-static inline bool isinf(double x) { return __builtin_isinf(x); }
-static inline bool isinf(float x) { return __builtin_isinf(x); }
-static inline bool isnan(double x) { return __builtin_isnan(x); }
-static inline bool isnan(float x) { return __builtin_isnan(x); }
-static inline bool isnormal(double x) { return __builtin_isnormal(x); }
-static inline bool isnormal(float x) { return __builtin_isnormal(x); }
-static inline double lgamma(double x) { return __builtin_lgamma(x); }
-static inline float lgamma(float x) { return __builtin_lgammaf(x); }
-static inline double log10(double x) { return __builtin_log10(x); }
-static inline float log10(float x) { return __builtin_log10f(x); }
-static inline double log(double x) { return __builtin_log(x); }
-static inline float log(float x) { return __builtin_logf(x); }
-static inline double pow(double x, double y) { return __builtin_pow(x, y); }
-static inline float pow(float x, float y) { return __builtin_powf(x, y); }
-static inline double round(double x) { return __builtin_round(x); }
-static inline float round(float x) { return __builtin_roundf(x); }
-static inline double sin(double x) { return __builtin_sin(x); }
-static inline float sin(float x) { return __builtin_sinf(x); }
-static inline double sinh(double x) { return __builtin_sinh(x); }
-static inline float sinh(float x) { return __builtin_sinhf(x); }
-static inline double sqrt(double x) { return __builtin_sqrt(x); }
-static inline float sqrt(float x) { return __builtin_sqrtf(x); }
-static inline double tan(double x) { return __builtin_tan(x); }
-static inline float tan(float x) { return __builtin_tanf(x); }
-static inline double tanh(double x) { return __builtin_tanh(x); }
-static inline float tanh(float x) { return __builtin_tanhf(x); }
+inline constexpr float acos(float x) { return __builtin_acosf(x); }
+inline constexpr float acosh(float x) { return __builtin_acoshf(x); }
+inline constexpr float asin(float x) { return __builtin_asinf(x); }
+inline constexpr float asinh(float x) { return __builtin_asinhf(x); }
+inline constexpr float atan(float x) { return __builtin_atanf(x); }
+inline constexpr float atan2(float y, float x) { return __builtin_atan2f(y, x); }
+inline constexpr float atanh(float x) { return __builtin_atanhf(x); }
+inline constexpr float ceil(float x) { return __builtin_ceilf(x); }
+inline constexpr float copysign(float x, float y) { return __builtin_copysignf(x, y); }
+inline constexpr float cos(float x) { return __builtin_cosf(x); }
+inline constexpr float cosh(float x) { return __builtin_coshf(x); }
+inline constexpr float exp(float x) { return __builtin_expf(x); }
+inline constexpr float expm1(float x) { return __builtin_expm1f(x); }
+inline constexpr float fabs(float x) { return __builtin_fabsf(x); }
+inline constexpr float floor(float x) { return __builtin_floorf(x); }
+inline constexpr float fmax(float x, float y) { return __builtin_fmaxf(x, y); }
+inline constexpr float fmod(float x, float y) { return __builtin_fmodf(x, y); }
+inline constexpr float frexp(float x, int *exp) { return __builtin_frexpf(x, exp); }
+inline constexpr float ldexp(float x, int exp) { return __builtin_ldexpf(x, exp); }
+inline constexpr float lgamma(float x) { return __builtin_lgammaf(x); }
+inline constexpr float lgamma_r(float x, int *signgamp) { return __builtin_lgammaf_r(x, signgamp); }
+inline constexpr float log1p(float x) { return __builtin_log1pf(x); }
+inline constexpr float log10(float x) { return __builtin_log10f(x); }
+inline constexpr float log(float x) { return __builtin_logf(x); }
+inline constexpr float modf(float x, float *iptr) { return __builtin_modff(x, iptr); }
+inline constexpr float nanf(const char *tagp) { return __builtin_nanf(tagp); }
+inline constexpr float nearbyint(float x) { return __builtin_nearbyintf(x); }
+inline constexpr float pow(float x, float y) { return __builtin_powf(x, y); }
+inline constexpr float round(float x) { return __builtin_roundf(x); }
+inline constexpr float scalbn(float x, int exp) { return __builtin_scalbnf(x, exp); }
+inline constexpr float sin(float x) { return __builtin_sinf(x); }
+inline constexpr float sinh(float x) { return __builtin_sinhf(x); }
+inline constexpr float sqrt(float x) { return __builtin_sqrtf(x); }
+inline constexpr float tan(float x) { return __builtin_tanf(x); }
+inline constexpr float tanh(float x) { return __builtin_tanhf(x); }
+inline constexpr float trunc(float x) { return __builtin_truncf(x); }
+
+inline constexpr int fpclassify(double x) { return __builtin_fpclassify(FP_NAN, FP_INFINITE, FP_NORMAL, FP_SUBNORMAL, FP_ZERO, x); }
+inline constexpr bool isfinite(double x) { return __builtin_isfinite(x); }
+inline constexpr bool isinf(double x) { return __builtin_isinf(x); }
+inline constexpr bool isnan(double x) { return __builtin_isnan(x); }
+inline constexpr bool isnormal(double x) { return __builtin_isnormal(x); }
+inline constexpr bool signbit(double x) { return __builtin_signbit(x); }
+
+inline constexpr double acos(double x) { return __builtin_acos(x); }
+inline constexpr double acosh(double x) { return __builtin_acosh(x); }
+inline constexpr double asin(double x) { return __builtin_asin(x); }
+inline constexpr double asinh(double x) { return __builtin_asinh(x); }
+inline constexpr double atan(double x) { return __builtin_atan(x); }
+inline constexpr double atan2(double y, double x) { return __builtin_atan2(y, x); }
+inline constexpr double atanh(double x) { return __builtin_atanh(x); }
+inline constexpr double ceil(double x) { return __builtin_ceil(x); }
+inline constexpr double copysign(double x, double y) { return __builtin_copysign(x, y); }
+inline constexpr double cos(double x) { return __builtin_cos(x); }
+inline constexpr double cosh(double x) { return __builtin_cosh(x); }
+inline constexpr double erf(double x) { return __builtin_erf(x); }
+inline constexpr double erfc(double x) { return __builtin_erfc(x); }
+inline constexpr double exp(double x) { return __builtin_exp(x); }
+inline constexpr double expm1(double x) { return __builtin_expm1(x); }
+inline constexpr double fabs(double x) { return __builtin_fabs(x); }
+inline constexpr double floor(double x) { return __builtin_floor(x); }
+inline constexpr double fmax(double x, double y) { return __builtin_fmax(x, y); }
+inline constexpr double fmod(double x, double y) { return __builtin_fmod(x, y); }
+inline constexpr double frexp(double x, int *exp) { return __builtin_frexp(x, exp); }
+inline constexpr double hypot(double x, double y) { return __builtin_hypot(x, y); }
+inline constexpr double ldexp(double x, int exp) { return __builtin_scalbn(x, exp); }
+inline constexpr double lgamma(double x) { return __builtin_lgamma(x); }
+inline constexpr double lgamma_r(double x, int *signgamp) { return __builtin_lgamma_r(x, signgamp); }
+inline constexpr double log(double x) { return __builtin_log(x); }
+inline constexpr double log1p(double x) { return __builtin_log1p(x); }
+inline constexpr double log10(double x) { return __builtin_log10(x); }
+inline constexpr double log2(double x) { return __builtin_log2(x); }
+inline constexpr double logb(double x) { return __builtin_logb(x); }
+inline constexpr double modf(double x, double *iptr) { return __builtin_modf(x, iptr); }
+inline constexpr double nan(const char *tagp) { return __builtin_nan(tagp); }
+inline constexpr double nearbyint(double x) { return __builtin_nearbyint(x); }
+inline constexpr double pow(double x, double y) { return __builtin_pow(x, y); }
+inline constexpr double rint(double x) { return __builtin_rint(x); }
+inline constexpr double round(double x) { return __builtin_round(x); }
+inline constexpr double scalb(double x, double exp) { return __builtin_scalb(x, exp); }
+inline constexpr double scalbn(double x, int exp) { return __builtin_scalbn(x, exp); }
+inline constexpr double sin(double x) { return __builtin_sin(x); }
+inline constexpr double sinh(double x) { return __builtin_sinh(x); }
+inline constexpr double sqrt(double x) { return __builtin_sqrt(x); }
+inline constexpr double tan(double x) { return __builtin_tan(x); }
+inline constexpr double tanh(double x) { return __builtin_tanh(x); }
+inline constexpr double tgamma(double x) { return __builtin_tgamma(x); }
+inline constexpr double trunc(double x) { return __builtin_trunc(x); }

 }

--- a/poincare/include/poincare/approximation_helper.h
+++ b/poincare/include/poincare/approximation_helper.h
@@ -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> NeglectRealOrImaginaryPartIfNeglectable(std::complex<T> result, std::complex<T> input1, std::complex<T> input2 = 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);
--- a/poincare/include/poincare/power.h
+++ b/poincare/include/poincare/power.h
@@ -34,6 +34,7 @@ public:
  int polynomialDegree(Context * context, const char * symbolName) const override;
  int getPolynomialCoefficients(Context * context, const char * symbolName, Expression coefficients[], ExpressionNode::SymbolicComputation symbolicComputation) const override;

+  template<typename T> static Complex<T> computeNotPrincipalRealRootOfRationalPow(const std::complex<T> c, T p, T q);
  template<typename T> static Complex<T> compute(const std::complex<T> c, const std::complex<T> d, Preferences::ComplexFormat complexFormat);

 private:
@@ -59,11 +60,12 @@ private:
  template<typename T> static MatrixComplex<T> computeOnMatrixAndComplex(const MatrixComplex<T> m, const std::complex<T> d, Preferences::ComplexFormat complexFormat);
  template<typename T> static MatrixComplex<T> computeOnMatrices(const MatrixComplex<T> m, const MatrixComplex<T> n, Preferences::ComplexFormat complexFormat);
  Evaluation<float> approximate(SinglePrecision p, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const override {
-    return ApproximationHelper::MapReduce<float>(this, context, complexFormat, angleUnit, compute<float>, computeOnComplexAndMatrix<float>, computeOnMatrixAndComplex<float>, computeOnMatrices<float>);
+    return templatedApproximate<float>(context, complexFormat, angleUnit);
  }
  Evaluation<double> approximate(DoublePrecision p, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const override {
-    return ApproximationHelper::MapReduce<double>(this, context, complexFormat, angleUnit, compute<double>, computeOnComplexAndMatrix<double>, computeOnMatrixAndComplex<double>, computeOnMatrices<double>);
+    return templatedApproximate<double>(context, complexFormat, angleUnit);
  }
+ template<typename T> Evaluation<T> templatedApproximate(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const;
 };

 class Power final : public Expression {
--- a/poincare/include/poincare/trigonometry.h
+++ b/poincare/include/poincare/trigonometry.h
@@ -23,9 +23,6 @@ public:
  static Expression table(const Expression e, ExpressionNode::Type type,  ExpressionNode::ReductionContext reductionContext); // , Function f, bool inverse
  template <typename T> static std::complex<T> ConvertToRadian(const std::complex<T> c, Preferences::AngleUnit angleUnit);
  template <typename T> static std::complex<T> ConvertRadianToAngleUnit(const std::complex<T> c, Preferences::AngleUnit angleUnit);
-  template <typename T> static std::complex<T> RoundToMeaningfulDigits(const std::complex<T> result, const std::complex<T> input);
-private:
-  template <typename T> static T RoundToMeaningfulDigits(T result, T input);
 };

 }
--- a/poincare/src/approximation_helper.cpp
+++ b/poincare/src/approximation_helper.cpp
@@ -10,9 +10,19 @@ extern "C" {

 namespace Poincare {

-template <typename T> T absMod(T a, T b) {
-  T result = std::fmod(std::fabs(a), b);
-  return result > b/2 ? b-result : result;
+template <typename T> bool isNegligeable(T x, T precision, T norm1, T norm2) {
+  T absX = std::fabs(x);
+  return absX <= 10.0*precision && absX/norm1 <= precision && absX/norm2 <= precision;
+}
+
+template < typename T> T minimalNonNullMagnitudeOfParts(std::complex<T> c) {
+  T absRealInput = std::fabs(c.real());
+  T absImagInput = std::fabs(c.imag());
+  // If the magnitude of one part is null, ignore it
+  if (absRealInput == 0.0 || (absImagInput > 0.0 && absImagInput < absRealInput)) {
+    return absImagInput;
+  }
+  return absRealInput;
 }

 static inline int absInt(int x) { return x < 0 ? -x : x; }
@@ -27,16 +37,26 @@ template <typename T> int ApproximationHelper::PositiveIntegerApproximationIfPos
  return absInt((int)scalar);
 }

-template <typename T> std::complex<T> ApproximationHelper::TruncateRealOrImaginaryPartAccordingToArgument(std::complex<T> c) {
-  T arg = std::arg(c);
-  T precision = 10*Expression::Epsilon<T>();
-  if (absMod<T>(arg, (T)M_PI) <= precision) {
-    c.imag(0);
+template <typename T> std::complex<T> ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(std::complex<T> result, std::complex<T> input1, std::complex<T> input2) {
+  /* Cheat: openbsd  functions (cos, sin, tan, cosh, acos, pow...) 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 acos(1) = 6E-17 or cos(π/2) = 4E-17, we round the result
+   * to its 1E-6 or 1E-14 precision when its ratio with the argument (π/2 in the
+   * example) is smaller than epsilon. This way, we have sin(π) ~ 0 and
+   * sin(1E-15)=1E-15.
+   * We can't do that for all evaluation as the user can operate on values as
+   * small as 1E-308 (in double) and most results still be correct. */
+  T magnitude1 = minimalNonNullMagnitudeOfParts(input1);
+  T magnitude2 = minimalNonNullMagnitudeOfParts(input2);
+  T precision = 10.0*Expression::Epsilon<T>();
+  if (isNegligeable(result.imag(), precision, magnitude1, magnitude2)) {
+    result.imag(0);
  }
-  if (absMod<T>(arg-(T)M_PI/2.0, (T)M_PI) <= precision) {
-    c.real(0);
+  if (isNegligeable(result.real(), precision, magnitude1, magnitude2)) {
+    result.real(0);
  }
-  return c;
+  return result;
 }

 template<typename T> Evaluation<T> ApproximationHelper::Map(const ExpressionNode * expression, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, ComplexCompute<T> compute) {
@@ -106,8 +126,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::NeglectRealOrImaginaryPartIfNeglectable<float>(std::complex<float>,std::complex<float>,std::complex<float>);
+template std::complex<double> Poincare::ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable<double>(std::complex<double>,std::complex<double>,std::complex<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);
--- a/poincare/src/arc_cosine.cpp
+++ b/poincare/src/arc_cosine.cpp
@@ -44,7 +44,7 @@ Complex<T> ArcCosineNode::computeOnComplex(const std::complex<T> c, Preferences:
      result.imag(-result.imag()); // other side of the cut
    }
  }
-  result = Trigonometry::RoundToMeaningfulDigits(result, c);
+  result = ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(result, c);
  return Complex<T>::Builder(Trigonometry::ConvertRadianToAngleUnit(result, angleUnit));
 }

--- a/poincare/src/arc_sine.cpp
+++ b/poincare/src/arc_sine.cpp
@@ -44,7 +44,7 @@ Complex<T> ArcSineNode::computeOnComplex(const std::complex<T> c, Preferences::C
      result.imag(-result.imag()); // other side of the cut
    }
  }
-  result = Trigonometry::RoundToMeaningfulDigits(result, c);
+  result = ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(result, c);
  return Complex<T>::Builder(Trigonometry::ConvertRadianToAngleUnit(result, angleUnit));
 }

--- a/poincare/src/arc_tangent.cpp
+++ b/poincare/src/arc_tangent.cpp
@@ -40,7 +40,7 @@ Complex<T> ArcTangentNode::computeOnComplex(const std::complex<T> c, Preferences
      result.real(-result.real()); // other side of the cut
    }
  }
-  result = Trigonometry::RoundToMeaningfulDigits(result, c);
+  result = ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(result, c);
  return Complex<T>::Builder(Trigonometry::ConvertRadianToAngleUnit(result, angleUnit));
 }

--- a/poincare/src/arithmetic.cpp
+++ b/poincare/src/arithmetic.cpp
@@ -7,6 +7,9 @@ Integer Arithmetic::LCM(const Integer & a, const Integer & b) {
  if (a.isZero() || b.isZero()) {
    return Integer(0);
  }
+  if (a.isEqualTo(b)) {
+    return a;
+  }
  Integer signResult = Integer::Division(Integer::Multiplication(a, b), GCD(a, b)).quotient;
  signResult.setNegative(false);
  return signResult;
@@ -16,7 +19,9 @@ Integer Arithmetic::GCD(const Integer & a, const Integer & b) {
  if (a.isOverflow() || b.isOverflow()) {
    return Integer::Overflow(false);
  }
-
+  if (a.isEqualTo(b)) {
+    return a;
+  }
  Integer i = a;
  Integer j = b;
  i.setNegative(false);
--- a/poincare/src/cosine.cpp
+++ b/poincare/src/cosine.cpp
@@ -19,7 +19,7 @@ template<typename T>
 Complex<T> CosineNode::computeOnComplex(const std::complex<T> c, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit) {
  std::complex<T> angleInput = Trigonometry::ConvertToRadian(c, angleUnit);
  std::complex<T> res = std::cos(angleInput);
-  return Complex<T>::Builder(Trigonometry::RoundToMeaningfulDigits(res, angleInput));
+  return Complex<T>::Builder(ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(res, angleInput));
 }

 Layout CosineNode::createLayout(Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const {
--- a/poincare/src/hyperbolic_arc_cosine.cpp
+++ b/poincare/src/hyperbolic_arc_cosine.cpp
@@ -23,7 +23,7 @@ Complex<T> HyperbolicArcCosineNode::computeOnComplex(const std::complex<T> c, Pr
   * on this cut. We followed the convention chosen by the lib c++ of llvm on
   * ]-inf+0i, 1+0i] (warning: atanh takes the other side of the cut values on
   * ]-inf-0i, 1-0i[).*/
-  return Complex<T>::Builder(Trigonometry::RoundToMeaningfulDigits(result, c));
+  return Complex<T>::Builder(ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(result, c));
 }

 template Complex<float> Poincare::HyperbolicArcCosineNode::computeOnComplex<float>(std::complex<float>, Preferences::ComplexFormat, Preferences::AngleUnit);
--- a/poincare/src/hyperbolic_arc_sine.cpp
+++ b/poincare/src/hyperbolic_arc_sine.cpp
@@ -27,7 +27,7 @@ Complex<T> HyperbolicArcSineNode::computeOnComplex(const std::complex<T> c, Pref
  if (c.real() == 0 && c.imag() < 1) {
    result.real(-result.real()); // other side of the cut
  }
-  return Complex<T>::Builder(Trigonometry::RoundToMeaningfulDigits(result, c));
+  return Complex<T>::Builder(ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(result, c));
 }

 template Complex<float> Poincare::HyperbolicArcSineNode::computeOnComplex<float>(std::complex<float>, Preferences::ComplexFormat, Preferences::AngleUnit);
--- a/poincare/src/hyperbolic_arc_tangent.cpp
+++ b/poincare/src/hyperbolic_arc_tangent.cpp
@@ -27,7 +27,7 @@ Complex<T> HyperbolicArcTangentNode::computeOnComplex(const std::complex<T> c, P
  if (c.imag() == 0 && c.real() > 1) {
    result.imag(-result.imag()); // other side of the cut
  }
-  return Complex<T>::Builder(Trigonometry::RoundToMeaningfulDigits(result, c));
+  return Complex<T>::Builder(ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(result, c));
 }

 template Complex<float> Poincare::HyperbolicArcTangentNode::computeOnComplex<float>(std::complex<float>, Preferences::ComplexFormat, Preferences::AngleUnit);
--- a/poincare/src/hyperbolic_cosine.cpp
+++ b/poincare/src/hyperbolic_cosine.cpp
@@ -16,7 +16,7 @@ int HyperbolicCosineNode::serialize(char * buffer, int bufferSize, Preferences::

 template<typename T>
 Complex<T> HyperbolicCosineNode::computeOnComplex(const std::complex<T> c, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit) {
-  return Complex<T>::Builder(Trigonometry::RoundToMeaningfulDigits(std::cosh(c), c));
+  return Complex<T>::Builder(ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(std::cosh(c), c));
 }

 template Complex<float> Poincare::HyperbolicCosineNode::computeOnComplex<float>(std::complex<float>, Preferences::ComplexFormat, Preferences::AngleUnit);
--- a/poincare/src/hyperbolic_sine.cpp
+++ b/poincare/src/hyperbolic_sine.cpp
@@ -16,7 +16,7 @@ int HyperbolicSineNode::serialize(char * buffer, int bufferSize, Preferences::Pr

 template<typename T>
 Complex<T> HyperbolicSineNode::computeOnComplex(const std::complex<T> c, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit) {
-  return Complex<T>::Builder(Trigonometry::RoundToMeaningfulDigits(std::sinh(c), c));
+  return Complex<T>::Builder(ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(std::sinh(c), c));
 }

 template Complex<float> Poincare::HyperbolicSineNode::computeOnComplex<float>(std::complex<float>, Preferences::ComplexFormat, Preferences::AngleUnit);
--- a/poincare/src/hyperbolic_tangent.cpp
+++ b/poincare/src/hyperbolic_tangent.cpp
@@ -16,7 +16,7 @@ int HyperbolicTangentNode::serialize(char * buffer, int bufferSize, Preferences:

 template<typename T>
 Complex<T> HyperbolicTangentNode::computeOnComplex(const std::complex<T> c, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit) {
-  return Complex<T>::Builder(Trigonometry::RoundToMeaningfulDigits(std::tanh(c), c));
+  return Complex<T>::Builder(ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(std::tanh(c), c));
 }

 template Complex<float> Poincare::HyperbolicTangentNode::computeOnComplex<float>(std::complex<float>, Preferences::ComplexFormat, Preferences::AngleUnit);
--- a/poincare/src/multiplication.cpp
+++ b/poincare/src/multiplication.cpp
@@ -607,8 +607,10 @@ Expression Multiplication::privateShallowReduce(ExpressionNode::ReductionContext
   * which also are multiplications themselves. */
  mergeMultiplicationChildrenInPlace();

+  Context * context = reductionContext.context();
+
  // Step 2: Sort the children
-  sortChildrenInPlace([](const ExpressionNode * e1, const ExpressionNode * e2, bool canBeInterrupted) { return ExpressionNode::SimplificationOrder(e1, e2, true, canBeInterrupted); }, reductionContext.context(), true);
+  sortChildrenInPlace([](const ExpressionNode * e1, const ExpressionNode * e2, bool canBeInterrupted) { return ExpressionNode::SimplificationOrder(e1, e2, true, canBeInterrupted); }, context, true);

  // Step 3: Handle matrices
  /* Thanks to the simplification order, all matrix children (if any) are the
@@ -683,7 +685,7 @@ Expression Multiplication::privateShallowReduce(ExpressionNode::ReductionContext
     * interval). */

    if (multiplicationChildIndex >= 0) {
-      if (childAtIndex(multiplicationChildIndex).deepIsMatrix(reductionContext.context())) {
+      if (childAtIndex(multiplicationChildIndex).deepIsMatrix(context)) {
        return *this;
      }
      removeChildInPlace(resultMatrix, resultMatrix.numberOfChildren());
@@ -696,7 +698,7 @@ Expression Multiplication::privateShallowReduce(ExpressionNode::ReductionContext
      }
    }
    replaceWithInPlace(resultMatrix);
-    return resultMatrix.shallowReduce(reductionContext.context());
+    return resultMatrix.shallowReduce(context);
  }

  /* Step 4: Gather like terms. For example, turn pi^2*pi^3 into pi^5. Thanks to
@@ -706,7 +708,7 @@ Expression Multiplication::privateShallowReduce(ExpressionNode::ReductionContext
  while (i < numberOfChildren()-1) {
    Expression oi = childAtIndex(i);
    Expression oi1 = childAtIndex(i+1);
-    if (oi.recursivelyMatches(Expression::IsRandom, reductionContext.context(), true)) {
+    if (oi.recursivelyMatches(Expression::IsRandom, context, true)) {
      // Do not factorize random or randint
      i++;
      continue;
@@ -754,7 +756,7 @@ Expression Multiplication::privateShallowReduce(ExpressionNode::ReductionContext
    /* Replacing sin/cos by tan factors may have mixed factors and factors are
     * guaranteed to be sorted (according ot SimplificationOrder) at the end of
     * shallowReduce */
-    sortChildrenInPlace([](const ExpressionNode * e1, const ExpressionNode * e2, bool canBeInterrupted) { return ExpressionNode::SimplificationOrder(e1, e2, true, canBeInterrupted); }, reductionContext.context(), true);
+    sortChildrenInPlace([](const ExpressionNode * e1, const ExpressionNode * e2, bool canBeInterrupted) { return ExpressionNode::SimplificationOrder(e1, e2, true, canBeInterrupted); }, context, true);
  }

  /* Step 6: We remove rational children that appeared in the middle of sorted
@@ -777,6 +779,14 @@ Expression Multiplication::privateShallowReduce(ExpressionNode::ReductionContext
      if (childAtIndex(0).isNumber()) {
        Number o0 = childAtIndex(0).convert<Rational>();
        Number m = Number::Multiplication(o0, static_cast<Number &>(o));
+        if ((IsInfinity(m, context) || m.isUndefined())
+            && !IsInfinity(o0, context) && !o0.isUndefined()
+            && !IsInfinity(o, context) && !o.isUndefined())
+        {
+          // Stop the reduction due to a multiplication overflow
+          SetInterruption(true);
+          return *this;
+        }
        replaceChildAtIndexInPlace(0, m);
        removeChildAtIndexInPlace(i);
      } else {
@@ -798,7 +808,7 @@ Expression Multiplication::privateShallowReduce(ExpressionNode::ReductionContext
    const Expression c = childAtIndex(0);
    if (c.type() == ExpressionNode::Type::Rational && static_cast<const Rational &>(c).isZero()) {
      // Check that other children don't match inf or unit
-      bool infiniteOrUnitFactor = recursivelyMatches([](const Expression e, Context * context) { return Expression::IsInfinity(e,context) || e.type() == ExpressionNode::Type::Unit; }, reductionContext.context());
+      bool infiniteOrUnitFactor = recursivelyMatches([](const Expression e, Context * context) { return Expression::IsInfinity(e,context) || e.type() == ExpressionNode::Type::Unit; }, context);
      if (!infiniteOrUnitFactor) {
        replaceWithInPlace(c);
        return c;
@@ -817,7 +827,7 @@ Expression Multiplication::privateShallowReduce(ExpressionNode::ReductionContext
   * reduce expressions such as (x+y)^(-1)*(x+y)(a+b).
   * If there is a random somewhere, do not expand. */
  Expression p = parent();
-  bool hasRandom = recursivelyMatches(Expression::IsRandom, reductionContext.context(), true);
+  bool hasRandom = recursivelyMatches(Expression::IsRandom, context, true);
  if (shouldExpand
      && (p.isUninitialized() || p.type() != ExpressionNode::Type::Multiplication)
      && !hasRandom)
@@ -844,7 +854,7 @@ Expression Multiplication::privateShallowReduce(ExpressionNode::ReductionContext
   * - All children are either real or ComplexCartesian (allChildrenAreReal == 0)
   *   We can bubble up ComplexCartesian nodes.
   * Do not simplify if there are randoms !*/
-  if (!hasRandom && allChildrenAreReal(reductionContext.context()) == 0) {
+  if (!hasRandom && allChildrenAreReal(context) == 0) {
    int nbChildren = numberOfChildren();
    int i = nbChildren-1;
    // Children are sorted so ComplexCartesian nodes are at the end
--- a/poincare/src/nth_root.cpp
+++ b/poincare/src/nth_root.cpp
@@ -46,18 +46,12 @@ Evaluation<T> NthRootNode::templatedApproximate(Context * context, Preferences::
    /* If the complexFormat is Real, we look for nthroot of form root(x,q) with
     * x real and q integer because they might have a real form which does not
     * correspond to the principale angle. */
-    if (complexFormat == Preferences::ComplexFormat::Real) {
+    if (complexFormat == Preferences::ComplexFormat::Real && indexc.imag() == 0.0 && std::round(indexc.real()) == indexc.real()) {
      // root(x, q) with q integer and x real
-      if (basec.imag() == 0.0 && indexc.imag() == 0.0 && std::round(indexc.real()) == indexc.real()) {
-        std::complex<T> absBasec = basec;
-        absBasec.real(std::fabs(absBasec.real()));
-        // compute root(|x|, q)
-        Complex<T> absBasePowIndex = PowerNode::compute(absBasec, std::complex<T>(1.0)/(indexc), complexFormat);
-        // q odd if (-1)^q = -1
-        if (std::pow((T)-1.0, (T)indexc.real()) < 0.0) {
-          return basec.real() < 0 ? Complex<T>::Builder(-absBasePowIndex.stdComplex()) : absBasePowIndex;
-        }
-      }
+      Complex<T> result = PowerNode::computeNotPrincipalRealRootOfRationalPow(basec, (T)1.0, indexc.real());
+       if (!result.isUndefined()) {
+         return result;
+       }
    }
    result = PowerNode::compute(basec, std::complex<T>(1.0)/(indexc), complexFormat);
  }
--- a/poincare/src/power.cpp
+++ b/poincare/src/power.cpp
@@ -127,6 +127,32 @@ bool PowerNode::childAtIndexNeedsUserParentheses(const Expression & child, int c

 // 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() == 0 && std::pow((T)-1.0, q) < 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() < 0 && std::pow((T)-1.0, p) < 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;
@@ -151,8 +177,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::NeglectRealOrImaginaryPartIfNeglectable(result, c, d));
 }

 // Layout
@@ -257,6 +288,50 @@ template<typename T> MatrixComplex<T> PowerNode::computeOnMatrices(const MatrixC
  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 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) {
@@ -896,20 +971,6 @@ Expression Power::shallowBeautify(ExpressionNode::ReductionContext reductionCont
    return result;
  }

-  /* Optional Step 3: if the ReductionTarget is the SystemForApproximation or
-   * SystemForAnalysis, turn a^(p/q) into (root(a, q))^p
-   * Indeed, root(a, q) can have a real root which is not the principale angle
-   * but that we want to return in real complex format. This special case is
-   * handled in NthRoot approximation but not in Power approximation. */
-  if (reductionContext.target() != ExpressionNode::ReductionTarget::User && childAtIndex(1).type() == ExpressionNode::Type::Rational) {
-    Integer p = childAtIndex(1).convert<Rational>().signedIntegerNumerator();
-    Integer q = childAtIndex(1).convert<Rational>().integerDenominator();
-    Expression nthRoot = q.isOne() ? childAtIndex(0) : NthRoot::Builder(childAtIndex(0), Rational::Builder(q));
-    Expression result = p.isOne() ? nthRoot : Power::Builder(nthRoot, Rational::Builder(p));
-    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));
--- a/poincare/src/sine.cpp
+++ b/poincare/src/sine.cpp
@@ -19,7 +19,7 @@ template<typename T>
 Complex<T> SineNode::computeOnComplex(const std::complex<T> c, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit) {
  std::complex<T> angleInput = Trigonometry::ConvertToRadian(c, angleUnit);
  std::complex<T> res = std::sin(angleInput);
-  return Complex<T>::Builder(Trigonometry::RoundToMeaningfulDigits(res, angleInput));
+  return Complex<T>::Builder(ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(res, angleInput));
 }

 Layout SineNode::createLayout(Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const {
--- a/poincare/src/square_root.cpp
+++ b/poincare/src/square_root.cpp
@@ -30,13 +30,7 @@ int SquareRootNode::serialize(char * buffer, int bufferSize, Preferences::PrintF
 template<typename T>
 Complex<T> SquareRootNode::computeOnComplex(const std::complex<T> c, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit) {
  std::complex<T> result = std::sqrt(c);
-  /* 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 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));
+  return Complex<T>::Builder(ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(result, std::complex<T>(std::log(std::abs(c)), std::arg(c))));
 }

 Expression SquareRootNode::shallowReduce(ReductionContext reductionContext) {
--- a/poincare/src/tangent.cpp
+++ b/poincare/src/tangent.cpp
@@ -30,7 +30,7 @@ template<typename T>
 Complex<T> TangentNode::computeOnComplex(const std::complex<T> c, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit) {
  std::complex<T> angleInput = Trigonometry::ConvertToRadian(c, angleUnit);
  std::complex<T> res = std::tan(angleInput);
-  return Complex<T>::Builder(Trigonometry::RoundToMeaningfulDigits(res, angleInput));
+  return Complex<T>::Builder(ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(res, angleInput));
 }

 Expression TangentNode::shallowReduce(ReductionContext reductionContext) {
--- a/poincare/src/trigonometry.cpp
+++ b/poincare/src/trigonometry.cpp
@@ -412,34 +412,9 @@ std::complex<T> Trigonometry::ConvertRadianToAngleUnit(const std::complex<T> c,
  return c;
 }

-template<typename T>
-T Trigonometry::RoundToMeaningfulDigits(T result, T input) {
-  /* Cheat: 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 acos(1) = 6E-17 or cos(π/2) = 4E-17, we round the result
-   * to its 1E-6 or 1E-14 precision when its ratio with the argument (π/2 in the
-   * example) is smaller than epsilon. This way, we have sin(π) ~ 0 and
-   * sin(1E-15)=1E-15.
-   * We can't do that for all evaluation as the user can operate on values as
-   * small as 1E-308 (in double) and most results still be correct. */
-  if (input == 0.0 || std::fabs(result/input) <= Expression::Epsilon<T>()) {
-     T precision = 10*Expression::Epsilon<T>();
-     result = std::round(result/precision)*precision;
-  }
-  return result;
-}
-
-template <typename T>
-std::complex<T> Trigonometry::RoundToMeaningfulDigits(const std::complex<T> result, const std::complex<T> input) {
-  return std::complex<T>(RoundToMeaningfulDigits(result.real(), std::abs(input)), RoundToMeaningfulDigits(result.imag(), std::abs(input)));
-}
-
 template std::complex<float> Trigonometry::ConvertToRadian<float>(std::complex<float>, Preferences::AngleUnit);
 template std::complex<double> Trigonometry::ConvertToRadian<double>(std::complex<double>, Preferences::AngleUnit);
 template std::complex<float> Trigonometry::ConvertRadianToAngleUnit<float>(std::complex<float>, Preferences::AngleUnit);
 template std::complex<double> Trigonometry::ConvertRadianToAngleUnit<double>(std::complex<double>, Preferences::AngleUnit);
-template std::complex<float> Trigonometry::RoundToMeaningfulDigits<float>(std::complex<float>, std::complex<float>);
-template std::complex<double> Trigonometry::RoundToMeaningfulDigits<double>(std::complex<double>, std::complex<double>);

 }
--- a/poincare/test/approximation.cpp
+++ b/poincare/test/approximation.cpp
@@ -149,6 +149,12 @@ 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);
+  assert_expression_simplifies_approximates_to<float>("3.5^2.0000001", "12.25");
+  assert_expression_simplifies_approximates_to<float>("3.7^2.0000001", "13.69");
 }

 QUIZ_CASE(poincare_approximation_subtraction) {
@@ -772,6 +778,13 @@ QUIZ_CASE(poincare_approximation_trigonometry_functions) {
  assert_expression_approximates_to<float>("atanh(𝐢-4)", "-0.238878+1.50862×𝐢", Radian, Cartesian, 6);
  assert_expression_approximates_to<float>("atanh(𝐢-4)", "-0.238878+1.50862×𝐢", Degree, Cartesian, 6);

+  // Check that the complex part is not neglected
+  assert_expression_approximates_to<double>("atanh(0.99999999999+1.0ᴇ-26×𝐢)", "13.01+5ᴇ-16×𝐢", Radian, Cartesian, 4);
+  assert_expression_approximates_to<double>("atanh(0.99999999999+1.0ᴇ-60×𝐢)", "13.01+5ᴇ-50×𝐢", Radian, Cartesian, 4);
+  assert_expression_approximates_to<double>("atanh(0.99999999999+1.0ᴇ-150×𝐢)", "13.01+5ᴇ-140×𝐢", Radian, Cartesian, 4);
+  assert_expression_approximates_to<double>("atanh(0.99999999999+1.0ᴇ-250×𝐢)", "13.01+5ᴇ-240×𝐢", Radian, Cartesian, 4);
+  assert_expression_approximates_to<double>("atanh(0.99999999999+1.0ᴇ-300×𝐢)", "13.01+5ᴇ-290×𝐢", Radian, Cartesian, 4);
+
  // WARNING: evaluate on branch cut can be multivalued
  assert_expression_approximates_to<double>("acos(2)", "1.3169578969248×𝐢", Radian);
  assert_expression_approximates_to<double>("acos(2)", "75.456129290217×𝐢", Degree);
@@ -831,11 +844,12 @@ QUIZ_CASE(poincare_approximation_complex_format) {
  assert_expression_approximates_to<double>("√(-1)", "unreal", Radian, Real);
  assert_expression_approximates_to<double>("√(-1)×√(-1)", "unreal", Radian, Real);
  assert_expression_approximates_to<double>("ln(-2)", "unreal", Radian, Real);
-  assert_expression_approximates_to<double>("(-8)^(1/3)", "unreal", Radian, Real); // Power always approximates to the principal root (even if unreal)
+  // Power/Root approximates to the first REAL root in Real mode
+  assert_expression_simplifies_approximates_to<double>("(-8)^(1/3)", "-2", Radian, Real); // Power have to be simplified first in order to spot the right form c^(p/q) with p, q integers
  assert_expression_approximates_to<double>("root(-8,3)", "-2", Radian, Real); // Root approximates to the first REAL root in Real mode
  assert_expression_approximates_to<double>("8^(1/3)", "2", Radian, Real);
-  assert_expression_approximates_to<float>("(-8)^(2/3)", "unreal", Radian, Real); // Power always approximates to the principal root (even if unreal)
-  assert_expression_approximates_to<float>("root(-8, 3)^2", "4", Radian, Real); // Root approximates to the first REAL root in Real mode
+  assert_expression_simplifies_approximates_to<float>("(-8)^(2/3)", "4", Radian, Real); // Power have to be simplified first (cf previous comment)
+  assert_expression_approximates_to<float>("root(-8, 3)^2", "4", Radian, Real);
  assert_expression_approximates_to<double>("root(-8,3)", "-2", Radian, Real);

  // Cartesian
@@ -924,6 +938,19 @@ QUIZ_CASE(poincare_approximation_mix) {
  assert_expression_approximates_to<double>("sin(3)2(4+2)", "1.6934400967184", Radian);
  assert_expression_approximates_to<float>("4/2×(2+3)", "10");
  assert_expression_approximates_to<double>("4/2×(2+3)", "10");
+
+  assert_expression_simplifies_and_approximates_to("1.0092^(20)", "1.2010050593402");
+  assert_expression_simplifies_and_approximates_to("1.0092^(50)×ln(3/2)", "0.6409373488899", Degree, Cartesian, 13);
+  assert_expression_simplifies_and_approximates_to("1.0092^(50)×ln(1.0092)", "1.447637354655ᴇ-2", Degree, Cartesian, 13);
+  assert_expression_approximates_to<double>("1.0092^(20)", "1.2010050593402");
+  assert_expression_approximates_to<double>("1.0092^(50)×ln(3/2)", "0.6409373488899", Degree, Cartesian, 13);
+  assert_expression_approximates_to<double>("1.0092^(50)×ln(1.0092)", "1.447637354655ᴇ-2", Degree, Cartesian, 13);
+  assert_expression_simplifies_approximates_to<double>("1.0092^(20)", "1.2010050593402");
+  assert_expression_simplifies_approximates_to<double>("1.0092^(50)×ln(3/2)", "0.6409373488899", Degree, Cartesian, 13);
+  //assert_expression_approximates_to<float>("1.0092^(20)", "1.201005"); TODO does not work
+  assert_expression_approximates_to<float>("1.0092^(50)×ln(3/2)", "0.6409366");
+  //assert_expression_simplifies_approximates_to<float>("1.0092^(20)", "1.2010050593402"); TODO does not work
+  //assert_expression_simplifies_approximates_to<float>("1.0092^(50)×ln(3/2)", "6.4093734888993ᴇ-1"); TODO does not work
 }


--- a/poincare/test/helper.cpp
+++ b/poincare/test/helper.cpp
@@ -108,6 +108,26 @@ void assert_expression_approximates_to(const char * expression, const char * app
    }, numberOfDigits);
 }

+void assert_expression_simplifies_and_approximates_to(const char * expression, const char * approximation, Preferences::AngleUnit angleUnit, Preferences::ComplexFormat complexFormat, int numberOfSignificantDigits) {
+  int numberOfDigits = numberOfSignificantDigits > 0 ? numberOfSignificantDigits : PrintFloat::k_numberOfStoredSignificantDigits;
+  assert_parsed_expression_process_to(expression, approximation, SystemForApproximation, complexFormat, angleUnit, ReplaceAllSymbolsWithDefinitionsOrUndefined, [](Expression e, ExpressionNode::ReductionContext reductionContext) {
+      Expression reduced;
+      Expression approximated;
+      e.simplifyAndApproximate(&reduced, &approximated, reductionContext.context(), reductionContext.complexFormat(), reductionContext.angleUnit(), reductionContext.symbolicComputation());
+      return approximated;
+    }, numberOfDigits);
+}
+
+template<typename T>
+void assert_expression_simplifies_approximates_to(const char * expression, const char * approximation, Preferences::AngleUnit angleUnit, Preferences::ComplexFormat complexFormat, int numberOfSignificantDigits) {
+  int numberOfDigits = sizeof(T) == sizeof(double) ? PrintFloat::k_numberOfStoredSignificantDigits : PrintFloat::k_numberOfPrintedSignificantDigits;
+  numberOfDigits = numberOfSignificantDigits > 0 ? numberOfSignificantDigits : numberOfDigits;
+  assert_parsed_expression_process_to(expression, approximation, SystemForApproximation, complexFormat, angleUnit, ReplaceAllSymbolsWithDefinitionsOrUndefined, [](Expression e, ExpressionNode::ReductionContext reductionContext) {
+      e = e.simplify(reductionContext);
+      return e.approximate<T>(reductionContext.context(), reductionContext.complexFormat(), reductionContext.angleUnit());
+    }, numberOfDigits);
+}
+
 void assert_layout_serialize_to(Poincare::Layout layout, const char * serialization) {
  constexpr int bufferSize = 255;
  char buffer[bufferSize];
@@ -122,3 +142,5 @@ void assert_expression_layouts_as(Poincare::Expression expression, Poincare::Lay

 template void assert_expression_approximates_to<float>(char const*, char const *, Poincare::Preferences::AngleUnit, Poincare::Preferences::ComplexFormat, int);
 template void assert_expression_approximates_to<double>(char const*, char const *,  Poincare::Preferences::AngleUnit, Poincare::Preferences::ComplexFormat, int);
+template void assert_expression_simplifies_approximates_to<float>(char const*, char const *, Poincare::Preferences::AngleUnit, Poincare::Preferences::ComplexFormat, int);
+template void assert_expression_simplifies_approximates_to<double>(char const*, char const *,  Poincare::Preferences::AngleUnit, Poincare::Preferences::ComplexFormat, int);
--- a/poincare/test/helper.h
+++ b/poincare/test/helper.h
@@ -43,6 +43,9 @@ void assert_parsed_expression_simplify_to(const char * expression, const char *

 template<typename T>
 void assert_expression_approximates_to(const char * expression, const char * approximation, Poincare::Preferences::AngleUnit angleUnit = Degree, Poincare::Preferences::ComplexFormat complexFormat = Cartesian, int numberOfSignificantDigits = -1);
+void assert_expression_simplifies_and_approximates_to(const char * expression, const char * approximation, Poincare::Preferences::AngleUnit angleUnit = Degree, Poincare::Preferences::ComplexFormat complexFormat = Cartesian, int numberOfSignificantDigits = -1);
+template<typename T>
+void assert_expression_simplifies_approximates_to(const char * expression, const char * approximation, Poincare::Preferences::AngleUnit angleUnit = Degree, Poincare::Preferences::ComplexFormat complexFormat = Cartesian, int numberOfSignificantDigits = -1);


 // Layout serializing
--- a/poincare/test/simplification.cpp
+++ b/poincare/test/simplification.cpp
@@ -1231,8 +1231,8 @@ QUIZ_CASE(poincare_simplification_reduction_target) {
  assert_parsed_expression_simplify_to("(1+x)/(1+x)", "1", User);

  // Apply rule x^(2/3) --> root(x,3)^2 for ReductionTarget = System
-  assert_parsed_expression_simplify_to("x^(2/3)", "root(x,3)^2", SystemForApproximation);
-  assert_parsed_expression_simplify_to("x^(2/3)", "root(x,3)^2", SystemForAnalysis);
+  assert_parsed_expression_simplify_to("x^(2/3)", "x^\u00122/3\u0013", SystemForApproximation);
+  assert_parsed_expression_simplify_to("x^(2/3)", "x^\u00122/3\u0013", SystemForAnalysis);
  assert_parsed_expression_simplify_to("x^(2/3)", "x^\u00122/3\u0013", User);
  assert_parsed_expression_simplify_to("x^(1/3)", "root(x,3)", SystemForApproximation);
  assert_parsed_expression_simplify_to("x^(1/3)", "root(x,3)", SystemForAnalysis);