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 /* 0.95 */
  27 #define THRESH  ((0xe666666666666666ULL>>(64-(LDBL_MANH_SIZE-1)))|LDBL_NBIT)
  28
  29 long double asinl(long double x)
  30 {
  31         union IEEEl2bits u;
  32         long double z,r,s;
  33         uint16_t expsign, expt;
  34
  35         u.e = x;
  36         expsign = u.xbits.expsign;
  37         expt = expsign & 0x7fff;
  38         if (expt >= 0x3fff) {   /* |x| >= 1 or nan */
  39                 if (expt == 0x3fff &&
  40                     ((u.bits.manh&~LDBL_NBIT)|u.bits.manl) == 0)
  41                         /* asin(+-1)=+-pi/2 with inexact */
  42                         return x*pio2_hi + 0x1p-120f;
  43                 return 0/(x-x);
  44         }
  45         if (expt < 0x3fff - 1) {  /* |x| < 0.5 */
  46                 if (expt < 0x3fff - 32) {  /* |x|<0x1p-32, asinl(x)=x */
  47                         /* return x with inexact if x!=0 */
  48                         FORCE_EVAL(x + 0x1p120f);
  49                         return x;
  50                 }
  51                 return x + x*__invtrigl_R(x*x);
  52         }
  53         /* 1 > |x| >= 0.5 */
  54         z = (1.0 - fabsl(x))*0.5;
  55         s = sqrtl(z);
  56         r = __invtrigl_R(z);
  57         if (u.bits.manh >= THRESH) { /* if |x| is close to 1 */
  58                 x = pio2_hi - (2*(s+s*r)-pio2_lo);
  59         } else {
  60                 long double f, c;
  61                 u.e = s;
  62                 u.bits.manl = 0;
  63                 f = u.e;
  64                 c = (z-f*f)/(s+f);
  65                 x = 0.5*pio2_hi-(2*s*r - (pio2_lo-2*c) - (0.5*pio2_hi-2*f));
  66         }
  67         if (expsign>>15)
  68                 return -x;
  69         return x;
  70 }
  71 #endif