Home | History | Annotate | Line # | Download | only in math 
      1 /* e_acoshl.c -- long double version of e_acosh.c.
      2  * Conversion to long double by Jakub Jelinek, jj (at) ultra.linux.cz.
      3  */
      4 
      5 /*
      6  * ====================================================
      7  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
      8  *
      9  * Developed at SunPro, a Sun Microsystems, Inc. business.
     10  * Permission to use, copy, modify, and distribute this
     11  * software is freely granted, provided that this notice
     12  * is preserved.
     13  * ====================================================
     14  */
     15 
     16 /* acoshq(x)
     17  * Method :
     18  *	Based on
     19  *		acoshl(x) = logq [ x + sqrtq(x*x-1) ]
     20  *	we have
     21  *		acoshl(x) := logq(x)+ln2,	if x is large; else
     22  *		acoshl(x) := logq(2x-1/(sqrtq(x*x-1)+x)) if x>2; else
     23  *		acoshl(x) := log1pq(t+sqrtq(2.0*t+t*t)); where t=x-1.
     24  *
     25  * Special cases:
     26  *	acoshl(x) is NaN with signal if x<1.
     27  *	acoshl(NaN) is NaN without signal.
     28  */
     29 
     30 #include "quadmath-imp.h"
     31 
     32 static const __float128
     33 one	= 1.0,
     34 ln2	= 0.6931471805599453094172321214581766Q;
     35 
     36 __float128
     37 acoshq(__float128 x)
     38 {
     39 	__float128 t;
     40 	uint64_t lx;
     41 	int64_t hx;
     42 	GET_FLT128_WORDS64(hx,lx,x);
     43 	if(hx<0x3fff000000000000LL) {		/* x < 1 */
     44 	    return (x-x)/(x-x);
     45 	} else if(hx >=0x4035000000000000LL) {	/* x > 2**54 */
     46 	    if(hx >=0x7fff000000000000LL) {	/* x is inf of NaN */
     47 		return x+x;
     48 	    } else
     49 		return logq(x)+ln2;	/* acoshl(huge)=logq(2x) */
     50 	} else if(((hx-0x3fff000000000000LL)|lx)==0) {
     51 	    return 0;			/* acosh(1) = 0 */
     52 	} else if (hx > 0x4000000000000000LL) {	/* 2**28 > x > 2 */
     53 	    t=x*x;
     54 	    return logq(2*x-one/(x+sqrtq(t-one)));
     55 	} else {			/* 1<x<2 */
     56 	    t = x-one;
     57 	    return log1pq(t+sqrtq(2*t+t*t));
     58 	}
     59 }
     60