fusl/src/math/scalbn.c

#include <math.h>
#include <stdint.h>

double scalbn(double x, int n) {
  union {
    double f;
    uint64_t i;
  } u;
  double_t y = x;

  if (n > 1023) {
    y *= 0x1p1023;
    n -= 1023;
    if (n > 1023) {
      y *= 0x1p1023;
      n -= 1023;
      if (n > 1023)
        n = 1023;
    }
  } else if (n < -1022) {
    y *= 0x1p-1022;
    n += 1022;
    if (n < -1022) {
      y *= 0x1p-1022;
      n += 1022;
      if (n < -1022)
        n = -1022;
    }
  }
  u.i = (uint64_t)(0x3ff + n) << 52;
  x = y * u.f;
  return x;
}