compiler/stdc/math/s_sincosl.c

   1 /* s_sincosl.c -- long double version of s_sincos.c
   2  *
   3  * Copyright (C) 2013 Elliot Saba
   4  * Developed at the University of Washington
   5  *
   6  * Permission to use, copy, modify, and distribute this
   7  * software is freely granted, provided that this notice
   8  * is preserved.
   9  * ====================================================
  10 */
  11
  12 #include <aros/debug.h>
  13
  14 #include <float.h>
  15 #include "math.h"
  16 #include "math_private.h"
  17 #include "k_sincosl.h"
  18
  19 #if defined(__x86_64__)
  20 #include "x86_64/ieeefp.h"
  21 #elif defined(__i386__)
  22 #include "i386/ieeefp.h"
  23 #endif
  24
  25 #include "math_private.h"
  26 #if LDBL_MANT_DIG == 64
  27 #include "../ld80/e_rem_pio2l.h"
  28 #elif LDBL_MANT_DIG == 113
  29 #include "../ld128/e_rem_pio2l.h"
  30 #else
  31 #error "Unsupported long double format"
  32 #endif
  33
  34 void
  35 sincosl(long double x, long double *sn, long double *cs)
  36 {
  37         union IEEEl2bits z;
  38         int e0;
  39         __unused int sgn;
  40         long double y[2];
  41
  42         z.e = x;
  43         sgn = z.bits.sign;
  44         z.bits.sign = 0;
  45
  46         ENTERV();
  47
  48         /* Optimize the case where x is already within range. */
  49         if (z.e < M_PI_4) {
  50                 /*
  51                  * If x = +-0 or x is a subnormal number, then sin(x) = x and
  52                  * cos(x) = 1.
  53                  */
  54                 if (z.bits.exp == 0) {
  55                         *sn = x;
  56                         *cs = 1;
  57                 } else
  58                         __kernel_sincosl(x, 0, 0, sn, cs);
  59                 RETURNV();
  60         }
  61
  62         /* If x = NaN or Inf, then sin(x) and cos(x) are NaN. */
  63         if (z.bits.exp == 32767) {
  64                 *sn = x - x;
  65                 *cs = x - x;
  66                 RETURNV();
  67         }
  68
  69         /* Range reduction. */
  70         e0 = __ieee754_rem_pio2l(x, y);
  71
  72         switch (e0 & 3) {
  73         case 0:
  74                 __kernel_sincosl(y[0], y[1], 1, sn, cs);
  75                 break;
  76         case 1:
  77                 __kernel_sincosl(y[0], y[1], 1, cs, sn);
  78                 *cs = -*cs;
  79                 break;
  80         case 2:
  81                 __kernel_sincosl(y[0], y[1], 1, sn, cs);
  82                 *sn = -*sn;
  83                 *cs = -*cs;
  84                 break;
  85         default:
  86                 __kernel_sincosl(y[0], y[1], 1, cs, sn);
  87                 *sn = -*sn;
  88         }
  89
  90         RETURNV();
  91 }