fpu_hyperb.c revision 1.10
1 /* $NetBSD: fpu_hyperb.c,v 1.10 2013/04/19 14:05:12 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.10 2013/04/19 14:05:12 isaki Exp $"); 61 62 #include "fpu_emulate.h" 63 64 /* 65 * fpu_hyperb.c: defines the following functions 66 * 67 * fpu_atanh(), fpu_cosh(), fpu_sinh(), and fpu_tanh() 68 */ 69 70 /* 71 * 1 1 + x 72 * atanh(x) = ---*log(-------) 73 * 2 1 - x 74 */ 75 struct fpn * 76 fpu_atanh(struct fpemu *fe) 77 { 78 struct fpn x; 79 struct fpn t; 80 struct fpn *r; 81 82 if (ISNAN(&fe->fe_f2)) 83 return &fe->fe_f2; 84 if (ISINF(&fe->fe_f2)) 85 return fpu_newnan(fe); 86 87 /* 88 * if x is +0/-0, 68000PRM.pdf says it returns +0/-0 but 89 * my real 68882 returns +0 for both. 90 */ 91 if (ISZERO(&fe->fe_f2)) { 92 fe->fe_f2.fp_sign = 0; 93 return &fe->fe_f2; 94 } 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 (;;) { 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 sin(), s1 := s0 * x^2 / (2n+1)2n 171 * for cos(), 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 >= FP_NMANT) 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 CPYFPN(&fe->fe_f2, r); 217 218 return &fe->fe_f2; 219 } 220 221 struct fpn * 222 fpu_sinh(struct fpemu *fe) 223 { 224 struct fpn s0; 225 struct fpn *r; 226 227 if (ISNAN(&fe->fe_f2)) 228 return &fe->fe_f2; 229 if (ISINF(&fe->fe_f2)) 230 return &fe->fe_f2; 231 232 CPYFPN(&s0, &fe->fe_f2); 233 r = __fpu_sinhcosh_taylor(fe, &s0, 2); 234 CPYFPN(&fe->fe_f2, r); 235 236 return &fe->fe_f2; 237 } 238 239 struct fpn * 240 fpu_tanh(struct fpemu *fe) 241 { 242 struct fpn x; 243 struct fpn s; 244 struct fpn *r; 245 int sign; 246 247 if (ISNAN(&fe->fe_f2)) 248 return &fe->fe_f2; 249 250 if (ISINF(&fe->fe_f2)) { 251 sign = fe->fe_f2.fp_sign; 252 fpu_const(&fe->fe_f2, FPU_CONST_1); 253 fe->fe_f2.fp_sign = sign; 254 return &fe->fe_f2; 255 } 256 257 CPYFPN(&x, &fe->fe_f2); 258 259 /* sinh(x) */ 260 CPYFPN(&fe->fe_f2, &x); 261 r = fpu_sinh(fe); 262 CPYFPN(&s, r); 263 264 /* cosh(x) */ 265 CPYFPN(&fe->fe_f2, &x); 266 r = fpu_cosh(fe); 267 CPYFPN(&fe->fe_f2, r); 268 269 CPYFPN(&fe->fe_f1, &s); 270 r = fpu_div(fe); 271 272 CPYFPN(&fe->fe_f2, r); 273 274 return &fe->fe_f2; 275 } 276
Indexes created Tue Sep 30 06:09:59 GMT 2025