scalar_io.hpp

Go to the documentation of this file.

#pragma once
 
#include "numeric_limits.hpp"
#include "scalar_concepts.hpp"
#include "types.hpp"
#if UNI20_HAS_MPREAL
#include "mpreal.hpp"
#endif
#if UNI20_HAS_MPREAL
#include "mpcomplex.hpp"
#endif
 
#include <cerrno>
#include <cmath>
#include <cstdlib>
#include <fmt/format.h>
#include <istream>
#include <stdexcept>
#include <string>
#include <string_view>
#include <type_traits>
 
namespace uni20
{
 
enum class real_format_notation
{
  general,
  fixed,
  scientific
};
 
struct scalar_format_options
{
    int precision = -1;
    real_format_notation notation = real_format_notation::general;
    bool normalize_negative_zero = true;
    std::string_view imaginary_unit = "i";
    bool exact_fractions = true;
};
 
namespace detail
{
#if UNI20_HAS_MPREAL
// Decimal rounding of an exact rational must not pass through a binary float:
// even an arbitrarily wide temporary can move a decimal halfway case off its tie.
inline std::string format_exact_decimal(exact_rational const& value, scalar_format_options const& options)
{
  auto q = value.native_handle();
  int const requested = options.precision;
  if (requested < 0)
  {
    if (requested == -1 && mpz_cmp_ui(mpq_denref(q), 1) == 0) return value.to_string();
    throw std::invalid_argument("exact decimal output requires an explicit digit count or fraction formatting");
  }
  if (requested == 0 && options.notation == real_format_notation::general)
    throw std::invalid_argument("general decimal output requires a positive digit count");
 
  mp_integer numerator, denominator, power, rounded, remainder;
  mpz_abs(numerator.value, mpq_numref(q));
  mpz_set(denominator.value, mpq_denref(q));
  bool const negative = mpq_sgn(q) < 0;
  long exponent = 0;
  bool const fixed = options.notation == real_format_notation::fixed;
  bool const general = options.notation == real_format_notation::general;
  if (!fixed && mpz_sgn(numerator.value) != 0)
  {
    // Decimal digit estimates can be one too large. Correct them by exact comparisons.
    exponent = static_cast<long>(mpz_sizeinbase(numerator.value, 10)) -
               static_cast<long>(mpz_sizeinbase(denominator.value, 10));
    auto compare_power = [&](long e) {
      mpz_ui_pow_ui(power.value, 10, static_cast<unsigned long>(e < 0 ? -e : e));
      if (e >= 0)
      {
        mpz_mul(remainder.value, denominator.value, power.value);
        return mpz_cmp(numerator.value, remainder.value);
      }
      mpz_mul(remainder.value, numerator.value, power.value);
      return mpz_cmp(remainder.value, denominator.value);
    };
    while (compare_power(exponent) < 0)
      --exponent;
    while (compare_power(exponent + 1) >= 0)
      ++exponent;
  }
  long const fractional = general ? static_cast<long>(requested) - 1 : requested;
  long const scale = fixed ? requested : fractional - exponent;
  mpz_ui_pow_ui(power.value, 10, static_cast<unsigned long>(scale < 0 ? -scale : scale));
  if (scale >= 0)
    mpz_mul(numerator.value, numerator.value, power.value);
  else
    mpz_mul(denominator.value, denominator.value, power.value);
  mpz_tdiv_qr(rounded.value, remainder.value, numerator.value, denominator.value);
  mpz_mul_2exp(remainder.value, remainder.value, 1);
  int const halfway = mpz_cmp(remainder.value, denominator.value);
  if (halfway > 0 || (halfway == 0 && mpz_odd_p(rounded.value))) mpz_add_ui(rounded.value, rounded.value, 1);
  std::string digits(mpz_sizeinbase(rounded.value, 10) + 2, '\0');
  mpz_get_str(digits.data(), 10, rounded.value);
  digits.resize(std::char_traits<char>::length(digits.c_str()));
  if (!fixed && digits.size() > static_cast<std::size_t>(fractional + 1))
  {
    digits.pop_back(); // A rounding carry, e.g. 9.995 -> 1.000e+01.
    ++exponent;
  }
  bool const scientific = !fixed && (!general || exponent < -4 || exponent >= requested);
  long const decimal_places = fixed ? requested : (scientific ? fractional : fractional - exponent);
  if (decimal_places > 0)
  {
    auto places = static_cast<std::size_t>(decimal_places);
    if (digits.size() <= places) digits.insert(0, places + 1 - digits.size(), '0');
    digits.insert(digits.size() - places, 1, '.');
    if (general)
    {
      while (digits.back() == '0')
        digits.pop_back();
      if (digits.back() == '.') digits.pop_back();
    }
  }
  if (scientific) digits += fmt::format("e{}{:02}", exponent < 0 ? '-' : '+', exponent < 0 ? -exponent : exponent);
  if (negative) digits.insert(0, 1, '-');
  return digits;
}
 
inline std::string format_mpreal(mpreal const& value, scalar_format_options const& options)
{
  if (value.is_exact())
    return options.exact_fractions ? value.to_string() : format_exact_decimal(value.exact_value(), options);
  auto const requested = options.precision < 0 ? mpfr_get_str_ndigits(10, value.precision().bit_count())
                                               : static_cast<std::size_t>(options.precision);
  if (requested > static_cast<std::size_t>(std::numeric_limits<int>::max()))
    throw std::length_error("mpreal: decimal formatting precision is too large");
  char const* format = options.notation == real_format_notation::fixed        ? "%.*RNf"
                       : options.notation == real_format_notation::scientific ? "%.*RNe"
                                                                              : "%.*RNg";
  char* buffer = nullptr;
  int const count = mpfr_asprintf(&buffer, format, static_cast<int>(requested), value.native_handle());
  std::unique_ptr<char, decltype(&mpfr_free_str)> owned(buffer, &mpfr_free_str);
  if (count < 0) throw std::runtime_error("mpreal: decimal formatting failed");
  return std::string(buffer, static_cast<std::size_t>(count));
}
#endif
 
template <typename T> [[nodiscard]] int effective_precision(scalar_format_options const& options)
{
  if (options.precision >= 0)
  {
    return options.precision;
  }
  return uni20::numeric_limits<std::remove_cv_t<T>>::max_digits10();
}
 
[[nodiscard]] inline char printf_notation(real_format_notation notation)
{
  switch (notation)
  {
    case real_format_notation::fixed:
      return 'f';
    case real_format_notation::scientific:
      return 'e';
    case real_format_notation::general:
      return 'g';
  }
  return 'g';
}
 
#if UNI20_HAS_MPREAL
inline bool scalar_signbit(mpreal const& value) { return uni20::signbit(value); }
#endif
 
template <typename T> [[nodiscard]] bool scalar_signbit(T value)
{
  return std::signbit(static_cast<long double>(value));
}
 
template <typename T> [[nodiscard]] std::string fmt_real(T value, scalar_format_options const& options)
{
  int const precision = effective_precision<T>(options);
  switch (options.notation)
  {
    case real_format_notation::fixed:
      return fmt::format("{:.{}f}", value, precision);
    case real_format_notation::scientific:
      return fmt::format("{:.{}e}", value, precision);
    case real_format_notation::general:
      return fmt::format("{:.{}g}", value, precision);
  }
  return fmt::format("{}", value);
}
 
#if UNI20_HAS_FLOAT128 && UNI20_FLOAT128_PROVIDER_MPLAPACK
[[nodiscard]] inline std::string format_float128(uni20::float128 value, scalar_format_options const& options)
{
#if defined(MPLAPACK_BINARY128_MODE) && (MPLAPACK_BINARY128_MODE == MPLAPACK_BINARY128_MODE_FLOAT128)
  int const precision = effective_precision<uni20::float128>(options);
  char const notation = printf_notation(options.notation);
  std::string const format = fmt::format("%.{}{}", precision, notation);
 
  char fixed_buffer[128] = {};
  int const fixed_count = ::strfromf128(fixed_buffer, sizeof(fixed_buffer), format.c_str(), value);
  if (fixed_count < 0)
  {
    throw std::runtime_error("failed to format uni20::float128 value");
  }
  if (static_cast<std::size_t>(fixed_count) < sizeof(fixed_buffer))
  {
    return std::string(fixed_buffer, static_cast<std::size_t>(fixed_count));
  }
 
  std::string dynamic_buffer(static_cast<std::size_t>(fixed_count) + 1, '\0');
  int const dynamic_count = ::strfromf128(dynamic_buffer.data(), dynamic_buffer.size(), format.c_str(), value);
  if (dynamic_count < 0)
  {
    throw std::runtime_error("failed to format uni20::float128 value");
  }
  dynamic_buffer.resize(static_cast<std::size_t>(dynamic_count));
  return dynamic_buffer;
#else
  // This diagnostic fallback does not require binary128 numeric_limits, which
  // older standard libraries may omit. It is not admitted by DataTableValue.
  (void)options;
  return fmt::format("{}", static_cast<long double>(value));
#endif
}
 
[[nodiscard]] inline uni20::float128 parse_float128(std::string_view text)
{
  std::string buffer(text);
  char* end = nullptr;
  errno = 0;
#if defined(MPLAPACK_BINARY128_MODE) && (MPLAPACK_BINARY128_MODE == MPLAPACK_BINARY128_MODE_FLOAT128)
  uni20::float128 const parsed = ::strtof128(buffer.c_str(), &end);
#else
  long double const parsed = std::strtold(buffer.c_str(), &end);
#endif
  if (end == buffer.c_str() || *end != '\0')
  {
    throw std::invalid_argument(fmt::format("invalid float128 value: {}", text));
  }
  if (errno == ERANGE)
  {
    throw std::out_of_range(fmt::format("float128 value out of range: {}", text));
  }
  return static_cast<uni20::float128>(parsed);
}
#endif
 
template <typename T> [[nodiscard]] T parse_builtin_real(std::string_view text)
{
  std::string buffer(text);
  char* end = nullptr;
  errno = 0;
 
  if constexpr (std::same_as<T, float>)
  {
    float const parsed = std::strtof(buffer.c_str(), &end);
    if (end == buffer.c_str() || *end != '\0')
    {
      throw std::invalid_argument(fmt::format("invalid real value: {}", text));
    }
    if (errno == ERANGE)
    {
      throw std::out_of_range(fmt::format("real value out of range: {}", text));
    }
    return parsed;
  }
  else if constexpr (std::same_as<T, double>)
  {
    double const parsed = std::strtod(buffer.c_str(), &end);
    if (end == buffer.c_str() || *end != '\0')
    {
      throw std::invalid_argument(fmt::format("invalid real value: {}", text));
    }
    if (errno == ERANGE)
    {
      throw std::out_of_range(fmt::format("real value out of range: {}", text));
    }
    return parsed;
  }
  else
  {
    long double const parsed = std::strtold(buffer.c_str(), &end);
    if (end == buffer.c_str() || *end != '\0')
    {
      throw std::invalid_argument(fmt::format("invalid real value: {}", text));
    }
    if (errno == ERANGE)
    {
      throw std::out_of_range(fmt::format("real value out of range: {}", text));
    }
    return static_cast<T>(parsed);
  }
}
} // namespace detail
 
template <Real T> [[nodiscard]] std::string format_real(T value, scalar_format_options const& options = {})
{
  using real_type = std::remove_cv_t<T>;
#if UNI20_HAS_MPREAL
  if constexpr (std::same_as<real_type, mpreal>)
  {
    if (options.normalize_negative_zero && value == 0) value = mpreal(value.precision());
    return detail::format_mpreal(value, options);
  }
  else
#endif
  {
    if (options.normalize_negative_zero && value == real_type{})
    {
      value = real_type{};
    }
 
#if UNI20_HAS_FLOAT128
    if constexpr (std::same_as<real_type, uni20::float128>)
    {
      return detail::format_float128(value, options);
    }
    else
#endif
    {
      return detail::fmt_real(value, options);
    }
  }
}
 
template <Complex T> [[nodiscard]] std::string format_complex(T const& value, scalar_format_options const& options = {})
{
  using real_type = make_real_t<T>;
  real_type imag = value.imag();
  bool const negative_imag = detail::scalar_signbit(imag) && !(options.normalize_negative_zero && imag == 0);
  if (negative_imag)
  {
    imag = -imag;
  }
  else if (options.normalize_negative_zero && imag == 0 && detail::scalar_signbit(imag))
  {
    imag = -imag;
  }
 
  return format_real(value.real(), options) + (negative_imag ? "-" : "+") + format_real(imag, options) +
         std::string(options.imaginary_unit);
}
 
template <RealOrComplex T>
[[nodiscard]] std::string format_scalar(T const& value, scalar_format_options const& options = {})
{
  if constexpr (Real<T>)
  {
    return format_real(value, options);
  }
  else
  {
    return format_complex(value, options);
  }
}
 
template <Real T>
#if UNI20_HAS_MPREAL
  requires(!std::same_as<std::remove_cv_t<T>, mpreal>)
#endif
[[nodiscard]] T parse_real(std::string_view text)
{
  using real_type = std::remove_cv_t<T>;
#if UNI20_HAS_FLOAT128
  if constexpr (std::same_as<real_type, uni20::float128>)
  {
    return detail::parse_float128(text);
  }
  else
#endif
  {
    return detail::parse_builtin_real<real_type>(text);
  }
}
 
#if UNI20_HAS_MPREAL
template <Real T>
  requires std::same_as<std::remove_cv_t<T>, mpreal>
[[nodiscard]] T parse_real(std::string_view text, Precision precision)
{
  return mpreal(text, precision);
}
#endif
 
template <Real T> std::istream& read_real(std::istream& stream, T& value)
{
  std::string token;
  stream >> token;
  if (!stream)
  {
    return stream;
  }
 
  try
  {
#if UNI20_HAS_MPREAL
    if constexpr (std::same_as<std::remove_cv_t<T>, mpreal>)
      value = parse_real<T>(token, value.precision());
    else
#endif
      value = parse_real<T>(token);
  }
  catch (...)
  {
    stream.setstate(std::ios::failbit);
  }
  return stream;
}
 
} // namespace uni20
 
#if UNI20_HAS_FLOAT128
template <> struct fmt::formatter<uni20::float128>
{
    constexpr auto parse(format_parse_context& ctx)
    {
      auto it = ctx.begin();
      if (it != ctx.end() && *it != '}')
      {
        throw fmt::format_error("uni20::float128 formatter currently supports only default '{}'");
      }
      return it;
    }
 
    template <typename FormatContext> auto format(uni20::float128 value, FormatContext& ctx) const
    {
      return fmt::format_to(ctx.out(), "{}", uni20::format_real(value));
    }
};
 
template <> struct fmt::formatter<uni20::complex<uni20::float128>>
{
    constexpr auto parse(format_parse_context& ctx)
    {
      auto it = ctx.begin();
      if (it != ctx.end() && *it != '}')
      {
        throw fmt::format_error("uni20::complex<uni20::float128> formatter currently supports only default '{}'");
      }
      return it;
    }
 
    template <typename FormatContext>
    auto format(uni20::complex<uni20::float128> const& value, FormatContext& ctx) const
    {
      return fmt::format_to(ctx.out(), "{}", uni20::format_complex(value));
    }
};
#endif