1/* tmul -- test file for mpc_mul. 2 3Copyright (C) 2002, 2005, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2020 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, 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 133static void 134bug20200206 (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 163static void 164timemul (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 221static void 222bug20221130 (void) 223{ 224 mpc_t b, c_conj, res, ref; 225 mpfr_prec_t prec; 226 mpc_init2 (b, 20); 227 mpc_init2 (c_conj, 20); 228 mpc_init2 (res, 2); 229 mpc_init2 (ref, 5); 230 mpc_set_str (b, "(0x1p+0 0x2p+0)", 16, MPC_RNDNN); 231 mpc_set_str (c_conj, "(-0xap+0 0x1.4p+4)", 16, MPC_RNDNN); 232 mpc_set_str (ref, "(-0x3.2p+4 0x0p+0)", 16, MPC_RNDNN); 233 for (prec = 5; prec <= 2000; prec++) 234 { 235 mpc_set_prec (res, prec); 236 mpc_mul (res, b, c_conj, MPC_RNDZZ); 237 if (mpc_cmp (res, ref) != 0 || mpfr_signbit (mpc_imagref (res))) 238 { 239 printf ("Error in bug20221130 for prec=%lu\n", prec); 240 mpfr_printf ("expected (%Ra %Ra)\n", mpc_realref (ref), mpc_imagref (ref)); 241 mpfr_printf ("got (%Ra %Ra)\n", mpc_realref (res), mpc_imagref (res)); 242 exit (1); 243 } 244 } 245 mpc_clear (b); 246 mpc_clear (c_conj); 247 mpc_clear (res); 248 mpc_clear (ref); 249} 250 251int 252main (void) 253{ 254 test_start (); 255 256#ifdef TIMING 257 timemul (); 258#endif 259 260 bug20200206 (); 261 bug20221130 (); 262 check_regular (); 263 264 data_check_template ("mul.dsc", "mul.dat"); 265 266 tgeneric_template ("mul.dsc", 2, 4096, 41, 1024); 267 268 test_end (); 269 270 return 0; 271} 272