Bug#153022: [powerpc libm] exp() in directed rounding modes gives wrong results
Jonathan Nieder wrote:
> debdiff attached.
Better debdiff attached.
diff -u eglibc-2.13/debian/changelog eglibc-2.13/debian/changelog
--- eglibc-2.13/debian/changelog
+++ eglibc-2.13/debian/changelog
@@ -1,3 +1,17 @@
+eglibc (2.13-38+test) local; urgency=low
+
+ * patches/any/cvs-exp-rounding-mode.diff: __ieee754_exp: save and
+ restore rounding mode and use round-to-nearest for all computations.
+ Closes: #153022
+ * patches/any/cvs-sin-accuracy.diff: fix incorrect rounding of sin.
+ * patches/any/cvs-sin-cos-tan-rounding.diff: __sin, __cos, tan: save
+ and restore rounding mode and use round-to-nearest for all
+ computations.
+ * patches/any/cvs-pow-rounding-mode.diff: __ieee754_pow: save and
+ restore rounding mode and use round-to-nearest for all computations.
+
+ -- Jonathan Nieder <jrnieder@gmail.com> Sun, 03 Mar 2013 23:41:19 -0800
+
eglibc (2.13-38) unstable; urgency=low
[ Adam Conrad ]
diff -u eglibc-2.13/debian/patches/series eglibc-2.13/debian/patches/series
--- eglibc-2.13/debian/patches/series
+++ eglibc-2.13/debian/patches/series
@@ -376,0 +377,4 @@
+any/cvs-exp-rounding-mode.diff
+any/cvs-sin-accuracy.diff
+any/cvs-sin-cos-tan-rounding.diff
+any/cvs-pow-rounding-mode.diff
only in patch2:
unchanged:
--- eglibc-2.13.orig/debian/patches/any/cvs-pow-rounding-mode.diff
+++ eglibc-2.13/debian/patches/any/cvs-pow-rounding-mode.diff
@@ -0,0 +1,329 @@
+2012-03-05 Joseph Myers <joseph@codesourcery.com>
+
+ [BZ #3976]
+ * sysdeps/ieee754/dbl-64/e_pow.c: Include <fenv.h>.
+ (__ieee754_pow): Save and restore rounding mode and use
+ round-to-nearest for main computations.
+ * math/libm-test.inc (pow_test_tonearest): New function.
+ (pow_test_towardzero): Likewise.
+ (pow_test_downward): Likewise.
+ (pow_test_upward): Likewise.
+ (main): Call the new functions.
+ * sysdeps/i386/fpu/libm-test-ulps: Update.
+ * sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
+---
+ math/libm-test.inc | 109 ++++++++++++++++++++++++++++++++++++++
+ sysdeps/i386/fpu/libm-test-ulps | 60 +++++++++++++++++++++
+ sysdeps/ieee754/dbl-64/e_pow.c | 14 ++++-
+ sysdeps/x86_64/fpu/libm-test-ulps | 48 +++++++++++++++++
+ 4 files changed, 229 insertions(+), 2 deletions(-)
+
+--- a/math/libm-test.inc
++++ b/math/libm-test.inc
+@@ -5175,6 +5175,111 @@ pow_test (void)
+ END (pow);
+ }
+
++
++static void
++pow_test_tonearest (void)
++{
++ int save_round_mode;
++ errno = 0;
++ FUNC(pow) (0, 0);
++ if (errno == ENOSYS)
++ /* Function not implemented. */
++ return;
++
++ START (pow_tonearest);
++
++ save_round_mode = fegetround ();
++
++ if (!fesetround (FE_TONEAREST))
++ {
++ TEST_ff_f (pow, 1.0625L, 1.125L, 1.070582293028761362162622578677070098674L);
++ TEST_ff_f (pow, 1.5L, 1.03125L, 1.519127098714743184071644334163037684948L);
++ }
++
++ fesetround (save_round_mode);
++
++ END (pow_tonearest);
++}
++
++
++static void
++pow_test_towardzero (void)
++{
++ int save_round_mode;
++ errno = 0;
++ FUNC(pow) (0, 0);
++ if (errno == ENOSYS)
++ /* Function not implemented. */
++ return;
++
++ START (pow_towardzero);
++
++ save_round_mode = fegetround ();
++
++ if (!fesetround (FE_TOWARDZERO))
++ {
++ TEST_ff_f (pow, 1.0625L, 1.125L, 1.070582293028761362162622578677070098674L);
++ TEST_ff_f (pow, 1.5L, 1.03125L, 1.519127098714743184071644334163037684948L);
++ }
++
++ fesetround (save_round_mode);
++
++ END (pow_towardzero);
++}
++
++
++static void
++pow_test_downward (void)
++{
++ int save_round_mode;
++ errno = 0;
++ FUNC(pow) (0, 0);
++ if (errno == ENOSYS)
++ /* Function not implemented. */
++ return;
++
++ START (pow_downward);
++
++ save_round_mode = fegetround ();
++
++ if (!fesetround (FE_DOWNWARD))
++ {
++ TEST_ff_f (pow, 1.0625L, 1.125L, 1.070582293028761362162622578677070098674L);
++ TEST_ff_f (pow, 1.5L, 1.03125L, 1.519127098714743184071644334163037684948L);
++ }
++
++ fesetround (save_round_mode);
++
++ END (pow_downward);
++}
++
++
++static void
++pow_test_upward (void)
++{
++ int save_round_mode;
++ errno = 0;
++ FUNC(pow) (0, 0);
++ if (errno == ENOSYS)
++ /* Function not implemented. */
++ return;
++
++ START (pow_upward);
++
++ save_round_mode = fegetround ();
++
++ if (!fesetround (FE_UPWARD))
++ {
++ TEST_ff_f (pow, 1.0625L, 1.125L, 1.070582293028761362162622578677070098674L);
++ TEST_ff_f (pow, 1.5L, 1.03125L, 1.519127098714743184071644334163037684948L);
++ }
++
++ fesetround (save_round_mode);
++
++ END (pow_upward);
++}
++
++
+ static void
+ remainder_test (void)
+ {
+@@ -6824,6 +6929,10 @@ main (int argc, char **argv)
+ fabs_test ();
+ hypot_test ();
+ pow_test ();
++ pow_test_tonearest ();
++ pow_test_towardzero ();
++ pow_test_downward ();
++ pow_test_upward ();
+ sqrt_test ();
+
+ /* Error and gamma functions: */
+--- a/sysdeps/i386/fpu/libm-test-ulps
++++ b/sysdeps/i386/fpu/libm-test-ulps
+@@ -844,6 +844,42 @@ ifloat: 1
+ ildouble: 1
+ ldouble: 1
+
++# pow_downward
++Test "pow_downward (1.0625, 1.125) == 1.070582293028761362162622578677070098674":
++double: 1
++idouble: 1
++ildouble: 1
++ldouble: 1
++Test "pow_downward (1.5, 1.03125) == 1.519127098714743184071644334163037684948":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++# pow_towardzero
++Test "pow_towardzero (1.0625, 1.125) == 1.070582293028761362162622578677070098674":
++double: 1
++idouble: 1
++ildouble: 1
++ldouble: 1
++Test "pow_towardzero (1.5, 1.03125) == 1.519127098714743184071644334163037684948":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++# pow_upward
++Test "pow_upward (1.0625, 1.125) == 1.070582293028761362162622578677070098674":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "pow_upward (1.5, 1.03125) == 1.519127098714743184071644334163037684948":
++double: 1
++idouble: 1
++ildouble: 1
++ldouble: 1
++
+ # sin_downward
+ Test "sin_downward (1) == 0.8414709848078965066525023216302989996226":
+ ildouble: 1
+@@ -1685,6 +1721,30 @@ float: 1
+ ifloat: 1
+ ildouble: 1
+ ldouble: 1
++
++Function: "pow_downward":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++Function: "pow_towardzero":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++Function: "pow_upward":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
+
+ Function: "sin_downward":
+ double: 1
+--- a/sysdeps/ieee754/dbl-64/e_pow.c
++++ b/sysdeps/ieee754/dbl-64/e_pow.c
+@@ -1,7 +1,7 @@
+ /*
+ * IBM Accurate Mathematical Library
+ * written by International Business Machines Corp.
+- * Copyright (C) 2001, 2002, 2004 Free Software Foundation
++ * Copyright (C) 2001, 2002, 2004, 2012 Free Software Foundation
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+@@ -42,6 +42,7 @@
+ #include "MathLib.h"
+ #include "upow.tbl"
+ #include "math_private.h"
++#include <fenv.h>
+
+
+ double __exp1(double x, double xx, double error);
+@@ -79,6 +80,12 @@ double __ieee754_pow(double x, double y)
+ (u.i[HIGH_HALF]==0 && u.i[LOW_HALF]!=0)) &&
+ /* 2^-1023< x<= 2^-1023 * 0x1.0000ffffffff */
+ (v.i[HIGH_HALF]&0x7fffffff) < 0x4ff00000) { /* if y<-1 or y>1 */
++ fenv_t env;
++ double retval;
++
++ feholdexcept (&env);
++ fesetround (FE_TONEAREST);
++
+ z = log1(x,&aa,&error); /* x^y =e^(y log (X)) */
+ t = y*134217729.0;
+ y1 = t - (t-y);
+@@ -92,7 +99,10 @@ double __ieee754_pow(double x, double y)
+ a2 = (a-a1)+aa;
+ error = error*ABS(y);
+ t = __exp1(a1,a2,1.9e16*error); /* return -10 or 0 if wasn't computed exactly */
+- return (t>0)?t:power1(x,y);
++ retval = (t>0)?t:power1(x,y);
++
++ feupdateenv (&env);
++ return retval;
+ }
+
+ if (x == 0) {
+--- a/sysdeps/x86_64/fpu/libm-test-ulps
++++ b/sysdeps/x86_64/fpu/libm-test-ulps
+@@ -867,6 +867,36 @@ Test "log1p (-0.25) == -0.28768207245178
+ float: 1
+ ifloat: 1
+
++# pow_downward
++Test "pow_downward (1.0625, 1.125) == 1.070582293028761362162622578677070098674":
++ildouble: 1
++ldouble: 1
++Test "pow_downward (1.5, 1.03125) == 1.519127098714743184071644334163037684948":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++# pow_towardzero
++Test "pow_towardzero (1.0625, 1.125) == 1.070582293028761362162622578677070098674":
++ildouble: 1
++ldouble: 1
++Test "pow_towardzero (1.5, 1.03125) == 1.519127098714743184071644334163037684948":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++# pow_upward
++Test "pow_upward (1.0625, 1.125) == 1.070582293028761362162622578677070098674":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "pow_upward (1.5, 1.03125) == 1.519127098714743184071644334163037684948":
++ildouble: 1
++ldouble: 1
++
+ # sin_downward
+ Test "sin_downward (1) == 0.8414709848078965066525023216302989996226":
+ ildouble: 1
+@@ -1644,6 +1674,24 @@ Function: "log1p":
+ float: 1
+ ifloat: 1
+
++Function: "pow_downward":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++Function: "pow_towardzero":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++Function: "pow_upward":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
+ Function: "sin_downward":
+ float: 1
+ ifloat: 1
only in patch2:
unchanged:
--- eglibc-2.13.orig/debian/patches/any/cvs-sin-cos-tan-rounding.diff
+++ eglibc-2.13/debian/patches/any/cvs-sin-cos-tan-rounding.diff
@@ -0,0 +1,2099 @@
+2012-03-02 Joseph Myers <joseph@codesourcery.com>
+
+ [BZ #3976]
+ * sysdeps/ieee754/dbl-64/s_sin.c: Include <fenv.h>
+ (__sin): Save and restore rounding mode and use round-to-nearest
+ for all computations.
+ (__cos): Save and restore rounding mode and use round-to-nearest
+ for all computations.
+ * sysdeps/ieee754/dbl-64/s_tan.c: Include "math_private.h" and
+ <fenv.h>.
+ (tan): Save and restore rounding mode and use round-to-nearest for
+ all computations.
+ * math/libm-test.inc (cos_test_tonearest): New function.
+ (cos_test_towardzero): Likewise.
+ (cos_test_downward): Likewise.
+ (cos_test_upward): Likewise.
+ (sin_test_tonearest): Likewise.
+ (sin_test_towardzero): Likewise.
+ (sin_test_downward): Likewise.
+ (sin_test_upward): Likewise.
+ (tan_test_tonearest): Likewise.
+ (tan_test_towardzero): Likewise.
+ (tan_test_downward): Likewise.
+ (tan_test_upward): Likewise.
+ (main): Call the new functions.
+ * sysdeps/i386/fpu/libm-test-ulps: Update.
+ * sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
+---
+ math/libm-test.inc | 421 +++++++++++++++++++++++++++++++++
+ sysdeps/i386/fpu/libm-test-ulps | 483 ++++++++++++++++++++++++++++++++++++++
+ sysdeps/ieee754/dbl-64/s_sin.c | 111 ++++++---
+ sysdeps/ieee754/dbl-64/s_tan.c | 146 +++++++-----
+ sysdeps/x86_64/fpu/libm-test-ulps | 378 +++++++++++++++++++++++++++++
+ 6 files changed, 1470 insertions(+), 97 deletions(-)
+
+--- a/math/libm-test.inc
++++ b/math/libm-test.inc
+@@ -2026,6 +2026,142 @@ cos_test (void)
+
+
+ static void
++cos_test_tonearest (void)
++{
++ int save_round_mode;
++ errno = 0;
++ FUNC(cos) (0);
++ if (errno == ENOSYS)
++ /* Function not implemented. */
++ return;
++
++ START (cos_tonearest);
++
++ save_round_mode = fegetround ();
++
++ if (!fesetround (FE_TONEAREST))
++ {
++ TEST_f_f (cos, 1, 0.5403023058681397174009366074429766037323L);
++ TEST_f_f (cos, 2, -0.4161468365471423869975682295007621897660L);
++ TEST_f_f (cos, 3, -0.9899924966004454572715727947312613023937L);
++ TEST_f_f (cos, 4, -0.6536436208636119146391681830977503814241L);
++ TEST_f_f (cos, 5, 0.2836621854632262644666391715135573083344L);
++ TEST_f_f (cos, 6, 0.9601702866503660205456522979229244054519L);
++ TEST_f_f (cos, 7, 0.7539022543433046381411975217191820122183L);
++ TEST_f_f (cos, 8, -0.1455000338086135258688413818311946826093L);
++ TEST_f_f (cos, 9, -0.9111302618846769883682947111811653112463L);
++ TEST_f_f (cos, 10, -0.8390715290764524522588639478240648345199L);
++ }
++
++ fesetround (save_round_mode);
++
++ END (cos_tonearest);
++}
++
++
++static void
++cos_test_towardzero (void)
++{
++ int save_round_mode;
++ errno = 0;
++ FUNC(cos) (0);
++ if (errno == ENOSYS)
++ /* Function not implemented. */
++ return;
++
++ START (cos_towardzero);
++
++ save_round_mode = fegetround ();
++
++ if (!fesetround (FE_TOWARDZERO))
++ {
++ TEST_f_f (cos, 1, 0.5403023058681397174009366074429766037323L);
++ TEST_f_f (cos, 2, -0.4161468365471423869975682295007621897660L);
++ TEST_f_f (cos, 3, -0.9899924966004454572715727947312613023937L);
++ TEST_f_f (cos, 4, -0.6536436208636119146391681830977503814241L);
++ TEST_f_f (cos, 5, 0.2836621854632262644666391715135573083344L);
++ TEST_f_f (cos, 6, 0.9601702866503660205456522979229244054519L);
++ TEST_f_f (cos, 7, 0.7539022543433046381411975217191820122183L);
++ TEST_f_f (cos, 8, -0.1455000338086135258688413818311946826093L);
++ TEST_f_f (cos, 9, -0.9111302618846769883682947111811653112463L);
++ TEST_f_f (cos, 10, -0.8390715290764524522588639478240648345199L);
++ }
++
++ fesetround (save_round_mode);
++
++ END (cos_towardzero);
++}
++
++
++static void
++cos_test_downward (void)
++{
++ int save_round_mode;
++ errno = 0;
++ FUNC(cos) (0);
++ if (errno == ENOSYS)
++ /* Function not implemented. */
++ return;
++
++ START (cos_downward);
++
++ save_round_mode = fegetround ();
++
++ if (!fesetround (FE_DOWNWARD))
++ {
++ TEST_f_f (cos, 1, 0.5403023058681397174009366074429766037323L);
++ TEST_f_f (cos, 2, -0.4161468365471423869975682295007621897660L);
++ TEST_f_f (cos, 3, -0.9899924966004454572715727947312613023937L);
++ TEST_f_f (cos, 4, -0.6536436208636119146391681830977503814241L);
++ TEST_f_f (cos, 5, 0.2836621854632262644666391715135573083344L);
++ TEST_f_f (cos, 6, 0.9601702866503660205456522979229244054519L);
++ TEST_f_f (cos, 7, 0.7539022543433046381411975217191820122183L);
++ TEST_f_f (cos, 8, -0.1455000338086135258688413818311946826093L);
++ TEST_f_f (cos, 9, -0.9111302618846769883682947111811653112463L);
++ TEST_f_f (cos, 10, -0.8390715290764524522588639478240648345199L);
++ }
++
++ fesetround (save_round_mode);
++
++ END (cos_downward);
++}
++
++
++static void
++cos_test_upward (void)
++{
++ int save_round_mode;
++ errno = 0;
++ FUNC(cos) (0);
++ if (errno == ENOSYS)
++ /* Function not implemented. */
++ return;
++
++ START (cos_upward);
++
++ save_round_mode = fegetround ();
++
++ if (!fesetround (FE_UPWARD))
++ {
++ TEST_f_f (cos, 1, 0.5403023058681397174009366074429766037323L);
++ TEST_f_f (cos, 2, -0.4161468365471423869975682295007621897660L);
++ TEST_f_f (cos, 3, -0.9899924966004454572715727947312613023937L);
++ TEST_f_f (cos, 4, -0.6536436208636119146391681830977503814241L);
++ TEST_f_f (cos, 5, 0.2836621854632262644666391715135573083344L);
++ TEST_f_f (cos, 6, 0.9601702866503660205456522979229244054519L);
++ TEST_f_f (cos, 7, 0.7539022543433046381411975217191820122183L);
++ TEST_f_f (cos, 8, -0.1455000338086135258688413818311946826093L);
++ TEST_f_f (cos, 9, -0.9111302618846769883682947111811653112463L);
++ TEST_f_f (cos, 10, -0.8390715290764524522588639478240648345199L);
++ }
++
++ fesetround (save_round_mode);
++
++ END (cos_upward);
++}
++
++
++static void
+ cosh_test (void)
+ {
+ errno = 0;
+@@ -5745,6 +5881,142 @@ sin_test (void)
+
+
+ static void
++sin_test_tonearest (void)
++{
++ int save_round_mode;
++ errno = 0;
++ FUNC(sin) (0);
++ if (errno == ENOSYS)
++ /* Function not implemented. */
++ return;
++
++ START (sin_tonearest);
++
++ save_round_mode = fegetround ();
++
++ if (!fesetround (FE_TONEAREST))
++ {
++ TEST_f_f (sin, 1, 0.8414709848078965066525023216302989996226L);
++ TEST_f_f (sin, 2, 0.9092974268256816953960198659117448427023L);
++ TEST_f_f (sin, 3, 0.1411200080598672221007448028081102798469L);
++ TEST_f_f (sin, 4, -0.7568024953079282513726390945118290941359L);
++ TEST_f_f (sin, 5, -0.9589242746631384688931544061559939733525L);
++ TEST_f_f (sin, 6, -0.2794154981989258728115554466118947596280L);
++ TEST_f_f (sin, 7, 0.6569865987187890903969990915936351779369L);
++ TEST_f_f (sin, 8, 0.9893582466233817778081235982452886721164L);
++ TEST_f_f (sin, 9, 0.4121184852417565697562725663524351793439L);
++ TEST_f_f (sin, 10, -0.5440211108893698134047476618513772816836L);
++ }
++
++ fesetround (save_round_mode);
++
++ END (sin_tonearest);
++}
++
++
++static void
++sin_test_towardzero (void)
++{
++ int save_round_mode;
++ errno = 0;
++ FUNC(sin) (0);
++ if (errno == ENOSYS)
++ /* Function not implemented. */
++ return;
++
++ START (sin_towardzero);
++
++ save_round_mode = fegetround ();
++
++ if (!fesetround (FE_TOWARDZERO))
++ {
++ TEST_f_f (sin, 1, 0.8414709848078965066525023216302989996226L);
++ TEST_f_f (sin, 2, 0.9092974268256816953960198659117448427023L);
++ TEST_f_f (sin, 3, 0.1411200080598672221007448028081102798469L);
++ TEST_f_f (sin, 4, -0.7568024953079282513726390945118290941359L);
++ TEST_f_f (sin, 5, -0.9589242746631384688931544061559939733525L);
++ TEST_f_f (sin, 6, -0.2794154981989258728115554466118947596280L);
++ TEST_f_f (sin, 7, 0.6569865987187890903969990915936351779369L);
++ TEST_f_f (sin, 8, 0.9893582466233817778081235982452886721164L);
++ TEST_f_f (sin, 9, 0.4121184852417565697562725663524351793439L);
++ TEST_f_f (sin, 10, -0.5440211108893698134047476618513772816836L);
++ }
++
++ fesetround (save_round_mode);
++
++ END (sin_towardzero);
++}
++
++
++static void
++sin_test_downward (void)
++{
++ int save_round_mode;
++ errno = 0;
++ FUNC(sin) (0);
++ if (errno == ENOSYS)
++ /* Function not implemented. */
++ return;
++
++ START (sin_downward);
++
++ save_round_mode = fegetround ();
++
++ if (!fesetround (FE_DOWNWARD))
++ {
++ TEST_f_f (sin, 1, 0.8414709848078965066525023216302989996226L);
++ TEST_f_f (sin, 2, 0.9092974268256816953960198659117448427023L);
++ TEST_f_f (sin, 3, 0.1411200080598672221007448028081102798469L);
++ TEST_f_f (sin, 4, -0.7568024953079282513726390945118290941359L);
++ TEST_f_f (sin, 5, -0.9589242746631384688931544061559939733525L);
++ TEST_f_f (sin, 6, -0.2794154981989258728115554466118947596280L);
++ TEST_f_f (sin, 7, 0.6569865987187890903969990915936351779369L);
++ TEST_f_f (sin, 8, 0.9893582466233817778081235982452886721164L);
++ TEST_f_f (sin, 9, 0.4121184852417565697562725663524351793439L);
++ TEST_f_f (sin, 10, -0.5440211108893698134047476618513772816836L);
++ }
++
++ fesetround (save_round_mode);
++
++ END (sin_downward);
++}
++
++
++static void
++sin_test_upward (void)
++{
++ int save_round_mode;
++ errno = 0;
++ FUNC(sin) (0);
++ if (errno == ENOSYS)
++ /* Function not implemented. */
++ return;
++
++ START (sin_upward);
++
++ save_round_mode = fegetround ();
++
++ if (!fesetround (FE_UPWARD))
++ {
++ TEST_f_f (sin, 1, 0.8414709848078965066525023216302989996226L);
++ TEST_f_f (sin, 2, 0.9092974268256816953960198659117448427023L);
++ TEST_f_f (sin, 3, 0.1411200080598672221007448028081102798469L);
++ TEST_f_f (sin, 4, -0.7568024953079282513726390945118290941359L);
++ TEST_f_f (sin, 5, -0.9589242746631384688931544061559939733525L);
++ TEST_f_f (sin, 6, -0.2794154981989258728115554466118947596280L);
++ TEST_f_f (sin, 7, 0.6569865987187890903969990915936351779369L);
++ TEST_f_f (sin, 8, 0.9893582466233817778081235982452886721164L);
++ TEST_f_f (sin, 9, 0.4121184852417565697562725663524351793439L);
++ TEST_f_f (sin, 10, -0.5440211108893698134047476618513772816836L);
++ }
++
++ fesetround (save_round_mode);
++
++ END (sin_upward);
++}
++
++
++static void
+ sincos_test (void)
+ {
+ FLOAT sin_res, cos_res;
+@@ -5865,6 +6137,143 @@ tan_test (void)
+ END (tan);
+ }
+
++
++static void
++tan_test_tonearest (void)
++{
++ int save_round_mode;
++ errno = 0;
++ FUNC(tan) (0);
++ if (errno == ENOSYS)
++ /* Function not implemented. */
++ return;
++
++ START (tan_tonearest);
++
++ save_round_mode = fegetround ();
++
++ if (!fesetround (FE_TONEAREST))
++ {
++ TEST_f_f (tan, 1, 1.5574077246549022305069748074583601730873L);
++ TEST_f_f (tan, 2, -2.1850398632615189916433061023136825434320L);
++ TEST_f_f (tan, 3, -0.1425465430742778052956354105339134932261L);
++ TEST_f_f (tan, 4, 1.1578212823495775831373424182673239231198L);
++ TEST_f_f (tan, 5, -3.3805150062465856369827058794473439087096L);
++ TEST_f_f (tan, 6, -0.2910061913847491570536995888681755428312L);
++ TEST_f_f (tan, 7, 0.8714479827243187364564508896003135663222L);
++ TEST_f_f (tan, 8, -6.7997114552203786999252627596086333648814L);
++ TEST_f_f (tan, 9, -0.4523156594418098405903708757987855343087L);
++ TEST_f_f (tan, 10, 0.6483608274590866712591249330098086768169L);
++ }
++
++ fesetround (save_round_mode);
++
++ END (tan_tonearest);
++}
++
++
++static void
++tan_test_towardzero (void)
++{
++ int save_round_mode;
++ errno = 0;
++ FUNC(tan) (0);
++ if (errno == ENOSYS)
++ /* Function not implemented. */
++ return;
++
++ START (tan_towardzero);
++
++ save_round_mode = fegetround ();
++
++ if (!fesetround (FE_TOWARDZERO))
++ {
++ TEST_f_f (tan, 1, 1.5574077246549022305069748074583601730873L);
++ TEST_f_f (tan, 2, -2.1850398632615189916433061023136825434320L);
++ TEST_f_f (tan, 3, -0.1425465430742778052956354105339134932261L);
++ TEST_f_f (tan, 4, 1.1578212823495775831373424182673239231198L);
++ TEST_f_f (tan, 5, -3.3805150062465856369827058794473439087096L);
++ TEST_f_f (tan, 6, -0.2910061913847491570536995888681755428312L);
++ TEST_f_f (tan, 7, 0.8714479827243187364564508896003135663222L);
++ TEST_f_f (tan, 8, -6.7997114552203786999252627596086333648814L);
++ TEST_f_f (tan, 9, -0.4523156594418098405903708757987855343087L);
++ TEST_f_f (tan, 10, 0.6483608274590866712591249330098086768169L);
++ }
++
++ fesetround (save_round_mode);
++
++ END (tan_towardzero);
++}
++
++
++static void
++tan_test_downward (void)
++{
++ int save_round_mode;
++ errno = 0;
++ FUNC(tan) (0);
++ if (errno == ENOSYS)
++ /* Function not implemented. */
++ return;
++
++ START (tan_downward);
++
++ save_round_mode = fegetround ();
++
++ if (!fesetround (FE_DOWNWARD))
++ {
++ TEST_f_f (tan, 1, 1.5574077246549022305069748074583601730873L);
++ TEST_f_f (tan, 2, -2.1850398632615189916433061023136825434320L);
++ TEST_f_f (tan, 3, -0.1425465430742778052956354105339134932261L);
++ TEST_f_f (tan, 4, 1.1578212823495775831373424182673239231198L);
++ TEST_f_f (tan, 5, -3.3805150062465856369827058794473439087096L);
++ TEST_f_f (tan, 6, -0.2910061913847491570536995888681755428312L);
++ TEST_f_f (tan, 7, 0.8714479827243187364564508896003135663222L);
++ TEST_f_f (tan, 8, -6.7997114552203786999252627596086333648814L);
++ TEST_f_f (tan, 9, -0.4523156594418098405903708757987855343087L);
++ TEST_f_f (tan, 10, 0.6483608274590866712591249330098086768169L);
++ }
++
++ fesetround (save_round_mode);
++
++ END (tan_downward);
++}
++
++
++static void
++tan_test_upward (void)
++{
++ int save_round_mode;
++ errno = 0;
++ FUNC(tan) (0);
++ if (errno == ENOSYS)
++ /* Function not implemented. */
++ return;
++
++ START (tan_upward);
++
++ save_round_mode = fegetround ();
++
++ if (!fesetround (FE_UPWARD))
++ {
++ TEST_f_f (tan, 1, 1.5574077246549022305069748074583601730873L);
++ TEST_f_f (tan, 2, -2.1850398632615189916433061023136825434320L);
++ TEST_f_f (tan, 3, -0.1425465430742778052956354105339134932261L);
++ TEST_f_f (tan, 4, 1.1578212823495775831373424182673239231198L);
++ TEST_f_f (tan, 5, -3.3805150062465856369827058794473439087096L);
++ TEST_f_f (tan, 6, -0.2910061913847491570536995888681755428312L);
++ TEST_f_f (tan, 7, 0.8714479827243187364564508896003135663222L);
++ TEST_f_f (tan, 8, -6.7997114552203786999252627596086333648814L);
++ TEST_f_f (tan, 9, -0.4523156594418098405903708757987855343087L);
++ TEST_f_f (tan, 10, 0.6483608274590866712591249330098086768169L);
++ }
++
++ fesetround (save_round_mode);
++
++ END (tan_upward);
++}
++
++
+ static void
+ tanh_test (void)
+ {
+@@ -6363,9 +6772,21 @@ main (int argc, char **argv)
+ atan_test ();
+ atan2_test ();
+ cos_test ();
++ cos_test_tonearest ();
++ cos_test_towardzero ();
++ cos_test_downward ();
++ cos_test_upward ();
+ sin_test ();
++ sin_test_tonearest ();
++ sin_test_towardzero ();
++ sin_test_downward ();
++ sin_test_upward ();
+ sincos_test ();
+ tan_test ();
++ tan_test_tonearest ();
++ tan_test_towardzero ();
++ tan_test_downward ();
++ tan_test_upward ();
+
+ /* Hyperbolic functions: */
+ acosh_test ();
+--- a/sysdeps/i386/fpu/libm-test-ulps
++++ b/sysdeps/i386/fpu/libm-test-ulps
+@@ -280,6 +280,130 @@ ifloat: 1
+ ildouble: 1
+ ldouble: 1
+
++# cos_downward
++Test "cos_downward (1) == 0.5403023058681397174009366074429766037323":
++double: 1
++idouble: 1
++ildouble: 1
++ldouble: 1
++Test "cos_downward (10) == -0.8390715290764524522588639478240648345199":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++Test "cos_downward (3) == -0.9899924966004454572715727947312613023937":
++double: 1
++idouble: 1
++Test "cos_downward (4) == -0.6536436208636119146391681830977503814241":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "cos_downward (5) == 0.2836621854632262644666391715135573083344":
++float: 1
++ifloat: 1
++Test "cos_downward (7) == 0.7539022543433046381411975217191820122183":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "cos_downward (8) == -0.1455000338086135258688413818311946826093":
++ildouble: 1
++ldouble: 1
++Test "cos_downward (9) == -0.9111302618846769883682947111811653112463":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++# cos_tonearest
++Test "cos_tonearest (8) == -0.1455000338086135258688413818311946826093":
++ildouble: 1
++ldouble: 1
++
++# cos_towardzero
++Test "cos_towardzero (1) == 0.5403023058681397174009366074429766037323":
++double: 1
++idouble: 1
++ildouble: 1
++ldouble: 1
++Test "cos_towardzero (10) == -0.8390715290764524522588639478240648345199":
++ildouble: 1
++ldouble: 1
++Test "cos_towardzero (2) == -0.4161468365471423869975682295007621897660":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "cos_towardzero (3) == -0.9899924966004454572715727947312613023937":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "cos_towardzero (4) == -0.6536436208636119146391681830977503814241":
++double: 1
++idouble: 1
++Test "cos_towardzero (5) == 0.2836621854632262644666391715135573083344":
++float: 1
++ifloat: 1
++Test "cos_towardzero (7) == 0.7539022543433046381411975217191820122183":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "cos_towardzero (8) == -0.1455000338086135258688413818311946826093":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++
++# cos_upward
++Test "cos_upward (1) == 0.5403023058681397174009366074429766037323":
++float: 1
++ifloat: 1
++Test "cos_upward (10) == -0.8390715290764524522588639478240648345199":
++ildouble: 1
++ldouble: 1
++Test "cos_upward (2) == -0.4161468365471423869975682295007621897660":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "cos_upward (3) == -0.9899924966004454572715727947312613023937":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "cos_upward (4) == -0.6536436208636119146391681830977503814241":
++double: 1
++idouble: 1
++Test "cos_upward (5) == 0.2836621854632262644666391715135573083344":
++double: 1
++idouble: 1
++ildouble: 1
++ldouble: 1
++Test "cos_upward (6) == 0.9601702866503660205456522979229244054519":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "cos_upward (7) == 0.7539022543433046381411975217191820122183":
++double: 1
++idouble: 1
++Test "cos_upward (8) == -0.1455000338086135258688413818311946826093":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++
+ # cosh
+ Test "cosh (0.75) == 1.29468328467684468784170818539018176":
+ ildouble: 1
+@@ -720,6 +844,135 @@ ifloat: 1
+ ildouble: 1
+ ldouble: 1
+
++# sin_downward
++Test "sin_downward (1) == 0.8414709848078965066525023216302989996226":
++ildouble: 1
++ldouble: 1
++Test "sin_downward (10) == -0.5440211108893698134047476618513772816836":
++double: 1
++idouble: 1
++ildouble: 1
++ldouble: 1
++Test "sin_downward (2) == 0.9092974268256816953960198659117448427023":
++double: 1
++idouble: 1
++Test "sin_downward (3) == 0.1411200080598672221007448028081102798469":
++ildouble: 1
++ldouble: 1
++Test "sin_downward (4) == -0.7568024953079282513726390945118290941359":
++double: 1
++idouble: 1
++ildouble: 1
++ldouble: 1
++Test "sin_downward (5) == -0.9589242746631384688931544061559939733525":
++double: 1
++idouble: 1
++ildouble: 1
++ldouble: 1
++Test "sin_downward (6) == -0.2794154981989258728115554466118947596280":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++Test "sin_downward (7) == 0.6569865987187890903969990915936351779369":
++ildouble: 1
++ldouble: 1
++Test "sin_downward (8) == 0.9893582466233817778081235982452886721164":
++double: 1
++idouble: 1
++ildouble: 1
++ldouble: 1
++Test "sin_downward (9) == 0.4121184852417565697562725663524351793439":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++# sin_tonearest
++Test "sin_tonearest (10) == -0.5440211108893698134047476618513772816836":
++ildouble: 1
++ldouble: 1
++Test "sin_tonearest (4) == -0.7568024953079282513726390945118290941359":
++ildouble: 1
++ldouble: 1
++
++# sin_towardzero
++Test "sin_towardzero (1) == 0.8414709848078965066525023216302989996226":
++ildouble: 1
++ldouble: 1
++Test "sin_towardzero (10) == -0.5440211108893698134047476618513772816836":
++float: 1
++ifloat: 1
++Test "sin_towardzero (2) == 0.9092974268256816953960198659117448427023":
++double: 1
++idouble: 1
++Test "sin_towardzero (3) == 0.1411200080598672221007448028081102798469":
++ildouble: 1
++ldouble: 1
++Test "sin_towardzero (4) == -0.7568024953079282513726390945118290941359":
++float: 1
++ifloat: 1
++Test "sin_towardzero (5) == -0.9589242746631384688931544061559939733525":
++float: 1
++ifloat: 1
++Test "sin_towardzero (6) == -0.2794154981989258728115554466118947596280":
++ildouble: 1
++ldouble: 1
++Test "sin_towardzero (7) == 0.6569865987187890903969990915936351779369":
++ildouble: 1
++ldouble: 1
++Test "sin_towardzero (8) == 0.9893582466233817778081235982452886721164":
++double: 1
++idouble: 1
++ildouble: 1
++ldouble: 1
++Test "sin_towardzero (9) == 0.4121184852417565697562725663524351793439":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++# sin_upward
++Test "sin_upward (1) == 0.8414709848078965066525023216302989996226":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++Test "sin_upward (10) == -0.5440211108893698134047476618513772816836":
++float: 1
++ifloat: 1
++Test "sin_upward (2) == 0.9092974268256816953960198659117448427023":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "sin_upward (3) == 0.1411200080598672221007448028081102798469":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++Test "sin_upward (4) == -0.7568024953079282513726390945118290941359":
++float: 1
++ifloat: 1
++Test "sin_upward (5) == -0.9589242746631384688931544061559939733525":
++float: 1
++ifloat: 1
++Test "sin_upward (6) == -0.2794154981989258728115554466118947596280":
++ildouble: 1
++ldouble: 1
++Test "sin_upward (7) == 0.6569865987187890903969990915936351779369":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++Test "sin_upward (8) == 0.9893582466233817778081235982452886721164":
++float: 1
++ifloat: 1
++
+ # sincos
+ Test "sincos (M_PI_6l*2.0, &sin_res, &cos_res) puts 0.5 in cos_res":
+ double: 1
+@@ -751,6 +1004,152 @@ Test "tan (pi/4) == 1":
+ double: 1
+ idouble: 1
+
++# tan_downward
++Test "tan_downward (1) == 1.5574077246549022305069748074583601730873":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++Test "tan_downward (10) == 0.6483608274590866712591249330098086768169":
++float: 1
++ifloat: 1
++Test "tan_downward (2) == -2.1850398632615189916433061023136825434320":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++Test "tan_downward (3) == -0.1425465430742778052956354105339134932261":
++double: 1
++idouble: 1
++ildouble: 1
++ldouble: 1
++Test "tan_downward (4) == 1.1578212823495775831373424182673239231198":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "tan_downward (6) == -0.2910061913847491570536995888681755428312":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++Test "tan_downward (7) == 0.8714479827243187364564508896003135663222":
++double: 1
++idouble: 1
++Test "tan_downward (8) == -6.7997114552203786999252627596086333648814":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "tan_downward (9) == -0.4523156594418098405903708757987855343087":
++float: 1
++ifloat: 1
++
++# tan_tonearest
++Test "tan_tonearest (1) == 1.5574077246549022305069748074583601730873":
++ildouble: 1
++ldouble: 1
++Test "tan_tonearest (2) == -2.1850398632615189916433061023136825434320":
++ildouble: 1
++ldouble: 1
++Test "tan_tonearest (6) == -0.2910061913847491570536995888681755428312":
++ildouble: 1
++ldouble: 1
++Test "tan_tonearest (8) == -6.7997114552203786999252627596086333648814":
++ildouble: 1
++ldouble: 1
++Test "tan_tonearest (9) == -0.4523156594418098405903708757987855343087":
++ildouble: 1
++ldouble: 1
++
++# tan_towardzero
++Test "tan_towardzero (1) == 1.5574077246549022305069748074583601730873":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++Test "tan_towardzero (10) == 0.6483608274590866712591249330098086768169":
++float: 1
++ifloat: 1
++Test "tan_towardzero (2) == -2.1850398632615189916433061023136825434320":
++ildouble: 1
++ldouble: 1
++Test "tan_towardzero (3) == -0.1425465430742778052956354105339134932261":
++float: 1
++ifloat: 1
++Test "tan_towardzero (4) == 1.1578212823495775831373424182673239231198":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "tan_towardzero (5) == -3.3805150062465856369827058794473439087096":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "tan_towardzero (6) == -0.2910061913847491570536995888681755428312":
++ildouble: 1
++ldouble: 1
++Test "tan_towardzero (7) == 0.8714479827243187364564508896003135663222":
++double: 1
++idouble: 1
++Test "tan_towardzero (8) == -6.7997114552203786999252627596086333648814":
++double: 1
++idouble: 1
++ildouble: 2
++ldouble: 2
++Test "tan_towardzero (9) == -0.4523156594418098405903708757987855343087":
++double: 1
++idouble: 1
++ildouble: 1
++ldouble: 1
++
++# tan_upward
++Test "tan_upward (1) == 1.5574077246549022305069748074583601730873":
++ildouble: 1
++ldouble: 1
++Test "tan_upward (10) == 0.6483608274590866712591249330098086768169":
++double: 1
++idouble: 1
++ildouble: 1
++ldouble: 1
++Test "tan_upward (2) == -2.1850398632615189916433061023136825434320":
++ildouble: 1
++ldouble: 1
++Test "tan_upward (3) == -0.1425465430742778052956354105339134932261":
++float: 1
++ifloat: 1
++Test "tan_upward (4) == 1.1578212823495775831373424182673239231198":
++double: 1
++idouble: 1
++Test "tan_upward (5) == -3.3805150062465856369827058794473439087096":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "tan_upward (6) == -0.2910061913847491570536995888681755428312":
++ildouble: 1
++ldouble: 1
++Test "tan_upward (7) == 0.8714479827243187364564508896003135663222":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "tan_upward (8) == -6.7997114552203786999252627596086333648814":
++double: 1
++idouble: 1
++ildouble: 2
++ldouble: 2
++Test "tan_upward (9) == -0.4523156594418098405903708757987855343087":
++double: 1
++idouble: 1
++ildouble: 1
++ldouble: 1
++
+ # tgamma
+ Test "tgamma (-0.5) == -2 sqrt (pi)":
+ double: 2
+@@ -1089,6 +1488,34 @@ ifloat: 1
+ ildouble: 1
+ ldouble: 1
+
++Function: "cos_downward":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++Function: "cos_tonearest":
++ildouble: 1
++ldouble: 1
++
++Function: "cos_towardzero":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++Function: "cos_upward":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
+ Function: "cosh":
+ ildouble: 1
+
+@@ -1259,6 +1686,34 @@ ifloat: 1
+ ildouble: 1
+ ldouble: 1
+
++Function: "sin_downward":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++Function: "sin_tonearest":
++ildouble: 1
++ldouble: 1
++
++Function: "sin_towardzero":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++Function: "sin_upward":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
+ Function: "sincos":
+ double: 1
+ float: 1
+@@ -1275,6 +1730,34 @@ Function: "tan":
+ double: 1
+ idouble: 1
+
++Function: "tan_downward":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++Function: "tan_tonearest":
++ildouble: 1
++ldouble: 1
++
++Function: "tan_towardzero":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 2
++ldouble: 2
++
++Function: "tan_upward":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 2
++ldouble: 2
++
+ Function: "tgamma":
+ double: 2
+ float: 1
+--- a/sysdeps/ieee754/dbl-64/s_sin.c
++++ b/sysdeps/ieee754/dbl-64/s_sin.c
+@@ -55,6 +55,7 @@
+ #include "MathLib.h"
+ #include "sincos.tbl"
+ #include "math_private.h"
++#include <fenv.h>
+
+ static const double
+ sn3 = -1.66666666666664880952546298448555E-01,
+@@ -97,12 +98,17 @@ double __sin(double x){
+ #if 0
+ int4 nn;
+ #endif
++ fenv_t env;
++ double retval = 0;
++
++ feholdexcept (&env);
++ fesetround (FE_TONEAREST);
+
+ u.x = x;
+ m = u.i[HIGH_HALF];
+ k = 0x7fffffff&m; /* no sign */
+ if (k < 0x3e500000) /* if x->0 =>sin(x)=x */
+- return x;
++ { retval = x; goto ret; }
+ /*---------------------------- 2^-26 < |x|< 0.25 ----------------------*/
+ else if (k < 0x3fd00000){
+ xx = x*x;
+@@ -110,7 +116,8 @@ double __sin(double x){
+ t = ((((s5.x*xx + s4.x)*xx + s3.x)*xx + s2.x)*xx + s1.x)*(xx*x);
+ res = x+t;
+ cor = (x-res)+t;
+- return (res == res + 1.07*cor)? res : slow(x);
++ retval = (res == res + 1.07*cor)? res : slow(x);
++ goto ret;
+ } /* else if (k < 0x3fd00000) */
+ /*---------------------------- 0.25<|x|< 0.855469---------------------- */
+ else if (k < 0x3feb6000) {
+@@ -127,7 +134,8 @@ double __sin(double x){
+ cor=(ssn+s*ccs-sn*c)+cs*s;
+ res=sn+cor;
+ cor=(sn-res)+cor;
+- return (res==res+1.096*cor)? res : slow1(x);
++ retval = (res==res+1.096*cor)? res : slow1(x);
++ goto ret;
+ } /* else if (k < 0x3feb6000) */
+
+ /*----------------------- 0.855469 <|x|<2.426265 ----------------------*/
+@@ -153,7 +161,8 @@ double __sin(double x){
+ cor=(ccs-s*ssn-cs*c)-sn*s;
+ res=cs+cor;
+ cor=(cs-res)+cor;
+- return (res==res+1.020*cor)? ((m>0)?res:-res) : slow2(x);
++ retval = (res==res+1.020*cor)? ((m>0)?res:-res) : slow2(x);
++ goto ret;
+ } /* else if (k < 0x400368fd) */
+
+ /*-------------------------- 2.426265<|x|< 105414350 ----------------------*/
+@@ -179,7 +188,8 @@ double __sin(double x){
+ res = a+t;
+ cor = (a-res)+t;
+ cor = (cor>0)? 1.02*cor+eps : 1.02*cor -eps;
+- return (res == res + cor)? res : sloww(a,da,x);
++ retval = (res == res + cor)? res : sloww(a,da,x);
++ goto ret;
+ }
+ else {
+ if (a>0)
+@@ -200,7 +210,8 @@ double __sin(double x){
+ res=sn+cor;
+ cor=(sn-res)+cor;
+ cor = (cor>0)? 1.035*cor+eps : 1.035*cor-eps;
+- return (res==res+cor)? ((m)?res:-res) : sloww1(a,da,x);
++ retval = (res==res+cor)? ((m)?res:-res) : sloww1(a,da,x);
++ goto ret;
+ }
+ break;
+
+@@ -222,7 +233,8 @@ double __sin(double x){
+ res=cs+cor;
+ cor=(cs-res)+cor;
+ cor = (cor>0)? 1.025*cor+eps : 1.025*cor-eps;
+- return (res==res+cor)? ((n&2)?-res:res) : sloww2(a,da,x,n);
++ retval = (res==res+cor)? ((n&2)?-res:res) : sloww2(a,da,x,n);
++ goto ret;
+
+ break;
+
+@@ -258,7 +270,8 @@ double __sin(double x){
+ res = a+t;
+ cor = (a-res)+t;
+ cor = (cor>0)? 1.02*cor+eps : 1.02*cor -eps;
+- return (res == res + cor)? res : bsloww(a,da,x,n);
++ retval = (res == res + cor)? res : bsloww(a,da,x,n);
++ goto ret;
+ }
+ else {
+ if (a>0) {m=1;t=a;db=da;}
+@@ -277,7 +290,8 @@ double __sin(double x){
+ res=sn+cor;
+ cor=(sn-res)+cor;
+ cor = (cor>0)? 1.035*cor+eps : 1.035*cor-eps;
+- return (res==res+cor)? ((m)?res:-res) : bsloww1(a,da,x,n);
++ retval = (res==res+cor)? ((m)?res:-res) : bsloww1(a,da,x,n);
++ goto ret;
+ }
+ break;
+
+@@ -299,7 +313,8 @@ double __sin(double x){
+ res=cs+cor;
+ cor=(cs-res)+cor;
+ cor = (cor>0)? 1.025*cor+eps : 1.025*cor-eps;
+- return (res==res+cor)? ((n&2)?-res:res) : bsloww2(a,da,x,n);
++ retval = (res==res+cor)? ((n&2)?-res:res) : bsloww2(a,da,x,n);
++ goto ret;
+
+ break;
+
+@@ -313,17 +328,20 @@ double __sin(double x){
+ n = __branred(x,&a,&da);
+ switch (n) {
+ case 0:
+- if (a*a < 0.01588) return bsloww(a,da,x,n);
+- else return bsloww1(a,da,x,n);
++ if (a*a < 0.01588) retval = bsloww(a,da,x,n);
++ else retval = bsloww1(a,da,x,n);
++ goto ret;
+ break;
+ case 2:
+- if (a*a < 0.01588) return bsloww(-a,-da,x,n);
+- else return bsloww1(-a,-da,x,n);
++ if (a*a < 0.01588) retval = bsloww(-a,-da,x,n);
++ else retval = bsloww1(-a,-da,x,n);
++ goto ret;
+ break;
+
+ case 1:
+ case 3:
+- return bsloww2(a,da,x,n);
++ retval = bsloww2(a,da,x,n);
++ goto ret;
+ break;
+ }
+
+@@ -333,9 +351,13 @@ double __sin(double x){
+ else {
+ if (k == 0x7ff00000 && u.i[LOW_HALF] == 0)
+ __set_errno (EDOM);
+- return x / x;
++ retval = x / x;
++ goto ret;
+ }
+- return 0; /* unreachable */
++
++ ret:
++ feupdateenv (&env);
++ return retval;
+ }
+
+
+@@ -350,11 +372,17 @@ double __cos(double x)
+ mynumber u,v;
+ int4 k,m,n;
+
++ fenv_t env;
++ double retval = 0;
++
++ feholdexcept (&env);
++ fesetround (FE_TONEAREST);
++
+ u.x = x;
+ m = u.i[HIGH_HALF];
+ k = 0x7fffffff&m;
+
+- if (k < 0x3e400000 ) return 1.0; /* |x|<2^-27 => cos(x)=1 */
++ if (k < 0x3e400000 ) { retval = 1.0; goto ret; } /* |x|<2^-27 => cos(x)=1 */
+
+ else if (k < 0x3feb6000 ) {/* 2^-27 < |x| < 0.855469 */
+ y=ABS(x);
+@@ -371,7 +399,8 @@ double __cos(double x)
+ cor=(ccs-s*ssn-cs*c)-sn*s;
+ res=cs+cor;
+ cor=(cs-res)+cor;
+- return (res==res+1.020*cor)? res : cslow2(x);
++ retval = (res==res+1.020*cor)? res : cslow2(x);
++ goto ret;
+
+ } /* else if (k < 0x3feb6000) */
+
+@@ -385,7 +414,8 @@ double __cos(double x)
+ res = a+t;
+ cor = (a-res)+t;
+ cor = (cor>0)? 1.02*cor+1.0e-31 : 1.02*cor -1.0e-31;
+- return (res == res + cor)? res : csloww(a,da,x);
++ retval = (res == res + cor)? res : csloww(a,da,x);
++ goto ret;
+ }
+ else {
+ if (a>0) {m=1;t=a;db=da;}
+@@ -404,7 +434,8 @@ double __cos(double x)
+ res=sn+cor;
+ cor=(sn-res)+cor;
+ cor = (cor>0)? 1.035*cor+1.0e-31 : 1.035*cor-1.0e-31;
+- return (res==res+cor)? ((m)?res:-res) : csloww1(a,da,x);
++ retval = (res==res+cor)? ((m)?res:-res) : csloww1(a,da,x);
++ goto ret;
+ }
+
+ } /* else if (k < 0x400368fd) */
+@@ -431,7 +462,8 @@ double __cos(double x)
+ res = a+t;
+ cor = (a-res)+t;
+ cor = (cor>0)? 1.02*cor+eps : 1.02*cor -eps;
+- return (res == res + cor)? res : csloww(a,da,x);
++ retval = (res == res + cor)? res : csloww(a,da,x);
++ goto ret;
+ }
+ else {
+ if (a>0) {m=1;t=a;db=da;}
+@@ -450,7 +482,8 @@ double __cos(double x)
+ res=sn+cor;
+ cor=(sn-res)+cor;
+ cor = (cor>0)? 1.035*cor+eps : 1.035*cor-eps;
+- return (res==res+cor)? ((m)?res:-res) : csloww1(a,da,x);
++ retval = (res==res+cor)? ((m)?res:-res) : csloww1(a,da,x);
++ goto ret;
+ }
+ break;
+
+@@ -471,7 +504,8 @@ double __cos(double x)
+ res=cs+cor;
+ cor=(cs-res)+cor;
+ cor = (cor>0)? 1.025*cor+eps : 1.025*cor-eps;
+- return (res==res+cor)? ((n)?-res:res) : csloww2(a,da,x,n);
++ retval = (res==res+cor)? ((n)?-res:res) : csloww2(a,da,x,n);
++ goto ret;
+
+ break;
+
+@@ -506,7 +540,8 @@ double __cos(double x)
+ res = a+t;
+ cor = (a-res)+t;
+ cor = (cor>0)? 1.02*cor+eps : 1.02*cor -eps;
+- return (res == res + cor)? res : bsloww(a,da,x,n);
++ retval = (res == res + cor)? res : bsloww(a,da,x,n);
++ goto ret;
+ }
+ else {
+ if (a>0) {m=1;t=a;db=da;}
+@@ -525,7 +560,8 @@ double __cos(double x)
+ res=sn+cor;
+ cor=(sn-res)+cor;
+ cor = (cor>0)? 1.035*cor+eps : 1.035*cor-eps;
+- return (res==res+cor)? ((m)?res:-res) : bsloww1(a,da,x,n);
++ retval = (res==res+cor)? ((m)?res:-res) : bsloww1(a,da,x,n);
++ goto ret;
+ }
+ break;
+
+@@ -546,7 +582,8 @@ double __cos(double x)
+ res=cs+cor;
+ cor=(cs-res)+cor;
+ cor = (cor>0)? 1.025*cor+eps : 1.025*cor-eps;
+- return (res==res+cor)? ((n)?-res:res) : bsloww2(a,da,x,n);
++ retval = (res==res+cor)? ((n)?-res:res) : bsloww2(a,da,x,n);
++ goto ret;
+ break;
+
+ }
+@@ -558,17 +595,20 @@ double __cos(double x)
+ n = __branred(x,&a,&da);
+ switch (n) {
+ case 1:
+- if (a*a < 0.01588) return bsloww(-a,-da,x,n);
+- else return bsloww1(-a,-da,x,n);
++ if (a*a < 0.01588) retval = bsloww(-a,-da,x,n);
++ else retval = bsloww1(-a,-da,x,n);
++ goto ret;
+ break;
+ case 3:
+- if (a*a < 0.01588) return bsloww(a,da,x,n);
+- else return bsloww1(a,da,x,n);
++ if (a*a < 0.01588) retval = bsloww(a,da,x,n);
++ else retval = bsloww1(a,da,x,n);
++ goto ret;
+ break;
+
+ case 0:
+ case 2:
+- return bsloww2(a,da,x,n);
++ retval = bsloww2(a,da,x,n);
++ goto ret;
+ break;
+ }
+
+@@ -580,10 +620,13 @@ double __cos(double x)
+ else {
+ if (k == 0x7ff00000 && u.i[LOW_HALF] == 0)
+ __set_errno (EDOM);
+- return x / x; /* |x| > 2^1024 */
++ retval = x / x; /* |x| > 2^1024 */
++ goto ret;
+ }
+- return 0;
+
++ ret:
++ feupdateenv (&env);
++ return retval;
+ }
+
+ /************************************************************************/
+--- a/sysdeps/ieee754/dbl-64/s_tan.c
++++ b/sysdeps/ieee754/dbl-64/s_tan.c
+@@ -40,6 +40,8 @@
+ #include "mpa.h"
+ #include "MathLib.h"
+ #include "math.h"
++#include "math_private.h"
++#include <fenv.h>
+
+ static double tanMp(double);
+ void __mptan(double, mp_no *, int);
+@@ -58,21 +60,28 @@ double tan(double x) {
+ mp_no mpy;
+ #endif
+
++ fenv_t env;
++ double retval;
++
+ int __branred(double, double *, double *);
+ int __mpranred(double, mp_no *, int);
+
++ feholdexcept (&env);
++ fesetround (FE_TONEAREST);
++
+ /* x=+-INF, x=NaN */
+ num.d = x; ux = num.i[HIGH_HALF];
+ if ((ux&0x7ff00000)==0x7ff00000) {
+ if ((ux&0x7fffffff)==0x7ff00000)
+ __set_errno (EDOM);
+- return x-x;
++ retval = x-x;
++ goto ret;
+ }
+
+ w=(x<ZERO) ? -x : x;
+
+ /* (I) The case abs(x) <= 1.259e-8 */
+- if (w<=g1.d) return x;
++ if (w<=g1.d) { retval = x; goto ret; }
+
+ /* (II) The case 1.259e-8 < abs(x) <= 0.0608 */
+ if (w<=g2.d) {
+@@ -80,7 +89,7 @@ double tan(double x) {
+ /* First stage */
+ x2 = x*x;
+ t2 = x*x2*(d3.d+x2*(d5.d+x2*(d7.d+x2*(d9.d+x2*d11.d))));
+- if ((y=x+(t2-u1.d*t2)) == x+(t2+u1.d*t2)) return y;
++ if ((y=x+(t2-u1.d*t2)) == x+(t2+u1.d*t2)) { retval = y; goto ret; }
+
+ /* Second stage */
+ c1 = x2*(a15.d+x2*(a17.d+x2*(a19.d+x2*(a21.d+x2*(a23.d+x2*(a25.d+
+@@ -100,8 +109,9 @@ double tan(double x) {
+ MUL2(x2,xx2,c2,cc2,c1,cc1,t1,t2,t3,t4,t5,t6,t7,t8)
+ MUL2(x ,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8)
+ ADD2(x ,zero.d,c2,cc2,c1,cc1,t1,t2)
+- if ((y=c1+(cc1-u2.d*c1)) == c1+(cc1+u2.d*c1)) return y;
+- return tanMp(x);
++ if ((y=c1+(cc1-u2.d*c1)) == c1+(cc1+u2.d*c1)) { retval = y; goto ret; }
++ retval = tanMp(x);
++ goto ret;
+ }
+
+ /* (III) The case 0.0608 < abs(x) <= 0.787 */
+@@ -112,10 +122,10 @@ double tan(double x) {
+ z = w-xfg[i][0].d; z2 = z*z; s = (x<ZERO) ? MONE : ONE;
+ pz = z+z*z2*(e0.d+z2*e1.d);
+ fi = xfg[i][1].d; gi = xfg[i][2].d; t2 = pz*(gi+fi)/(gi-pz);
+- if ((y=fi+(t2-fi*u3.d))==fi+(t2+fi*u3.d)) return (s*y);
++ if ((y=fi+(t2-fi*u3.d))==fi+(t2+fi*u3.d)) { retval = (s*y); goto ret; }
+ t3 = (t2<ZERO) ? -t2 : t2;
+ t4 = fi*ua3.d+t3*ub3.d;
+- if ((y=fi+(t2-t4))==fi+(t2+t4)) return (s*y);
++ if ((y=fi+(t2-t4))==fi+(t2+t4)) { retval = (s*y); goto ret; }
+
+ /* Second stage */
+ ffi = xfg[i][3].d;
+@@ -133,8 +143,9 @@ double tan(double x) {
+ SUB2(one.d,zero.d,c3,cc3,c1,cc1,t1,t2)
+ DIV2(c2,cc2,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
+
+- if ((y=c3+(cc3-u4.d*c3))==c3+(cc3+u4.d*c3)) return (s*y);
+- return tanMp(x);
++ if ((y=c3+(cc3-u4.d*c3))==c3+(cc3+u4.d*c3)) { retval = (s*y); goto ret; }
++ retval = tanMp(x);
++ goto ret;
+ }
+
+ /* (---) The case 0.787 < abs(x) <= 25 */
+@@ -152,20 +163,20 @@ double tan(double x) {
+ else {ya= a; yya= da; sy= ONE;}
+
+ /* (IV),(V) The case 0.787 < abs(x) <= 25, abs(y) <= 1e-7 */
+- if (ya<=gy1.d) return tanMp(x);
++ if (ya<=gy1.d) { retval = tanMp(x); goto ret; }
+
+ /* (VI) The case 0.787 < abs(x) <= 25, 1e-7 < abs(y) <= 0.0608 */
+ if (ya<=gy2.d) {
+ a2 = a*a;
+ t2 = da+a*a2*(d3.d+a2*(d5.d+a2*(d7.d+a2*(d9.d+a2*d11.d))));
+ if (n) {
+- /* First stage -cot */
+- EADD(a,t2,b,db)
+- DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
+- if ((y=c+(dc-u6.d*c))==c+(dc+u6.d*c)) return (-y); }
++ /* First stage -cot */
++ EADD(a,t2,b,db)
++ DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
++ if ((y=c+(dc-u6.d*c))==c+(dc+u6.d*c)) { retval = (-y); goto ret; } }
+ else {
+- /* First stage tan */
+- if ((y=a+(t2-u5.d*a))==a+(t2+u5.d*a)) return y; }
++ /* First stage tan */
++ if ((y=a+(t2-u5.d*a))==a+(t2+u5.d*a)) { retval = y; goto ret; } }
+ /* Second stage */
+ /* Range reduction by algorithm ii */
+ t = (x*hpinv.d + toint.d);
+@@ -201,13 +212,14 @@ double tan(double x) {
+ ADD2(a ,da ,c2,cc2,c1,cc1,t1,t2)
+
+ if (n) {
+- /* Second stage -cot */
+- DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
+- if ((y=c2+(cc2-u8.d*c2)) == c2+(cc2+u8.d*c2)) return (-y); }
++ /* Second stage -cot */
++ DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
++ if ((y=c2+(cc2-u8.d*c2)) == c2+(cc2+u8.d*c2)) { retval = (-y); goto ret; } }
+ else {
+- /* Second stage tan */
+- if ((y=c1+(cc1-u7.d*c1)) == c1+(cc1+u7.d*c1)) return y; }
+- return tanMp(x);
++ /* Second stage tan */
++ if ((y=c1+(cc1-u7.d*c1)) == c1+(cc1+u7.d*c1)) { retval = y; goto ret; } }
++ retval = tanMp(x);
++ goto ret;
+ }
+
+ /* (VII) The case 0.787 < abs(x) <= 25, 0.0608 < abs(y) <= 0.787 */
+@@ -221,17 +233,17 @@ double tan(double x) {
+ if (n) {
+ /* -cot */
+ t2 = pz*(fi+gi)/(fi+pz);
+- if ((y=gi-(t2-gi*u10.d))==gi-(t2+gi*u10.d)) return (-sy*y);
++ if ((y=gi-(t2-gi*u10.d))==gi-(t2+gi*u10.d)) { retval = (-sy*y); goto ret; }
+ t3 = (t2<ZERO) ? -t2 : t2;
+ t4 = gi*ua10.d+t3*ub10.d;
+- if ((y=gi-(t2-t4))==gi-(t2+t4)) return (-sy*y); }
++ if ((y=gi-(t2-t4))==gi-(t2+t4)) { retval = (-sy*y); goto ret; } }
+ else {
+ /* tan */
+ t2 = pz*(gi+fi)/(gi-pz);
+- if ((y=fi+(t2-fi*u9.d))==fi+(t2+fi*u9.d)) return (sy*y);
++ if ((y=fi+(t2-fi*u9.d))==fi+(t2+fi*u9.d)) { retval = (sy*y); goto ret; }
+ t3 = (t2<ZERO) ? -t2 : t2;
+ t4 = fi*ua9.d+t3*ub9.d;
+- if ((y=fi+(t2-t4))==fi+(t2+t4)) return (sy*y); }
++ if ((y=fi+(t2-t4))==fi+(t2+t4)) { retval = (sy*y); goto ret; } }
+
+ /* Second stage */
+ ffi = xfg[i][3].d;
+@@ -252,13 +264,14 @@ double tan(double x) {
+ if (n) {
+ /* -cot */
+ DIV2(c1,cc1,c2,cc2,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
+- if ((y=c3+(cc3-u12.d*c3))==c3+(cc3+u12.d*c3)) return (-sy*y); }
++ if ((y=c3+(cc3-u12.d*c3))==c3+(cc3+u12.d*c3)) { retval = (-sy*y); goto ret; } }
+ else {
+ /* tan */
+ DIV2(c2,cc2,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
+- if ((y=c3+(cc3-u11.d*c3))==c3+(cc3+u11.d*c3)) return (sy*y); }
++ if ((y=c3+(cc3-u11.d*c3))==c3+(cc3+u11.d*c3)) { retval = (sy*y); goto ret; } }
+
+- return tanMp(x);
++ retval = tanMp(x);
++ goto ret;
+ }
+
+ /* (---) The case 25 < abs(x) <= 1e8 */
+@@ -280,20 +293,20 @@ double tan(double x) {
+ else {ya= a; yya= da; sy= ONE;}
+
+ /* (+++) The case 25 < abs(x) <= 1e8, abs(y) <= 1e-7 */
+- if (ya<=gy1.d) return tanMp(x);
++ if (ya<=gy1.d) { retval = tanMp(x); goto ret; }
+
+ /* (VIII) The case 25 < abs(x) <= 1e8, 1e-7 < abs(y) <= 0.0608 */
+ if (ya<=gy2.d) {
+ a2 = a*a;
+ t2 = da+a*a2*(d3.d+a2*(d5.d+a2*(d7.d+a2*(d9.d+a2*d11.d))));
+ if (n) {
+- /* First stage -cot */
+- EADD(a,t2,b,db)
+- DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
+- if ((y=c+(dc-u14.d*c))==c+(dc+u14.d*c)) return (-y); }
++ /* First stage -cot */
++ EADD(a,t2,b,db)
++ DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
++ if ((y=c+(dc-u14.d*c))==c+(dc+u14.d*c)) { retval = (-y); goto ret; } }
+ else {
+- /* First stage tan */
+- if ((y=a+(t2-u13.d*a))==a+(t2+u13.d*a)) return y; }
++ /* First stage tan */
++ if ((y=a+(t2-u13.d*a))==a+(t2+u13.d*a)) { retval = y; goto ret; } }
+
+ /* Second stage */
+ MUL2(a,da,a,da,x2,xx2,t1,t2,t3,t4,t5,t6,t7,t8)
+@@ -315,13 +328,14 @@ double tan(double x) {
+ ADD2(a ,da ,c2,cc2,c1,cc1,t1,t2)
+
+ if (n) {
+- /* Second stage -cot */
+- DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
+- if ((y=c2+(cc2-u16.d*c2)) == c2+(cc2+u16.d*c2)) return (-y); }
++ /* Second stage -cot */
++ DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
++ if ((y=c2+(cc2-u16.d*c2)) == c2+(cc2+u16.d*c2)) { retval = (-y); goto ret; } }
+ else {
+- /* Second stage tan */
+- if ((y=c1+(cc1-u15.d*c1)) == c1+(cc1+u15.d*c1)) return (y); }
+- return tanMp(x);
++ /* Second stage tan */
++ if ((y=c1+(cc1-u15.d*c1)) == c1+(cc1+u15.d*c1)) { retval = (y); goto ret; } }
++ retval = tanMp(x);
++ goto ret;
+ }
+
+ /* (IX) The case 25 < abs(x) <= 1e8, 0.0608 < abs(y) <= 0.787 */
+@@ -334,17 +348,17 @@ double tan(double x) {
+ if (n) {
+ /* -cot */
+ t2 = pz*(fi+gi)/(fi+pz);
+- if ((y=gi-(t2-gi*u18.d))==gi-(t2+gi*u18.d)) return (-sy*y);
++ if ((y=gi-(t2-gi*u18.d))==gi-(t2+gi*u18.d)) { retval = (-sy*y); goto ret; }
+ t3 = (t2<ZERO) ? -t2 : t2;
+ t4 = gi*ua18.d+t3*ub18.d;
+- if ((y=gi-(t2-t4))==gi-(t2+t4)) return (-sy*y); }
++ if ((y=gi-(t2-t4))==gi-(t2+t4)) { retval = (-sy*y); goto ret; } }
+ else {
+ /* tan */
+ t2 = pz*(gi+fi)/(gi-pz);
+- if ((y=fi+(t2-fi*u17.d))==fi+(t2+fi*u17.d)) return (sy*y);
++ if ((y=fi+(t2-fi*u17.d))==fi+(t2+fi*u17.d)) { retval = (sy*y); goto ret; }
+ t3 = (t2<ZERO) ? -t2 : t2;
+ t4 = fi*ua17.d+t3*ub17.d;
+- if ((y=fi+(t2-t4))==fi+(t2+t4)) return (sy*y); }
++ if ((y=fi+(t2-t4))==fi+(t2+t4)) { retval = (sy*y); goto ret; } }
+
+ /* Second stage */
+ ffi = xfg[i][3].d;
+@@ -365,12 +379,13 @@ double tan(double x) {
+ if (n) {
+ /* -cot */
+ DIV2(c1,cc1,c2,cc2,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
+- if ((y=c3+(cc3-u20.d*c3))==c3+(cc3+u20.d*c3)) return (-sy*y); }
++ if ((y=c3+(cc3-u20.d*c3))==c3+(cc3+u20.d*c3)) { retval = (-sy*y); goto ret; } }
+ else {
+ /* tan */
+ DIV2(c2,cc2,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
+- if ((y=c3+(cc3-u19.d*c3))==c3+(cc3+u19.d*c3)) return (sy*y); }
+- return tanMp(x);
++ if ((y=c3+(cc3-u19.d*c3))==c3+(cc3+u19.d*c3)) { retval = (sy*y); goto ret; } }
++ retval = tanMp(x);
++ goto ret;
+ }
+
+ /* (---) The case 1e8 < abs(x) < 2**1024 */
+@@ -381,7 +396,7 @@ double tan(double x) {
+ else {ya= a; yya= da; sy= ONE;}
+
+ /* (+++) The case 1e8 < abs(x) < 2**1024, abs(y) <= 1e-7 */
+- if (ya<=gy1.d) return tanMp(x);
++ if (ya<=gy1.d) { retval = tanMp(x); goto ret; }
+
+ /* (X) The case 1e8 < abs(x) < 2**1024, 1e-7 < abs(y) <= 0.0608 */
+ if (ya<=gy2.d) {
+@@ -391,10 +406,10 @@ double tan(double x) {
+ /* First stage -cot */
+ EADD(a,t2,b,db)
+ DIV2(one.d,zero.d,b,db,c,dc,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
+- if ((y=c+(dc-u22.d*c))==c+(dc+u22.d*c)) return (-y); }
++ if ((y=c+(dc-u22.d*c))==c+(dc+u22.d*c)) { retval = (-y); goto ret; } }
+ else {
+ /* First stage tan */
+- if ((y=a+(t2-u21.d*a))==a+(t2+u21.d*a)) return y; }
++ if ((y=a+(t2-u21.d*a))==a+(t2+u21.d*a)) { retval = y; goto ret; } }
+
+ /* Second stage */
+ /* Reduction by algorithm iv */
+@@ -423,11 +438,12 @@ double tan(double x) {
+ if (n) {
+ /* Second stage -cot */
+ DIV2(one.d,zero.d,c1,cc1,c2,cc2,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
+- if ((y=c2+(cc2-u24.d*c2)) == c2+(cc2+u24.d*c2)) return (-y); }
++ if ((y=c2+(cc2-u24.d*c2)) == c2+(cc2+u24.d*c2)) { retval = (-y); goto ret; } }
+ else {
+ /* Second stage tan */
+- if ((y=c1+(cc1-u23.d*c1)) == c1+(cc1+u23.d*c1)) return y; }
+- return tanMp(x);
++ if ((y=c1+(cc1-u23.d*c1)) == c1+(cc1+u23.d*c1)) { retval = y; goto ret; } }
++ retval = tanMp(x);
++ goto ret;
+ }
+
+ /* (XI) The case 1e8 < abs(x) < 2**1024, 0.0608 < abs(y) <= 0.787 */
+@@ -440,17 +456,17 @@ double tan(double x) {
+ if (n) {
+ /* -cot */
+ t2 = pz*(fi+gi)/(fi+pz);
+- if ((y=gi-(t2-gi*u26.d))==gi-(t2+gi*u26.d)) return (-sy*y);
++ if ((y=gi-(t2-gi*u26.d))==gi-(t2+gi*u26.d)) { retval = (-sy*y); goto ret; }
+ t3 = (t2<ZERO) ? -t2 : t2;
+ t4 = gi*ua26.d+t3*ub26.d;
+- if ((y=gi-(t2-t4))==gi-(t2+t4)) return (-sy*y); }
++ if ((y=gi-(t2-t4))==gi-(t2+t4)) { retval = (-sy*y); goto ret; } }
+ else {
+ /* tan */
+ t2 = pz*(gi+fi)/(gi-pz);
+- if ((y=fi+(t2-fi*u25.d))==fi+(t2+fi*u25.d)) return (sy*y);
++ if ((y=fi+(t2-fi*u25.d))==fi+(t2+fi*u25.d)) { retval = (sy*y); goto ret; }
+ t3 = (t2<ZERO) ? -t2 : t2;
+ t4 = fi*ua25.d+t3*ub25.d;
+- if ((y=fi+(t2-t4))==fi+(t2+t4)) return (sy*y); }
++ if ((y=fi+(t2-t4))==fi+(t2+t4)) { retval = (sy*y); goto ret; } }
+
+ /* Second stage */
+ ffi = xfg[i][3].d;
+@@ -471,15 +487,19 @@ double tan(double x) {
+ if (n) {
+ /* -cot */
+ DIV2(c1,cc1,c2,cc2,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
+- if ((y=c3+(cc3-u28.d*c3))==c3+(cc3+u28.d*c3)) return (-sy*y); }
++ if ((y=c3+(cc3-u28.d*c3))==c3+(cc3+u28.d*c3)) { retval = (-sy*y); goto ret; } }
+ else {
+ /* tan */
+ DIV2(c2,cc2,c1,cc1,c3,cc3,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10)
+- if ((y=c3+(cc3-u27.d*c3))==c3+(cc3+u27.d*c3)) return (sy*y); }
+- return tanMp(x);
++ if ((y=c3+(cc3-u27.d*c3))==c3+(cc3+u27.d*c3)) { retval = (sy*y); goto ret; } }
++ retval = tanMp(x);
++ goto ret;
++
++ ret:
++ feupdateenv (&env);
++ return retval;
+ }
+
+-
+ /* multiple precision stage */
+ /* Convert x to multi precision number,compute tan(x) by mptan() routine */
+ /* and converts result back to double */
+--- a/sysdeps/x86_64/fpu/libm-test-ulps
++++ b/sysdeps/x86_64/fpu/libm-test-ulps
+@@ -327,6 +327,107 @@ Test "cos (0.80190127184058835) == 0.695
+ double: 1
+ idouble: 1
+
++# cos_downward
++Test "cos_downward (1) == 0.5403023058681397174009366074429766037323":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "cos_downward (2) == -0.4161468365471423869975682295007621897660":
++float: 1
++ifloat: 1
++Test "cos_downward (3) == -0.9899924966004454572715727947312613023937":
++float: 1
++ifloat: 1
++Test "cos_downward (4) == -0.6536436208636119146391681830977503814241":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "cos_downward (5) == 0.2836621854632262644666391715135573083344":
++float: 1
++ifloat: 1
++Test "cos_downward (7) == 0.7539022543433046381411975217191820122183":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "cos_downward (8) == -0.1455000338086135258688413818311946826093":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "cos_downward (9) == -0.9111302618846769883682947111811653112463":
++ildouble: 1
++ldouble: 1
++
++# cos_tonearest
++Test "cos_tonearest (7) == 0.7539022543433046381411975217191820122183":
++float: 1
++ifloat: 1
++Test "cos_tonearest (8) == -0.1455000338086135258688413818311946826093":
++ildouble: 1
++ldouble: 1
++Test "cos_tonearest (9) == -0.9111302618846769883682947111811653112463":
++ildouble: 1
++ldouble: 1
++
++# cos_towardzero
++Test "cos_towardzero (1) == 0.5403023058681397174009366074429766037323":
++ildouble: 1
++ldouble: 1
++Test "cos_towardzero (10) == -0.8390715290764524522588639478240648345199":
++ildouble: 1
++ldouble: 1
++Test "cos_towardzero (2) == -0.4161468365471423869975682295007621897660":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "cos_towardzero (3) == -0.9899924966004454572715727947312613023937":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "cos_towardzero (5) == 0.2836621854632262644666391715135573083344":
++float: 1
++ifloat: 1
++Test "cos_towardzero (7) == 0.7539022543433046381411975217191820122183":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "cos_towardzero (8) == -0.1455000338086135258688413818311946826093":
++float: 1
++ifloat: 1
++
++# cos_upward
++Test "cos_upward (10) == -0.8390715290764524522588639478240648345199":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "cos_upward (2) == -0.4161468365471423869975682295007621897660":
++ildouble: 1
++ldouble: 1
++Test "cos_upward (3) == -0.9899924966004454572715727947312613023937":
++ildouble: 1
++ldouble: 1
++Test "cos_upward (5) == 0.2836621854632262644666391715135573083344":
++ildouble: 1
++ldouble: 1
++Test "cos_upward (6) == 0.9601702866503660205456522979229244054519":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "cos_upward (7) == 0.7539022543433046381411975217191820122183":
++float: 1
++ifloat: 1
++Test "cos_upward (9) == -0.9111302618846769883682947111811653112463":
++float: 2
++ifloat: 2
++
+ # cpow
+ Test "Real part of: cpow (0.75 + 1.25 i, 0.0 + 1.0 i) == 0.331825439177608832276067945276730566 + 0.131338600281188544930936345230903032 i":
+ float: 1
+@@ -766,6 +867,107 @@ Test "log1p (-0.25) == -0.28768207245178
+ float: 1
+ ifloat: 1
+
++# sin_downward
++Test "sin_downward (1) == 0.8414709848078965066525023216302989996226":
++ildouble: 1
++ldouble: 1
++Test "sin_downward (10) == -0.5440211108893698134047476618513772816836":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "sin_downward (3) == 0.1411200080598672221007448028081102798469":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "sin_downward (4) == -0.7568024953079282513726390945118290941359":
++ildouble: 1
++ldouble: 1
++Test "sin_downward (5) == -0.9589242746631384688931544061559939733525":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "sin_downward (6) == -0.2794154981989258728115554466118947596280":
++float: 1
++ifloat: 1
++Test "sin_downward (7) == 0.6569865987187890903969990915936351779369":
++ildouble: 1
++ldouble: 1
++Test "sin_downward (8) == 0.9893582466233817778081235982452886721164":
++ildouble: 1
++ldouble: 1
++Test "sin_downward (9) == 0.4121184852417565697562725663524351793439":
++ildouble: 1
++ldouble: 1
++
++# sin_tonearest
++Test "sin_tonearest (1) == 0.8414709848078965066525023216302989996226":
++float: 1
++ifloat: 1
++Test "sin_tonearest (10) == -0.5440211108893698134047476618513772816836":
++ildouble: 1
++ldouble: 1
++Test "sin_tonearest (4) == -0.7568024953079282513726390945118290941359":
++ildouble: 1
++ldouble: 1
++Test "sin_tonearest (9) == 0.4121184852417565697562725663524351793439":
++ildouble: 1
++ldouble: 1
++
++# sin_towardzero
++Test "sin_towardzero (1) == 0.8414709848078965066525023216302989996226":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "sin_towardzero (10) == -0.5440211108893698134047476618513772816836":
++float: 1
++ifloat: 1
++Test "sin_towardzero (3) == 0.1411200080598672221007448028081102798469":
++ildouble: 1
++ldouble: 1
++Test "sin_towardzero (4) == -0.7568024953079282513726390945118290941359":
++float: 1
++ifloat: 1
++Test "sin_towardzero (5) == -0.9589242746631384688931544061559939733525":
++float: 1
++ifloat: 1
++Test "sin_towardzero (6) == -0.2794154981989258728115554466118947596280":
++ildouble: 1
++ldouble: 1
++Test "sin_towardzero (7) == 0.6569865987187890903969990915936351779369":
++ildouble: 1
++ldouble: 1
++Test "sin_towardzero (8) == 0.9893582466233817778081235982452886721164":
++ildouble: 1
++ldouble: 1
++Test "sin_towardzero (9) == 0.4121184852417565697562725663524351793439":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++# sin_upward
++Test "sin_upward (1) == 0.8414709848078965066525023216302989996226":
++float: 1
++ifloat: 1
++Test "sin_upward (2) == 0.9092974268256816953960198659117448427023":
++float: 2
++ifloat: 2
++ildouble: 1
++ldouble: 1
++Test "sin_upward (4) == -0.7568024953079282513726390945118290941359":
++float: 1
++ifloat: 1
++Test "sin_upward (6) == -0.2794154981989258728115554466118947596280":
++ildouble: 1
++ldouble: 1
++Test "sin_upward (9) == 0.4121184852417565697562725663524351793439":
++float: 1
++ifloat: 1
++
+ # sincos
+ Test "sincos (M_PI_6l*2.0, &sin_res, &cos_res) puts 0.5 in cos_res":
+ double: 1
+@@ -798,6 +1000,112 @@ Test "tan (pi/4) == 1":
+ double: 1
+ idouble: 1
+
++# tan_downward
++Test "tan_downward (1) == 1.5574077246549022305069748074583601730873":
++float: 1
++ifloat: 1
++Test "tan_downward (10) == 0.6483608274590866712591249330098086768169":
++float: 1
++ifloat: 1
++Test "tan_downward (2) == -2.1850398632615189916433061023136825434320":
++float: 1
++ifloat: 1
++Test "tan_downward (4) == 1.1578212823495775831373424182673239231198":
++ildouble: 1
++ldouble: 1
++Test "tan_downward (5) == -3.3805150062465856369827058794473439087096":
++ildouble: 1
++ldouble: 1
++Test "tan_downward (6) == -0.2910061913847491570536995888681755428312":
++float: 1
++ifloat: 1
++Test "tan_downward (8) == -6.7997114552203786999252627596086333648814":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "tan_downward (9) == -0.4523156594418098405903708757987855343087":
++float: 1
++ifloat: 1
++
++# tan_tonearest
++Test "tan_tonearest (6) == -0.2910061913847491570536995888681755428312":
++ildouble: 1
++ldouble: 1
++Test "tan_tonearest (8) == -6.7997114552203786999252627596086333648814":
++ildouble: 1
++ldouble: 1
++Test "tan_tonearest (9) == -0.4523156594418098405903708757987855343087":
++ildouble: 1
++ldouble: 1
++
++# tan_towardzero
++Test "tan_towardzero (10) == 0.6483608274590866712591249330098086768169":
++float: 1
++ifloat: 1
++Test "tan_towardzero (2) == -2.1850398632615189916433061023136825434320":
++ildouble: 1
++ldouble: 1
++Test "tan_towardzero (3) == -0.1425465430742778052956354105339134932261":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "tan_towardzero (4) == 1.1578212823495775831373424182673239231198":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "tan_towardzero (5) == -3.3805150062465856369827058794473439087096":
++float: 1
++ifloat: 1
++Test "tan_towardzero (6) == -0.2910061913847491570536995888681755428312":
++ildouble: 1
++ldouble: 1
++Test "tan_towardzero (8) == -6.7997114552203786999252627596086333648814":
++ildouble: 2
++ldouble: 2
++Test "tan_towardzero (9) == -0.4523156594418098405903708757987855343087":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++# tan_upward
++Test "tan_upward (1) == 1.5574077246549022305069748074583601730873":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "tan_upward (10) == 0.6483608274590866712591249330098086768169":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "tan_upward (2) == -2.1850398632615189916433061023136825434320":
++ildouble: 1
++ldouble: 1
++Test "tan_upward (3) == -0.1425465430742778052956354105339134932261":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++Test "tan_upward (5) == -3.3805150062465856369827058794473439087096":
++float: 1
++ifloat: 1
++Test "tan_upward (6) == -0.2910061913847491570536995888681755428312":
++ildouble: 1
++ldouble: 1
++Test "tan_upward (7) == 0.8714479827243187364564508896003135663222":
++ildouble: 1
++ldouble: 1
++Test "tan_upward (8) == -6.7997114552203786999252627596086333648814":
++ildouble: 2
++ldouble: 2
++Test "tan_upward (9) == -0.4523156594418098405903708757987855343087":
++ildouble: 1
++ldouble: 1
++
+ # tgamma
+ Test "tgamma (-0.5) == -2 sqrt (pi)":
+ double: 1
+@@ -1144,6 +1452,30 @@ ifloat: 1
+ ildouble: 1
+ ldouble: 1
+
++Function: "cos_downward":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++Function: "cos_tonearest":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++Function: "cos_towardzero":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++Function: "cos_upward":
++float: 2
++ifloat: 2
++ildouble: 1
++ldouble: 1
++
+ Function: Real part of "cpow":
+ double: 2
+ float: 5
+@@ -1312,6 +1644,30 @@ Function: "log1p":
+ float: 1
+ ifloat: 1
+
++Function: "sin_downward":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++Function: "sin_tonearest":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++Function: "sin_towardzero":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++Function: "sin_upward":
++float: 2
++ifloat: 2
++ildouble: 1
++ldouble: 1
++
+ Function: "sincos":
+ double: 1
+ float: 1
+@@ -1324,6 +1680,28 @@ Function: "tan":
+ double: 1
+ idouble: 1
+
++Function: "tan_downward":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++Function: "tan_tonearest":
++ildouble: 1
++ldouble: 1
++
++Function: "tan_towardzero":
++float: 1
++ifloat: 1
++ildouble: 2
++ldouble: 2
++
++Function: "tan_upward":
++float: 1
++ifloat: 1
++ildouble: 2
++ldouble: 2
++
+ Function: "tgamma":
+ double: 1
+ float: 1
only in patch2:
unchanged:
--- eglibc-2.13.orig/debian/patches/any/cvs-sin-accuracy.diff
+++ eglibc-2.13/debian/patches/any/cvs-sin-accuracy.diff
@@ -0,0 +1,33 @@
+2011-07-03 Andreas Jaeger <aj@suse.de>
+
+ [BZ #10709]
+ * sysdeps/ieee754/dbl-64/s_sin.c (__sin): Fix incorrect rounding
+ of sin. Patch suggested by Paul Zimmermann <zimmerma+gcc@loria.fr>.
+ * math/libm-test.inc (sin_test): Add test case.
+---
+ ChangeLog | 7 +++++++
+ math/libm-test.inc | 1 +
+ sysdeps/ieee754/dbl-64/s_sin.c | 2 +-
+ 3 files changed, 9 insertions(+), 1 deletion(-)
+
+--- a/math/libm-test.inc
++++ b/math/libm-test.inc
+@@ -5736,6 +5736,7 @@ sin_test (void)
+
+ #ifdef TEST_DOUBLE
+ TEST_f_f (sin, 0.80190127184058835, 0.71867942238767868);
++ TEST_f_f (sin, 2.522464e-1, 2.4957989804940911e-1);
+ #endif
+
+ END (sin);
+--- a/sysdeps/ieee754/dbl-64/s_sin.c
++++ b/sysdeps/ieee754/dbl-64/s_sin.c
+@@ -127,7 +127,7 @@ double __sin(double x){
+ cor=(ssn+s*ccs-sn*c)+cs*s;
+ res=sn+cor;
+ cor=(sn-res)+cor;
+- return (res==res+1.025*cor)? res : slow1(x);
++ return (res==res+1.096*cor)? res : slow1(x);
+ } /* else if (k < 0x3feb6000) */
+
+ /*----------------------- 0.855469 <|x|<2.426265 ----------------------*/
only in patch2:
unchanged:
--- eglibc-2.13.orig/debian/patches/any/cvs-exp-rounding-mode.diff
+++ eglibc-2.13/debian/patches/any/cvs-exp-rounding-mode.diff
@@ -0,0 +1,392 @@
+2012-03-02 Joseph Myers <joseph@codesourcery.com>
+
+ [BZ #3976]
+ * sysdeps/ieee754/dbl-64/e_exp.c: Include <fenv.h>.
+ (__ieee754_exp): Save and restore rounding mode and use
+ round-to-nearest for all computations.
+ * math/libm-test.inc (exp_test_tonearest): New function.
+ (exp_test_towardzero): Likewise.
+ (exp_test_downward): Likewise.
+ (exp_test_upward): Likewise.
+ (main): Call the new functions.
+ * sysdeps/i386/fpu/libm-test-ulps: Update.
+ * sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
+---
+ ChangeLog | 14 +++++
+ math/libm-test.inc | 112 ++++++++++++++++++++++++++++++++++++++
+ sysdeps/i386/fpu/libm-test-ulps | 67 +++++++++++++++++++++++
+ sysdeps/ieee754/dbl-64/e_exp.c | 38 ++++++++-----
+ sysdeps/x86_64/fpu/libm-test-ulps | 51 +++++++++++++++++
+ 5 files changed, 268 insertions(+), 14 deletions(-)
+
+--- a/math/libm-test.inc
++++ b/math/libm-test.inc
+@@ -2540,6 +2540,114 @@ exp_test (void)
+
+
+ static void
++exp_test_tonearest (void)
++{
++ int save_round_mode;
++ errno = 0;
++ FUNC(exp) (0);
++ if (errno == ENOSYS)
++ /* Function not implemented. */
++ return;
++
++ START (exp_tonearest);
++
++ save_round_mode = fegetround ();
++
++ if (!fesetround (FE_TONEAREST))
++ {
++ TEST_f_f (exp, 1, M_El);
++ TEST_f_f (exp, 2, M_E2l);
++ TEST_f_f (exp, 3, M_E3l);
++ }
++
++ fesetround (save_round_mode);
++
++ END (exp_tonearest);
++}
++
++
++static void
++exp_test_towardzero (void)
++{
++ int save_round_mode;
++ errno = 0;
++ FUNC(exp) (0);
++ if (errno == ENOSYS)
++ /* Function not implemented. */
++ return;
++
++ START (exp_towardzero);
++
++ save_round_mode = fegetround ();
++
++ if (!fesetround (FE_TOWARDZERO))
++ {
++ TEST_f_f (exp, 1, M_El);
++ TEST_f_f (exp, 2, M_E2l);
++ TEST_f_f (exp, 3, M_E3l);
++ }
++
++ fesetround (save_round_mode);
++
++ END (exp_towardzero);
++}
++
++
++static void
++exp_test_downward (void)
++{
++ int save_round_mode;
++ errno = 0;
++ FUNC(exp) (0);
++ if (errno == ENOSYS)
++ /* Function not implemented. */
++ return;
++
++ START (exp_downward);
++
++ save_round_mode = fegetround ();
++
++ if (!fesetround (FE_DOWNWARD))
++ {
++ TEST_f_f (exp, 1, M_El);
++ TEST_f_f (exp, 2, M_E2l);
++ TEST_f_f (exp, 3, M_E3l);
++ }
++
++ fesetround (save_round_mode);
++
++ END (exp_downward);
++}
++
++
++static void
++exp_test_upward (void)
++{
++ int save_round_mode;
++ errno = 0;
++ FUNC(exp) (0);
++ if (errno == ENOSYS)
++ /* Function not implemented. */
++ return;
++
++ START (exp_upward);
++
++ save_round_mode = fegetround ();
++
++ if (!fesetround (FE_UPWARD))
++ {
++ TEST_f_f (exp, 1, M_El);
++ TEST_f_f (exp, 2, M_E2l);
++ TEST_f_f (exp, 3, M_E3l);
++ }
++
++ fesetround (save_round_mode);
++
++ END (exp_upward);
++}
++
++
++static void
+ exp10_test (void)
+ {
+ errno = 0;
+@@ -6268,6 +6376,10 @@ main (int argc, char **argv)
+
+ /* Exponential and logarithmic functions: */
+ exp_test ();
++ exp_test_tonearest ();
++ exp_test_towardzero ();
++ exp_test_downward ();
++ exp_test_upward ();
+ exp10_test ();
+ exp2_test ();
+ expm1_test ();
+--- a/sysdeps/i386/fpu/libm-test-ulps
++++ b/sysdeps/i386/fpu/libm-test-ulps
+@@ -453,6 +453,51 @@ Test "exp10 (3) == 1000":
+ ildouble: 8
+ ldouble: 8
+
++# exp_downward
++Test "exp_downward (1) == e":
++ildouble: 1
++ldouble: 1
++Test "exp_downward (2) == e^2":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 2
++ldouble: 2
++Test "exp_downward (3) == e^3":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++# exp_towardzero
++Test "exp_towardzero (1) == e":
++ildouble: 1
++ldouble: 1
++Test "exp_towardzero (2) == e^2":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 2
++ldouble: 2
++Test "exp_towardzero (3) == e^3":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++# exp_upward
++Test "exp_upward (1) == e":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++
+ # expm1
+ Test "expm1 (1) == M_El - 1.0":
+ ildouble: 1
+@@ -1138,6 +1183,28 @@ Function: "exp10":
+ ildouble: 8
+ ldouble: 8
+
++Function: "exp_downward":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 2
++ldouble: 2
++
++Function: "exp_towardzero":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++ildouble: 2
++ldouble: 2
++
++Function: "exp_upward":
++double: 1
++float: 1
++idouble: 1
++ifloat: 1
++
+ Function: "expm1":
+ ildouble: 1
+
+--- a/sysdeps/ieee754/dbl-64/e_exp.c
++++ b/sysdeps/ieee754/dbl-64/e_exp.c
+@@ -1,7 +1,7 @@
+ /*
+ * IBM Accurate Mathematical Library
+ * written by International Business Machines Corp.
+- * Copyright (C) 2001 Free Software Foundation
++ * Copyright (C) 2001, 2012 Free Software Foundation
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+@@ -39,6 +39,7 @@
+ #include "MathLib.h"
+ #include "uexp.tbl"
+ #include "math_private.h"
++#include <fenv.h>
+
+ double __slowexp(double);
+
+@@ -53,6 +54,11 @@ double __ieee754_exp(double x) {
+ int4 k;
+ #endif
+ int4 i,j,m,n,ex;
++ fenv_t env;
++ double retval;
++
++ feholdexcept (&env);
++ fesetround (FE_TONEAREST);
+
+ junk1.x = x;
+ m = junk1.i[HIGH_HALF];
+@@ -85,18 +91,19 @@ double __ieee754_exp(double x) {
+ rem=(bet + bet*eps)+al*eps;
+ res = al + rem;
+ cor = (al - res) + rem;
+- if (res == (res+cor*err_0)) return res*binexp.x;
+- else return __slowexp(x); /*if error is over bound */
++ if (res == (res+cor*err_0)) { retval = res*binexp.x; goto ret; }
++ else { retval = __slowexp(x); goto ret; } /*if error is over bound */
+ }
+
+- if (n <= smallint) return 1.0;
++ if (n <= smallint) { retval = 1.0; goto ret; }
+
+ if (n >= badint) {
+- if (n > infint) return(x+x); /* x is NaN */
+- if (n < infint) return ( (x>0) ? (hhuge*hhuge) : (tiny*tiny) );
++ if (n > infint) { retval = x+x; goto ret; } /* x is NaN */
++ if (n < infint) { retval = (x>0) ? (hhuge*hhuge) : (tiny*tiny); goto ret; }
+ /* x is finite, cause either overflow or underflow */
+- if (junk1.i[LOW_HALF] != 0) return (x+x); /* x is NaN */
+- return ((x>0)?inf.x:zero ); /* |x| = inf; return either inf or 0 */
++ if (junk1.i[LOW_HALF] != 0) { retval = x+x; goto ret; } /* x is NaN */
++ retval = (x>0)?inf.x:zero; /* |x| = inf; return either inf or 0 */
++ goto ret;
+ }
+
+ y = x*log2e.x + three51.x;
+@@ -121,8 +128,8 @@ double __ieee754_exp(double x) {
+ if (res < 1.0) {res+=res; cor+=cor; ex-=1;}
+ if (ex >=-1022) {
+ binexp.i[HIGH_HALF] = (1023+ex)<<20;
+- if (res == (res+cor*err_0)) return res*binexp.x;
+- else return __slowexp(x); /*if error is over bound */
++ if (res == (res+cor*err_0)) { retval = res*binexp.x; goto ret; }
++ else { retval = __slowexp(x); goto ret; } /*if error is over bound */
+ }
+ ex = -(1022+ex);
+ binexp.i[HIGH_HALF] = (1023-ex)<<20;
+@@ -135,15 +142,19 @@ double __ieee754_exp(double x) {
+ cor = (t-res)+y;
+ if (res == (res + eps*cor))
+ { binexp.i[HIGH_HALF] = 0x00100000;
+- return (res-1.0)*binexp.x;
++ retval = (res-1.0)*binexp.x;
++ goto ret;
+ }
+- else return __slowexp(x); /* if error is over bound */
++ else { retval = __slowexp(x); goto ret; } /* if error is over bound */
+ }
+ else {
+ binexp.i[HIGH_HALF] =(junk1.i[LOW_HALF]+767)<<20;
+- if (res == (res+cor*err_0)) return res*binexp.x*t256.x;
+- else return __slowexp(x);
++ if (res == (res+cor*err_0)) { retval = res*binexp.x*t256.x; goto ret; }
++ else { retval = __slowexp(x); goto ret; }
+ }
++ ret:
++ feupdateenv (&env);
++ return retval;
+ }
+
+ /************************************************************************/
+--- a/sysdeps/x86_64/fpu/libm-test-ulps
++++ b/sysdeps/x86_64/fpu/libm-test-ulps
+@@ -502,6 +502,41 @@ ifloat: 2
+ double: 6
+ idouble: 6
+
++# exp_downward
++Test "exp_downward (1) == e":
++ildouble: 1
++ldouble: 1
++Test "exp_downward (2) == e^2":
++float: 1
++ifloat: 1
++ildouble: 2
++ldouble: 2
++Test "exp_downward (3) == e^3":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++# exp_towardzero
++Test "exp_towardzero (1) == e":
++ildouble: 1
++ldouble: 1
++Test "exp_towardzero (2) == e^2":
++float: 1
++ifloat: 1
++ildouble: 2
++ldouble: 2
++Test "exp_towardzero (3) == e^3":
++float: 1
++ifloat: 1
++ildouble: 1
++ldouble: 1
++
++# exp_upward
++Test "exp_upward (1) == e":
++float: 1
++ifloat: 1
++
+ # expm1
+ Test "expm1 (0.75) == 1.11700001661267466854536981983709561":
+ double: 1
+@@ -1203,6 +1238,22 @@ ifloat: 2
+ double: 6
+ idouble: 6
+
++Function: "exp_downward":
++float: 1
++ifloat: 1
++ildouble: 2
++ldouble: 2
++
++Function: "exp_towardzero":
++float: 1
++ifloat: 1
++ildouble: 2
++ldouble: 2
++
++Function: "exp_upward":
++float: 1
++ifloat: 1
++
+ Function: "expm1":
+ double: 1
+ float: 1
Reply to: