text: //          Copyright Naoki Shibata 2010 - 2019.
  // Distributed under the Boost Software License, Version 1.0.
  //    (See accompanying file LICENSE.txt or copy at
  //          http://www.boost.org/LICENSE_1_0.txt)
  
  // Always use -ffp-contract=off option to compile SLEEF.
  
  #include <stdio.h>
  
  #include <assert.h>
  #include <stdint.h>
  #include <limits.h>
  
  #include "misc.h"
  
  #ifdef DORENAME
  #include "rename.h"
  #endif
  
  #if (defined(_MSC_VER))
  #pragma fp_contract (off)
  #endif
  
  static INLINE CONST long double mlal(long double x, long double y, long double z) { return x * y + z; }
  static INLINE CONST long double xrintl(long double x) { return x < 0 ? (int)(x - 0.5) : (int)(x + 0.5); }
  static INLINE CONST int64_t xceill(long double x) { return (int64_t)x + (x < 0 ? 0 : 1); }
  static INLINE CONST long double xtruncl(long double x) { return (long double)(int)x; }
  
  static INLINE CONST int xisnanl(long double x) { return x != x; }
  static INLINE CONST int xisinfl(long double x) { return x == SLEEF_INFINITYl || x == -SLEEF_INFINITYl; }
  static INLINE CONST int xisminfl(long double x) { return x == -SLEEF_INFINITYl; }
  static INLINE CONST int xispinfl(long double x) { return x == SLEEF_INFINITYl; }
  
  static INLINE CONST long double xfabsl(long double x) { return x >= 0 ? x : -x; }
  
  //
  
  #ifndef NDEBUG
  static int checkfp(long double x) {
    if (xisinfl(x) || xisnanl(x)) return 1;
    return 0;
  }
  #endif
  
  static INLINE CONST long double upperl(long double d) {
    union {
      long double ld;
      uint32_t u[4];
    } cnv;
  
    cnv.ld = d;
    cnv.u[0] = 0;
    return cnv.ld;
  }
  
  static INLINE CONST Sleef_longdouble2 dl(long double h, long double l) {
    Sleef_longdouble2 ret;
    ret.x = h; ret.y = l;
    return ret;
  }
  
  static INLINE CONST Sleef_longdouble2 dlnormalize_l2_l2(Sleef_longdouble2 t) {
    Sleef_longdouble2 s;
  
    s.x = t.x + t.y;
    s.y = t.x - s.x + t.y;
  
    return s;
  }
  
  static INLINE CONST Sleef_longdouble2 dlscale_l2_l2_l(Sleef_longdouble2 d, long double s) {
    Sleef_longdouble2 r;
  
    r.x = d.x * s;
    r.y = d.y * s;
  
    return r;
  }
  
  static INLINE CONST Sleef_longdouble2 dlneg_l2_l2(Sleef_longdouble2 d) {
    Sleef_longdouble2 r;
  
    r.x = -d.x;
    r.y = -d.y;
  
    return r;
  }
  
  static INLINE CONST Sleef_longdouble2 dladd_l2_l_l(long double x, long double y) {
    // |x| >= |y|
  
    Sleef_longdouble2 r;
  
  #ifndef NDEBUG
    if (!(checkfp(x) || checkfp(y) || xfabsl(x) >= xfabsl(y))) {
      fprintf(stderr, "[dladd_l2_l_l : %Lg, %Lg]\n", x, y);
      fflush(stderr);
    }
  #endif
  
    r.x = x + y;
    r.y = x - r.x + y;
  
    return r;
  }
  
  static INLINE CONST Sleef_longdouble2 dladd2_l2_l_l(long double x, long double y) {
    Sleef_longdouble2 r;
  
    r.x = x + y;
    long double v = r.x - x;
    r.y = (x - (r.x - v)) + (y - v);
  
    return r;
  }
  
  static INLINE CONST Sleef_longdouble2 dladd_l2_l2_l(Sleef_longdouble2 x, long double y) {
    // |x| >= |y|
  
    Sleef_longdouble2 r;
  
  #ifndef NDEBUG
    if (!(checkfp(x.x) || checkfp(y) || xfabsl(x.x) >= xfabsl(y))) {
      fprintf(stderr, "[dladd_l2_l2_l : %Lg %Lg]\n", x.x, y);
      fflush(stderr);
    }
  #endif
  
    r.x = x.x + y;
    r.y = x.x - r.x + y + x.y;
  
    return r;
  }
  
  static INLINE CONST Sleef_longdouble2 dladd2_l2_l2_l(Sleef_longdouble2 x, long double y) {
    // |x| >= |y|
  
    Sleef_longdouble2 r;
  
    r.x  = x.x + y;
    long double v = r.x - x.x;
    r.y = (x.x - (r.x - v)) + (y - v);
    r.y += x.y;
  
    return r;
  }
  
  static INLINE CONST Sleef_longdouble2 dladd_l2_l_l2(long double x, Sleef_longdouble2 y) {
    // |x| >= |y|
  
    Sleef_longdouble2 r;
  
  #ifndef NDEBUG
    if (!(checkfp(x) || checkfp(y.x) || xfabsl(x) >= xfabsl(y.x))) {
      fprintf(stderr, "[dladd_l2_l_l2 : %Lg %Lg]\n", x, y.x);
      fflush(stderr);
    }
  #endif
  
    r.x = x + y.x;
    r.y = x - r.x + y.x + y.y;
  
    return r;
  }
  
  static INLINE CONST Sleef_longdouble2 dladd2_l2_l_l2(long double x, Sleef_longdouble2 y) {
    Sleef_longdouble2 r;
  
    r.x  = x + y.x;
    long double v = r.x - x;
    r.y = (x - (r.x - v)) + (y.x - v) + y.y;
  
    return r;
  }
  
  static INLINE CONST Sleef_longdouble2 dladd_l2_l2_l2(Sleef_longdouble2 x, Sleef_longdouble2 y) {
    // |x| >= |y|
  
    Sleef_longdouble2 r;
  
  #ifndef NDEBUG
    if (!(checkfp(x.x) || checkfp(y.x) || xfabsl(x.x) >= xfabsl(y.x))) {
      fprintf(stderr, "[dladd_l2_l2_l2 : %Lg %Lg]\n", x.x, y.x);
      fflush(stderr);
    }
  #endif
  
    r.x = x.x + y.x;
    r.y = x.x - r.x + y.x + x.y + y.y;
  
    return r;
  }
  
  static INLINE CONST Sleef_longdouble2 dladd2_l2_l2_l2(Sleef_longdouble2 x, Sleef_longdouble2 y) {
    Sleef_longdouble2 r;
  
    r.x  = x.x + y.x;
    long double v = r.x - x.x;
    r.y = (x.x - (r.x - v)) + (y.x - v);
    r.y += x.y + y.y;
  
    return r;
  }
  
  static INLINE CONST Sleef_longdouble2 dlsub_l2_l2_l2(Sleef_longdouble2 x, Sleef_longdouble2 y) {
    // |x| >= |y|
  
    Sleef_longdouble2 r;
  
  #ifndef NDEBUG
    if (!(checkfp(x.x) || checkfp(y.x) || xfabsl(x.x) >= xfabsl(y.x))) {
      fprintf(stderr, "[dlsub_l2_l2_l2 : %Lg %Lg]\n", x.x, y.x);
      fflush(stderr);
    }
  #endif
  
    r.x = x.x - y.x;
    r.y = x.x - r.x - y.x + x.y - y.y;
  
    return r;
  }
  
  static INLINE CONST Sleef_longdouble2 dldiv_l2_l2_l2(Sleef_longdouble2 n, Sleef_longdouble2 d) {
    long double t = 1.0 / d.x;
    long double dh  = upperl(d.x), dl  = d.x - dh;
    long double th  = upperl(t  ), tl  = t   - th;
    long double nhh = upperl(n.x), nhl = n.x - nhh;
  
    Sleef_longdouble2 q;
  
    q.x = n.x * t;
  
    long double u = -q.x + nhh * th + nhh * tl + nhl * th + nhl * tl +
      q.x * (1 - dh * th - dh * tl - dl * th - dl * tl);
  
    q.y = t * (n.y - q.x * d.y) + u;
  
    return q;
  }
  
  static INLINE CONST Sleef_longdouble2 dlmul_l2_l_l(long double x, long double y) {
    long double xh = upperl(x), xl = x - xh;
    long double yh = upperl(y), yl = y - yh;
    Sleef_longdouble2 r;
  
    r.x = x * y;
    r.y = xh * yh - r.x + xl * yh + xh * yl + xl * yl;
  
    return r;
  }
  
  static INLINE CONST Sleef_longdouble2 dlmul_l2_l2_l(Sleef_longdouble2 x, long double y) {
    long double xh = upperl(x.x), xl = x.x - xh;
    long double yh = upperl(y  ), yl = y   - yh;
    Sleef_longdouble2 r;
  
    r.x = x.x * y;
    r.y = xh * yh - r.x + xl * yh + xh * yl + xl * yl + x.y * y;
  
    return r;
  }
  
  static INLINE CONST Sleef_longdouble2 dlmul_l2_l2_l2(Sleef_longdouble2 x, Sleef_longdouble2 y) {
    long double xh = upperl(x.x), xl = x.x - xh;
    long double yh = upperl(y.x), yl = y.x - yh;
    Sleef_longdouble2 r;
  
    r.x = x.x * y.x;
    r.y = xh * yh - r.x + xl * yh + xh * yl + xl * yl + x.x * y.y + x.y * y.x;
  
    return r;
  }
  
  static INLINE CONST Sleef_longdouble2 dlsqu_l2_l2(Sleef_longdouble2 x) {
    long double xh = upperl(x.x), xl = x.x - xh;
    Sleef_longdouble2 r;
  
    r.x = x.x * x.x;
    r.y = xh * xh - r.x + (xh + xh) * xl + xl * xl + x.x * (x.y + x.y);
  
    return r;
  }
  
  static INLINE CONST Sleef_longdouble2 dlrec_l2_l(long double d) {
    long double t = 1.0 / d;
    long double dh = upperl(d), dl = d - dh;
    long double th = upperl(t), tl = t - th;
    Sleef_longdouble2 q;
  
    q.x = t;
    q.y = t * (1 - dh * th - dh * tl - dl * th - dl * tl);
  
    return q;
  }
  
  static INLINE CONST Sleef_longdouble2 dlrec_l2_l2(Sleef_longdouble2 d) {
    long double t = 1.0 / d.x;
    long double dh = upperl(d.x), dl = d.x - dh;
    long double th = upperl(t  ), tl = t   - th;
    Sleef_longdouble2 q;
  
    q.x = t;
    q.y = t * (1 - dh * th - dh * tl - dl * th - dl * tl - d.y * t);
  
    return q;
  }
  
  /*
  static INLINE CONST Sleef_longdouble2 dlsqrt_l2_l2(Sleef_longdouble2 d) {
    long double t = sqrt(d.x + d.y);
    return dlscale_l2_l2_l(dlmul_l2_l2_l2(dladd2_l2_l2_l2(d, dlmul_l2_l_l(t, t)), dlrec_l2_l(t)), 0.5);
  }
  */
  
  //
  
  EXPORT CONST Sleef_longdouble2 xsincospil_u05(long double d) {
    long double u, s, t;
    Sleef_longdouble2 r, x, s2;
  
    u = d * 4;
    int64_t q = xceill(u) & ~(int64_t)1;
  
    s = u - (long double)q;
    t = s;
    s = s * s;
    s2 = dlmul_l2_l_l(t, t);
    
    //
  
    u = 4.59265607313529833157632e-17L;
    u = mlal(u, s, -2.04096140520547829627419e-14L);
    u = mlal(u, s, 6.94845264320316515640316e-12L);
    u = mlal(u, s, -1.75724767308629210422023e-09L);
    u = mlal(u, s, 3.13361689037693212744991e-07L);
    u = mlal(u, s, -3.65762041821772284521155e-05L);
    u = mlal(u, s, 0.00249039457019272015784594L);
    x = dladd2_l2_l_l2(u * s, dl(-0.0807455121882807817044873L, -2.40179063154839769223037e-21L));
    x = dladd2_l2_l2_l2(dlmul_l2_l2_l2(s2, x), dl(0.785398163397448309628202L, -1.25420305812534448752181e-20L));
    
    x = dlmul_l2_l2_l(x, t);
    r.x = x.x + x.y;
  
    //
  
    u = -2.00423964577657539380734e-18L;
    u = mlal(u, s, 1.00185574457758689324113e-15L);
    u = mlal(u, s, -3.89807283423502620989528e-13L);
    u = mlal(u, s, 1.15011591257563133685341e-10L);
    u = mlal(u, s, -2.461136950493305818105e-08L);
    u = mlal(u, s, 3.59086044859150791782134e-06L);
    u = mlal(u, s, -0.00032599188692739001335938L);
    x = dladd2_l2_l_l2(u * s, dl(0.0158543442438155008529635L, -6.97556143018517384674258e-22L));
    x = dladd2_l2_l2_l2(dlmul_l2_l2_l2(s2, x), dl(-0.308425137534042456829379L, -9.19882299434302978226668e-21L));
    
    x = dladd2_l2_l2_l(dlmul_l2_l2_l2(x, s2), 1);
    r.y = x.x + x.y;
  
    //
    
    if ((q & 2) != 0) { s = r.y; r.y = r.x; r.x = s; }
    if ((q & 4) != 0) { r.x = -r.x; }
    if (((q+2) & 4) != 0) { r.y = -r.y; }
  
    if (xisinfl(d)) { r.x = r.y = SLEEF_NAN; }
    if (!xisinfl(d) && xfabsl(d) > TRIGRANGEMAX3) { r.x = r.y = 0; }
  
    return r;
  }
  
  EXPORT CONST Sleef_longdouble2 xsincospil_u35(long double d) {
    long double u, s, t;
    Sleef_longdouble2 r;
  
    u = d * 4;
    int64_t q = xceill(u) & ~(int64_t)1;
  
    s = u - (long double)q;
    t = s;
    s = s * s;
    
    //
  
    u = -0.2023275819380976135024e-13L;
    u = mlal(u, s, +0.6948176964255957574946e-11L);
    u = mlal(u, s, -0.1757247450021535880723e-8L);
    u = mlal(u, s, +0.3133616889379195970541e-6L);
    u = mlal(u, s, -0.3657620418215300856408e-4L);
    u = mlal(u, s, +0.2490394570192717262476e-2L);
    u = mlal(u, s, -0.8074551218828078160284e-1L);
    u = mlal(u, s, +0.7853981633974483096282e+0L);
    
    r.x = u * t;
  
    //
  
    u = +0.9933418221428971922705e-15L;
    u = mlal(u, s, -0.3897923064055824005357e-12L);
    u = mlal(u, s, +0.1150115771521792692066e-9L);
    u = mlal(u, s, -0.2461136949725905367314e-7L);
    u = mlal(u, s, +0.3590860448589084195081e-5L);
    u = mlal(u, s, -0.3259918869273895914840e-3L);
    u = mlal(u, s, +0.1585434424381550079706e-1L);
    u = mlal(u, s, -0.3084251375340424568294e+0L);
    u = mlal(u, s, 1.0L);
  
    r.y = u;
  
    //
    
    if ((q & 2) != 0) { s = r.y; r.y = r.x; r.x = s; }
    if ((q & 4) != 0) { r.x = -r.x; }
    if (((q+2) & 4) != 0) { r.y = -r.y; }
  
    if (xisinfl(d)) { r.x = r.y = SLEEF_NAN; }
    if (!xisinfl(d) && xfabsl(d) > TRIGRANGEMAX3) { r.x = r.y = 0; }
  
    return r;
  }
  