mdl_state_decimal.c

Go to the documentation of this file.

#include <stdint.h>
#include <string.h>
 
#include "mdl_state_internal.h"
#include "ra8_attributes.h"

typedef enum : uint8_t {
  k_state_big_words = 48, 
} mdl_state_big_limit_t;

typedef enum : int32_t {
  k_state_decimal_prime            = 5,     
  k_state_decimal_scale_max        = 400,   
  k_state_divide_bits              = 64,    
  k_state_binary64_fraction_bits   = 52,    
  k_state_binary64_precision_bits  = 53,    
  k_state_binary64_sign_shift      = 63,    
  k_state_binary64_exponent_bias   = 1023,  
  k_state_binary64_exponent_max    = 1023,  
  k_state_binary64_exponent_min    = -1022, 
  k_state_binary64_subnormal_scale = 1074,  
} mdl_state_decimal_limit_t;

typedef struct {
  uint32_t word[k_state_big_words]; 
  uint8_t  used;                    
} mdl_state_big_t;

RA8_INTERNAL static void internal_mdl_state_big_init(mdl_state_big_t* big, uint64_t value)
{
  memset(big, 0, sizeof(*big));
  big->word[0] = (uint32_t)value;
  big->word[1] = (uint32_t)(value >> 32U);
  big->used    = (uint8_t)((big->word[1] != 0U) ? 2U : 1U);
}

RA8_INTERNAL static void internal_mdl_state_big_trim(mdl_state_big_t* big)
{
  while ((big->used > 1U) && (big->word[big->used - 1U] == 0U)) {
    --big->used;
  }
}

RA8_INTERNAL static uint16_t internal_mdl_state_big_bits(const mdl_state_big_t* big)
{
  if ((big->used == 1U) && (big->word[0] == 0U)) {
    return 0U;
  }
  uint32_t high = big->word[big->used - 1U];
  uint8_t  bits = 0U;
  while (high != 0U) {
    high >>= 1U;
    ++bits;
  }
  return (uint16_t)(((uint16_t)(big->used - 1U) * 32U) + bits);
}

RA8_INTERNAL static bool internal_mdl_state_big_mul(mdl_state_big_t* big, uint32_t factor)
{
  uint64_t carry = 0U;
  for (uint8_t i = 0U; i < big->used; ++i) {
    const uint64_t product = ((uint64_t)big->word[i] * factor) + carry;
    big->word[i]           = (uint32_t)product;
    carry                  = product >> 32U;
  }
  if (carry != 0U) {
    if (big->used >= (uint8_t)k_state_big_words) {
      return false;
    }
    big->word[big->used] = (uint32_t)carry;
    ++big->used;
  }
  return true;
}

RA8_INTERNAL static bool
internal_mdl_state_big_shift(const mdl_state_big_t* src, uint16_t shift, mdl_state_big_t* dst)
{
  const uint16_t words = shift / 32U;
  const uint8_t  bits  = (uint8_t)(shift % 32U);
  if (((uint16_t)src->used + words + (bits != 0U ? 1U : 0U)) > k_state_big_words) {
    return false;
  }
  memset(dst, 0, sizeof(*dst));
  uint64_t carry = 0U;
  for (uint8_t i = 0U; i < src->used; ++i) {
    const uint64_t value           = ((uint64_t)src->word[i] << bits) | carry;
    dst->word[(uint16_t)i + words] = (uint32_t)value;
    carry                          = value >> 32U;
  }
  dst->used = (uint8_t)((uint16_t)src->used + words);
  if (carry != 0U) {
    dst->word[dst->used] = (uint32_t)carry;
    ++dst->used;
  }
  internal_mdl_state_big_trim(dst);
  return true;
}

RA8_INTERNAL static int internal_mdl_state_big_compare(const mdl_state_big_t* lhs,
                                                       const mdl_state_big_t* rhs)
{
  if (lhs->used != rhs->used) {
    return lhs->used > rhs->used ? 1 : -1;
  }
  for (uint8_t i = lhs->used; i > 0U; --i) {
    if (lhs->word[i - 1U] != rhs->word[i - 1U]) {
      return lhs->word[i - 1U] > rhs->word[i - 1U] ? 1 : -1;
    }
  }
  return 0;
}

RA8_INTERNAL static void internal_mdl_state_big_subtract(mdl_state_big_t*       lhs,
                                                         const mdl_state_big_t* rhs)
{
  uint64_t borrow = 0U;
  for (uint8_t i = 0U; i < lhs->used; ++i) {
    const uint64_t subtrahend = (uint64_t)(i < rhs->used ? rhs->word[i] : 0U) + borrow;
    const uint64_t current    = lhs->word[i];
    lhs->word[i]              = (uint32_t)(current - subtrahend);
    borrow                    = current < subtrahend ? 1U : 0U;
  }
  internal_mdl_state_big_trim(lhs);
}

RA8_INTERNAL static int internal_mdl_state_compare_power(const mdl_state_big_t* numerator,
                                                         const mdl_state_big_t* denominator,
                                                         int32_t                binary_scale,
                                                         int32_t                exponent)
{
  mdl_state_big_t shifted = {};
  const int32_t   shift   = binary_scale - exponent;
  if (shift >= 0) {
    if (!internal_mdl_state_big_shift(numerator, (uint16_t)shift, &shifted)) {
      return 1;
    }
    return internal_mdl_state_big_compare(&shifted, denominator);
  }
  if (!internal_mdl_state_big_shift(denominator, (uint16_t)-shift, &shifted)) {
    return -1;
  }
  return internal_mdl_state_big_compare(numerator, &shifted);
}

RA8_INTERNAL static bool internal_mdl_state_divide(const mdl_state_big_t* numerator,
                                                   const mdl_state_big_t* denominator,
                                                   int32_t                binary_shift,
                                                   uint64_t*              out)
{
  mdl_state_big_t work    = {};
  mdl_state_big_t divisor = {};
  if (binary_shift >= 0) {
    if (!internal_mdl_state_big_shift(numerator, (uint16_t)binary_shift, &work)) {
      return false;
    }
    divisor = *denominator;
  } else {
    work = *numerator;
    if (!internal_mdl_state_big_shift(denominator, (uint16_t)-binary_shift, &divisor)) {
      return false;
    }
  }
  uint64_t      quotient = 0U;
  const int32_t top =
    (int32_t)internal_mdl_state_big_bits(&work) - (int32_t)internal_mdl_state_big_bits(&divisor);
  if (top >= (int32_t)k_state_divide_bits) {
    return false;
  }
  for (int32_t bit = top; bit >= 0; --bit) {
    mdl_state_big_t trial = {};
    if (!internal_mdl_state_big_shift(&divisor, (uint16_t)bit, &trial)) {
      return false;
    }
    if (internal_mdl_state_big_compare(&work, &trial) >= 0) {
      internal_mdl_state_big_subtract(&work, &trial);
      quotient |= UINT64_C(1) << (uint32_t)bit;
    }
  }
  mdl_state_big_t twice = {};
  if (!internal_mdl_state_big_shift(&work, 1U, &twice)) {
    return false;
  }
  const int comparison = internal_mdl_state_big_compare(&twice, &divisor);
  if ((comparison > 0) || ((comparison == 0) && ((quotient & 1U) != 0U))) {
    ++quotient;
  }
  *out = quotient;
  return true;
}

RA8_INTERNAL static bool internal_mdl_state_build_rational(uint64_t         mantissa,
                                                           int32_t          decimal_scale,
                                                           mdl_state_big_t* numerator,
                                                           mdl_state_big_t* denominator)
{
  internal_mdl_state_big_init(numerator, mantissa);
  internal_mdl_state_big_init(denominator, 1U);
  const uint32_t   count  = (uint32_t)(decimal_scale < 0 ? -decimal_scale : decimal_scale);
  mdl_state_big_t* factor = decimal_scale < 0 ? denominator : numerator;
  for (uint32_t i = 0U; i < count; ++i) {
    if (!internal_mdl_state_big_mul(factor, (uint32_t)k_state_decimal_prime)) {
      return false;
    }
  }
  return true;
}

RA8_INTERNAL static bool internal_mdl_state_encode(uint64_t quotient,
                                                   int32_t  exponent,
                                                   bool     negative,
                                                   bool     normal,
                                                   double*  out)
{
  uint64_t bits = negative ? (UINT64_C(1) << (uint32_t)k_state_binary64_sign_shift) : 0U;
  if (normal) {
    if (quotient == (UINT64_C(1) << (uint32_t)k_state_binary64_precision_bits)) {
      quotient >>= 1U;
      ++exponent;
    }
    if ((exponent > (int32_t)k_state_binary64_exponent_max) ||
        (quotient < (UINT64_C(1) << (uint32_t)k_state_binary64_fraction_bits))) {
      return false;
    }
    bits |= (uint64_t)(exponent + (int32_t)k_state_binary64_exponent_bias)
            << (uint32_t)k_state_binary64_fraction_bits;
    bits |= quotient - (UINT64_C(1) << (uint32_t)k_state_binary64_fraction_bits);
  } else {
    if ((quotient == 0U) ||
        (quotient > (UINT64_C(1) << (uint32_t)k_state_binary64_fraction_bits))) {
      return false;
    }
    if (quotient == (UINT64_C(1) << (uint32_t)k_state_binary64_fraction_bits)) {
      bits |= UINT64_C(1) << (uint32_t)k_state_binary64_fraction_bits;
    } else {
      bits |= quotient;
    }
  }
  memcpy(out, &bits, sizeof(bits));
  return true;
}
 
RA8_PRIV bool priv_mdl_state_decimal_to_binary64(uint64_t mantissa,
                                                 int32_t  decimal_scale,
                                                 bool     negative,
                                                 double*  out)
{
  if ((out == nullptr) || (decimal_scale < -(int32_t)k_state_decimal_scale_max) ||
      (decimal_scale > (int32_t)k_state_decimal_scale_max)) {
    return false;
  }
  if (mantissa == 0U) {
    const uint64_t bits = negative ? (UINT64_C(1) << (uint32_t)k_state_binary64_sign_shift) : 0U;
    memcpy(out, &bits, sizeof(bits));
    return true;
  }
  mdl_state_big_t numerator   = {};
  mdl_state_big_t denominator = {};
  if (!internal_mdl_state_build_rational(mantissa, decimal_scale, &numerator, &denominator)) {
    return false;
  }
  int32_t exponent = (int32_t)internal_mdl_state_big_bits(&numerator) -
                     (int32_t)internal_mdl_state_big_bits(&denominator) + decimal_scale;
  if (internal_mdl_state_compare_power(&numerator, &denominator, decimal_scale, exponent) < 0) {
    --exponent;
  }
  if (exponent > (int32_t)k_state_binary64_exponent_max) {
    return false;
  }
  const bool    normal   = exponent >= (int32_t)k_state_binary64_exponent_min;
  const int32_t shift    = normal
                             ? (decimal_scale + (int32_t)k_state_binary64_fraction_bits - exponent)
                             : (decimal_scale + (int32_t)k_state_binary64_subnormal_scale);
  uint64_t      quotient = 0U;
  if (!internal_mdl_state_divide(&numerator, &denominator, shift, &quotient)) {
    return false;
  }
  return internal_mdl_state_encode(quotient, exponent, negative, normal, out);
}