writer_float_math.hpp

#pragma once
 
template <typename Float>
struct Writer::DecimalScale
{
  int exponent = 0;  
  Float scale = 1;   
};
 
template <typename Float>
struct Writer::ScientificDigits
{
  Float digits =
      0;  
  Float scale = 1;   
  int exponent = 0;  
};
 
template <typename Float>
Float Writer::Power10(int exponent)
{
  Float result = 1;
  Float base = 10;
  unsigned int remaining = static_cast<unsigned int>(exponent < 0 ? -exponent : exponent);
 
  while (remaining != 0)
  {
    if ((remaining & 1U) != 0U)
    {
      if (exponent < 0)
      {
        result /= base;
      }
      else
      {
        result *= base;
      }
    }
 
    remaining >>= 1U;
    if (remaining != 0U)
    {
      base *= base;
    }
  }
 
  return result;
}
 
template <typename Float>
Float Writer::RoundDecimal(Float value, uint8_t precision)
{
  Float scale = Power10<Float>(static_cast<int>(precision));
  Float scaled = value * scale;
  if (!std::isfinite(scaled))
  {
    // Scaling overflowed: propagate as infinity so callers' isfinite() guard
    // rejects the result cleanly instead of silently returning the original
    // value and later producing OUT_OF_RANGE. Input is expected to be a
    // non-negative magnitude, so plain infinity is sufficient.
    return std::numeric_limits<Float>::infinity();
  }
 
  return std::nearbyint(scaled) / scale;
}
 
template <typename Float>
Writer::ScientificDigits<Float> Writer::RoundScientificDigits(Float value,
                                                              uint8_t precision)
{
  ScientificDigits<Float> result{};
  auto normalized = NormalizeDecimal(value);
  result.exponent = (value == 0) ? 0 : normalized.exponent;
  Float decimal_scale = Power10<Float>(result.exponent);
  Float mantissa = (value == 0) ? 0 : value / decimal_scale;
  result.scale = Power10<Float>(static_cast<int>(precision));
  result.digits = std::nearbyint(mantissa * result.scale);
  if (result.digits >= static_cast<Float>(10) * result.scale)
  {
    result.digits /= 10;
    ++result.exponent;
  }
 
  return result;
}
 
template <typename Float>
Writer::DecimalScale<Float> Writer::NormalizeDecimal(Float value)
{
  DecimalScale<Float> normalized{};
  if (value == 0)
  {
    return normalized;
  }
 
  int binary_exponent = 0;
  std::frexp(value, &binary_exponent);
  constexpr Float log10_of_2 = static_cast<Float>(0.30102999566398119521373889472449L);
  normalized.exponent =
      static_cast<int>(static_cast<Float>(binary_exponent - 1) * log10_of_2);
  normalized.scale = Power10<Float>(normalized.exponent);
 
  Float scaled = value / normalized.scale;
  while (scaled < 1)
  {
    normalized.scale /= 10;
    --normalized.exponent;
    scaled *= 10;
  }
  while (scaled >= 10)
  {
    normalized.scale *= 10;
    ++normalized.exponent;
    scaled /= 10;
  }
 
  return normalized;
}
 
template <typename Float>
uint8_t Writer::ExtractDigit(Float& value, Float scale)
{
  Float scaled = value / scale;
  // Bias to correct floating-point rounding when the true digit value is an
  // integer but division leaves it just below (e.g. 1.9999999... instead of 2).
  // 1e-12 is effective for double (epsilon ~2.2e-16) but effectively zero for
  // float (epsilon ~1.2e-7). For float, digit extraction may have up to 1-ULP
  // error in the last digit; this is a known limitation of the approach.
  // Do NOT increase this bias to fix float: values like 1.999999f have a
  // legitimate float representation ~9.5e-7 below 2.0, so any bias large
  // enough to "fix" float rounding would also cause false carry on such values.
  auto digit = static_cast<int>(scaled + static_cast<Float>(1e-12L));
  if (digit < 0)
  {
    digit = 0;
  }
  else if (digit > 9)
  {
    digit = 9;
  }
 
  value -= static_cast<Float>(digit) * scale;
  // Zero-clamp: clear tiny negative residuals that are FP rounding artifacts.
  // 1e-9 is effective for double but coarse for long double; for float it is
  // large enough to clear genuine residuals without clamping real fractional
  // remainders (float residuals after digit extraction are < epsilon * scale).
  Float epsilon = scale * static_cast<Float>(1e-9L);
  if (value < 0 && value > -epsilon)
  {
    value = 0;
  }
 
  return static_cast<uint8_t>(digit);
}