nsz Git - musl/blob - src/math/acosl.c

   1 /* origin: FreeBSD /usr/src/lib/msun/src/e_acosl.c */
   2 /*
   3  * ====================================================
   4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   5  *
   6  * Developed at SunSoft, a Sun Microsystems, Inc. business.
   7  * Permission to use, copy, modify, and distribute this
   8  * software is freely granted, provided that this notice
   9  * is preserved.
  10  * ====================================================
  11  */
  12 /*
  13  * See comments in acos.c.
  14  * Converted to long double by David Schultz <das@FreeBSD.ORG>.
  15  */
  16
  17 #include "libm.h"
  18
  19 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
  20 long double acosl(long double x)
  21 {
  22         return acos(x);
  23 }
  24 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
  25 #include "__invtrigl.h"
  26 #if LDBL_MANT_DIG == 64
  27 #define CLEARBOTTOM(u) (u.i.m &= -1ULL << 32)
  28 #elif LDBL_MANT_DIG == 113
  29 #define CLEARBOTTOM(u) (u.i.lo = 0)
  30 #endif
  31
  32 long double acosl(long double x)
  33 {
  34         union ldshape u = {x};
  35         long double z, s, c, f;
  36         uint16_t e = u.i.se & 0x7fff;
  37
  38         /* |x| >= 1 or nan */
  39         if (e >= 0x3fff) {
  40                 if (x == 1)
  41                         return 0;
  42                 if (x == -1)
  43                         return 2*pio2_hi + 0x1p-120f;
  44                 return 0/(x-x);
  45         }
  46         /* |x| < 0.5 */
  47         if (e < 0x3fff - 1) {
  48                 if (e < 0x3fff - LDBL_MANT_DIG - 1)
  49                         return pio2_hi + 0x1p-120f;
  50                 return pio2_hi - (__invtrigl_R(x*x)*x - pio2_lo + x);
  51         }
  52         /* x < -0.5 */
  53         if (u.i.se >> 15) {
  54                 z = (1 + x)*0.5;
  55                 s = sqrtl(z);
  56                 return 2*(pio2_hi - (__invtrigl_R(z)*s - pio2_lo + s));
  57         }
  58         /* x > 0.5 */
  59         z = (1 - x)*0.5;
  60         s = sqrtl(z);
  61         u.f = s;
  62         CLEARBOTTOM(u);
  63         f = u.f;
  64         c = (z - f*f)/(s + f);
  65         return 2*(__invtrigl_R(z)*s + c + f);
  66 }
  67 #endif