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