mirror of
https://github.com/vale981/arb
synced 2025-03-05 09:21:38 -05:00
don't use the three-point complex multiplication by default, as it introduces too much error in some cases
This commit is contained in:
Fredrik Johansson
parent
45ad953f69
commit
8bb4f5a6a9
5 changed files with 331 additions and 8 deletions
|
|
@ -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)
|
||||
|
|
|
|||
2
fmpcb.h
2
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)
|
||||
{
|
||||
|
|
|
|||
15
fmpcb/mul.c
15
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
|
||||
}
|
||||
|
||||
|
|
|
|||
91
fmpcb/mul_alt.c
Normal file
91
fmpcb/mul_alt.c
Normal file
|
|
@ -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
|
||||
}
|
||||
|
||||
221
fmpcb/test/t-mul_alt.c
Normal file
221
fmpcb/test/t-mul_alt.c
Normal file
|
|
@ -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;
|
||||
}
|
||||
Loading…
Add table
Reference in a new issue