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 * 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 * 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 * 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 * 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 * 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
Indexes created Thu Sep 25 16:09:42 GMT 2025