diff --git a/include/tgmath.h b/include/tgmath.h
index 832b052b5e106c4b2de2cd20e11dbcf53625c6f5..25ad0071d4e8cfe592e7891f02a7f77f7a8152e4 100644
--- a/include/tgmath.h
+++ b/include/tgmath.h
@@ -27,34 +27,67 @@ sizeof(double) == sizeof(long double)
 /* return type */
 
 #ifdef __GNUC__
+/*
+the conditional expression (that selects the right function
+for evaluation) should be casted to the return type of the
+selected function (otherwise the result type is determined
+by the conversion rules applied to all the function return
+types so it is long double or long double _Complex)
+
+this cannot be done in c99 (c11 has _Generic that would help
+but is not yet widely supported), so the typeof extension of
+gcc is used and that ?: has special semantics when one of
+the operands is a null pointer constant
+
+unfortunately the standard is overly strict about the
+definition of null pointer constants (current gcc does not
+follow the standard but clang does) so we have to use a hack
+to be able to tell integer and floating-point types apart:
+both gcc and clang evaluate sizeof(void) to 1 even though it
+is a constraint violation, we only use this on clang for now
+*/
+/* predicates that evaluate to integer constant expressions */
+#ifdef __clang__
+/* TODO: __c_IS_FP(1.0) is a constraint violation */
+#define __c_IS_FP(x) (!(sizeof*(__typeof__(0?(int*)0:(void*)__IS_FP(x)))0-1))
+#else
+#define __c_IS_FP(x) __IS_FP(x)
+#endif
+#define __c_IS_CX(x) (__c_IS_FP(x) && sizeof(x) == sizeof((x)+I))
+#define __c_FLTCX(x) (__c_IS_CX(x) && sizeof(x) == sizeof(float complex))
+#define __c_DBLCX(x) (__c_IS_CX(x) && sizeof(x) == sizeof(double complex))
+#define __c_LDBLCX(x) (__c_IS_CX(x) && sizeof(x) == sizeof(long double complex) && sizeof(long double) != sizeof(double))
+/* if c then t else void */
+#define __type1(c,t) __typeof__(*(0?(t*)0:(void*)!(c)))
+/* if c then t1 else t2 */
+#define __type2(c,t1,t2) __typeof__(*(0?(__type1(c,t1)*)0:(__type1(!(c),t2)*)0))
 /* cast to double when x is integral, otherwise use typeof(x) */
-#define __RETCAST(x) (__typeof__(*( \
-	0 ? (__typeof__(0 ? (double *)0 : (void *)__IS_FP(x)))0 : \
-	    (__typeof__(0 ? (__typeof__(x) *)0 : (void *)!__IS_FP(x)))0 )))
-/* 2 args case, consider complex types (for cpow) */
-#define __RETCAST_2(x, y) (__typeof__(*( \
-	0 ? (__typeof__(0 ? (double *)0 : \
-		(void *)!((!__IS_FP(x) || !__IS_FP(y)) && __FLT((x)+(y)+1.0f))))0 : \
-	0 ? (__typeof__(0 ? (double complex *)0 : \
-		(void *)!((!__IS_FP(x) || !__IS_FP(y)) && __FLTCX((x)+(y)))))0 : \
-	    (__typeof__(0 ? (__typeof__((x)+(y)) *)0 : \
-		(void *)((!__IS_FP(x) || !__IS_FP(y)) && (__FLT((x)+(y)+1.0f) || __FLTCX((x)+(y))))))0 )))
-/* 3 args case, don't consider complex types (fma only) */
-#define __RETCAST_3(x, y, z) (__typeof__(*( \
-	0 ? (__typeof__(0 ? (double *)0 : \
-		(void *)!((!__IS_FP(x) || !__IS_FP(y) || !__IS_FP(z)) && __FLT((x)+(y)+(z)+1.0f))))0 : \
-	    (__typeof__(0 ? (__typeof__((x)+(y)) *)0 : \
-		(void *)((!__IS_FP(x) || !__IS_FP(y) || !__IS_FP(z)) && __FLT((x)+(y)+(z)+1.0f))))0 )))
+#define __RETCAST(x) ( \
+	__type2(__c_IS_FP(x), __typeof__(x), double))
+/* 2 args case, should work for complex types (cpow) */
+#define __RETCAST_2(x, y) ( \
+	__type2(__c_IS_FP(x) && __c_IS_FP(y), \
+		__typeof__((x)+(y)), \
+		__typeof__((x)+(y)+1.0)))
+/* 3 args case (fma only) */
+#define __RETCAST_3(x, y, z) ( \
+	__type2(__c_IS_FP(x) && __c_IS_FP(y) && __c_IS_FP(z), \
+		__typeof__((x)+(y)+(z)), \
+		__typeof__((x)+(y)+(z)+1.0)))
 /* drop complex from the type of x */
-#define __TO_REAL(x) *( \
-	0 ? (__typeof__(0 ? (double *)0 : (void *)!__DBLCX(x)))0 : \
-	0 ? (__typeof__(0 ? (float *)0 : (void *)!__FLTCX(x)))0 : \
-	0 ? (__typeof__(0 ? (long double *)0 : (void *)!__LDBLCX(x)))0 : \
-	    (__typeof__(0 ? (__typeof__(x) *)0 : (void *)__IS_CX(x)))0 )
+/* TODO: wrong when sizeof(long double)==sizeof(double) */
+#define __RETCAST_REAL(x) (  \
+	__type2(sizeof(__RETCAST(x)0+I)==sizeof(float complex), float, \
+	__type2(sizeof(__RETCAST(x)0+I)==sizeof(double complex), double, \
+		long double)))
+/* add complex to the type of x */
+#define __RETCAST_CX(x) (__typeof__(__RETCAST(x)0+I))
 #else
 #define __RETCAST(x)
 #define __RETCAST_2(x, y)
 #define __RETCAST_3(x, y, z)
+#define __RETCAST_REAL(x)
+#define __RETCAST_CX(x)
 #endif
 
 /* function selection */
@@ -76,12 +109,12 @@ sizeof(double) == sizeof(long double)
 	__LDBL((x)+(y)) ? fun ## l (x, y) : \
 	fun(x, y) ))
 
-#define __tg_complex(fun, x) (__RETCAST((x)+I)( \
+#define __tg_complex(fun, x) (__RETCAST_CX(x)( \
 	__FLTCX((x)+I) && __IS_FP(x) ? fun ## f (x) : \
 	__LDBLCX((x)+I) ? fun ## l (x) : \
 	fun(x) ))
 
-#define __tg_complex_retreal(fun, x) (__RETCAST(__TO_REAL(x))( \
+#define __tg_complex_retreal(fun, x) (__RETCAST_REAL(x)( \
 	__FLTCX((x)+I) && __IS_FP(x) ? fun ## f (x) : \
 	__LDBLCX((x)+I) ? fun ## l (x) : \
 	fun(x) ))
@@ -115,7 +148,7 @@ sizeof(double) == sizeof(long double)
 	__LDBL((x)+(y)) ? powl(x, y) : \
 	pow(x, y) ))
 
-#define __tg_real_complex_fabs(x) (__RETCAST(__TO_REAL(x))( \
+#define __tg_real_complex_fabs(x) (__RETCAST_REAL(x)( \
 	__FLTCX(x) ? cabsf(x) : \
 	__DBLCX(x) ? cabs(x) : \
 	__LDBLCX(x) ? cabsl(x) : \