14d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle/* from: FreeBSD: head/lib/msun/src/e_acosh.c 176451 2008-02-22 02:30:36Z das */ 24d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle 34d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle/* @(#)s_asinh.c 5.1 93/09/24 */ 44d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle/* 54d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle * ==================================================== 64d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 74d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle * 84d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle * Developed at SunPro, 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#include <sys/cdefs.h> 164d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle__FBSDID("$FreeBSD$"); 174d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle 184d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle/* 194d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle * See s_asinh.c for complete comments. 204d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle * 214d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle * Converted to long double by David Schultz <das@FreeBSD.ORG> and 224d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle * Bruce D. Evans. 234d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle */ 244d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle 254d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#include <float.h> 264d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#ifdef __i386__ 274d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#include <ieeefp.h> 284d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#endif 294d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle 304d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#include "fpmath.h" 314d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#include "math.h" 324d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#include "math_private.h" 334d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle 344d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle/* EXP_LARGE is the threshold above which we use asinh(x) ~= log(2x). */ 354d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle/* EXP_TINY is the threshold below which we use asinh(x) ~= x. */ 364d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#if LDBL_MANT_DIG == 64 374d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#define EXP_LARGE 34 384d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#define EXP_TINY -34 394d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#elif LDBL_MANT_DIG == 113 404d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#define EXP_LARGE 58 414d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#define EXP_TINY -58 424d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#else 434d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#error "Unsupported long double format" 444d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#endif 454d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle 464d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#if LDBL_MAX_EXP != 0x4000 474d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle/* We also require the usual expsign encoding. */ 484d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#error "Unsupported long double format" 494d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#endif 504d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle 514d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#define BIAS (LDBL_MAX_EXP - 1) 524d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle 534d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravlestatic const double 544d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravleone = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ 554d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravlehuge= 1.00000000000000000000e+300; 564d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle 574d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#if LDBL_MANT_DIG == 64 584d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravlestatic const union IEEEl2bits 594d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravleu_ln2 = LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L); 604d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#define ln2 u_ln2.e 614d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#elif LDBL_MANT_DIG == 113 624d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravlestatic const long double 634d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravleln2 = 6.93147180559945309417232121458176568e-1L; /* 0x162e42fefa39ef35793c7673007e6.0p-113 */ 644d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#else 654d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#error "Unsupported long double format" 664d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle#endif 674d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle 684d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravlelong double 694d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravleasinhl(long double x) 704d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle{ 714d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle long double t, w; 724d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle uint16_t hx, ix; 734d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle 744d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle ENTERI(); 754d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle GET_LDBL_EXPSIGN(hx, x); 764d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle ix = hx & 0x7fff; 774d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle if (ix >= 0x7fff) RETURNI(x+x); /* x is inf, NaN or misnormal */ 784d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle if (ix < BIAS + EXP_TINY) { /* |x| < TINY, or misnormal */ 794d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle if (huge + x > one) RETURNI(x); /* return x inexact except 0 */ 804d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle } 814d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle if (ix >= BIAS + EXP_LARGE) { /* |x| >= LARGE, or misnormal */ 824d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle w = logl(fabsl(x))+ln2; 834d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle } else if (ix >= 0x4000) { /* LARGE > |x| >= 2.0, or misnormal */ 844d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle t = fabsl(x); 854d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle w = logl(2.0*t+one/(sqrtl(x*x+one)+t)); 864d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle } else { /* 2.0 > |x| >= TINY, or misnormal */ 874d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle t = x*x; 884d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle w =log1pl(fabsl(x)+t/(one+sqrtl(one+t))); 894d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle } 904d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle RETURNI((hx & 0x8000) == 0 ? w : -w); 914d77c1151c40010d137e4a2fa8629bff4bea72b0Calin Juravle} 92