1 /* tmul -- test file for mpc_mul.
2
3 Copyright (C) 2002, 2005, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2020 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
cmpmul(mpc_srcptr x,mpc_srcptr y,mpc_rnd_t rnd)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
testmul(long a,long b,long c,long d,mpfr_prec_t prec,mpc_rnd_t rnd)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
check_regular(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, MPFR_RNDN);
106 mpfr_set_str (mpc_imagref (x), "0x90p+244", 16, MPFR_RNDN);
107 mpfr_set_str (mpc_realref (y), "0xECp-146", 16, MPFR_RNDN);
108 mpfr_set_str (mpc_imagref (y), "0xACp-471", 16, MPFR_RNDN);
109 cmpmul (x, y, MPC_RNDNN);
110 mpfr_set_str (mpc_realref (x), "0xB4p+733", 16, MPFR_RNDN);
111 mpfr_set_str (mpc_imagref (x), "0x90p+244", 16, MPFR_RNDN);
112 mpfr_set_str (mpc_realref (y), "0xACp-471", 16, MPFR_RNDN);
113 mpfr_set_str (mpc_imagref (y), "-0xECp-146", 16, MPFR_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 static void
bug20200206(void)134 bug20200206 (void)
135 {
136 mpfr_exp_t emin = mpfr_get_emin ();
137 mpc_t x, y, z;
138
139 mpfr_set_emin (-1073);
140 mpc_init2 (x, 53);
141 mpc_init2 (y, 53);
142 mpc_init2 (z, 53);
143 mpfr_set_d (mpc_realref (x), -6.0344722345057644e-272, MPFR_RNDN);
144 mpfr_set_d (mpc_imagref (x), -4.8536770224196353e-204, MPFR_RNDN);
145 mpfr_set_d (mpc_realref (y), 1.3834775731431992e-246, MPFR_RNDN);
146 mpfr_set_d (mpc_imagref (y), 2.9246270396940562e-124, MPFR_RNDN);
147 mpc_mul (z, x, y, MPC_RNDNN);
148 if (mpfr_regular_p (mpc_realref (z)) &&
149 mpfr_get_exp (mpc_realref (z)) < -1073)
150 {
151 printf ("Error, mpc_mul returns an out-of-range exponent:\n");
152 mpfr_dump (mpc_realref (z));
153 printf ("Bug most probably in MPFR, please upgrade to MPFR 4.1.0 or later\n");
154 exit (1);
155 }
156 mpc_clear (x);
157 mpc_clear (y);
158 mpc_clear (z);
159 mpfr_set_emin (emin);
160 }
161
162 #ifdef TIMING
163 static void
timemul(void)164 timemul (void)
165 {
166 /* measures the time needed with different precisions for naive and */
167 /* Karatsuba multiplication */
168
169 mpc_t x, y, z;
170 unsigned long int i, j;
171 const unsigned long int tests = 10000;
172 struct tms time_old, time_new;
173 double passed1, passed2;
174
175 mpc_init (x);
176 mpc_init (y);
177 mpc_init_set_ui_ui (z, 1, 0, MPC_RNDNN);
178
179 for (i = 1; i < 50; i++)
180 {
181 mpc_set_prec (x, i * BITS_PER_MP_LIMB);
182 mpc_set_prec (y, i * BITS_PER_MP_LIMB);
183 mpc_set_prec (z, i * BITS_PER_MP_LIMB);
184 test_default_random (x, -1, 1, 128, 25);
185 test_default_random (y, -1, 1, 128, 25);
186
187 times (&time_old);
188 for (j = 0; j < tests; j++)
189 mpc_mul_naive (z, x, y, MPC_RNDNN);
190 times (&time_new);
191 passed1 = ((double) (time_new.tms_utime - time_old.tms_utime)) / 100;
192
193 times (&time_old);
194 for (j = 0; j < tests; j++)
195 mpc_mul_karatsuba (z, x, y, MPC_RNDNN);
196 times (&time_new);
197 passed2 = ((double) (time_new.tms_utime - time_old.tms_utime)) / 100;
198
199 printf ("Time for %3li limbs naive/Karatsuba: %5.2f %5.2f\n", i,
200 passed1, passed2);
201 }
202
203 mpc_clear (x);
204 mpc_clear (y);
205 mpc_clear (z);
206 }
207 #endif
208
209 #define MPC_FUNCTION_CALL \
210 P[0].mpc_inex = mpc_mul (P[1].mpc, P[2].mpc, P[3].mpc, P[4].mpc_rnd)
211 #define MPC_FUNCTION_CALL_SYMMETRIC \
212 P[0].mpc_inex = mpc_mul (P[1].mpc, P[3].mpc, P[2].mpc, P[4].mpc_rnd)
213 #define MPC_FUNCTION_CALL_REUSE_OP1 \
214 P[0].mpc_inex = mpc_mul (P[1].mpc, P[1].mpc, P[3].mpc, P[4].mpc_rnd)
215 #define MPC_FUNCTION_CALL_REUSE_OP2 \
216 P[0].mpc_inex = mpc_mul (P[1].mpc, P[2].mpc, P[1].mpc, P[4].mpc_rnd)
217
218 #include "data_check.tpl"
219 #include "tgeneric.tpl"
220
221 int
main(void)222 main (void)
223 {
224 test_start ();
225
226 #ifdef TIMING
227 timemul ();
228 #endif
229
230 bug20200206 ();
231 check_regular ();
232
233 data_check_template ("mul.dsc", "mul.dat");
234
235 tgeneric_template ("mul.dsc", 2, 4096, 41, 1024);
236
237 test_end ();
238
239 return 0;
240 }
241