/* tsqr -- test file for mpc_sqr.

Copyright (C) 2002, 2005, 2008, 2010, 2011 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>

#include "mpc-tests.h"

static void

cmpsqr (mpc_srcptr x, mpc_rnd_t rnd)

   /* computes the square of x with the specific function or by simple     */

   /* multiplication using the rounding mode rnd and compares the results  */

   /* and return values.                                                   */

   /* In our current test suite, the real and imaginary parts of x have    */

   /* the same precision, and we use this precision also for the result.   */

   /* Furthermore, we check whether computing the square in the same       */

   /* place yields the same result.                                        */

   /* We also compute the result with four times the precision and check   */

   /* whether the rounding is correct. Error reports in this part of the   */

   /* algorithm might still be wrong, though, since there are two          */

   /* consecutive roundings.                                               */

{

  mpc_t z, t, u;

  int   inexact_z, inexact_t;

  mpc_init2 (z, MPC_MAX_PREC (x));

  mpc_init2 (t, MPC_MAX_PREC (x));

  mpc_init2 (u, 4 * MPC_MAX_PREC (x));

  inexact_z = mpc_sqr (z, x, rnd);

  inexact_t = mpc_mul (t, x, x, rnd);

  if (mpc_cmp (z, t))

    {

      fprintf (stderr, "sqr and mul differ for rnd=(%s,%s) \nx=",

               mpfr_print_rnd_mode(MPC_RND_RE(rnd)),

               mpfr_print_rnd_mode(MPC_RND_IM(rnd)));

      mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);

      fprintf (stderr, "\nmpc_sqr gives ");

      mpc_out_str (stderr, 2, 0, z, MPC_RNDNN);

      fprintf (stderr, "\nmpc_mul gives ");

      mpc_out_str (stderr, 2, 0, t, MPC_RNDNN);

      fprintf (stderr, "\n");

      exit (1);

    }

  if (inexact_z != inexact_t)

    {

      fprintf (stderr, "The return values of sqr and mul differ for rnd=(%s,%s) \nx=  ",

               mpfr_print_rnd_mode(MPC_RND_RE(rnd)),

               mpfr_print_rnd_mode(MPC_RND_IM(rnd)));

      mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);

      fprintf (stderr, "\nx^2=");

      mpc_out_str (stderr, 2, 0, z, MPC_RNDNN);

      fprintf (stderr, "\nmpc_sqr gives %i", inexact_z);

      fprintf (stderr, "\nmpc_mul gives %i", inexact_t);

      fprintf (stderr, "\n");

      exit (1);

    }

  mpc_set (t, x, MPC_RNDNN);

  inexact_t = mpc_sqr (t, t, rnd);

  if (mpc_cmp (z, t))

    {

      fprintf (stderr, "sqr and sqr in place differ for rnd=(%s,%s) \nx=",

               mpfr_print_rnd_mode(MPC_RND_RE(rnd)),

               mpfr_print_rnd_mode(MPC_RND_IM(rnd)));

      mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);

      fprintf (stderr, "\nmpc_sqr          gives ");

      mpc_out_str (stderr, 2, 0, z, MPC_RNDNN);

      fprintf (stderr, "\nmpc_sqr in place gives ");

      mpc_out_str (stderr, 2, 0, t, MPC_RNDNN);

      fprintf (stderr, "\n");

      exit (1);

    }

  if (inexact_z != inexact_t)

    {

      fprintf (stderr, "The return values of sqr and sqr in place differ for rnd=(%s,%s) \nx=  ",

               mpfr_print_rnd_mode(MPC_RND_RE(rnd)),

               mpfr_print_rnd_mode(MPC_RND_IM(rnd)));

      mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);

      fprintf (stderr, "\nx^2=");

      mpc_out_str (stderr, 2, 0, z, MPC_RNDNN);

      fprintf (stderr, "\nmpc_sqr          gives %i", inexact_z);

      fprintf (stderr, "\nmpc_sqr in place gives %i", inexact_t);

      fprintf (stderr, "\n");

      exit (1);

    }

  mpc_sqr (u, x, rnd);

  mpc_set (t, u, rnd);

  if (mpc_cmp (z, t))

    {

      fprintf (stderr, "rounding in sqr might be incorrect for rnd=(%s,%s) \nx=",

               mpfr_print_rnd_mode(MPC_RND_RE(rnd)),

               mpfr_print_rnd_mode(MPC_RND_IM(rnd)));

      mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);

      fprintf (stderr, "\nmpc_sqr                     gives ");

      mpc_out_str (stderr, 2, 0, z, MPC_RNDNN);

      fprintf (stderr, "\nmpc_sqr quadruple precision gives ");

      mpc_out_str (stderr, 2, 0, u, MPC_RNDNN);

      fprintf (stderr, "\nand is rounded to                 ");

      mpc_out_str (stderr, 2, 0, t, MPC_RNDNN);

      fprintf (stderr, "\n");

      exit (1);

    }

  mpc_clear (z);

  mpc_clear (t);

  mpc_clear (u);

}

static void

testsqr (long a, long b, mpfr_prec_t prec, mpc_rnd_t rnd)

{

  mpc_t x;

  mpc_init2 (x, prec);

  mpc_set_si_si (x, a, b, rnd);

  cmpsqr (x, rnd);

  mpc_clear (x);

}

static void

reuse_bug (void)

{

  mpc_t z1;

  /* reuse bug found by Paul Zimmermann 20081021 */

  mpc_init2 (z1, 2);

  /* RE (z1^2) overflows, IM(z^2) = -0 */

  mpfr_set_str (mpc_realref (z1), "0.11", 2, GMP_RNDN);

  mpfr_mul_2si (mpc_realref (z1), mpc_realref (z1), mpfr_get_emax (), GMP_RNDN);

  mpfr_set_ui (mpc_imagref (z1), 0, GMP_RNDN);

  mpc_conj (z1, z1, MPC_RNDNN);

  mpc_sqr (z1, z1, MPC_RNDNN);

  if (!mpfr_inf_p (mpc_realref (z1)) || mpfr_signbit (mpc_realref (z1))

      ||!mpfr_zero_p (mpc_imagref (z1)) || !mpfr_signbit (mpc_imagref (z1)))

    {

      printf ("Error: Regression, bug 20081021 reproduced\n");

      MPC_OUT (z1);

      exit (1);

    }

  mpc_clear (z1);

}

int

main (void)

{

  DECL_FUNC (CC, f, mpc_sqr);

  test_start ();

  testsqr (247, -65, 8, 24);

  testsqr (5, -896, 3, 2);

  testsqr (-3, -512, 2, 16);

  testsqr (266013312, 121990769, 27, 0);

  testsqr (170, 9, 8, 0);

  testsqr (768, 85, 8, 16);

  testsqr (145, 1816, 8, 24);

  testsqr (0, 1816, 8, 24);

  testsqr (145, 0, 8, 24);

  data_check (f, "sqr.dat");

  tgeneric (f, 2, 1024, 1, 0);

  reuse_bug ();

  test_end ();

  return 0;

}