From 8bb4f5a6a9b8ad0adf2d148d24d0d4e68f33a65b Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Sat, 10 Nov 2012 18:43:43 +0100 Subject: [PATCH] don't use the three-point complex multiplication by default, as it introduces too much error in some cases --- doc/source/fmpcb.rst | 10 ++ fmpcb.h | 2 + fmpcb/mul.c | 15 ++- fmpcb/mul_alt.c | 91 +++++++++++++++++ fmpcb/test/t-mul_alt.c | 221 +++++++++++++++++++++++++++++++++++++++++ 5 files changed, 331 insertions(+), 8 deletions(-) create mode 100644 fmpcb/mul_alt.c create mode 100644 fmpcb/test/t-mul_alt.c diff --git a/doc/source/fmpcb.rst b/doc/source/fmpcb.rst index d84998a2..e09f06d2 100644 --- a/doc/source/fmpcb.rst +++ b/doc/source/fmpcb.rst @@ -166,10 +166,20 @@ Arithmetic Sets *z* to the product of *x* and *y*. If at least one part of *x* or *y* is zero, the operations is reduced to two real multiplications. + +.. function:: void fmpcb_mul_alt(fmpcb_t z, const fmpcb_t x, const fmpcb_t y, long prec) + + Sets *z* to the product of *x* and *y*. If at least one part of + *x* or *y* is zero, the operations is reduced to two real multiplications. Otherwise, letting `x = a + bi`, `y = c + di`, `z = e + fi`, we use the formula `e = ac - bd`, `f = (a+b)(c+d) - ac - bd`, which requires three real multiplications instead of four. + The drawback of this algorithm is that the numerical stability is much + worse than for the default algorithm. In particular, if one operand + has a large error and the other a small error, the output error will + be about twice that of the large input error, rather than about the same. + .. function:: void fmpcb_addmul(fmpcb_t z, const fmpcb_t x, const fmpcb_t y, long prec) .. function:: void fmpcb_submul(fmpcb_t z, const fmpcb_t x, const fmpcb_t y, long prec) diff --git a/fmpcb.h b/fmpcb.h index 0b30b8a0..434e2156 100644 --- a/fmpcb.h +++ b/fmpcb.h @@ -259,6 +259,8 @@ fmpcb_mul_onei(fmpcb_t z, const fmpcb_t x) void fmpcb_mul(fmpcb_t z, const fmpcb_t x, const fmpcb_t y, long prec); +void fmpcb_mul_alt(fmpcb_t z, const fmpcb_t x, const fmpcb_t y, long prec); + static __inline__ void fmpcb_addmul(fmpcb_t z, const fmpcb_t x, const fmpcb_t y, long prec) { diff --git a/fmpcb/mul.c b/fmpcb/mul.c index 9ac6407f..d865f67b 100644 --- a/fmpcb/mul.c +++ b/fmpcb/mul.c @@ -59,26 +59,26 @@ fmpcb_mul(fmpcb_t z, const fmpcb_t x, const fmpcb_t y, long prec) } else { - fmprb_t t, u, v; + fmprb_t t, u, v, w; fmprb_init(t); fmprb_init(u); fmprb_init(v); - - fmprb_add(t, a, b, prec); - fmprb_add(u, c, d, prec); - fmprb_mul(v, t, u, prec); + fmprb_init(w); fmprb_mul(t, a, c, prec); fmprb_mul(u, b, d, prec); + fmprb_mul(v, a, d, prec); + fmprb_mul(w, b, c, prec); + fmprb_sub(e, t, u, prec); - fmprb_sub(f, v, t, prec); - fmprb_sub(f, f, u, prec); + fmprb_add(f, v, w, prec); fmprb_clear(t); fmprb_clear(u); fmprb_clear(v); + fmprb_clear(w); } #undef a @@ -88,4 +88,3 @@ fmpcb_mul(fmpcb_t z, const fmpcb_t x, const fmpcb_t y, long prec) #undef e #undef f } - diff --git a/fmpcb/mul_alt.c b/fmpcb/mul_alt.c new file mode 100644 index 00000000..cd39258b --- /dev/null +++ b/fmpcb/mul_alt.c @@ -0,0 +1,91 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB 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 General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2012 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb.h" + +void +fmpcb_mul_alt(fmpcb_t z, const fmpcb_t x, const fmpcb_t y, long prec) +{ +#define a fmpcb_realref(x) +#define b fmpcb_imagref(x) +#define c fmpcb_realref(y) +#define d fmpcb_imagref(y) +#define e fmpcb_realref(z) +#define f fmpcb_imagref(z) + + if (fmprb_is_zero(b)) + { + fmprb_mul(f, d, a, prec); + fmprb_mul(e, c, a, prec); + } + else if (fmprb_is_zero(d)) + { + fmprb_mul(f, b, c, prec); + fmprb_mul(e, a, c, prec); + } + else if (fmprb_is_zero(a)) + { + fmprb_mul(e, c, b, prec); + fmprb_mul(f, d, b, prec); + fmpcb_mul_onei(z, z); + } + else if (fmprb_is_zero(c)) + { + fmprb_mul(e, a, d, prec); + fmprb_mul(f, b, d, prec); + fmpcb_mul_onei(z, z); + } + else + { + fmprb_t t, u, v; + + fmprb_init(t); + fmprb_init(u); + fmprb_init(v); + + fmprb_add(t, a, b, prec); + fmprb_add(u, c, d, prec); + fmprb_mul(v, t, u, prec); + + fmprb_mul(t, a, c, prec); + fmprb_mul(u, b, d, prec); + + fmprb_sub(e, t, u, prec); + fmprb_sub(f, v, t, prec); + fmprb_sub(f, f, u, prec); + + fmprb_clear(t); + fmprb_clear(u); + fmprb_clear(v); + } + +#undef a +#undef b +#undef c +#undef d +#undef e +#undef f +} + diff --git a/fmpcb/test/t-mul_alt.c b/fmpcb/test/t-mul_alt.c new file mode 100644 index 00000000..57ed5747 --- /dev/null +++ b/fmpcb/test/t-mul_alt.c @@ -0,0 +1,221 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB 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 General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2012 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("mul_alt...."); + fflush(stdout); + + flint_randinit(state); + + /* test aliasing of c and a */ + for (iter = 0; iter < 100000; iter++) + { + fmpcb_t a, b, c; + long prec; + + fmpcb_init(a); + fmpcb_init(b); + fmpcb_init(c); + + fmpcb_randtest(a, state, 1 + n_randint(state, 200), 10); + fmpcb_randtest(b, state, 1 + n_randint(state, 200), 10); + fmpcb_randtest(c, state, 1 + n_randint(state, 200), 10); + + prec = 2 + n_randint(state, 200); + + fmpcb_mul_alt(c, a, b, prec); + fmpcb_mul_alt(a, a, b, prec); + + if (!fmpcb_equal(a, c)) + { + printf("FAIL: aliasing c, a\n\n"); + printf("a = "); fmpcb_print(a); printf("\n\n"); + printf("b = "); fmpcb_print(b); printf("\n\n"); + printf("c = "); fmpcb_print(c); printf("\n\n"); + abort(); + } + + fmpcb_clear(a); + fmpcb_clear(b); + fmpcb_clear(c); + } + + /* test aliasing of c and b */ + for (iter = 0; iter < 100000; iter++) + { + fmpcb_t a, b, c; + long prec; + + fmpcb_init(a); + fmpcb_init(b); + fmpcb_init(c); + + fmpcb_randtest(a, state, 1 + n_randint(state, 200), 10); + fmpcb_randtest(b, state, 1 + n_randint(state, 200), 10); + fmpcb_randtest(c, state, 1 + n_randint(state, 200), 10); + + prec = 2 + n_randint(state, 200); + + fmpcb_mul_alt(c, a, b, prec); + fmpcb_mul_alt(b, a, b, prec); + + if (!fmpcb_equal(b, c)) + { + printf("FAIL: aliasing b, a\n\n"); + printf("a = "); fmpcb_print(a); printf("\n\n"); + printf("b = "); fmpcb_print(b); printf("\n\n"); + printf("c = "); fmpcb_print(c); printf("\n\n"); + abort(); + } + + fmpcb_clear(a); + fmpcb_clear(b); + fmpcb_clear(c); + } + + /* test aliasing a, a */ + for (iter = 0; iter < 100000; iter++) + { + fmpcb_t a, b, c, d; + long prec; + + fmpcb_init(a); + fmpcb_init(b); + fmpcb_init(c); + fmpcb_init(d); + + fmpcb_randtest(a, state, 1 + n_randint(state, 200), 10); + fmpcb_randtest(b, state, 1 + n_randint(state, 200), 10); + fmpcb_randtest(c, state, 1 + n_randint(state, 200), 10); + + prec = 2 + n_randint(state, 200); + + fmpcb_set(b, a); + fmpcb_mul_alt(c, a, a, prec); + fmpcb_mul_alt(d, a, b, prec); + + if (!fmpcb_equal(c, d)) + { + printf("FAIL: aliasing a, a\n\n"); + printf("a = "); fmpcb_print(a); printf("\n\n"); + printf("b = "); fmpcb_print(b); printf("\n\n"); + printf("c = "); fmpcb_print(c); printf("\n\n"); + printf("d = "); fmpcb_print(d); printf("\n\n"); + abort(); + } + + fmpcb_clear(a); + fmpcb_clear(b); + fmpcb_clear(c); + fmpcb_clear(d); + } + + /* test aliasing a, a, a */ + for (iter = 0; iter < 100000; iter++) + { + fmpcb_t a, b, c; + long prec; + + fmpcb_init(a); + fmpcb_init(b); + fmpcb_init(c); + + fmpcb_randtest(a, state, 1 + n_randint(state, 200), 10); + fmpcb_randtest(b, state, 1 + n_randint(state, 200), 10); + fmpcb_randtest(c, state, 1 + n_randint(state, 200), 10); + + prec = 2 + n_randint(state, 200); + + fmpcb_set(b, a); + fmpcb_mul_alt(c, a, b, prec); + fmpcb_mul_alt(a, a, a, prec); + + if (!fmpcb_equal(a, c)) + { + printf("FAIL: aliasing a, a, a\n\n"); + printf("a = "); fmpcb_print(a); printf("\n\n"); + printf("b = "); fmpcb_print(b); printf("\n\n"); + printf("c = "); fmpcb_print(c); printf("\n\n"); + abort(); + } + + fmpcb_clear(a); + fmpcb_clear(b); + fmpcb_clear(c); + } + + /* test a*(b+c) = a*b + a*c */ + for (iter = 0; iter < 100000; iter++) + { + fmpcb_t a, b, c, d, e, f; + + fmpcb_init(a); + fmpcb_init(b); + fmpcb_init(c); + fmpcb_init(d); + fmpcb_init(e); + fmpcb_init(f); + + fmpcb_randtest(a, state, 1 + n_randint(state, 200), 10); + fmpcb_randtest(b, state, 1 + n_randint(state, 200), 10); + fmpcb_randtest(c, state, 1 + n_randint(state, 200), 10); + + fmpcb_add(d, b, c, 2 + n_randint(state, 200)); + fmpcb_mul_alt(e, a, d, 2 + n_randint(state, 200)); + + fmpcb_mul_alt(d, a, b, 2 + n_randint(state, 200)); + fmpcb_mul_alt(f, a, c, 2 + n_randint(state, 200)); + fmpcb_add(f, d, f, 2 + n_randint(state, 200)); + + if (!fmpcb_overlaps(e, f)) + { + printf("FAIL: a*(b+c) = a*b + a*c\n\n"); + printf("a = "); fmpcb_print(a); printf("\n\n"); + printf("b = "); fmpcb_print(b); printf("\n\n"); + printf("c = "); fmpcb_print(c); printf("\n\n"); + printf("e = "); fmpcb_print(e); printf("\n\n"); + printf("f = "); fmpcb_print(f); printf("\n\n"); + abort(); + } + + fmpcb_clear(a); + fmpcb_clear(b); + fmpcb_clear(c); + fmpcb_clear(d); + fmpcb_clear(e); + fmpcb_clear(f); + } + + flint_randclear(state); + _fmpz_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +}