Attachment #243657 for bug #272765

View | Details | Raw Unified | Return to bug 272765
Collapse All | Expand All




#endif
#endif

/* Write a pari script to test things externally. */
#ifdef DOPRINT
#include <stdio.h>

#ifndef DOPRINT_SWIZZLE
#define	DOPRINT_SWIZZLE		0
#endif

#ifdef DOPRINT_LD80

#define	DOPRINT_START(xp) do {						\
	uint64_t __lx;							\
	uint16_t __hx;							\
									\
	/* Hack to give more-problematic args. */			\
	EXTRACT_LDBL80_WORDS(__hx, __lx, *xp);				\
	__lx ^= DOPRINT_SWIZZLE;					\
	INSERT_LDBL80_WORDS(*xp, __hx, __lx);				\
	printf("x = %.21Lg; ", (long double)*xp);			\
} while (0)
#define	DOPRINT_END1(v)							\
	printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
#define	DOPRINT_END2(hi, lo)						\
	printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n",		\
	    (long double)(hi), (long double)(lo))

#elif defined(DOPRINT_D64)

#define	DOPRINT_START(xp) do {						\
	uint32_t __hx, __lx;						\
									\
	EXTRACT_WORDS(__hx, __lx, *xp);					\
	__lx ^= DOPRINT_SWIZZLE;					\
	INSERT_WORDS(*xp, __hx, __lx);					\
	printf("x = %.21Lg; ", (long double)*xp);			\
} while (0)
#define	DOPRINT_END1(v)							\
	printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
#define	DOPRINT_END2(hi, lo)						\
	printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n",		\
	    (long double)(hi), (long double)(lo))

#elif defined(DOPRINT_F32)

#define	DOPRINT_START(xp) do {						\
	uint32_t __hx;							\
									\
	GET_FLOAT_WORD(__hx, *xp);					\
	__hx ^= DOPRINT_SWIZZLE;					\
	SET_FLOAT_WORD(*xp, __hx);					\
	printf("x = %.21Lg; ", (long double)*xp);			\
} while (0)
#define	DOPRINT_END1(v)							\
	printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
#define	DOPRINT_END2(hi, lo)						\
	printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n",		\
	    (long double)(hi), (long double)(lo))

#else /* !DOPRINT_LD80 && !DOPRINT_D64 (LD128 only) */

#ifndef DOPRINT_SWIZZLE_HIGH
#define	DOPRINT_SWIZZLE_HIGH	0
#endif

#define	DOPRINT_START(xp) do {						\
	uint64_t __lx, __llx;						\
	uint16_t __hx;							\
									\
	EXTRACT_LDBL128_WORDS(__hx, __lx, __llx, *xp);			\
	__llx ^= DOPRINT_SWIZZLE;					\
	__lx ^= DOPRINT_SWIZZLE_HIGH;					\
	INSERT_LDBL128_WORDS(*xp, __hx, __lx, __llx);			\
	printf("x = %.36Lg; ", (long double)*xp);					\
} while (0)
#define	DOPRINT_END1(v)							\
	printf("y = %.36Lg; z = 0; show(x, y, z);\n", (long double)(v))
#define	DOPRINT_END2(hi, lo)						\
	printf("y = %.36Lg; z = %.36Lg; show(x, y, z);\n",		\
	    (long double)(hi), (long double)(lo))

#endif /* DOPRINT_LD80 */

#else /* !DOPRINT */
#define	DOPRINT_START(xp)
#define	DOPRINT_END1(v)
#define	DOPRINT_END2(hi, lo)
#endif /* DOPRINT */

#define	RETURNP(x) do {			\
	DOPRINT_END1(x);		\
	RETURNF(x);			\
} while (0)
#define	RETURNPI(x) do {		\
	DOPRINT_END1(x);		\
	RETURNI(x);			\
} while (0)
#define	RETURN2P(x, y) do {		\
	DOPRINT_END2((x), (y));		\
	RETURNF((x) + (y));		\
} while (0)
#define	RETURN2PI(x, y) do {		\
	DOPRINT_END2((x), (y));		\
	RETURNI((x) + (y));		\
} while (0)
#ifdef STRUCT_RETURN
#define	RETURNSP(rp) do {		\
	if (!(rp)->lo_set)		\
		RETURNF((rp)->hi);	\
	RETURNF((rp)->hi + (rp)->lo);	\
} while (0)
#define	RETURNSPI(rp) do {		\
	if (!(rp)->lo_set)		\
		RETURNI((rp)->hi);	\
	RETURNI((rp)->hi + (rp)->lo);	\
} while (0)
#endif

#define	SUM2P(x, y) ({			\
	const __typeof (x) __x = (x);	\
	const __typeof (y) __y = (y);	\
					\
	DOPRINT_END2(__x, __y);		\
	__x + __y;			\
})





	int k;
	uint16_t hx, ix;

	DOPRINT_START(&x);

	/* Filter out exceptional cases. */
	u.e = x;
	hx = u.xbits.expsign;

	if (ix >= BIAS + 13) {		/* |x| >= 8192 or x is NaN */
		if (ix == BIAS + LDBL_MAX_EXP) {
			if (hx & 0x8000)  /* x is -Inf, -NaN or unsupported */
				RETURNF(-1 / x);
			RETURNF(x + x);	/* x is +Inf, +NaN or unsupported */
		}
		if (x > o_threshold)
			RETURNF(huge * huge);
		if (x < u_threshold)
			RETURNF(tiny * tiny);
	} else if (ix < BIAS - 75) {	/* |x| < 0x1p-75 (includes pseudos) */
		RETURNF(1 + x);		/* 1 with inexact iff x != 0 */
	}

	ENTERI();

	int k, n, n2;
	uint16_t hx, ix;

	DOPRINT_START(&x);

	/* Filter out exceptional cases. */
	u.e = x;
	hx = u.xbits.expsign;

	if (ix >= BIAS + 6) {		/* |x| >= 64 or x is NaN */
		if (ix == BIAS + LDBL_MAX_EXP) {
			if (hx & 0x8000)  /* x is -Inf, -NaN or unsupported */
				RETURNF(-1 / x - 1);
			RETURNF(x + x);	/* x is +Inf, +NaN or unsupported */
		}
		if (x > o_threshold)
			RETURNF(huge * huge);
		/*
		 * expm1l() never underflows, but it must avoid
		 * unrepresentable large negative exponents.  We used a

		 * in the same way as large ones here.
		 */
		if (hx & 0x8000)	/* x <= -64 */
			RETURNF(tiny - 1);	/* good for x < -65ln2 - eps */
	}

	ENTERI();

	if (T1 < x && x < T2) {
		if (ix < BIAS - 74) {	/* |x| < 0x1p-74 (includes pseudos) */
			/* x (rounded) with inexact if x != 0: */
			RETURNI(x == 0 ? x :
			    (0x1p100 * x + fabsl(x)) * 0x1p-100);
		}


		hx2_hi = x_hi * x_hi / 2;
		hx2_lo = x_lo * (x + x_hi) / 2;
		if (ix >= BIAS - 7)
			RETURNI((hx2_hi + x_hi) + (hx2_lo + x_lo + q));
		else
			RETURNI(x + (hx2_lo + q + hx2_hi));
	}

	/* Reduce x to (k*ln2 + endpoint[n2] + r1 + r2). */




	int i, k;
	int16_t ax, hx;

	DOPRINT_START(&x);
	EXTRACT_LDBL80_WORDS(hx, lx, x);
	if (hx < 0x3fff) {		/* x < 1, or x neg NaN */
		ax = hx & 0x7fff;
		if (ax >= 0x3fff) {	/* x <= -1, or x neg NaN */
			if (ax == 0x3fff && lx == 0x8000000000000000ULL)
				RETURNF(-1 / zero);	/* log1p(-1) = -Inf */
			/* log1p(x < 1, or x [pseudo-]NaN) = qNaN: */
			RETURNF((x - x) / (x - x));
		}
		if (ax <= 0x3fbe) {	/* |x| < 2**-64 */
			if ((int)x == 0)
				RETURNF(x);	/* x with inexact if x != 0 */
		}
		f_hi = 1;
		f_lo = x;
	} else if (hx >= 0x7fff) {	/* x +Inf or non-neg NaN */
		RETURNF(x + x);		/* log1p(Inf or NaN) = Inf or qNaN */
					/* log1p(pseudo-Inf) = qNaN */
					/* log1p(pseudo-NaN) = qNaN */
					/* log1p(unnormal) = qNaN */

#endif

	_3sumF(val_hi, val_lo, F_hi(i) + dk * ln2_hi);
	RETURNI(val_hi + val_lo);
}

#ifdef STRUCT_RETURN

	struct ld r;

	ENTERI();
	DOPRINT_START(&x);
	k_logl(x, &r);
	RETURNSPI(&r);
}

	long double hi, lo;

	ENTERI();
	DOPRINT_START(&x);
	k_logl(x, &r);
	if (!r.lo_set)
		RETURNI(r.hi);
	_2sumF(r.hi, r.lo);
	hi = (float)r.hi;
	lo = r.lo + (r.hi - hi);
	RETURNI(invln10_hi * hi + 
	    (invln10_lo_plus_hi * lo + invln10_lo * hi));
}

long double

	long double hi, lo;

	ENTERI();
	DOPRINT_START(&x);
	k_logl(x, &r);
	if (!r.lo_set)
		RETURNI(r.hi);
	_2sumF(r.hi, r.lo);
	hi = (float)r.hi;
	lo = r.lo + (r.hi - hi);
	RETURNI(invln2_hi * hi +
	    (invln2_lo_plus_hi * lo + invln2_lo * hi));
}

#endif /* STRUCT_RETURN */




	int k;
	uint16_t hx, ix;

	DOPRINT_START(&x);

	/* Filter out exceptional cases. */
	u.e = x;
	hx = u.xbits.expsign;

	if (ix >= BIAS + 13) {		/* |x| >= 8192 or x is NaN */
		if (ix == BIAS + LDBL_MAX_EXP) {
			if (hx & 0x8000)  /* x is -Inf or -NaN */
				RETURNF(-1 / x);
			RETURNF(x + x);	/* x is +Inf or +NaN */
		}
		if (x > o_threshold)
			RETURNF(huge * huge);
		if (x < u_threshold)
			RETURNF(tiny * tiny);
	} else if (ix < BIAS - 114) {	/* |x| < 0x1p-114 */
		RETURNF(1 + x);		/* 1 with inexact iff x != 0 */
	}

	ENTERI();

	int k, n, n2;
	uint16_t hx, ix;

	DOPRINT_START(&x);

	/* Filter out exceptional cases. */
	u.e = x;
	hx = u.xbits.expsign;

	if (ix >= BIAS + 7) {		/* |x| >= 128 or x is NaN */
		if (ix == BIAS + LDBL_MAX_EXP) {
			if (hx & 0x8000)  /* x is -Inf or -NaN */
				RETURNF(-1 / x - 1);
			RETURNF(x + x);	/* x is +Inf or +NaN */
		}
		if (x > o_threshold)
			RETURNF(huge * huge);
		/*
		 * expm1l() never underflows, but it must avoid
		 * unrepresentable large negative exponents.  We used a

		 * in the same way as large ones here.
		 */
		if (hx & 0x8000)	/* x <= -128 */
			RETURNF(tiny - 1);	/* good for x < -114ln2 - eps */
	}

	ENTERI();

		if (x < T3) {
			if (ix < BIAS - 113) {	/* |x| < 0x1p-113 */
				/* x (rounded) with inexact if x != 0: */
				RETURNI(x == 0 ? x :
				    (0x1p200 * x + fabsl(x)) * 0x1p-200);
			}
			q = x * x2 * C3 + x2 * x2 * (C4 + x * (C5 + x * (C6 +

		hx2_hi = x_hi * x_hi / 2;
		hx2_lo = x_lo * (x + x_hi) / 2;
		if (ix >= BIAS - 7)
			RETURNI((hx2_hi + x_hi) + (hx2_lo + x_lo + q));
		else
			RETURNI(x + (hx2_lo + q + hx2_hi));
	}

	/* Reduce x to (k*ln2 + endpoint[n2] + r1 + r2). */




	int i, k;
	int16_t ax, hx;

	DOPRINT_START(&x);
	EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
	if (hx < 0x3fff) {		/* x < 1, or x neg NaN */
		ax = hx & 0x7fff;
		if (ax >= 0x3fff) {	/* x <= -1, or x neg NaN */
			if (ax == 0x3fff && (lx | llx) == 0)
				RETURNF(-1 / zero);	/* log1p(-1) = -Inf */
			/* log1p(x < 1, or x NaN) = qNaN: */
			RETURNF((x - x) / (x - x));
		}
		if (ax <= 0x3f8d) {	/* |x| < 2**-113 */
			if ((int)x == 0)
				RETURNF(x);	/* x with inexact if x != 0 */
		}
		f_hi = 1;
		f_lo = x;
	} else if (hx >= 0x7fff) {	/* x +Inf or non-neg NaN */
		RETURNF(x + x);		/* log1p(Inf or NaN) = Inf or qNaN */
	} else if (hx < 0x40e1) {	/* 1 <= x < 2**226 */
		f_hi = x;
		f_lo = 1;

#endif

	_3sumF(val_hi, val_lo, F_hi(i) + dk * ln2_hi);
	RETURNI(val_hi + val_lo);
}

#ifdef STRUCT_RETURN

	struct ld r;

	ENTERI();
	DOPRINT_START(&x);
	k_logl(x, &r);
	RETURNSPI(&r);
}

	long double hi, lo;

	ENTERI();
	DOPRINT_START(&x);
	k_logl(x, &r);
	if (!r.lo_set)
		RETURNI(r.hi);
	_2sumF(r.hi, r.lo);
	hi = (float)r.hi;
	lo = r.lo + (r.hi - hi);
	RETURNI(invln10_hi * hi + (invln10_lo_plus_hi * lo + invln10_lo * hi));
	    invln10_lo_plus_hi * lo + invln10_lo * hi);
}

long double

	long double hi, lo;

	ENTERI();
	DOPRINT_START(&x);
	k_logl(x, &r);
	if (!r.lo_set)
		RETURNI(r.hi);
	_2sumF(r.hi, r.lo);
	hi = (float)r.hi;
	lo = r.lo + (r.hi - hi);
	RETURNI(invln2_hi * hi + (invln2_lo_plus_hi * lo + invln2_lo * hi));
	    invln2_lo_plus_hi * lo + invln2_lo * hi);
}

#endif /* STRUCT_RETURN */

Return to bug 272765

Lines 698-824 irintl(long double x) Link Here

(-)b/lib/msun/src/math_private.h (-110 / +5 lines)
698	#endif	698	#endif
699	#endif	699	#endif
700		700
701	/* Write a pari script to test things externally. */
702	#ifdef DOPRINT
703	#include <stdio.h>
704
705	#ifndef DOPRINT_SWIZZLE
706	#define DOPRINT_SWIZZLE 0
707	#endif
708
709	#ifdef DOPRINT_LD80
710
711	#define DOPRINT_START(xp) do { \
712	uint64_t __lx; \
713	uint16_t __hx; \
714	\
715	/* Hack to give more-problematic args. */ \
716	EXTRACT_LDBL80_WORDS(__hx, __lx, *xp); \
717	__lx ^= DOPRINT_SWIZZLE; \
718	INSERT_LDBL80_WORDS(*xp, __hx, __lx); \
719	printf("x = %.21Lg; ", (long double)*xp); \
720	} while (0)
721	#define DOPRINT_END1(v) \
722	printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
723	#define DOPRINT_END2(hi, lo) \
724	printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \
725	(long double)(hi), (long double)(lo))
726
727	#elif defined(DOPRINT_D64)
728
729	#define DOPRINT_START(xp) do { \
730	uint32_t __hx, __lx; \
731	\
732	EXTRACT_WORDS(__hx, __lx, *xp); \
733	__lx ^= DOPRINT_SWIZZLE; \
734	INSERT_WORDS(*xp, __hx, __lx); \
735	printf("x = %.21Lg; ", (long double)*xp); \
736	} while (0)
737	#define DOPRINT_END1(v) \
738	printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
739	#define DOPRINT_END2(hi, lo) \
740	printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \
741	(long double)(hi), (long double)(lo))
742
743	#elif defined(DOPRINT_F32)
744
745	#define DOPRINT_START(xp) do { \
746	uint32_t __hx; \
747	\
748	GET_FLOAT_WORD(__hx, *xp); \
749	__hx ^= DOPRINT_SWIZZLE; \
750	SET_FLOAT_WORD(*xp, __hx); \
751	printf("x = %.21Lg; ", (long double)*xp); \
752	} while (0)
753	#define DOPRINT_END1(v) \
754	printf("y = %.21Lg; z = 0; show(x, y, z);\n", (long double)(v))
755	#define DOPRINT_END2(hi, lo) \
756	printf("y = %.21Lg; z = %.21Lg; show(x, y, z);\n", \
757	(long double)(hi), (long double)(lo))
758
759	#else /* !DOPRINT_LD80 && !DOPRINT_D64 (LD128 only) */
760
761	#ifndef DOPRINT_SWIZZLE_HIGH
762	#define DOPRINT_SWIZZLE_HIGH 0
763	#endif
764
765	#define DOPRINT_START(xp) do { \
766	uint64_t __lx, __llx; \
767	uint16_t __hx; \
768	\
769	EXTRACT_LDBL128_WORDS(__hx, __lx, __llx, *xp); \
770	__llx ^= DOPRINT_SWIZZLE; \
771	__lx ^= DOPRINT_SWIZZLE_HIGH; \
772	INSERT_LDBL128_WORDS(*xp, __hx, __lx, __llx); \
773	printf("x = %.36Lg; ", (long double)*xp); \
774	} while (0)
775	#define DOPRINT_END1(v) \
776	printf("y = %.36Lg; z = 0; show(x, y, z);\n", (long double)(v))
777	#define DOPRINT_END2(hi, lo) \
778	printf("y = %.36Lg; z = %.36Lg; show(x, y, z);\n", \
779	(long double)(hi), (long double)(lo))
780
781	#endif /* DOPRINT_LD80 */
782
783	#else /* !DOPRINT */
784	#define DOPRINT_START(xp)
785	#define DOPRINT_END1(v)
786	#define DOPRINT_END2(hi, lo)
787	#endif /* DOPRINT */
788
789	#define RETURNP(x) do { \
790	DOPRINT_END1(x); \
791	RETURNF(x); \
792	} while (0)
793	#define RETURNPI(x) do { \
794	DOPRINT_END1(x); \
795	RETURNI(x); \
796	} while (0)
797	#define RETURN2P(x, y) do { \
798	DOPRINT_END2((x), (y)); \
799	RETURNF((x) + (y)); \
800	} while (0)
801	#define RETURN2PI(x, y) do { \
802	DOPRINT_END2((x), (y)); \
803	RETURNI((x) + (y)); \
804	} while (0)
805	#ifdef STRUCT_RETURN	701	#ifdef STRUCT_RETURN
806	#define RETURNSP(rp) do { \	702	#define RETURNSP(rp) do { \
807	if (!(rp)->lo_set) \	703	if (!(rp)->lo_set) \
808	RETURNP((rp)->hi); \	704	RETURNF((rp)->hi); \
809	RETURN2P((rp)->hi, (rp)->lo); \	705	RETURNF((rp)->hi + (rp)->lo); \
810	} while (0)	706	} while (0)
811	#define RETURNSPI(rp) do { \	707	#define RETURNSPI(rp) do { \
812	if (!(rp)->lo_set) \	708	if (!(rp)->lo_set) \
813	RETURNPI((rp)->hi); \	709	RETURNI((rp)->hi); \
814	RETURN2PI((rp)->hi, (rp)->lo); \	710	RETURNI((rp)->hi + (rp)->lo); \
815	} while (0)	711	} while (0)
816	#endif	712	#endif
		713
817	#define SUM2P(x, y) ({ \	714	#define SUM2P(x, y) ({ \
818	const __typeof (x) __x = (x); \	715	const __typeof (x) __x = (x); \
819	const __typeof (y) __y = (y); \	716	const __typeof (y) __y = (y); \
820	\
821	DOPRINT_END2(__x, __y); \
822	__x + __y; \	717	__x + __y; \
823	})	718	})
824		719

Lines 77-84 expl(long double x) Link Here

(-)b/lib/msun/ld80/s_expl.c (-16 / +12 lines)
77	int k;	77	int k;
78	uint16_t hx, ix;	78	uint16_t hx, ix;
79		79
80	DOPRINT_START(&x);
81
82	/* Filter out exceptional cases. */	80	/* Filter out exceptional cases. */
83	u.e = x;	81	u.e = x;
84	hx = u.xbits.expsign;	82	hx = u.xbits.expsign;
Lines 86-100 expl(long double x) Link Here
86	if (ix >= BIAS + 13) { /* \|x\| >= 8192 or x is NaN */	84	if (ix >= BIAS + 13) { /* \|x\| >= 8192 or x is NaN */
87	if (ix == BIAS + LDBL_MAX_EXP) {	85	if (ix == BIAS + LDBL_MAX_EXP) {
88	if (hx & 0x8000) /* x is -Inf, -NaN or unsupported */	86	if (hx & 0x8000) /* x is -Inf, -NaN or unsupported */
89	RETURNP(-1 / x);	87	RETURNF(-1 / x);
90	RETURNP(x + x); /* x is +Inf, +NaN or unsupported */	88	RETURNF(x + x); /* x is +Inf, +NaN or unsupported */
91	}	89	}
92	if (x > o_threshold)	90	if (x > o_threshold)
93	RETURNP(huge * huge);	91	RETURNF(huge * huge);
94	if (x < u_threshold)	92	if (x < u_threshold)
95	RETURNP(tiny * tiny);	93	RETURNF(tiny * tiny);
96	} else if (ix < BIAS - 75) { /* \|x\| < 0x1p-75 (includes pseudos) */	94	} else if (ix < BIAS - 75) { /* \|x\| < 0x1p-75 (includes pseudos) */
97	RETURN2P(1, x); /* 1 with inexact iff x != 0 */	95	RETURNF(1 + x); /* 1 with inexact iff x != 0 */
98	}	96	}
99		97
100	ENTERI();	98	ENTERI();
Lines 168-175 expm1l(long double x) Link Here
168	int k, n, n2;	166	int k, n, n2;
169	uint16_t hx, ix;	167	uint16_t hx, ix;
170		168
171	DOPRINT_START(&x);
172
173	/* Filter out exceptional cases. */	169	/* Filter out exceptional cases. */
174	u.e = x;	170	u.e = x;
175	hx = u.xbits.expsign;	171	hx = u.xbits.expsign;
Lines 177-187 expm1l(long double x) Link Here
177	if (ix >= BIAS + 6) { /* \|x\| >= 64 or x is NaN */	173	if (ix >= BIAS + 6) { /* \|x\| >= 64 or x is NaN */
178	if (ix == BIAS + LDBL_MAX_EXP) {	174	if (ix == BIAS + LDBL_MAX_EXP) {
179	if (hx & 0x8000) /* x is -Inf, -NaN or unsupported */	175	if (hx & 0x8000) /* x is -Inf, -NaN or unsupported */
180	RETURNP(-1 / x - 1);	176	RETURNF(-1 / x - 1);
181	RETURNP(x + x); /* x is +Inf, +NaN or unsupported */	177	RETURNF(x + x); /* x is +Inf, +NaN or unsupported */
182	}	178	}
183	if (x > o_threshold)	179	if (x > o_threshold)
184	RETURNP(huge * huge);	180	RETURNF(huge * huge);
185	/*	181	/*
186	* expm1l() never underflows, but it must avoid	182	* expm1l() never underflows, but it must avoid
187	* unrepresentable large negative exponents. We used a	183	* unrepresentable large negative exponents. We used a
Lines 190-196 expm1l(long double x) Link Here
190	* in the same way as large ones here.	186	* in the same way as large ones here.
191	*/	187	*/
192	if (hx & 0x8000) /* x <= -64 */	188	if (hx & 0x8000) /* x <= -64 */
193	RETURN2P(tiny, -1); /* good for x < -65ln2 - eps */	189	RETURNF(tiny - 1); /* good for x < -65ln2 - eps */
194	}	190	}
195		191
196	ENTERI();	192	ENTERI();
Lines 198-204 expm1l(long double x) Link Here
198	if (T1 < x && x < T2) {	194	if (T1 < x && x < T2) {
199	if (ix < BIAS - 74) { /* \|x\| < 0x1p-74 (includes pseudos) */	195	if (ix < BIAS - 74) { /* \|x\| < 0x1p-74 (includes pseudos) */
200	/* x (rounded) with inexact if x != 0: */	196	/* x (rounded) with inexact if x != 0: */
201	RETURNPI(x == 0 ? x :	197	RETURNI(x == 0 ? x :
202	(0x1p100 * x + fabsl(x)) * 0x1p-100);	198	(0x1p100 * x + fabsl(x)) * 0x1p-100);
203	}	199	}
204		200
Lines 219-227 expm1l(long double x) Link Here
219	hx2_hi = x_hi * x_hi / 2;	215	hx2_hi = x_hi * x_hi / 2;
220	hx2_lo = x_lo * (x + x_hi) / 2;	216	hx2_lo = x_lo * (x + x_hi) / 2;
221	if (ix >= BIAS - 7)	217	if (ix >= BIAS - 7)
222	RETURN2PI(hx2_hi + x_hi, hx2_lo + x_lo + q);	218	RETURNI((hx2_hi + x_hi) + (hx2_lo + x_lo + q));
223	else	219	else
224	RETURN2PI(x, hx2_lo + q + hx2_hi);	220	RETURNI(x + (hx2_lo + q + hx2_hi));
225	}	221	}
226		222
227	/* Reduce x to (kln2 + endpoint[n2] + r1 + r2). /	223	/* Reduce x to (kln2 + endpoint[n2] + r1 + r2). /

Lines 560-583 log1pl(long double x) Link Here

(-)b/lib/msun/ld80/s_logl.c (-15 / +11 lines)
560	int i, k;	560	int i, k;
561	int16_t ax, hx;	561	int16_t ax, hx;
562		562
563	DOPRINT_START(&x);
564	EXTRACT_LDBL80_WORDS(hx, lx, x);	563	EXTRACT_LDBL80_WORDS(hx, lx, x);
565	if (hx < 0x3fff) { /* x < 1, or x neg NaN */	564	if (hx < 0x3fff) { /* x < 1, or x neg NaN */
566	ax = hx & 0x7fff;	565	ax = hx & 0x7fff;
567	if (ax >= 0x3fff) { /* x <= -1, or x neg NaN */	566	if (ax >= 0x3fff) { /* x <= -1, or x neg NaN */
568	if (ax == 0x3fff && lx == 0x8000000000000000ULL)	567	if (ax == 0x3fff && lx == 0x8000000000000000ULL)
569	RETURNP(-1 / zero); /* log1p(-1) = -Inf */	568	RETURNF(-1 / zero); /* log1p(-1) = -Inf */
570	/* log1p(x < 1, or x [pseudo-]NaN) = qNaN: */	569	/* log1p(x < 1, or x [pseudo-]NaN) = qNaN: */
571	RETURNP((x - x) / (x - x));	570	RETURNF((x - x) / (x - x));
572	}	571	}
573	if (ax <= 0x3fbe) { /* \|x\| < 2*-64 /	572	if (ax <= 0x3fbe) { /* \|x\| < 2*-64 /
574	if ((int)x == 0)	573	if ((int)x == 0)
575	RETURNP(x); /* x with inexact if x != 0 */	574	RETURNF(x); /* x with inexact if x != 0 */
576	}	575	}
577	f_hi = 1;	576	f_hi = 1;
578	f_lo = x;	577	f_lo = x;
579	} else if (hx >= 0x7fff) { /* x +Inf or non-neg NaN */	578	} else if (hx >= 0x7fff) { /* x +Inf or non-neg NaN */
580	RETURNP(x + x); /* log1p(Inf or NaN) = Inf or qNaN */	579	RETURNF(x + x); /* log1p(Inf or NaN) = Inf or qNaN */
581	/* log1p(pseudo-Inf) = qNaN */	580	/* log1p(pseudo-Inf) = qNaN */
582	/* log1p(pseudo-NaN) = qNaN */	581	/* log1p(pseudo-NaN) = qNaN */
583	/* log1p(unnormal) = qNaN */	582	/* log1p(unnormal) = qNaN */
Lines 658-664 log1pl(long double x) Link Here
658	#endif	657	#endif
659		658
660	_3sumF(val_hi, val_lo, F_hi(i) + dk * ln2_hi);	659	_3sumF(val_hi, val_lo, F_hi(i) + dk * ln2_hi);
661	RETURN2PI(val_hi, val_lo);	660	RETURNI(val_hi + val_lo);
662	}	661	}
663		662
664	#ifdef STRUCT_RETURN	663	#ifdef STRUCT_RETURN
Lines 669-675 logl(long double x) Link Here
669	struct ld r;	668	struct ld r;
670		669
671	ENTERI();	670	ENTERI();
672	DOPRINT_START(&x);
673	k_logl(x, &r);	671	k_logl(x, &r);
674	RETURNSPI(&r);	672	RETURNSPI(&r);
675	}	673	}
Lines 690-704 log10l(long double x) Link Here
690	long double hi, lo;	688	long double hi, lo;
691		689
692	ENTERI();	690	ENTERI();
693	DOPRINT_START(&x);
694	k_logl(x, &r);	691	k_logl(x, &r);
695	if (!r.lo_set)	692	if (!r.lo_set)
696	RETURNPI(r.hi);	693	RETURNI(r.hi);
697	_2sumF(r.hi, r.lo);	694	_2sumF(r.hi, r.lo);
698	hi = (float)r.hi;	695	hi = (float)r.hi;
699	lo = r.lo + (r.hi - hi);	696	lo = r.lo + (r.hi - hi);
700	RETURN2PI(invln10_hi * hi,	697	RETURNI(invln10_hi * hi +
701	invln10_lo_plus_hi * lo + invln10_lo * hi);	698	(invln10_lo_plus_hi * lo + invln10_lo * hi));
702	}	699	}
703		700
704	long double	701	long double
Lines 708-722 log2l(long double x) Link Here
708	long double hi, lo;	705	long double hi, lo;
709		706
710	ENTERI();	707	ENTERI();
711	DOPRINT_START(&x);
712	k_logl(x, &r);	708	k_logl(x, &r);
713	if (!r.lo_set)	709	if (!r.lo_set)
714	RETURNPI(r.hi);	710	RETURNI(r.hi);
715	_2sumF(r.hi, r.lo);	711	_2sumF(r.hi, r.lo);
716	hi = (float)r.hi;	712	hi = (float)r.hi;
717	lo = r.lo + (r.hi - hi);	713	lo = r.lo + (r.hi - hi);
718	RETURN2PI(invln2_hi * hi,	714	RETURNI(invln2_hi * hi +
719	invln2_lo_plus_hi * lo + invln2_lo * hi);	715	(invln2_lo_plus_hi * lo + invln2_lo * hi));
720	}	716	}
721		717
722	#endif /* STRUCT_RETURN */	718	#endif /* STRUCT_RETURN */

Lines 65-72 expl(long double x) Link Here

(-)b/lib/msun/ld128/s_expl.c (-16 / +12 lines)
65	int k;	65	int k;
66	uint16_t hx, ix;	66	uint16_t hx, ix;
67		67
68	DOPRINT_START(&x);
69
70	/* Filter out exceptional cases. */	68	/* Filter out exceptional cases. */
71	u.e = x;	69	u.e = x;
72	hx = u.xbits.expsign;	70	hx = u.xbits.expsign;
Lines 74-88 expl(long double x) Link Here
74	if (ix >= BIAS + 13) { /* \|x\| >= 8192 or x is NaN */	72	if (ix >= BIAS + 13) { /* \|x\| >= 8192 or x is NaN */
75	if (ix == BIAS + LDBL_MAX_EXP) {	73	if (ix == BIAS + LDBL_MAX_EXP) {
76	if (hx & 0x8000) /* x is -Inf or -NaN */	74	if (hx & 0x8000) /* x is -Inf or -NaN */
77	RETURNP(-1 / x);	75	RETURNF(-1 / x);
78	RETURNP(x + x); /* x is +Inf or +NaN */	76	RETURNF(x + x); /* x is +Inf or +NaN */
79	}	77	}
80	if (x > o_threshold)	78	if (x > o_threshold)
81	RETURNP(huge * huge);	79	RETURNF(huge * huge);
82	if (x < u_threshold)	80	if (x < u_threshold)
83	RETURNP(tiny * tiny);	81	RETURNF(tiny * tiny);
84	} else if (ix < BIAS - 114) { /* \|x\| < 0x1p-114 */	82	} else if (ix < BIAS - 114) { /* \|x\| < 0x1p-114 */
85	RETURN2P(1, x); /* 1 with inexact iff x != 0 */	83	RETURNF(1 + x); /* 1 with inexact iff x != 0 */
86	}	84	}
87		85
88	ENTERI();	86	ENTERI();
Lines 210-217 expm1l(long double x) Link Here
210	int k, n, n2;	208	int k, n, n2;
211	uint16_t hx, ix;	209	uint16_t hx, ix;
212		210
213	DOPRINT_START(&x);
214
215	/* Filter out exceptional cases. */	211	/* Filter out exceptional cases. */
216	u.e = x;	212	u.e = x;
217	hx = u.xbits.expsign;	213	hx = u.xbits.expsign;
Lines 219-229 expm1l(long double x) Link Here
219	if (ix >= BIAS + 7) { /* \|x\| >= 128 or x is NaN */	215	if (ix >= BIAS + 7) { /* \|x\| >= 128 or x is NaN */
220	if (ix == BIAS + LDBL_MAX_EXP) {	216	if (ix == BIAS + LDBL_MAX_EXP) {
221	if (hx & 0x8000) /* x is -Inf or -NaN */	217	if (hx & 0x8000) /* x is -Inf or -NaN */
222	RETURNP(-1 / x - 1);	218	RETURNF(-1 / x - 1);
223	RETURNP(x + x); /* x is +Inf or +NaN */	219	RETURNF(x + x); /* x is +Inf or +NaN */
224	}	220	}
225	if (x > o_threshold)	221	if (x > o_threshold)
226	RETURNP(huge * huge);	222	RETURNF(huge * huge);
227	/*	223	/*
228	* expm1l() never underflows, but it must avoid	224	* expm1l() never underflows, but it must avoid
229	* unrepresentable large negative exponents. We used a	225	* unrepresentable large negative exponents. We used a
Lines 232-238 expm1l(long double x) Link Here
232	* in the same way as large ones here.	228	* in the same way as large ones here.
233	*/	229	*/
234	if (hx & 0x8000) /* x <= -128 */	230	if (hx & 0x8000) /* x <= -128 */
235	RETURN2P(tiny, -1); /* good for x < -114ln2 - eps */	231	RETURNF(tiny - 1); /* good for x < -114ln2 - eps */
236	}	232	}
237		233
238	ENTERI();	234	ENTERI();
Lines 244-250 expm1l(long double x) Link Here
244	if (x < T3) {	240	if (x < T3) {
245	if (ix < BIAS - 113) { /* \|x\| < 0x1p-113 */	241	if (ix < BIAS - 113) { /* \|x\| < 0x1p-113 */
246	/* x (rounded) with inexact if x != 0: */	242	/* x (rounded) with inexact if x != 0: */
247	RETURNPI(x == 0 ? x :	243	RETURNI(x == 0 ? x :
248	(0x1p200 * x + fabsl(x)) * 0x1p-200);	244	(0x1p200 * x + fabsl(x)) * 0x1p-200);
249	}	245	}
250	q = x * x2 * C3 + x2 * x2 * (C4 + x * (C5 + x * (C6 +	246	q = x * x2 * C3 + x2 * x2 * (C4 + x * (C5 + x * (C6 +
Lines 265-273 expm1l(long double x) Link Here
265	hx2_hi = x_hi * x_hi / 2;	261	hx2_hi = x_hi * x_hi / 2;
266	hx2_lo = x_lo * (x + x_hi) / 2;	262	hx2_lo = x_lo * (x + x_hi) / 2;
267	if (ix >= BIAS - 7)	263	if (ix >= BIAS - 7)
268	RETURN2PI(hx2_hi + x_hi, hx2_lo + x_lo + q);	264	RETURNI((hx2_hi + x_hi) + (hx2_lo + x_lo + q));
269	else	265	else
270	RETURN2PI(x, hx2_lo + q + hx2_hi);	266	RETURNI(x + (hx2_lo + q + hx2_hi));
271	}	267	}
272		268
273	/* Reduce x to (kln2 + endpoint[n2] + r1 + r2). /	269	/* Reduce x to (kln2 + endpoint[n2] + r1 + r2). /

Lines 573-596 log1pl(long double x) Link Here

(-)b/lib/msun/ld128/s_logl.c (-15 / +9 lines)
573	int i, k;	573	int i, k;
574	int16_t ax, hx;	574	int16_t ax, hx;
575		575
576	DOPRINT_START(&x);
577	EXTRACT_LDBL128_WORDS(hx, lx, llx, x);	576	EXTRACT_LDBL128_WORDS(hx, lx, llx, x);
578	if (hx < 0x3fff) { /* x < 1, or x neg NaN */	577	if (hx < 0x3fff) { /* x < 1, or x neg NaN */
579	ax = hx & 0x7fff;	578	ax = hx & 0x7fff;
580	if (ax >= 0x3fff) { /* x <= -1, or x neg NaN */	579	if (ax >= 0x3fff) { /* x <= -1, or x neg NaN */
581	if (ax == 0x3fff && (lx \| llx) == 0)	580	if (ax == 0x3fff && (lx \| llx) == 0)
582	RETURNP(-1 / zero); /* log1p(-1) = -Inf */	581	RETURNF(-1 / zero); /* log1p(-1) = -Inf */
583	/* log1p(x < 1, or x NaN) = qNaN: */	582	/* log1p(x < 1, or x NaN) = qNaN: */
584	RETURNP((x - x) / (x - x));	583	RETURNF((x - x) / (x - x));
585	}	584	}
586	if (ax <= 0x3f8d) { /* \|x\| < 2*-113 /	585	if (ax <= 0x3f8d) { /* \|x\| < 2*-113 /
587	if ((int)x == 0)	586	if ((int)x == 0)
588	RETURNP(x); /* x with inexact if x != 0 */	587	RETURNF(x); /* x with inexact if x != 0 */
589	}	588	}
590	f_hi = 1;	589	f_hi = 1;
591	f_lo = x;	590	f_lo = x;
592	} else if (hx >= 0x7fff) { /* x +Inf or non-neg NaN */	591	} else if (hx >= 0x7fff) { /* x +Inf or non-neg NaN */
593	RETURNP(x + x); /* log1p(Inf or NaN) = Inf or qNaN */	592	RETURNF(x + x); /* log1p(Inf or NaN) = Inf or qNaN */
594	} else if (hx < 0x40e1) { /* 1 <= x < 2*226 /	593	} else if (hx < 0x40e1) { /* 1 <= x < 2*226 /
595	f_hi = x;	594	f_hi = x;
596	f_lo = 1;	595	f_lo = 1;
Lines 669-675 log1pl(long double x) Link Here
669	#endif	668	#endif
670		669
671	_3sumF(val_hi, val_lo, F_hi(i) + dk * ln2_hi);	670	_3sumF(val_hi, val_lo, F_hi(i) + dk * ln2_hi);
672	RETURN2PI(val_hi, val_lo);	671	RETURNI(val_hi + val_lo);
673	}	672	}
674		673
675	#ifdef STRUCT_RETURN	674	#ifdef STRUCT_RETURN
Lines 680-686 logl(long double x) Link Here
680	struct ld r;	679	struct ld r;
681		680
682	ENTERI();	681	ENTERI();
683	DOPRINT_START(&x);
684	k_logl(x, &r);	682	k_logl(x, &r);
685	RETURNSPI(&r);	683	RETURNSPI(&r);
686	}	684	}
Lines 708-722 log10l(long double x) Link Here
708	long double hi, lo;	706	long double hi, lo;
709		707
710	ENTERI();	708	ENTERI();
711	DOPRINT_START(&x);
712	k_logl(x, &r);	709	k_logl(x, &r);
713	if (!r.lo_set)	710	if (!r.lo_set)
714	RETURNPI(r.hi);	711	RETURNI(r.hi);
715	_2sumF(r.hi, r.lo);	712	_2sumF(r.hi, r.lo);
716	hi = (float)r.hi;	713	hi = (float)r.hi;
717	lo = r.lo + (r.hi - hi);	714	lo = r.lo + (r.hi - hi);
718	RETURN2PI(invln10_hi * hi,	715	RETURNI(invln10_hi * hi + (invln10_lo_plus_hi * lo + invln10_lo * hi));
719	invln10_lo_plus_hi * lo + invln10_lo * hi);
720	}	716	}
721		717
722	long double	718	long double
Lines 726-740 log2l(long double x) Link Here
726	long double hi, lo;	722	long double hi, lo;
727		723
728	ENTERI();	724	ENTERI();
729	DOPRINT_START(&x);
730	k_logl(x, &r);	725	k_logl(x, &r);
731	if (!r.lo_set)	726	if (!r.lo_set)
732	RETURNPI(r.hi);	727	RETURNI(r.hi);
733	_2sumF(r.hi, r.lo);	728	_2sumF(r.hi, r.lo);
734	hi = (float)r.hi;	729	hi = (float)r.hi;
735	lo = r.lo + (r.hi - hi);	730	lo = r.lo + (r.hi - hi);
736	RETURN2PI(invln2_hi * hi,	731	RETURNI(invln2_hi * hi + (invln2_lo_plus_hi * lo + invln2_lo * hi));
737	invln2_lo_plus_hi * lo + invln2_lo * hi);
738	}	732	}
739		733
740	#endif /* STRUCT_RETURN */	734	#endif /* STRUCT_RETURN */