tests/tsqr.c

   1 /* tsqr -- test file for mpc_sqr.
   2
   3 Copyright (C) 2002, 2005, 2008, 2010, 2011 INRIA
   4
   5 This file is part of GNU MPC.
   6
   7 GNU MPC is free software; you can redistribute it and/or modify it under
   8 the terms of the GNU Lesser General Public License as published by the
   9 Free Software Foundation; either version 3 of the License, or (at your
  10 option) any later version.
  11
  12 GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
  13 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  14 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
  15 more details.
  16
  17 You should have received a copy of the GNU Lesser General Public License
  18 along with this program. If not, see http://www.gnu.org/licenses/ .
  19 */
  20
  21 #include <stdlib.h>
  22 #include "mpc-tests.h"
  23
  24 static void
  25 cmpsqr (mpc_srcptr x, mpc_rnd_t rnd)
  26    /* computes the square of x with the specific function or by simple     */
  27    /* multiplication using the rounding mode rnd and compares the results  */
  28    /* and return values.                                                   */
  29    /* In our current test suite, the real and imaginary parts of x have    */
  30    /* the same precision, and we use this precision also for the result.   */
  31    /* Furthermore, we check whether computing the square in the same       */
  32    /* place yields the same result.                                        */
  33    /* We also compute the result with four times the precision and check   */
  34    /* whether the rounding is correct. Error reports in this part of the   */
  35    /* algorithm might still be wrong, though, since there are two          */
  36    /* consecutive roundings.                                               */
  37 {
  38   mpc_t z, t, u;
  39   int   inexact_z, inexact_t;
  40
  41   mpc_init2 (z, MPC_MAX_PREC (x));
  42   mpc_init2 (t, MPC_MAX_PREC (x));
  43   mpc_init2 (u, 4 * MPC_MAX_PREC (x));
  44
  45   inexact_z = mpc_sqr (z, x, rnd);
  46   inexact_t = mpc_mul (t, x, x, rnd);
  47
  48   if (mpc_cmp (z, t))
  49     {
  50       fprintf (stderr, "sqr and mul differ for rnd=(%s,%s) \nx=",
  51                mpfr_print_rnd_mode(MPC_RND_RE(rnd)),
  52                mpfr_print_rnd_mode(MPC_RND_IM(rnd)));
  53       mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);
  54       fprintf (stderr, "\nmpc_sqr gives ");
  55       mpc_out_str (stderr, 2, 0, z, MPC_RNDNN);
  56       fprintf (stderr, "\nmpc_mul gives ");
  57       mpc_out_str (stderr, 2, 0, t, MPC_RNDNN);
  58       fprintf (stderr, "\n");
  59       exit (1);
  60     }
  61   if (inexact_z != inexact_t)
  62     {
  63       fprintf (stderr, "The return values of sqr and mul differ for rnd=(%s,%s) \nx=  ",
  64                mpfr_print_rnd_mode(MPC_RND_RE(rnd)),
  65                mpfr_print_rnd_mode(MPC_RND_IM(rnd)));
  66       mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);
  67       fprintf (stderr, "\nx^2=");
  68       mpc_out_str (stderr, 2, 0, z, MPC_RNDNN);
  69       fprintf (stderr, "\nmpc_sqr gives %i", inexact_z);
  70       fprintf (stderr, "\nmpc_mul gives %i", inexact_t);
  71       fprintf (stderr, "\n");
  72       exit (1);
  73     }
  74
  75   mpc_set (t, x, MPC_RNDNN);
  76   inexact_t = mpc_sqr (t, t, rnd);
  77   if (mpc_cmp (z, t))
  78     {
  79       fprintf (stderr, "sqr and sqr in place differ for rnd=(%s,%s) \nx=",
  80                mpfr_print_rnd_mode(MPC_RND_RE(rnd)),
  81                mpfr_print_rnd_mode(MPC_RND_IM(rnd)));
  82       mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);
  83       fprintf (stderr, "\nmpc_sqr          gives ");
  84       mpc_out_str (stderr, 2, 0, z, MPC_RNDNN);
  85       fprintf (stderr, "\nmpc_sqr in place gives ");
  86       mpc_out_str (stderr, 2, 0, t, MPC_RNDNN);
  87       fprintf (stderr, "\n");
  88       exit (1);
  89     }
  90   if (inexact_z != inexact_t)
  91     {
  92       fprintf (stderr, "The return values of sqr and sqr in place differ for rnd=(%s,%s) \nx=  ",
  93                mpfr_print_rnd_mode(MPC_RND_RE(rnd)),
  94                mpfr_print_rnd_mode(MPC_RND_IM(rnd)));
  95       mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);
  96       fprintf (stderr, "\nx^2=");
  97       mpc_out_str (stderr, 2, 0, z, MPC_RNDNN);
  98       fprintf (stderr, "\nmpc_sqr          gives %i", inexact_z);
  99       fprintf (stderr, "\nmpc_sqr in place gives %i", inexact_t);
 100       fprintf (stderr, "\n");
 101       exit (1);
 102     }
 103
 104   mpc_sqr (u, x, rnd);
 105   mpc_set (t, u, rnd);
 106   if (mpc_cmp (z, t))
 107     {
 108       fprintf (stderr, "rounding in sqr might be incorrect for rnd=(%s,%s) \nx=",
 109                mpfr_print_rnd_mode(MPC_RND_RE(rnd)),
 110                mpfr_print_rnd_mode(MPC_RND_IM(rnd)));
 111       mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);
 112       fprintf (stderr, "\nmpc_sqr                     gives ");
 113       mpc_out_str (stderr, 2, 0, z, MPC_RNDNN);
 114       fprintf (stderr, "\nmpc_sqr quadruple precision gives ");
 115       mpc_out_str (stderr, 2, 0, u, MPC_RNDNN);
 116       fprintf (stderr, "\nand is rounded to                 ");
 117       mpc_out_str (stderr, 2, 0, t, MPC_RNDNN);
 118       fprintf (stderr, "\n");
 119       exit (1);
 120     }
 121
 122   mpc_clear (z);
 123   mpc_clear (t);
 124   mpc_clear (u);
 125 }
 126
 127
 128 static void
 129 testsqr (long a, long b, mpfr_prec_t prec, mpc_rnd_t rnd)
 130 {
 131   mpc_t x;
 132
 133   mpc_init2 (x, prec);
 134
 135   mpc_set_si_si (x, a, b, rnd);
 136
 137   cmpsqr (x, rnd);
 138
 139   mpc_clear (x);
 140 }
 141
 142
 143 static void
 144 reuse_bug (void)
 145 {
 146   mpc_t z1;
 147
 148   /* reuse bug found by Paul Zimmermann 20081021 */
 149   mpc_init2 (z1, 2);
 150   /* RE (z1^2) overflows, IM(z^2) = -0 */
 151   mpfr_set_str (mpc_realref (z1), "0.11", 2, GMP_RNDN);
 152   mpfr_mul_2si (mpc_realref (z1), mpc_realref (z1), mpfr_get_emax (), GMP_RNDN);
 153   mpfr_set_ui (mpc_imagref (z1), 0, GMP_RNDN);
 154   mpc_conj (z1, z1, MPC_RNDNN);
 155   mpc_sqr (z1, z1, MPC_RNDNN);
 156   if (!mpfr_inf_p (mpc_realref (z1)) || mpfr_signbit (mpc_realref (z1))
 157       ||!mpfr_zero_p (mpc_imagref (z1)) || !mpfr_signbit (mpc_imagref (z1)))
 158     {
 159       printf ("Error: Regression, bug 20081021 reproduced\n");
 160       MPC_OUT (z1);
 161       exit (1);
 162     }
 163
 164   mpc_clear (z1);
 165 }
 166
 167 int
 168 main (void)
 169 {
 170   DECL_FUNC (CC, f, mpc_sqr);
 171   test_start ();
 172
 173   testsqr (247, -65, 8, 24);
 174   testsqr (5, -896, 3, 2);
 175   testsqr (-3, -512, 2, 16);
 176   testsqr (266013312, 121990769, 27, 0);
 177   testsqr (170, 9, 8, 0);
 178   testsqr (768, 85, 8, 16);
 179   testsqr (145, 1816, 8, 24);
 180   testsqr (0, 1816, 8, 24);
 181   testsqr (145, 0, 8, 24);
 182
 183   data_check (f, "sqr.dat");
 184   tgeneric (f, 2, 1024, 1, 0);
 185
 186   reuse_bug ();
 187
 188   test_end ();
 189
 190   return 0;
 191 }