Slow Inverse Square Root

sisr.R

## there are simpler ways to do this in R. See
## for example https://github.com/itchyny/fastinvsqrt/blob/main/src/r/fastinvsqrt.r
## contains a reasonable way to do this unreasonable thing.
## What follows is an unreasonable way to do an unreasonable thing.

doubleTofloatBits <- function(x) {
  # Convert double to character in scientific notation
  x_char <- format(x, scientific = TRUE)
  
  # Parse the character representation
  parts <- strsplit(x_char, "e")[[1]]
  mantissa <- as.numeric(parts[1])
  exponent <- as.integer(parts[2])
  
  # Normalize the mantissa to be within [1, 2) as per IEEE 754
  while (abs(mantissa) >= 2) {
    mantissa <- mantissa / 2
    exponent <- exponent + 1
  }
  while (abs(mantissa) < 1 && mantissa != 0) {
    mantissa <- mantissa * 2
    exponent <- exponent - 1
  }
  
  # Calculate sign, exponent, and fraction bits
  sign_bit <- as.integer(x < 0)
  exp_bits <- bitwShiftL(exponent + 127, 23)
  frac_bits <- bitwAnd(as.integer((abs(mantissa) - 1) * 2^23), 0x7FFFFF)
  
  # Combine all bits
  float_bits <- bitwOr(bitwOr(bitwShiftL(sign_bit, 31), exp_bits), frac_bits)
  
  # Return raw vector
  packBits(intToBits(float_bits), type="raw")
}

floatBitsToDouble <- function(raw_bits) {
  # Convert raw vector to unsigned 32-bit integer
  bits <- sum(as.integer(raw_bits) * 256^(0:3))
  
  # Extract sign, exponent, and fraction
  sign <- ifelse(bits >= 2^31, -1, 1)
  exponent <- bitwAnd(bitwShiftR(bits, 23), 0xFF)
  fraction <- bitwAnd(bits, 0x7FFFFF)
  
  if (exponent == 0) {
    # Denormalized number (including zero)
    unbiased_exp <- -126
    significand <- fraction / 2^23
  } else if (exponent == 255) {
    # Infinity or NaN
    return(if (fraction == 0) sign * Inf else NaN)
  } else {
    # Normalized number
    unbiased_exp <- exponent - 127
    significand <- 1 + fraction / 2^23
  }
  
  # Return reconstructed value: (-1)^sign * 2^unbiased_exp * significand
  sign * 2^unbiased_exp * significand
}

slow_inv_sqrt <- function(input) {
  # Hard-code 0x5f3759df
  magic_int <- 1597463007L
  # Convert input to integer for bit manipulation
  input_int <- packBits(rawToBits(doubleTofloatBits(input)), type="integer")
  # Perform subtraction and right shift by 1 bit
  result_int <- magic_int - bitwShiftR(input_int, 1)
  # Convert result back to raw
  result_raw <- packBits(intToBits(result_int), type="raw")
  
  # Convert raw result back to double
  floatBitsToDouble(result_raw)
}