Attachment #243812 for bug #272742

Lines 688-693 irintl(long double x) Link Here

(-)b/lib/msun/src/math_private.h (+53 lines)
688	}	688	}
689	#endif	689	#endif
690		690
		691	/*
		692	* The following are fast floor macros for 0 <= \|x\| < 0x1p(N-1), where
		693	* N is the precision of the type of x. These macros are used in the
		694	* half-cycle trignometric functions (e.g., sinpi(x)).
		695	*/
		696	#define FFLOORF(x, j0, ix) do { \
		697	(j0) = (((ix) >> 23) & 0xff) - 0x7f; \
		698	(ix) &= ~(0x007fffff >> (j0)); \
		699	SET_FLOAT_WORD((x), (ix)); \
		700	} while (0)
		701
		702	#define FFLOOR(x, j0, ix, lx) do { \
		703	(j0) = (((ix) >> 20) & 0x7ff) - 0x3ff; \
		704	if ((j0) < 20) { \
		705	(ix) &= ~(0x000fffff >> (j0)); \
		706	(lx) = 0; \
		707	} else { \
		708	(lx) &= ~((uint32_t)0xffffffff >> ((j0) - 20)); \
		709	} \
		710	INSERT_WORDS((x), (ix), (lx)); \
		711	} while (0)
		712
		713	#define FFLOORL80(x, j0, ix, lx) do { \
		714	j0 = ix - 0x3fff + 1; \
		715	if ((j0) < 32) { \
		716	(lx) = ((lx) >> 32) << 32; \
		717	(lx) &= ~((((lx) << 32)-1) >> (j0)); \
		718	} else { \
		719	uint64_t _m; \
		720	_m = (uint64_t)-1 >> (j0); \
		721	if ((lx) & _m) (lx) &= ~_m; \
		722	} \
		723	INSERT_LDBL80_WORDS((x), (ix), (lx)); \
		724	} while (0)
		725
		726	#define FFLOORL128(x, ai, ar) do { \
		727	union IEEEl2bits u; \
		728	uint64_t m; \
		729	int e; \
		730	u.e = (x); \
		731	e = u.bits.exp - 16383; \
		732	if (e < 48) { \
		733	m = ((1llu << 49) - 1) >> (e + 1); \
		734	u.bits.manh &= ~m; \
		735	u.bits.manl = 0; \
		736	} else { \
		737	m = (uint64_t)-1 >> (e - 48); \
		738	u.bits.manl &= ~m; \
		739	} \
		740	(ai) = u.e; \
		741	(ar) = (x) - (ai); \
		742	} while (0)
		743
691	#ifdef DEBUG	744	#ifdef DEBUG
692	#if defined(__amd64__) \|\| defined(__i386__)	745	#if defined(__amd64__) \|\| defined(__i386__)
693	#define breakpoint() asm("int $3")	746	#define breakpoint() asm("int $3")

Lines 1-5 Link Here

(-)b/lib/msun/src/s_cospif.c (-9 / +7 lines)
1	/*-	1	/*-
2	* Copyright (c) 2017 Steven G. Kargl	2	* Copyright (c) 2017,2023 Steven G. Kargl
3	* All rights reserved.	3	* All rights reserved.
4	*	4	*
5	* Redistribution and use in source and binary forms, with or without	5	* Redistribution and use in source and binary forms, with or without
Lines 71-82 cospif(float x) Link Here
71	return (c);	71	return (c);
72	}	72	}
73		73
74	if (ix < 0x4b000000) { /* 1 <= \|x\| < 0x1p23 */	74	if (ix < 0x4b000000) { /* 1 <= \|x\| < 0x1p23 */
75	/* Determine integer part of ax. */	75	FFLOORF(x, j0, ix); /* Integer part of ax. */
76	j0 = ((ix >> 23) & 0xff) - 0x7f;
77	ix &= ~(0x007fffff >> j0);
78	SET_FLOAT_WORD(x, ix);
79
80	ax -= x;	76	ax -= x;
81	GET_FLOAT_WORD(ix, ax);	77	GET_FLOAT_WORD(ix, ax);
82		78
Lines 103-109 cospif(float x) Link Here
103	return (vzero / vzero);	99	return (vzero / vzero);
104		100
105	/*	101	/*
106	* \|x\| >= 0x1p23 is always an even integer, so return 1.	102	* For 0x1p23 <= \|x\| < 0x1p24 need to determine if x is an even
		103	* or odd integer to return +1 or -1.
		104	* For \|x\| >= 0x1p24, it is always an even integer, so return 1.
107	*/	105	*/
108	return (1);	106	return (ix < 0x4b800000 ? ((ix & 1) ? -1 : 1) : 1);
109	}	107	}




/*-
 * Copyright (c) 2017,2023 Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without

		return ((hx & 0x80000000) ? -s : s);
	}

	if (ix < 0x4b000000) {		/* 1 <= |x| < 0x1p23 */
		FFLOORF(x, j0, ix);	/* Integer part of ax. */
		j0 = ((ix >> 23) & 0xff) - 0x7f;
		ix &= ~(0x007fffff >> j0);
		SET_FLOAT_WORD(x, ix);

		ax -= x;
		GET_FLOAT_WORD(ix, ax);





/*-
 * Copyright (c) 2017,2023 Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without

float
tanpif(float x)
{
	float ax, hi, lo, odd, t;
	uint32_t hx, ix, j0;

	GET_FLOAT_WORD(hx, x);

			}
			t = __kernel_tanpif(ax);
		} else if (ix == 0x3f000000)
			t = 1 / vzero;
		else
			t = - __kernel_tanpif(1 - ax);
		return ((hx & 0x80000000) ? -t : t);
	}

	if (ix < 0x4b000000) {		/* 1 <= |x| < 0x1p23 */
		FFLOORF(x, j0, ix);	/* Integer part of ax. */
		odd = (uint32_t)x & 1 ? -1 : 1;
		ix &= ~(0x007fffff >> j0);
		SET_FLOAT_WORD(x, ix);

		ax -= x;
		GET_FLOAT_WORD(ix, ax);

		if (ix < 0x3f000000)		/* |x| < 0.5 */
			t = ix == 0 ? copysignf(0, odd) : __kernel_tanpif(ax);
		else if (ix == 0x3f000000)
			t = odd / vzero;
		else
			t = - __kernel_tanpif(1 - ax);
		return ((hx & 0x80000000) ? -t : t);

		return (vzero / vzero);

	/*
	 * For 0x1p23 <= |x| < 0x1p24 need to determine if x is an even
	 * or odd integer to set t = +0 or -0.
	 * For |x| >= 0x1p24, it is always an even integer, so t = 0.
	 */
	t = ix >= 0x4b800000 ? 0 : (copysignf(0, (ix & 1) ? -1 : 1));
	return ((hx & 0x80000000) ? -t : t);
}




/*-
 * Copyright (c) 2017, 2023 Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without

	}

	if (ix < 0x43300000) {		/* 1 <= |x| < 0x1p52 */
		FFLOOR(x, j0, ix, lx);	/* Integer part of ax. */
		j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
		if (j0 < 20) {
			ix &= ~(0x000fffff >> j0);
			lx = 0;
		} else {
			lx &= ~((uint32_t)0xffffffff >> (j0 - 20));
		}
		INSERT_WORDS(x, ix, lx);

		ax -= x;
		EXTRACT_WORDS(ix, lx, ax);


		if (ix < 0x3fe00000) {		/* |x| < 0.5 */
			if (ix < 0x3fd00000)	/* |x| < 0.25 */
				c = ix == 0 ? 1 : __kernel_cospi(ax);

		return (vzero / vzero);

	/*
	 * For 0x1p52 <= |x| < 0x1p53 need to determine if x is an even
	 * or odd integer to return +1 or -1.
	 * For |x| >= 0x1p53, it is always an even integer, so return 1.
	 */
	return (ix < 0x43400000 ? ((lx & 1) ? -1 : 1) : 1);
}

#if LDBL_MANT_DIG == 53




/*-
 * Copyright (c) 2017, 2023 Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without

	}

	if (ix < 0x43300000) {			/* 1 <= |x| < 0x1p52 */
		FFLOOR(x, j0, ix, lx);	/* Integer part of ax. */
		j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
		if (j0 < 20) {
			ix &= ~(0x000fffff >> j0);
			lx = 0;
		} else {
			lx &= ~((uint32_t)0xffffffff >> (j0 - 20));
		}
		INSERT_WORDS(x, ix, lx);

		ax -= x;
		EXTRACT_WORDS(ix, lx, ax);


Lines 1-5 Link Here

(-)b/lib/msun/src/s_tanpi.c (-21 / +20 lines)
1	/*-	1	/*-
2	* Copyright (c) 2017 Steven G. Kargl	2	* Copyright (c) 2017, 2023 Steven G. Kargl
3	* All rights reserved.	3	* All rights reserved.
4	*	4	*
5	* Redistribution and use in source and binary forms, with or without	5	* Redistribution and use in source and binary forms, with or without
Lines 56-66 Link Here
56	* 5. Special cases:	56	* 5. Special cases:
57	*	57	*
58	* tanpi(+-0) = +-0	58	* tanpi(+-0) = +-0
59	* tanpi(+-n) = +-0, for positive integers n.	59	* tanpi(n) = +0 for positive even and negative odd integer n.
		60	* tanpi(n) = -0 for positive odd and negative even integer n.
60	* tanpi(+-n+1/4) = +-1, for positive integers n.	61	* tanpi(+-n+1/4) = +-1, for positive integers n.
61	* tanpi(+-n+1/2) = NaN, for positive integers n.	62	* tanpi(n+1/2) = +inf and raises the FE_DIVBYZERO exception for
62	* tanpi(+-inf) = NaN. Raises the "invalid" floating-point exception.	63	* even integers n.
63	* tanpi(nan) = NaN. Raises the "invalid" floating-point exception.	64	* tanpi(n+1/2) = -inf and raises the FE_DIVBYZERO exception for
		65	* odd integers n.
		66	* tanpi(+-inf) = NaN and raises the FE_INVALID exception.
		67	* tanpi(nan) = NaN and raises the FE_INVALID exception.
64	*/	68	*/
65		69
66	#include <float.h>	70	#include <float.h>
Lines 106-112 volatile static const double vzero = 0; Link Here
106	double	110	double
107	tanpi(double x)	111	tanpi(double x)
108	{	112	{
109	double ax, hi, lo, t;	113	double ax, hi, lo, odd, t;
110	uint32_t hx, ix, j0, lx;	114	uint32_t hx, ix, j0, lx;
111		115
112	EXTRACT_WORDS(hx, lx, x);	116	EXTRACT_WORDS(hx, lx, x);
Lines 132-161 tanpi(double x) Link Here
132	}	136	}
133	t = __kernel_tanpi(ax);	137	t = __kernel_tanpi(ax);
134	} else if (ax == 0.5)	138	} else if (ax == 0.5)
135	return ((ax - ax) / (ax - ax));	139	t = 1 / vzero;
136	else	140	else
137	t = - __kernel_tanpi(1 - ax);	141	t = - __kernel_tanpi(1 - ax);
138	return ((hx & 0x80000000) ? -t : t);	142	return ((hx & 0x80000000) ? -t : t);
139	}	143	}
140		144
141	if (ix < 0x43300000) { /* 1 <= \|x\| < 0x1p52 */	145	if (ix < 0x43300000) { /* 1 <= \|x\| < 0x1p52 */
142	/* Determine integer part of ax. */	146	FFLOOR(x, j0, ix, lx); /* Integer part of ax. */
143	j0 = ((ix >> 20) & 0x7ff) - 0x3ff;	147	odd = (uint64_t)x & 1 ? -1 : 1;
144	if (j0 < 20) {
145	ix &= ~(0x000fffff >> j0);
146	lx = 0;
147	} else {
148	lx &= ~(((uint32_t)(0xffffffff)) >> (j0 - 20));
149	}
150	INSERT_WORDS(x,ix,lx);
151
152	ax -= x;	148	ax -= x;
153	EXTRACT_WORDS(ix, lx, ax);	149	EXTRACT_WORDS(ix, lx, ax);
154		150
155	if (ix < 0x3fe00000) /* \|x\| < 0.5 */	151	if (ix < 0x3fe00000) /* \|x\| < 0.5 */
156	t = ax == 0 ? 0 : __kernel_tanpi(ax);	152	t = ix == 0 ? copysign(0, odd) : __kernel_tanpi(ax);
157	else if (ax == 0.5)	153	else if (ax == 0.5)
158	return ((ax - ax) / (ax - ax));	154	t = odd / vzero;
159	else	155	else
160	t = - __kernel_tanpi(1 - ax);	156	t = - __kernel_tanpi(1 - ax);
161		157
Lines 167-175 tanpi(double x) Link Here
167	return (vzero / vzero);	163	return (vzero / vzero);
168		164
169	/*	165	/*
170	* \|x\| >= 0x1p52 is always an integer, so return +-0.	166	* For 0x1p52 <= \|x\| < 0x1p53 need to determine if x is an even
		167	* or odd integer to set t = +0 or -0.
		168	* For \|x\| >= 0x1p54, it is always an even integer, so t = 0.
171	*/	169	*/
172	return (copysign(0, x));	170	t = ix >= 0x43400000 ? 0 : (copysign(0, (lx & 1) ? -1 : 1));
		171	return ((hx & 0x80000000) ? -t : t);
173	}	172	}
174		173
175	#if LDBL_MANT_DIG == 53	174	#if LDBL_MANT_DIG == 53




/*-
 * Copyright (c) 2017, 2023 Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without

		RETURNI(c);
	}

	if (ix < 0x403e) {			/* 1 <= |x| < 0x1p63 */
		FFLOORL80(x, j0, ix, lx);	/* Integer part of ax. */
		j0 = ix - 0x3fff + 1;
		if (j0 < 32) {
			lx = (lx >> 32) << 32;
			lx &= ~(((lx << 32)-1) >> j0);
		} else {
			m = (uint64_t)-1 >> (j0 + 1);
			if (lx & m) lx &= ~m;
		}
		INSERT_LDBL80_WORDS(x, ix, lx);

		ax -= x;
		EXTRACT_LDBL80_WORDS(ix, lx, ax);


		RETURNI(vzero / vzero);

	/*
	 * For 0x1p63 <= |x| < 0x1p64 need to determine if x is an even
	 * or odd integer to return t = +1 or -1.
	 * For |x| >= 0x1p64, it is always an even integer, so t = 1.
	 */
	RETURNI(ix >= 0x403f ? 1 : ((lx & 1) ? -1 : 1));
}

Lines 1-5 Link Here

(-)b/lib/msun/ld80/s_sinpil.c (-13 / +3 lines)
1	/*-	1	/*-
2	* Copyright (c) 2017 Steven G. Kargl	2	* Copyright (c) 2017, 2023 Steven G. Kargl
3	* All rights reserved.	3	* All rights reserved.
4	*	4	*
5	* Redistribution and use in source and binary forms, with or without	5	* Redistribution and use in source and binary forms, with or without
Lines 88-105 sinpil(long double x) Link Here
88	RETURNI((hx & 0x8000) ? -s : s);	88	RETURNI((hx & 0x8000) ? -s : s);
89	}	89	}
90		90
91	if (ix < 0x403e) { /* 1 <= \|x\| < 0x1p63 */	91	if (ix < 0x403e) { /* 1 <= \|x\| < 0x1p63 */
92	/* Determine integer part of ax. */	92	FFLOORL80(x, j0, ix, lx); /* Integer part of ax. */
93	j0 = ix - 0x3fff + 1;
94	if (j0 < 32) {
95	lx = (lx >> 32) << 32;
96	lx &= ~(((lx << 32)-1) >> j0);
97	} else {
98	m = (uint64_t)-1 >> (j0 + 1);
99	if (lx & m) lx &= ~m;
100	}
101	INSERT_LDBL80_WORDS(x, ix, lx);
102
103	ax -= x;	93	ax -= x;
104	EXTRACT_LDBL80_WORDS(ix, lx, ax);	94	EXTRACT_LDBL80_WORDS(ix, lx, ax);
105		95




/*-
 * Copyright (c) 2017, 2023 Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without

long double
tanpil(long double x)
{
	long double ax, hi, lo, odd, t;
	uint64_t lx, m;
	uint32_t j0;
	uint16_t hx, ix;

			}
			t = __kernel_tanpil(ax);
		} else if (ax == 0.5)
			t = 1 / vzero;
		else
			t = -__kernel_tanpil(1 - ax);
		RETURNI((hx & 0x8000) ? -t : t);
	}

	if (ix < 0x403e) {			/* 1 <= |x| < 0x1p63 */
		FFLOORL80(x, j0, ix, lx);	/* Integer part of ax. */
		odd = (uint64_t)x & 1 ? -1 : 1;
		if (j0 < 32) {
			lx = (lx >> 32) << 32;
			lx &= ~(((lx << 32)-1) >> j0);
		} else {
			m = (uint64_t)-1 >> (j0 + 1);
			if (lx & m) lx &= ~m;
		}
		INSERT_LDBL80_WORDS(x, ix, lx);

		ax -= x;
		EXTRACT_LDBL80_WORDS(ix, lx, ax);

		if (ix < 0x3ffe)		/* |x| < 0.5 */
			t = ix == 0 ? copysignl(0, odd) : __kernel_tanpil(ax);
		else if (ax == 0.5L)
			t = odd / vzero;
		else
			t = -__kernel_tanpil(1 - ax);
		RETURNI((hx & 0x8000) ? -t : t);

		RETURNI(vzero / vzero);

	/*
	 * For 0x1p63 <= |x| < 0x1p64 need to determine if x is an even
	 * or odd integer to set t = +0 or -0.
	 * For |x| >= 0x1p64, it is always an even integer, so t = 0.
	 */
	t = ix >= 0x403f ? 0 : (copysignl(0, (lx & 1) ? -1 : 1));
	RETURNI((hx & 0x8000) ? -t : t);
}




/*-
 * Copyright (c) 2017-2023 Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without

 * See ../src/s_cospi.c for implementation details.
 */

#include "fpmath.h"
#include "math.h"
#include "math_private.h"


long double
cospil(long double x)
{
	long double ai, ar, ax, c;
	uint32_t ix;

	ax = fabsl(x);


	}

	if (ax < 0x1p112) {
		/* Split ax = ai + ar with 0 <= ar < 1. */
		FFLOORL128(ax, ai, ar);
		ax -= xf;                               /* Remainder */
		if (ax < 0) {
			ax += 1;
			xf -= 1;
		}

		if (ar < 0.5) {
			if (ar < 0.25)
				c = ar == 0 ? 1 : __kernel_cospil(ar);
			else
				c = __kernel_sinpil(0.5 - ar);
		} else {
			if (ar < 0.75) {
				if (ar == 0.5)
					return (0);
				c = -__kernel_sinpil(ar - 0.5);
			} else
				c = -__kernel_cospil(1 - ar);
		}
		return (fmodl(ai, 2.L) == 0 ? c : -c);
		if (xf > 0x1p64)
			xf -= 0x1p64;
		if (xf > 0x1p32)
			xf -= 0x1p32;
		ix = (uint32_t)xf;
		return (ix & 1 ? -c : c);
	}

	if (isinf(x) || isnan(x))
		return (vzero / vzero);

	/*
	 * For |x| >= 0x1p113, it is always an even integer, so return 1.
	 */
	if (ax >= 0x1p113)
		return (1);
	/*
	 * For 0x1p112 <= |x| < 0x1p113 need to determine if x is an even
	 * or odd integer to return 1 or -1.
	 */

	return (fmodl(ax, 2.L) == 0 ? 1 : -1);
}




/*-
 * Copyright (c) 2017-2023 Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without

 * See ../src/s_sinpi.c for implementation details.
 */

#include "fpmath.h"
#include "math.h"
#include "math_private.h"


long double
sinpil(long double x)
{
	long double ai, ar, ax, hi, lo, s, xhi, xlo;
	uint32_t ix;

	ax = fabsl(x);


	}

	if (ax < 0x1p112) {
		/* Split ax = ai + ar with 0 <= ar < 1. */
		FFLOORL128(ax, ai, ar);
		ax -= xf;                               /* Remainder */
		if (ax < 0) {
			ax += 1;
			xf -= 1;
		}

		if (ar == 0) {
			s = 0;
		} else {
			if (ar < 0.5) {
				if (ar <= 0.25)
					s = __kernel_sinpil(ar);
				else
					s = __kernel_cospil(0.5 - ar);
			} else {
				if (ar < 0.75)
					s = __kernel_cospil(ar - 0.5);
				else
					s = __kernel_sinpil(1 - ar);
			}

			s = fmodl(ai, 2.L) == 0 ? s : -s;
				xf -= 0x1p64;
			if (xf > 0x1p32)
				xf -= 0x1p32;
			ix = (uint32_t)xf;
			if (ix & 1) s = -s;
		}
		return (x < 0 ? -s : s);
	}




/*-
 * Copyright (c) 2017-2023 Steven G. Kargl
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without

 * See ../src/s_tanpi.c for implementation details.
 */

#include "fpmath.h"
#include "math.h"
#include "math_private.h"


long double
tanpil(long double x)
{
	long double ai, ar, ax, hi, lo, t;
	double odd;

	ax = fabsl(x);


			}
			t = __kernel_tanpil(ax);
		} else if (ax == 0.5)
			t = 1 / vzero;
		else
			t = -__kernel_tanpil(1 - ax);
		return (x < 0 ? -t : t);
	}

	if (ax < 0x1p112) {
		/* Split ax = ai + ar with 0 <= ar < 1. */
		FFLOORL128(ax, ai, ar);
		odd = fmodl(ai, 2.L) == 0 ? 1 : -1;
		if (ar < 0.5)
			t = ar == 0 ? copysign(0., odd) : __kernel_tanpil(ar);
		else if (ar == 0.5)
			t = odd / vzero;

		if (ax < 0.5)
			t = ax == 0 ? 0 : __kernel_tanpil(ax);
		else if (ax == 0.5)
			return ((ax - ax) / (ax - ax));
		else
			t = -__kernel_tanpil(1 - ar);
		return (x < 0 ? -t : t);
	}


		return (vzero / vzero);

	/*
	 * For 0x1p112 <= |x| < 0x1p113 need to determine if x is an even
	 * or odd integer to set t = +0 or -0.
	 * For |x| >= 0x1p113, it is always an even integer, so t = 0.
	 */
	t = fmodl(ax,2.L) == 0  ? 0 : copysign(0., -1.);
	return (copysignl(t, x));
}

Return to bug 272742

Lines 1-5 Link Here

(-)b/lib/msun/src/s_sinpif.c (-7 / +3 lines)
1	/*-	1	/*-
2	* Copyright (c) 2017 Steven G. Kargl	2	* Copyright (c) 2017,2023 Steven G. Kargl
3	* All rights reserved.	3	* All rights reserved.
4	*	4	*
5	* Redistribution and use in source and binary forms, with or without	5	* Redistribution and use in source and binary forms, with or without
Lines 81-92 sinpif(float x) Link Here
81	return ((hx & 0x80000000) ? -s : s);	81	return ((hx & 0x80000000) ? -s : s);
82	}	82	}
83		83
84	if (ix < 0x4b000000) { /* 1 <= \|x\| < 0x1p23 */	84	if (ix < 0x4b000000) { /* 1 <= \|x\| < 0x1p23 */
85	/* Determine integer part of ax. */	85	FFLOORF(x, j0, ix); /* Integer part of ax. */
86	j0 = ((ix >> 23) & 0xff) - 0x7f;
87	ix &= ~(0x007fffff >> j0);
88	SET_FLOAT_WORD(x, ix);
89
90	ax -= x;	86	ax -= x;
91	GET_FLOAT_WORD(ix, ax);	87	GET_FLOAT_WORD(ix, ax);
92		88

Lines 1-5 Link Here

(-)b/lib/msun/src/s_tanpif.c (-12 / +12 lines)
1	/*-	1	/*-
2	* Copyright (c) 2017 Steven G. Kargl	2	* Copyright (c) 2017,2023 Steven G. Kargl
3	* All rights reserved.	3	* All rights reserved.
4	*	4	*
5	* Redistribution and use in source and binary forms, with or without	5	* Redistribution and use in source and binary forms, with or without
Lines 58-64 volatile static const float vzero = 0; Link Here
58	float	58	float
59	tanpif(float x)	59	tanpif(float x)
60	{	60	{
61	float ax, hi, lo, t;	61	float ax, hi, lo, odd, t;
62	uint32_t hx, ix, j0;	62	uint32_t hx, ix, j0;
63		63
64	GET_FLOAT_WORD(hx, x);	64	GET_FLOAT_WORD(hx, x);
Lines 79-103 tanpif(float x) Link Here
79	}	79	}
80	t = __kernel_tanpif(ax);	80	t = __kernel_tanpif(ax);
81	} else if (ix == 0x3f000000)	81	} else if (ix == 0x3f000000)
82	return ((ax - ax) / (ax - ax));	82	t = 1 / vzero;
83	else	83	else
84	t = - __kernel_tanpif(1 - ax);	84	t = - __kernel_tanpif(1 - ax);
85	return ((hx & 0x80000000) ? -t : t);	85	return ((hx & 0x80000000) ? -t : t);
86	}	86	}
87		87
88	if (ix < 0x4b000000) { /* 1 <= \|x\| < 0x1p23 */	88	if (ix < 0x4b000000) { /* 1 <= \|x\| < 0x1p23 */
89	/* Determine integer part of ax. */	89	FFLOORF(x, j0, ix); /* Integer part of ax. */
90	j0 = ((ix >> 23) & 0xff) - 0x7f;	90	odd = (uint32_t)x & 1 ? -1 : 1;
91	ix &= ~(0x007fffff >> j0);
92	SET_FLOAT_WORD(x, ix);
93
94	ax -= x;	91	ax -= x;
95	GET_FLOAT_WORD(ix, ax);	92	GET_FLOAT_WORD(ix, ax);
96		93
97	if (ix < 0x3f000000) /* \|x\| < 0.5 */	94	if (ix < 0x3f000000) /* \|x\| < 0.5 */
98	t = ix == 0 ? 0 : __kernel_tanpif(ax);	95	t = ix == 0 ? copysignf(0, odd) : __kernel_tanpif(ax);
99	else if (ix == 0x3f000000)	96	else if (ix == 0x3f000000)
100	return ((ax - ax) / (ax - ax));	97	t = odd / vzero;
101	else	98	else
102	t = - __kernel_tanpif(1 - ax);	99	t = - __kernel_tanpif(1 - ax);
103	return ((hx & 0x80000000) ? -t : t);	100	return ((hx & 0x80000000) ? -t : t);
Lines 108-114 tanpif(float x) Link Here
108	return (vzero / vzero);	105	return (vzero / vzero);
109		106
110	/*	107	/*
111	* \|x\| >= 0x1p23 is always an integer, so return +-0.	108	* For 0x1p23 <= \|x\| < 0x1p24 need to determine if x is an even
		109	* or odd integer to set t = +0 or -0.
		110	* For \|x\| >= 0x1p24, it is always an even integer, so t = 0.
112	*/	111	*/
113	return (copysignf(0, x));	112	t = ix >= 0x4b800000 ? 0 : (copysignf(0, (ix & 1) ? -1 : 1));
		113	return ((hx & 0x80000000) ? -t : t);
114	}	114	}

Lines 1-5 Link Here

(-)b/lib/msun/src/s_cospi.c (-14 / +6 lines)
1	/*-	1	/*-
2	* Copyright (c) 2017 Steven G. Kargl	2	* Copyright (c) 2017, 2023 Steven G. Kargl
3	* All rights reserved.	3	* All rights reserved.
4	*	4	*
5	* Redistribution and use in source and binary forms, with or without	5	* Redistribution and use in source and binary forms, with or without
Lines 104-123 cospi(double x) Link Here
104	}	104	}
105		105
106	if (ix < 0x43300000) { /* 1 <= \|x\| < 0x1p52 */	106	if (ix < 0x43300000) { /* 1 <= \|x\| < 0x1p52 */
107	/* Determine integer part of ax. */	107	FFLOOR(x, j0, ix, lx); /* Integer part of ax. */
108	j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
109	if (j0 < 20) {
110	ix &= ~(0x000fffff >> j0);
111	lx = 0;
112	} else {
113	lx &= ~((uint32_t)0xffffffff >> (j0 - 20));
114	}
115	INSERT_WORDS(x, ix, lx);
116
117	ax -= x;	108	ax -= x;
118	EXTRACT_WORDS(ix, lx, ax);	109	EXTRACT_WORDS(ix, lx, ax);
119		110
120
121	if (ix < 0x3fe00000) { /* \|x\| < 0.5 */	111	if (ix < 0x3fe00000) { /* \|x\| < 0.5 */
122	if (ix < 0x3fd00000) /* \|x\| < 0.25 */	112	if (ix < 0x3fd00000) /* \|x\| < 0.25 */
123	c = ix == 0 ? 1 : __kernel_cospi(ax);	113	c = ix == 0 ? 1 : __kernel_cospi(ax);
Lines 143-151 cospi(double x) Link Here
143	return (vzero / vzero);	133	return (vzero / vzero);
144		134
145	/*	135	/*
146	* \|x\| >= 0x1p52 is always an even integer, so return 1.	136	* For 0x1p52 <= \|x\| < 0x1p53 need to determine if x is an even
		137	* or odd integer to return +1 or -1.
		138	* For \|x\| >= 0x1p53, it is always an even integer, so return 1.
147	*/	139	*/
148	return (1);	140	return (ix < 0x43400000 ? ((lx & 1) ? -1 : 1) : 1);
149	}	141	}
150		142
151	#if LDBL_MANT_DIG == 53	143	#if LDBL_MANT_DIG == 53

Lines 1-5 Link Here

(-)b/lib/msun/src/s_sinpi.c (-11 / +2 lines)
1	/*-	1	/*-
2	* Copyright (c) 2017 Steven G. Kargl	2	* Copyright (c) 2017, 2023 Steven G. Kargl
3	* All rights reserved.	3	* All rights reserved.
4	*	4	*
5	* Redistribution and use in source and binary forms, with or without	5	* Redistribution and use in source and binary forms, with or without
Lines 118-133 sinpi(double x) Link Here
118	}	118	}
119		119
120	if (ix < 0x43300000) { /* 1 <= \|x\| < 0x1p52 */	120	if (ix < 0x43300000) { /* 1 <= \|x\| < 0x1p52 */
121	/* Determine integer part of ax. */	121	FFLOOR(x, j0, ix, lx); /* Integer part of ax. */
122	j0 = ((ix >> 20) & 0x7ff) - 0x3ff;
123	if (j0 < 20) {
124	ix &= ~(0x000fffff >> j0);
125	lx = 0;
126	} else {
127	lx &= ~((uint32_t)0xffffffff >> (j0 - 20));
128	}
129	INSERT_WORDS(x, ix, lx);
130
131	ax -= x;	122	ax -= x;
132	EXTRACT_WORDS(ix, lx, ax);	123	EXTRACT_WORDS(ix, lx, ax);
133		124

Lines 1-5 Link Here

(-)b/lib/msun/ld80/s_cospil.c (-15 / +7 lines)
1	/*-	1	/*-
2	* Copyright (c) 2017 Steven G. Kargl	2	* Copyright (c) 2017, 2023 Steven G. Kargl
3	* All rights reserved.	3	* All rights reserved.
4	*	4	*
5	* Redistribution and use in source and binary forms, with or without	5	* Redistribution and use in source and binary forms, with or without
Lines 80-97 cospil(long double x) Link Here
80	RETURNI(c);	80	RETURNI(c);
81	}	81	}
82		82
83	if (ix < 0x403e) { /* 1 <= \|x\| < 0x1p63 */	83	if (ix < 0x403e) { /* 1 <= \|x\| < 0x1p63 */
84	/* Determine integer part of ax. */	84	FFLOORL80(x, j0, ix, lx); /* Integer part of ax. */
85	j0 = ix - 0x3fff + 1;
86	if (j0 < 32) {
87	lx = (lx >> 32) << 32;
88	lx &= ~(((lx << 32)-1) >> j0);
89	} else {
90	m = (uint64_t)-1 >> (j0 + 1);
91	if (lx & m) lx &= ~m;
92	}
93	INSERT_LDBL80_WORDS(x, ix, lx);
94
95	ax -= x;	85	ax -= x;
96	EXTRACT_LDBL80_WORDS(ix, lx, ax);	86	EXTRACT_LDBL80_WORDS(ix, lx, ax);
97		87
Lines 123-129 cospil(long double x) Link Here
123	RETURNI(vzero / vzero);	113	RETURNI(vzero / vzero);
124		114
125	/*	115	/*
126	* \|x\| >= 0x1p63 is always an even integer, so return 1.	116	* For 0x1p63 <= \|x\| < 0x1p64 need to determine if x is an even
		117	* or odd integer to return t = +1 or -1.
		118	* For \|x\| >= 0x1p64, it is always an even integer, so t = 1.
127	*/	119	*/
128	RETURNI(1);	120	RETURNI(ix >= 0x403f ? 1 : ((lx & 1) ? -1 : 1));
129	}	121	}

Lines 1-5 Link Here

(-)b/lib/msun/ld80/s_tanpil.c (-20 / +14 lines)
1	/*-	1	/*-
2	* Copyright (c) 2017 Steven G. Kargl	2	* Copyright (c) 2017, 2023 Steven G. Kargl
3	* All rights reserved.	3	* All rights reserved.
4	*	4	*
5	* Redistribution and use in source and binary forms, with or without	5	* Redistribution and use in source and binary forms, with or without
Lines 72-78 volatile static const double vzero = 0; Link Here
72	long double	72	long double
73	tanpil(long double x)	73	tanpil(long double x)
74	{	74	{
75	long double ax, hi, lo, t;	75	long double ax, hi, lo, odd, t;
76	uint64_t lx, m;	76	uint64_t lx, m;
77	uint32_t j0;	77	uint32_t j0;
78	uint16_t hx, ix;	78	uint16_t hx, ix;
Lines 98-128 tanpil(long double x) Link Here
98	}	98	}
99	t = __kernel_tanpil(ax);	99	t = __kernel_tanpil(ax);
100	} else if (ax == 0.5)	100	} else if (ax == 0.5)
101	RETURNI((ax - ax) / (ax - ax));	101	t = 1 / vzero;
102	else	102	else
103	t = -__kernel_tanpil(1 - ax);	103	t = -__kernel_tanpil(1 - ax);
104	RETURNI((hx & 0x8000) ? -t : t);	104	RETURNI((hx & 0x8000) ? -t : t);
105	}	105	}
106		106
107	if (ix < 0x403e) { /* 1 <= \|x\| < 0x1p63 */	107	if (ix < 0x403e) { /* 1 <= \|x\| < 0x1p63 */
108	/* Determine integer part of ax. */	108	FFLOORL80(x, j0, ix, lx); /* Integer part of ax. */
109	j0 = ix - 0x3fff + 1;	109	odd = (uint64_t)x & 1 ? -1 : 1;
110	if (j0 < 32) {
111	lx = (lx >> 32) << 32;
112	lx &= ~(((lx << 32)-1) >> j0);
113	} else {
114	m = (uint64_t)-1 >> (j0 + 1);
115	if (lx & m) lx &= ~m;
116	}
117	INSERT_LDBL80_WORDS(x, ix, lx);
118
119	ax -= x;	110	ax -= x;
120	EXTRACT_LDBL80_WORDS(ix, lx, ax);	111	EXTRACT_LDBL80_WORDS(ix, lx, ax);
121		112
122	if (ix < 0x3ffe) /* \|x\| < 0.5 */	113	if (ix < 0x3ffe) /* \|x\| < 0.5 */
123	t = ax == 0 ? 0 : __kernel_tanpil(ax);	114	t = ix == 0 ? copysignl(0, odd) : __kernel_tanpil(ax);
124	else if (ax == 0.5)	115	else if (ax == 0.5L)
125	RETURNI((ax - ax) / (ax - ax));	116	t = odd / vzero;
126	else	117	else
127	t = -__kernel_tanpil(1 - ax);	118	t = -__kernel_tanpil(1 - ax);
128	RETURNI((hx & 0x8000) ? -t : t);	119	RETURNI((hx & 0x8000) ? -t : t);
Lines 133-139 tanpil(long double x) Link Here
133	RETURNI(vzero / vzero);	124	RETURNI(vzero / vzero);
134		125
135	/*	126	/*
136	* \|x\| >= 0x1p63 is always an integer, so return +-0.	127	* For 0x1p63 <= \|x\| < 0x1p64 need to determine if x is an even
		128	* or odd integer to set t = +0 or -0.
		129	* For \|x\| >= 0x1p64, it is always an even integer, so t = 0.
137	*/	130	*/
138	RETURNI(copysignl(0, x));	131	t = ix >= 0x403f ? 0 : (copysignl(0, (lx & 1) ? -1 : 1));
		132	RETURNI((hx & 0x8000) ? -t : t);
139	}	133	}

Lines 1-5 Link Here

(-)b/lib/msun/ld128/s_cospil.c (-27 / +23 lines)
1	/*-	1	/*-
2	* Copyright (c) 2017-2021 Steven G. Kargl	2	* Copyright (c) 2017-2023 Steven G. Kargl
3	* All rights reserved.	3	* All rights reserved.
4	*	4	*
5	* Redistribution and use in source and binary forms, with or without	5	* Redistribution and use in source and binary forms, with or without
Lines 28-33 Link Here
28	* See ../src/s_cospi.c for implementation details.	28	* See ../src/s_cospi.c for implementation details.
29	*/	29	*/
30		30
		31	#include "fpmath.h"
31	#include "math.h"	32	#include "math.h"
32	#include "math_private.h"	33	#include "math_private.h"
33		34
Lines 46-53 volatile static const double vzero = 0; Link Here
46	long double	47	long double
47	cospil(long double x)	48	cospil(long double x)
48	{	49	{
49	long double ax, c, xf;	50	long double ai, ar, ax, c;
50	uint32_t ix;
51		51
52	ax = fabsl(x);	52	ax = fabsl(x);
53		53
Lines 72-112 cospil(long double x) Link Here
72	}	72	}
73		73
74	if (ax < 0x1p112) {	74	if (ax < 0x1p112) {
75	/* Split x = n + r with 0 <= r < 1. */	75	/* Split ax = ai + ar with 0 <= ar < 1. */
76	xf = (ax + 0x1p112L) - 0x1p112L; /* Integer part */	76	FFLOORL128(ax, ai, ar);
77	ax -= xf; /* Remainder */
78	if (ax < 0) {
79	ax += 1;
80	xf -= 1;
81	}
82		77
83	if (ax < 0.5) {	78	if (ar < 0.5) {
84	if (ax < 0.25)	79	if (ar < 0.25)
85	c = ax == 0 ? 1 : __kernel_cospil(ax);	80	c = ar == 0 ? 1 : __kernel_cospil(ar);
86	else	81	else
87	c = __kernel_sinpil(0.5 - ax);	82	c = __kernel_sinpil(0.5 - ar);
88	} else {	83	} else {
89	if (ax < 0.75) {	84	if (ar < 0.75) {
90	if (ax == 0.5)	85	if (ar == 0.5)
91	return (0);	86	return (0);
92	c = -__kernel_sinpil(ax - 0.5);	87	c = -__kernel_sinpil(ar - 0.5);
93	} else	88	} else
94	c = -__kernel_cospil(1 - ax);	89	c = -__kernel_cospil(1 - ar);
95	}	90	}
96		91	return (fmodl(ai, 2.L) == 0 ? c : -c);
97	if (xf > 0x1p64)
98	xf -= 0x1p64;
99	if (xf > 0x1p32)
100	xf -= 0x1p32;
101	ix = (uint32_t)xf;
102	return (ix & 1 ? -c : c);
103	}	92	}
104		93
105	if (isinf(x) \|\| isnan(x))	94	if (isinf(x) \|\| isnan(x))
106	return (vzero / vzero);	95	return (vzero / vzero);
107		96
108	/*	97	/*
109	* \|x\| >= 0x1p112 is always an even integer, so return 1.	98	* For \|x\| >= 0x1p113, it is always an even integer, so return 1.
110	*/	99	*/
111	return (1);	100	if (ax >= 0x1p113)
		101	return (1);
		102	/*
		103	* For 0x1p112 <= \|x\| < 0x1p113 need to determine if x is an even
		104	* or odd integer to return 1 or -1.
		105	*/
		106
		107	return (fmodl(ax, 2.L) == 0 ? 1 : -1);
112	}	108	}

Lines 1-5 Link Here

(-)b/lib/msun/ld128/s_sinpil.c (-24 / +14 lines)
1	/*-	1	/*-
2	* Copyright (c) 2017-2021 Steven G. Kargl	2	* Copyright (c) 2017-2023 Steven G. Kargl
3	* All rights reserved.	3	* All rights reserved.
4	*	4	*
5	* Redistribution and use in source and binary forms, with or without	5	* Redistribution and use in source and binary forms, with or without
Lines 28-33 Link Here
28	* See ../src/s_sinpi.c for implementation details.	28	* See ../src/s_sinpi.c for implementation details.
29	*/	29	*/
30		30
		31	#include "fpmath.h"
31	#include "math.h"	32	#include "math.h"
32	#include "math_private.h"	33	#include "math_private.h"
33		34
Lines 46-53 volatile static const double vzero = 0; Link Here
46	long double	47	long double
47	sinpil(long double x)	48	sinpil(long double x)
48	{	49	{
49	long double ax, hi, lo, s, xf, xhi, xlo;	50	long double ai, ar, ax, hi, lo, s, xhi, xlo;
50	uint32_t ix;
51		51
52	ax = fabsl(x);	52	ax = fabsl(x);
53		53
Lines 78-112 sinpil(long double x) Link Here
78	}	78	}
79		79
80	if (ax < 0x1p112) {	80	if (ax < 0x1p112) {
81	/* Split x = n + r with 0 <= r < 1. */	81	/* Split ax = ai + ar with 0 <= ar < 1. */
82	xf = (ax + 0x1p112L) - 0x1p112L; /* Integer part */	82	FFLOORL128(ax, ai, ar);
83	ax -= xf; /* Remainder */
84	if (ax < 0) {
85	ax += 1;
86	xf -= 1;
87	}
88		83
89	if (ax == 0) {	84	if (ar == 0) {
90	s = 0;	85	s = 0;
91	} else {	86	} else {
92	if (ax < 0.5) {	87	if (ar < 0.5) {
93	if (ax <= 0.25)	88	if (ar <= 0.25)
94	s = __kernel_sinpil(ax);	89	s = __kernel_sinpil(ar);
95	else	90	else
96	s = __kernel_cospil(0.5 - ax);	91	s = __kernel_cospil(0.5 - ar);
97	} else {	92	} else {
98	if (ax < 0.75)	93	if (ar < 0.75)
99	s = __kernel_cospil(ax - 0.5);	94	s = __kernel_cospil(ar - 0.5);
100	else	95	else
101	s = __kernel_sinpil(1 - ax);	96	s = __kernel_sinpil(1 - ar);
102	}	97	}
103		98
104	if (xf > 0x1p64)	99	s = fmodl(ai, 2.L) == 0 ? s : -s;
105	xf -= 0x1p64;
106	if (xf > 0x1p32)
107	xf -= 0x1p32;
108	ix = (uint32_t)xf;
109	if (ix & 1) s = -s;
110	}	100	}
111	return (x < 0 ? -s : s);	101	return (x < 0 ? -s : s);
112	}	102	}

Lines 1-5 Link Here

(-)b/lib/msun/ld128/s_tanpil.c (-20 / +19 lines)
1	/*-	1	/*-
2	* Copyright (c) 2017-2021 Steven G. Kargl	2	* Copyright (c) 2017-2023 Steven G. Kargl
3	* All rights reserved.	3	* All rights reserved.
4	*	4	*
5	* Redistribution and use in source and binary forms, with or without	5	* Redistribution and use in source and binary forms, with or without
Lines 28-33 Link Here
28	* See ../src/s_tanpi.c for implementation details.	28	* See ../src/s_tanpi.c for implementation details.
29	*/	29	*/
30		30
		31	#include "fpmath.h"
31	#include "math.h"	32	#include "math.h"
32	#include "math_private.h"	33	#include "math_private.h"
33		34
Lines 69-76 volatile static const double vzero = 0; Link Here
69	long double	70	long double
70	tanpil(long double x)	71	tanpil(long double x)
71	{	72	{
72	long double ax, hi, lo, xf, t;	73	long double ai, ar, ax, hi, lo, t;
73	uint32_t ix;	74	double odd;
74		75
75	ax = fabsl(x);	76	ax = fabsl(x);
76		77
Lines 88-114 tanpil(long double x) Link Here
88	}	89	}
89	t = __kernel_tanpil(ax);	90	t = __kernel_tanpil(ax);
90	} else if (ax == 0.5)	91	} else if (ax == 0.5)
91	return ((ax - ax) / (ax - ax));	92	t = 1 / vzero;
92	else	93	else
93	t = -__kernel_tanpil(1 - ax);	94	t = -__kernel_tanpil(1 - ax);
94	return (x < 0 ? -t : t);	95	return (x < 0 ? -t : t);
95	}	96	}
96		97
97	if (ix < 0x1p112) {	98	if (ax < 0x1p112) {
98	/* Split x = n + r with 0 <= r < 1. */	99	/* Split ax = ai + ar with 0 <= ar < 1. */
99	xf = (ax + 0x1p112L) - 0x1p112L; /* Integer part */	100	FFLOORL128(ax, ai, ar);
100	ax -= xf; /* Remainder */	101	odd = fmodl(ai, 2.L) == 0 ? 1 : -1;
101	if (ax < 0) {	102	if (ar < 0.5)
102	ax += 1;	103	t = ar == 0 ? copysign(0., odd) : __kernel_tanpil(ar);
103	xf -= 1;	104	else if (ar == 0.5)
104	}	105	t = odd / vzero;
105
106	if (ax < 0.5)
107	t = ax == 0 ? 0 : __kernel_tanpil(ax);
108	else if (ax == 0.5)
109	return ((ax - ax) / (ax - ax));
110	else	106	else
111	t = -__kernel_tanpil(1 - ax);	107	t = -__kernel_tanpil(1 - ar);
112	return (x < 0 ? -t : t);	108	return (x < 0 ? -t : t);
113	}	109	}
114		110
Lines 117-123 tanpil(long double x) Link Here
117	return (vzero / vzero);	113	return (vzero / vzero);
118		114
119	/*	115	/*
120	* \|x\| >= 0x1p112 is always an integer, so return +-0.	116	* For 0x1p112 <= \|x\| < 0x1p113 need to determine if x is an even
		117	* or odd integer to set t = +0 or -0.
		118	* For \|x\| >= 0x1p113, it is always an even integer, so t = 0.
121	*/	119	*/
122	return (copysignl(0, x));	120	t = fmodl(ax,2.L) == 0 ? 0 : copysign(0., -1.);
		121	return (copysignl(t, x));
123	}	122	}