[Gmp-commit] /var/hg/gmp: 5 new changesets
mercurial at gmplib.org
mercurial at gmplib.org
Sun Oct 28 15:51:29 CET 2012
details: /var/hg/gmp/rev/b30d41e3a372
changeset: 15094:b30d41e3a372
user: Torbjorn Granlund <tege at gmplib.org>
date: Sun Oct 28 13:44:03 2012 +0100
description:
Whitespace cleanup.
details: /var/hg/gmp/rev/5d842e6151fd
changeset: 15095:5d842e6151fd
user: Torbjorn Granlund <tege at gmplib.org>
date: Sun Oct 28 15:48:24 2012 +0100
description:
* mpn/generic/brootinv.c: New file, code from overhauled binv_root.
* mpn/generic/bsqrtinv.c: New file, code from overhauled binv_sqroot.
* mpn/generic/bsqrt.c: New file.
details: /var/hg/gmp/rev/e8dec0678826
changeset: 15096:e8dec0678826
user: Torbjorn Granlund <tege at gmplib.org>
date: Sun Oct 28 15:49:17 2012 +0100
description:
Declare new functions.
details: /var/hg/gmp/rev/161eb8599bd0
changeset: 15097:161eb8599bd0
user: Torbjorn Granlund <tege at gmplib.org>
date: Sun Oct 28 15:50:54 2012 +0100
description:
Overhaul; (binv_root, binv_sqroot): Remove.
details: /var/hg/gmp/rev/7a309cd32a28
changeset: 15098:7a309cd32a28
user: Torbjorn Granlund <tege at gmplib.org>
date: Sun Oct 28 15:51:08 2012 +0100
description:
*** empty log message ***
diffstat:
ChangeLog | 8 +
configure.in | 2 +-
gmp-impl.h | 9 +
mpn/generic/broot.c | 14 +-
mpn/generic/brootinv.c | 87 ++++++++++++++
mpn/generic/bsqrt.c | 37 ++++++
mpn/generic/bsqrtinv.c | 94 +++++++++++++++
mpn/generic/perfpow.c | 289 +++++++++++++-----------------------------------
8 files changed, 320 insertions(+), 220 deletions(-)
diffs (truncated from 795 to 300 lines):
diff -r 015a9fbab1ef -r 7a309cd32a28 ChangeLog
--- a/ChangeLog Sun Oct 28 13:21:34 2012 +0100
+++ b/ChangeLog Sun Oct 28 15:51:08 2012 +0100
@@ -1,5 +1,13 @@
2012-10-28 Torbjorn Granlund <tege at gmplib.org>
+ * configure.in (gmp_mpn_functions): Add new files.
+ * gmp-impl.h: Declare new functions.
+ * mpn/generic/perfpow.c: Overhaul.
+ (binv_root, binv_sqroot): Remove.
+ * mpn/generic/brootinv.c: New file, code from overhauled binv_root.
+ * mpn/generic/bsqrtinv.c: New file, code from overhauled binv_sqroot.
+ * mpn/generic/bsqrt.c: New file.
+
* tests/mpn/t-broot.c: Add a forgotten TMP_MARK.
2012-10-28 Niels Möller <nisse at lysator.liu.se>
diff -r 015a9fbab1ef -r 7a309cd32a28 configure.in
--- a/configure.in Sun Oct 28 13:21:34 2012 +0100
+++ b/configure.in Sun Oct 28 15:51:08 2012 +0100
@@ -2679,7 +2679,7 @@
sbpi1_bdiv_q sbpi1_bdiv_qr \
dcpi1_bdiv_q dcpi1_bdiv_qr \
mu_bdiv_q mu_bdiv_qr \
- bdiv_q bdiv_qr broot \
+ bdiv_q bdiv_qr broot brootinv bsqrt bsqrtinv \
divexact bdiv_dbm1c redc_1 redc_2 redc_n powm powlo powm_sec \
trialdiv remove \
and_n andn_n nand_n ior_n iorn_n nior_n xor_n xnor_n \
diff -r 015a9fbab1ef -r 7a309cd32a28 gmp-impl.h
--- a/gmp-impl.h Sun Oct 28 13:21:34 2012 +0100
+++ b/gmp-impl.h Sun Oct 28 15:51:08 2012 +0100
@@ -1632,6 +1632,15 @@
#define mpn_broot __MPN(broot)
__GMP_DECLSPEC void mpn_broot (mp_ptr, mp_srcptr, mp_size_t, mp_limb_t);
+#define mpn_brootinv __MPN(brootinv)
+__GMP_DECLSPEC void mpn_brootinv (mp_ptr, mp_srcptr, mp_bitcnt_t, mp_limb_t, mp_ptr);
+
+#define mpn_bsqrt __MPN(bsqrt)
+__GMP_DECLSPEC void mpn_bsqrt (mp_ptr, mp_srcptr, mp_bitcnt_t, mp_ptr);
+
+#define mpn_bsqrtinv __MPN(bsqrtinv)
+__GMP_DECLSPEC int mpn_bsqrtinv (mp_ptr, mp_srcptr, mp_bitcnt_t, mp_ptr);
+
#if defined (_CRAY)
#define MPN_COPY_INCR(dst, src, n) \
do { \
diff -r 015a9fbab1ef -r 7a309cd32a28 mpn/generic/broot.c
--- a/mpn/generic/broot.c Sun Oct 28 13:21:34 2012 +0100
+++ b/mpn/generic/broot.c Sun Oct 28 15:51:08 2012 +0100
@@ -58,7 +58,7 @@
static void
mpn_broot_invm1 (mp_ptr rp, mp_srcptr ap, mp_size_t n, mp_limb_t k)
{
- mp_size_t sizes[GMP_LIMB_BITS * 2];
+ mp_size_t sizes[GMP_LIMB_BITS * 2];
mp_ptr akm1, tp, rnp, ep, scratch;
mp_limb_t a0, r0, km1, kinv;
mp_size_t rn;
@@ -72,7 +72,7 @@
ASSERT (k >= 3);
TMP_MARK;
-
+
akm1 = TMP_ALLOC_LIMBS (4*n);
tp = akm1 + n;
@@ -89,8 +89,8 @@
/* 4 bits: a^{1/k - 1} (mod 16):
- a % 8
- 1 3 5 7
+ a % 8
+ 1 3 5 7
k%4 +-------
1 |1 1 1 1
3 |1 9 9 1
@@ -119,7 +119,7 @@
}
/* FIXME: Special case for two limb iteration. */
- rnp = TMP_ALLOC_LIMBS (2*n);
+ rnp = TMP_ALLOC_LIMBS (2*n);
ep = rnp + n;
/* FIXME: Possible to this on the fly with some bit fiddling. */
@@ -133,7 +133,7 @@
/* Compute x^k. What's the best way to handle the doubled
precision? */
MPN_ZERO (rp + rn, sizes[i] - rn);
- mpn_powlo (rnp, rp, &k, 1, sizes[i], tp);
+ mpn_powlo (rnp, rp, &k, 1, sizes[i], tp);
/* Multiply by a^{k-1}. Can use wraparound; low part is
000...01. */
@@ -180,5 +180,5 @@
mpn_broot_invm1 (tp, ap, n, k);
mpn_mullo_n (rp, tp, ap, n);
- TMP_FREE;
+ TMP_FREE;
}
diff -r 015a9fbab1ef -r 7a309cd32a28 mpn/generic/brootinv.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/mpn/generic/brootinv.c Sun Oct 28 15:51:08 2012 +0100
@@ -0,0 +1,87 @@
+/* mpn_brootinv, compute r such that r^k * y = 1 (mod 2^b).
+
+ Contributed to the GNU project by Martin Boij (as part of perfpow.c).
+
+Copyright 2009, 2010, 2012 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"
+
+/* Compute r such that r^k * y = 1 (mod 2^b).
+
+ Iterates
+ r' <-- k^{-1} ((k+1) r - r^{k+1} y) (mod 2^b)
+ using Hensel lifting, each time doubling the number of known bits in r.
+
+ Works just for odd k. Else the Hensel lifting degenerates.
+
+ FIXME:
+ (1) Simplify to do precision book-keeping in limbs rather than bits.
+
+ (2) Rewrite iteration as
+ r' <-- r - k^{-1} r (r^k y - 1)
+ and take advantage of the zero low part of r^k y - 1.
+
+ (3) Use wrap-around trick.
+
+ (4) Use a small table to get starting value.
+
+ If we prefer to get y^{1/k} rather than y^{-1/k}, we could also switch to
+ the iteration
+
+ r' <-- r - k^{-1} r (r^k y^{k-1} - 1)
+
+ which converges to y^{1/n - 1}, and we then get y^{1/n} by a single mullo.
+*/
+void
+mpn_brootinv (mp_ptr rp, mp_srcptr yp, mp_bitcnt_t bnb, mp_limb_t k, mp_ptr tp)
+{
+ mp_ptr tp2, tp3;
+ mp_limb_t kinv, k2;
+ mp_size_t bn, order[GMP_LIMB_BITS + 1];
+ int i, d;
+
+ ASSERT (bnb > 0);
+ ASSERT ((k & 1) != 0);
+
+ bn = 1 + (bnb - 1) / GMP_LIMB_BITS;
+
+ tp2 = tp + bn;
+ tp3 = tp + 2 * bn;
+ k2 = k + 1;
+
+ binvert_limb (kinv, k);
+
+ d = 0;
+ for (; bnb != 1; bnb = (bnb + 1) >> 1)
+ order[d++] = 1 + (bnb - 1) / GMP_LIMB_BITS;
+
+ rp[0] = 1;
+ for (i = d - 1; i >= 0; i--)
+ {
+ bn = order[i];
+
+ mpn_mul_1 (tp, rp, bn, k2);
+
+ mpn_powlo (tp2, rp, &k2, 1, bn, tp3);
+ mpn_mullo_n (rp, yp, tp2, bn);
+
+ mpn_sub_n (tp2, tp, rp, bn);
+ mpn_pi1_bdiv_q_1 (rp, tp2, bn, k, kinv, 0);
+ }
+}
diff -r 015a9fbab1ef -r 7a309cd32a28 mpn/generic/bsqrt.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/mpn/generic/bsqrt.c Sun Oct 28 15:51:08 2012 +0100
@@ -0,0 +1,37 @@
+/* mpn_bsqrt, a^{1/2} (mod 2^n).
+
+Copyright 2009, 2010, 2012 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
+mpn_bsqrt (mp_ptr rp, mp_srcptr ap, mp_bitcnt_t nb, mp_ptr tp)
+{
+ mp_ptr sp;
+ mp_size_t n;
+
+ ASSERT (nb > 0);
+
+ n = nb / GMP_NUMB_BITS;
+ sp = tp + n;
+
+ mpn_bsqrtinv (sp, ap, nb, tp);
+ mpn_mullo_n (rp, sp, ap, n);
+}
diff -r 015a9fbab1ef -r 7a309cd32a28 mpn/generic/bsqrtinv.c
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/mpn/generic/bsqrtinv.c Sun Oct 28 15:51:08 2012 +0100
@@ -0,0 +1,94 @@
+/* mpn_bsqrtinv, compute r such that r^2 * y = 1 (mod 2^{b+1}).
+
+ Contributed to the GNU project by Martin Boij (as part of perfpow.c).
+
+Copyright 2009, 2010, 2012 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"
+
+/* Compute r such that r^2 * y = 1 (mod 2^{b+1}).
+ Return non-zero if such an integer r exists.
+
+ Iterates
+ r' <-- (3r - r^3 y) / 2
+ using Hensel lifting. Since we divide by two, the Hensel lifting is
+ somewhat degenerates. Therefore, we lift from 2^b to 2^{b+1}-1.
+
+ FIXME:
+ (1) Simplify to do precision book-keeping in limbs rather than bits.
+
+ (2) Rewrite iteration as
+ r' <-- r - r (r^2 y - 1) / 2
+ and take advantage of zero low part of r^2 y - 1.
+
+ (3) Use wrap-around trick.
+
+ (4) Use a small table to get starting value.
+*/
+int
+mpn_bsqrtinv (mp_ptr rp, mp_srcptr yp, mp_bitcnt_t bnb, mp_ptr tp)
+{
+ mp_ptr tp2, tp3;
+ mp_limb_t k;
+ mp_size_t bn, order[GMP_LIMB_BITS + 1];
+ int i, d;
+
+ ASSERT (bnb > 0);
+
+ bn = 1 + bnb / GMP_LIMB_BITS;
+
+ tp2 = tp + bn;
+ tp3 = tp + 2 * bn;
+ k = 3;
+
More information about the gmp-commit
mailing list