From 615ea353a5dc933e973fb82d76c58159cff88a15 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Wed, 7 Nov 2012 11:48:22 +0100 Subject: [PATCH] some more work on the fmpcb module; test code for multiplication --- doc/source/fmpz_holonomic.rst | 2 +- doc/source/history.rst | 21 ++++ fmpcb.h | 80 +++--------- fmpcb/mul.c | 91 ++++++++++++++ fmpcb/randtest.c | 34 ++++++ fmpcb/test/t-mul.c | 221 ++++++++++++++++++++++++++++++++++ 6 files changed, 388 insertions(+), 61 deletions(-) create mode 100644 fmpcb/mul.c create mode 100644 fmpcb/randtest.c create mode 100644 fmpcb/test/t-mul.c diff --git a/doc/source/fmpz_holonomic.rst b/doc/source/fmpz_holonomic.rst index 9370f210..baa28cf8 100644 --- a/doc/source/fmpz_holonomic.rst +++ b/doc/source/fmpz_holonomic.rst @@ -325,7 +325,7 @@ Sequence evaluation .. function:: void fmpz_holonomic_get_nth_fmpq(fmpq_t res, const fmpz_holonomic_t op, const fmpq * initial, long n0, long n) -.. function:: mp_limb_t fmpz_holonomic_get_nth_nmod(const fmpz_holonomic_t op, mp_srcptr initial, long n0, long n, nmod_t mod) +.. function:: mp_limb_t fmpz_holonomic_get_nth_nmod(const fmpz_holonomic_t op, mp_srcptr initial, ulong n0, ulong n, nmod_t mod) Computes element `c(n)` in the sequence annihilated by the difference operator *op*, given the diff --git a/doc/source/history.rst b/doc/source/history.rst index f02aa1a8..5e5d74a5 100644 --- a/doc/source/history.rst +++ b/doc/source/history.rst @@ -4,6 +4,27 @@ History and changes For more details, view the detailed commit log in the git repository https://github.com/fredrik-johansson/arb +* version 0.3-git + + * converted documentation to sphinx + * initial support for complex numbers and polynomials + + * root isolation for complex polynomials + + * new module fmpz_holonomic for functions/sequences + defined by linear differential/difference equations + with polynomial coefficients + + * functions for creating various special sequences and functions + * some closure properties for sequences + * Taylor series expansion for differential equations + * computing the nth entry of a sequence using binary splitting + * computing the nth entry mod p using fast multipoint evaluation + + * generic code for evaluating hypergeometric series + * matrix powering + * various other helper functions + * 2012-09-29 - version 0.2 * code for computing the gamma function (Karatsuba, Stirling's series) diff --git a/fmpcb.h b/fmpcb.h index b645d903..f7cf0b8c 100644 --- a/fmpcb.h +++ b/fmpcb.h @@ -125,6 +125,13 @@ fmpcb_swap(fmpcb_t z, fmpcb_t x) fmprb_swap(fmpcb_imagref(z), fmpcb_imagref(x)); } +static __inline__ int +fmpcb_equal(const fmpcb_t x, const fmpcb_t y) +{ + return fmprb_equal(fmpcb_realref(x), fmpcb_realref(y)) && + fmprb_equal(fmpcb_imagref(x), fmpcb_imagref(y)); +} + static __inline__ int fmpcb_overlaps(const fmpcb_t x, const fmpcb_t y) { @@ -199,66 +206,7 @@ fmpcb_mul_onei(fmpcb_t z, const fmpcb_t x) } } - -static __inline__ void -fmpcb_mul(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) - -/* - TODO: use these optimizations, but fix aliasing issues - - if (fmprb_is_zero(b)) - { - fmpcb_mul_fmprb(z, y, a, prec); - } - else if (fmprb_is_zero(d)) - { - fmpcb_mul_fmprb(z, x, c, prec); - } - else if (fmprb_is_zero(a)) - { - fmpcb_mul_fmprb(z, y, b, prec); - fmpcb_mul_i(z, z); - } - else if (fmprb_is_zero(c)) - { - fmpcb_mul_fmprb(z, x, d, prec); - fmpcb_mul_i(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(fmpcb_realref(z), t, u, prec); - fmprb_sub(fmpcb_imagref(z), v, t, prec); - fmprb_sub(fmpcb_imagref(z), fmpcb_imagref(z), u, prec); - - fmprb_clear(t); - fmprb_clear(u); - fmprb_clear(v); - } - -#undef a -#undef b -#undef c -#undef d -} +void fmpcb_mul(fmpcb_t z, const fmpcb_t x, const fmpcb_t y, long prec); static __inline__ void fmpcb_inv(fmpcb_t z, const fmpcb_t x, long prec) @@ -300,6 +248,18 @@ _fmpcb_vec_set(fmpcb_struct * res, const fmpcb_struct * vec, long len) fmpcb_set(res + i, vec + i); } +static __inline__ void +fmpcb_print(const fmpcb_t x) +{ + printf("("); + fmprb_print(fmpcb_realref(x)); + printf(", "); + fmprb_print(fmpcb_imagref(x)); + printf(")"); +} + void fmpcb_printd(const fmpcb_t z, long digits); +void fmpcb_randtest(fmpcb_t z, flint_rand_t state, long prec, long mag_bits); + #endif diff --git a/fmpcb/mul.c b/fmpcb/mul.c new file mode 100644 index 00000000..9ac6407f --- /dev/null +++ b/fmpcb/mul.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(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/randtest.c b/fmpcb/randtest.c new file mode 100644 index 00000000..2a6bf1f2 --- /dev/null +++ b/fmpcb/randtest.c @@ -0,0 +1,34 @@ +/*============================================================================= + + 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_randtest(fmpcb_t z, flint_rand_t state, long prec, long mag_bits) +{ + fmprb_randtest(fmpcb_realref(z), state, prec, mag_bits); + fmprb_randtest(fmpcb_imagref(z), state, prec, mag_bits); +} + diff --git a/fmpcb/test/t-mul.c b/fmpcb/test/t-mul.c new file mode 100644 index 00000000..ec59f6d2 --- /dev/null +++ b/fmpcb/test/t-mul.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...."); + 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(c, a, b, prec); + fmpcb_mul(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(c, a, b, prec); + fmpcb_mul(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(c, a, a, prec); + fmpcb_mul(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(c, a, b, prec); + fmpcb_mul(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(e, a, d, 2 + n_randint(state, 200)); + + fmpcb_mul(d, a, b, 2 + n_randint(state, 200)); + fmpcb_mul(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; +}