1/* tmul -- test file for mpc_mul.
2
3Copyright (C) 2002, 2005, 2008, 2009, 2010, 2011, 2012 INRIA
4
5This file is part of GNU MPC.
6
7GNU MPC is free software; you can redistribute it and/or modify it under
8the terms of the GNU Lesser General Public License as published by the
9Free Software Foundation; either version 3 of the License, or (at your
10option) any later version.
11
12GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
13WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15more details.
16
17You should have received a copy of the GNU Lesser General Public License
18along 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
27static void
28cmpmul (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
65static void
66testmul (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
83static void
84check_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
135static void
136timemul (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
182int
183main (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}
202