From 1995781a9f859d84899322984891d199285f2f42 Mon Sep 17 00:00:00 2001 From: Arthur Camouseigt Date: Tue, 4 Aug 2020 10:22:35 +0200 Subject: [PATCH] [Poincare/IEEE754] Changed methods to use std function Implemented std::nextafter to replace hand made one. Methods next and previous are no longer making the difference between -0 and +0 Change-Id: I42e1a073623b70656d9df954694803840cf3088c --- liba/Makefile | 2 + liba/include/math.h | 4 ++ liba/src/external/openbsd/s_nextafter.c | 71 ++++++++++++++++++++++++ liba/src/external/openbsd/s_nextafterf.c | 62 +++++++++++++++++++++ libaxx/include/cmath | 4 ++ poincare/Makefile | 1 - poincare/include/poincare/ieee754.h | 28 ---------- poincare/src/solver.cpp | 8 ++- poincare/test/ieee754.cpp | 40 ------------- 9 files changed, 149 insertions(+), 71 deletions(-) create mode 100644 liba/src/external/openbsd/s_nextafter.c create mode 100644 liba/src/external/openbsd/s_nextafterf.c delete mode 100644 poincare/test/ieee754.cpp diff --git a/liba/Makefile b/liba/Makefile index 3639ceb91..639bf196d 100644 --- a/liba/Makefile +++ b/liba/Makefile @@ -70,6 +70,8 @@ liba_src += $(addprefix liba/src/external/openbsd/, \ s_logbf.c \ s_modf.c \ s_modff.c \ + s_nextafter.c \ + s_nextafterf.c \ s_rint.c \ s_roundf.c \ s_scalbnf.c \ diff --git a/liba/include/math.h b/liba/include/math.h index 3cb95bcb5..4c05d5bcb 100644 --- a/liba/include/math.h +++ b/liba/include/math.h @@ -70,6 +70,7 @@ float log10f(float x); float logf(float x); float modff(float x, float *iptr); float nearbyintf(float x); +float nextafterf(float from, float to); float powf(float x, float y); float roundf(float x); float scalbnf(float x, int exp); @@ -110,6 +111,7 @@ double log2(double x); double logb(double x); double modf(double x, double *iptr); double nearbyint(double x); +double nextafter(double from, double to); double pow(double x, double y); double rint(double x); double round(double x); @@ -156,6 +158,7 @@ extern int signgam; #define modff(x, iptr) __builtin_modff(x, iptr) #define nanf(tagp) __builtin_nanf(tagp) #define nearbyintf(x) __builtin_nearbyintf(x) +#define nextafterf(from, to) __builtin_nextafterf(from, to) #define powf(x, y) __builtin_powf(x, y) #define roundf(x) __builtin_roundf(x) #define scalbnf(x, exp) __builtin_scalbnf(x, exp) @@ -196,6 +199,7 @@ extern int signgam; #define modf(x, iptr) __builtin_modf(x, iptr) #define nan(tagp) __builtin_nan(tagp) #define nearbyint(x) __builtin_nearbyint(x) +#define nextafter(from, to) __builtin_nextafter(from, to) #define pow(x, y) __builtin_pow(x, y) #define rint(x) __builtin_rint(x) #define round(x) __builtin_round(x) diff --git a/liba/src/external/openbsd/s_nextafter.c b/liba/src/external/openbsd/s_nextafter.c new file mode 100644 index 000000000..199926cd8 --- /dev/null +++ b/liba/src/external/openbsd/s_nextafter.c @@ -0,0 +1,71 @@ +/* @(#)s_nextafter.c 5.1 93/09/24 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* IEEE functions + * nextafter(x,y) + * return the next machine floating-point number of x in the + * direction toward y. + * Special cases: + */ + +#include "math.h" +#include "math_private.h" + +double +nextafter(double x, double y) +{ + int32_t hx,hy,ix,iy; + u_int32_t lx,ly; + + EXTRACT_WORDS(hx,lx,x); + EXTRACT_WORDS(hy,ly,y); + ix = hx&0x7fffffff; /* |x| */ + iy = hy&0x7fffffff; /* |y| */ + + if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) || /* x is nan */ + ((iy>=0x7ff00000)&&((iy-0x7ff00000)|ly)!=0)) /* y is nan */ + return x+y; + if(x==y) return x; /* x=y, return x */ + if((ix|lx)==0) { /* x == 0 */ + INSERT_WORDS(x,hy&0x80000000,1); /* return +-minsubnormal */ + y = x*x; + if(y==x) return y; else return x; /* raise underflow flag */ + } + if(hx>=0) { /* x > 0 */ + if(hx>hy||((hx==hy)&&(lx>ly))) { /* x > y, x -= ulp */ + if(lx==0) hx -= 1; + lx -= 1; + } else { /* x < y, x += ulp */ + lx += 1; + if(lx==0) hx += 1; + } + } else { /* x < 0 */ + if(hy>=0||hx>hy||((hx==hy)&&(lx>ly))){/* x < y, x -= ulp */ + if(lx==0) hx -= 1; + lx -= 1; + } else { /* x > y, x += ulp */ + lx += 1; + if(lx==0) hx += 1; + } + } + hy = hx&0x7ff00000; + if(hy>=0x7ff00000) return x+x; /* overflow */ + if(hy<0x00100000) { /* underflow */ + y = x*x; + if(y!=x) { /* raise underflow flag */ + INSERT_WORDS(y,hx,lx); + return y; + } + } + INSERT_WORDS(x,hx,lx); + return x; +} diff --git a/liba/src/external/openbsd/s_nextafterf.c b/liba/src/external/openbsd/s_nextafterf.c new file mode 100644 index 000000000..e5c3d9f9f --- /dev/null +++ b/liba/src/external/openbsd/s_nextafterf.c @@ -0,0 +1,62 @@ +/* s_nextafterf.c -- float version of s_nextafter.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include "math.h" +#include "math_private.h" + +float +nextafterf(float x, float y) +{ + int32_t hx,hy,ix,iy; + + GET_FLOAT_WORD(hx,x); + GET_FLOAT_WORD(hy,y); + ix = hx&0x7fffffff; /* |x| */ + iy = hy&0x7fffffff; /* |y| */ + + if((ix>0x7f800000) || /* x is nan */ + (iy>0x7f800000)) /* y is nan */ + return x+y; + if(x==y) return x; /* x=y, return x */ + if(ix==0) { /* x == 0 */ + SET_FLOAT_WORD(x,(hy&0x80000000)|1);/* return +-minsubnormal */ + y = x*x; + if(y==x) return y; else return x; /* raise underflow flag */ + } + if(hx>=0) { /* x > 0 */ + if(hx>hy) { /* x > y, x -= ulp */ + hx -= 1; + } else { /* x < y, x += ulp */ + hx += 1; + } + } else { /* x < 0 */ + if(hy>=0||hx>hy){ /* x < y, x -= ulp */ + hx -= 1; + } else { /* x > y, x += ulp */ + hx += 1; + } + } + hy = hx&0x7f800000; + if(hy>=0x7f800000) return x+x; /* overflow */ + if(hy<0x00800000) { /* underflow */ + y = x*x; + if(y!=x) { /* raise underflow flag */ + SET_FLOAT_WORD(y,hx); + return y; + } + } + SET_FLOAT_WORD(x,hx); + return x; +} \ No newline at end of file diff --git a/libaxx/include/cmath b/libaxx/include/cmath index 9fdad6b7a..cc42319e7 100644 --- a/libaxx/include/cmath +++ b/libaxx/include/cmath @@ -39,6 +39,7 @@ #undef modff #undef nanf #undef nearbyintf +#undef nextafter #undef powf #undef roundf #undef scalbnf @@ -81,6 +82,7 @@ #undef modf #undef nan #undef nearbyint +#undef nextafter #undef pow #undef rint #undef round @@ -134,6 +136,7 @@ inline constexpr float logb(float x) { return __builtin_logbf(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 nextafter(float from, float to) { return __builtin_nextafterf(from, to); } 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); } @@ -183,6 +186,7 @@ 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 nextafter(double from, double to) { return __builtin_nextafter(from, to); } 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); } diff --git a/poincare/Makefile b/poincare/Makefile index 7ffd237cb..c289935cc 100644 --- a/poincare/Makefile +++ b/poincare/Makefile @@ -178,7 +178,6 @@ tests_src += $(addprefix poincare/test/,\ function_solver.cpp\ helper.cpp\ helpers.cpp\ - ieee754.cpp\ integer.cpp\ layout.cpp\ layout_cursor.cpp\ diff --git a/poincare/include/poincare/ieee754.h b/poincare/include/poincare/ieee754.h index c1b3c0bcc..accda9d33 100644 --- a/poincare/include/poincare/ieee754.h +++ b/poincare/include/poincare/ieee754.h @@ -73,40 +73,12 @@ public: } return exponentBase10; } - static T next(T f) { - return nextOrPrevious(f, true); - } - - static T previous(T f) { - return nextOrPrevious(f, false); - } private: union uint_float { uint64_t ui; T f; }; - static T nextOrPrevious(T f, bool isNext) { - if (std::isinf(f) || std::isnan(f)) { - return f; - } - uint_float u; - u.ui = 0; - u.f = f; - uint64_t oneBitOnSignBit = (uint64_t)1 << (k_exponentNbBits + k_mantissaNbBits); - if ((isNext && (u.ui & oneBitOnSignBit) > 0) // next: Negative float - || (!isNext && (u.ui & oneBitOnSignBit) == 0)) { // previous: Positive float - if ((isNext && u.ui == oneBitOnSignBit) // next: -0.0 - || (!isNext && u.ui == 0.0)) { // previous: 0.0 - u.ui = isNext ? 0 : oneBitOnSignBit; - } else { - u.ui -= 1; - } - } else { // next: Positive float, previous: Negative float - u.ui += 1; - } - return u.f; - } constexpr static size_t k_signNbBits = 1; constexpr static size_t k_exponentNbBits = sizeof(T) == sizeof(float) ? 8 : 11; diff --git a/poincare/src/solver.cpp b/poincare/src/solver.cpp index a5ea72bc3..37d8ec976 100644 --- a/poincare/src/solver.cpp +++ b/poincare/src/solver.cpp @@ -196,10 +196,14 @@ Coordinate2D Solver::IncreasingFunctionRoot(double ax, double bx, double * representable double between min and max strictly. If there is, we choose * it instead, otherwise, we reached the most precise result possible. */ if (currentAbscissa == min) { - currentAbscissa = IEEE754::next(min); + if (currentAbscissa != -INFINITY) { + currentAbscissa = std::nextafter(currentAbscissa, INFINITY); + } } if (currentAbscissa == max) { - currentAbscissa = IEEE754::previous(max); + if (currentAbscissa != INFINITY) { + currentAbscissa = std::nextafter(currentAbscissa, -INFINITY); + } } if (currentAbscissa == min || currentAbscissa == max) { break; diff --git a/poincare/test/ieee754.cpp b/poincare/test/ieee754.cpp deleted file mode 100644 index 85bdee9d4..000000000 --- a/poincare/test/ieee754.cpp +++ /dev/null @@ -1,40 +0,0 @@ -#include -#include -#include -#include -#include "helper.h" - -using namespace Poincare; - -template -void assert_next_and_previous_IEEE754_is(T a, T b) { - T next = IEEE754::next(a); - T previous = IEEE754::previous(b); - quiz_assert((std::isnan(next) && std::isnan(b)) || next == b); - quiz_assert((std::isnan(previous) && std::isnan(a)) || previous == a); -} - -QUIZ_CASE(ieee754_next_and_previous) { - assert_next_and_previous_IEEE754_is(0.0f, 1.4E-45f); - assert_next_and_previous_IEEE754_is(INFINITY, INFINITY); - assert_next_and_previous_IEEE754_is(NAN, NAN); - assert_next_and_previous_IEEE754_is(2.2566837E-10f, 2.2566839E-10f); - assert_next_and_previous_IEEE754_is(-INFINITY, -INFINITY); - assert_next_and_previous_IEEE754_is(-0.0f, 0.0f); - assert_next_and_previous_IEEE754_is(-1.4E-45f, -0.0f); - assert_next_and_previous_IEEE754_is(-3.4359738E10f, -3.43597363E10f); - quiz_assert(IEEE754::next(3.4028235E38f) == INFINITY); - quiz_assert(IEEE754::previous(-3.4028235E38f) == -INFINITY); - - assert_next_and_previous_IEEE754_is(0.0f, 4.94065645841246544176568792868E-324); - assert_next_and_previous_IEEE754_is(INFINITY, INFINITY); - assert_next_and_previous_IEEE754_is(NAN, NAN); - assert_next_and_previous_IEEE754_is(1.936766735060658315512927142E-282, 1.93676673506065869777770528092E-282); - assert_next_and_previous_IEEE754_is(-INFINITY, -INFINITY); - assert_next_and_previous_IEEE754_is(-0.0, 0.0); - assert_next_and_previous_IEEE754_is(-4.94065645841246544176568792868E-324, -0.0); - assert_next_and_previous_IEEE754_is(-1.38737906372912431085182213247E201, -1.38737906372912403890916981028E201); - quiz_assert(IEEE754::next(1.79769313486231570814527423732E308) == INFINITY); - quiz_assert(IEEE754::previous(-1.79769313486231570814527423732E308) == -INFINITY); - -}