contrib/openbsd_libm/src/e_acosl.c

   1 /*      $OpenBSD: e_acosl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $    */
   2 /* @(#)e_acos.c 1.3 95/01/18 */
   3 /* FreeBSD: head/lib/msun/src/e_acos.c 176451 2008-02-22 02:30:36Z das */
   4 /*
   5  * ====================================================
   6  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   7  *
   8  * Developed at SunSoft, a Sun Microsystems, Inc. business.
   9  * Permission to use, copy, modify, and distribute this
  10  * software is freely granted, provided that this notice
  11  * is preserved.
  12  * ====================================================
  13  */
  14
  15 /*
  16  * See comments in e_acos.c.
  17  * Converted to long double by David Schultz <das@FreeBSD.ORG>.
  18  * Adapted for OpenBSD by Martynas Venckus <martynas@openbsd.org>.
  19  */
  20
  21 #include <float.h>
  22 #include <math.h>
  23
  24 #include "invtrig.h"
  25 #include "math_private.h"
  26
  27 #ifdef EXT_IMPLICIT_NBIT
  28 #define LDBL_NBIT       0
  29 #else /* EXT_IMPLICIT_NBIT */
  30 #define LDBL_NBIT       0x80000000
  31 #endif /* EXT_IMPLICIT_NBIT */
  32
  33 static const long double
  34 one=  1.00000000000000000000e+00;
  35
  36 #if defined(__amd64__) || defined(__i386__)
  37 /* XXX Work around the fact that gcc truncates long double constants on i386 */
  38 static volatile double
  39 pi1 =  3.14159265358979311600e+00,      /*  0x1.921fb54442d18p+1  */
  40 pi2 =  1.22514845490862001043e-16;      /*  0x1.1a80000000000p-53 */
  41 #define pi      ((long double)pi1 + pi2)
  42 #else
  43 static const long double
  44 pi =  3.14159265358979323846264338327950280e+00L;
  45 #endif
  46
  47 long double
  48 acosl(long double x)
  49 {
  50         union {
  51                 long double e;
  52                 struct ieee_ext bits;
  53         } u;
  54         long double z,p,q,r,w,s,c,df;
  55         int16_t expsign, expt;
  56         u.e = x;
  57         expsign = (u.bits.ext_sign << 15) | u.bits.ext_exp;
  58         expt = expsign & 0x7fff;
  59         if(expt >= BIAS) {      /* |x| >= 1 */
  60             if(expt==BIAS && ((u.bits.ext_frach&~LDBL_NBIT)
  61 #ifdef EXT_FRACHMBITS
  62                 | u.bits.ext_frachm
  63 #endif /* EXT_FRACHMBITS */
  64 #ifdef EXT_FRACLMBITS
  65                 | u.bits.ext_fraclm
  66 #endif /* EXT_FRACLMBITS */
  67                 | u.bits.ext_fracl)==0) {
  68                 if (expsign>0) return 0.0;      /* acos(1) = 0  */
  69                 else return pi+2.0*pio2_lo;     /* acos(-1)= pi */
  70             }
  71             return (x-x)/(x-x);         /* acos(|x|>1) is NaN */
  72         }
  73         if(expt<BIAS-1) {       /* |x| < 0.5 */
  74             if(expt<ACOS_CONST) return pio2_hi+pio2_lo;/*x tiny: acosl=pi/2*/
  75             z = x*x;
  76             p = P(z);
  77             q = Q(z);
  78             r = p/q;
  79             return pio2_hi - (x - (pio2_lo-x*r));
  80         } else  if (expsign<0) {        /* x < -0.5 */
  81             z = (one+x)*0.5;
  82             p = P(z);
  83             q = Q(z);
  84             s = sqrtl(z);
  85             r = p/q;
  86             w = r*s-pio2_lo;
  87             return pi - 2.0*(s+w);
  88         } else {                        /* x > 0.5 */
  89             z = (one-x)*0.5;
  90             s = sqrtl(z);
  91             u.e = s;
  92             u.bits.ext_fracl = 0;
  93 #ifdef EXT_FRACLMBITS
  94             u.bits.ext_fraclm = 0;
  95 #endif /* EXT_FRACLMBITS */
  96             df = u.e;
  97             c  = (z-df*df)/(s+df);
  98             p = P(z);
  99             q = Q(z);
 100             r = p/q;
 101             w = r*s+c;
 102             return 2.0*(df+w);
 103         }
 104 }