fpu_trig.c revision 1.10
1 /* $NetBSD: fpu_trig.c,v 1.10 2013/04/18 13:40:25 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_trig.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_trig.c,v 1.10 2013/04/18 13:40:25 isaki Exp $"); 61 62 #include "fpu_emulate.h" 63 64 static inline struct fpn *fpu_cos_halfpi(struct fpemu *); 65 static inline struct fpn *fpu_sin_halfpi(struct fpemu *); 66 67 struct fpn * 68 fpu_acos(struct fpemu *fe) 69 { 70 /* stub */ 71 return &fe->fe_f2; 72 } 73 74 struct fpn * 75 fpu_asin(struct fpemu *fe) 76 { 77 /* stub */ 78 return &fe->fe_f2; 79 } 80 81 struct fpn * 82 fpu_atan(struct fpemu *fe) 83 { 84 /* stub */ 85 return &fe->fe_f2; 86 } 87 88 /* 89 * taylor expansion used by sin(), cos(), sinh(), cosh(). 90 * hyperb is for sinh(), cosh(). 91 */ 92 struct fpn * 93 fpu_sincos_taylor(struct fpemu *fe, struct fpn *s0, uint32_t f, int hyperb) 94 { 95 struct fpn res; 96 struct fpn x2; 97 struct fpn *s1; 98 struct fpn *r; 99 int sign; 100 uint32_t k; 101 102 /* x2 := x * x */ 103 CPYFPN(&fe->fe_f1, &fe->fe_f2); 104 r = fpu_mul(fe); 105 CPYFPN(&x2, r); 106 107 /* res := s0 */ 108 CPYFPN(&res, s0); 109 110 sign = 1; /* sign := (-1)^n */ 111 112 for (;;) { 113 /* (f1 :=) s0 * x^2 */ 114 CPYFPN(&fe->fe_f1, s0); 115 CPYFPN(&fe->fe_f2, &x2); 116 r = fpu_mul(fe); 117 CPYFPN(&fe->fe_f1, r); 118 119 /* 120 * for sin(), s1 := s0 * x^2 / (2n+1)2n 121 * for cos(), s1 := s0 * x^2 / 2n(2n-1) 122 */ 123 k = f * (f + 1); 124 fpu_explode(fe, &fe->fe_f2, FTYPE_LNG, &k); 125 s1 = fpu_div(fe); 126 127 /* break if s1 is enough small */ 128 if (ISZERO(s1)) 129 break; 130 if (res.fp_exp - s1->fp_exp >= FP_NMANT) 131 break; 132 133 /* s0 := s1 for next loop */ 134 CPYFPN(s0, s1); 135 136 CPYFPN(&fe->fe_f2, s1); 137 if (!hyperb) { 138 /* (-1)^n */ 139 fe->fe_f2.fp_sign = sign; 140 } 141 142 /* res += s1 */ 143 CPYFPN(&fe->fe_f1, &res); 144 r = fpu_add(fe); 145 CPYFPN(&res, r); 146 147 f += 2; 148 sign ^= 1; 149 } 150 151 CPYFPN(&fe->fe_f2, &res); 152 return &fe->fe_f2; 153 } 154 155 /* 156 * inf (-1)^n 2n 157 * cos(x) = \sum { ------ * x } 158 * n=0 2n! 159 * 160 * x^2 x^4 x^6 161 * = 1 - --- + --- - --- + ... 162 * 2! 4! 6! 163 * 164 * = s0 - s1 + s2 - s3 + ... 165 * 166 * x*x x*x x*x 167 * s0 = 1, s1 = ---*s0, s2 = ---*s1, s3 = ---*s2, ... 168 * 1*2 3*4 5*6 169 * 170 * here 0 <= x < pi/2 171 */ 172 static inline struct fpn * 173 fpu_cos_halfpi(struct fpemu *fe) 174 { 175 struct fpn s0; 176 177 /* s0 := 1 */ 178 fpu_const(&s0, FPU_CONST_1); 179 180 return fpu_sincos_taylor(fe, &s0, 1, 0); 181 } 182 183 /* 184 * inf (-1)^n 2n+1 185 * sin(x) = \sum { ------- * x } 186 * n=0 (2n+1)! 187 * 188 * x^3 x^5 x^7 189 * = x - --- + --- - --- + ... 190 * 3! 5! 7! 191 * 192 * = s0 - s1 + s2 - s3 + ... 193 * 194 * x*x x*x x*x 195 * s0 = x, s1 = ---*s0, s2 = ---*s1, s3 = ---*s2, ... 196 * 2*3 4*5 6*7 197 * 198 * here 0 <= x < pi/2 199 */ 200 static inline struct fpn * 201 fpu_sin_halfpi(struct fpemu *fe) 202 { 203 struct fpn s0; 204 205 /* s0 := x */ 206 CPYFPN(&s0, &fe->fe_f2); 207 208 return fpu_sincos_taylor(fe, &s0, 2, 0); 209 } 210 211 /* 212 * cos(x): 213 * 214 * if (x < 0) { 215 * x = abs(x); 216 * } 217 * if (x > 2*pi) { 218 * x %= 2*pi; 219 * } 220 * if (x > pi) { 221 * x -= pi; 222 * sign inverse; 223 * } 224 * if (x > pi/2) { 225 * y = sin(x - pi/2); 226 * sign inverse; 227 * } else { 228 * y = cos(x); 229 * } 230 * if (sign) { 231 * y = -y; 232 * } 233 */ 234 struct fpn * 235 fpu_cos(struct fpemu *fe) 236 { 237 struct fpn x; 238 struct fpn p; 239 struct fpn *r; 240 int sign; 241 242 if (ISNAN(&fe->fe_f2)) 243 return &fe->fe_f2; 244 if (ISINF(&fe->fe_f2)) 245 return fpu_newnan(fe); 246 247 CPYFPN(&x, &fe->fe_f2); 248 249 /* x = abs(input) */ 250 x.fp_sign = 0; 251 sign = 0; 252 253 /* p <- 2*pi */ 254 fpu_const(&p, FPU_CONST_PI); 255 p.fp_exp++; 256 257 /* 258 * if (x > 2*pi*N) 259 * cos(x) is cos(x - 2*pi*N) 260 */ 261 CPYFPN(&fe->fe_f1, &x); 262 CPYFPN(&fe->fe_f2, &p); 263 r = fpu_cmp(fe); 264 if (r->fp_sign == 0) { 265 CPYFPN(&fe->fe_f1, &x); 266 CPYFPN(&fe->fe_f2, &p); 267 r = fpu_mod(fe); 268 CPYFPN(&x, r); 269 } 270 271 /* p <- pi */ 272 p.fp_exp--; 273 274 /* 275 * if (x > pi) 276 * cos(x) is -cos(x - pi) 277 */ 278 CPYFPN(&fe->fe_f1, &x); 279 CPYFPN(&fe->fe_f2, &p); 280 fe->fe_f2.fp_sign = 1; 281 r = fpu_add(fe); 282 if (r->fp_sign == 0) { 283 CPYFPN(&x, r); 284 sign ^= 1; 285 } 286 287 /* p <- pi/2 */ 288 p.fp_exp--; 289 290 /* 291 * if (x > pi/2) 292 * cos(x) is -sin(x - pi/2) 293 * else 294 * cos(x) 295 */ 296 CPYFPN(&fe->fe_f1, &x); 297 CPYFPN(&fe->fe_f2, &p); 298 fe->fe_f2.fp_sign = 1; 299 r = fpu_add(fe); 300 if (r->fp_sign == 0) { 301 CPYFPN(&fe->fe_f2, r); 302 r = fpu_sin_halfpi(fe); 303 sign ^= 1; 304 } else { 305 CPYFPN(&fe->fe_f2, &x); 306 r = fpu_cos_halfpi(fe); 307 } 308 309 CPYFPN(&fe->fe_f2, r); 310 fe->fe_f2.fp_sign = sign; 311 312 return &fe->fe_f2; 313 } 314 315 /* 316 * sin(x): 317 * 318 * if (x < 0) { 319 * x = abs(x); 320 * sign = 1; 321 * } 322 * if (x > 2*pi) { 323 * x %= 2*pi; 324 * } 325 * if (x > pi) { 326 * x -= pi; 327 * sign inverse; 328 * } 329 * if (x > pi/2) { 330 * y = cos(x - pi/2); 331 * } else { 332 * y = sin(x); 333 * } 334 * if (sign) { 335 * y = -y; 336 * } 337 */ 338 struct fpn * 339 fpu_sin(struct fpemu *fe) 340 { 341 struct fpn x; 342 struct fpn p; 343 struct fpn *r; 344 int sign; 345 346 if (ISNAN(&fe->fe_f2)) 347 return &fe->fe_f2; 348 if (ISINF(&fe->fe_f2)) 349 return fpu_newnan(fe); 350 351 CPYFPN(&x, &fe->fe_f2); 352 353 /* x = abs(input) */ 354 sign = x.fp_sign; 355 x.fp_sign = 0; 356 357 /* p <- 2*pi */ 358 fpu_const(&p, FPU_CONST_PI); 359 p.fp_exp++; 360 361 /* 362 * if (x > 2*pi*N) 363 * sin(x) is sin(x - 2*pi*N) 364 */ 365 CPYFPN(&fe->fe_f1, &x); 366 CPYFPN(&fe->fe_f2, &p); 367 r = fpu_cmp(fe); 368 if (r->fp_sign == 0) { 369 CPYFPN(&fe->fe_f1, &x); 370 CPYFPN(&fe->fe_f2, &p); 371 r = fpu_mod(fe); 372 CPYFPN(&x, r); 373 } 374 375 /* p <- pi */ 376 p.fp_exp--; 377 378 /* 379 * if (x > pi) 380 * sin(x) is -sin(x - pi) 381 */ 382 CPYFPN(&fe->fe_f1, &x); 383 CPYFPN(&fe->fe_f2, &p); 384 fe->fe_f2.fp_sign = 1; 385 r = fpu_add(fe); 386 if (r->fp_sign == 0) { 387 CPYFPN(&x, r); 388 sign ^= 1; 389 } 390 391 /* p <- pi/2 */ 392 p.fp_exp--; 393 394 /* 395 * if (x > pi/2) 396 * sin(x) is cos(x - pi/2) 397 * else 398 * sin(x) 399 */ 400 CPYFPN(&fe->fe_f1, &x); 401 CPYFPN(&fe->fe_f2, &p); 402 fe->fe_f2.fp_sign = 1; 403 r = fpu_add(fe); 404 if (r->fp_sign == 0) { 405 CPYFPN(&fe->fe_f2, r); 406 r = fpu_cos_halfpi(fe); 407 } else { 408 CPYFPN(&fe->fe_f2, &x); 409 r = fpu_sin_halfpi(fe); 410 } 411 412 CPYFPN(&fe->fe_f2, r); 413 fe->fe_f2.fp_sign = sign; 414 415 return &fe->fe_f2; 416 } 417 418 /* 419 * tan(x) = sin(x) / cos(x) 420 */ 421 struct fpn * 422 fpu_tan(struct fpemu *fe) 423 { 424 struct fpn x; 425 struct fpn s; 426 struct fpn *r; 427 428 if (ISNAN(&fe->fe_f2)) 429 return &fe->fe_f2; 430 if (ISINF(&fe->fe_f2)) 431 return fpu_newnan(fe); 432 433 CPYFPN(&x, &fe->fe_f2); 434 435 /* sin(x) */ 436 CPYFPN(&fe->fe_f2, &x); 437 r = fpu_sin(fe); 438 CPYFPN(&s, r); 439 440 /* cos(x) */ 441 CPYFPN(&fe->fe_f2, &x); 442 r = fpu_cos(fe); 443 CPYFPN(&fe->fe_f2, r); 444 445 CPYFPN(&fe->fe_f1, &s); 446 r = fpu_div(fe); 447 448 CPYFPN(&fe->fe_f2, r); 449 450 return &fe->fe_f2; 451 } 452 453 struct fpn * 454 fpu_sincos(struct fpemu *fe, int regc) 455 { 456 struct fpn x; 457 struct fpn *r; 458 459 CPYFPN(&x, &fe->fe_f2); 460 461 /* cos(x) */ 462 r = fpu_cos(fe); 463 fpu_implode(fe, r, FTYPE_EXT, &fe->fe_fpframe->fpf_regs[regc]); 464 465 /* sin(x) */ 466 CPYFPN(&fe->fe_f2, &x); 467 r = fpu_sin(fe); 468 return r; 469 } 470
Indexes created Tue Sep 30 06:09:59 GMT 2025