e_remainder.c 1.9 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485
  1. /*
  2. * ====================================================
  3. * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  4. *
  5. * Developed at SunPro, a Sun Microsystems, Inc. business.
  6. * Permission to use, copy, modify, and distribute this
  7. * software is freely granted, provided that this notice
  8. * is preserved.
  9. * ====================================================
  10. */
  11. /* __ieee754_remainder(x,p)
  12. * Return :
  13. * returns x REM p = x - [x/p]*p as if in infinite
  14. * precise arithmetic, where [x/p] is the (infinite bit)
  15. * integer nearest x/p (in half way case choose the even one).
  16. * Method :
  17. * Based on fmod() return x-[x/p]chopped*p exactlp.
  18. */
  19. #include "math.h"
  20. #include "math_private.h"
  21. static const double zero = 0.0;
  22. double __ieee754_remainder(double x, double p)
  23. {
  24. int32_t hx,hp;
  25. u_int32_t sx,lx,lp;
  26. double p_half;
  27. EXTRACT_WORDS(hx,lx,x);
  28. EXTRACT_WORDS(hp,lp,p);
  29. sx = hx&0x80000000;
  30. hp &= 0x7fffffff;
  31. hx &= 0x7fffffff;
  32. /* purge off exception values */
  33. if((hp|lp)==0) return (x*p)/(x*p); /* p = 0 */
  34. if((hx>=0x7ff00000)|| /* x not finite */
  35. ((hp>=0x7ff00000)&& /* p is NaN */
  36. (((hp-0x7ff00000)|lp)!=0)))
  37. return (x*p)/(x*p);
  38. if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p); /* now x < 2p */
  39. if (((hx-hp)|(lx-lp))==0) return zero*x;
  40. x = fabs(x);
  41. p = fabs(p);
  42. if (hp<0x00200000) {
  43. if(x+x>p) {
  44. x-=p;
  45. if(x+x>=p) x -= p;
  46. }
  47. } else {
  48. p_half = 0.5*p;
  49. if(x>p_half) {
  50. x-=p;
  51. if(x>=p_half) x -= p;
  52. }
  53. }
  54. GET_HIGH_WORD(hx,x);
  55. SET_HIGH_WORD(x,hx^sx);
  56. return x;
  57. }
  58. /*
  59. * wrapper remainder(x,p)
  60. */
  61. #ifndef _IEEE_LIBM
  62. double remainder(double x, double y)
  63. {
  64. double z = __ieee754_remainder(x, y);
  65. if (_LIB_VERSION == _IEEE_ || isnan(y))
  66. return z;
  67. if (y == 0.0)
  68. return __kernel_standard(x, y, 28); /* remainder(x,0) */
  69. return z;
  70. }
  71. strong_alias(remainder, drem)
  72. #else
  73. strong_alias(__ieee754_remainder, remainder)
  74. strong_alias(__ieee754_remainder, drem)
  75. #endif
  76. libm_hidden_def(remainder)