![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to fpu_exp.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Isaki's NoNo m68k/m88k emulator] / nono / fpe|
Return to fpu_exp.c CVS log | Up to [Isaki's NoNo m68k/m88k emulator] / nono / fpe |
1.1.1.2 ! root 1: /* $NetBSD: fpu_exp.c,v 1.11 2017/01/15 11:56:11 isaki Exp $ */
1.1 root 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 "fpu_emulate.h"
35:
36: /* The number of items to terminate the Taylor expansion */
37: #define MAX_ITEMS (2000)
38:
39: /*
40: * fpu_exp.c: defines fpu_etox(), fpu_etoxm1(), fpu_tentox(), and fpu_twotox();
41: */
42:
43: /*
44: * x^2 x^3 x^4
45: * exp(x) = 1 + x + --- + --- + --- + ...
46: * 2! 3! 4!
47: */
48: static struct fpn *
49: fpu_etox_taylor(struct fpemu *fe)
50: {
51: struct fpn res;
52: struct fpn x;
53: struct fpn s0;
54: struct fpn *s1;
55: struct fpn *r;
56: uint32_t k;
57:
58: CPYFPN(&x, &fe->fe_f2);
59: CPYFPN(&s0, &fe->fe_f2);
60:
61: /* res := 1 + x */
62: fpu_const(&fe->fe_f1, FPU_CONST_1);
63: r = fpu_add(fe);
64: CPYFPN(&res, r);
65:
66: k = 2;
67: for (; k < MAX_ITEMS; k++) {
68: /* s1 = s0 * x / k */
69: CPYFPN(&fe->fe_f1, &s0);
70: CPYFPN(&fe->fe_f2, &x);
71: r = fpu_mul(fe);
72:
73: CPYFPN(&fe->fe_f1, r);
74: fpu_explode(fe, &fe->fe_f2, FTYPE_LNG, &k);
75: s1 = fpu_div(fe);
76:
77: /* break if s1 is enough small */
78: if (ISZERO(s1))
79: break;
80: if (res.fp_exp - s1->fp_exp >= EXT_FRACBITS)
81: break;
82:
83: /* s0 := s1 for next loop */
84: CPYFPN(&s0, s1);
85:
86: /* res += s1 */
87: CPYFPN(&fe->fe_f2, s1);
88: CPYFPN(&fe->fe_f1, &res);
89: r = fpu_add(fe);
90: CPYFPN(&res, r);
91: }
92:
93: CPYFPN(&fe->fe_f2, &res);
94: return &fe->fe_f2;
95: }
96:
97: /*
98: * exp(x) = 2^k * exp(r) with k = round(x / ln2) and r = x - k * ln2
99: *
100: * Algorithm partially taken from libm, where exp(r) is approximated by a
101: * rational function of r. We use the Taylor expansion instead.
102: */
103: struct fpn *
104: fpu_etox(struct fpemu *fe)
105: {
106: struct fpn x, *fp;
107: int k;
108:
109: if (ISNAN(&fe->fe_f2))
110: return &fe->fe_f2;
111: if (ISINF(&fe->fe_f2)) {
112: if (fe->fe_f2.fp_sign)
113: fpu_const(&fe->fe_f2, FPU_CONST_0);
114: return &fe->fe_f2;
115: }
116: if (ISZERO(&fe->fe_f2)) {
117: fpu_const(&fe->fe_f2, FPU_CONST_1);
118: return &fe->fe_f2;
119: }
120:
121: /*
122: * return inf if x >= 2^14
123: * return +0 if x <= -2^14
124: */
125: if (fe->fe_f2.fp_exp >= 14) {
126: if (fe->fe_f2.fp_sign) {
127: fe->fe_f2.fp_class = FPC_ZERO;
128: fe->fe_f2.fp_sign = 0;
129: } else {
130: fe->fe_f2.fp_class = FPC_INF;
131: }
132: return &fe->fe_f2;
133: }
134:
135: CPYFPN(&x, &fe->fe_f2);
136:
137: /* k = round(x / ln2) */
138: CPYFPN(&fe->fe_f1, &fe->fe_f2);
139: fpu_const(&fe->fe_f2, FPU_CONST_LN_2);
140: fp = fpu_div(fe);
141: CPYFPN(&fe->fe_f2, fp);
142: fp = fpu_int(fe);
143: if (ISZERO(fp)) {
144: /* k = 0 */
145: CPYFPN(&fe->fe_f2, &x);
146: fp = fpu_etox_taylor(fe);
147: return fp;
148: }
149: /* extract k as integer format from fpn format */
150: k = fp->fp_mant[0] >> (FP_LG - fp->fp_exp);
151: if (fp->fp_sign)
152: k *= -1;
153:
154: /* exp(r) = exp(x - k * ln2) */
155: CPYFPN(&fe->fe_f1, fp);
156: fpu_const(&fe->fe_f2, FPU_CONST_LN_2);
157: fp = fpu_mul(fe);
158: fp->fp_sign = !fp->fp_sign;
159: CPYFPN(&fe->fe_f1, fp);
160: CPYFPN(&fe->fe_f2, &x);
161: fp = fpu_add(fe);
162: CPYFPN(&fe->fe_f2, fp);
163: fp = fpu_etox_taylor(fe);
164:
165: /* 2^k */
166: fp->fp_exp += k;
167:
168: return fp;
169: }
170:
171: /*
172: * exp(x) - 1
173: */
174: struct fpn *
175: fpu_etoxm1(struct fpemu *fe)
176: {
177: struct fpn *fp;
178:
179: /* if x is +0/-0, return +0/-0 */
180: if (ISZERO(&fe->fe_f2))
181: return &fe->fe_f2;
182:
183: fp = fpu_etox(fe);
184:
185: CPYFPN(&fe->fe_f1, fp);
186: /* build a 1.0 */
187: fp = fpu_const(&fe->fe_f2, FPU_CONST_1);
188: fe->fe_f2.fp_sign = !fe->fe_f2.fp_sign;
189: /* fp = f2 - 1.0 */
190: fp = fpu_add(fe);
191:
192: return fp;
193: }
194:
195: /*
196: * 10^x = exp(x * ln10)
197: */
198: struct fpn *
199: fpu_tentox(struct fpemu *fe)
200: {
201: struct fpn *fp;
202:
203: /* build a ln10 */
204: fp = fpu_const(&fe->fe_f1, FPU_CONST_LN_10);
205: /* fp = ln10 * f2 */
206: fp = fpu_mul(fe);
207:
208: /* copy the result to the src opr */
209: CPYFPN(&fe->fe_f2, fp);
210:
211: return fpu_etox(fe);
212: }
213:
214: /*
215: * 2^x = exp(x * ln2)
216: */
217: struct fpn *
218: fpu_twotox(struct fpemu *fe)
219: {
220: struct fpn *fp;
221:
222: /* build a ln2 */
223: fp = fpu_const(&fe->fe_f1, FPU_CONST_LN_2);
224: /* fp = ln2 * f2 */
225: fp = fpu_mul(fe);
226:
227: /* copy the result to the src opr */
228: CPYFPN(&fe->fe_f2, fp);
229:
230: return fpu_etox(fe);
231: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.