logh.c
2.21 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
// logh - software implementation of function "half logh(half x)"
#include <stdint.h>
#include "extmath.h"
//#if !(SOFT_FOPS_HALF & EXTMATH_FPOP_MASK_CF2I)
// #error NO HW FSQRTH implementation
//#endif
// software implementation based on glibc/openlibm s_logf.c
static const half
ln2_hi = 6.9313812256e-01H, /* 0x3f317180 */
ln2_lo = 9.0580006145e-06H, /* 0x3717f7d1 */
two12 = 4096.0H, /* 0x6C00 */ /* float: two25 = 3.355443200e+07, -> 0x4c000000 */
/* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */
Lg1 = 0xaaaaaa.0p-24H, /* 0.66666662693 */
Lg2 = 0xccce13.0p-25H, /* 0.40000972152 */
Lg3 = 0x91e9ee.0p-25H, /* 0.28498786688 */
Lg4 = 0xf89e26.0p-26H; /* 0.24279078841 */
static const half
zero = 0.0H;
half logh(half x)
{
half hfsq,f,s,z,R,w,t1,t2,dk;
int16_t k,ix,i,j;
GET_HALF_WORD(ix,x);
k=0;
if (ix < 0x0400) { /* x < 2**-15 */
if ((ix&0x7fff)==0)
return -two12/zero; /* log(+-0)=-inf */
if (ix<0)
return (x-x)/zero; /* log(-#) = NaN */
k -= 12; x *= two12; /* subnormal number, scale up x */
GET_HALF_WORD(ix,x);
}
if (ix >= 0x7c00)
return x+x;
k += (ix>>10)-15;
ix &= 0x03ff;
i = (ix+0x0257)&0x0400; // float: i = (ix+(0x95f64<<3))&0x00800000;
SET_HALF_WORD(x,ix|(i^0x3C00)); /* normalize x or x/2 */
k += (i>>10);
f = x-(half)1.0H;
if((0x03ff&(0x8000+ix))<0xc000) { /* -2**-9 <= f < 2**-9 */
if(f==zero) {
if(k==0) {
return zero;
} else {
dk=(half)k;
return dk*ln2_hi+dk*ln2_lo;
}
}
R = f*f*((half)0.5H - (half) 0.33333333333333333H * f);
if(k==0)
return f-R;
else {
dk=(half)k;
return dk*ln2_hi-((R-dk*ln2_lo)-f);
}
}
s = f/((half)2.0H+f);
dk = (half)k;
z = s*s;
i = ix-0x30A; /* float: i = ix-(0x6147a<<3); */
w = z*z;
j = 0x35C-ix; /* float: j = (0x6b851<<3)-ix; */
t1= w*(Lg2+w*Lg4);
t2= z*(Lg1+w*Lg3);
i |= j;
R = t2+t1;
if(i>0) {
hfsq=(half)0.5H * f*f;
if(k==0)
return f-(hfsq-s*(hfsq+R));
else
return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
} else {
if(k==0)
return f-s*(f-R);
else
return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
}
}