APFloat Constant Generator

#include <mpfr.h>

#include <bit>
#include <cinttypes>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <string>
#include <chrono>

using std::string;

string string_val(const char *str, mpfr_prec_t Precision, mpfr_rnd_t RND) {
  mpfr_t value;
  mpfr_init2(value, Precision);
  mpfr_strtofr(value, str, nullptr, 10, RND);

  // Scale to [1, 2) (exact), remembering the binary exponent
  long exp = mpfr_get_exp(value) - 1;
  if (mpfr_mul_2si(value, value, -exp, MPFR_RNDN) != 0)
    printf("Rounding occured for [%s]\n", str);

  // 1 leading hex digit + enough fraction digits for the remaining (Precision - 1) bits
  size_t ndigits = 1 + (Precision - 1 + 3) / 4;
  mpfr_exp_t hex_exp;
  char *digits = mpfr_get_str(nullptr, &hex_exp, 16, ndigits, value, MPFR_RNDN);
  char buf[999];
  snprintf(buf, sizeof(buf), "0x%c.%sp%+ld", digits[0], digits + 1, exp);

  mpfr_free_str(digits);
  mpfr_clear(value);
  return buf;
}

// Rounds the decimal in `str` to the nearest double-double under RND.
//
// hi is always RN(x): a canonical double-double requires RN(hi + lo) == hi, so
// hi is fixed by x and only lo is free. The value increases with lo, so a
// directed mode on the value is the same directed mode on the residual, except
// for the two sign-relative modes, which have to be resolved against the sign
// of x first.
//
// The pair spans more than the 106 bits of a fixed-width significand, because
// lo carries its own exponent; that is why rounding through a flat 106-bit
// format loses the last bit of lo.
void double_double_val(const char *str, mpfr_rnd_t RND, double &hi,
                       double &lo) {
  mpfr_prec_t Precision = 2112;
  mpfr_t value, hi_exact, residual;
  mpfr_init2(value, Precision);
  mpfr_init2(hi_exact, Precision);
  mpfr_init2(residual, Precision);

  mpfr_strtofr(value, str, nullptr, 10, MPFR_RNDN);

  hi = mpfr_get_d(value, MPFR_RNDN);

  // residual = value - hi, exact at this precision
  mpfr_set_d(hi_exact, hi, MPFR_RNDN);
  if (mpfr_sub(residual, value, hi_exact, MPFR_RNDN) != 0)
    printf("Rounding occured computing the residual for [%s]\n", str);

  // Toward/away from zero refer to the value, not the residual, and the
  // residual takes either sign.
  mpfr_rnd_t ResidualRND = RND;
  if (RND == MPFR_RNDZ)
    ResidualRND = mpfr_sgn(value) < 0 ? MPFR_RNDU : MPFR_RNDD;
  else if (RND == MPFR_RNDA)
    ResidualRND = mpfr_sgn(value) < 0 ? MPFR_RNDD : MPFR_RNDU;

  lo = mpfr_get_d(residual, ResidualRND);

  // APFloat wants a canonical double-double: adding lo back must not change hi.
  if (hi + lo != hi)
    printf("Not a canonical double-double for [%s]\n", str);

  mpfr_clear(residual);
  mpfr_clear(hi_exact);
  mpfr_clear(value);
}

uint64_t f64_to_u64(double x) {
  return std::bit_cast<uint64_t>(x);
}

// These literals are not sufficient for ieee_binary256
static const struct { const char *label; const char *str; } constants[] = {
  { "e",          "2.7182818284590452353602874713526624977572470936999595749669676277..." },
  { "log2e",      "1.4426950408889634073599246810018921374266459541529859341354494069..." },
  { "log10e",     "0.4342944819032518276511289189166050822943970058036665661144537831..." },
  { "pi",         "3.1415926535897932384626433832795028841971693993751058209749445923..." },
  { "inv_pi",     "0.3183098861837906715377675267450287240689192914809128974953346881..." },
  { "inv_sqrtpi", "0.5641895835477562869480794515607725858440506293289988568440857217..." },
  { "ln2",        "0.6931471805599453094172321214581765680755001343602552541206800094..." },
  { "ln10",       "2.3025850929940456840179914546843642076011014886287729760333279009..." },
  { "sqrt2",      "1.4142135623730950488016887242096980785696718753769480731766797379..." },
  { "sqrt3",      "1.7320508075688772935274463415058723669428052538103806280558069794..." },
  { "inv_sqrt3",  "0.5773502691896257645091487805019574556476017512701268760186023264..." },
  { "egamma",     "0.5772156649015328606065120900824024310421593359399235988057672348..." },
  { "phi",        "1.6180339887498948482045868343656381177203091798057628621354486227..." },
};

// struct { const char *name; mpfr_rnd_t rnd; } modes[] = {
//   { "rmNearestTiesToEven", MPFR_RNDN },
//   { "rmTowardPositive",    MPFR_RNDU },
//   { "rmTowardNegative",    MPFR_RNDD },
//   { "rmTowardZero",        MPFR_RNDZ },
// };

void print_literals(mpfr_prec_t Precision) {
  for (auto &c : constants) {
    string Down = string_val(c.str, Precision, MPFR_RNDD);
    string Near = string_val(c.str, Precision, MPFR_RNDN);
    string Up = string_val(c.str, Precision, MPFR_RNDU);
    printf("    {APFloat::MathConstant::%s, \"%s\",\n", c.label, c.label);
    printf("     \"%s\",\n", Down.c_str());
    printf("     \"%s\",\n", Near.c_str());
    printf("     \"%s\"},\n", Up.c_str());
  }
  printf("\n");
}

void print_PPCDoubleDouble_literals() {
  for (auto &c : constants) {
    double Down_hi, Down_lo;
    double Near_hi, Near_lo;
    double Up_hi, Up_lo;
    double_double_val(c.str, MPFR_RNDD, Down_hi, Down_lo);
    double_double_val(c.str, MPFR_RNDN, Near_hi, Near_lo);
    double_double_val(c.str, MPFR_RNDU, Up_hi, Up_lo);
    printf("    {APFloat::MathConstant::%s, \"%s\",\n", c.label, c.label);
    printf("     \"0x%" PRIx64 "u, 0x%" PRIx64 "u,\n", f64_to_u64(Down_hi), f64_to_u64(Down_lo));
    printf("     \"0x%" PRIx64 "u, 0x%" PRIx64 "u,\n", f64_to_u64(Near_hi), f64_to_u64(Near_lo));
    printf("     \"0x%" PRIx64 "u, 0x%" PRIx64 "u},\n", f64_to_u64(Up_hi), f64_to_u64(Up_lo));
  }
  printf("\n");
}

int main(void) {

#if 1
  mpfr_prec_t Precision = 53;
  print_literals(Precision);
#else
  print_PPCDoubleDouble_literals();
#endif

  return 0;
}