/* Test sipe.h against GNU MPFR
 *
 * Usage: test-sipe <emax> [ <prec_max> ]
 *
 * All values whose exponent of the integral significand is <= emax
 * in absolute value will be tested (including 0). All precisions
 * from precision 2 (to prec_max, if provided) will be tested.
 * FIXME: The condition on the exponent seems wrong as one starts with
 * 2^(-emax), so that the quantum exponent is not -emax; moreover, it
 * seems that the highest quantum exponent is emax - 1. Check r63219,
 * and change the code or his comment.
 *
 * Copyright (c) 2011-2022 Vincent Lefevre <vincent@vinc17.net>
 * Arenaire/AriC projects, INRIA / LIP / ENS-Lyon, France.
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 3
 * of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, see https://www.gnu.org/licenses/
 * or write to the Free Software Foundation, Inc., 51 Franklin St,
 * Fifth Floor, Boston, MA 02110-1301, USA.
 */

#define SVNID "$Id: test-sipe.c 146217 2022-02-28 17:03:13Z vinc17/cventin $"

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <mpfr.h>

#ifndef SIPE_USE_FMA
#define SIPE_USE_FMA 1
#endif

#include "sipe.h"

#define STRINGIZE(S) #S
#define MAKE_STR(S) STRINGIZE(S)

#if SIPE_TRUNC
#define RND MPFR_RNDZ
#else
#define RND MPFR_RNDN
#endif

static mpfr_rnd_t rnd = RND;

static long emax;

/* Syntax check of SIPE_ROUND (during compilation). */
static void check_round_macro (int prec)
{
#ifdef SIPE_FLOAT
  sipe_t r = 1;
#else
  sipe_t r = { 1, 0 };
#endif

  if (1)
    SIPE_ROUND (r, prec);
  else
    SIPE_ROUND (r, prec);

  (void) r;
}

static sipe_t first (mpfr_ptr u, int prec)
{
  sipe_t r = sipe_set_si_2exp (1, -emax, prec);
#ifdef SIPE_FLOAT
  sipe_t s = ldexp (1.0, -emax);
#else
  sipe_t s = { 1, -emax };
  SIPE_ROUND (s, prec);
#endif

  if (! sipe_eq (r, s))
    {
      putchar ('\n');
      fflush (stdout);
      fprintf (stderr, "test-sipe: error for sipe_set_si_2exp on (1,%ld)\n",
               -emax);
      exit (1);
    }

  mpfr_set_ui_2exp (u, 1, -emax, MPFR_RNDN);
  return r;
}

static inline void s2mpfr (mpfr_ptr sm, sipe_t rs, int prec)
{
  mpfr_set_si_2exp (sm, sipe_intm (rs, prec),
                    sipe_qexp (rs, prec), MPFR_RNDN);
}

#define LOOP(U,BODY)                                             \
  do                                                             \
    {                                                            \
      U##s = first (U##m, prec);                                 \
      while (1)                                                  \
        {                                                        \
          for (U##n = ! sipe_zero_p (U##s); U##n >= 0; U##n--)   \
            { BODY;                                              \
              mpfr_neg (U##m, U##m, MPFR_RNDN);                  \
              U##s = sipe_neg (U##s, prec); }                    \
          if (sipe_zero_p (U##s))                                \
            break;                                               \
          mpfr_nextabove (U##m);                                 \
          U##s = sipe_nextabove (U##s, prec);                    \
          if (sipe_qexp (U##s, prec) < emax)                     \
            continue;                                            \
          mpfr_set_zero (U##m, 0);                               \
          U##s = sipe_set_si (0, prec);                          \
        }                                                        \
    }                                                            \
  while (0)

int main (int argc, char **argv)
{
  char *end;
  int prec;
  long prec_max = SIPE_PREC_MAX;
  long ctr_nextabove = 0, ctr_nextbelow = 0, ctr_nexttozero = 0,
    ctr_neg = 0, ctr_abs = 0, ctr_eqne0 = 0, ctr_legt0 = 0, ctr_gelt0 = 0,
    ctr_add = 0, ctr_sub = 0, ctr_mul = 0, ctr_fma = 0, ctr_fms = 0,
    ctr_2mul = 0, ctr_add_si = 0, ctr_sub_si = 0, ctr_mul_si = 0,
    ctr_eqne = 0, ctr_legt = 0, ctr_gelt = 0, ctr_min = 0, ctr_max = 0;

#include "version-info.h"

  if (argc < 2 || argc > 3)
    {
      fprintf (stderr, "Usage: test-sipe <emax> [ <prec_max> ]\n");
      exit (1);
    }

  emax = strtol (argv[1], &end, 10);
  if (*end != '\0' || emax < 1 || emax > INT_MAX)
    {
      fprintf (stderr, "test-sipe: bad emax value\n");
      exit (1);
    }

  if (argc == 3)
    {
      prec_max = strtol (argv[2], &end, 10);
      if (*end != '\0' || prec_max < 2 || prec_max > SIPE_PREC_MAX)
        {
          fprintf (stderr, "test-sipe: bad prec_max value\n");
          exit (1);
        }
    }

  for (prec = 2; prec <= prec_max; prec++)
    {
      mpfr_t xm, ym, zm, rm, sm;
      sipe_t xs, ys, zs, rs, ss;
      int xn, yn;
#if SIPE_USE_FMA
      int zn;
#endif
      int bm, bsf, bst;
      sipe_int_t i, imax;
      clock_t t;

      printf ("Testing precision %d with emax = %ld...", prec, emax);
      fflush (stdout);
      t = clock ();
      check_round_macro (prec);
      imax = (sipe_int_t) 1 << prec;
      mpfr_inits2 (prec, xm, ym, zm, rm, sm, (mpfr_ptr) 0);
      LOOP (x,
        {
          if (! sipe_zero_p (xs))
            {
              ++ctr_nextabove;
              rs = sipe_nextabove (xs, prec);
              s2mpfr (sm, rs, prec);
              mpfr_set (rm, xm, MPFR_RNDN);
              mpfr_nextabove (rm);
              if (! mpfr_equal_p (rm, sm))
                {
                  putchar ('\n');
                  fflush (stdout);
                  mpfr_fprintf
                    (stderr,
                     "test-sipe: error for nextabove on %Rb\n"
                     "expected %Rb\n"
                     "     got %Rb\n",
                     xm, rm, sm);
                  exit (1);
                }
              ++ctr_nextbelow;
              rs = sipe_nextbelow (xs, prec);
              s2mpfr (sm, rs, prec);
              mpfr_set (rm, xm, MPFR_RNDN);
              mpfr_nextbelow (rm);
              if (! mpfr_equal_p (rm, sm))
                {
                  putchar ('\n');
                  fflush (stdout);
                  mpfr_fprintf
                    (stderr,
                     "test-sipe: error for nextbelow on %Rb\n"
                     "expected %Rb\n"
                     "     got %Rb\n",
                     xm, rm, sm);
                  exit (1);
                }
              ++ctr_nexttozero;
              rs = sipe_nexttozero (xs, prec);
              s2mpfr (sm, rs, prec);
              mpfr_set (rm, xm, MPFR_RNDN);
              mpfr_set_zero (ym, 1);
              mpfr_nexttoward (rm, ym);
              if (! mpfr_equal_p (rm, sm))
                {
                  putchar ('\n');
                  fflush (stdout);
                  mpfr_fprintf
                    (stderr,
                     "test-sipe: error for nexttozero on %Rb\n"
                     "expected %Rb\n"
                     "     got %Rb\n",
                     xm, rm, sm);
                  exit (1);
                }
            }  /* ! sipe_zero_p (xs) */
          ++ctr_neg;
          rs = sipe_neg (xs, prec);
          s2mpfr (sm, rs, prec);
          mpfr_neg (rm, xm, rnd);
          if (! mpfr_equal_p (rm, sm))
            {
              putchar ('\n');
              fflush (stdout);
              mpfr_fprintf
                (stderr,
                 "test-sipe: error for neg on %Rb\n"
                 "expected %Rb\n"
                 "     got %Rb\n",
                 xm, rm, sm);
              exit (1);
            }
          ++ctr_abs;
          rs = sipe_abs (xs, prec);
          s2mpfr (sm, rs, prec);
          mpfr_abs (rm, xm, rnd);
          if (! mpfr_equal_p (rm, sm))
            {
              putchar ('\n');
              fflush (stdout);
              mpfr_fprintf
                (stderr,
                 "test-sipe: error for abs on %Rb\n"
                 "expected %Rb\n"
                 "     got %Rb\n",
                 xm, rm, sm);
              exit (1);
            }
          ++ctr_eqne0;
          bsf = sipe_eq0 (xs);
          bst = sipe_ne0 (xs);
          bm = ! mpfr_zero_p (xm);
          if (bsf == bm || bst != bm || sipe_zero_p (xs) != bsf)
            {
              putchar ('\n');
              fflush (stdout);
              mpfr_fprintf
                (stderr,
                 "test-sipe: error for sipe_eq0/ne0 on %Rb\n"
                 "(%d / %d / %d)\n", xm, bm, bsf, bst);
              exit (1);
            }
          ++ctr_legt0;
          bsf = sipe_le0 (xs);
          bst = sipe_gt0 (xs);
          bm = mpfr_sgn (xm) > 0;
          if (bsf == bm || bst != bm)
            {
              putchar ('\n');
              fflush (stdout);
              mpfr_fprintf
                (stderr,
                 "test-sipe: error for sipe_le0/gt0 on %Rb\n"
                 "(%d / %d / %d)\n", xm, bm, bsf, bst);
              exit (1);
            }
          ++ctr_gelt0;
          bsf = sipe_ge0 (xs);
          bst = sipe_lt0 (xs);
          bm = mpfr_sgn (xm) < 0;
          if (bsf == bm || bst != bm)
            {
              putchar ('\n');
              fflush (stdout);
              mpfr_fprintf
                (stderr,
                 "test-sipe: error for sipe_ge0/lt0 on %Rb\n"
                 "(%d / %d / %d)\n", xm, bm, bsf, bst);
              exit (1);
            }
          LOOP (y,
            {
              ++ctr_add;
              rs = sipe_add (xs, ys, prec);
              s2mpfr (sm, rs, prec);
              mpfr_add (rm, xm, ym, rnd);
              if (! mpfr_equal_p (rm, sm))
                {
                  putchar ('\n');
                  fflush (stdout);
                  mpfr_fprintf
                    (stderr,
                     "test-sipe: error for add on %Rb / %Rb\n"
                     "expected %Rb\n"
                     "     got %Rb\n",
                     xm, ym, rm, sm);
                  exit (1);
                }
              ++ctr_sub;
              rs = sipe_sub (xs, ys, prec);
              s2mpfr (sm, rs, prec);
              mpfr_sub (rm, xm, ym, rnd);
              if (! mpfr_equal_p (rm, sm))
                {
                  putchar ('\n');
                  fflush (stdout);
                  mpfr_fprintf
                    (stderr,
                     "test-sipe: error for sub on %Rb / %Rb\n"
                     "expected %Rb\n"
                     "     got %Rb\n",
                     xm, ym, rm, sm);
                  exit (1);
                }
              ++ctr_mul;
              rs = sipe_mul (xs, ys, prec);
              s2mpfr (sm, rs, prec);
              mpfr_mul (rm, xm, ym, rnd);
              if (! mpfr_equal_p (rm, sm))
                {
                  putchar ('\n');
                  fflush (stdout);
                  mpfr_fprintf
                    (stderr,
                     "test-sipe: error for mul on %Rb / %Rb\n"
                     "expected %Rb\n"
                     "     got %Rb\n",
                     xm, ym, rm, sm);
                  exit (1);
                }
              ++ctr_2mul;
              SIPE_2MUL (zs, ss, xs, ys, prec);
              s2mpfr (sm, ss, prec);
              mpfr_fms (rm, xm, ym, rm, rnd);
              if (! sipe_eq (rs, zs) || ! mpfr_equal_p (rm, sm))
                {
                  putchar ('\n');
                  fflush (stdout);
                  mpfr_fprintf
                    (stderr, "test-sipe: error for SIPE_2MUL on %Rb / %Rb\n"
                     "expected (", xm, ym);
                  sipe_outbin (stderr, rs, prec);
                  mpfr_fprintf (stderr, ",%Rb)\n     got (", rm);
                  sipe_outbin (stderr, zs, prec);
                  mpfr_fprintf (stderr, ",%Rb)\n", sm);
                  exit (1);
                }
#if SIPE_USE_FMA
              LOOP (z,
                {
                  ++ctr_fma;
                  rs = sipe_fma (xs, ys, zs, prec);
                  s2mpfr (sm, rs, prec);
                  mpfr_fma (rm, xm, ym, zm, rnd);
                  if (! mpfr_equal_p (rm, sm))
                    {
                      putchar ('\n');
                      fflush (stdout);
                      mpfr_fprintf
                        (stderr,
                         "test-sipe: error for fma on %Rb / %Rb / %Rb\n"
                         "expected %Rb\n"
                         "     got %Rb\n",
                         xm, ym, zm, rm, sm);
                      exit (1);
                    }
                  ++ctr_fms;
                  rs = sipe_fms (xs, ys, zs, prec);
                  s2mpfr (sm, rs, prec);
                  mpfr_fms (rm, xm, ym, zm, rnd);
                  if (! mpfr_equal_p (rm, sm))
                    {
                      putchar ('\n');
                      fflush (stdout);
                      mpfr_fprintf
                        (stderr,
                         "test-sipe: error for fms on %Rb / %Rb / %Rb\n"
                         "expected %Rb\n"
                         "     got %Rb\n",
                         xm, ym, zm, rm, sm);
                      exit (1);
                    }
                } );  /* LOOP (z, ...) */
#endif
              if (prec <= 3)
                {
                  ++ctr_eqne;
                  bsf = sipe_eq (xs, ys);
                  bst = sipe_ne (xs, ys);
                  bm = ! mpfr_equal_p (xm, ym);
                  if (bsf == bm || bst != bm)
                    {
                      putchar ('\n');
                      fflush (stdout);
                      mpfr_fprintf
                        (stderr,
                         "test-sipe: error for sipe_eq/ne on %Rb / %Rb\n"
                         "(%d / %d / %d)\n", xm, ym, bm, bsf, bst);
                      exit (1);
                    }
                  ++ctr_legt;
                  bsf = sipe_le (xs, ys);
                  bst = sipe_gt (xs, ys);
                  bm = ! mpfr_lessequal_p (xm, ym);
                  if (bsf == bm || bst != bm)
                    {
                      putchar ('\n');
                      fflush (stdout);
                      mpfr_fprintf
                        (stderr,
                         "test-sipe: error for sipe_le/gt on %Rb / %Rb\n"
                         "(%d / %d / %d)\n", xm, ym, bm, bsf, bst);
                      exit (1);
                    }
                  ++ctr_gelt;
                  bsf = sipe_ge (xs, ys);
                  bst = sipe_lt (xs, ys);
                  bm = ! mpfr_greaterequal_p (xm, ym);
                  if (bsf == bm || bst != bm)
                    {
                      putchar ('\n');
                      fflush (stdout);
                      mpfr_fprintf
                        (stderr,
                         "test-sipe: error for sipe_ge/lt on %Rb / %Rb\n"
                         "(%d / %d / %d)\n", xm, ym, bm, bsf, bst);
                      exit (1);
                    }
                  ++ctr_min;
                  rs = sipe_min (xs, ys);
                  s2mpfr (sm, rs, prec);
                  mpfr_min (rm, xm, ym, rnd);
                  if (! mpfr_equal_p (rm, sm))
                    {
                      putchar ('\n');
                      fflush (stdout);
                      mpfr_fprintf
                        (stderr,
                         "test-sipe: error for min on %Rb / %Rb\n"
                         "expected %Rb\n"
                         "     got %Rb\n",
                         xm, ym, rm, sm);
                      exit (1);
                    }
                  ++ctr_max;
                  rs = sipe_max (xs, ys);
                  s2mpfr (sm, rs, prec);
                  mpfr_max (rm, xm, ym, rnd);
                  if (! mpfr_equal_p (rm, sm))
                    {
                      putchar ('\n');
                      fflush (stdout);
                      mpfr_fprintf
                        (stderr,
                         "test-sipe: error for max on %Rb / %Rb\n"
                         "expected %Rb\n"
                         "     got %Rb\n",
                         xm, ym, rm, sm);
                      exit (1);
                    }
                }  /* prec <= 3 */
            } );  /* LOOP (y, ...) */
          for (i = -imax; i <= imax; i++)
            {
              ++ctr_add_si;
              rs = sipe_add_si (xs, i, prec);
              s2mpfr (sm, rs, prec);
              mpfr_add_si (rm, xm, i, rnd);
              if (! mpfr_equal_p (rm, sm))
                {
                  putchar ('\n');
                  fflush (stdout);
                  mpfr_fprintf
                    (stderr,
                     "test-sipe: error for add_si on %Rb / %d\n"
                     "expected %Rb\n"
                     "     got %Rb\n",
                     xm, i, rm, sm);
                  exit (1);
                }
              ++ctr_sub_si;
              rs = sipe_sub_si (xs, i, prec);
              s2mpfr (sm, rs, prec);
              mpfr_sub_si (rm, xm, i, rnd);
              if (! mpfr_equal_p (rm, sm))
                {
                  putchar ('\n');
                  fflush (stdout);
                  mpfr_fprintf
                    (stderr,
                     "test-sipe: error for sub_si on %Rb / %d\n"
                     "expected %Rb\n"
                     "     got %Rb\n",
                     xm, i, rm, sm);
                  exit (1);
                }
              ++ctr_mul_si;
              rs = sipe_mul_si (xs, i, prec);
              s2mpfr (sm, rs, prec);
              mpfr_mul_si (rm, xm, i, rnd);
              if (! mpfr_equal_p (rm, sm))
                {
                  putchar ('\n');
                  fflush (stdout);
                  mpfr_fprintf
                    (stderr,
                     "test-sipe: error for mul_si on %Rb / %d\n"
                     "expected %Rb\n"
                     "     got %Rb\n",
                     xm, i, rm, sm);
                  exit (1);
                }
            }
        } );  /* LOOP (x, ...) */
      mpfr_clears (xm, ym, zm, rm, sm, (mpfr_ptr) 0);
      t = clock () - t;
      printf (" %.2f s\n", (double) t / CLOCKS_PER_SEC);
      fflush (stdout);
    }

  printf ("Number of tests:\n");
  printf ("%12ld sipe_nextabove\n", ctr_nextabove);
  printf ("%12ld sipe_nextbelow\n", ctr_nextbelow);
  printf ("%12ld sipe_nexttozero\n", ctr_nexttozero);
  printf ("%12ld sipe_neg\n", ctr_neg);
  printf ("%12ld sipe_abs\n", ctr_abs);
  printf ("%12ld sipe_eqne0\n", ctr_eqne0);
  printf ("%12ld sipe_legt0\n", ctr_legt0);
  printf ("%12ld sipe_gelt0\n", ctr_gelt0);
  printf ("%12ld sipe_add\n", ctr_add);
  printf ("%12ld sipe_sub\n", ctr_sub);
  printf ("%12ld sipe_mul\n", ctr_mul);
  printf ("%12ld SIPE_2MUL\n", ctr_2mul);
  printf ("%12ld sipe_fma\n", ctr_fma);
  printf ("%12ld sipe_fms\n", ctr_fms);
  printf ("%12ld sipe_eqne\n", ctr_eqne);
  printf ("%12ld sipe_legt\n", ctr_legt);
  printf ("%12ld sipe_gelt\n", ctr_gelt);
  printf ("%12ld sipe_min\n", ctr_min);
  printf ("%12ld sipe_max\n", ctr_max);
  printf ("%12ld sipe_add_si\n", ctr_add_si);
  printf ("%12ld sipe_sub_si\n", ctr_sub_si);
  printf ("%12ld sipe_mul_si\n", ctr_mul_si);

  return 0;
}