Browse Source

math: add exception handling functionality

According to standards SVID and SYSV.
Modified lgamma calling in case when 'signgam' variable should not be used.

Signed-off-by: Waldemar Brodkorb <wbx@uclibc-ng.org>
Sergey Cherkashin 4 years ago
parent
commit
ea38f4d89c
100 changed files with 4375 additions and 433 deletions
  1. 25 1
      include/math.h
  2. 5 2
      libc/sysdeps/linux/common/bits/mathcalls.h
  3. 3 7
      libc/sysdeps/linux/powerpc/bits/mathdef.h
  4. 0 11
      libc/sysdeps/linux/powerpc/bits/wordsize.h
  5. 34 3
      libm/Makefile.in
  6. 1 4
      libm/e_acos.c
  7. 1 4
      libm/e_acosh.c
  8. 0 3
      libm/e_asin.c
  9. 1 4
      libm/e_atan2.c
  10. 0 3
      libm/e_atanh.c
  11. 1 4
      libm/e_cosh.c
  12. 0 3
      libm/e_exp.c
  13. 0 2
      libm/e_exp10.c
  14. 0 3
      libm/e_fmod.c
  15. 1 4
      libm/e_hypot.c
  16. 0 4
      libm/e_j0.c
  17. 0 4
      libm/e_j1.c
  18. 0 4
      libm/e_jn.c
  19. 2 28
      libm/e_lgamma_r.c
  20. 0 3
      libm/e_log.c
  21. 0 3
      libm/e_log10.c
  22. 0 2
      libm/e_log2.c
  23. 0 3
      libm/e_pow.c
  24. 0 4
      libm/e_remainder.c
  25. 0 5
      libm/e_scalb.c
  26. 0 3
      libm/e_sinh.c
  27. 0 3
      libm/e_sqrt.c
  28. 32 157
      libm/float_wrappers.c
  29. 973 0
      libm/k_standard.c
  30. 31 0
      libm/k_standardf.c
  31. 107 0
      libm/k_standardl.c
  32. 32 139
      libm/ldouble_wrappers.c
  33. 10 0
      libm/math_private.h
  34. 40 0
      libm/w_acos.c
  35. 38 0
      libm/w_acosf.c
  36. 34 0
      libm/w_acosh.c
  37. 33 0
      libm/w_acoshf.c
  38. 34 0
      libm/w_acoshl.c
  39. 42 0
      libm/w_acosl.c
  40. 39 0
      libm/w_asin.c
  41. 38 0
      libm/w_asinf.c
  42. 41 0
      libm/w_asinl.c
  43. 46 0
      libm/w_atan2.c
  44. 46 0
      libm/w_atan2f.c
  45. 44 0
      libm/w_atan2l.c
  46. 37 0
      libm/w_atanh.c
  47. 36 0
      libm/w_atanhf.c
  48. 38 0
      libm/w_atanhl.c
  49. 34 0
      libm/w_cosh.c
  50. 37 0
      libm/w_coshf.c
  51. 40 0
      libm/w_coshl.c
  52. 38 0
      libm/w_exp.c
  53. 43 0
      libm/w_exp10.c
  54. 42 0
      libm/w_exp10f.c
  55. 44 0
      libm/w_exp10l.c
  56. 32 13
      libm/w_exp2.c
  57. 22 0
      libm/w_exp2f.c
  58. 24 0
      libm/w_exp2l.c
  59. 38 0
      libm/w_expf.c
  60. 39 0
      libm/w_expl.c
  61. 34 0
      libm/w_fmod.c
  62. 33 0
      libm/w_fmodf.c
  63. 35 0
      libm/w_fmodl.c
  64. 34 0
      libm/w_hypot.c
  65. 37 0
      libm/w_hypotf.c
  66. 40 0
      libm/w_hypotl.c
  67. 63 0
      libm/w_j0.c
  68. 68 0
      libm/w_j0f.c
  69. 68 0
      libm/w_j0l.c
  70. 64 0
      libm/w_j1.c
  71. 68 0
      libm/w_j1f.c
  72. 68 0
      libm/w_j1l.c
  73. 63 0
      libm/w_jn.c
  74. 66 0
      libm/w_jnf.c
  75. 95 0
      libm/w_jnl.c
  76. 37 0
      libm/w_lgamma.c
  77. 38 0
      libm/w_lgamma_r.c
  78. 30 0
      libm/w_lgammaf.c
  79. 41 0
      libm/w_lgammaf_r.c
  80. 39 0
      libm/w_lgammal.c
  81. 44 0
      libm/w_lgammal_r.c
  82. 46 0
      libm/w_log.c
  83. 46 0
      libm/w_log10.c
  84. 45 0
      libm/w_log10f.c
  85. 45 0
      libm/w_log10l.c
  86. 46 0
      libm/w_log2.c
  87. 45 0
      libm/w_log2f.c
  88. 47 0
      libm/w_log2l.c
  89. 45 0
      libm/w_logf.c
  90. 47 0
      libm/w_logl.c
  91. 78 0
      libm/w_pow.c
  92. 76 0
      libm/w_powf.c
  93. 78 0
      libm/w_powl.c
  94. 35 0
      libm/w_remainder.c
  95. 34 0
      libm/w_remainderf.c
  96. 36 0
      libm/w_remainderl.c
  97. 86 0
      libm/w_scalb.c
  98. 80 0
      libm/w_scalbf.c
  99. 83 0
      libm/w_scalbl.c
  100. 34 0
      libm/w_sinh.c

+ 25 - 1
include/math.h

@@ -192,6 +192,30 @@ extern long double __REDIRECT_NTH (nexttowardl, (long double __x, long double __
 #if defined __USE_MISC || defined __USE_XOPEN
 /* This variable is used by `gamma' and `lgamma'.  */
 extern int signgam;
+#else
+/* This is used when standart of libm
+ * should prevent lgamma(x) of modifying signgam variable.
+ */
+# define lgammaf(arg) \
+({ \
+int local_signgam = 0; \
+float result = lgammaf_r((float)arg, &local_signgam); \
+result; \
+}) 
+
+# define lgamma(arg) \
+({ \
+int local_signgam = 0; \
+double result = lgamma_r(arg, &local_signgam); \
+result; \
+}) 
+
+# define lgammal(arg) \
+({ \
+int local_signgam = 0; \
+long double result = lgammal_r(arg, &local_signgam); \
+result; \
+})
 #endif
 
 
@@ -361,7 +385,7 @@ struct exception
 # ifdef __cplusplus
 extern int matherr (struct __exception *__exc) throw ();
 # else
-extern int matherr (struct exception *__exc);
+extern int __attribute__ ((weak)) matherr (struct exception *__exc);
 # endif
 
 # define X_TLOSS	1.41484755040568800000e+16

+ 5 - 2
libc/sysdeps/linux/common/bits/mathcalls.h

@@ -300,14 +300,17 @@ __END_NAMESPACE_C99
 __MATHCALLI (gamma,, (_Mdouble_))
 #endif
 
-#ifdef __USE_MISC
+//#ifdef __USE_MISC
+/* To provide compatibility with finite-math-only in C99 and C11 
+ * standerts function lgamma_r should be declared, when defined __USE_MISC.
+ */
 /* Reentrant version of lgamma.  This function uses the global variable
    `signgam'.  The reentrant version instead takes a pointer and stores
    the value through it.  */
 __MATHCALL (lgamma,_r, (_Mdouble_, int *__signgamp))
 /* __MATHCALLI does not work here, probably due to ,_r, */
 libm_hidden_proto(lgamma_r)
-#endif
+//#endif
 
 
 #if defined __USE_MISC || defined __USE_XOPEN_EXTENDED || defined __USE_ISOC99

+ 3 - 7
libc/sysdeps/linux/powerpc/bits/mathdef.h

@@ -44,10 +44,6 @@ typedef double double_t;	/* `double' expressions are evaluated as
 
 #endif	/* ISO C99 */
 
-#ifndef __NO_LONG_DOUBLE_MATH
-/* Signal that we do not really have a `long double'.  The disables the
-   declaration of all the `long double' function variants.  */
-# if !defined __UCLIBC_HAS_LONG_DOUBLE_MATH__
-#  define __NO_LONG_DOUBLE_MATH	1
-# endif
-#endif  /* __NO_LONG_DOUBLE_MATH */
+#if !defined __NO_LONG_DOUBLE_MATH && !defined __UCLIBC_HAS_LONG_DOUBLE_MATH__
+# define __NO_LONG_DOUBLE_MATH	1
+#endif

+ 0 - 11
libc/sysdeps/linux/powerpc/bits/wordsize.h

@@ -1,14 +1,3 @@
 /* Determine the wordsize from the preprocessor defines.  */
 
 #define __WORDSIZE	32
-
-#if !defined __NO_LONG_DOUBLE_MATH && !defined __LONG_DOUBLE_MATH_OPTIONAL
-
-/* Signal the glibc ABI didn't used to have a `long double'.
-   The changes all the `long double' function variants to be redirects
-   to the double functions.  */
-# define __LONG_DOUBLE_MATH_OPTIONAL   1
-# ifndef __LONG_DOUBLE_128__
-#  define __NO_LONG_DOUBLE_MATH        1
-# endif
-#endif

+ 34 - 3
libm/Makefile.in

@@ -64,8 +64,37 @@ libm_CSRC := \
 	s_fpclassify.c s_fpclassifyf.c s_signbit.c s_signbitf.c \
 	s_isnan.c s_isnanf.c s_isinf.c s_isinff.c s_finitef.c \
 	s_fdim.c s_fma.c s_fmax.c s_fmin.c \
-	s_remquo.c w_exp2.c \
-	cexp.c sincos.c
+	s_remquo.c \
+	cexp.c sincos.c	\
+	w_acos.c w_acosf.c w_acosl.c \
+	w_asin.c w_asinf.c w_asinl.c \
+	w_atan2.c w_atan2f.c w_atan2l.c \
+	w_hypot.c w_hypotf.c w_hypotl.c \
+	w_cosh.c w_coshf.c w_coshl.c \
+	w_exp.c w_expf.c w_expl.c \
+	w_exp2.c w_exp2f.c w_exp2l.c \
+	w_exp10.c w_exp10f.c w_exp10l.c \
+	w_j0.c w_j0f.c w_j0l.c \
+	w_j1.c w_j1f.c w_j1l.c \
+	w_jn.c w_jnf.c w_jnl.c \
+	w_lgamma_r.c w_lgammaf_r.c w_lgammal_r.c \
+	w_lgamma.c w_lgammaf.c w_lgammal.c \
+	w_tgamma.c w_tgammaf.c w_tgammal.c \
+	w_log.c w_logf.c w_logl.c \
+	w_log2.c w_log2f.c w_log2l.c \
+	w_log10.c w_log10f.c w_log10l.c \
+	w_pow.c w_powf.c w_powl.c \
+	w_sinh.c w_sinhf.c w_sinhl.c \
+	w_fmod.c w_fmodf.c w_fmodl.c \
+	w_sqrt.c w_sqrtf.c w_sqrtl.c \
+	w_remainder.c w_remainderf.c w_remainderl.c \
+	w_acosh.c w_acoshf.c w_acoshl.c \
+	w_atanh.c w_atanhf.c w_atanhl.c \
+	w_scalb.c w_scalbf.c w_scalbl.c 
+
+ifeq ($(UCLIBC_HAS_FENV),y)
+	libm_CSRC += k_standard.c k_standardf.c k_standardl.c
+endif
 
 FL_MOBJ := \
 	acosf.o \
@@ -230,7 +259,9 @@ libm_CSRC := \
 	s_expm1.c s_scalbn.c s_copysign.c e_acos.c e_asin.c e_atan2.c \
 	k_cos.c e_cosh.c e_exp.c e_fmod.c e_log.c e_log10.c e_pow.c \
 	k_sin.c e_sinh.c e_sqrt.c k_tan.c e_rem_pio2.c k_rem_pio2.c \
-	s_finite.c e_exp10.c
+	s_finite.c e_exp10.c 	\
+	matherr_wrapers.c k_standart.c 
+
 # We'll add sqrtf to avoid problems with libstdc++
 FL_MOBJ := sqrtf.o
 endif

+ 1 - 4
libm/e_acos.c

@@ -94,7 +94,4 @@ double __ieee754_acos(double x)
 	    w = r*s+c;
 	    return 2.0*(df+w);
 	}
-}
-
-strong_alias(__ieee754_acos, acos)
-libm_hidden_def(acos)
+}

+ 1 - 4
libm/e_acosh.c

@@ -50,9 +50,6 @@ double __ieee754_acosh(double x)
 	    return __ieee754_log(2.0*x-one/(x+__ieee754_sqrt(t-one)));
 	} else {			/* 1<x<2 */
 	    t = x-one;
-	    return log1p(t+sqrt(2.0*t+t*t));
+	    return log1p(t+__ieee754_sqrt(2.0*t+t*t));
 	}
 }
-
-strong_alias(__ieee754_acosh, acosh)
-libm_hidden_def(acosh)

+ 0 - 3
libm/e_asin.c

@@ -104,6 +104,3 @@ double __ieee754_asin(double x)
 	}
 	if(hx>0) return t; else return -t;
 }
-
-strong_alias(__ieee754_asin, asin)
-libm_hidden_def(asin)

+ 1 - 4
libm/e_atan2.c

@@ -113,7 +113,4 @@ double __ieee754_atan2(double y, double x)
 	    default: /* case 3 */
 	    	    return  (z-pi_lo)-pi;/* atan(-,-) */
 	}
-}
-
-strong_alias(__ieee754_atan2, atan2)
-libm_hidden_def(atan2)
+}

+ 0 - 3
libm/e_atanh.c

@@ -54,6 +54,3 @@ double __ieee754_atanh(double x)
 	    t = 0.5*log1p((x+x)/(one-x));
 	if(hx>=0) return t; else return -t;
 }
-
-strong_alias(__ieee754_atanh, atanh)
-libm_hidden_def(atanh)

+ 1 - 4
libm/e_cosh.c

@@ -76,7 +76,4 @@ double __ieee754_cosh(double x)
 
     /* |x| > overflowthresold, cosh(x) overflow */
 	return huge*huge;
-}
-
-strong_alias(__ieee754_cosh, cosh)
-libm_hidden_def(cosh)
+}

+ 0 - 3
libm/e_exp.c

@@ -155,6 +155,3 @@ double __ieee754_exp(double x)	/* default IEEE double exp */
 	    return y*twom1000;
 	}
 }
-
-strong_alias(__ieee754_exp, exp)
-libm_hidden_def(exp)

+ 0 - 2
libm/e_exp10.c

@@ -29,5 +29,3 @@ double __ieee754_exp10 (double arg)
        replaced sometime (soon?).  */
     return __ieee754_exp (M_LN10 * arg);
 }
-strong_alias(__ieee754_exp10, exp10)
-libm_hidden_def(exp10)

+ 0 - 3
libm/e_fmod.c

@@ -124,6 +124,3 @@ double __ieee754_fmod(double x, double y)
 	}
 	return x;		/* exact output */
 }
-
-strong_alias(__ieee754_fmod, fmod)
-libm_hidden_def(fmod)

+ 1 - 4
libm/e_hypot.c

@@ -115,7 +115,4 @@ double __ieee754_hypot(double x, double y)
 	    SET_HIGH_WORD(t1,high+(k<<20));
 	    return t1*w;
 	} else return w;
-}
-
-strong_alias(__ieee754_hypot, hypot)
-libm_hidden_def(hypot)
+}

+ 0 - 4
libm/e_j0.c

@@ -123,8 +123,6 @@ double __ieee754_j0(double x)
 	}
 }
 
-strong_alias(__ieee754_j0, j0)
-
 static const double
 u00  = -7.38042951086872317523e-02, /* 0xBFB2E4D6, 0x99CBD01F */
 u01  =  1.76666452509181115538e-01, /* 0x3FC69D01, 0x9DE9E3FC */
@@ -190,8 +188,6 @@ double __ieee754_y0(double x)
 	return(u/v + tpi*(__ieee754_j0(x)*__ieee754_log(x)));
 }
 
-strong_alias(__ieee754_y0, y0)
-
 /* The asymptotic expansions of pzero is
  *	1 - 9/128 s^2 + 11025/98304 s^4 - ...,	where s = 1/x.
  * For x >= 2, We approximate pzero by

+ 0 - 4
libm/e_j1.c

@@ -118,8 +118,6 @@ double __ieee754_j1(double x)
 	return(x*0.5+r/s);
 }
 
-strong_alias(__ieee754_j1, j1)
-
 static const double U0[5] = {
  -1.96057090646238940668e-01, /* 0xBFC91866, 0x143CBC8A */
   5.04438716639811282616e-02, /* 0x3FA9D3C7, 0x76292CD1 */
@@ -183,8 +181,6 @@ double __ieee754_y1(double x)
         return(x*(u/v) + tpi*(__ieee754_j1(x)*__ieee754_log(x)-one/x));
 }
 
-strong_alias(__ieee754_y1, y1)
-
 /* For x >= 8, the asymptotic expansions of pone is
  *	1 + 15/128 s^2 - 4725/2^15 s^4 - ...,	where s = 1/x.
  * We approximate pone by

+ 0 - 4
libm/e_jn.c

@@ -200,8 +200,6 @@ double __ieee754_jn(int n, double x)
 	if(sgn==1) return -b; else return b;
 }
 
-strong_alias(__ieee754_jn, jn)
-
 double __ieee754_yn(int n, double x)
 {
 	int32_t i,hx,ix,lx;
@@ -258,5 +256,3 @@ double __ieee754_yn(int n, double x)
 	}
 	if(sign>0) return b; else return -b;
 }
-
-strong_alias(__ieee754_yn, yn)

+ 2 - 28
libm/e_lgamma_r.c

@@ -295,36 +295,11 @@ double __ieee754_lgamma_r(double x, int *signgamp)
 	return r;
 }
 
-strong_alias(__ieee754_lgamma_r, lgamma_r)
-libm_hidden_def(lgamma_r)
-
-/* __ieee754_lgamma(x)
- * Return the logarithm of the Gamma function of x.
- */
-double __ieee754_lgamma(double x)
-{
-	return __ieee754_lgamma_r(x, &signgam);
-}
-
-strong_alias(__ieee754_lgamma, lgamma);
-libm_hidden_def(lgamma)
-
-
-/* NB: gamma function is an old name for lgamma.
- * It is deprecated.
- * Some C math libraries redefine it as a "true gamma", i.e.,
- * not a ln(|Gamma(x)|) but just Gamma(x), but standards
- * introduced tgamma name for that.
- */
-strong_alias(__ieee754_lgamma_r, gamma_r)
-strong_alias(__ieee754_lgamma, gamma)
-libm_hidden_def(gamma)
-
 
 /* double tgamma(double x)
  * Return the Gamma function of x.
  */
-double tgamma(double x)
+double __ieee754_tgamma(double x)
 {
 	int sign_of_gamma;
 	int32_t hx;
@@ -351,5 +326,4 @@ double tgamma(double x)
 
 	x = exp(lgamma_r(x, &sign_of_gamma));
 	return sign_of_gamma >= 0 ? x : -x;
-}
-libm_hidden_def(tgamma)
+}

+ 0 - 3
libm/e_log.c

@@ -127,6 +127,3 @@ double __ieee754_log(double x)
 		     return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
 	}
 }
-
-strong_alias(__ieee754_log, log)
-libm_hidden_def(log)

+ 0 - 3
libm/e_log10.c

@@ -78,6 +78,3 @@ double __ieee754_log10(double x)
 	z  = y*log10_2lo + ivln10*__ieee754_log(x);
 	return  z+y*log10_2hi;
 }
-
-strong_alias(__ieee754_log10, log10)
-libm_hidden_def(log10)

+ 0 - 2
libm/e_log2.c

@@ -115,5 +115,3 @@ double __ieee754_log2(double x)
 	    return dk-((s*(f-R))-f)/ln2;
 	}
 }
-strong_alias(__ieee754_log2,log2)
-libm_hidden_def(log2)

+ 0 - 3
libm/e_pow.c

@@ -299,6 +299,3 @@ double __ieee754_pow(double x, double y)
 	else SET_HIGH_WORD(z,j);
 	return s*z;
 }
-
-strong_alias(__ieee754_pow, pow)
-libm_hidden_def(pow)

+ 0 - 4
libm/e_remainder.c

@@ -63,7 +63,3 @@ double __ieee754_remainder(double x, double p)
 	SET_HIGH_WORD(x,hx^sx);
 	return x;
 }
-
-strong_alias(__ieee754_remainder, remainder)
-strong_alias(__ieee754_remainder, drem)
-libm_hidden_def(remainder)

+ 0 - 5
libm/e_scalb.c

@@ -31,8 +31,3 @@ double __ieee754_scalb(double x, double fn)
 	if (-fn > 65000.0) return scalbn(x,-65000);
 	return scalbn(x,(int)fn);
 }
-
-#if defined __UCLIBC_SUSV3_LEGACY__
-strong_alias(__ieee754_scalb, scalb)
-libm_hidden_def(scalb)
-#endif /* UCLIBC_SUSV3_LEGACY */

+ 0 - 3
libm/e_sinh.c

@@ -70,6 +70,3 @@ double __ieee754_sinh(double x)
     /* |x| > overflowthresold, sinh(x) overflow */
 	return x*shuge;
 }
-
-strong_alias(__ieee754_sinh, sinh)
-libm_hidden_def(sinh)

+ 0 - 3
libm/e_sqrt.c

@@ -180,9 +180,6 @@ double __ieee754_sqrt(double x)
 	return z;
 }
 
-strong_alias(__ieee754_sqrt, sqrt)
-libm_hidden_def(sqrt)
-
 
 /*
 Other methods  (use floating-point arithmetic)

+ 32 - 157
libm/float_wrappers.c

@@ -37,94 +37,78 @@ long long func##f (float x) \
 	return func((double)x); \
 }
 
-#ifndef __DO_XSI_MATH__
-# undef L_j0f	/* float j0f(float x); */
-# undef L_j1f	/* float j1f(float x); */
-# undef L_jnf	/* float jnf(int n, float x); */
-# undef L_y0f	/* float y0f(float x); */
-# undef L_y1f	/* float y1f(float x); */
-# undef L_ynf	/* float ynf(int n, float x); */
-#endif
-
 /* Implement the following, as defined by SuSv3 */
 #if 0
-float       acosf(float);
-float       acoshf(float);
-float       asinf(float);
 float       asinhf(float);
-float       atan2f(float, float);
 float       atanf(float);
-float       atanhf(float);
 float       cargf(float complex);
 float       cbrtf(float);
 float       ceilf(float);
 float       copysignf(float, float);
 float       cosf(float);
-float       coshf(float);
 float       erfcf(float);
 float       erff(float);
-float       exp2f(float);
-float       expf(float);
 float       expm1f(float);
 float       fabsf(float);
 float       floorf(float);
-float       fmodf(float, float);
 float       frexpf(float value, int *);
-float       hypotf(float, float);
 int         ilogbf(float);
 float       ldexpf(float, int);
-float       lgammaf(float);
 long long   llroundf(float);
-float       log10f(float);
 float       log1pf(float);
-float       log2f(float);
 float       logbf(float);
-float       logf(float);
 long        lroundf(float);
 float       modff(float, float *);
-float       powf(float, float);
-float       remainderf(float, float);
 float       rintf(float);
 float       roundf(float);
 float       scalbnf(float, int);
 float       sinf(float);
-float       sinhf(float);
-float       sqrtf(float);
 float       tanf(float);
 float       tanhf(float);
 #endif
 
-#ifdef L_acosf
-WRAPPER1(acos)
-#endif
-
-#ifdef L_acoshf
-WRAPPER1(acosh)
-#endif
-
-#ifdef L_asinf
-WRAPPER1(asin)
+/* The following functions implemented as wrappers
+ * in separate files (w_funcf.c)
+ */
+#if 0
+float       acosf(float);
+float       acoshf(float);
+float       asinf(float);
+float       atan2f(float, float);
+float       atanhf(float);
+float       coshf(float);
+float       exp2f(float);
+float       expf(float);
+float       fmodf(float, float);
+float       hypotf(float, float);
+float       lgammaf(float);
+float       log10f(float);
+float       log2f(float);
+float       logf(float);
+float       powf(float, float);
+float       remainderf(float, float);
+float       sinhf(float);
+float       sqrtf(float);
+float		j0f(float x);
+float		j1f(float x);
+float		jnf(int n, float x);
+float		y0f(float x);
+float		y1f(float x);
+float		ynf(int n, float x);
+float 		tgammaf(float x);
+float 		scalbf(float x, float fn);
+float 		gammaf(float x);
+float	 	scalbl(float x, float fn);
 #endif
 
 #ifdef L_asinhf
 WRAPPER1(asinh)
 #endif
 
-#ifdef L_atan2f
-float atan2f (float x, float y)
-{
-	return (float) atan2( (double)x, (double)y );
-}
-#endif
-
 #ifdef L_atanf
 WRAPPER1(atan)
 #endif
 
-#ifdef L_atanhf
-WRAPPER1(atanh)
-#endif
-
 #ifdef L_cargf
 float cargf (float complex x)
 {
@@ -152,10 +136,6 @@ WRAPPER1(cos)
 libm_hidden_def(cosf)
 #endif
 
-#ifdef L_coshf
-WRAPPER1(cosh)
-#endif
-
 #ifdef L_erfcf
 WRAPPER1(erfc)
 #endif
@@ -164,14 +144,6 @@ WRAPPER1(erfc)
 WRAPPER1(erf)
 #endif
 
-#ifdef L_exp2f
-WRAPPER1(exp2)
-#endif
-
-#ifdef L_expf
-WRAPPER1(exp)
-#endif
-
 #ifdef L_expm1f
 WRAPPER1(expm1)
 #endif
@@ -212,13 +184,6 @@ float fminf (float x, float y)
 }
 #endif
 
-#ifdef L_fmodf
-float fmodf (float x, float y)
-{
-	return (float) fmod( (double)x, (double)y );
-}
-#endif
-
 #ifdef L_frexpf
 float frexpf (float x, int *_exp)
 {
@@ -226,32 +191,10 @@ float frexpf (float x, int *_exp)
 }
 #endif
 
-#ifdef L_hypotf
-float hypotf (float x, float y)
-{
-	return (float) hypot( (double)x, (double)y );
-}
-#endif
-
 #ifdef L_ilogbf
 int_WRAPPER1(ilogb)
 #endif
 
-#ifdef L_j0f
-WRAPPER1(j0)
-#endif
-
-#ifdef L_j1f
-WRAPPER1(j1)
-#endif
-
-#ifdef L_jnf
-float jnf(int n, float x)
-{
-	return (float) jn(n, (double)x);
-}
-#endif
-
 #ifdef L_ldexpf
 float ldexpf (float x, int _exp)
 {
@@ -259,10 +202,6 @@ float ldexpf (float x, int _exp)
 }
 #endif
 
-#ifdef L_lgammaf
-WRAPPER1(lgamma)
-#endif
-
 #ifdef L_llrintf
 long_long_WRAPPER1(llrint)
 #endif
@@ -271,26 +210,14 @@ long_long_WRAPPER1(llrint)
 long_long_WRAPPER1(llround)
 #endif
 
-#ifdef L_log10f
-WRAPPER1(log10)
-#endif
-
 #ifdef L_log1pf
 WRAPPER1(log1p)
 #endif
 
-#ifdef L_log2f
-WRAPPER1(log2)
-#endif
-
 #ifdef L_logbf
 WRAPPER1(logb)
 #endif
 
-#ifdef L_logf
-WRAPPER1(log)
-#endif
-
 #ifdef L_lrintf
 long_WRAPPER1(lrint)
 #endif
@@ -320,20 +247,6 @@ float nexttowardf (float x, long double y)
 }
 #endif
 
-#ifdef L_powf
-float powf (float x, float y)
-{
-	return (float) pow( (double)x, (double)y );
-}
-#endif
-
-#ifdef L_remainderf
-float remainderf (float x, float y)
-{
-	return (float) remainder( (double)x, (double)y );
-}
-#endif
-
 #ifdef L_remquof
 float remquof (float x, float y, int *quo)
 {
@@ -368,14 +281,6 @@ WRAPPER1(sin)
 libm_hidden_def(sinf)
 #endif
 
-#ifdef L_sinhf
-WRAPPER1(sinh)
-#endif
-
-#ifdef L_sqrtf
-WRAPPER1(sqrt)
-#endif
-
 #ifdef L_tanf
 WRAPPER1(tan)
 #endif
@@ -384,40 +289,10 @@ WRAPPER1(tan)
 WRAPPER1(tanh)
 #endif
 
-#ifdef L_tgammaf
-WRAPPER1(tgamma)
-#endif
-
 #ifdef L_truncf
 WRAPPER1(trunc)
 #endif
 
-#if defined L_scalbf && defined __UCLIBC_SUSV3_LEGACY__
-float scalbf (float x, float y)
-{
-	return (float) scalb( (double)x, (double)y );
-}
-#endif
-
-#ifdef L_gammaf
-WRAPPER1(gamma)
-#endif
-
 #ifdef L_significandf
 WRAPPER1(significand)
 #endif
-
-#ifdef L_y0f
-WRAPPER1(y0)
-#endif
-
-#ifdef L_y1f
-WRAPPER1(y1)
-#endif
-
-#ifdef L_ynf
-float ynf(int n, float x)
-{
-	return (float) yn(n, (double)x);
-}
-#endif

+ 973 - 0
libm/k_standard.c

@@ -0,0 +1,973 @@
+/* @(#)k_standard.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#if defined(LIBM_SCCS) && !defined(lint)
+static char rcsid[] = "$NetBSD: k_standard.c,v 1.6 1995/05/10 20:46:35 jtc Exp $";
+#endif
+
+#include <math.h>
+#include "math_private.h"
+#include <errno.h>
+
+#include <assert.h>
+
+#ifndef _USE_WRITE
+#include <stdio.h>			/* fputs(), stderr */
+#define	WRITE2(u,v)	fputs(u, stderr)
+#else	/* !defined(_USE_WRITE) */
+#include <unistd.h>			/* write */
+#define	WRITE2(u,v)	write(2, u, v)
+#undef fflush
+#endif	/* !defined(_USE_WRITE) */
+
+/* XXX gcc versions until now don't delay the 0.0/0.0 division until
+   runtime but produce NaN at compile time.  This is wrong since the
+   exceptions are not set correctly.  */
+#if 0
+static const double zero = 0.0;	/* used as const */
+#else
+static double zero = 0.0;	/* used as const */
+#endif
+
+/*
+ * Standard conformance (non-IEEE) on exception cases.
+ * Mapping:
+ *	1 -- acos(|x|>1)
+ *	2 -- asin(|x|>1)
+ *	3 -- atan2(+-0,+-0)
+ *	4 -- hypot overflow
+ *	5 -- cosh overflow
+ *	6 -- exp overflow
+ *	7 -- exp underflow
+ *	8 -- y0(0)
+ *	9 -- y0(-ve)
+ *	10-- y1(0)
+ *	11-- y1(-ve)
+ *	12-- yn(0)
+ *	13-- yn(-ve)
+ *	14-- lgamma(finite) overflow
+ *	15-- lgamma(-integer)
+ *	16-- log(0)
+ *	17-- log(x<0)
+ *	18-- log10(0)
+ *	19-- log10(x<0)
+ *	20-- pow(0.0,0.0)
+ *	21-- pow(x,y) overflow
+ *	22-- pow(x,y) underflow
+ *	23-- pow(0,negative)
+ *	24-- pow(neg,non-integral)
+ *	25-- sinh(finite) overflow
+ *	26-- sqrt(negative)
+ *	27-- fmod(x,0)
+ *	28-- remainder(x,0)
+ *	29-- acosh(x<1)
+ *	30-- atanh(|x|>1)
+ *	31-- atanh(|x|=1)
+ *	32-- scalb overflow
+ *	33-- scalb underflow
+ *	34-- j0(|x|>X_TLOSS)
+ *	35-- y0(x>X_TLOSS)
+ *	36-- j1(|x|>X_TLOSS)
+ *	37-- y1(x>X_TLOSS)
+ *	38-- jn(|x|>X_TLOSS, n)
+ *	39-- yn(x>X_TLOSS, n)
+ *	40-- tgamma(finite) overflow
+ *	41-- tgamma(-integer)
+ *	42-- pow(NaN,0.0)
+ *	43-- +0**neg
+ *	44-- exp2 overflow	done
+ *	45-- exp2 underflow	done
+ *	46-- exp10 overflow
+ *	47-- exp10 underflow
+ *	48-- log2(0)
+ *	49-- log2(x<0)
+ *	50-- tgamma(+-0)
+ */
+
+
+double
+__kernel_standard(double x, double y, int type)
+{
+	struct exception exc;
+#ifndef HUGE_VAL	/* this is the only routine that uses HUGE_VAL */
+#define HUGE_VAL inf
+	double inf = 0.0;
+
+	SET_HIGH_WORD(inf,0x7ff00000);	/* set inf to infinite */
+#endif
+
+	/* The SVID struct exception uses a field "char *name;".  */
+#define CSTR(func) ((char *) (type < 100				\
+			      ? func					\
+			      : (type < 200 ? func "f" : func "l")))
+
+#ifdef _USE_WRITE
+	(void) fflush(stdout);
+#endif
+	exc.arg1 = x;
+	exc.arg2 = y;
+	switch(type) {
+	    case 1:
+	    case 101:
+	    case 201:
+		/* acos(|x|>1) */
+		exc.type = DOMAIN;
+		exc.name = CSTR ("acos");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = HUGE;
+		else
+		  exc.retval = NAN;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (EDOM);
+		else if (!matherr(&exc)) {
+		  if(_LIB_VERSION == _SVID_) {
+		    (void) WRITE2("acos: DOMAIN error\n", 19);
+		  }
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 2:
+	    case 102:
+	    case 202:
+		/* asin(|x|>1) */
+		exc.type = DOMAIN;
+		exc.name = CSTR ("asin");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = HUGE;
+		else
+		  exc.retval = NAN;
+		if(_LIB_VERSION == _POSIX_)
+		  __set_errno (EDOM);
+		else if (!matherr(&exc)) {
+		  if(_LIB_VERSION == _SVID_) {
+			(void) WRITE2("asin: DOMAIN error\n", 19);
+		  }
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 3:
+	    case 103:
+	    case 203:
+		/* atan2(+-0,+-0) */
+		exc.arg1 = y;
+		exc.arg2 = x;
+		exc.type = DOMAIN;
+		exc.name = CSTR ("atan2");
+		assert (_LIB_VERSION == _SVID_);
+		exc.retval = HUGE;
+		if(_LIB_VERSION == _POSIX_)
+		  __set_errno (EDOM);
+		else if (!matherr(&exc)) {
+		  if(_LIB_VERSION == _SVID_) {
+			(void) WRITE2("atan2: DOMAIN error\n", 20);
+		      }
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 4:
+	    case 104:
+	    case 204:
+		/* hypot(finite,finite) overflow */
+		exc.type = OVERFLOW;
+		exc.name = CSTR ("hypot");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = HUGE;
+		else
+		  exc.retval = HUGE_VAL;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+			__set_errno (ERANGE);
+		}
+		break;
+	    case 5:
+	    case 105:
+	    case 205:
+		/* cosh(finite) overflow */
+		exc.type = OVERFLOW;
+		exc.name = CSTR ("cosh");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = HUGE;
+		else
+		  exc.retval = HUGE_VAL;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+			__set_errno (ERANGE);
+		}
+		break;
+	    case 6:
+	    case 106:
+	    case 206:
+		/* exp(finite) overflow */
+		exc.type = OVERFLOW;
+		exc.name = CSTR ("exp");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = HUGE;
+		else
+		  exc.retval = HUGE_VAL;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+			__set_errno (ERANGE);
+		}
+		break;
+	    case 7:
+	    case 107:
+	    case 207:
+		/* exp(finite) underflow */
+		exc.type = UNDERFLOW;
+		exc.name = CSTR ("exp");
+		exc.retval = zero;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+			__set_errno (ERANGE);
+		}
+		break;
+	    case 8:
+	    case 108:
+	    case 208:
+		/* y0(0) = -inf */
+		exc.type = DOMAIN;	/* should be SING for IEEE */
+		exc.name = CSTR ("y0");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = -HUGE;
+		else
+		  exc.retval = -HUGE_VAL;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+		  if (_LIB_VERSION == _SVID_) {
+			(void) WRITE2("y0: DOMAIN error\n", 17);
+		      }
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 9:
+	    case 109:
+	    case 209:
+		/* y0(x<0) = NaN */
+		exc.type = DOMAIN;
+		exc.name = CSTR ("y0");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = -HUGE;
+		else
+		  exc.retval = NAN;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (EDOM);
+		else if (!matherr(&exc)) {
+		  if (_LIB_VERSION == _SVID_) {
+			(void) WRITE2("y0: DOMAIN error\n", 17);
+		      }
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 10:
+	    case 110:
+	    case 210:
+		/* y1(0) = -inf */
+		exc.type = DOMAIN;	/* should be SING for IEEE */
+		exc.name = CSTR ("y1");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = -HUGE;
+		else
+		  exc.retval = -HUGE_VAL;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+		  if (_LIB_VERSION == _SVID_) {
+			(void) WRITE2("y1: DOMAIN error\n", 17);
+		      }
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 11:
+	    case 111:
+	    case 211:
+		/* y1(x<0) = NaN */
+		exc.type = DOMAIN;
+		exc.name = CSTR ("y1");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = -HUGE;
+		else
+		  exc.retval = NAN;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (EDOM);
+		else if (!matherr(&exc)) {
+		  if (_LIB_VERSION == _SVID_) {
+			(void) WRITE2("y1: DOMAIN error\n", 17);
+		      }
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 12:
+	    case 112:
+	    case 212:
+		/* yn(n,0) = -inf */
+		exc.type = DOMAIN;	/* should be SING for IEEE */
+		exc.name = CSTR ("yn");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = -HUGE;
+		else
+		  exc.retval = ((x < 0 && ((int) x & 1) != 0)
+				? HUGE_VAL
+				: -HUGE_VAL);
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+		  if (_LIB_VERSION == _SVID_) {
+			(void) WRITE2("yn: DOMAIN error\n", 17);
+		      }
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 13:
+	    case 113:
+	    case 213:
+		/* yn(x<0) = NaN */
+		exc.type = DOMAIN;
+		exc.name = CSTR ("yn");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = -HUGE;
+		else
+		  exc.retval = NAN;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (EDOM);
+		else if (!matherr(&exc)) {
+		  if (_LIB_VERSION == _SVID_) {
+			(void) WRITE2("yn: DOMAIN error\n", 17);
+		      }
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 14:
+	    case 114:
+	    case 214:
+		/* lgamma(finite) overflow */
+		exc.type = OVERFLOW;
+		exc.name = CSTR ("lgamma");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = HUGE;
+		else
+		  exc.retval = HUGE_VAL;
+		if (_LIB_VERSION == _POSIX_)
+			__set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+			__set_errno (ERANGE);
+		}
+		break;
+	    case 15:
+	    case 115:
+	    case 215:
+		/* lgamma(-integer) or lgamma(0) */
+		exc.type = SING;
+		exc.name = CSTR ("lgamma");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = HUGE;
+		else
+		  exc.retval = HUGE_VAL;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+		  if (_LIB_VERSION == _SVID_) {
+			(void) WRITE2("lgamma: SING error\n", 19);
+		      }
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 16:
+	    case 116:
+	    case 216:
+		/* log(0) */
+		exc.type = SING;
+		exc.name = CSTR ("log");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = -HUGE;
+		else
+		  exc.retval = -HUGE_VAL;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+		  if (_LIB_VERSION == _SVID_) {
+			(void) WRITE2("log: SING error\n", 16);
+		      }
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 17:
+	    case 117:
+	    case 217:
+		/* log(x<0) */
+		exc.type = DOMAIN;
+		exc.name = CSTR ("log");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = -HUGE;
+		else
+		  exc.retval = NAN;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (EDOM);
+		else if (!matherr(&exc)) {
+		  if (_LIB_VERSION == _SVID_) {
+			(void) WRITE2("log: DOMAIN error\n", 18);
+		      }
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 18:
+	    case 118:
+	    case 218:
+		/* log10(0) */
+		exc.type = SING;
+		exc.name = CSTR ("log10");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = -HUGE;
+		else
+		  exc.retval = -HUGE_VAL;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+		  if (_LIB_VERSION == _SVID_) {
+			(void) WRITE2("log10: SING error\n", 18);
+		      }
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 19:
+	    case 119:
+	    case 219:
+		/* log10(x<0) */
+		exc.type = DOMAIN;
+		exc.name = CSTR ("log10");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = -HUGE;
+		else
+		  exc.retval = NAN;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (EDOM);
+		else if (!matherr(&exc)) {
+		  if (_LIB_VERSION == _SVID_) {
+			(void) WRITE2("log10: DOMAIN error\n", 20);
+		      }
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 20:
+	    case 120:
+	    case 220:
+		/* pow(0.0,0.0) */
+		/* error only if _LIB_VERSION == _SVID_ */
+		exc.type = DOMAIN;
+		exc.name = CSTR ("pow");
+		exc.retval = zero;
+		if (_LIB_VERSION != _SVID_) exc.retval = 1.0;
+		else if (!matherr(&exc)) {
+			(void) WRITE2("pow(0,0): DOMAIN error\n", 23);
+			__set_errno (EDOM);
+		}
+		break;
+	    case 21:
+	    case 121:
+	    case 221:
+		/* pow(x,y) overflow */
+		exc.type = OVERFLOW;
+		exc.name = CSTR ("pow");
+		if (_LIB_VERSION == _SVID_) {
+		  exc.retval = HUGE;
+		  y *= 0.5;
+		  if(x<zero&&rint(y)!=y) exc.retval = -HUGE;
+		} else {
+		  exc.retval = HUGE_VAL;
+		  y *= 0.5;
+		  if(x<zero&&rint(y)!=y) exc.retval = -HUGE_VAL;
+		}
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+			__set_errno (ERANGE);
+		}
+		break;
+	    case 22:
+	    case 122:
+	    case 222:
+		/* pow(x,y) underflow */
+		exc.type = UNDERFLOW;
+		exc.name = CSTR ("pow");
+		exc.retval =  zero;
+		y *= 0.5;
+		if (x < zero && rint (y) != y)
+		  exc.retval = -zero;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+			__set_errno (ERANGE);
+		}
+		break;
+	    case 23:
+	    case 123:
+	    case 223:
+		/* -0**neg */
+		exc.type = DOMAIN;
+		exc.name = CSTR ("pow");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = zero;
+		else
+		  exc.retval = -HUGE_VAL;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+		  if (_LIB_VERSION == _SVID_) {
+			(void) WRITE2("pow(0,neg): DOMAIN error\n", 25);
+		      }
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 43:
+	    case 143:
+	    case 243:
+		/* +0**neg */
+		exc.type = DOMAIN;
+		exc.name = CSTR ("pow");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = zero;
+		else
+		  exc.retval = HUGE_VAL;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+		  if (_LIB_VERSION == _SVID_) {
+			(void) WRITE2("pow(0,neg): DOMAIN error\n", 25);
+		      }
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 24:
+	    case 124:
+	    case 224:
+		/* neg**non-integral */
+		exc.type = DOMAIN;
+		exc.name = CSTR ("pow");
+		if (_LIB_VERSION == _SVID_)
+		    exc.retval = zero;
+		else
+		    exc.retval = zero/zero;	/* X/Open allow NaN */
+		if (_LIB_VERSION == _POSIX_)
+		   __set_errno (EDOM);
+		else if (!matherr(&exc)) {
+		  if (_LIB_VERSION == _SVID_) {
+			(void) WRITE2("neg**non-integral: DOMAIN error\n", 32);
+		      }
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 25:
+	    case 125:
+	    case 225:
+		/* sinh(finite) overflow */
+		exc.type = OVERFLOW;
+		exc.name = CSTR ("sinh");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = ( (x>zero) ? HUGE : -HUGE);
+		else
+		  exc.retval = ( (x>zero) ? HUGE_VAL : -HUGE_VAL);
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+			__set_errno (ERANGE);
+		}
+		break;
+	    case 26:
+	    case 126:
+	    case 226:
+		/* sqrt(x<0) */
+		exc.type = DOMAIN;
+		exc.name = CSTR ("sqrt");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = zero;
+		else
+		  exc.retval = zero/zero;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (EDOM);
+		else if (!matherr(&exc)) {
+		  if (_LIB_VERSION == _SVID_) {
+			(void) WRITE2("sqrt: DOMAIN error\n", 19);
+		      }
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 27:
+	    case 127:
+	    case 227:
+		/* fmod(x,0) */
+		exc.type = DOMAIN;
+		exc.name = CSTR ("fmod");
+		if (_LIB_VERSION == _SVID_)
+		    exc.retval = x;
+		else
+		    exc.retval = zero/zero;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (EDOM);
+		else if (!matherr(&exc)) {
+		  if (_LIB_VERSION == _SVID_) {
+		    (void) WRITE2("fmod:  DOMAIN error\n", 20);
+		  }
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 28:
+	    case 128:
+	    case 228:
+		/* remainder(x,0) */
+		exc.type = DOMAIN;
+		exc.name = CSTR ("remainder");
+		exc.retval = zero/zero;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (EDOM);
+		else if (!matherr(&exc)) {
+		  if (_LIB_VERSION == _SVID_) {
+		    (void) WRITE2("remainder: DOMAIN error\n", 24);
+		  }
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 29:
+	    case 129:
+	    case 229:
+		/* acosh(x<1) */
+		exc.type = DOMAIN;
+		exc.name = CSTR ("acosh");
+		exc.retval = zero/zero;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (EDOM);
+		else if (!matherr(&exc)) {
+		  if (_LIB_VERSION == _SVID_) {
+		    (void) WRITE2("acosh: DOMAIN error\n", 20);
+		  }
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 30:
+	    case 130:
+	    case 230:
+		/* atanh(|x|>1) */
+		exc.type = DOMAIN;
+		exc.name = CSTR ("atanh");
+		exc.retval = zero/zero;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (EDOM);
+		else if (!matherr(&exc)) {
+		  if (_LIB_VERSION == _SVID_) {
+		    (void) WRITE2("atanh: DOMAIN error\n", 20);
+		  }
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 31:
+	    case 131:
+	    case 231:
+		/* atanh(|x|=1) */
+		exc.type = SING;
+		exc.name = CSTR ("atanh");
+		exc.retval = x/zero;	/* sign(x)*inf */
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+		  if (_LIB_VERSION == _SVID_) {
+		    (void) WRITE2("atanh: SING error\n", 18);
+		  }
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 32:
+	    case 132:
+	    case 232:
+		/* scalb overflow; SVID also returns +-HUGE_VAL */
+		exc.type = OVERFLOW;
+		exc.name = CSTR ("scalb");
+		exc.retval = x > zero ? HUGE_VAL : -HUGE_VAL;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+			__set_errno (ERANGE);
+		}
+		break;
+	    case 33:
+	    case 133:
+	    case 233:
+		/* scalb underflow */
+		exc.type = UNDERFLOW;
+		exc.name = CSTR ("scalb");
+		exc.retval = copysign(zero,x);
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+			__set_errno (ERANGE);
+		}
+		break;
+	    case 34:
+	    case 134:
+	    case 234:
+		/* j0(|x|>X_TLOSS) */
+		exc.type = TLOSS;
+		exc.name = CSTR ("j0");
+		exc.retval = zero;
+		if (_LIB_VERSION == _POSIX_)
+			__set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+			if (_LIB_VERSION == _SVID_) {
+				(void) WRITE2(exc.name, 2);
+				(void) WRITE2(": TLOSS error\n", 14);
+			}
+			__set_errno (ERANGE);
+		}
+		break;
+	    case 35:
+	    case 135:
+	    case 235:
+		/* y0(x>X_TLOSS) */
+		exc.type = TLOSS;
+		exc.name = CSTR ("y0");
+		exc.retval = zero;
+		if (_LIB_VERSION == _POSIX_)
+			__set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+			if (_LIB_VERSION == _SVID_) {
+				(void) WRITE2(exc.name, 2);
+				(void) WRITE2(": TLOSS error\n", 14);
+			}
+			__set_errno (ERANGE);
+		}
+		break;
+	    case 36:
+	    case 136:
+	    case 236:
+		/* j1(|x|>X_TLOSS) */
+		exc.type = TLOSS;
+		exc.name = CSTR ("j1");
+		exc.retval = zero;
+		if (_LIB_VERSION == _POSIX_)
+			__set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+			if (_LIB_VERSION == _SVID_) {
+				(void) WRITE2(exc.name, 2);
+				(void) WRITE2(": TLOSS error\n", 14);
+			}
+			__set_errno (ERANGE);
+		}
+		break;
+	    case 37:
+	    case 137:
+	    case 237:
+		/* y1(x>X_TLOSS) */
+		exc.type = TLOSS;
+		exc.name = CSTR ("y1");
+		exc.retval = zero;
+		if (_LIB_VERSION == _POSIX_)
+			__set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+			if (_LIB_VERSION == _SVID_) {
+				(void) WRITE2(exc.name, 2);
+				(void) WRITE2(": TLOSS error\n", 14);
+			}
+			__set_errno (ERANGE);
+		}
+		break;
+	    case 38:
+	    case 138:
+	    case 238:
+		/* jn(|x|>X_TLOSS) */
+		exc.type = TLOSS;
+		exc.name = CSTR ("jn");
+		exc.retval = zero;
+		if (_LIB_VERSION == _POSIX_)
+			__set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+			if (_LIB_VERSION == _SVID_) {
+				(void) WRITE2(exc.name, 2);
+				(void) WRITE2(": TLOSS error\n", 14);
+			}
+			__set_errno (ERANGE);
+		}
+		break;
+	    case 39:
+	    case 139:
+	    case 239:
+		/* yn(x>X_TLOSS) */
+		exc.type = TLOSS;
+		exc.name = CSTR ("yn");
+		exc.retval = zero;
+		if (_LIB_VERSION == _POSIX_)
+			__set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+			if (_LIB_VERSION == _SVID_) {
+				(void) WRITE2(exc.name, 2);
+				(void) WRITE2(": TLOSS error\n", 14);
+			}
+			__set_errno (ERANGE);
+		}
+		break;
+	    case 40:
+	    case 140:
+	    case 240:
+		/* tgamma(finite) overflow */
+		exc.type = OVERFLOW;
+		exc.name = CSTR ("tgamma");
+		exc.retval = copysign (HUGE_VAL, x);
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+		  __set_errno (ERANGE);
+		}
+		break;
+	    case 41:
+	    case 141:
+	    case 241:
+		/* tgamma(-integer) */
+		exc.type = SING;
+		exc.name = CSTR ("tgamma");
+		exc.retval = NAN;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (EDOM);
+		else if (!matherr(&exc)) {
+		  if (_LIB_VERSION == _SVID_) {
+			(void) WRITE2("tgamma: SING error\n", 18);
+			exc.retval = HUGE_VAL;
+		      }
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 42:
+	    case 142:
+	    case 242:
+		/* pow(NaN,0.0) */
+		/* error only if _LIB_VERSION == _SVID_ & _XOPEN_ */
+		exc.type = DOMAIN;
+		exc.name = CSTR ("pow");
+		exc.retval = x;
+		if (_LIB_VERSION == _IEEE_ ||
+		    _LIB_VERSION == _POSIX_) exc.retval = 1.0;
+		else if (!matherr(&exc)) {
+			__set_errno (EDOM);
+		}
+		break;
+
+	    case 44:
+	    case 144:
+	    case 244:
+		/* exp(finite) overflow */
+		exc.type = OVERFLOW;
+		exc.name = CSTR ("exp2");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = HUGE;
+		else
+		  exc.retval = HUGE_VAL;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+			__set_errno (ERANGE);
+		}
+		break;
+	    case 45:
+	    case 145:
+	    case 245:
+		/* exp(finite) underflow */
+		exc.type = UNDERFLOW;
+		exc.name = CSTR ("exp2");
+		exc.retval = zero;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+			__set_errno (ERANGE);
+		}
+		break;
+
+	    case 46:
+	    case 146:
+	    case 246:
+		/* exp(finite) overflow */
+		exc.type = OVERFLOW;
+		exc.name = CSTR ("exp10");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = HUGE;
+		else
+		  exc.retval = HUGE_VAL;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+			__set_errno (ERANGE);
+		}
+		break;
+	    case 47:
+	    case 147:
+	    case 247:
+		/* exp(finite) underflow */
+		exc.type = UNDERFLOW;
+		exc.name = CSTR ("exp10");
+		exc.retval = zero;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+			__set_errno (ERANGE);
+		}
+		break;
+	    case 48:
+	    case 148:
+	    case 248:
+		/* log2(0) */
+		exc.type = SING;
+		exc.name = CSTR ("log2");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = -HUGE;
+		else
+		  exc.retval = -HUGE_VAL;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 49:
+	    case 149:
+	    case 249:
+		/* log2(x<0) */
+		exc.type = DOMAIN;
+		exc.name = CSTR ("log2");
+		if (_LIB_VERSION == _SVID_)
+		  exc.retval = -HUGE;
+		else
+		  exc.retval = NAN;
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (EDOM);
+		else if (!matherr(&exc)) {
+		  __set_errno (EDOM);
+		}
+		break;
+	    case 50:
+	    case 150:
+	    case 250:
+		/* tgamma(+-0) */
+		exc.type = SING;
+		exc.name = CSTR ("tgamma");
+		exc.retval = copysign (HUGE_VAL, x);
+		if (_LIB_VERSION == _POSIX_)
+		  __set_errno (ERANGE);
+		else if (!matherr(&exc)) {
+		  if (_LIB_VERSION == _SVID_)
+		    (void) WRITE2("tgamma: SING error\n", 18);
+		  __set_errno (ERANGE);
+		}
+		break;
+
+		/* #### Last used is 50/150/250 ### */
+	}
+	return exc.retval;
+}

+ 31 - 0
libm/k_standardf.c

@@ -0,0 +1,31 @@
+/* Implement __kernel_standard_f.
+   Copyright (C) 2011-2016 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
+
+#include <math.h>
+#include "math_private.h"
+
+
+/* Handle errors for a libm function as specified by TYPE (see
+   comments in k_standard.c for details), with arguments X and Y,
+   returning the appropriate return value for that function.  */
+
+float
+__kernel_standard_f (float x, float y, int type)
+{
+  return __kernel_standard (x, y, type);
+}

+ 107 - 0
libm/k_standardl.c

@@ -0,0 +1,107 @@
+/* Implement __kernel_standard_l.
+   Copyright (C) 2012-2016 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.
+
+   Parts based on k_standard.c from fdlibm: */
+
+/* @(#)k_standard.c 5.1 93/09/24 */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include <math.h>
+#include "math_private.h"
+#include <fenv.h>
+#include <float.h>
+#include <errno.h>
+
+
+static double zero = 0.0;
+
+/* Handle errors for a libm function as specified by TYPE (see
+   comments in k_standard.c for details), with arguments X and Y,
+   returning the appropriate return value for that function.  */
+
+long double
+__kernel_standard_l (long double x, long double y, int type)
+{
+  double dx, dy;
+  struct exception exc;
+  fenv_t env;
+
+  feholdexcept (&env);
+  dx = x;
+  dy = y;
+  math_force_eval (dx);
+  math_force_eval (dy);
+  fesetenv (&env);
+
+  switch (type)
+    {
+    case 221:
+      /* powl (x, y) overflow.  */
+      exc.arg1 = dx;
+      exc.arg2 = dy;
+      exc.type = OVERFLOW;
+      exc.name = (char *) "powl";
+      if (_LIB_VERSION == _SVID_)
+	{
+	  exc.retval = HUGE;
+	  y *= 0.5;
+	  if (x < zero && rint (y) != y)
+	    exc.retval = -HUGE;
+	}
+      else
+	{
+	  exc.retval = HUGE_VAL;
+	  y *= 0.5;
+	  if (x < zero && rint (y) != y)
+	    exc.retval = -HUGE_VAL;
+	}
+      if (_LIB_VERSION == _POSIX_)
+	__set_errno (ERANGE);
+      else if (!matherr (&exc))
+	__set_errno (ERANGE);
+      return exc.retval;
+
+    case 222:
+      /* powl (x, y) underflow.  */
+      exc.arg1 = dx;
+      exc.arg2 = dy;
+      exc.type = UNDERFLOW;
+      exc.name = (char *) "powl";
+      exc.retval = zero;
+      y *= 0.5;
+      if (x < zero && rint (y) != y)
+	exc.retval = -zero;
+      if (_LIB_VERSION == _POSIX_)
+	__set_errno (ERANGE);
+      else if (!matherr (&exc))
+	__set_errno (ERANGE);
+      return exc.retval;
+
+    default:
+      return __kernel_standard (dx, dy, type);
+    }
+}

+ 32 - 139
libm/ldouble_wrappers.c

@@ -42,34 +42,17 @@ long long func##l(long double x) \
 	return func((double) x); \
 }
 
-#ifndef __DO_XSI_MATH__
-# undef L_j0l  /* long double j0l(long double x); */
-# undef L_j1l  /* long double j1l(long double x); */
-# undef L_jnl  /* long double jnl(int n, long double x); */
-# undef L_y0l  /* long double y0l(long double x); */
-# undef L_y1l  /* long double y1l(long double x); */
-# undef L_ynl  /* long double ynl(int n, long double x); */
-#endif
-
 /* Implement the following, as defined by SuSv3 */
 #if 0
-long double acoshl(long double);
-long double acosl(long double);
 long double asinhl(long double);
-long double asinl(long double);
-long double atan2l(long double, long double);
-long double atanhl(long double);
 long double atanl(long double);
 long double cargl(long double complex);
 long double cbrtl(long double);
 long double ceill(long double);
 long double copysignl(long double, long double);
-long double coshl(long double);
 long double cosl(long double);
 long double erfcl(long double);
 long double erfl(long double);
-long double exp2l(long double);
-long double expl(long double);
 long double expm1l(long double);
 long double fabsl(long double);
 long double fdiml(long double, long double);
@@ -77,65 +60,69 @@ long double floorl(long double);
 long double fmal(long double, long double, long double);
 long double fmaxl(long double, long double);
 long double fminl(long double, long double);
-long double fmodl(long double, long double);
 long double frexpl(long double value, int *);
-long double hypotl(long double, long double);
 int         ilogbl(long double);
 long double ldexpl(long double, int);
-long double lgammal(long double);
 long long   llrintl(long double);
 long long   llroundl(long double);
-long double log10l(long double);
 long double log1pl(long double);
-long double log2l(long double);
 long double logbl(long double);
-long double logl(long double);
 long        lrintl(long double);
 long        lroundl(long double);
 long double modfl(long double, long double *);
 long double nearbyintl(long double);
 long double nextafterl(long double, long double);
 long double nexttowardl(long double, long double);
-long double powl(long double, long double);
-long double remainderl(long double, long double);
 long double remquol(long double, long double, int *);
 long double rintl(long double);
 long double roundl(long double);
 long double scalblnl(long double, long);
 long double scalbnl(long double, int);
-long double sinhl(long double);
 long double sinl(long double);
-long double sqrtl(long double);
 long double tanhl(long double);
 long double tanl(long double);
-long double tgammal(long double);
 long double truncl(long double);
 #endif
 
-#ifdef L_acoshl
-WRAPPER1(acosh)
-#endif
+/* The following functions implemented as wrappers
+ * in separate files (w_funcl.c)
+ */
+#if 0
+long double 	acosl(long double);
+long double 	acoshl(long double);
+long double 	asinl(long double);
+long double 	atan2l(long double, long double);
+long double 	atanhl(long double);
+long double 	coshl(long double);
+long double 	exp2l(long double);
+long double 	expl(long double);
+long double 	fmodl(long double, long double);
+long double 	hypotl(long double, long double);
+long double 	lgammal(long double);
+long double 	log10l(long double);
+long double 	log2l(long double);
+long double 	logl(long double);
+long double 	powl(long double, long double);
+long double 	remainderl(long double, long double);
+long double 	sinhl(long double);
+long double 	sqrtl(long double);
+long double		j0l(long double x);
+long double		j1l(long double x);
+long double		jnl(int n, long double x);
+long double		y0l(long double x);
+long double		y1l(long double x);
+long double		ynl(int n, long double x);
+long double 	tgammal(long double x);
+long double 	scalbl(long double x, long double fn);
+long double 	gammal(long double x);
+long double 	scalbl(long double x, long double fn);
 
-#ifdef L_acosl
-WRAPPER1(acos)
 #endif
 
 #ifdef L_asinhl
 WRAPPER1(asinh)
 #endif
 
-#ifdef L_asinl
-WRAPPER1(asin)
-#endif
-
-#ifdef L_atan2l
-WRAPPER2(atan2)
-#endif
-
-#ifdef L_atanhl
-WRAPPER1(atanh)
-#endif
-
 #ifdef L_atanl
 WRAPPER1(atan)
 #endif
@@ -159,10 +146,6 @@ WRAPPER1(ceil)
 WRAPPER2(copysign)
 #endif
 
-#ifdef L_coshl
-WRAPPER1(cosh)
-#endif
-
 #ifdef L_cosl
 WRAPPER1(cos)
 libm_hidden_def(cosl)
@@ -176,15 +159,6 @@ WRAPPER1(erfc)
 WRAPPER1(erf)
 #endif
 
-#ifdef L_exp2l
-WRAPPER1(exp2)
-#endif
-
-#ifdef L_expl
-WRAPPER1(exp)
-libm_hidden_def(expl)
-#endif
-
 #ifdef L_expm1l
 WRAPPER1(expm1)
 #endif
@@ -216,10 +190,6 @@ WRAPPER2(fmax)
 WRAPPER2(fmin)
 #endif
 
-#ifdef L_fmodl
-WRAPPER2(fmod)
-#endif
-
 #ifdef L_frexpl
 long double frexpl (long double x, int *ex)
 {
@@ -227,34 +197,10 @@ long double frexpl (long double x, int *ex)
 }
 #endif
 
-#ifdef L_gammal
-WRAPPER1(gamma)
-#endif
-
-#ifdef L_hypotl
-WRAPPER2(hypot)
-libm_hidden_def(hypotl)
-#endif
-
 #ifdef L_ilogbl
 int_WRAPPER1(ilogb)
 #endif
 
-#ifdef L_j0l
-	WRAPPER1(j0)
-#endif
-
-#ifdef L_j1l
-	WRAPPER1(j1)
-#endif
-
-#ifdef L_jnl
-long double jnl(int n, long double x)
-{
-	return (long double) jn(n, (double)x);
-}
-#endif
-
 #ifdef L_ldexpl
 long double ldexpl (long double x, int ex)
 {
@@ -262,10 +208,6 @@ long double ldexpl (long double x, int ex)
 }
 #endif
 
-#ifdef L_lgammal
-WRAPPER1(lgamma)
-#endif
-
 #ifdef L_llrintl
 long_long_WRAPPER1(llrint)
 #endif
@@ -274,26 +216,14 @@ long_long_WRAPPER1(llrint)
 long_long_WRAPPER1(llround)
 #endif
 
-#ifdef L_log10l
-WRAPPER1(log10)
-#endif
-
 #ifdef L_log1pl
 WRAPPER1(log1p)
 #endif
 
-#ifdef L_log2l
-WRAPPER1(log2)
-#endif
-
 #ifdef L_logbl
 WRAPPER1(logb)
 #endif
 
-#ifdef L_logl
-WRAPPER1(log)
-#endif
-
 #ifdef L_lrintl
 long_WRAPPER1(lrint)
 #endif
@@ -332,14 +262,6 @@ long double nexttowardl(long double x, long double y)
 #endif
 #endif
 
-#ifdef L_powl
-WRAPPER2(pow)
-#endif
-
-#ifdef L_remainderl
-WRAPPER2(remainder)
-#endif
-
 #ifdef L_remquol
 long double remquol (long double x, long double y, int *quo)
 {
@@ -369,21 +291,11 @@ long double scalbnl (long double x, int ex)
 }
 #endif
 
-/* scalb is an obsolete function */
-
-#ifdef L_sinhl
-WRAPPER1(sinh)
-#endif
-
 #ifdef L_sinl
 WRAPPER1(sin)
 libm_hidden_def(sinl)
 #endif
 
-#ifdef L_sqrtl
-WRAPPER1(sqrt)
-#endif
-
 #ifdef L_tanhl
 WRAPPER1(tanh)
 #endif
@@ -392,10 +304,6 @@ WRAPPER1(tanh)
 WRAPPER1(tan)
 #endif
 
-#ifdef L_tgammal
-WRAPPER1(tgamma)
-#endif
-
 #ifdef L_truncl
 WRAPPER1(trunc)
 #endif
@@ -404,21 +312,6 @@ WRAPPER1(trunc)
 WRAPPER1(significand)
 #endif
 
-#ifdef L_y0l
-WRAPPER1(y0)
-#endif
-
-#ifdef L_y1l
-WRAPPER1(y1)
-#endif
-
-#ifdef L_ynl
-long double ynl(int n, long double x)
-{
-	return (long double) yn(n, (double)x);
-}
-#endif
-
 
 #if defined __DO_C99_MATH__ && !defined __NO_LONG_DOUBLE_MATH
 

+ 10 - 0
libm/math_private.h

@@ -184,6 +184,16 @@ extern double __kernel_sin (double,double,int) attribute_hidden;
 extern double __kernel_cos (double,double) attribute_hidden;
 extern double __kernel_tan (double,double,int) attribute_hidden;
 extern int    __kernel_rem_pio2 (double*,double*,int,int,int,const int*) attribute_hidden;
+extern double __kernel_standard(double x, double y, int type) attribute_hidden;
+extern float  __kernel_standard_f (float, float, int) attribute_hidden;
+#ifndef __NO_LONG_DOUBLE_MATH
+extern long double __kernel_standard_l (long double, long double, int) attribute_hidden;
+#endif
+/* wrappers functions for internal use */
+extern float __lgammaf_r (float, int*);
+extern double __lgamma_r (double, int*);
+extern long double __lgammal_r(long double, int*);
+extern double __ieee754_tgamma(double);
 
 /*
  * math_opt_barrier(x): safely load x, even if it was manipulated

+ 40 - 0
libm/w_acos.c

@@ -0,0 +1,40 @@
+/* Copyright (C) 2011-2016 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
+
+#include <math.h>
+#include "math_private.h"
+#if __UCLIBC_HAS_FENV__
+#include <fenv.h>
+#endif
+
+/* wrapper acos */
+double
+acos (double x)
+{
+#if __UCLIBC_HAS_FENV__
+  if (__builtin_expect (isgreater (fabs (x), 1.0), 0)
+      && _LIB_VERSION != _IEEE_)
+    {
+      /* acos(|x|>1) */
+      feraiseexcept (FE_INVALID);
+      return __kernel_standard (x, x, 1);
+    }
+#endif
+  return __ieee754_acos (x);
+}
+libm_hidden_def(acos)

+ 38 - 0
libm/w_acosf.c

@@ -0,0 +1,38 @@
+/* Copyright (C) 2011-2016 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
+
+#include <math.h>
+#include "math_private.h"
+#if __UCLIBC_HAS_FENV__
+#include <fenv.h>
+#endif
+/* wrapper acosf */
+float
+acosf (float x)
+{
+#if __UCLIBC_HAS_FENV__
+  if (__builtin_expect (isgreater (fabsf (x), 1.0f), 0)
+      && _LIB_VERSION != _IEEE_)
+    {
+      /* acos(|x|>1) */
+      feraiseexcept (FE_INVALID);
+      return __kernel_standard_f (x, x, 101);
+    }
+#endif
+  return (float) __ieee754_acos ((float) x);
+}

+ 34 - 0
libm/w_acosh.c

@@ -0,0 +1,34 @@
+/* Copyright (C) 2011-2016 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
+
+#include <math.h>
+#include "math_private.h"
+
+
+/* wrapper acosh */
+double
+acosh (double x)
+{
+#if __UCLIBC_HAS_FENV__
+  if (__builtin_expect (isless (x,  1.0), 0) && _LIB_VERSION != _IEEE_)
+    /* acosh(x<1) */
+    return __kernel_standard (x, x, 29);
+#endif
+  return __ieee754_acosh (x);
+}
+libm_hidden_def(acosh)

+ 33 - 0
libm/w_acoshf.c

@@ -0,0 +1,33 @@
+/* Copyright (C) 2011-2016 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
+
+#include <math.h>
+#include "math_private.h"
+
+
+/* wrapper acoshf */
+float
+acoshf (float x)
+{
+#if __UCLIBC_HAS_FENV__
+  if (__builtin_expect (isless (x, 1.0f), 0) && _LIB_VERSION != _IEEE_)
+    /* acosh(x<1) */
+    return __kernel_standard_f (x, x, 129);
+#endif
+    return (float) __ieee754_acosh ((double) x);
+}

+ 34 - 0
libm/w_acoshl.c

@@ -0,0 +1,34 @@
+/* Copyright (C) 2011-2016 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
+
+#include <math.h>
+#include "math_private.h"
+
+#if !defined __NO_LONG_DOUBLE_MATH
+/* wrapper acosl */
+long double
+acoshl (long double x)
+{
+#if __UCLIBC_HAS_FENV__
+  if (__builtin_expect (isless (x, 1.0L), 0) && _LIB_VERSION != _IEEE_)
+    /* acosh(x<1) */
+    return __kernel_standard_l (x, x, 229);
+#endif
+  return (long double) __ieee754_acosh ((double) x);
+}
+#endif

+ 42 - 0
libm/w_acosl.c

@@ -0,0 +1,42 @@
+/* Copyright (C) 2011-2016 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
+
+
+#include <math.h>
+#include "math_private.h"
+#if __UCLIBC_HAS_FENV__
+#include <fenv.h>
+#endif
+
+#if !defined __NO_LONG_DOUBLE_MATH
+/* wrapper acosl */
+long double
+acosl (long double x)
+{
+#if __UCLIBC_HAS_FENV__
+  if (__builtin_expect (isgreater (fabsl (x), 1.0L), 0)
+      && _LIB_VERSION != _IEEE_)
+    {
+      /* acos(|x|>1) */
+      feraiseexcept (FE_INVALID);
+      return __kernel_standard_l (x, x, 201);
+    }
+#endif
+  return (long double) __ieee754_acos ((double) x);
+}
+#endif

+ 39 - 0
libm/w_asin.c

@@ -0,0 +1,39 @@
+/* Copyright (C) 2011-2016 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
+
+#include <math.h>
+#include "math_private.h"
+#if __UCLIBC_HAS_FENV__
+#include <fenv.h>
+#endif
+/* wrapper asin */
+double
+asin (double x)
+{
+#if __UCLIBC_HAS_FENV__
+  if (__builtin_expect (isgreater (fabs (x),  1.0), 0)
+      && _LIB_VERSION != _IEEE_)
+    {
+      /* asin(|x|>1) */
+      feraiseexcept (FE_INVALID);
+      return __kernel_standard (x, x, 2);
+    }
+#endif
+  return __ieee754_asin (x);
+}
+libm_hidden_def(asin)

+ 38 - 0
libm/w_asinf.c

@@ -0,0 +1,38 @@
+/* Copyright (C) 2011-2016 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
+
+#include <math.h>
+#include "math_private.h"
+#if __UCLIBC_HAS_FENV__
+#include <fenv.h>
+#endif
+/* wrapper asinf */
+float
+asinf (float x)
+{
+#if __UCLIBC_HAS_FENV__
+  if (__builtin_expect (isgreater (fabsf (x), 1.0f), 0)
+      && _LIB_VERSION != _IEEE_)
+    {
+      /* asin(|x|>1) */
+      feraiseexcept (FE_INVALID);
+      return __kernel_standard_f (x, x, 102);
+    }
+#endif
+  return (float) __ieee754_asin ((double)x);
+}

+ 41 - 0
libm/w_asinl.c

@@ -0,0 +1,41 @@
+/* Copyright (C) 2011-2016 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
+
+#include <math.h>
+#include "math_private.h"
+#if __UCLIBC_HAS_FENV__
+#include <fenv.h>
+#endif
+
+#if !defined __NO_LONG_DOUBLE_MATH
+/* wrapper asinl */
+long double
+asinl (long double x)
+{
+# if __UCLIBC_HAS_FENV__
+  if (__builtin_expect (isgreater (fabsl (x), 1.0L), 0)
+      && _LIB_VERSION != _IEEE_)
+    {
+      /* asin(|x|>1) */
+      feraiseexcept (FE_INVALID);
+      return __kernel_standard_l (x, x, 202);
+    }
+# endif
+  return (long double)__ieee754_asin ((double)x);
+}
+#endif

+ 46 - 0
libm/w_atan2.c

@@ -0,0 +1,46 @@
+/* Copyright (C) 2011-2016 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
+
+/*
+ * wrapper atan2(y,x)
+ */
+
+#include <math.h>
+#include "math_private.h"
+#if __UCLIBC_HAS_FENV__
+#include <errno.h>
+#endif
+
+double
+atan2 (double y, double x)
+{
+#if __UCLIBC_HAS_FENV__
+  double z;
+
+  if (__builtin_expect (x == 0.0 && y == 0.0, 0) && _LIB_VERSION == _SVID_)
+    return __kernel_standard (y, x, 3); /* atan2(+-0,+-0) */
+
+  z = __ieee754_atan2 (y, x);
+  if (__builtin_expect (z == 0.0 && y != 0.0 && isfinite (x),0))
+    __set_errno (ERANGE);
+  return z;
+#else
+  return __ieee754_atan2 (y, x);
+#endif
+}
+libm_hidden_def(atan2)

+ 46 - 0
libm/w_atan2f.c

@@ -0,0 +1,46 @@
+/* Copyright (C) 2011-2016 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
+
+/*
+ * wrapper atan2f(y,x)
+ */
+
+
+#include <math.h>
+#include "math_private.h"
+#if __UCLIBC_HAS_FENV__
+#include <errno.h>
+#endif
+
+float
+atan2f (float y, float x)
+{
+#if __UCLIBC_HAS_FENV__
+  float z;
+
+  if (__builtin_expect (x == 0.0f && y == 0.0f, 0) && _LIB_VERSION == _SVID_)
+    return __kernel_standard_f (y, x, 103); /* atan2(+-0,+-0) */
+
+  z = (float) __ieee754_atan2 ((double) y, (double) x);
+  if (__builtin_expect (z == 0.0f && y != 0.0f && isfinite (x), 0))
+    __set_errno (ERANGE);
+  return z;
+#else
+  return (float) __ieee754_atan2 ((double) y, (double) x);
+#endif
+}

+ 44 - 0
libm/w_atan2l.c

@@ -0,0 +1,44 @@
+/* Copyright (C) 2011-2016 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
+
+#include <math.h>
+#include "math_private.h"
+#if __UCLIBC_HAS_FENV__
+#include <errno.h>
+#endif
+
+#if !defined __NO_LONG_DOUBLE_MATH
+/* wrapper atan2l(y,x) */
+long double
+atan2l (long double y, long double x)
+{
+# if __UCLIBC_HAS_FENV__
+  long double z;
+
+  if (__builtin_expect (x == 0.0L && y == 0.0L, 0) && _LIB_VERSION == _SVID_)
+    return __kernel_standard_l (y, x, 203); /* atan2(+-0,+-0) */
+
+  z = (long double) __ieee754_atan2 ((double)y,(double) x);
+  if (__builtin_expect (z == 0.0L && y != 0.0L && isfinite (x),0))
+    __set_errno (ERANGE);
+  return z;
+# else
+  return (long double) __ieee754_atan2 ((double)y,(double) x);
+# endif /* __UCLIBC_HAS_FENV__ */
+}
+#endif /* __NO_LONG_DOUBLE_MATH */

+ 37 - 0
libm/w_atanh.c

@@ -0,0 +1,37 @@
+/* Copyright (C) 2011-2016 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
+
+#include <math.h>
+#include "math_private.h"
+
+
+/* wrapper atanh */
+double
+atanh (double x)
+{
+#if __UCLIBC_HAS_FENV__
+  if (__builtin_expect (isgreaterequal (fabs (x), 1.0), 0)
+      && _LIB_VERSION != _IEEE_)
+    return __kernel_standard (x, x,
+			      fabs (x) > 1.0
+			      ? 30		/* atanh(|x|>1) */
+			      : 31);		/* atanh(|x|==1) */
+#endif
+  return __ieee754_atanh (x);
+}
+libm_hidden_def(atanh)

+ 36 - 0
libm/w_atanhf.c

@@ -0,0 +1,36 @@
+/* Copyright (C) 2011-2016 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
+
+#include <math.h>
+#include "math_private.h"
+
+
+/* wrapper atanhf */
+float
+atanhf (float x)
+{
+#if __UCLIBC_HAS_FENV__
+  if (__builtin_expect (isgreaterequal (fabsf (x), 1.0f), 0)
+      && _LIB_VERSION != _IEEE_)
+    return __kernel_standard_f (x, x,
+				fabsf (x) > 1.0f
+				? 130		/* atanh(|x|>1) */
+				: 131);		/* atanh(|x|==1) */
+#endif
+  return (float) __ieee754_atanh ((double) x);
+}

+ 38 - 0
libm/w_atanhl.c

@@ -0,0 +1,38 @@
+/* Copyright (C) 2011-2016 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
+
+#include <math.h>
+#include "math_private.h"
+