/* ibm2ieee - Converts a number from IBM 370 single precision floating
 point format to IEEE 754 single precision format. For normalized
numbers, the IBM format has greater range but less precision than the
IEEE format. Numbers within the overlapping range are converted
exactly. Numbers which are too large are converted to IEEE Infinity
with the correct sign. Numbers which are too small are converted to
IEEE denormalized numbers with a potential loss of precision (including
complete loss of precision which results in zero with the correct
sign). When precision is lost, rounding is toward zero (because it's
fast and easy -- if someone really wants round to nearest it shouldn't
be TOO difficult). */

/*
 * https://bytes.com/topic/c/answers/221981-c-code-converting-ibm-370-floating-point-ieee-754-a
 */

// [[Rcpp::export]]
void ibm2ieee(void *to, const void *from, int len)
{
  unsigned fr; /* fraction */
  int exp; /* exponent */
  int sgn; /* sign */

  for (; len-- > 0; to = (char *)to + 4, from = (char *)from + 4) {
    /* split into sign, exponent, and fraction */
    /* eliminate conversion -- expect little endian */
    fr = *(long *)from; /* pick up value */
    sgn = fr >> 31; /* save sign */
    fr <<= 1; /* shift sign out */
    exp = fr >> 25; /* save exponent */
    fr <<= 7; /* shift exponent out */

    if (fr == 0) { /* short-circuit for zero */
    exp = 0;
      goto done;
    }

    /* adjust exponent from base 16 offset 64 radix point before first digit
    to base 2 offset 127 radix point after first digit */
    /* (exp - 64) * 4 + 127 - 1 == exp * 4 - 256 + 126 == (exp << 2) - 130 */
    exp = (exp << 2) - 130;

    /* (re)normalize */
    while (fr < 0x80000000) { /* 3 times max for normalized input */
    --exp;
      fr <<= 1;
    }

    if (exp <= 0) { /* underflow */
    if (exp < -24) /* complete underflow - return properly signed zero */
    fr = 0;
    else /* partial underflow - return denormalized number */
    fr >>= -exp;
    exp = 0;
    } else if (exp >= 255) { /* overflow - return infinity */
    fr = 0;
      exp = 255;
    } else { /* just a plain old number - remove the assumed high bit */
    fr <<= 1;
    }

    done:
      /* put the pieces back together and return it */
      *(unsigned *)to = (fr >> 9) | (exp << 23) | (sgn << 31);

  }
}