fpu_hyperb.c revision 1.13
1 /* $NetBSD: fpu_hyperb.c,v 1.13 2013/04/20 04:55:44 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_hyperb.c 10/24/95 32 */ 33 34 /* 35 * Copyright (c) 2011 Tetsuya Isaki. All rights reserved. 36 * 37 * Redistribution and use in source and binary forms, with or without 38 * modification, are permitted provided that the following conditions 39 * are met: 40 * 1. Redistributions of source code must retain the above copyright 41 * notice, this list of conditions and the following disclaimer. 42 * 2. Redistributions in binary form must reproduce the above copyright 43 * notice, this list of conditions and the following disclaimer in the 44 * documentation and/or other materials provided with the distribution. 45 * 46 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR 47 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 48 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 49 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, 50 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, 51 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 52 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED 53 * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, 54 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 55 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 56 * SUCH DAMAGE. 57 */ 58 59 #include <sys/cdefs.h> 60 __KERNEL_RCSID(0, "$NetBSD: fpu_hyperb.c,v 1.13 2013/04/20 04:55:44 isaki Exp $"); 61 62 #include <machine/ieee.h> 63 64 #include "fpu_emulate.h" 65 66 /* The number of items to terminate the Taylor expansion */ 67 #define MAX_ITEMS (2000) 68 69 /* 70 * fpu_hyperb.c: defines the following functions 71 * 72 * fpu_atanh(), fpu_cosh(), fpu_sinh(), and fpu_tanh() 73 */ 74 75 /* 76 * 1 1 + x 77 * atanh(x) = ---*log(-------) 78 * 2 1 - x 79 */ 80 struct fpn * 81 fpu_atanh(struct fpemu *fe) 82 { 83 struct fpn x; 84 struct fpn t; 85 struct fpn *r; 86 87 if (ISNAN(&fe->fe_f2)) 88 return &fe->fe_f2; 89 if (ISINF(&fe->fe_f2)) 90 return fpu_newnan(fe); 91 92 /* 93 * if x is +0/-0, 68000PRM.pdf says it returns +0/-0 but 94 * my real 68882 returns +0 for both. 95 */ 96 if (ISZERO(&fe->fe_f2)) { 97 fe->fe_f2.fp_sign = 0; 98 return &fe->fe_f2; 99 } 100 101 /* 102 * -INF if x == -1 103 * +INF if x == 1 104 */ 105 r = &fe->fe_f2; 106 if (r->fp_exp == 0 && r->fp_mant[0] == FP_1 && 107 r->fp_mant[1] == 0 && r->fp_mant[2] == 0) { 108 r->fp_class = FPC_INF; 109 return r; 110 } 111 112 /* NAN if |x| > 1 */ 113 if (fe->fe_f2.fp_exp >= 0) 114 return fpu_newnan(fe); 115 116 CPYFPN(&x, &fe->fe_f2); 117 118 /* t := 1 - x */ 119 fpu_const(&fe->fe_f1, FPU_CONST_1); 120 fe->fe_f2.fp_sign = !fe->fe_f2.fp_sign; 121 r = fpu_add(fe); 122 CPYFPN(&t, r); 123 124 /* r := 1 + x */ 125 fpu_const(&fe->fe_f1, FPU_CONST_1); 126 CPYFPN(&fe->fe_f2, &x); 127 r = fpu_add(fe); 128 129 /* (1-x)/(1+x) */ 130 CPYFPN(&fe->fe_f1, r); 131 CPYFPN(&fe->fe_f2, &t); 132 r = fpu_div(fe); 133 134 /* log((1-x)/(1+x)) */ 135 CPYFPN(&fe->fe_f2, r); 136 r = fpu_logn(fe); 137 138 /* r /= 2 */ 139 r->fp_exp--; 140 141 return r; 142 } 143 144 /* 145 * taylor expansion used by sinh(), cosh(). 146 */ 147 static struct fpn * 148 __fpu_sinhcosh_taylor(struct fpemu *fe, struct fpn *s0, uint32_t f) 149 { 150 struct fpn res; 151 struct fpn x2; 152 struct fpn *s1; 153 struct fpn *r; 154 int sign; 155 uint32_t k; 156 157 /* x2 := x * x */ 158 CPYFPN(&fe->fe_f1, &fe->fe_f2); 159 r = fpu_mul(fe); 160 CPYFPN(&x2, r); 161 162 /* res := s0 */ 163 CPYFPN(&res, s0); 164 165 sign = 1; /* sign := (-1)^n */ 166 167 for (; f < (2 * MAX_ITEMS); ) { 168 /* (f1 :=) s0 * x^2 */ 169 CPYFPN(&fe->fe_f1, s0); 170 CPYFPN(&fe->fe_f2, &x2); 171 r = fpu_mul(fe); 172 CPYFPN(&fe->fe_f1, r); 173 174 /* 175 * for sinh(), s1 := s0 * x^2 / (2n+1)2n 176 * for cosh(), s1 := s0 * x^2 / 2n(2n-1) 177 */ 178 k = f * (f + 1); 179 fpu_explode(fe, &fe->fe_f2, FTYPE_LNG, &k); 180 s1 = fpu_div(fe); 181 182 /* break if s1 is enough small */ 183 if (ISZERO(s1)) 184 break; 185 if (res.fp_exp - s1->fp_exp >= EXT_FRACBITS) 186 break; 187 188 /* s0 := s1 for next loop */ 189 CPYFPN(s0, s1); 190 191 /* res += s1 */ 192 CPYFPN(&fe->fe_f2, s1); 193 CPYFPN(&fe->fe_f1, &res); 194 r = fpu_add(fe); 195 CPYFPN(&res, r); 196 197 f += 2; 198 sign ^= 1; 199 } 200 201 CPYFPN(&fe->fe_f2, &res); 202 return &fe->fe_f2; 203 } 204 205 struct fpn * 206 fpu_cosh(struct fpemu *fe) 207 { 208 struct fpn s0; 209 struct fpn *r; 210 211 if (ISNAN(&fe->fe_f2)) 212 return &fe->fe_f2; 213 214 if (ISINF(&fe->fe_f2)) { 215 fe->fe_f2.fp_sign = 0; 216 return &fe->fe_f2; 217 } 218 219 fpu_const(&s0, FPU_CONST_1); 220 r = __fpu_sinhcosh_taylor(fe, &s0, 1); 221 CPYFPN(&fe->fe_f2, r); 222 223 return &fe->fe_f2; 224 } 225 226 struct fpn * 227 fpu_sinh(struct fpemu *fe) 228 { 229 struct fpn s0; 230 struct fpn *r; 231 232 if (ISNAN(&fe->fe_f2)) 233 return &fe->fe_f2; 234 if (ISINF(&fe->fe_f2)) 235 return &fe->fe_f2; 236 237 CPYFPN(&s0, &fe->fe_f2); 238 r = __fpu_sinhcosh_taylor(fe, &s0, 2); 239 CPYFPN(&fe->fe_f2, r); 240 241 return &fe->fe_f2; 242 } 243 244 struct fpn * 245 fpu_tanh(struct fpemu *fe) 246 { 247 struct fpn x; 248 struct fpn s; 249 struct fpn *r; 250 int sign; 251 252 if (ISNAN(&fe->fe_f2)) 253 return &fe->fe_f2; 254 255 if (ISINF(&fe->fe_f2)) { 256 sign = fe->fe_f2.fp_sign; 257 fpu_const(&fe->fe_f2, FPU_CONST_1); 258 fe->fe_f2.fp_sign = sign; 259 return &fe->fe_f2; 260 } 261 262 CPYFPN(&x, &fe->fe_f2); 263 264 /* sinh(x) */ 265 CPYFPN(&fe->fe_f2, &x); 266 r = fpu_sinh(fe); 267 CPYFPN(&s, r); 268 269 /* cosh(x) */ 270 CPYFPN(&fe->fe_f2, &x); 271 r = fpu_cosh(fe); 272 CPYFPN(&fe->fe_f2, r); 273 274 CPYFPN(&fe->fe_f1, &s); 275 r = fpu_div(fe); 276 277 CPYFPN(&fe->fe_f2, r); 278 279 return &fe->fe_f2; 280 } 281
Indexes created Sun Sep 28 07:09:56 GMT 2025