outline details from C18 annex F

49 files changed, 581 insertions, 14 deletions

diff --git a/src/math/acos.c b/src/math/acos.c
index 2f740d03..9a356100 100644
--- a/src/math/acos.c
+++ b/src/math/acos.c
@@ -1,12 +1,15 @@
 # define TGSOURCE "acos.c"
 #include <math.h>
+#include "fenv.h"
 #include "errno.h"
 #include "_tgmath.h"
 
 /** arc cosine **/
 TYPE TGFN(acos)(TYPE x)
 {
-	if (x < -1 || x > 1) {
+	if (TGFN(fabs)(x) > 1) {
+		feraiseexcept(FE_INVALID);
+		return NAN;
 		errno = EDOM;	/* ARGUMENT(x) not in the range [-1, +1] */
 		return TGHUGE;
 	}
@@ -17,6 +20,10 @@ TYPE TGFN(acos)(TYPE x)
 		return TGHUGE;
 	}
 
+	if (x == 1.0) {
+		return 0.0;
+	}
+
 	/* RETURN_SUCCESS(a value in range `[0, PI()]'); */
 	return (TYPE)0.0;
 }
diff --git a/src/math/acosh.c b/src/math/acosh.c
index aeb15816..f557c875 100644
--- a/src/math/acosh.c
+++ b/src/math/acosh.c
@@ -1,9 +1,23 @@
 # define TGSOURCE "acosh.c"
 #include "_tgmath.h"
 #include <math.h>
+#include "fenv.h"
 
 TYPE TGFN(acosh)(TYPE x)
 {
+	if (x == 1.0) {
+		return 0.0;
+	}
+
+	if (x < 1.0) {
+		feraiseexcept(FE_INVALID);
+		return NAN;
+	}
+
+	if (fpclassify(x) == FP_INFINITE && !signbit(x)) {
+		return x;
+	}
+
 	return x;
 }
 
diff --git a/src/math/asin.c b/src/math/asin.c
index 92dae8c8..a94338f5 100644
--- a/src/math/asin.c
+++ b/src/math/asin.c
@@ -6,11 +6,15 @@
 /** arc sine **/
 TYPE TGFN(asin)(TYPE x)
 {
-	if (x < -1.0 || x > 1.0) {
+	if (TGFN(fabs)(x) > 1.0) {
 		errno = EDOM;	/* ARGUMENT(x) is not in the range [-1, +1] */
 		return TGHUGE;
 	}
 
+	if (fpclassify(x) == FP_ZERO) {
+		return x;
+	}
+
 	if (0) {
 		errno = ERANGE;	/* The result cannot be represented */
 		/* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */
diff --git a/src/math/asinh.c b/src/math/asinh.c
index 421f574d..9f4589e6 100644
--- a/src/math/asinh.c
+++ b/src/math/asinh.c
@@ -4,6 +4,12 @@
 
 TYPE TGFN(asinh)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:		return x;
+	case FP_INFINITE:	return x;
+	default:		break;
+	}
+
 	return x;
 }
 
diff --git a/src/math/atan.c b/src/math/atan.c
index 141367cb..6d9eee63 100644
--- a/src/math/atan.c
+++ b/src/math/atan.c
@@ -3,6 +3,8 @@
 #include "errno.h"
 #include "_tgmath.h"
 
+#include "M_PI_2.c"
+
 /** arc tangent **/
 TYPE TGFN(atan)(TYPE x)
 {
@@ -11,6 +13,12 @@ TYPE TGFN(atan)(TYPE x)
 	TYPE power = 1.0;
 	int i;
 
+	switch (fpclassify(x)) {
+	case FP_ZERO:		return x;
+	case FP_INFINITE:	return signbit(x) ? - M_PI_2 : M_PI_2;
+	default:		break;
+	}
+
 	if (0) {
 		errno = ERANGE; /* The result cannot be represented */
 		/* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */
diff --git a/src/math/atan2.c b/src/math/atan2.c
index c4c56a63..3b88b617 100644
--- a/src/math/atan2.c
+++ b/src/math/atan2.c
@@ -4,11 +4,60 @@
 #include "errno.h"
 #include "nonstd/assert.h"
 
+#include "M_PI.c"
+#include "M_PI_2.c"
+
 /** arc tangent **/
 TYPE TGFN(atan2)(TYPE y, TYPE x)
 {
+	int classy = fpclassify(y);
+	int classx = fpclassify(x);
 	ASSERT_NONZERO(x);
 
+	if (classy == FP_ZERO && classx == 0) {
+		if (signbit(x)) {
+			return copysign(M_PI, y);
+		} else {
+			return y;
+		}
+	}
+
+	if (classy == FP_ZERO) {
+		if (x < 0) {
+			return copysign(M_PI, y);
+		} else if (x > 0) {
+			return y;
+		}
+	}
+
+	if (classx == FP_ZERO) {
+		if (y < 0) {
+			return - M_PI_2;
+		} else if (y > 0) {
+			return M_PI_2;
+		}
+	}
+
+	if (classx == FP_INFINITE) {
+		if (classy == FP_INFINITE) {
+			if (signbit(x)) {
+				return copysign(3.0 * M_PI / 4.0, y);
+			} else {
+				return copysign(M_PI / 4.0, y);
+			}
+		} else if (y > 0) {
+			if (signbit(x)) {
+				return copysign(M_PI, y);
+			} else {
+				return copysign(0.0, y);
+			}
+		}
+	}
+
+	if (classy == FP_INFINITE && classx != FP_INFINITE) {
+		return copysign(M_PI_2, y);
+	}
+
 	if (y == 0 && x == 0) {
 		errno = EDOM; /* ARGUMENT(y) and ARGUMENT(x) are both LITERAL(0)) */
 		return TGHUGE;
diff --git a/src/math/atanh.c b/src/math/atanh.c
index eef1cfc6..bfd4e57b 100644
--- a/src/math/atanh.c
+++ b/src/math/atanh.c
@@ -1,9 +1,24 @@
 # define TGSOURCE "atanh.c"
 #include "_tgmath.h"
 #include <math.h>
+#include "fenv.h"
 
 TYPE TGFN(atanh)(TYPE x)
 {
+	if (fpclassify(x) == FP_ZERO) {
+		return x;
+	}
+
+	if (TGFN(fabs)(x) == 1.0) {
+		feraiseexcept(FE_DIVBYZERO);
+		return copysign(INFINITY, x);
+	}
+
+	if (TGFN(fabs)(x) > 1.0) {
+		feraiseexcept(FE_INVALID);
+		return NAN;
+	}
+
 	return x;
 }
 
diff --git a/src/math/cbrt.c b/src/math/cbrt.c
index 156093ba..74c48d52 100644
--- a/src/math/cbrt.c
+++ b/src/math/cbrt.c
@@ -4,6 +4,11 @@
 
 TYPE TGFN(cbrt)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:		return x;
+	case FP_INFINITE:	return x;
+	}
+
 	return x;
 }
 
diff --git a/src/math/ceil.c b/src/math/ceil.c
index 7af395ac..d90ba15f 100644
--- a/src/math/ceil.c
+++ b/src/math/ceil.c
@@ -6,6 +6,12 @@
 /** round up to nearest integer **/
 TYPE TGFN(ceil)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:		return x;
+	case FP_INFINITE:	return x;
+	default:		break;
+	}
+
 	if (0) {
 		errno = ERANGE; /* The result cannot be represented */
 		/* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */
diff --git a/src/math/copysign.c b/src/math/copysign.c
index a462348b..2ace3c1d 100644
--- a/src/math/copysign.c
+++ b/src/math/copysign.c
@@ -4,7 +4,16 @@
 
 TYPE TGFN(copysign)(TYPE x, TYPE y)
 {
-	return x - y;
+	if (isnan(x)) {
+		if (signbit(y)) {
+			/* return -NaN; */
+		} else {
+			/* return NaN; */
+		}
+	}
+
+	x = TGFN(fabs)(x);
+	return signbit(y) ? -x : x;
 }
 
 /*
diff --git a/src/math/cos.c b/src/math/cos.c
index d10629a3..86a1c995 100644
--- a/src/math/cos.c
+++ b/src/math/cos.c
@@ -1,6 +1,7 @@
 # define TGSOURCE "cos.c"
 #include <math.h>
 #include "_tgmath.h"
+#include "fenv.h"
 #include "errno.h"
 
 /** cosine **/
@@ -12,6 +13,12 @@ TYPE TGFN(cos)(TYPE x)
 	TYPE cosine = 1.0;
 	TYPE power = 1.0;
 
+	switch (fpclassify(x)) {
+	case FP_ZERO:		return 1.0;
+	case FP_INFINITE:	feraiseexcept(FE_INVALID); return NAN;
+	default:		break;
+	}
+
 	if (0) {
 		errno = ERANGE; /* The result cannot be represented */
 		/* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */
diff --git a/src/math/cosh.c b/src/math/cosh.c
index 9065cd03..0f981e36 100644
--- a/src/math/cosh.c
+++ b/src/math/cosh.c
@@ -6,6 +6,12 @@
 /** hyperbolic cosine **/
 TYPE TGFN(cosh)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:		return 1.0;
+	case FP_INFINITE:	return TGFN(fabs)(x);
+	default:		break;
+	}
+
 	if (0) {
 		errno = ERANGE; /* The magnitude of ARGUMENT(x) is too large */
 		/* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */
diff --git a/src/math/erf.c b/src/math/erf.c
index 33059480..2115982b 100644
--- a/src/math/erf.c
+++ b/src/math/erf.c
@@ -4,6 +4,12 @@
 
 TYPE TGFN(erf)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:		return x;
+	case FP_INFINITE:	return signbit(x) ? -1.0 : 1.0;
+	default:		break;
+	}
+
 	return x;
 }
 
diff --git a/src/math/erfc.c b/src/math/erfc.c
index 84579e10..b937590d 100644
--- a/src/math/erfc.c
+++ b/src/math/erfc.c
@@ -4,6 +4,10 @@
 
 TYPE TGFN(erfc)(TYPE x)
 {
+	if (fpclassify(x) == FP_INFINITE) {
+		return signbit(x) ? 2.0 : 0.0;
+	}
+
 	return x;
 }
 
diff --git a/src/math/exp.c b/src/math/exp.c
index f244a60f..412697b2 100644
--- a/src/math/exp.c
+++ b/src/math/exp.c
@@ -12,6 +12,12 @@ TYPE TGFN(exp)(TYPE x)
 	TYPE exponent = 1.0;
 	TYPE power = 1.0;
 
+	switch (fpclassify(x)) {
+	case FP_ZERO:		return 1.0;
+	case FP_INFINITE:	return signbit(x) ? 0.0 : x;
+	default:		break;
+	}
+
 	if (0) {
 		errno = ERANGE; /* The magnitude of ARGUMENT(x) is too large */
 		/* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */
diff --git a/src/math/exp2.c b/src/math/exp2.c
index 599634e6..16ff5940 100644
--- a/src/math/exp2.c
+++ b/src/math/exp2.c
@@ -4,6 +4,12 @@
 
 TYPE TGFN(exp2)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:		return 1.0;
+	case FP_INFINITE:	return signbit(x) ? 0.0 : x;
+	default:		break;
+	}
+
 	return x;
 }
 
diff --git a/src/math/expm1.c b/src/math/expm1.c
index 152c6c1e..40684b05 100644
--- a/src/math/expm1.c
+++ b/src/math/expm1.c
@@ -4,6 +4,12 @@
 
 TYPE TGFN(expm1)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:		return x;
+	case FP_INFINITE:	return signbit(x) ? -1.0 : x;
+	default:		break;
+	}
+
 	return x;
 }
 
diff --git a/src/math/fabs.c b/src/math/fabs.c
index 9664b702..5fc42750 100644
--- a/src/math/fabs.c
+++ b/src/math/fabs.c
@@ -6,6 +6,12 @@
 /** absolute value **/
 TYPE TGFN(fabs)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:		return 0.0;
+	case FP_INFINITE:	return INFINITY;
+	default:		break;
+	}
+
 	if (0) {
 		errno = ERANGE; /* The result cannot be represented */
 		/* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */
@@ -13,7 +19,7 @@ TYPE TGFN(fabs)(TYPE x)
 	}
 
 	/* RETURN_SUCCESS(ABS(ARGUMENT(x))); */
-	return x < 0.0 ? -x : x;
+	return signbit(x) ? -x : x;
 }
 
 /***
diff --git a/src/math/fdim.c b/src/math/fdim.c
index c39c6753..4f6b3896 100644
--- a/src/math/fdim.c
+++ b/src/math/fdim.c
@@ -4,7 +4,11 @@
 
 TYPE TGFN(fdim)(TYPE x, TYPE y)
 {
-	return x - y;
+	if (x > y) {
+		return x - y;
+	}
+
+	return 0;
 }
 
 /*
diff --git a/src/math/floor.c b/src/math/floor.c
index baa7d3a6..41350fe8 100644
--- a/src/math/floor.c
+++ b/src/math/floor.c
@@ -6,6 +6,12 @@
 /** round down to nearest integer **/
 TYPE TGFN(floor)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:		return x;
+	case FP_INFINITE:	return x;
+	default:		break;
+	}
+
 	if (0) {
 		errno = ERANGE; /* The result cannot be represented */
 		/* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */
diff --git a/src/math/fma.c b/src/math/fma.c
index 18b059d1..845afea4 100644
--- a/src/math/fma.c
+++ b/src/math/fma.c
@@ -1,10 +1,32 @@
 # define TGSOURCE "fma.c"
 #include "_tgmath.h"
 #include <math.h>
+#include "fenv.h"
 
 TYPE TGFN(fma)(TYPE x, TYPE y, TYPE z)
 {
-	return x - y - z;
+	int classx = fpclassify(x);
+	int classy = fpclassify(y);
+	/* int classz = fpclassify(z); */
+
+	if (classx == FP_INFINITE && classy == FP_ZERO) {
+		feraiseexcept(FE_INVALID);
+		return NAN;
+	}
+
+	if (classx == FP_ZERO && classy == FP_INFINITE) {
+		feraiseexcept(FE_INVALID);
+		return NAN;
+	}
+
+	/*
+	if (x * y == -z && fpclassify(x*y) == classz == FP_INFINITE) {
+		feraiseexcept(FE_INVALID);
+		return NAN;
+	}
+	*/
+
+	return x * y + z;
 }
 
 /*
diff --git a/src/math/fmax.c b/src/math/fmax.c
index dc354ca3..7e5d8a61 100644
--- a/src/math/fmax.c
+++ b/src/math/fmax.c
@@ -4,7 +4,13 @@
 
 TYPE TGFN(fmax)(TYPE x, TYPE y)
 {
-	return x - y;
+	if (isnan(x) && !isnan(y)) {
+		return y;
+	} else if (!isnan(x) && isnan(y)) {
+		return x;
+	}
+
+	return x > y ? x : y;
 }
 
 /*
diff --git a/src/math/fmin.c b/src/math/fmin.c
index 822f893b..4408f432 100644
--- a/src/math/fmin.c
+++ b/src/math/fmin.c
@@ -4,7 +4,13 @@
 
 TYPE TGFN(fmin)(TYPE x, TYPE y)
 {
-	return x - y;
+	if (isnan(x) && !isnan(y)) {
+		return y;
+	} else if (!isnan(x) && isnan(y)) {
+		return x;
+	}
+
+	return x < y ? x : y;
 }
 
 /*
diff --git a/src/math/fmod.c b/src/math/fmod.c
index d8254d40..4370830f 100644
--- a/src/math/fmod.c
+++ b/src/math/fmod.c
@@ -2,10 +2,27 @@
 #include <math.h>
 #include "_tgmath.h"
 #include "errno.h"
+#include "fenv.h"
 
 /** floating-point remainder **/
 TYPE TGFN(fmod)(TYPE x, TYPE y)
 {
+	int classx = fpclassify(x);
+	int classy = fpclassify(y);
+
+	if (classx == FP_ZERO && classy != FP_ZERO) {
+		return x;
+	}
+
+	if (classx == FP_INFINITE && classy == FP_ZERO) {
+		feraiseexcept(FE_INVALID);
+		return NAN;
+	}
+
+	if (classx != FP_INFINITE && classy == FP_INFINITE) {
+		return x;
+	}
+
 	if (y == 0) {
 		return 0;
 	}
diff --git a/src/math/frexp.c b/src/math/frexp.c
index 40bb59c8..b60dadc5 100644
--- a/src/math/frexp.c
+++ b/src/math/frexp.c
@@ -6,6 +6,13 @@
 /** extract mantissa and exponent **/
 TYPE TGFN(frexp)(TYPE value, int *exp)
 {
+	switch (fpclassify(value)) {
+	case FP_ZERO:		*exp = 0; return value;
+	case FP_INFINITE:	*exp = /* unspecified */0; return value;
+	case FP_NAN:		*exp = /* unspecified */0; return value;
+	default:		break;
+	}
+
 	if (0) {
 		errno = ERANGE; /* The result cannot be represented */
 		/* RETURN(CONSTANT(HUGE_VAL), A range error occurred); */
diff --git a/src/math/hypot.c b/src/math/hypot.c
index 6ab2ecf7..ce902c59 100644
--- a/src/math/hypot.c
+++ b/src/math/hypot.c
@@ -4,6 +4,14 @@
 
 TYPE TGFN(hypot)(TYPE x, TYPE y)
 {
+	if (fpclassify(y) == FP_ZERO) {
+		return fabs(x);
+	}
+
+	if (fpclassify(x) == FP_INFINITE) {
+		return INFINITY;
+	}
+
 	return x - y;
 }
 
diff --git a/src/math/ilogb.c b/src/math/ilogb.c
index 0d826e04..d1295815 100644
--- a/src/math/ilogb.c
+++ b/src/math/ilogb.c
@@ -1,9 +1,18 @@
 # define TGSOURCE "ilogb.c"
 #include "_tgmath.h"
 #include <math.h>
+#include "limits.h"
+#include "fenv.h"
 
-TYPE TGFN(ilogb)(TYPE x)
+int TGFN(ilogb)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:		feraiseexcept(FE_INVALID); return FP_ILOGB0;
+	case FP_INFINITE:	feraiseexcept(FE_INVALID); return INT_MAX;
+	case FP_NAN:		feraiseexcept(FE_INVALID); return FP_ILOGBNAN;
+	default:		break;
+	}
+
 	return x;
 }
 
diff --git a/src/math/ldexp.c b/src/math/ldexp.c
index ab143982..33e2f381 100644
--- a/src/math/ldexp.c
+++ b/src/math/ldexp.c
@@ -6,6 +6,16 @@
 /** multiply by a power of 2 **/
 TYPE TGFN(ldexp)(TYPE x, int exp)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:           return x;
+	case FP_INFINITE:       return x;
+	default:                break;
+	}
+
+	if (exp == 0) {
+		return x;
+	}
+
 	if (0) {
 		errno = ERANGE; /* The result cannot be represented */
 		/* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */
diff --git a/src/math/lgamma.c b/src/math/lgamma.c
index 3612c825..3f0799c2 100644
--- a/src/math/lgamma.c
+++ b/src/math/lgamma.c
@@ -1,9 +1,25 @@
 # define TGSOURCE "lgamma.c"
 #include "_tgmath.h"
 #include <math.h>
+#include "fenv.h"
 
 TYPE TGFN(lgamma)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:		feraiseexcept(FE_DIVBYZERO); return INFINITY;
+	case FP_INFINITE:	return INFINITY;
+	default:		break;
+	}
+
+	if (x == 1.0 || x == 2.0) {
+		return 0.0;
+	}
+
+	if (signbit(x)) {
+		feraiseexcept(FE_DIVBYZERO);
+		return INFINITY;
+	}
+
 	return x;
 }
 
diff --git a/src/math/llrint.c b/src/math/llrint.c
index a967f651..299a8258 100644
--- a/src/math/llrint.c
+++ b/src/math/llrint.c
@@ -1,9 +1,23 @@
 # define TGSOURCE "llrint.c"
 #include "_tgmath.h"
 #include <math.h>
+#include "limits.h"
+#include "fenv.h"
 
 long long int TGFN(llrint)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:
+		return 0;
+
+	case FP_INFINITE:
+		feraiseexcept(FE_INVALID);
+		return signbit(x) ? LLONG_MIN : LLONG_MAX;
+
+	default:
+		break;
+	}
+
 	return x;
 }
 
diff --git a/src/math/log.c b/src/math/log.c
index ee809cf5..abc29e4e 100644
--- a/src/math/log.c
+++ b/src/math/log.c
@@ -2,11 +2,24 @@
 #include <math.h>
 #include "_tgmath.h"
 #include "errno.h"
+#include "fenv.h"
 
 /** natural logarithm **/
 TYPE TGFN(log)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:		feraiseexcept(FE_DIVBYZERO); return - INFINITY;
+	case FP_INFINITE:	if (!signbit(x)) { return x; }
+	default:		break;
+	}
+
+	if (x == 1.0) {
+		return 0.0;
+	}
+
 	if (x < 0) {
+		feraiseexcept(FE_INVALID);
+		return NAN;
 		errno = EDOM; /* ARGUMENT(x) is negative */
 		return TGHUGE;
 	}
diff --git a/src/math/log10.c b/src/math/log10.c
index cea1d45b..e504df44 100644
--- a/src/math/log10.c
+++ b/src/math/log10.c
@@ -2,11 +2,24 @@
 #include <math.h>
 #include "_tgmath.h"
 #include "errno.h"
+#include "fenv.h"
 
 /** base-10 logarithm **/
 TYPE TGFN(log10)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:		feraiseexcept(FE_INVALID); return - INFINITY;
+	case FP_INFINITE:	if (!signbit(x)) { return x; } break;
+	default:		break;
+	}
+
+	if (x == 1.0) {
+		return 0.0;
+	}
+
 	if (x < 0) {
+		feraiseexcept(FE_INVALID);
+		return NAN;
 		errno = EDOM; /* ARGUMENT(x) is negative */
 		return TGHUGE;
 	}
diff --git a/src/math/log1p.c b/src/math/log1p.c
index 30e4a8d1..5f3e2684 100644
--- a/src/math/log1p.c
+++ b/src/math/log1p.c
@@ -1,9 +1,26 @@
 # define TGSOURCE "log1p.c"
 #include "_tgmath.h"
 #include <math.h>
+#include "fenv.h"
 
 TYPE TGFN(log1p)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:		return x;
+	case FP_INFINITE:	if (!signbit(x)) { return x; } break;
+	default:		break;
+	}
+
+	if (x == -1) {
+		feraiseexcept(FE_DIVBYZERO);
+		return - INFINITY;
+	}
+
+	if (x < -1) {
+		feraiseexcept(FE_INVALID);
+		return NAN;
+	}
+
 	return x;
 }
 
diff --git a/src/math/log2.c b/src/math/log2.c
index c438d7ce..55569065 100644
--- a/src/math/log2.c
+++ b/src/math/log2.c
@@ -1,9 +1,25 @@
 #define TGSOURCE "log2.c"
 #include "_tgmath.h"
 #include <math.h>
+#include "fenv.h"
 
 TYPE TGFN(log2)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:		feraiseexcept(FE_DIVBYZERO); return - INFINITY;
+	case FP_INFINITE:	if (!signbit(x)) { return x; } break;
+	default:		break;
+	}
+
+	if (x == 1.0) {
+		return 0.0;
+	}
+
+	if (x < 0.0) {
+		feraiseexcept(FE_INVALID);
+		return NAN;
+	}
+
 	return x;
 }
 
diff --git a/src/math/logb.c b/src/math/logb.c
index d0bbfaca..8dd29371 100644
--- a/src/math/logb.c
+++ b/src/math/logb.c
@@ -1,9 +1,16 @@
 # define TGSOURCE "logb.c"
 #include "_tgmath.h"
 #include <math.h>
+#include "fenv.h"
 
 TYPE TGFN(logb)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:		feraiseexcept(FE_DIVBYZERO); return - INFINITY;
+	case FP_INFINITE:	return x;
+	default:		break;
+	}
+
 	return x;
 }
 
diff --git a/src/math/lrint.c b/src/math/lrint.c
index 16eaf7be..3ce1dd66 100644
--- a/src/math/lrint.c
+++ b/src/math/lrint.c
@@ -1,9 +1,23 @@
 # define TGSOURCE "lrint.c"
 #include "_tgmath.h"
 #include <math.h>
+#include "limits.h"
+#include "fenv.h"
 
 long int TGFN(lrint)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:
+		return 0;
+
+	case FP_INFINITE:
+		feraiseexcept(FE_INVALID);
+		return signbit(x) ? LONG_MIN : LONG_MAX;
+
+	default:
+		break;
+	}
+
 	return x;
 }
 
diff --git a/src/math/modf.c b/src/math/modf.c
index fc09fc2d..63db6094 100644
--- a/src/math/modf.c
+++ b/src/math/modf.c
@@ -6,6 +6,12 @@
 /** decompose floating-point numbers **/
 TYPE TGFN(modf)(TYPE value, TYPE *iptr)
 {
+	switch (fpclassify(value)) {
+	case FP_INFINITE:	*iptr = value; return copysign(0.0, value);
+	case FP_NAN:		*iptr = value; return value;
+	default:		break;
+	}
+
 	if (0) {
 		errno = ERANGE; /* The result cannot be represented */
 		/* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */
@@ -13,7 +19,6 @@ TYPE TGFN(modf)(TYPE value, TYPE *iptr)
 	}
 
 	/* RETURN_SUCCESS(the signed fractional part of ARGUMENT(value)); */
-	(void)iptr;
 	return value;
 }
 
diff --git a/src/math/nan.c b/src/math/nan.c
index bd0cbbc4..c0874e5f 100644
--- a/src/math/nan.c
+++ b/src/math/nan.c
@@ -1,11 +1,22 @@
 # define TGSOURCE "nan.c"
 #include "_tgmath.h"
 #include <math.h>
+#include "string.h"
+#include "stdlib.h"
+
+#define strtodf strtof
+#define strtodl strtold
 
 TYPE TGFN(nan)(const char *tagp)
 {
-	(void)tagp;
-	return 0.0;
+	if (tagp) {
+		char ncharseq[strlen(tagp) + 6];
+		strcpy(ncharseq, "NAN(");
+		strcat(ncharseq, tagp);
+		strcat(ncharseq, ")");
+		return TGFN(strtod)(ncharseq, (char**)NULL);
+	}
+	return TGFN(strtod)("NAN", (char**)NULL);
 }
 
 /*
diff --git a/src/math/nearbyint.c b/src/math/nearbyint.c
index d5d2cc21..ac9417db 100644
--- a/src/math/nearbyint.c
+++ b/src/math/nearbyint.c
@@ -4,6 +4,12 @@
 
 TYPE TGFN(nearbyint)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:		return x;
+	case FP_INFINITE:	return x;
+	default:		break;
+	}
+
 	return x;
 }
 
diff --git a/src/math/pow.c b/src/math/pow.c
index c0d613e5..8a070005 100644
--- a/src/math/pow.c
+++ b/src/math/pow.c
@@ -2,10 +2,91 @@
 #include <math.h>
 #include "_tgmath.h"
 #include "errno.h"
+#include "fenv.h"
 
 /** exponentiation **/
 TYPE TGFN(pow)(TYPE x, TYPE y)
 {
+	int classx = fpclassify(x);
+	int classy = fpclassify(y);
+
+	if (classx == FP_ZERO) {
+		if (classy == FP_INFINITE && signbit(y)) {
+			feraiseexcept(FE_DIVBYZERO);
+			return INFINITY;
+		}
+
+		if (y < 0) {
+			feraiseexcept(FE_DIVBYZERO);
+			/* if (y == odd integer) { */
+				return signbit(x) ? - INFINITY : INFINITY;
+			/* } else { */
+				return INFINITY;
+			/* } */
+		}
+
+		if (y > 0) {
+			/* if (y == odd integer) { */
+				return x;
+			/* } else { */
+				return 0.0;
+			/* } */
+		}
+	}
+
+	if (x == -1.0 && classy == FP_INFINITE) {
+		return 1;
+	}
+
+	if (x == 1.0) {
+		return 1.0;
+	}
+
+	if (classy == FP_ZERO) {
+		return 1.0;
+	}
+
+	if (x < 0 /* && y != integer */) {
+		feraiseexcept(FE_INVALID);
+		return NAN;
+	}
+
+	if (classy == FP_INFINITE) {
+		int sb = signbit(y);
+		TYPE ab = TGFN(fabs)(x);
+		if ((sb && ab < 1.0) || (!sb && ab > 1.0)) {
+			return INFINITY;
+		} else {
+			return 0.0;
+		}
+	}
+
+	if (classx == FP_INFINITE) {
+		if (!signbit(x)) {
+			if (y < 0) {
+				return 0.0;
+			} else if (y > 0) {
+				return INFINITY;
+			}
+		}
+
+		if (y < 0) {
+			/* if (y == odd integer ) { */
+				return -0.0;
+			/* } else { */
+				return 0.0;
+			/* } */
+		}
+
+		if (y > 0) {
+			/* if (y == odd integer ) { */
+				return - INFINITY;
+			/* } else { */
+				return INFINITY;
+			/* } */
+		}
+	}
+
 	if (x < 0 /* && !isintegral(y) */) {
 		errno = EDOM; /* ARGUMENT(x) is negative and ARGUMENT(y) is not an integer */
 		return TGHUGE;
diff --git a/src/math/rint.c b/src/math/rint.c
index a1868cbb..a97c4709 100644
--- a/src/math/rint.c
+++ b/src/math/rint.c
@@ -4,6 +4,12 @@
 
 TYPE TGFN(rint)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:		return x;
+	case FP_INFINITE:	return x;
+	default:		break;
+	}
+
 	return x;
 }
 
diff --git a/src/math/round.c b/src/math/round.c
index c154bdbf..5eedd578 100644
--- a/src/math/round.c
+++ b/src/math/round.c
@@ -1,10 +1,33 @@
 # define TGSOURCE "round.c"
 #include "_tgmath.h"
 #include <math.h>
+#include "fenv.h"
+
+#ifndef __GNUC__
+#pragma STDC FENV_ACCESS ON
+#endif
 
 TYPE TGFN(round)(TYPE x)
 {
-	return x;
+	switch (fpclassify(x)) {
+	case FP_ZERO:		return x;
+	case FP_INFINITE:	return x;
+	default:		break;
+	}
+
+	fenv_t save_env;
+	feholdexcept(&save_env);
+
+	TYPE ret = TGFN(rint)(x);
+
+	if (fetestexcept(FE_INEXACT)) {
+		fesetround(FE_TOWARDZERO);
+		ret = TGFN(rint)(TGFN(copysign)(0.5 + TGFN(fabs)(x), x));
+	}
+
+	feupdateenv(&save_env);
+
+	return ret;
 }
 
 /*
diff --git a/src/math/scalbn.c b/src/math/scalbn.c
index bf5b4aae..f01843ce 100644
--- a/src/math/scalbn.c
+++ b/src/math/scalbn.c
@@ -4,7 +4,17 @@
 
 TYPE TGFN(scalbn)(TYPE x, int n)
 {
-	return x - n;
+	switch (fpclassify(x)) {
+	case FP_ZERO:		return x;
+	case FP_INFINITE:	return x;
+	default:		break;
+	}
+
+	if (n == 0) {
+		return x;
+	}
+
+	return TGFN(ldexp)(x, n);
 }
 
 /*
diff --git a/src/math/sin.c b/src/math/sin.c
index b753cecb..e93a663b 100644
--- a/src/math/sin.c
+++ b/src/math/sin.c
@@ -2,6 +2,7 @@
 #include <math.h>
 #include "_tgmath.h"
 #include "errno.h"
+#include "fenv.h"
 
 /** sine **/
 TYPE TGFN(sin)(TYPE x)
@@ -12,6 +13,12 @@ TYPE TGFN(sin)(TYPE x)
 	TYPE sine = 0.0;
 	TYPE power = 1.0;
 
+	switch (fpclassify(x)) {
+	case FP_ZERO:		return x;
+	case FP_INFINITE:	feraiseexcept(FE_INVALID); return NAN;
+	default:		break;
+	}
+
 	if (0) {
 		errno = ERANGE; /* The result cannot be represented */
 		/* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */
diff --git a/src/math/sinh.c b/src/math/sinh.c
index e5942de6..2e8d498c 100644
--- a/src/math/sinh.c
+++ b/src/math/sinh.c
@@ -6,6 +6,12 @@
 /** hyperbolic sine **/
 TYPE TGFN(sinh)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:		return x;
+	case FP_INFINITE:	return x;
+	default:		break;
+	}
+
 	if (0) {
 		errno = ERANGE; /* The magnitude of ARGUMENT(x) is too large */
 		/* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */
diff --git a/src/math/tan.c b/src/math/tan.c
index 2e60bd2f..061e6cdd 100644
--- a/src/math/tan.c
+++ b/src/math/tan.c
@@ -2,10 +2,17 @@
 #include <math.h>
 #include "_tgmath.h"
 #include "errno.h"
+#include "fenv.h"
 
 /** tangent **/
 TYPE TGFN(tan)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:		return x;
+	case FP_INFINITE:	feraiseexcept(FE_INVALID); return NAN;
+	default:		break;
+	}
+
 	if (0) {
 		errno = ERANGE; /* The result cannot be represented */
 		/* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */
diff --git a/src/math/tanh.c b/src/math/tanh.c
index d64be283..cf5bf028 100644
--- a/src/math/tanh.c
+++ b/src/math/tanh.c
@@ -6,6 +6,12 @@
 /** hyperbolic tangent **/
 TYPE TGFN(tanh)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:		return x;
+	case FP_INFINITE:	return signbit(x) ? -1.0 : 1.0;
+	default:		break;
+	}
+
 	if (0) {
 		errno = ERANGE; /* The result cannot be represented */
 		/* RETURN_FAILURE(CONSTANT(HUGE_VAL), A range error occurred); */
diff --git a/src/math/tgamma.c b/src/math/tgamma.c
index b80584ac..63e25a5b 100644
--- a/src/math/tgamma.c
+++ b/src/math/tgamma.c
@@ -1,9 +1,27 @@
 # define TGSOURCE "tgamma.c"
 #include "_tgmath.h"
 #include <math.h>
+#include "fenv.h"
 
 TYPE TGFN(tgamma)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:
+		feraiseexcept(FE_DIVBYZERO);
+		return signbit(x) ? - INFINITY : INFINITY;
+
+	case FP_INFINITE:
+		if (signbit(x)) {
+			feraiseexcept(FE_INVALID);
+			return NAN;
+		} else {
+			return INFINITY;
+		}
+
+	default:
+		break;
+	}
+
 	return x;
 }
 
diff --git a/src/math/trunc.c b/src/math/trunc.c
index 4b1081a6..4ee27109 100644
--- a/src/math/trunc.c
+++ b/src/math/trunc.c
@@ -4,6 +4,12 @@
 
 TYPE TGFN(trunc)(TYPE x)
 {
+	switch (fpclassify(x)) {
+	case FP_ZERO:		return x;
+	case FP_INFINITE:	return x;
+	default:		break;
+	}
+
 	return x;
 }