1 /* $NetBSD: fpu_exp.c,v 1.11 2017/01/15 11:56:11 isaki Exp $ */
2
3 /*
4 * Copyright (c) 1995 Ken Nakata
5 * All rights reserved.
6 *
7 * Redistribution and use in source and binary forms, with or without
8 * modification, are permitted provided that the following conditions
9 * are met:
10 * 1. Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 * 2. Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the distribution.
15 * 3. Neither the name of the author nor the names of its contributors
16 * may be used to endorse or promote products derived from this software
17 * without specific prior written permission.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
20 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
23 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
24 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
25 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
26 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
27 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
28 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
29 * SUCH DAMAGE.
30 *
31 * @(#)fpu_exp.c 10/24/95
32 */
33
34 #include <sys/cdefs.h>
35 __KERNEL_RCSID(0, "$NetBSD: fpu_exp.c,v 1.11 2017/01/15 11:56:11 isaki Exp $");
36
37 #include <machine/ieee.h>
38
39 #include "fpu_emulate.h"
40
41 /* The number of items to terminate the Taylor expansion */
42 #define MAX_ITEMS (2000)
43
44 /*
45 * fpu_exp.c: defines fpu_etox(), fpu_etoxm1(), fpu_tentox(), and fpu_twotox();
46 */
47
48 /*
49 * x^2 x^3 x^4
50 * exp(x) = 1 + x + --- + --- + --- + ...
51 * 2! 3! 4!
52 */
53 static struct fpn *
fpu_etox_taylor(struct fpemu * fe)54 fpu_etox_taylor(struct fpemu *fe)
55 {
56 struct fpn res;
57 struct fpn x;
58 struct fpn s0;
59 struct fpn *s1;
60 struct fpn *r;
61 uint32_t k;
62
63 CPYFPN(&x, &fe->fe_f2);
64 CPYFPN(&s0, &fe->fe_f2);
65
66 /* res := 1 + x */
67 fpu_const(&fe->fe_f1, FPU_CONST_1);
68 r = fpu_add(fe);
69 CPYFPN(&res, r);
70
71 k = 2;
72 for (; k < MAX_ITEMS; k++) {
73 /* s1 = s0 * x / k */
74 CPYFPN(&fe->fe_f1, &s0);
75 CPYFPN(&fe->fe_f2, &x);
76 r = fpu_mul(fe);
77
78 CPYFPN(&fe->fe_f1, r);
79 fpu_explode(fe, &fe->fe_f2, FTYPE_LNG, &k);
80 s1 = fpu_div(fe);
81
82 /* break if s1 is enough small */
83 if (ISZERO(s1))
84 break;
85 if (res.fp_exp - s1->fp_exp >= EXT_FRACBITS)
86 break;
87
88 /* s0 := s1 for next loop */
89 CPYFPN(&s0, s1);
90
91 /* res += s1 */
92 CPYFPN(&fe->fe_f2, s1);
93 CPYFPN(&fe->fe_f1, &res);
94 r = fpu_add(fe);
95 CPYFPN(&res, r);
96 }
97
98 CPYFPN(&fe->fe_f2, &res);
99 return &fe->fe_f2;
100 }
101
102 /*
103 * exp(x) = 2^k * exp(r) with k = round(x / ln2) and r = x - k * ln2
104 *
105 * Algorithm partially taken from libm, where exp(r) is approximated by a
106 * rational function of r. We use the Taylor expansion instead.
107 */
108 struct fpn *
fpu_etox(struct fpemu * fe)109 fpu_etox(struct fpemu *fe)
110 {
111 struct fpn x, *fp;
112 int k;
113
114 if (ISNAN(&fe->fe_f2))
115 return &fe->fe_f2;
116 if (ISINF(&fe->fe_f2)) {
117 if (fe->fe_f2.fp_sign)
118 fpu_const(&fe->fe_f2, FPU_CONST_0);
119 return &fe->fe_f2;
120 }
121
122 /*
123 * return inf if x >= 2^14
124 * return +0 if x <= -2^14
125 */
126 if (fe->fe_f2.fp_exp >= 14) {
127 if (fe->fe_f2.fp_sign) {
128 fe->fe_f2.fp_class = FPC_ZERO;
129 fe->fe_f2.fp_sign = 0;
130 } else {
131 fe->fe_f2.fp_class = FPC_INF;
132 }
133 return &fe->fe_f2;
134 }
135
136 CPYFPN(&x, &fe->fe_f2);
137
138 /* k = round(x / ln2) */
139 CPYFPN(&fe->fe_f1, &fe->fe_f2);
140 fpu_const(&fe->fe_f2, FPU_CONST_LN_2);
141 fp = fpu_div(fe);
142 CPYFPN(&fe->fe_f2, fp);
143 fp = fpu_int(fe);
144 if (ISZERO(fp)) {
145 /* k = 0 */
146 CPYFPN(&fe->fe_f2, &x);
147 fp = fpu_etox_taylor(fe);
148 return fp;
149 }
150 /* extract k as integer format from fpn format */
151 k = fp->fp_mant[0] >> (FP_LG - fp->fp_exp);
152 if (fp->fp_sign)
153 k *= -1;
154
155 /* exp(r) = exp(x - k * ln2) */
156 CPYFPN(&fe->fe_f1, fp);
157 fpu_const(&fe->fe_f2, FPU_CONST_LN_2);
158 fp = fpu_mul(fe);
159 fp->fp_sign = !fp->fp_sign;
160 CPYFPN(&fe->fe_f1, fp);
161 CPYFPN(&fe->fe_f2, &x);
162 fp = fpu_add(fe);
163 CPYFPN(&fe->fe_f2, fp);
164 fp = fpu_etox_taylor(fe);
165
166 /* 2^k */
167 fp->fp_exp += k;
168
169 return fp;
170 }
171
172 /*
173 * exp(x) - 1
174 */
175 struct fpn *
fpu_etoxm1(struct fpemu * fe)176 fpu_etoxm1(struct fpemu *fe)
177 {
178 struct fpn *fp;
179
180 fp = fpu_etox(fe);
181
182 CPYFPN(&fe->fe_f1, fp);
183 /* build a 1.0 */
184 fp = fpu_const(&fe->fe_f2, FPU_CONST_1);
185 fe->fe_f2.fp_sign = !fe->fe_f2.fp_sign;
186 /* fp = f2 - 1.0 */
187 fp = fpu_add(fe);
188
189 return fp;
190 }
191
192 /*
193 * 10^x = exp(x * ln10)
194 */
195 struct fpn *
fpu_tentox(struct fpemu * fe)196 fpu_tentox(struct fpemu *fe)
197 {
198 struct fpn *fp;
199
200 /* build a ln10 */
201 fp = fpu_const(&fe->fe_f1, FPU_CONST_LN_10);
202 /* fp = ln10 * f2 */
203 fp = fpu_mul(fe);
204
205 /* copy the result to the src opr */
206 CPYFPN(&fe->fe_f2, fp);
207
208 return fpu_etox(fe);
209 }
210
211 /*
212 * 2^x = exp(x * ln2)
213 */
214 struct fpn *
fpu_twotox(struct fpemu * fe)215 fpu_twotox(struct fpemu *fe)
216 {
217 struct fpn *fp;
218
219 /* build a ln2 */
220 fp = fpu_const(&fe->fe_f1, FPU_CONST_LN_2);
221 /* fp = ln2 * f2 */
222 fp = fpu_mul(fe);
223
224 /* copy the result to the src opr */
225 CPYFPN(&fe->fe_f2, fp);
226
227 return fpu_etox(fe);
228 }
229