1 /* $NetBSD: fpu_hyperb.c,v 1.16 2013/10/11 03:37:08 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.16 2013/10/11 03:37:08 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 /* if x is +0/-0, return +0/-0 */ 93 if (ISZERO(&fe->fe_f2)) 94 return &fe->fe_f2; 95 96 /* 97 * -INF if x == -1 98 * +INF if x == 1 99 */ 100 r = &fe->fe_f2; 101 if (r->fp_exp == 0 && r->fp_mant[0] == FP_1 && 102 r->fp_mant[1] == 0 && r->fp_mant[2] == 0) { 103 r->fp_class = FPC_INF; 104 return r; 105 } 106 107 /* NAN if |x| > 1 */ 108 if (fe->fe_f2.fp_exp >= 0) 109 return fpu_newnan(fe); 110 111 CPYFPN(&x, &fe->fe_f2); 112 113 /* t := 1 - x */ 114 fpu_const(&fe->fe_f1, FPU_CONST_1); 115 fe->fe_f2.fp_sign = !fe->fe_f2.fp_sign; 116 r = fpu_add(fe); 117 CPYFPN(&t, r); 118 119 /* r := 1 + x */ 120 fpu_const(&fe->fe_f1, FPU_CONST_1); 121 CPYFPN(&fe->fe_f2, &x); 122 r = fpu_add(fe); 123 124 /* (1-x)/(1+x) */ 125 CPYFPN(&fe->fe_f1, r); 126 CPYFPN(&fe->fe_f2, &t); 127 r = fpu_div(fe); 128 129 /* log((1-x)/(1+x)) */ 130 CPYFPN(&fe->fe_f2, r); 131 r = fpu_logn(fe); 132 133 /* r /= 2 */ 134 r->fp_exp--; 135 136 return r; 137 } 138 139 /* 140 * taylor expansion used by sinh(), cosh(). 141 */ 142 static struct fpn * 143 __fpu_sinhcosh_taylor(struct fpemu *fe, struct fpn *s0, uint32_t f) 144 { 145 struct fpn res; 146 struct fpn x2; 147 struct fpn *s1; 148 struct fpn *r; 149 int sign; 150 uint32_t k; 151 152 /* x2 := x * x */ 153 CPYFPN(&fe->fe_f1, &fe->fe_f2); 154 r = fpu_mul(fe); 155 CPYFPN(&x2, r); 156 157 /* res := s0 */ 158 CPYFPN(&res, s0); 159 160 sign = 1; /* sign := (-1)^n */ 161 162 for (; f < (2 * MAX_ITEMS); ) { 163 /* (f1 :=) s0 * x^2 */ 164 CPYFPN(&fe->fe_f1, s0); 165 CPYFPN(&fe->fe_f2, &x2); 166 r = fpu_mul(fe); 167 CPYFPN(&fe->fe_f1, r); 168 169 /* 170 * for sinh(), s1 := s0 * x^2 / (2n+1)2n 171 * for cosh(), s1 := s0 * x^2 / 2n(2n-1) 172 */ 173 k = f * (f + 1); 174 fpu_explode(fe, &fe->fe_f2, FTYPE_LNG, &k); 175 s1 = fpu_div(fe); 176 177 /* break if s1 is enough small */ 178 if (ISZERO(s1)) 179 break; 180 if (res.fp_exp - s1->fp_exp >= EXT_FRACBITS) 181 break; 182 183 /* s0 := s1 for next loop */ 184 CPYFPN(s0, s1); 185 186 /* res += s1 */ 187 CPYFPN(&fe->fe_f2, s1); 188 CPYFPN(&fe->fe_f1, &res); 189 r = fpu_add(fe); 190 CPYFPN(&res, r); 191 192 f += 2; 193 sign ^= 1; 194 } 195 196 CPYFPN(&fe->fe_f2, &res); 197 return &fe->fe_f2; 198 } 199 200 struct fpn * 201 fpu_cosh(struct fpemu *fe) 202 { 203 struct fpn s0; 204 struct fpn *r; 205 206 if (ISNAN(&fe->fe_f2)) 207 return &fe->fe_f2; 208 209 if (ISINF(&fe->fe_f2)) { 210 fe->fe_f2.fp_sign = 0; 211 return &fe->fe_f2; 212 } 213 214 fpu_const(&s0, FPU_CONST_1); 215 r = __fpu_sinhcosh_taylor(fe, &s0, 1); 216 217 return r; 218 } 219 220 struct fpn * 221 fpu_sinh(struct fpemu *fe) 222 { 223 struct fpn s0; 224 struct fpn *r; 225 226 if (ISNAN(&fe->fe_f2)) 227 return &fe->fe_f2; 228 if (ISINF(&fe->fe_f2)) 229 return &fe->fe_f2; 230 231 /* if x is +0/-0, return +0/-0 */ 232 if (ISZERO(&fe->fe_f2)) 233 return &fe->fe_f2; 234 235 CPYFPN(&s0, &fe->fe_f2); 236 r = __fpu_sinhcosh_taylor(fe, &s0, 2); 237 238 return r; 239 } 240 241 struct fpn * 242 fpu_tanh(struct fpemu *fe) 243 { 244 struct fpn x; 245 struct fpn s; 246 struct fpn *r; 247 int sign; 248 249 if (ISNAN(&fe->fe_f2)) 250 return &fe->fe_f2; 251 252 /* if x is +0/-0, return +0/-0 */ 253 if (ISZERO(&fe->fe_f2)) 254 return &fe->fe_f2; 255 256 if (ISINF(&fe->fe_f2)) { 257 sign = fe->fe_f2.fp_sign; 258 fpu_const(&fe->fe_f2, FPU_CONST_1); 259 fe->fe_f2.fp_sign = sign; 260 return &fe->fe_f2; 261 } 262 263 CPYFPN(&x, &fe->fe_f2); 264 265 /* sinh(x) */ 266 CPYFPN(&fe->fe_f2, &x); 267 r = fpu_sinh(fe); 268 CPYFPN(&s, r); 269 270 /* cosh(x) */ 271 CPYFPN(&fe->fe_f2, &x); 272 r = fpu_cosh(fe); 273 CPYFPN(&fe->fe_f2, r); 274 275 CPYFPN(&fe->fe_f1, &s); 276 r = fpu_div(fe); 277 278 return r; 279 } 280