// SPDX-FileCopyrightText: The Eigen Authors
// SPDX-License-Identifier: MPL-2.0
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <limits>
#include <mpfr.h>
int main() {
  const float pairs[][2] = {{0.86628270149230957f, -0.24927544593811035f},
                            {0.27852725982666016f, 0.021227836608886719f},
                            {-0.39414405822753906f, 0.74842500686645508f},
                            {-0.19227957725524902f, 0.89961457252502441f},
                            {-0.77999615669250488f, 0.41786742210388184f},
                            {0.81306219100952148f, 0.16147041320800781f},
                            {-0.35762596130371094f, -0.36932516098022461f},
                            {-0.35938382148742676f, 0.54462099075317383f},
                            {0.78412365913391113f, -0.25450515747070312f},
                            {0.64735150337219238f, -0.17155909538269043f},
                            {-0.020654678344726562f, -0.22532916069030762f},
                            {-0.82065844535827637f, -0.90036630630493164f},
                            {0.084037542343139648f, 0.73479199409484863f},
                            {-0.67060613632202148f, 0.38916158676147461f},
                            {0.31683993339538574f, 0.56197762489318848f},
                            {0.87025904655456543f, 0.6010749340057373f}};
  mpfr_t x, y, product, sum, scale;
  mpfr_inits2(256, x, y, product, sum, scale, (mpfr_ptr)0);
  mpfr_set_zero(sum, 0);
  mpfr_set_zero(scale, 0);
  for (const auto &pair : pairs) {
    mpfr_set_flt(x, pair[0], MPFR_RNDN);
    mpfr_set_flt(y, pair[1], MPFR_RNDN);
    mpfr_mul(product, x, y, MPFR_RNDN);
    mpfr_add(sum, sum, product, MPFR_RNDN);
    mpfr_abs(product, product, MPFR_RNDN);
    mpfr_add(scale, scale, product, MPFR_RNDN);
  }
  const double reference = mpfr_get_d(sum, MPFR_RNDN),
               magnitude = mpfr_get_d(scale, MPFR_RNDN);
  const double eigen = -0.0014439225196838379,
               openblas = -0.0014438522048294544;
  const double tolerance =
      64 * std::sqrt(16.0) * std::numeric_limits<float>::epsilon();
  mpfr_printf("MPFR 256-bit dot: %.50Rg\n", sum);
  std::printf("sum absolute products: %.17g\nEigen absolute error: "
              "%.17g\nOpenBLAS absolute error: %.17g\npair difference: "
              "%.17g\nold allowance: %.17g\ninput-scaled allowance: %.17g\n",
              magnitude, std::abs(eigen - reference),
              std::abs(openblas - reference), std::abs(eigen - openblas),
              tolerance * std::max(std::abs(eigen), std::abs(openblas)),
              64 * std::numeric_limits<float>::epsilon() * magnitude);
  const bool old_rejects =
      std::abs(eigen - openblas) >
      tolerance * std::max(std::abs(eigen), std::abs(openblas));
  const bool both_accurate =
      std::abs(eigen - reference) <
          std::numeric_limits<float>::epsilon() * magnitude &&
      std::abs(openblas - reference) <
          std::numeric_limits<float>::epsilon() * magnitude;
  mpfr_clears(x, y, product, sum, scale, (mpfr_ptr)0);
  return old_rejects && both_accurate ? 0 : 1;
}