overby.me
this repo has no description
1
2//
3// asinh
4//
5// by Ian Ollmann
6//
7// Based on algorithms from MathLib v3
8//
9// Copyright � 2005, Apple Computer Inc. All Rights Reserved.
10//
11
12#include <math.h>
13
14double asinh( double x )
15{
16 static const double ln2 = 0x1.62e42fefa39ef358p-1; //ln(2)
17
18 if( x != x ) return x + x;
19
20 long double fabsx = __builtin_fabs( x );
21
22 if( fabsx < 0x1.0p-27 ) //sqrt( negative epsilon )
23 {
24 if( x == 0.0L )
25 return x;
26
27 fabsx *= 0x1.0p55;
28 fabsx -= 0x1.0p-1022;
29 fabsx *= 0x1.0p-55;
30 }
31 else if( fabsx <= 4.0 / 3.0 )
32 {
33 double r = 1.0 / fabsx;
34
35 fabsx = log1p( fabsx + fabsx / ( r + sqrt( 1 + r * r)) );
36 }
37 else if( fabsx <= 0x1.0p27 ) //1/sqrt( negative epsilon )
38 {
39 fabsx = log( fabsx + fabsx + 1.0 / (fabsx + sqrt( 1.0 + fabsx * fabsx)) );
40 }
41 else
42 fabsx = log( fabsx ) + ln2;
43
44 if( x < 0 )
45 fabsx = -fabsx;
46
47 return fabsx;
48}
49
50
51long double asinhl( long double x )
52{
53 static const long double ln2 = 0x1.62e42fefa39ef358p-1L; //ln(2)
54
55 if( x != x ) return x + x;
56
57 long double fabsx = __builtin_fabsl( x );
58
59 if( fabsx < 0x1.0p-32 ) //sqrt( negative epsilon )
60 {
61 if( x == 0.0L )
62 return x;
63
64 fabsx *= 0x1.0p65L;
65 fabsx -= 0x1.0p-16382L;
66 fabsx *= 0x1.0p-65L;
67 }
68 else if( fabsx <= 4.0L / 3.0L )
69 {
70 long double r = 1.0L / fabsx;
71
72 fabsx = log1pl( fabsx + fabsx / ( r + sqrtl( 1 + r * r)) );
73 }
74 else if( fabsx <= 0x1.0p32 ) //1/sqrt( negative epsilon )
75 {
76 fabsx = logl( fabsx + fabsx + 1.0L / (fabsx + sqrtl( 1.0L + fabsx * fabsx)) );
77 }
78 else
79 fabsx = logl( fabsx ) + ln2;
80
81 if( x < 0 )
82 fabsx = -fabsx;
83
84 return fabsx;
85}