src/math/asinl.c

   1 /* origin: FreeBSD /usr/src/lib/msun/src/e_asinl.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 asin.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 asinl(long double x)
  21 {
  22         return asin(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 CLOSETO1(u) (u.i.m>>56 >= 0xf7)
  28 #define CLEARBOTTOM(u) (u.i.m &= -1ULL << 32)
  29 #elif LDBL_MANT_DIG == 113
  30 #define CLOSETO1(u) (u.i.top >= 0xee00)
  31 #define CLEARBOTTOM(u) (u.i.lo = 0)
  32 #endif
  33
  34 long double asinl(long double x)
  35 {
  36         union ldshape u = {x};
  37         long double z, r, s;
  38         uint16_t e = u.i.se & 0x7fff;
  39         int sign = u.i.se >> 15;
  40
  41         if (e >= 0x3fff) {   /* |x| >= 1 or nan */
  42                 /* asin(+-1)=+-pi/2 with inexact */
  43                 if (x == 1 || x == -1)
  44                         return x*pio2_hi + 0x1p-120f;
  45                 return 0/(x-x);
  46         }
  47         if (e < 0x3fff - 1) {  /* |x| < 0.5 */
  48                 if (e < 0x3fff - (LDBL_MANT_DIG+1)/2) {
  49                         /* return x with inexact if x!=0 */
  50                         FORCE_EVAL(x + 0x1p120f);
  51                         return x;
  52                 }
  53                 return x + x*__invtrigl_R(x*x);
  54         }
  55         /* 1 > |x| >= 0.5 */
  56         z = (1.0 - fabsl(x))*0.5;
  57         s = sqrtl(z);
  58         r = __invtrigl_R(z);
  59         if (CLOSETO1(u)) {
  60                 x = pio2_hi - (2*(s+s*r)-pio2_lo);
  61         } else {
  62                 long double f, c;
  63                 u.f = s;
  64                 CLEARBOTTOM(u);
  65                 f = u.f;
  66                 c = (z - f*f)/(s + f);
  67                 x = 0.5*pio2_hi-(2*s*r - (pio2_lo-2*c) - (0.5*pio2_hi-2*f));
  68         }
  69         return sign ? -x : x;
  70 }
  71 #endif