xref: /src/external/lgpl3/mpc/dist/tests/tmul.c

/* tmul -- test file for mpc_mul.

Copyright (C) 2002, 2005, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2020 INRIA

This file is part of GNU MPC.

GNU MPC is free software; you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License as published by the
Free Software Foundation; either version 3 of the License, or (at your
option) any later version.

GNU MPC 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 Lesser General Public License for
more details.

You should have received a copy of the GNU Lesser General Public License
along with this program. If not, see http://www.gnu.org/licenses/ .
*/

#include <stdlib.h>
#ifdef TIMING
#include <sys/times.h>
#endif
#include "mpc-tests.h"

static void
cmpmul (mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
   /* computes the product of x and y with the naive and Karatsuba methods */
   /* using the rounding mode rnd and compares the results and return      */
   /* values.                                                              */
   /* In our current test suite, the real and imaginary parts of x and y   */
   /* all have the same precision, and we use this precision also for the  */
   /* result.                                                              */
{
   mpc_t z, t;
   int   inex_z, inex_t;

   mpc_init2 (z, MPC_MAX_PREC (x));
   mpc_init2 (t, MPC_MAX_PREC (x));

   inex_z = mpc_mul_naive (z, x, y, rnd);
   inex_t = mpc_mul_karatsuba (t, x, y, rnd);

   if (mpc_cmp (z, t) != 0 || inex_z != inex_t) {
      fprintf (stderr, "mul_naive and mul_karatsuba differ for rnd=(%s,%s)\n",
               mpfr_print_rnd_mode(MPC_RND_RE(rnd)),
               mpfr_print_rnd_mode(MPC_RND_IM(rnd)));
      MPC_OUT (x);
      MPC_OUT (y);
      MPC_OUT (z);
      MPC_OUT (t);
      if (inex_z != inex_t) {
         fprintf (stderr, "inex_re (z): %s\n", MPC_INEX_STR (inex_z));
         fprintf (stderr, "inex_re (t): %s\n", MPC_INEX_STR (inex_t));
      }
      exit (1);
   }

   mpc_clear (z);
   mpc_clear (t);
}


static void
testmul (long a, long b, long c, long d, mpfr_prec_t prec, mpc_rnd_t rnd)
{
  mpc_t x, y;

  mpc_init2 (x, prec);
  mpc_init2 (y, prec);

  mpc_set_si_si (x, a, b, rnd);
  mpc_set_si_si (y, c, d, rnd);

  cmpmul (x, y, rnd);

  mpc_clear (x);
  mpc_clear (y);
}


static void
check_regular (void)
{
  mpc_t x, y;
  int rnd_re, rnd_im;
  mpfr_prec_t prec;

  testmul (247, -65, -223, 416, 8, 24);
  testmul (5, -896, 5, -32, 3, 2);
  testmul (-3, -512, -1, -1, 2, 16);
  testmul (266013312, 121990769, 110585572, 116491059, 27, 0);
  testmul (170, 9, 450, 251, 8, 0);
  testmul (768, 85, 169, 440, 8, 16);
  testmul (145, 1816, 848, 169, 8, 24);

  mpc_init2 (x, 1000);
  mpc_init2 (y, 1000);

  /* Bug 20081114: mpc_mul_karatsuba returned wrong inexact value for
     imaginary part */
  mpc_set_prec (x, 7);
  mpc_set_prec (y, 7);
  mpfr_set_str (mpc_realref (x), "0xB4p+733", 16, MPFR_RNDN);
  mpfr_set_str (mpc_imagref (x), "0x90p+244", 16, MPFR_RNDN);
  mpfr_set_str (mpc_realref (y), "0xECp-146", 16, MPFR_RNDN);
  mpfr_set_str (mpc_imagref (y), "0xACp-471", 16, MPFR_RNDN);
  cmpmul (x, y, MPC_RNDNN);
  mpfr_set_str (mpc_realref (x), "0xB4p+733", 16, MPFR_RNDN);
  mpfr_set_str (mpc_imagref (x), "0x90p+244", 16, MPFR_RNDN);
  mpfr_set_str (mpc_realref (y), "0xACp-471", 16, MPFR_RNDN);
  mpfr_set_str (mpc_imagref (y), "-0xECp-146", 16, MPFR_RNDN);
  cmpmul (x, y, MPC_RNDNN);

  for (prec = 2; prec < 1000; prec = (mpfr_prec_t) (prec * 1.1 + 1))
    {
      mpc_set_prec (x, prec);
      mpc_set_prec (y, prec);

      test_default_random (x, -1024, 1024, 128, 0);
      test_default_random (y, -1024, 1024, 128, 0);

      for (rnd_re = 0; rnd_re < 4; rnd_re ++)
        for (rnd_im = 0; rnd_im < 4; rnd_im ++)
          cmpmul (x, y, MPC_RND (rnd_re, rnd_im));
    }

  mpc_clear (x);
  mpc_clear (y);
}

static void
bug20200206 (void)
{
  mpfr_exp_t emin = mpfr_get_emin ();
  mpc_t x, y, z;

  mpfr_set_emin (-1073);
  mpc_init2 (x, 53);
  mpc_init2 (y, 53);
  mpc_init2 (z, 53);
  mpfr_set_d (mpc_realref (x), -6.0344722345057644e-272, MPFR_RNDN);
  mpfr_set_d (mpc_imagref (x), -4.8536770224196353e-204, MPFR_RNDN);
  mpfr_set_d (mpc_realref (y), 1.3834775731431992e-246, MPFR_RNDN);
  mpfr_set_d (mpc_imagref (y), 2.9246270396940562e-124, MPFR_RNDN);
  mpc_mul (z, x, y, MPC_RNDNN);
  if (mpfr_regular_p (mpc_realref (z)) &&
      mpfr_get_exp (mpc_realref (z)) < -1073)
    {
      printf ("Error, mpc_mul returns an out-of-range exponent:\n");
      mpfr_dump (mpc_realref (z));
      printf ("Bug most probably in MPFR, please upgrade to MPFR 4.1.0 or later\n");
      exit (1);
    }
  mpc_clear (x);
  mpc_clear (y);
  mpc_clear (z);
  mpfr_set_emin (emin);
}

#ifdef TIMING
static void
timemul (void)
{
  /* measures the time needed with different precisions for naive and */
  /* Karatsuba multiplication                                         */

  mpc_t             x, y, z;
  unsigned long int i, j;
  const unsigned long int tests = 10000;
  struct tms        time_old, time_new;
  double            passed1, passed2;

  mpc_init (x);
  mpc_init (y);
  mpc_init_set_ui_ui (z, 1, 0, MPC_RNDNN);

  for (i = 1; i < 50; i++)
    {
      mpc_set_prec (x, i * BITS_PER_MP_LIMB);
      mpc_set_prec (y, i * BITS_PER_MP_LIMB);
      mpc_set_prec (z, i * BITS_PER_MP_LIMB);
      test_default_random (x, -1, 1, 128, 25);
      test_default_random (y, -1, 1, 128, 25);

      times (&time_old);
      for (j = 0; j < tests; j++)
        mpc_mul_naive (z, x, y, MPC_RNDNN);
      times (&time_new);
      passed1 = ((double) (time_new.tms_utime - time_old.tms_utime)) / 100;

      times (&time_old);
      for (j = 0; j < tests; j++)
        mpc_mul_karatsuba (z, x, y, MPC_RNDNN);
      times (&time_new);
      passed2 = ((double) (time_new.tms_utime - time_old.tms_utime)) / 100;

      printf ("Time for %3li limbs naive/Karatsuba: %5.2f %5.2f\n", i,
              passed1, passed2);
    }

  mpc_clear (x);
  mpc_clear (y);
  mpc_clear (z);
}
#endif

#define MPC_FUNCTION_CALL                                               \
  P[0].mpc_inex = mpc_mul (P[1].mpc, P[2].mpc, P[3].mpc, P[4].mpc_rnd)
#define MPC_FUNCTION_CALL_SYMMETRIC                                     \
  P[0].mpc_inex = mpc_mul (P[1].mpc, P[3].mpc, P[2].mpc, P[4].mpc_rnd)
#define MPC_FUNCTION_CALL_REUSE_OP1                                     \
  P[0].mpc_inex = mpc_mul (P[1].mpc, P[1].mpc, P[3].mpc, P[4].mpc_rnd)
#define MPC_FUNCTION_CALL_REUSE_OP2                                     \
  P[0].mpc_inex = mpc_mul (P[1].mpc, P[2].mpc, P[1].mpc, P[4].mpc_rnd)

#include "data_check.tpl"
#include "tgeneric.tpl"

static void
bug20221130 (void)
{
  mpc_t b, c_conj, res, ref;
  mpfr_prec_t prec;
  mpc_init2 (b, 20);
  mpc_init2 (c_conj, 20);
  mpc_init2 (res, 2);
  mpc_init2 (ref, 5);
  mpc_set_str (b, "(0x1p+0 0x2p+0)", 16, MPC_RNDNN);
  mpc_set_str (c_conj, "(-0xap+0 0x1.4p+4)", 16, MPC_RNDNN);
  mpc_set_str (ref, "(-0x3.2p+4 0x0p+0)", 16, MPC_RNDNN);
  for (prec = 5; prec <= 2000; prec++)
  {
    mpc_set_prec (res, prec);
    mpc_mul (res, b, c_conj, MPC_RNDZZ);
    if (mpc_cmp (res, ref) != 0 || mpfr_signbit (mpc_imagref (res)))
    {
      printf ("Error in bug20221130 for prec=%lu\n", prec);
      mpfr_printf ("expected (%Ra %Ra)\n", mpc_realref (ref), mpc_imagref (ref));
      mpfr_printf ("got      (%Ra %Ra)\n", mpc_realref (res), mpc_imagref (res));
      exit (1);
    }
  }
  mpc_clear (b);
  mpc_clear (c_conj);
  mpc_clear (res);
  mpc_clear (ref);
}

int
main (void)
{
  test_start ();

#ifdef TIMING
  timemul ();
#endif

  bug20200206 ();
  bug20221130 ();
  check_regular ();

  data_check_template ("mul.dsc", "mul.dat");

  tgeneric_template ("mul.dsc", 2, 4096, 41, 1024);

  test_end ();

  return 0;
}