diff --git a/apps/probability/distribution/distribution.cpp b/apps/probability/distribution/distribution.cpp index ac0c8f517..aae3631d0 100644 --- 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; } diff --git a/apps/probability/test/distributions.cpp b/apps/probability/test/distributions.cpp index 4cb4b5c84..20bd9d9f9 100644 --- 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); } diff --git a/build/targets.all.mak b/build/targets.all.mak index 989c4fb7b..e6d989add 100644 --- 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 diff --git a/liba/include/math.h b/liba/include/math.h index 7fef24005..3cb95bcb5 100644 --- 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 diff --git a/liba/test/math.c b/liba/test/math.c index 64f874df8..c09fc30f2 100644 --- a/liba/test/math.c +++ b/liba/test/math.c @@ -1,210 +1,97 @@ // This file tests that each math fuction links. #include -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); } diff --git a/libaxx/include/cmath b/libaxx/include/cmath index e3c7835e2..6a3e3894e 100644 --- a/libaxx/include/cmath +++ b/libaxx/include/cmath @@ -3,11 +3,13 @@ #include +#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); } } diff --git a/poincare/include/poincare/approximation_helper.h b/poincare/include/poincare/approximation_helper.h index 169143faf..9b0b8ce17 100644 --- a/poincare/include/poincare/approximation_helper.h +++ b/poincare/include/poincare/approximation_helper.h @@ -10,7 +10,7 @@ namespace Poincare { namespace ApproximationHelper { template int PositiveIntegerApproximationIfPossible(const ExpressionNode * expression, bool * isUndefined, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit); - template std::complex TruncateRealOrImaginaryPartAccordingToArgument(std::complex c); + template std::complex NeglectRealOrImaginaryPartIfNeglectable(std::complex result, std::complex input1, std::complex input2 = 1.0); template using ComplexCompute = Complex(*)(const std::complex, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit); template Evaluation Map(const ExpressionNode * expression, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, ComplexCompute compute); diff --git a/poincare/include/poincare/power.h b/poincare/include/poincare/power.h index c5860e77e..77f2dda77 100644 --- 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 static Complex computeNotPrincipalRealRootOfRationalPow(const std::complex c, T p, T q); template static Complex compute(const std::complex c, const std::complex d, Preferences::ComplexFormat complexFormat); private: @@ -59,11 +60,12 @@ private: template static MatrixComplex computeOnMatrixAndComplex(const MatrixComplex m, const std::complex d, Preferences::ComplexFormat complexFormat); template static MatrixComplex computeOnMatrices(const MatrixComplex m, const MatrixComplex n, Preferences::ComplexFormat complexFormat); Evaluation approximate(SinglePrecision p, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const override { - return ApproximationHelper::MapReduce(this, context, complexFormat, angleUnit, compute, computeOnComplexAndMatrix, computeOnMatrixAndComplex, computeOnMatrices); + return templatedApproximate(context, complexFormat, angleUnit); } Evaluation approximate(DoublePrecision p, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const override { - return ApproximationHelper::MapReduce(this, context, complexFormat, angleUnit, compute, computeOnComplexAndMatrix, computeOnMatrixAndComplex, computeOnMatrices); + return templatedApproximate(context, complexFormat, angleUnit); } + template Evaluation templatedApproximate(Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit) const; }; class Power final : public Expression { diff --git a/poincare/include/poincare/trigonometry.h b/poincare/include/poincare/trigonometry.h index 4b0e3fd6e..0eac9d5c2 100644 --- 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 static std::complex ConvertToRadian(const std::complex c, Preferences::AngleUnit angleUnit); template static std::complex ConvertRadianToAngleUnit(const std::complex c, Preferences::AngleUnit angleUnit); - template static std::complex RoundToMeaningfulDigits(const std::complex result, const std::complex input); -private: - template static T RoundToMeaningfulDigits(T result, T input); }; } diff --git a/poincare/src/approximation_helper.cpp b/poincare/src/approximation_helper.cpp index 32f1f88ad..61c5f8d20 100644 --- a/poincare/src/approximation_helper.cpp +++ b/poincare/src/approximation_helper.cpp @@ -10,9 +10,19 @@ extern "C" { namespace Poincare { -template T absMod(T a, T b) { - T result = std::fmod(std::fabs(a), b); - return result > b/2 ? b-result : result; +template 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 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 int ApproximationHelper::PositiveIntegerApproximationIfPos return absInt((int)scalar); } -template std::complex ApproximationHelper::TruncateRealOrImaginaryPartAccordingToArgument(std::complex c) { - T arg = std::arg(c); - T precision = 10*Expression::Epsilon(); - if (absMod(arg, (T)M_PI) <= precision) { - c.imag(0); +template std::complex ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(std::complex result, std::complex input1, std::complex 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(); + if (isNegligeable(result.imag(), precision, magnitude1, magnitude2)) { + result.imag(0); } - if (absMod(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 Evaluation ApproximationHelper::Map(const ExpressionNode * expression, Context * context, Preferences::ComplexFormat complexFormat, Preferences::AngleUnit angleUnit, ComplexCompute compute) { @@ -106,8 +126,8 @@ template MatrixComplex ApproximationHelper::ElementWiseOnComplexM template int Poincare::ApproximationHelper::PositiveIntegerApproximationIfPossible(Poincare::ExpressionNode const*, bool*, Poincare::Context*, Poincare::Preferences::ComplexFormat, Poincare::Preferences::AngleUnit); template int Poincare::ApproximationHelper::PositiveIntegerApproximationIfPossible(Poincare::ExpressionNode const*, bool*, Poincare::Context*, Poincare::Preferences::ComplexFormat, Poincare::Preferences::AngleUnit); -template std::complex Poincare::ApproximationHelper::TruncateRealOrImaginaryPartAccordingToArgument(std::complex); -template std::complex Poincare::ApproximationHelper::TruncateRealOrImaginaryPartAccordingToArgument(std::complex); +template std::complex Poincare::ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(std::complex,std::complex,std::complex); +template std::complex Poincare::ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(std::complex,std::complex,std::complex); template Poincare::Evaluation Poincare::ApproximationHelper::Map(const Poincare::ExpressionNode * expression, Poincare::Context * context, Poincare::Preferences::ComplexFormat, Poincare::Preferences::AngleUnit angleUnit, Poincare::ApproximationHelper::ComplexCompute compute); template Poincare::Evaluation Poincare::ApproximationHelper::Map(const Poincare::ExpressionNode * expression, Poincare::Context * context, Poincare::Preferences::ComplexFormat, Poincare::Preferences::AngleUnit angleUnit, Poincare::ApproximationHelper::ComplexCompute compute); template Poincare::Evaluation Poincare::ApproximationHelper::MapReduce(const Poincare::ExpressionNode * expression, Poincare::Context * context, Poincare::Preferences::ComplexFormat, Poincare::Preferences::AngleUnit angleUnit, Poincare::ApproximationHelper::ComplexAndComplexReduction computeOnComplexes, Poincare::ApproximationHelper::ComplexAndMatrixReduction computeOnComplexAndMatrix, Poincare::ApproximationHelper::MatrixAndComplexReduction computeOnMatrixAndComplex, Poincare::ApproximationHelper::MatrixAndMatrixReduction computeOnMatrices); diff --git a/poincare/src/arc_cosine.cpp b/poincare/src/arc_cosine.cpp index 353c4b974..9f27e4f5d 100644 --- a/poincare/src/arc_cosine.cpp +++ b/poincare/src/arc_cosine.cpp @@ -44,7 +44,7 @@ Complex ArcCosineNode::computeOnComplex(const std::complex c, Preferences: result.imag(-result.imag()); // other side of the cut } } - result = Trigonometry::RoundToMeaningfulDigits(result, c); + result = ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(result, c); return Complex::Builder(Trigonometry::ConvertRadianToAngleUnit(result, angleUnit)); } diff --git a/poincare/src/arc_sine.cpp b/poincare/src/arc_sine.cpp index 2ee8b0fd1..4d7100ac8 100644 --- a/poincare/src/arc_sine.cpp +++ b/poincare/src/arc_sine.cpp @@ -44,7 +44,7 @@ Complex ArcSineNode::computeOnComplex(const std::complex c, Preferences::C result.imag(-result.imag()); // other side of the cut } } - result = Trigonometry::RoundToMeaningfulDigits(result, c); + result = ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(result, c); return Complex::Builder(Trigonometry::ConvertRadianToAngleUnit(result, angleUnit)); } diff --git a/poincare/src/arc_tangent.cpp b/poincare/src/arc_tangent.cpp index 2f8152411..19bb310bb 100644 --- a/poincare/src/arc_tangent.cpp +++ b/poincare/src/arc_tangent.cpp @@ -40,7 +40,7 @@ Complex ArcTangentNode::computeOnComplex(const std::complex c, Preferences result.real(-result.real()); // other side of the cut } } - result = Trigonometry::RoundToMeaningfulDigits(result, c); + result = ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(result, c); return Complex::Builder(Trigonometry::ConvertRadianToAngleUnit(result, angleUnit)); } diff --git a/poincare/src/arithmetic.cpp b/poincare/src/arithmetic.cpp index f4563db79..a7b50149b 100644 --- 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); diff --git a/poincare/src/cosine.cpp b/poincare/src/cosine.cpp index 91414e60c..2e216ce8e 100644 --- a/poincare/src/cosine.cpp +++ b/poincare/src/cosine.cpp @@ -19,7 +19,7 @@ template Complex CosineNode::computeOnComplex(const std::complex c, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit) { std::complex angleInput = Trigonometry::ConvertToRadian(c, angleUnit); std::complex res = std::cos(angleInput); - return Complex::Builder(Trigonometry::RoundToMeaningfulDigits(res, angleInput)); + return Complex::Builder(ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(res, angleInput)); } Layout CosineNode::createLayout(Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const { diff --git a/poincare/src/hyperbolic_arc_cosine.cpp b/poincare/src/hyperbolic_arc_cosine.cpp index 8ea38dd6c..eb6f4b668 100644 --- a/poincare/src/hyperbolic_arc_cosine.cpp +++ b/poincare/src/hyperbolic_arc_cosine.cpp @@ -23,7 +23,7 @@ Complex HyperbolicArcCosineNode::computeOnComplex(const std::complex 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::Builder(Trigonometry::RoundToMeaningfulDigits(result, c)); + return Complex::Builder(ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(result, c)); } template Complex Poincare::HyperbolicArcCosineNode::computeOnComplex(std::complex, Preferences::ComplexFormat, Preferences::AngleUnit); diff --git a/poincare/src/hyperbolic_arc_sine.cpp b/poincare/src/hyperbolic_arc_sine.cpp index 7ee51048c..3952a4f29 100644 --- a/poincare/src/hyperbolic_arc_sine.cpp +++ b/poincare/src/hyperbolic_arc_sine.cpp @@ -27,7 +27,7 @@ Complex HyperbolicArcSineNode::computeOnComplex(const std::complex c, Pref if (c.real() == 0 && c.imag() < 1) { result.real(-result.real()); // other side of the cut } - return Complex::Builder(Trigonometry::RoundToMeaningfulDigits(result, c)); + return Complex::Builder(ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(result, c)); } template Complex Poincare::HyperbolicArcSineNode::computeOnComplex(std::complex, Preferences::ComplexFormat, Preferences::AngleUnit); diff --git a/poincare/src/hyperbolic_arc_tangent.cpp b/poincare/src/hyperbolic_arc_tangent.cpp index 13c63bba5..09453d22a 100644 --- a/poincare/src/hyperbolic_arc_tangent.cpp +++ b/poincare/src/hyperbolic_arc_tangent.cpp @@ -27,7 +27,7 @@ Complex HyperbolicArcTangentNode::computeOnComplex(const std::complex c, P if (c.imag() == 0 && c.real() > 1) { result.imag(-result.imag()); // other side of the cut } - return Complex::Builder(Trigonometry::RoundToMeaningfulDigits(result, c)); + return Complex::Builder(ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(result, c)); } template Complex Poincare::HyperbolicArcTangentNode::computeOnComplex(std::complex, Preferences::ComplexFormat, Preferences::AngleUnit); diff --git a/poincare/src/hyperbolic_cosine.cpp b/poincare/src/hyperbolic_cosine.cpp index c722ec209..c963ee4ad 100644 --- 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 Complex HyperbolicCosineNode::computeOnComplex(const std::complex c, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit) { - return Complex::Builder(Trigonometry::RoundToMeaningfulDigits(std::cosh(c), c)); + return Complex::Builder(ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(std::cosh(c), c)); } template Complex Poincare::HyperbolicCosineNode::computeOnComplex(std::complex, Preferences::ComplexFormat, Preferences::AngleUnit); diff --git a/poincare/src/hyperbolic_sine.cpp b/poincare/src/hyperbolic_sine.cpp index 232bd8fcb..237804747 100644 --- 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 Complex HyperbolicSineNode::computeOnComplex(const std::complex c, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit) { - return Complex::Builder(Trigonometry::RoundToMeaningfulDigits(std::sinh(c), c)); + return Complex::Builder(ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(std::sinh(c), c)); } template Complex Poincare::HyperbolicSineNode::computeOnComplex(std::complex, Preferences::ComplexFormat, Preferences::AngleUnit); diff --git a/poincare/src/hyperbolic_tangent.cpp b/poincare/src/hyperbolic_tangent.cpp index 7c4ad545d..48a5e8830 100644 --- 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 Complex HyperbolicTangentNode::computeOnComplex(const std::complex c, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit) { - return Complex::Builder(Trigonometry::RoundToMeaningfulDigits(std::tanh(c), c)); + return Complex::Builder(ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(std::tanh(c), c)); } template Complex Poincare::HyperbolicTangentNode::computeOnComplex(std::complex, Preferences::ComplexFormat, Preferences::AngleUnit); diff --git a/poincare/src/multiplication.cpp b/poincare/src/multiplication.cpp index 5513bf17e..f0767a8c5 100644 --- 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(); Number m = Number::Multiplication(o0, static_cast(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(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 diff --git a/poincare/src/nth_root.cpp b/poincare/src/nth_root.cpp index 74ee32f52..db4748a47 100644 --- a/poincare/src/nth_root.cpp +++ b/poincare/src/nth_root.cpp @@ -46,18 +46,12 @@ Evaluation 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 absBasec = basec; - absBasec.real(std::fabs(absBasec.real())); - // compute root(|x|, q) - Complex absBasePowIndex = PowerNode::compute(absBasec, std::complex(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::Builder(-absBasePowIndex.stdComplex()) : absBasePowIndex; - } - } + Complex result = PowerNode::computeNotPrincipalRealRootOfRationalPow(basec, (T)1.0, indexc.real()); + if (!result.isUndefined()) { + return result; + } } result = PowerNode::compute(basec, std::complex(1.0)/(indexc), complexFormat); } diff --git a/poincare/src/power.cpp b/poincare/src/power.cpp index 495aa697a..c7362383e 100644 --- a/poincare/src/power.cpp +++ b/poincare/src/power.cpp @@ -127,6 +127,32 @@ bool PowerNode::childAtIndexNeedsUserParentheses(const Expression & child, int c // Private +template +Complex PowerNode::computeNotPrincipalRealRootOfRationalPow(const std::complex 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 absc = c; + absc.real(std::fabs(absc.real())); + // compute |c|^(p/q) which is a real + Complex absCPowD = PowerNode::compute(absc, std::complex(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::Builder(-absCPowD.stdComplex()) : absCPowD; + } + return Complex::Undefined(); +} + template Complex PowerNode::compute(const std::complex c, const std::complex d, Preferences::ComplexFormat complexFormat) { std::complex result; @@ -151,8 +177,13 @@ Complex PowerNode::compute(const std::complex c, const std::complex 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::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::Builder(ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(result, c, d)); } // Layout @@ -257,6 +288,50 @@ template MatrixComplex PowerNode::computeOnMatrices(const MatrixC return MatrixComplex::Undefined(); } +template Evaluation 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 base = childAtIndex(0)->approximate(T(), context, complexFormat, angleUnit); + if (base.type() != EvaluationNode::Type::Complex) { + goto defaultApproximation; + } + std::complex c = static_cast &>(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(childAtIndex(1)); + p = r->signedNumerator().approximate(); + q = r->denominator().approximate(); + } + /* 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(childAtIndex(1)->childAtIndex(0)); + const RationalNode * qRat = static_cast(childAtIndex(1)->childAtIndex(1)); + if (!pRat->denominator().isOne() || !qRat->denominator().isOne()) { + goto defaultApproximation; + } + p = pRat->signedNumerator().approximate(); + q = qRat->signedNumerator().approximate(); + } + /* 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 result = computeNotPrincipalRealRootOfRationalPow(c, p, q); + if (!result.isUndefined()) { + return result; + } + } +defaultApproximation: + return ApproximationHelper::MapReduce(this, context, complexFormat, angleUnit, compute, computeOnComplexAndMatrix, computeOnMatrixAndComplex, computeOnMatrices); +} + // 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().signedIntegerNumerator(); - Integer q = childAtIndex(1).convert().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::Builder(1.0)); diff --git a/poincare/src/sine.cpp b/poincare/src/sine.cpp index a4701e0c6..3668d8a3a 100644 --- a/poincare/src/sine.cpp +++ b/poincare/src/sine.cpp @@ -19,7 +19,7 @@ template Complex SineNode::computeOnComplex(const std::complex c, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit) { std::complex angleInput = Trigonometry::ConvertToRadian(c, angleUnit); std::complex res = std::sin(angleInput); - return Complex::Builder(Trigonometry::RoundToMeaningfulDigits(res, angleInput)); + return Complex::Builder(ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(res, angleInput)); } Layout SineNode::createLayout(Preferences::PrintFloatMode floatDisplayMode, int numberOfSignificantDigits) const { diff --git a/poincare/src/square_root.cpp b/poincare/src/square_root.cpp index ed0954476..0c8361d1e 100644 --- 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 Complex SquareRootNode::computeOnComplex(const std::complex c, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit) { std::complex 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::Builder(ApproximationHelper::TruncateRealOrImaginaryPartAccordingToArgument(result)); + return Complex::Builder(ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(result, std::complex(std::log(std::abs(c)), std::arg(c)))); } Expression SquareRootNode::shallowReduce(ReductionContext reductionContext) { diff --git a/poincare/src/tangent.cpp b/poincare/src/tangent.cpp index 0d44b0a5a..e140a4992 100644 --- a/poincare/src/tangent.cpp +++ b/poincare/src/tangent.cpp @@ -30,7 +30,7 @@ template Complex TangentNode::computeOnComplex(const std::complex c, Preferences::ComplexFormat, Preferences::AngleUnit angleUnit) { std::complex angleInput = Trigonometry::ConvertToRadian(c, angleUnit); std::complex res = std::tan(angleInput); - return Complex::Builder(Trigonometry::RoundToMeaningfulDigits(res, angleInput)); + return Complex::Builder(ApproximationHelper::NeglectRealOrImaginaryPartIfNeglectable(res, angleInput)); } Expression TangentNode::shallowReduce(ReductionContext reductionContext) { diff --git a/poincare/src/trigonometry.cpp b/poincare/src/trigonometry.cpp index dd802f97b..168f7e49d 100644 --- a/poincare/src/trigonometry.cpp +++ b/poincare/src/trigonometry.cpp @@ -412,34 +412,9 @@ std::complex Trigonometry::ConvertRadianToAngleUnit(const std::complex c, return c; } -template -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 precision = 10*Expression::Epsilon(); - result = std::round(result/precision)*precision; - } - return result; -} - -template -std::complex Trigonometry::RoundToMeaningfulDigits(const std::complex result, const std::complex input) { - return std::complex(RoundToMeaningfulDigits(result.real(), std::abs(input)), RoundToMeaningfulDigits(result.imag(), std::abs(input))); -} - template std::complex Trigonometry::ConvertToRadian(std::complex, Preferences::AngleUnit); template std::complex Trigonometry::ConvertToRadian(std::complex, Preferences::AngleUnit); template std::complex Trigonometry::ConvertRadianToAngleUnit(std::complex, Preferences::AngleUnit); template std::complex Trigonometry::ConvertRadianToAngleUnit(std::complex, Preferences::AngleUnit); -template std::complex Trigonometry::RoundToMeaningfulDigits(std::complex, std::complex); -template std::complex Trigonometry::RoundToMeaningfulDigits(std::complex, std::complex); } diff --git a/poincare/test/approximation.cpp b/poincare/test/approximation.cpp index 89bc3b8ac..723027d6e 100644 --- a/poincare/test/approximation.cpp +++ b/poincare/test/approximation.cpp @@ -149,6 +149,12 @@ QUIZ_CASE(poincare_approximation_power) { assert_expression_approximates_to_scalar("2^3", 8.0f); assert_expression_approximates_to_scalar("(3+𝐢)^(4+𝐢)", NAN); assert_expression_approximates_to_scalar("[[1,2][3,4]]^2", NAN); + + + assert_expression_approximates_to("(-10)^0.00000001", "unreal", Radian, Real); + assert_expression_approximates_to("(-10)^0.00000001", "1+3.141593ᴇ-8×𝐢", Radian, Cartesian); + assert_expression_simplifies_approximates_to("3.5^2.0000001", "12.25"); + assert_expression_simplifies_approximates_to("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("atanh(𝐢-4)", "-0.238878+1.50862×𝐢", Radian, Cartesian, 6); assert_expression_approximates_to("atanh(𝐢-4)", "-0.238878+1.50862×𝐢", Degree, Cartesian, 6); + // Check that the complex part is not neglected + assert_expression_approximates_to("atanh(0.99999999999+1.0ᴇ-26×𝐢)", "13.01+5ᴇ-16×𝐢", Radian, Cartesian, 4); + assert_expression_approximates_to("atanh(0.99999999999+1.0ᴇ-60×𝐢)", "13.01+5ᴇ-50×𝐢", Radian, Cartesian, 4); + assert_expression_approximates_to("atanh(0.99999999999+1.0ᴇ-150×𝐢)", "13.01+5ᴇ-140×𝐢", Radian, Cartesian, 4); + assert_expression_approximates_to("atanh(0.99999999999+1.0ᴇ-250×𝐢)", "13.01+5ᴇ-240×𝐢", Radian, Cartesian, 4); + assert_expression_approximates_to("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("acos(2)", "1.3169578969248×𝐢", Radian); assert_expression_approximates_to("acos(2)", "75.456129290217×𝐢", Degree); @@ -831,11 +844,12 @@ QUIZ_CASE(poincare_approximation_complex_format) { assert_expression_approximates_to("√(-1)", "unreal", Radian, Real); assert_expression_approximates_to("√(-1)×√(-1)", "unreal", Radian, Real); assert_expression_approximates_to("ln(-2)", "unreal", Radian, Real); - assert_expression_approximates_to("(-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("(-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("root(-8,3)", "-2", Radian, Real); // Root approximates to the first REAL root in Real mode assert_expression_approximates_to("8^(1/3)", "2", Radian, Real); - assert_expression_approximates_to("(-8)^(2/3)", "unreal", Radian, Real); // Power always approximates to the principal root (even if unreal) - assert_expression_approximates_to("root(-8, 3)^2", "4", Radian, Real); // Root approximates to the first REAL root in Real mode + assert_expression_simplifies_approximates_to("(-8)^(2/3)", "4", Radian, Real); // Power have to be simplified first (cf previous comment) + assert_expression_approximates_to("root(-8, 3)^2", "4", Radian, Real); assert_expression_approximates_to("root(-8,3)", "-2", Radian, Real); // Cartesian @@ -924,6 +938,19 @@ QUIZ_CASE(poincare_approximation_mix) { assert_expression_approximates_to("sin(3)2(4+2)", "1.6934400967184", Radian); assert_expression_approximates_to("4/2×(2+3)", "10"); assert_expression_approximates_to("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("1.0092^(20)", "1.2010050593402"); + assert_expression_approximates_to("1.0092^(50)×ln(3/2)", "0.6409373488899", Degree, Cartesian, 13); + assert_expression_approximates_to("1.0092^(50)×ln(1.0092)", "1.447637354655ᴇ-2", Degree, Cartesian, 13); + assert_expression_simplifies_approximates_to("1.0092^(20)", "1.2010050593402"); + assert_expression_simplifies_approximates_to("1.0092^(50)×ln(3/2)", "0.6409373488899", Degree, Cartesian, 13); + //assert_expression_approximates_to("1.0092^(20)", "1.201005"); TODO does not work + assert_expression_approximates_to("1.0092^(50)×ln(3/2)", "0.6409366"); + //assert_expression_simplifies_approximates_to("1.0092^(20)", "1.2010050593402"); TODO does not work + //assert_expression_simplifies_approximates_to("1.0092^(50)×ln(3/2)", "6.4093734888993ᴇ-1"); TODO does not work } diff --git a/poincare/test/helper.cpp b/poincare/test/helper.cpp index ecd631b65..67d035b41 100644 --- 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 +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(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(char const*, char const *, Poincare::Preferences::AngleUnit, Poincare::Preferences::ComplexFormat, int); template void assert_expression_approximates_to(char const*, char const *, Poincare::Preferences::AngleUnit, Poincare::Preferences::ComplexFormat, int); +template void assert_expression_simplifies_approximates_to(char const*, char const *, Poincare::Preferences::AngleUnit, Poincare::Preferences::ComplexFormat, int); +template void assert_expression_simplifies_approximates_to(char const*, char const *, Poincare::Preferences::AngleUnit, Poincare::Preferences::ComplexFormat, int); diff --git a/poincare/test/helper.h b/poincare/test/helper.h index cdac8d201..5e7f5942e 100644 --- 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 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 +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 diff --git a/poincare/test/simplification.cpp b/poincare/test/simplification.cpp index 0dd49cc10..cd5bd1c04 100644 --- 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);