e_log10.c 2.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101
  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_log10(x)
  12. * Return the base 10 logarithm of x
  13. *
  14. * Method :
  15. * Let log10_2hi = leading 40 bits of log10(2) and
  16. * log10_2lo = log10(2) - log10_2hi,
  17. * ivln10 = 1/log(10) rounded.
  18. * Then
  19. * n = ilogb(x),
  20. * if(n<0) n = n+1;
  21. * x = scalbn(x,-n);
  22. * log10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x))
  23. *
  24. * Note 1:
  25. * To guarantee log10(10**n)=n, where 10**n is normal, the rounding
  26. * mode must set to Round-to-Nearest.
  27. * Note 2:
  28. * [1/log(10)] rounded to 53 bits has error .198 ulps;
  29. * log10 is monotonic at all binary break points.
  30. *
  31. * Special cases:
  32. * log10(x) is NaN with signal if x < 0;
  33. * log10(+INF) is +INF with no signal; log10(0) is -INF with signal;
  34. * log10(NaN) is that NaN with no signal;
  35. * log10(10**N) = N for N=0,1,...,22.
  36. *
  37. * Constants:
  38. * The hexadecimal values are the intended ones for the following constants.
  39. * The decimal values may be used, provided that the compiler will convert
  40. * from decimal to binary accurately enough to produce the hexadecimal values
  41. * shown.
  42. */
  43. #include "math.h"
  44. #include "math_private.h"
  45. static const double
  46. two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
  47. ivln10 = 4.34294481903251816668e-01, /* 0x3FDBCB7B, 0x1526E50E */
  48. log10_2hi = 3.01029995663611771306e-01, /* 0x3FD34413, 0x509F6000 */
  49. log10_2lo = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
  50. static const double zero = 0.0;
  51. double attribute_hidden __ieee754_log10(double x)
  52. {
  53. double y,z;
  54. int32_t i,k,hx;
  55. u_int32_t lx;
  56. EXTRACT_WORDS(hx,lx,x);
  57. k=0;
  58. if (hx < 0x00100000) { /* x < 2**-1022 */
  59. if (((hx&0x7fffffff)|lx)==0)
  60. return -two54/zero; /* log(+-0)=-inf */
  61. if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
  62. k -= 54; x *= two54; /* subnormal number, scale up x */
  63. GET_HIGH_WORD(hx,x);
  64. }
  65. if (hx >= 0x7ff00000) return x+x;
  66. k += (hx>>20)-1023;
  67. i = ((u_int32_t)k&0x80000000)>>31;
  68. hx = (hx&0x000fffff)|((0x3ff-i)<<20);
  69. y = (double)(k+i);
  70. SET_HIGH_WORD(x,hx);
  71. z = y*log10_2lo + ivln10*__ieee754_log(x);
  72. return z+y*log10_2hi;
  73. }
  74. /*
  75. * wrapper log10(X)
  76. */
  77. #ifndef _IEEE_LIBM
  78. double log10(double x)
  79. {
  80. double z = __ieee754_log10(x);
  81. if (_LIB_VERSION == _IEEE_ || isnan(x))
  82. return z;
  83. if (x <= 0.0) {
  84. if(x == 0.0)
  85. return __kernel_standard(x, x, 18); /* log10(0) */
  86. return __kernel_standard(x, x, 19); /* log10(x<0) */
  87. }
  88. return z;
  89. }
  90. #else
  91. strong_alias(__ieee754_log10, log10)
  92. #endif
  93. libm_hidden_def(log10)