fpu_exp.c revision 1.8.16.1
1 /* $NetBSD: fpu_exp.c,v 1.8.16.1 2017/01/07 08:56:20 pgoyette 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.8.16.1 2017/01/07 08:56:20 pgoyette 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 j, 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 CPYFPN(&x, &fe->fe_f2); 123 124 /* k = round(x / ln2) */ 125 CPYFPN(&fe->fe_f1, &fe->fe_f2); 126 fpu_const(&fe->fe_f2, FPU_CONST_LN_2); 127 fp = fpu_div(fe); 128 CPYFPN(&fe->fe_f2, fp); 129 fp = fpu_int(fe); 130 if (ISZERO(fp)) { 131 /* k = 0 */ 132 CPYFPN(&fe->fe_f2, &x); 133 fp = fpu_etox_taylor(fe); 134 return fp; 135 } 136 /* extract k as integer format from fpn format */ 137 j = FP_LG - fp->fp_exp; 138 if (j < 0) { 139 if (fp->fp_sign) { 140 fp->fp_class = FPC_ZERO; /* k < -2^18 */ 141 fp->fp_sign = 0; 142 } else { 143 fp->fp_class = FPC_INF; /* k > 2^18 */ 144 } 145 return fp; 146 } 147 k = fp->fp_mant[0] >> j; 148 if (fp->fp_sign) 149 k *= -1; 150 151 /* exp(r) = exp(x - k * ln2) */ 152 CPYFPN(&fe->fe_f1, fp); 153 fpu_const(&fe->fe_f2, FPU_CONST_LN_2); 154 fp = fpu_mul(fe); 155 fp->fp_sign = !fp->fp_sign; 156 CPYFPN(&fe->fe_f1, fp); 157 CPYFPN(&fe->fe_f2, &x); 158 fp = fpu_add(fe); 159 CPYFPN(&fe->fe_f2, fp); 160 fp = fpu_etox_taylor(fe); 161 162 /* 2^k */ 163 fp->fp_exp += k; 164 165 return fp; 166 } 167 168 /* 169 * exp(x) - 1 170 */ 171 struct fpn * 172 fpu_etoxm1(struct fpemu *fe) 173 { 174 struct fpn *fp; 175 176 fp = fpu_etox(fe); 177 178 CPYFPN(&fe->fe_f1, fp); 179 /* build a 1.0 */ 180 fp = fpu_const(&fe->fe_f2, FPU_CONST_1); 181 fe->fe_f2.fp_sign = !fe->fe_f2.fp_sign; 182 /* fp = f2 - 1.0 */ 183 fp = fpu_add(fe); 184 185 return fp; 186 } 187 188 /* 189 * 10^x = exp(x * ln10) 190 */ 191 struct fpn * 192 fpu_tentox(struct fpemu *fe) 193 { 194 struct fpn *fp; 195 196 /* build a ln10 */ 197 fp = fpu_const(&fe->fe_f1, FPU_CONST_LN_10); 198 /* fp = ln10 * f2 */ 199 fp = fpu_mul(fe); 200 201 /* copy the result to the src opr */ 202 CPYFPN(&fe->fe_f2, fp); 203 204 return fpu_etox(fe); 205 } 206 207 /* 208 * 2^x = exp(x * ln2) 209 */ 210 struct fpn * 211 fpu_twotox(struct fpemu *fe) 212 { 213 struct fpn *fp; 214 215 /* build a ln2 */ 216 fp = fpu_const(&fe->fe_f1, FPU_CONST_LN_2); 217 /* fp = ln2 * f2 */ 218 fp = fpu_mul(fe); 219 220 /* copy the result to the src opr */ 221 CPYFPN(&fe->fe_f2, fp); 222 223 return fpu_etox(fe); 224 } 225
Indexes created Mon Oct 27 08:10:08 GMT 2025