tests/tmul.c

   1 /* tmul -- test file for mpc_mul.
   2
   3 Copyright (C) 2002, 2005, 2008, 2009, 2010, 2011, 2012 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 #ifdef TIMING
  23 #include <sys/times.h>
  24 #endif
  25 #include "mpc-tests.h"
  26
  27 static void
  28 cmpmul (mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
  29    /* computes the product of x and y with the naive and Karatsuba methods */
  30    /* using the rounding mode rnd and compares the results and return      */
  31    /* values.                                                              */
  32    /* In our current test suite, the real and imaginary parts of x and y   */
  33    /* all have the same precision, and we use this precision also for the  */
  34    /* result.                                                              */
  35 {
  36    mpc_t z, t;
  37    int   inex_z, inex_t;
  38
  39    mpc_init2 (z, MPC_MAX_PREC (x));
  40    mpc_init2 (t, MPC_MAX_PREC (x));
  41
  42    inex_z = mpc_mul_naive (z, x, y, rnd);
  43    inex_t = mpc_mul_karatsuba (t, x, y, rnd);
  44
  45    if (mpc_cmp (z, t) != 0 || inex_z != inex_t) {
  46       fprintf (stderr, "mul_naive and mul_karatsuba differ for rnd=(%s,%s)\n",
  47                mpfr_print_rnd_mode(MPC_RND_RE(rnd)),
  48                mpfr_print_rnd_mode(MPC_RND_IM(rnd)));
  49       MPC_OUT (x);
  50       MPC_OUT (y);
  51       MPC_OUT (z);
  52       MPC_OUT (t);
  53       if (inex_z != inex_t) {
  54          fprintf (stderr, "inex_re (z): %s\n", MPC_INEX_STR (inex_z));
  55          fprintf (stderr, "inex_re (t): %s\n", MPC_INEX_STR (inex_t));
  56       }
  57       exit (1);
  58    }
  59
  60    mpc_clear (z);
  61    mpc_clear (t);
  62 }
  63
  64
  65 static void
  66 testmul (long a, long b, long c, long d, mpfr_prec_t prec, mpc_rnd_t rnd)
  67 {
  68   mpc_t x, y;
  69
  70   mpc_init2 (x, prec);
  71   mpc_init2 (y, prec);
  72
  73   mpc_set_si_si (x, a, b, rnd);
  74   mpc_set_si_si (y, c, d, rnd);
  75
  76   cmpmul (x, y, rnd);
  77
  78   mpc_clear (x);
  79   mpc_clear (y);
  80 }
  81
  82
  83 static void
  84 check_regular (void)
  85 {
  86   mpc_t x, y;
  87   int rnd_re, rnd_im;
  88   mpfr_prec_t prec;
  89
  90   testmul (247, -65, -223, 416, 8, 24);
  91   testmul (5, -896, 5, -32, 3, 2);
  92   testmul (-3, -512, -1, -1, 2, 16);
  93   testmul (266013312, 121990769, 110585572, 116491059, 27, 0);
  94   testmul (170, 9, 450, 251, 8, 0);
  95   testmul (768, 85, 169, 440, 8, 16);
  96   testmul (145, 1816, 848, 169, 8, 24);
  97
  98   mpc_init2 (x, 1000);
  99   mpc_init2 (y, 1000);
 100
 101   /* Bug 20081114: mpc_mul_karatsuba returned wrong inexact value for
 102      imaginary part */
 103   mpc_set_prec (x, 7);
 104   mpc_set_prec (y, 7);
 105   mpfr_set_str (mpc_realref (x), "0xB4p+733", 16, GMP_RNDN);
 106   mpfr_set_str (mpc_imagref (x), "0x90p+244", 16, GMP_RNDN);
 107   mpfr_set_str (mpc_realref (y), "0xECp-146", 16, GMP_RNDN);
 108   mpfr_set_str (mpc_imagref (y), "0xACp-471", 16, GMP_RNDN);
 109   cmpmul (x, y, MPC_RNDNN);
 110   mpfr_set_str (mpc_realref (x), "0xB4p+733", 16, GMP_RNDN);
 111   mpfr_set_str (mpc_imagref (x), "0x90p+244", 16, GMP_RNDN);
 112   mpfr_set_str (mpc_realref (y), "0xACp-471", 16, GMP_RNDN);
 113   mpfr_set_str (mpc_imagref (y), "-0xECp-146", 16, GMP_RNDN);
 114   cmpmul (x, y, MPC_RNDNN);
 115
 116   for (prec = 2; prec < 1000; prec = (mpfr_prec_t) (prec * 1.1 + 1))
 117     {
 118       mpc_set_prec (x, prec);
 119       mpc_set_prec (y, prec);
 120
 121       test_default_random (x, -1024, 1024, 128, 0);
 122       test_default_random (y, -1024, 1024, 128, 0);
 123
 124       for (rnd_re = 0; rnd_re < 4; rnd_re ++)
 125         for (rnd_im = 0; rnd_im < 4; rnd_im ++)
 126           cmpmul (x, y, MPC_RND (rnd_re, rnd_im));
 127     }
 128
 129   mpc_clear (x);
 130   mpc_clear (y);
 131 }
 132
 133
 134 #ifdef TIMING
 135 static void
 136 timemul (void)
 137 {
 138   /* measures the time needed with different precisions for naive and */
 139   /* Karatsuba multiplication                                         */
 140
 141   mpc_t             x, y, z;
 142   unsigned long int i, j;
 143   const unsigned long int tests = 10000;
 144   struct tms        time_old, time_new;
 145   double            passed1, passed2;
 146
 147   mpc_init (x);
 148   mpc_init (y);
 149   mpc_init_set_ui_ui (z, 1, 0, MPC_RNDNN);
 150
 151   for (i = 1; i < 50; i++)
 152     {
 153       mpc_set_prec (x, i * BITS_PER_MP_LIMB);
 154       mpc_set_prec (y, i * BITS_PER_MP_LIMB);
 155       mpc_set_prec (z, i * BITS_PER_MP_LIMB);
 156       test_default_random (x, -1, 1, 128, 25);
 157       test_default_random (y, -1, 1, 128, 25);
 158
 159       times (&time_old);
 160       for (j = 0; j < tests; j++)
 161         mpc_mul_naive (z, x, y, MPC_RNDNN);
 162       times (&time_new);
 163       passed1 = ((double) (time_new.tms_utime - time_old.tms_utime)) / 100;
 164
 165       times (&time_old);
 166       for (j = 0; j < tests; j++)
 167         mpc_mul_karatsuba (z, x, y, MPC_RNDNN);
 168       times (&time_new);
 169       passed2 = ((double) (time_new.tms_utime - time_old.tms_utime)) / 100;
 170
 171       printf ("Time for %3li limbs naive/Karatsuba: %5.2f %5.2f\n", i,
 172               passed1, passed2);
 173     }
 174
 175   mpc_clear (x);
 176   mpc_clear (y);
 177   mpc_clear (z);
 178 }
 179 #endif
 180
 181
 182 int
 183 main (void)
 184 {
 185   DECL_FUNC (C_CC, f, mpc_mul);
 186   f.properties = FUNC_PROP_SYMETRIC;
 187
 188   test_start ();
 189
 190 #ifdef TIMING
 191   timemul ();
 192 #endif
 193
 194   check_regular ();
 195
 196   data_check (f, "mul.dat");
 197   tgeneric (f, 2, 4096, 41, 100);
 198
 199   test_end ();
 200   return 0;
 201 }