14d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle/* from: FreeBSD: head/lib/msun/src/e_acosh.c 176451 2008-02-22 02:30:36Z das */ 24d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle 34d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle/* @(#)e_acosh.c 1.3 95/01/18 */ 44d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle/* 54d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle * ==================================================== 64d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 74d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle * 84d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle * Developed at SunSoft, a Sun Microsystems, Inc. business. 94d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle * Permission to use, copy, modify, and distribute this 104d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle * software is freely granted, provided that this notice 114d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle * is preserved. 124d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle * ==================================================== 134d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle * 144d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle */ 154d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle 164d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#include <sys/cdefs.h> 174d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle__FBSDID("$FreeBSD$"); 184d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle 194d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle/* 204d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle * See e_acosh.c for complete comments. 214d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle * 224d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle * Converted to long double by David Schultz <das@FreeBSD.ORG> and 234d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle * Bruce D. Evans. 244d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle */ 254d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle 264d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#include <float.h> 274d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#ifdef __i386__ 284d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#include <ieeefp.h> 294d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#endif 304d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle 314d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#include "fpmath.h" 324d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#include "math.h" 334d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#include "math_private.h" 344d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle 354d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle/* EXP_LARGE is the threshold above which we use acosh(x) ~= log(2x). */ 364d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#if LDBL_MANT_DIG == 64 374d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#define EXP_LARGE 34 384d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#elif LDBL_MANT_DIG == 113 394d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#define EXP_LARGE 58 404d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#else 414d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#error "Unsupported long double format" 424d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#endif 434d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle 444d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#if LDBL_MAX_EXP != 0x4000 454d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle/* We also require the usual expsign encoding. */ 464d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#error "Unsupported long double format" 474d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#endif 484d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle 494d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#define BIAS (LDBL_MAX_EXP - 1) 504d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle 514d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravlestatic const double 524d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravleone = 1.0; 534d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle 544d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#if LDBL_MANT_DIG == 64 554d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravlestatic const union IEEEl2bits 564d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravleu_ln2 = LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L); 574d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#define ln2 u_ln2.e 584d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#elif LDBL_MANT_DIG == 113 594d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravlestatic const long double 604d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravleln2 = 6.93147180559945309417232121458176568e-1L; /* 0x162e42fefa39ef35793c7673007e6.0p-113 */ 614d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#else 624d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#error "Unsupported long double format" 634d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#endif 644d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle 654d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravlelong double 664d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravleacoshl(long double x) 674d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle{ 684d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle long double t; 694d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle int16_t hx; 704d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle 714d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle ENTERI(); 724d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle GET_LDBL_EXPSIGN(hx, x); 734d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle if (hx < 0x3fff) { /* x < 1, or misnormal */ 744d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle RETURNI((x-x)/(x-x)); 754d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle } else if (hx >= BIAS + EXP_LARGE) { /* x >= LARGE */ 764d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle if (hx >= 0x7fff) { /* x is inf, NaN or misnormal */ 774d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle RETURNI(x+x); 784d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle } else 794d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle RETURNI(logl(x)+ln2); /* acosh(huge)=log(2x), or misnormal */ 804d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle } else if (hx == 0x3fff && x == 1) { 814d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle RETURNI(0.0); /* acosh(1) = 0 */ 824d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle } else if (hx >= 0x4000) { /* LARGE > x >= 2, or misnormal */ 834d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle t=x*x; 844d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle RETURNI(logl(2.0*x-one/(x+sqrtl(t-one)))); 854d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle } else { /* 1<x<2 */ 864d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle t = x-one; 874d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle RETURNI(log1pl(t+sqrtl(2.0*t+t*t))); 884d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle } 894d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle} 90