github.com/afumu/libc@v0.0.6/musl/src/math/asinl.c (about)

     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