summaryrefslogtreecommitdiff
path: root/tmul.c
diff options
context:
space:
mode:
Diffstat (limited to 'tmul.c')
-rw-r--r--tmul.c334
1 files changed, 334 insertions, 0 deletions
diff --git a/tmul.c b/tmul.c
new file mode 100644
index 0000000..4970779
--- /dev/null
+++ b/tmul.c
@@ -0,0 +1,334 @@
+/* tmul -- test file for mpc_mul.
+
+Copyright (C) 2002 Andreas Enge, Paul Zimmermann
+
+This file is part of the MPC Library.
+
+The MPC Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The MPC Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the MPC Library; see the file COPYING.LIB. If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <sys/times.h>
+#include "gmp.h"
+#include "mpfr.h"
+#include "mpc.h"
+#include "mpc-impl.h"
+
+void cmpmul _PROTO((mpc_srcptr, mpc_srcptr, mp_rnd_t));
+void testmul _PROTO((long, long, long, long, mp_prec_t, mp_rnd_t));
+void special _PROTO((void));
+void timemul _PROTO((void));
+
+
+void cmpmul (mpc_srcptr x, mpc_srcptr y, mp_rnd_t rnd)
+ /* computes the product of x and y with the naive and Karatsuba methods */
+ /* using the rounding mode rnd and compares the results and return */
+ /* values. */
+ /* In our current test suite, the real and imaginary parts of x and y */
+ /* all have the same precision, and we use this precision also for the */
+ /* result. */
+ /* Furthermore, we check whether the multiplication with one of the */
+ /* input arguments being also the output variable yields the same */
+ /* result. */
+ /* We also compute the result with four times the precision and check */
+ /* whether the rounding is correct. Error reports in this part of the */
+ /* algorithm might still be wrong, though, since there are two */
+ /* consecutive roundings. */
+{
+ mpc_t z, t, u;
+ int inexact_z, inexact_t;
+
+ mpc_init2 (z, MPC_MAX_PREC (x));
+ mpc_init2 (t, MPC_MAX_PREC (x));
+ mpc_init2 (u, 4 * MPC_MAX_PREC (x));
+
+ inexact_z = mpc_mul_naive (z, x, y, rnd);
+ inexact_t = mpc_mul_karatsuba (t, x, y, rnd);
+
+ if (mpc_cmp (z, t))
+ {
+ fprintf (stderr, "mul and mul2 differ for rnd=(%s,%s) \nx=",
+ mpfr_print_rnd_mode(MPC_RND_RE(rnd)),
+ mpfr_print_rnd_mode(MPC_RND_IM(rnd)));
+ mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);
+ fprintf (stderr, "\nand y=");
+ mpc_out_str (stderr, 2, 0, y, MPC_RNDNN);
+ fprintf (stderr, "\nmpc_mul_naive gives ");
+ mpc_out_str (stderr, 2, 0, z, MPC_RNDNN);
+ fprintf (stderr, "\nmpc_mul_karatsuba gives ");
+ mpc_out_str (stderr, 2, 0, t, MPC_RNDNN);
+ fprintf (stderr, "\n");
+ exit (1);
+ }
+ if (inexact_z != inexact_t)
+ {
+ fprintf (stderr, "The return values of mul and mul2 differ for rnd=(%s,%s) \nx=",
+ mpfr_print_rnd_mode(MPC_RND_RE(rnd)),
+ mpfr_print_rnd_mode(MPC_RND_IM(rnd)));
+ mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);
+ fprintf (stderr, "\nand y=");
+ mpc_out_str (stderr, 2, 0, y, MPC_RNDNN);
+ fprintf (stderr, "\nand x*y=");
+ mpc_out_str (stderr, 2, 0, z, MPC_RNDNN);
+ fprintf (stderr, "\nmpc_mul_naive gives %i", inexact_z);
+ fprintf (stderr, "\nmpc_mul_karatsuba gives %i", inexact_t);
+ fprintf (stderr, "\n");
+ exit (1);
+ }
+
+ mpc_set (t, x, MPC_RNDNN);
+ inexact_t = mpc_mul (t, t, y, rnd);
+ if (mpc_cmp (z, t))
+ {
+ fprintf (stderr, "mul and mul with the first variable in place differ for rnd=(%s,%s) \nx=",
+ mpfr_print_rnd_mode(MPC_RND_RE(rnd)),
+ mpfr_print_rnd_mode(MPC_RND_IM(rnd)));
+ mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);
+ fprintf (stderr, "\nand y=");
+ mpc_out_str (stderr, 2, 0, y, MPC_RNDNN);
+ fprintf (stderr, "\nmpc_mul gives ");
+ mpc_out_str (stderr, 2, 0, z, MPC_RNDNN);
+ fprintf (stderr, "\nmpc_mul in place gives ");
+ mpc_out_str (stderr, 2, 0, t, MPC_RNDNN);
+ fprintf (stderr, "\n");
+ exit (1);
+ }
+ if (inexact_z != inexact_t)
+ {
+ fprintf (stderr, "The return values of mul and mul with the first variable in place differ for rnd=(%s,%s) \nx=",
+ mpfr_print_rnd_mode(MPC_RND_RE(rnd)),
+ mpfr_print_rnd_mode(MPC_RND_IM(rnd)));
+ mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);
+ fprintf (stderr, "\nand y=");
+ mpc_out_str (stderr, 2, 0, y, MPC_RNDNN);
+ fprintf (stderr, "\nand x*y=");
+ mpc_out_str (stderr, 2, 0, z, MPC_RNDNN);
+ fprintf (stderr, "\nmpc_mul gives %i", inexact_z);
+ fprintf (stderr, "\nmpc_mul in place gives %i", inexact_t);
+ fprintf (stderr, "\n");
+ exit (1);
+ }
+
+ mpc_set (t, y, MPC_RNDNN);
+ inexact_t = mpc_mul (t, x, t, rnd);
+ if (mpc_cmp (z, t))
+ {
+ fprintf (stderr, "mul and mul with the second variable in place differ for rnd=(%s,%s) \nx=",
+ mpfr_print_rnd_mode(MPC_RND_RE(rnd)),
+ mpfr_print_rnd_mode(MPC_RND_IM(rnd)));
+ mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);
+ fprintf (stderr, "\nand y=");
+ mpc_out_str (stderr, 2, 0, y, MPC_RNDNN);
+ fprintf (stderr, "\nmpc_mul gives ");
+ mpc_out_str (stderr, 2, 0, z, MPC_RNDNN);
+ fprintf (stderr, "\nmpc_mul in place gives ");
+ mpc_out_str (stderr, 2, 0, t, MPC_RNDNN);
+ fprintf (stderr, "\n");
+ exit (1);
+ }
+ if (inexact_z != inexact_t)
+ {
+ fprintf (stderr, "The return values of mul and mul with the second variable in place differ for rnd=(%s,%s) \nx=",
+ mpfr_print_rnd_mode(MPC_RND_RE(rnd)),
+ mpfr_print_rnd_mode(MPC_RND_IM(rnd)));
+ mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);
+ fprintf (stderr, "\nand y=");
+ mpc_out_str (stderr, 2, 0, y, MPC_RNDNN);
+ fprintf (stderr, "\nand x*y=");
+ mpc_out_str (stderr, 2, 0, z, MPC_RNDNN);
+ fprintf (stderr, "\nmpc_mul gives %i", inexact_z);
+ fprintf (stderr, "\nmpc_mul in place gives %i", inexact_t);
+ fprintf (stderr, "\n");
+ exit (1);
+ }
+
+ mpc_mul (u, x, y, rnd);
+ mpc_set (t, u, rnd);
+ if (mpc_cmp (z, t))
+ {
+ fprintf (stderr, "rounding in mul might be incorrect for rnd=(%s,%s) \nx=",
+ mpfr_print_rnd_mode(MPC_RND_RE(rnd)),
+ mpfr_print_rnd_mode(MPC_RND_IM(rnd)));
+ mpc_out_str (stderr, 2, 0, x, MPC_RNDNN);
+ fprintf (stderr, "\nand y=");
+ mpc_out_str (stderr, 2, 0, y, MPC_RNDNN);
+ fprintf (stderr, "\nmpc_mul gives ");
+ mpc_out_str (stderr, 2, 0, z, MPC_RNDNN);
+ fprintf (stderr, "\nmpc_mul quadruple precision gives ");
+ mpc_out_str (stderr, 2, 0, u, MPC_RNDNN);
+ fprintf (stderr, "\nand is rounded to ");
+ mpc_out_str (stderr, 2, 0, t, MPC_RNDNN);
+ fprintf (stderr, "\n");
+ exit (1);
+ }
+
+ mpc_clear (z);
+ mpc_clear (t);
+ mpc_clear (u);
+}
+
+
+void
+testmul (long a, long b, long c, long d, mp_prec_t prec, mp_rnd_t rnd)
+{
+ mpc_t x, y;
+
+ mpc_init2 (x, prec);
+ mpc_init2 (y, prec);
+
+ mpc_set_si_si (x, a, b, rnd);
+ mpc_set_si_si (y, c, d, rnd);
+
+ cmpmul (x, y, rnd);
+
+ mpc_clear (x);
+ mpc_clear (y);
+}
+
+
+void
+special ()
+{
+ mpc_t x, y, z, t;
+ int inexact;
+
+ mpc_init (x);
+ mpc_init (y);
+ mpc_init (z);
+ mpc_init (t);
+
+ mpc_set_prec (x, 8);
+ mpc_set_prec (y, 8);
+ mpc_set_prec (z, 8);
+ mpc_set_si_si (x, 4, 3, MPC_RNDNN);
+ mpc_set_si_si (y, 1, -2, MPC_RNDNN);
+ inexact = mpc_mul (z, x, y, MPC_RNDNN);
+ if (MPC_INEX_RE(inexact) || MPC_INEX_IM(inexact))
+ {
+ fprintf (stderr, "Error: (4+3*I)*(1-2*I) should be exact with prec=8\n");
+ exit (1);
+ }
+
+ mpc_set_prec (x, 27);
+ mpc_set_prec (y, 27);
+ mpfr_set_str_raw (MPC_RE(x), "1.11111011011000010101000000e-2");
+ mpfr_set_str_raw (MPC_IM(x), "1.11010001010110111001110001e-3");
+ mpfr_set_str_raw (MPC_RE(y), "1.10100101110110011011100100e-1");
+ mpfr_set_str_raw (MPC_IM(y), "1.10111100011000001100110011e-1");
+ cmpmul (x, y, 0);
+
+ mpc_clear (x);
+ mpc_clear (y);
+}
+
+
+void
+timemul ()
+{
+ /* measures the time needed with different precisions for naive and */
+ /* Karatsuba multiplication */
+
+ mpc_t x, y, z;
+ unsigned long int i, j;
+ const unsigned long int tests = 10000;
+ struct tms time_old, time_new;
+ double passed1, passed2;
+
+ mpc_init (x);
+ mpc_init (y);
+ mpc_init_set_ui_ui (z, 1, 0, MPC_RNDNN);
+
+ for (i = 1; i < 50; i++)
+ {
+ mpc_set_prec (x, i * BITS_PER_MP_LIMB);
+ mpc_set_prec (y, i * BITS_PER_MP_LIMB);
+ mpc_set_prec (z, i * BITS_PER_MP_LIMB);
+ mpc_random (x);
+ mpc_random (y);
+
+ times (&time_old);
+ for (j = 0; j < tests; j++)
+ mpc_mul_naive (z, x, y, MPC_RNDNN);
+ times (&time_new);
+ passed1 = ((double) (time_new.tms_utime - time_old.tms_utime)) / 100;
+
+ times (&time_old);
+ for (j = 0; j < tests; j++)
+ mpc_mul_karatsuba (z, x, y, MPC_RNDNN);
+ times (&time_new);
+ passed2 = ((double) (time_new.tms_utime - time_old.tms_utime)) / 100;
+
+ printf ("Time for %3li limbs naive/Karatsuba: %5.2f %5.2f\n", i,
+ passed1, passed2);
+ }
+
+ mpc_clear (x);
+ mpc_clear (y);
+ mpc_clear (z);
+}
+
+
+int
+main()
+{
+ mpc_t x, y;
+ mp_rnd_t rnd_re, rnd_im;
+ mp_prec_t prec;
+ int i;
+
+/*
+ timemul ();
+*/
+
+ special ();
+
+ testmul (247, -65, -223, 416, 8, 24);
+ testmul (5, -896, 5, -32, 3, 2);
+ testmul (-3, -512, -1, -1, 2, 16);
+ testmul (266013312, 121990769, 110585572, 116491059, 27, 0);
+ testmul (170, 9, 450, 251, 8, 0);
+ testmul (768, 85, 169, 440, 8, 16);
+ testmul (145, 1816, 848, 169, 8, 24);
+ testmul (0, 1816, 848, 169, 8, 24);
+ testmul (145, 0, 848, 169, 8, 24);
+ testmul (145, 1816, 0, 169, 8, 24);
+ testmul (145, 1816, 848, 0, 8, 24);
+
+ mpc_init (x);
+ mpc_init (y);
+
+ for (prec = 2; prec < 1000; prec++) {
+
+ mpc_set_prec (x, prec);
+ mpc_set_prec (y, prec);
+
+ for (i = 0; i < 1000/prec; i++)
+ {
+
+ mpc_random (x);
+ mpc_random (y);
+
+ for (rnd_re = 0; rnd_re < 4; rnd_re ++)
+ for (rnd_im = 0; rnd_im < 4; rnd_im ++)
+ cmpmul (x, y, RNDC(rnd_re, rnd_im));
+ }
+ }
+
+ mpc_clear (x);
+ mpc_clear (y);
+
+ return 0;
+}
View Lorry Controller Status