1
0

sinl.c 2.3 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182
  1. /* origin: FreeBSD /usr/src/lib/msun/src/s_sinl.c */
  2. /*-
  3. * Copyright (c) 2007 Steven G. Kargl
  4. * All rights reserved.
  5. *
  6. * Redistribution and use in source and binary forms, with or without
  7. * modification, are permitted provided that the following conditions
  8. * are met:
  9. * 1. Redistributions of source code must retain the above copyright
  10. * notice unmodified, this list of conditions, and the following
  11. * 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. *
  16. * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
  17. * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
  18. * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
  19. * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
  20. * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
  21. * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  22. * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  23. * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  24. * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
  25. * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  26. */
  27. #include "libm.h"
  28. #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
  29. long double sinl(long double x)
  30. {
  31. return sin(x);
  32. }
  33. #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
  34. long double sinl(long double x)
  35. {
  36. union IEEEl2bits z;
  37. int e0, s;
  38. long double y[2];
  39. long double hi, lo;
  40. z.e = x;
  41. s = z.bits.sign;
  42. z.bits.sign = 0;
  43. /* If x = +-0 or x is a subnormal number, then sin(x) = x */
  44. if (z.bits.exp == 0)
  45. return x;
  46. /* If x = NaN or Inf, then sin(x) = NaN. */
  47. if (z.bits.exp == 32767)
  48. return (x - x) / (x - x);
  49. /* Optimize the case where x is already within range. */
  50. if (z.e < M_PI_4) {
  51. hi = __sinl(z.e, 0, 0);
  52. return s ? -hi : hi;
  53. }
  54. e0 = __rem_pio2l(x, y);
  55. hi = y[0];
  56. lo = y[1];
  57. switch (e0 & 3) {
  58. case 0:
  59. hi = __sinl(hi, lo, 1);
  60. break;
  61. case 1:
  62. hi = __cosl(hi, lo);
  63. break;
  64. case 2:
  65. hi = - __sinl(hi, lo, 1);
  66. break;
  67. case 3:
  68. hi = - __cosl(hi, lo);
  69. break;
  70. }
  71. return hi;
  72. }
  73. #endif