@SuperDicq @steph 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
<https://www.gnu.org/licenses/>. */

static void
sincosx_mpn (mp1 si, mp1 co, mp1 xx, mp1 ix)
{
int i;
mp2 s[4], c[4];
mp1 tmp, x;

if (ix == NULL)
{
memset (si, 0, sizeof (mp1));
memset (co, 0, sizeof (mp1));
co[SZ-1] = 1;
memcpy (x, xx, sizeof (mp1));
}
else
mpn_sub_n (x, xx, ix, SZ);

for (i = 0; i < 1 << N; i++)
{
#define add_shift_mulh(d,x,s1,s2,sh,n) \
do { \
if (s2 != NULL) { \
if (sh > 0) { \
assert (sh < mpbpl); \
mpn_lshift (tmp, s1, SZ, sh); \
if (n) \
mpn_sub_n (tmp,tmp,s2+FRAC/mpbpl,SZ); \
else \
mpn_add_n (tmp,tmp,s2+FRAC/mpbpl,SZ); \
} else { \
if (n) \
mpn_sub_n (tmp,s1,s2+FRAC/mpbpl,SZ); \
else \
mpn_add_n (tmp,s1,s2+FRAC/mpbpl,SZ); \
} \
mpn_mul_n(d,tmp,x,SZ); \
} else \
mpn_mul_n(d,s1,x,SZ); \
assert(N+sh < mpbpl); \
if (N+sh > 0) mpn_rshift(d,d,2*SZ,N+sh); \
} while(0)
#define summ(d,ss,s,n) \
do { \
mpn_add_n(tmp,s[1]+FRAC/mpbpl,s[2]+FRAC/mpbpl,SZ); \
mpn_lshift(tmp,tmp,SZ,1); \
mpn_add_n(tmp,tmp,s[0]+FRAC/mpbpl,SZ); \
mpn_add_n(tmp,tmp,s[3]+FRAC/mpbpl,SZ); \
mpn_divmod_1(tmp,tmp,SZ,6); \
if (n) \
mpn_sub_n (d,ss,tmp,SZ); \
else \
mpn_add_n (d,ss,tmp,SZ); \
} while (0)

add_shift_mulh (s[0], x, co, NULL, 0, 0); /* s0 = h * c; */
add_shift_mulh (c[0], x, si, NULL, 0, 0); /* c0 = h * s; */
add_shift_mulh (s[1], x, co, c[0], 1, 1); /* s1 = h * (c - c0/2); */
add_shift_mulh (c[1], x, si, s[0], 1, 0); /* c1 = h * (s + s0/2); */
add_shift_mulh (s[2], x, co, c[1], 1, 1); /* s2 = h * (c - c1/2); */
add_shift_mulh (c[2], x, si, s[1], 1, 0); /* c2 = h * (s + s1/2); */
add_shift_mulh (s[3], x, co, c[2], 0, 1); /* s3 = h * (c - c2); */
add_shift_mulh (c[3], x, si, s[2], 0, 0); /* c3 = h * (s + s2); */
summ (si, si, s, 0); /* s = s + (s0+2*s1+2*s2+s3)/6; */
summ (co, co, c, 1); /* c = c - (c0+2*c1+2*c2+c3)/6; */
}
#undef add_shift_mulh
#undef summ
}