Extending the mpn interface
Niels Möller
nisse at lysator.liu.se
Sat Feb 23 17:45:43 CET 2013
nisse at lysator.liu.se (Niels Möller) writes:
> Another revision, this time as a patch. As you can see, I renamed things
> a bit more. Compile tested only...
New patch below, including some testcases (no tests for mpz_roinit_n and
MPZ_ROINIT_N though). I'd like to check this in now. Any objections?
While writing the tests, I discovered some additional things I miss. See
tests/mpz/t-limbs.c, the alt_mul function. IT would be useful to have a
function to allocate a new limb array (using gmp:s alloc function), and
use that later to replace the limb array in an mpz.
Useful for inplace operations, where the underlying mpn function can't
work in-place. Then it's desirable to pass the old limb array as input,
and the newly allocated array as output, and replace when done.
Regards,
/Niels
diff -r 32f62aca5666 Makefile.am
--- a/Makefile.am Thu Feb 21 11:03:04 2013 +0100
+++ b/Makefile.am Sat Feb 23 17:38:59 2013 +0100
@@ -174,7 +174,9 @@ MPZ_OBJECTS = mpz/abs$U.lo mpz/add$U.lo
mpz/ior$U.lo mpz/iset$U.lo mpz/iset_d$U.lo mpz/iset_si$U.lo \
mpz/iset_str$U.lo mpz/iset_ui$U.lo mpz/jacobi$U.lo mpz/kronsz$U.lo \
mpz/kronuz$U.lo mpz/kronzs$U.lo mpz/kronzu$U.lo \
- mpz/lcm$U.lo mpz/lcm_ui$U.lo mpz/lucnum_ui$U.lo mpz/lucnum2_ui$U.lo \
+ mpz/lcm$U.lo mpz/lcm_ui$U.lo mpz/limbs_finish$U.lo \
+ mpz/limbs_modify$U.lo mpz/limbs_read$U.lo mpz/limbs_write.lo \
+ mpz/lucnum_ui$U.lo mpz/lucnum2_ui$U.lo \
mpz/millerrabin$U.lo mpz/mod$U.lo mpz/mul$U.lo mpz/mul_2exp$U.lo \
mpz/mul_si$U.lo mpz/mul_ui$U.lo \
mpz/n_pow_ui$U.lo mpz/neg$U.lo mpz/nextprime$U.lo \
@@ -182,7 +184,7 @@ MPZ_OBJECTS = mpz/abs$U.lo mpz/add$U.lo
mpz/popcount$U.lo mpz/pow_ui$U.lo mpz/powm$U.lo mpz/powm_sec$U.lo \
mpz/powm_ui$U.lo mpz/primorial_ui$U.lo \
mpz/pprime_p$U.lo mpz/random$U.lo mpz/random2$U.lo \
- mpz/realloc$U.lo mpz/realloc2$U.lo mpz/remove$U.lo \
+ mpz/realloc$U.lo mpz/realloc2$U.lo mpz/remove$U.lo mpz/roinit_n$U.lo \
mpz/root$U.lo mpz/rootrem$U.lo mpz/rrandomb$U.lo mpz/scan0$U.lo \
mpz/scan1$U.lo mpz/set$U.lo mpz/set_d$U.lo mpz/set_f$U.lo \
mpz/set_q$U.lo mpz/set_si$U.lo mpz/set_str$U.lo mpz/set_ui$U.lo \
diff -r 32f62aca5666 gmp-h.in
--- a/gmp-h.in Thu Feb 21 11:03:04 2013 +0100
+++ b/gmp-h.in Sat Feb 23 17:38:59 2013 +0100
@@ -1117,6 +1117,22 @@ using std::FILE;
#define mpz_eor __gmpz_xor
__GMP_DECLSPEC void mpz_xor (mpz_ptr, mpz_srcptr, mpz_srcptr);
+#define mpz_limbs_read __gmpz_limbs_read
+__GMP_DECLSPEC mp_srcptr mpz_limbs_read (mpz_srcptr);
+
+#define mpz_limbs_write __gmpz_limbs_write
+__GMP_DECLSPEC mp_ptr mpz_limbs_write (mpz_ptr, mp_size_t);
+
+#define mpz_limbs_modify __gmpz_limbs_modify
+__GMP_DECLSPEC mp_ptr mpz_limbs_modify (mpz_ptr, mp_size_t);
+
+#define mpz_limbs_finish __gmpz_limbs_finish
+__GMP_DECLSPEC void mpz_limbs_finish (mpz_ptr, mp_size_t);
+
+#define mpz_roinit_n __gmpz_roinit_n
+__GMP_DECLSPEC mpz_srcptr mpz_roinit_n (mpz_ptr, mp_srcptr, mp_size_t);
+
+#define MPZ_ROINIT_N(xp, xs) {{0, (xs),(xp) }}
/**************** Rational (i.e. Q) routines. ****************/
diff -r 32f62aca5666 mpz/Makefile.am
--- a/mpz/Makefile.am Thu Feb 21 11:03:04 2013 +0100
+++ b/mpz/Makefile.am Sat Feb 23 17:38:59 2013 +0100
@@ -43,12 +43,13 @@ libmpz_la_SOURCES = aors.h aors_ui.h fit
import.c init.c init2.c inits.c inp_raw.c inp_str.c \
invert.c ior.c iset.c iset_d.c iset_si.c iset_str.c iset_ui.c \
jacobi.c kronsz.c kronuz.c kronzs.c kronzu.c \
- lcm.c lcm_ui.c lucnum_ui.c lucnum2_ui.c mfac_uiui.c millerrabin.c \
+ lcm.c lcm_ui.c limbs_read.c limbs_write.c limbs_modify.c limbs_finish.c \
+ lucnum_ui.c lucnum2_ui.c mfac_uiui.c millerrabin.c \
mod.c mul.c mul_2exp.c mul_si.c mul_ui.c n_pow_ui.c neg.c nextprime.c \
oddfac_1.c \
out_raw.c out_str.c perfpow.c perfsqr.c popcount.c pow_ui.c powm.c \
powm_sec.c powm_ui.c pprime_p.c prodlimbs.c primorial_ui.c random.c random2.c \
- realloc.c realloc2.c remove.c root.c rootrem.c rrandomb.c \
+ realloc.c realloc2.c remove.c roinit_n.c root.c rootrem.c rrandomb.c \
scan0.c scan1.c set.c set_d.c set_f.c set_q.c set_si.c set_str.c \
set_ui.c setbit.c size.c sizeinbase.c sqrt.c sqrtrem.c sub.c sub_ui.c \
swap.c tdiv_ui.c tdiv_q.c tdiv_q_2exp.c tdiv_q_ui.c tdiv_qr.c \
diff -r 32f62aca5666 mpz/limbs_finish.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/mpz/limbs_finish.c Sat Feb 23 17:38:59 2013 +0100
@@ -0,0 +1,29 @@
+/* mpz_finish_limbs -- Update mpz after writing to the limb array.
+
+Copyright 2013 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP 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 3 of the License, or (at your
+option) any later version.
+
+The GNU MP 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 GNU MP Library. If not, see http://www.gnu.org/licenses/. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+
+void
+mpz_limbs_finish (mpz_ptr x, mp_size_t xs)
+{
+ mp_size_t xn = ABS(xs);
+ MPN_NORMALIZE (PTR (x), xn);
+ SIZ (x) = xs < 0 ? -xn : xn;
+}
diff -r 32f62aca5666 mpz/limbs_modify.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/mpz/limbs_modify.c Sat Feb 23 17:38:59 2013 +0100
@@ -0,0 +1,28 @@
+/* mpz_limbs_modify -- Read-and-modify access to the mpn-style limb array.
+
+Copyright 2013 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP 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 3 of the License, or (at your
+option) any later version.
+
+The GNU MP 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 GNU MP Library. If not, see http://www.gnu.org/licenses/. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+
+mp_ptr
+mpz_limbs_modify (mpz_ptr x, mp_size_t n)
+{
+ ASSERT (n > 0);
+ return MPZ_REALLOC (x, n);
+}
diff -r 32f62aca5666 mpz/limbs_read.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/mpz/limbs_read.c Sat Feb 23 17:38:59 2013 +0100
@@ -0,0 +1,27 @@
+/* mpz_limbs_read -- Read access to the mpn-style limb array.
+
+Copyright 2013 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP 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 3 of the License, or (at your
+option) any later version.
+
+The GNU MP 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 GNU MP Library. If not, see http://www.gnu.org/licenses/. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+
+mp_srcptr
+mpz_limbs_read (mpz_srcptr x)
+{
+ return PTR(x);
+}
diff -r 32f62aca5666 mpz/limbs_write.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/mpz/limbs_write.c Sat Feb 23 17:38:59 2013 +0100
@@ -0,0 +1,28 @@
+/* mpz_limbs_write -- Write access to the mpn-style limb array.
+
+Copyright 2013 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP 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 3 of the License, or (at your
+option) any later version.
+
+The GNU MP 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 GNU MP Library. If not, see http://www.gnu.org/licenses/. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+
+mp_ptr
+mpz_limbs_write (mpz_ptr x, mp_size_t n)
+{
+ ASSERT (n > 0);
+ return MPZ_NEWALLOC (x, n);
+}
diff -r 32f62aca5666 mpz/roinit_n.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/mpz/roinit_n.c Sat Feb 23 17:38:59 2013 +0100
@@ -0,0 +1,30 @@
+/* mpz_roinit_n -- Initialize mpz with read-only limb array.
+
+Copyright 2013 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP 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 3 of the License, or (at your
+option) any later version.
+
+The GNU MP 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 GNU MP Library. If not, see http://www.gnu.org/licenses/. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+
+mpz_srcptr
+mpz_roinit_n (mpz_ptr x, mp_srcptr xp, mp_size_t xs)
+{
+ ALLOC (x) = 0;
+ SIZ (x) = xs;
+ PTR (x) = xp;
+ return x;
+}
diff -r 32f62aca5666 tests/mpz/Makefile.am
--- a/tests/mpz/Makefile.am Thu Feb 21 11:03:04 2013 +0100
+++ b/tests/mpz/Makefile.am Sat Feb 23 17:38:59 2013 +0100
@@ -30,7 +30,7 @@ check_PROGRAMS = t-addsub t-cmp t-mul t-
t-fac_ui t-mfac_uiui t-primorial_ui t-fib_ui t-lucnum_ui t-scan t-fits \
t-divis t-divis_2exp t-cong t-cong_2exp t-sizeinbase t-set_str \
t-aorsmul t-cmp_d t-cmp_si t-hamdist t-oddeven t-popcount t-set_f \
- t-io_raw t-import t-export t-pprime_p t-nextprime t-remove
+ t-io_raw t-import t-export t-pprime_p t-nextprime t-remove t-limbs
TESTS = $(check_PROGRAMS)
diff -r 32f62aca5666 tests/mpz/t-limbs.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tests/mpz/t-limbs.c Sat Feb 23 17:38:59 2013 +0100
@@ -0,0 +1,171 @@
+/* Test mpz_limbs_* functions
+
+Copyright 2013 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library test suite.
+
+The GNU MP Library test suite 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 3 of the License,
+or (at your option) any later version.
+
+The GNU MP Library test suite 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
+the GNU MP Library test suite. If not, see http://www.gnu.org/licenses/. */
+
+#include <stdlib.h>
+#include <stdio.h>
+
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "tests.h"
+
+#define COUNT 100
+#define BITSIZE 500
+
+/* Like mpz_add. For simplicity, support positive inputs only. */
+static void
+alt_add (mpz_ptr r, mpz_srcptr a, mpz_srcptr b)
+{
+ mp_size_t an = mpz_size (a);
+ mp_size_t bn = mpz_size (b);
+ mp_ptr rp;
+
+ ASSERT (an > 0);
+ ASSERT (bn > 0);
+ if (an < bn)
+ {
+ MP_SIZE_T_SWAP (an, bn);
+ MPZ_SRCPTR_SWAP (a, b);
+ }
+ rp = mpz_limbs_modify (r, an + 1);
+ rp[an] = mpn_add (rp, mpz_limbs_read (a), an, mpz_limbs_read (b), bn);
+ mpz_limbs_finish (r, an + 1);
+}
+
+static void
+check_funcs (const char *name,
+ void (*f)(mpz_ptr, mpz_srcptr, mpz_srcptr),
+ void (*ref_f)(mpz_ptr, mpz_srcptr, mpz_srcptr),
+ mpz_srcptr a, mpz_srcptr b)
+{
+ mpz_t r, ref;
+ mpz_inits (r, ref, NULL);
+
+ ref_f (ref, a, b);
+ MPZ_CHECK_FORMAT (ref);
+ f (r, a, b);
+ MPZ_CHECK_FORMAT (r);
+
+ if (mpz_cmp (r, ref) != 0)
+ {
+ printf ("%s failed, abits %u, bbits %u\n",
+ (unsigned) mpz_sizeinbase (a, 2),
+ (unsigned) mpz_sizeinbase (b, 2));
+ gmp_printf ("a = %Zx\n", a);
+ gmp_printf ("b = %Zx\n", b);
+ gmp_printf ("r = %Zx (bad)\n", r);
+ gmp_printf ("ref = %Zx\n", ref);
+ }
+ mpz_clears (r, ref, NULL);
+}
+
+static void
+check_add (void)
+{
+ gmp_randstate_ptr rands = RANDS;
+ mpz_t bs, a, b;
+ unsigned i;
+ mpz_inits (bs, a, b, NULL);
+ for (i = 0; i < COUNT; i++)
+ {
+ mpz_rrandomb (bs, rands, 32);
+ mpz_rrandomb (a, rands, 1 + mpz_get_ui (bs) % BITSIZE);
+ mpz_rrandomb (bs, rands, 32);
+ mpz_rrandomb (b, rands, 1 + mpz_get_ui (bs) % BITSIZE);
+
+ check_funcs ("add", alt_add, mpz_add, a, b);
+ }
+ mpz_clears (bs, a, b, NULL);
+}
+
+static void
+alt_mul (mpz_ptr r, mpz_srcptr a, mpz_srcptr b)
+{
+ mp_size_t an = mpz_size (a);
+ mp_size_t bn = mpz_size (b);
+ mp_srcptr ap, bp;
+ mp_ptr rp;
+ TMP_DECL;
+
+ TMP_MARK;
+
+ ASSERT (an > 0);
+ ASSERT (bn > 0);
+ if (an < bn)
+ {
+ MP_SIZE_T_SWAP (an, bn);
+ MPZ_SRCPTR_SWAP (a, b);
+ }
+ /* NOTE: This copying seems unnecessary; better to allocate new
+ result area, and free the old area when done. */
+ if (r == a)
+ {
+ mp_ptr tp = TMP_ALLOC_LIMBS (an);
+ MPN_COPY (tp, mpz_limbs_read (a), an);
+ ap = tp;
+ bp = (a == b) ? ap : mpz_limbs_read (b);
+ }
+ else if (r == b)
+ {
+ mp_ptr tp = TMP_ALLOC_LIMBS (bn);
+ MPN_COPY (tp, mpz_limbs_read (b), bn);
+ bp = tp;
+ ap = mpz_limbs_read (a);
+ }
+ else
+ {
+ ap = mpz_limbs_read (a);
+ bp = mpz_limbs_read (b);
+ }
+ mpn_mul (mpz_limbs_write (r, an + bn),
+ ap, an, bp, bn);
+
+ mpz_limbs_finish (r, an + bn);
+}
+
+void
+check_mul (void)
+{
+ gmp_randstate_ptr rands = RANDS;
+ mpz_t bs, a, b, r, ref;
+ unsigned i;
+ mpz_inits (bs, a, b, NULL);
+ for (i = 0; i < COUNT; i++)
+ {
+ mpz_rrandomb (bs, rands, 32);
+ mpz_rrandomb (a, rands, 1 + mpz_get_ui (bs) % BITSIZE);
+ mpz_rrandomb (bs, rands, 32);
+ mpz_rrandomb (b, rands, 1 + mpz_get_ui (bs) % BITSIZE);
+
+ check_funcs ("mul", alt_mul, mpz_mul, a, b);
+ }
+ mpz_clears (bs, a, b, NULL);
+}
+
+int
+main (int argc, char *argv[])
+{
+ tests_start ();
+ tests_end ();
+
+ check_add ();
+ check_mul ();
+
+ return 0;
+
+}
--
Niels Möller. PGP-encrypted email is preferred. Keyid C0B98E26.
Internet email is subject to wholesale government surveillance.
More information about the gmp-devel
mailing list