Dear all,
I am new to this list. I am a regular developer of the SageMath project
[1] in which we use mpz_t and mpq_t as our based type for integers and
rationals. We often have to compare integers with rationals and I came
to the conclusion that it would be convenient to have a GMP function
int mpq_cmp_z(mpq_t, mpz_t)
See the down stream ticket [2]. As we ship mpir by default, I made a
concrete proposal to them [3] and they kindly asked me to make the
proposal directly to GMP since they do want to keep their interface
consistent with the GMP project.
If you think that it would be a good addition, you might have a look in
the attached patch that was tested on the development version of GMP.
Best,
Vincent
[1] http://www.sagemath.org
[2] http://trac.sagemath.org/ticket/18304
[3] https://github.com/wbhart/mpir/pull/157
diff -r e5a96e5eb06d Makefile.am
--- a/Makefile.am Fri Aug 14 17:38:58 2015 +0200
+++ b/Makefile.am Sat Aug 15 12:35:02 2015 +0200
@@ -214,7 +214,7 @@
MPQ_OBJECTS = mpq/abs$U.lo mpq/aors$U.lo \
mpq/canonicalize$U.lo mpq/clear$U.lo mpq/clears$U.lo \
- mpq/cmp$U.lo mpq/cmp_si$U.lo mpq/cmp_ui$U.lo mpq/div$U.lo \
+ mpq/cmp$U.lo mpq/cmp_si$U.lo mpq/cmp_ui$U.lo mpq/cmp_z$U.lo mpq/div$U.lo \
mpq/get_d$U.lo mpq/get_den$U.lo mpq/get_num$U.lo mpq/get_str$U.lo \
mpq/init$U.lo mpq/inits$U.lo mpq/inp_str$U.lo mpq/inv$U.lo \
mpq/md_2exp$U.lo mpq/mul$U.lo mpq/neg$U.lo mpq/out_str$U.lo \
diff -r e5a96e5eb06d mpq/cmp_z.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/mpq/cmp_z.c Sat Aug 15 12:35:02 2015 +0200
@@ -0,0 +1,97 @@
+/* mpq_cmp_z(u,v) -- Compare U, V. Return positive, zero, or negative
+ based on if U > V, U == V, or U < V.
+
+Copyright 2015 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 2.1 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; see the file COPYING.LIB. If not, write to
+the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+MA 02110-1301, USA. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+
+int
+mpq_cmp_z(mpq_srcptr op, mpz_srcptr v)
+{
+ mp_size_t num_size = SIZ(NUM(op));
+ mp_size_t den_size = SIZ(DEN(op));
+ mp_size_t v_size = SIZ(v);
+ mp_size_t tmp_size;
+ mp_size_t num_sign;
+ mp_ptr tmp_ptr;
+ int cc;
+ TMP_DECL;
+
+ ASSERT (den_size > 0);
+
+ /* 0. Check signs */
+ if (num_size == 0)
+ return -v_size;
+ if (v_size == 0)
+ return num_size;
+ if ((num_size ^ v_size) < 0) /* different signs */
+ return num_size;
+ if (!mpz_cmp_ui(mpq_denref(op), 1)) /* denominator 1 */
+ return mpz_cmp(mpq_numref(op), v);
+
+ /* 1. Check to see if we can tell which operand is larger by just looking at the number of limbs */
+ num_sign = num_size;
+ num_size = ABS(num_size);
+ v_size = ABS(v_size);
+ tmp_size = v_size + den_size;
+
+ /* den*v is either size tmp_size or tmp_size-1 */
+ if(num_size > tmp_size)
+ return num_sign;
+ if(tmp_size > num_size + 1)
+ return -num_sign;
+
+ /* 2 . Same, but compare the number of significant bits. */
+ {
+ int cnt1, cnt2;
+ unsigned long int bits1, bits2;
+
+ count_leading_zeros(cnt1, PTR(NUM(op))[num_size-1]);
+ bits1 = num_size * GMP_NUMB_BITS - cnt1 + GMP_NAIL_BITS;
+
+ count_leading_zeros(cnt1, PTR(v)[v_size-1]);
+ count_leading_zeros(cnt2, PTR(DEN(op))[den_size-1]);
+ bits2 = tmp_size * GMP_NUMB_BITS - cnt1 - cnt2 + 2 * GMP_NAIL_BITS;
+
+ if (bits1 > bits2 + 1)
+ return num_sign;
+ if (bits2 > bits1 + 1)
+ return -num_sign;
+ }
+
+ /* 3. Finally, cross multiply and compare. */
+ TMP_MARK;
+ tmp_ptr = (mp_ptr) TMP_ALLOC(num_size * GMP_LIMB_BYTES);
+
+ if( v_size >= den_size)
+ tmp_size -= 0 == mpn_mul(tmp_ptr,
+ v->_mp_d, v_size,
+ op->_mp_den._mp_d, den_size);
+ else
+ tmp_size -= 0 == mpn_mul(tmp_ptr,
+ op->_mp_den._mp_d, den_size,
+ v->_mp_d, v_size);
+ cc = num_size - tmp_size != 0
+ ? num_size - tmp_size : mpn_cmp(op->_mp_num._mp_d, tmp_ptr, tmp_size);
+ TMP_FREE;
+ return num_sign < 0 ? -cc: cc;
+}
diff -r e5a96e5eb06d tests/mpq/Makefile.am
--- a/tests/mpq/Makefile.am Fri Aug 14 17:38:58 2015 +0200
+++ b/tests/mpq/Makefile.am Sat Aug 15 12:35:02 2015 +0200
@@ -21,7 +21,7 @@
AM_CPPFLAGS = -I$(top_srcdir) -I$(top_srcdir)/tests
LDADD = $(top_builddir)/tests/libtests.la $(top_builddir)/libgmp.la
-check_PROGRAMS = t-aors t-cmp t-cmp_ui t-cmp_si t-equal t-get_d t-get_str \
+check_PROGRAMS = t-aors t-cmp t-cmp_ui t-cmp_si t-cmp_z t-equal t-get_d t-get_str \
t-inp_str t-inv t-md_2exp t-set_f t-set_str io reuse
TESTS = $(check_PROGRAMS)
diff -r e5a96e5eb06d tests/mpq/t-cmp_z.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/tests/mpq/t-cmp_z.c Sat Aug 15 12:35:02 2015 +0200
@@ -0,0 +1,91 @@
+/* Test mpq_cmp_z.
+ *
+Copyright 2015 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 https://www.gnu.org/licenses/. */
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "tests.h"
+
+#define SGN(x) ((x) < 0 ? -1 : (x) > 0 ? 1 : 0)
+
+int ref_mpq_cmp_z(mpq_t a, mpz_t b)
+{
+ mpz_t bi;
+ int cc;
+ mpz_init(bi);
+
+ mpz_mul(bi, b, DEN(a));
+ cc = mpz_cmp(NUM(a), bi);
+ mpz_clear(bi);
+ return cc;
+}
+
+#ifndef SIZE
+#define SIZE 8 /* increasing this lowers the probability of finding an error */
+#endif
+
+int
+main (int argc, char **argv)
+{
+ mpq_t a;
+ mpz_t b;
+ mp_size_t size;
+ int reps = 10000;
+ int i;
+ int cc, ccref;
+
+ tests_start ();
+
+ if (argc == 2)
+ reps = atoi (argv[1]);
+
+ mpq_init (a);
+ mpz_init (b);
+
+ for (i = 0; i < reps; i++)
+ {
+ size = urandom () % SIZE - SIZE/2;
+ mpz_random2 (NUM (a), size);
+ do
+ {
+ size = urandom () % SIZE - SIZE/2;
+ mpz_random2 (DEN (a), size);
+ }
+ while (mpz_cmp_ui (DEN (a), 0) == 0);
+
+ size = urandom () % SIZE - SIZE/2;
+ mpz_random2 (b, size);
+
+ mpq_canonicalize (a);
+
+ ccref = ref_mpq_cmp_z (a, b);
+ cc = mpq_cmp_z (a, b);
+
+ if (SGN (ccref) != SGN (cc))
+ abort ();
+ }
+
+ mpq_clear (a);
+ mpz_clear (b);
+
+ tests_end ();
+ exit (0);
+}
_______________________________________________
gmp-devel mailing list
gmp-devel@gmplib.org
https://gmplib.org/mailman/listinfo/gmp-devel
