![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to fpu_trig.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Isaki's NoNo m68k/m88k emulator] / nono / fpe|
Return to fpu_trig.c CVS log | Up to [Isaki's NoNo m68k/m88k emulator] / nono / fpe |
1.1 ! root 1: /* $NetBSD: fpu_trig.c,v 1.17 2016/08/06 00:58:55 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 "fpu_emulate.h" ! 60: ! 61: /* ! 62: * arccos(x) = pi/2 - arcsin(x) ! 63: */ ! 64: struct fpn * ! 65: fpu_acos(struct fpemu *fe) ! 66: { ! 67: struct fpn *r; ! 68: ! 69: if (ISNAN(&fe->fe_f2)) ! 70: return &fe->fe_f2; ! 71: if (ISINF(&fe->fe_f2)) ! 72: return fpu_newnan(fe); ! 73: ! 74: r = fpu_asin(fe); ! 75: if (ISNAN(r)) { ! 76: return r; ! 77: } ! 78: CPYFPN(&fe->fe_f2, r); ! 79: ! 80: /* pi/2 - asin(x) */ ! 81: fpu_const(&fe->fe_f1, FPU_CONST_PI); ! 82: fe->fe_f1.fp_exp--; ! 83: fe->fe_f2.fp_sign = !fe->fe_f2.fp_sign; ! 84: r = fpu_add(fe); ! 85: ! 86: /* ! 87: * small hack ! 88: */ ! 89: r->fp_sign = 0; ! 90: ! 91: return r; ! 92: } ! 93: ! 94: /* ! 95: * x ! 96: * arcsin(x) = arctan(---------------) ! 97: * sqrt(1 - x^2) ! 98: */ ! 99: struct fpn * ! 100: fpu_asin(struct fpemu *fe) ! 101: { ! 102: struct fpn x; ! 103: struct fpn *r; ! 104: ! 105: if (ISNAN(&fe->fe_f2)) ! 106: return &fe->fe_f2; ! 107: if (ISZERO(&fe->fe_f2)) ! 108: return &fe->fe_f2; ! 109: ! 110: if (ISINF(&fe->fe_f2)) ! 111: return fpu_newnan(fe); ! 112: ! 113: /* ! 114: * +pi/2 if x == 1 ! 115: * -pi/2 if x == -1 ! 116: */ ! 117: r = &fe->fe_f2; ! 118: if (r->fp_exp == 0 && r->fp_mant[0] == FP_1 && ! 119: r->fp_mant[1] == 0 && r->fp_mant[2] == 0) { ! 120: ! 121: /* f1 <- pi/2 */ ! 122: fpu_const(&fe->fe_f1, FPU_CONST_PI); ! 123: fe->fe_f1.fp_exp--; ! 124: ! 125: fe->fe_f1.fp_sign = r->fp_sign; ! 126: return &fe->fe_f1; ! 127: } ! 128: ! 129: /* NAN if |x| > 1 */ ! 130: if (fe->fe_f2.fp_exp >= 0) ! 131: return fpu_newnan(fe); ! 132: ! 133: CPYFPN(&x, &fe->fe_f2); ! 134: ! 135: /* x^2 */ ! 136: CPYFPN(&fe->fe_f1, &fe->fe_f2); ! 137: r = fpu_mul(fe); ! 138: ! 139: /* 1 - x^2 */ ! 140: CPYFPN(&fe->fe_f2, r); ! 141: fe->fe_f2.fp_sign = 1; ! 142: fpu_const(&fe->fe_f1, FPU_CONST_1); ! 143: r = fpu_add(fe); ! 144: ! 145: /* sqrt(1-x^2) */ ! 146: CPYFPN(&fe->fe_f2, r); ! 147: r = fpu_sqrt(fe); ! 148: ! 149: /* x/sqrt */ ! 150: CPYFPN(&fe->fe_f2, r); ! 151: CPYFPN(&fe->fe_f1, &x); ! 152: r = fpu_div(fe); ! 153: ! 154: /* arctan */ ! 155: CPYFPN(&fe->fe_f2, r); ! 156: return fpu_atan(fe); ! 157: } ! 158: ! 159: /* ! 160: * arctan(x): ! 161: * ! 162: * if (x < 0) { ! 163: * x = abs(x); ! 164: * sign = 1; ! 165: * } ! 166: * y = arctan(x); ! 167: * if (sign) { ! 168: * y = -y; ! 169: * } ! 170: */ ! 171: struct fpn * ! 172: fpu_atan(struct fpemu *fe) ! 173: { ! 174: struct fpn a; ! 175: struct fpn x; ! 176: struct fpn v; ! 177: ! 178: if (ISNAN(&fe->fe_f2)) ! 179: return &fe->fe_f2; ! 180: if (ISZERO(&fe->fe_f2)) ! 181: return &fe->fe_f2; ! 182: ! 183: CPYFPN(&a, &fe->fe_f2); ! 184: ! 185: if (ISINF(&fe->fe_f2)) { ! 186: /* f2 <- pi/2 */ ! 187: fpu_const(&fe->fe_f2, FPU_CONST_PI); ! 188: fe->fe_f2.fp_exp--; ! 189: ! 190: fe->fe_f2.fp_sign = a.fp_sign; ! 191: return &fe->fe_f2; ! 192: } ! 193: ! 194: fpu_const(&x, FPU_CONST_1); ! 195: fpu_const(&fe->fe_f2, FPU_CONST_0); ! 196: CPYFPN(&v, &fe->fe_f2); ! 197: fpu_cordit1(fe, &x, &a, &fe->fe_f2, &v); ! 198: ! 199: return &fe->fe_f2; ! 200: } ! 201: ! 202: ! 203: /* ! 204: * fe_f1 := sin(in) ! 205: * fe_f2 := cos(in) ! 206: */ ! 207: static void ! 208: __fpu_sincos_cordic(struct fpemu *fe, const struct fpn *in) ! 209: { ! 210: struct fpn a; ! 211: struct fpn v; ! 212: ! 213: CPYFPN(&a, in); ! 214: fpu_const(&fe->fe_f1, FPU_CONST_0); ! 215: CPYFPN(&fe->fe_f2, &fpu_cordic_inv_gain1); ! 216: fpu_const(&v, FPU_CONST_1); ! 217: v.fp_sign = 1; ! 218: fpu_cordit1(fe, &fe->fe_f2, &fe->fe_f1, &a, &v); ! 219: } ! 220: ! 221: /* ! 222: * cos(x): ! 223: * ! 224: * if (x < 0) { ! 225: * x = abs(x); ! 226: * } ! 227: * if (x >= 2*pi) { ! 228: * x %= 2*pi; ! 229: * } ! 230: * if (x >= pi) { ! 231: * x -= pi; ! 232: * sign inverse; ! 233: * } ! 234: * if (x >= pi/2) { ! 235: * y = sin(x - pi/2); ! 236: * sign inverse; ! 237: * } else { ! 238: * y = cos(x); ! 239: * } ! 240: * if (sign) { ! 241: * y = -y; ! 242: * } ! 243: */ ! 244: struct fpn * ! 245: fpu_cos(struct fpemu *fe) ! 246: { ! 247: struct fpn x; ! 248: struct fpn p; ! 249: struct fpn *r; ! 250: int sign; ! 251: ! 252: if (ISNAN(&fe->fe_f2)) ! 253: return &fe->fe_f2; ! 254: if (ISINF(&fe->fe_f2)) ! 255: return fpu_newnan(fe); ! 256: ! 257: /* if x is +0/-0, return 1 */ ! 258: if (ISZERO(&fe->fe_f2)) { ! 259: fpu_const(&fe->fe_f2, FPU_CONST_1); ! 260: return &fe->fe_f2; ! 261: } ! 262: ! 263: /* x = abs(input) */ ! 264: sign = 0; ! 265: CPYFPN(&x, &fe->fe_f2); ! 266: x.fp_sign = 0; ! 267: ! 268: /* p <- 2*pi */ ! 269: fpu_const(&p, FPU_CONST_PI); ! 270: p.fp_exp++; ! 271: ! 272: /* ! 273: * if (x >= 2*pi) ! 274: * cos(x) is cos(x % 2*pi) ! 275: */ ! 276: CPYFPN(&fe->fe_f1, &x); ! 277: CPYFPN(&fe->fe_f2, &p); ! 278: r = fpu_cmp(fe); ! 279: if (r->fp_sign == 0) { ! 280: CPYFPN(&fe->fe_f1, &x); ! 281: CPYFPN(&fe->fe_f2, &p); ! 282: r = fpu_modrem(fe, 1); ! 283: CPYFPN(&x, r); ! 284: } ! 285: ! 286: /* p <- pi */ ! 287: p.fp_exp--; ! 288: ! 289: /* ! 290: * if (x >= pi) ! 291: * cos(x) is -cos(x - pi) ! 292: */ ! 293: CPYFPN(&fe->fe_f1, &x); ! 294: CPYFPN(&fe->fe_f2, &p); ! 295: fe->fe_f2.fp_sign = 1; ! 296: r = fpu_add(fe); ! 297: if (r->fp_sign == 0) { ! 298: CPYFPN(&x, r); ! 299: sign ^= 1; ! 300: } ! 301: ! 302: /* p <- pi/2 */ ! 303: p.fp_exp--; ! 304: ! 305: /* ! 306: * if (x >= pi/2) ! 307: * cos(x) is -sin(x - pi/2) ! 308: * else ! 309: * cos(x) ! 310: */ ! 311: CPYFPN(&fe->fe_f1, &x); ! 312: CPYFPN(&fe->fe_f2, &p); ! 313: fe->fe_f2.fp_sign = 1; ! 314: r = fpu_add(fe); ! 315: if (r->fp_sign == 0) { ! 316: __fpu_sincos_cordic(fe, r); ! 317: r = &fe->fe_f1; ! 318: sign ^= 1; ! 319: } else { ! 320: __fpu_sincos_cordic(fe, &x); ! 321: r = &fe->fe_f2; ! 322: } ! 323: r->fp_sign = sign; ! 324: return r; ! 325: } ! 326: ! 327: /* ! 328: * sin(x): ! 329: * ! 330: * if (x < 0) { ! 331: * x = abs(x); ! 332: * sign = 1; ! 333: * } ! 334: * if (x >= 2*pi) { ! 335: * x %= 2*pi; ! 336: * } ! 337: * if (x >= pi) { ! 338: * x -= pi; ! 339: * sign inverse; ! 340: * } ! 341: * if (x >= pi/2) { ! 342: * y = cos(x - pi/2); ! 343: * } else { ! 344: * y = sin(x); ! 345: * } ! 346: * if (sign) { ! 347: * y = -y; ! 348: * } ! 349: */ ! 350: struct fpn * ! 351: fpu_sin(struct fpemu *fe) ! 352: { ! 353: struct fpn x; ! 354: struct fpn p; ! 355: struct fpn *r; ! 356: int sign; ! 357: ! 358: if (ISNAN(&fe->fe_f2)) ! 359: return &fe->fe_f2; ! 360: if (ISINF(&fe->fe_f2)) ! 361: return fpu_newnan(fe); ! 362: ! 363: /* if x is +0/-0, return +0/-0 */ ! 364: if (ISZERO(&fe->fe_f2)) ! 365: return &fe->fe_f2; ! 366: ! 367: /* x = abs(input) */ ! 368: sign = fe->fe_f2.fp_sign; ! 369: CPYFPN(&x, &fe->fe_f2); ! 370: x.fp_sign = 0; ! 371: ! 372: /* p <- 2*pi */ ! 373: fpu_const(&p, FPU_CONST_PI); ! 374: p.fp_exp++; ! 375: ! 376: /* ! 377: * if (x >= 2*pi) ! 378: * sin(x) is sin(x % 2*pi) ! 379: */ ! 380: CPYFPN(&fe->fe_f1, &x); ! 381: CPYFPN(&fe->fe_f2, &p); ! 382: r = fpu_cmp(fe); ! 383: if (r->fp_sign == 0) { ! 384: CPYFPN(&fe->fe_f1, &x); ! 385: CPYFPN(&fe->fe_f2, &p); ! 386: r = fpu_modrem(fe, 1); ! 387: CPYFPN(&x, r); ! 388: } ! 389: ! 390: /* p <- pi */ ! 391: p.fp_exp--; ! 392: ! 393: /* ! 394: * if (x >= pi) ! 395: * sin(x) is -sin(x - pi) ! 396: */ ! 397: CPYFPN(&fe->fe_f1, &x); ! 398: CPYFPN(&fe->fe_f2, &p); ! 399: fe->fe_f2.fp_sign = 1; ! 400: r = fpu_add(fe); ! 401: if (r->fp_sign == 0) { ! 402: CPYFPN(&x, r); ! 403: sign ^= 1; ! 404: } ! 405: ! 406: /* if x is small enough, sin(x) is x */ ! 407: /* XXX the border value is not well considered... */ ! 408: if (r->fp_exp < -25) { ! 409: r->fp_sign = sign; ! 410: r->fp_sticky = 1; ! 411: return r; ! 412: } ! 413: ! 414: /* p <- pi/2 */ ! 415: p.fp_exp--; ! 416: ! 417: /* ! 418: * if (x >= pi/2) ! 419: * sin(x) is cos(x - pi/2) ! 420: * else ! 421: * sin(x) ! 422: */ ! 423: CPYFPN(&fe->fe_f1, &x); ! 424: CPYFPN(&fe->fe_f2, &p); ! 425: fe->fe_f2.fp_sign = 1; ! 426: r = fpu_add(fe); ! 427: if (r->fp_sign == 0) { ! 428: __fpu_sincos_cordic(fe, r); ! 429: r = &fe->fe_f2; ! 430: } else { ! 431: __fpu_sincos_cordic(fe, &x); ! 432: r = &fe->fe_f1; ! 433: } ! 434: r->fp_sign = sign; ! 435: return r; ! 436: } ! 437: ! 438: /* ! 439: * tan(x): ! 440: * ! 441: * if (x < 0) { ! 442: * x = abs(x); ! 443: * sign = 1; ! 444: * } ! 445: * if (x >= pi) { ! 446: * x %= pi; ! 447: * } ! 448: * if (x >= pi/2) { ! 449: * y = tan(pi - x); ! 450: * sign inverse; ! 451: * } else { ! 452: * y = tan(x); ! 453: * } ! 454: */ ! 455: struct fpn * ! 456: fpu_tan(struct fpemu *fe) ! 457: { ! 458: struct fpn x; ! 459: struct fpn p; ! 460: struct fpn *r; ! 461: int sign; ! 462: ! 463: if (ISNAN(&fe->fe_f2)) ! 464: return &fe->fe_f2; ! 465: if (ISINF(&fe->fe_f2)) ! 466: return fpu_newnan(fe); ! 467: ! 468: /* if x is +0/-0, return +0/-0 */ ! 469: if (ISZERO(&fe->fe_f2)) ! 470: return &fe->fe_f2; ! 471: ! 472: /* x = abs(input) */ ! 473: sign = fe->fe_f2.fp_sign; ! 474: CPYFPN(&x, &fe->fe_f2); ! 475: x.fp_sign = 0; ! 476: ! 477: /* p <- pi */ ! 478: fpu_const(&p, FPU_CONST_PI); ! 479: ! 480: /* ! 481: * if (x >= pi) ! 482: * tan(x) is tan(x % pi) ! 483: */ ! 484: CPYFPN(&fe->fe_f1, &x); ! 485: CPYFPN(&fe->fe_f2, &p); ! 486: r = fpu_cmp(fe); ! 487: if (r->fp_sign == 0) { ! 488: CPYFPN(&fe->fe_f1, &x); ! 489: CPYFPN(&fe->fe_f2, &p); ! 490: r = fpu_modrem(fe, 1); ! 491: CPYFPN(&x, r); ! 492: } ! 493: ! 494: /* ! 495: * if (x >= pi/2) ! 496: * tan(x) is -tan(pi - x) ! 497: * else ! 498: * tan(x) ! 499: */ ! 500: CPYFPN(&fe->fe_f1, &x); ! 501: CPYFPN(&fe->fe_f2, &p); ! 502: fe->fe_f2.fp_exp--; ! 503: r = fpu_cmp(fe); ! 504: if (r->fp_sign == 0) { ! 505: CPYFPN(&fe->fe_f1, &p); ! 506: CPYFPN(&fe->fe_f2, &x); ! 507: fe->fe_f2.fp_sign = 1; ! 508: r = fpu_add(fe); ! 509: sign ^= 1; ! 510: } else { ! 511: r = &x; ! 512: } ! 513: ! 514: /* tan(x) = sin(x)/cos(x) */ ! 515: __fpu_sincos_cordic(fe, r); ! 516: r = fpu_div(fe); ! 517: ! 518: if (sign) { ! 519: r->fp_sign = !r->fp_sign; ! 520: } ! 521: return r; ! 522: } ! 523: ! 524: struct fpn * ! 525: fpu_sincos(struct fpemu *fe, int regc) ! 526: { ! 527: struct fpn x; ! 528: struct fpn p; ! 529: struct fpn *r; ! 530: ! 531: if (ISNAN(&fe->fe_f2)) { ! 532: fpu_implode(fe, &fe->fe_f2, FTYPE_EXT, ! 533: &fe->fe_fpframe->fpf_regs[regc * 3]); ! 534: return &fe->fe_f2; ! 535: } ! 536: if (ISINF(&fe->fe_f2)) { ! 537: r = fpu_newnan(fe); ! 538: fpu_implode(fe, r, FTYPE_EXT, ! 539: &fe->fe_fpframe->fpf_regs[regc * 3]); ! 540: return r; ! 541: } ! 542: ! 543: /* if x is +0/-0, cos(x) is 1 and sin(x) is +0/-0 */ ! 544: if (ISZERO(&fe->fe_f2)) { ! 545: fpu_const(&fe->fe_f1, FPU_CONST_1); ! 546: fpu_implode(fe, &fe->fe_f1, FTYPE_EXT, ! 547: &fe->fe_fpframe->fpf_regs[regc * 3]); ! 548: return &fe->fe_f2; ! 549: } ! 550: ! 551: /* p <- 2*pi */ ! 552: fpu_const(&p, FPU_CONST_PI); ! 553: p.fp_exp++; ! 554: ! 555: /* ! 556: * if (x >= 2*pi) ! 557: * x %= 2*pi ! 558: */ ! 559: CPYFPN(&x, &fe->fe_f2); ! 560: CPYFPN(&fe->fe_f1, &x); ! 561: CPYFPN(&fe->fe_f2, &p); ! 562: r = fpu_cmp(fe); ! 563: if (r->fp_sign == 0) { ! 564: CPYFPN(&fe->fe_f1, &x); ! 565: CPYFPN(&fe->fe_f2, &p); ! 566: r = fpu_modrem(fe, 1); ! 567: CPYFPN(&x, r); ! 568: } ! 569: ! 570: __fpu_sincos_cordic(fe, &x); ! 571: ! 572: /* cos(x) */ ! 573: fpu_round_prec(fe, &fe->fe_f2); ! 574: fpu_implode(fe, &fe->fe_f2, FTYPE_EXT, &fe->fe_fpframe->fpf_regs[regc * 3]); ! 575: ! 576: /* sin(x) */ ! 577: return &fe->fe_f1; ! 578: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.