fpu_trig.c revision 1.9
1 /* $NetBSD: fpu_trig.c,v 1.9 2013/04/11 13:27:11 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.9 2013/04/11 13:27:11 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 r = fpu_cmp(fe); 281 if (r->fp_sign == 0) { 282 CPYFPN(&fe->fe_f1, &x); 283 CPYFPN(&fe->fe_f2, &p); 284 fe->fe_f2.fp_sign = 1; 285 r = fpu_add(fe); 286 CPYFPN(&x, r); 287 288 sign ^= 1; 289 } 290 291 /* p <- pi/2 */ 292 p.fp_exp--; 293 294 /* 295 * if (x > pi/2) 296 * cos(x) is -sin(x - pi/2) 297 * else 298 * cos(x) 299 */ 300 CPYFPN(&fe->fe_f1, &x); 301 CPYFPN(&fe->fe_f2, &p); 302 r = fpu_cmp(fe); 303 if (r->fp_sign == 0) { 304 CPYFPN(&fe->fe_f1, &x); 305 CPYFPN(&fe->fe_f2, &p); 306 fe->fe_f2.fp_sign = 1; 307 r = fpu_add(fe); 308 309 CPYFPN(&fe->fe_f2, r); 310 r = fpu_sin_halfpi(fe); 311 sign ^= 1; 312 } else { 313 CPYFPN(&fe->fe_f2, &x); 314 r = fpu_cos_halfpi(fe); 315 } 316 317 CPYFPN(&fe->fe_f2, r); 318 fe->fe_f2.fp_sign = sign; 319 320 return &fe->fe_f2; 321 } 322 323 /* 324 * sin(x): 325 * 326 * if (x < 0) { 327 * x = abs(x); 328 * sign = 1; 329 * } 330 * if (x > 2*pi) { 331 * x %= 2*pi; 332 * } 333 * if (x > pi) { 334 * x -= pi; 335 * sign inverse; 336 * } 337 * if (x > pi/2) { 338 * y = cos(x - pi/2); 339 * } else { 340 * y = sin(x); 341 * } 342 * if (sign) { 343 * y = -y; 344 * } 345 */ 346 struct fpn * 347 fpu_sin(struct fpemu *fe) 348 { 349 struct fpn x; 350 struct fpn p; 351 struct fpn *r; 352 int sign; 353 354 if (ISNAN(&fe->fe_f2)) 355 return &fe->fe_f2; 356 if (ISINF(&fe->fe_f2)) 357 return fpu_newnan(fe); 358 359 CPYFPN(&x, &fe->fe_f2); 360 361 /* x = abs(input) */ 362 sign = x.fp_sign; 363 x.fp_sign = 0; 364 365 /* p <- 2*pi */ 366 fpu_const(&p, FPU_CONST_PI); 367 p.fp_exp++; 368 369 /* 370 * if (x > 2*pi*N) 371 * sin(x) is sin(x - 2*pi*N) 372 */ 373 CPYFPN(&fe->fe_f1, &x); 374 CPYFPN(&fe->fe_f2, &p); 375 r = fpu_cmp(fe); 376 if (r->fp_sign == 0) { 377 CPYFPN(&fe->fe_f1, &x); 378 CPYFPN(&fe->fe_f2, &p); 379 r = fpu_mod(fe); 380 CPYFPN(&x, r); 381 } 382 383 /* p <- pi */ 384 p.fp_exp--; 385 386 /* 387 * if (x > pi) 388 * sin(x) is -sin(x - pi) 389 */ 390 CPYFPN(&fe->fe_f1, &x); 391 CPYFPN(&fe->fe_f2, &p); 392 r = fpu_cmp(fe); 393 if (r->fp_sign == 0) { 394 CPYFPN(&fe->fe_f1, &x); 395 CPYFPN(&fe->fe_f2, &p); 396 fe->fe_f2.fp_sign = 1; 397 r = fpu_add(fe); 398 CPYFPN(&x, r); 399 400 sign ^= 1; 401 } 402 403 /* p <- pi/2 */ 404 p.fp_exp--; 405 406 /* 407 * if (x > pi/2) 408 * sin(x) is cos(x - pi/2) 409 * else 410 * sin(x) 411 */ 412 CPYFPN(&fe->fe_f1, &x); 413 CPYFPN(&fe->fe_f2, &p); 414 r = fpu_cmp(fe); 415 if (r->fp_sign == 0) { 416 CPYFPN(&fe->fe_f1, &x); 417 CPYFPN(&fe->fe_f2, &p); 418 fe->fe_f2.fp_sign = 1; 419 r = fpu_add(fe); 420 421 CPYFPN(&fe->fe_f2, r); 422 r = fpu_cos_halfpi(fe); 423 } else { 424 CPYFPN(&fe->fe_f2, &x); 425 r = fpu_sin_halfpi(fe); 426 } 427 428 CPYFPN(&fe->fe_f2, r); 429 fe->fe_f2.fp_sign = sign; 430 431 return &fe->fe_f2; 432 } 433 434 /* 435 * tan(x) = sin(x) / cos(x) 436 */ 437 struct fpn * 438 fpu_tan(struct fpemu *fe) 439 { 440 struct fpn x; 441 struct fpn s; 442 struct fpn *r; 443 444 if (ISNAN(&fe->fe_f2)) 445 return &fe->fe_f2; 446 if (ISINF(&fe->fe_f2)) 447 return fpu_newnan(fe); 448 449 CPYFPN(&x, &fe->fe_f2); 450 451 /* sin(x) */ 452 CPYFPN(&fe->fe_f2, &x); 453 r = fpu_sin(fe); 454 CPYFPN(&s, r); 455 456 /* cos(x) */ 457 CPYFPN(&fe->fe_f2, &x); 458 r = fpu_cos(fe); 459 CPYFPN(&fe->fe_f2, r); 460 461 CPYFPN(&fe->fe_f1, &s); 462 r = fpu_div(fe); 463 464 CPYFPN(&fe->fe_f2, r); 465 466 return &fe->fe_f2; 467 } 468 469 struct fpn * 470 fpu_sincos(struct fpemu *fe, int regc) 471 { 472 struct fpn x; 473 struct fpn *r; 474 475 CPYFPN(&x, &fe->fe_f2); 476 477 /* cos(x) */ 478 r = fpu_cos(fe); 479 fpu_implode(fe, r, FTYPE_EXT, &fe->fe_fpframe->fpf_regs[regc]); 480 481 /* sin(x) */ 482 CPYFPN(&fe->fe_f2, &x); 483 r = fpu_sin(fe); 484 return r; 485 } 486
Indexes created Mon Sep 29 21:09:56 GMT 2025