scalbn.c 576 B

123456789101112131415161718192021222324252627282930313233
  1. #include <math.h>
  2. #include <stdint.h>
  3. double scalbn(double x, int n)
  4. {
  5. union {double f; uint64_t i;} u;
  6. double_t y = x;
  7. if (n > 1023) {
  8. y *= 0x1p1023;
  9. n -= 1023;
  10. if (n > 1023) {
  11. y *= 0x1p1023;
  12. n -= 1023;
  13. if (n > 1023)
  14. n = 1023;
  15. }
  16. } else if (n < -1022) {
  17. /* make sure final n < -53 to avoid double
  18. rounding in the subnormal range */
  19. y *= 0x1p-1022 * 0x1p53;
  20. n += 1022 - 53;
  21. if (n < -1022) {
  22. y *= 0x1p-1022 * 0x1p53;
  23. n += 1022 - 53;
  24. if (n < -1022)
  25. n = -1022;
  26. }
  27. }
  28. u.i = (uint64_t)(0x3ff+n)<<52;
  29. x = y * u.f;
  30. return x;
  31. }