[Gmp-commit] /var/hg/gmp: 2 new changesets

mercurial at gmplib.org mercurial at gmplib.org
Sat Oct 27 13:21:25 CEST 2012

```details:   /var/hg/gmp/rev/5d6dd71e0a0a
changeset: 15090:5d6dd71e0a0a
user:      Torbjorn Granlund <tege at gmplib.org>
date:      Sat Oct 27 11:43:28 2012 +0200
description:
Whitespace cleanup.

details:   /var/hg/gmp/rev/0b25e4cfd719
changeset: 15091:0b25e4cfd719
user:      Torbjorn Granlund <tege at gmplib.org>
date:      Sat Oct 27 13:21:09 2012 +0200
description:

diffstat:

mpn/generic/perfpow.c |  56 +++++++++++++++++++++++++++-----------------------
1 files changed, 30 insertions(+), 26 deletions(-)

diffs (87 lines):

diff -r 755243e9118f -r 0b25e4cfd719 mpn/generic/perfpow.c
--- a/mpn/generic/perfpow.c	Sat Oct 27 10:50:42 2012 +0200
+++ b/mpn/generic/perfpow.c	Sat Oct 27 13:21:09 2012 +0200
@@ -106,33 +106,31 @@
return ans;
}

-/*
-   Computes rp such that rp^k * yp = 1 (mod 2^b).
-   Algorithm:
-       Apply Hensel lifting repeatedly, each time
-       doubling (approx.) the number of known bits in rp.
+/* Compute r such that r^k * y = 1 (mod 2^b).

-   Uses the iteration
+   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.

-     r' <-- k^{-1} ((k+1) r - r^{k+1} y) (mod 2^b)
+   Works just for odd k.  Else the Hensel lifting degenerates.

-   FIXME: I think we get *exactly* doubled precision in each
-   iteration. Could be simplified to do precision book-keeping in
-   limbs rather than bits.
-
-   FIXME: Rewrite iteration as
+   FIXME:
+     (1) Simplify to do precision book-keeping in limbs rather than bits.

-     r' <-- r - k^{-1} r (r^k y - 1)
+     (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.

-   and take advantage of the zero low part of r^k y - 1.
+     (3) Use wrap-around trick.

-   If we prefer to get y^{1/k} rather than y^{-1/k}, we could also
-   switch to the iteration
+     (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.
+   which converges to y^{1/n - 1}, and we then get y^{1/n} by a single mullo.
*/
static void
binv_root (mp_ptr rp, mp_srcptr yp,
@@ -171,19 +169,25 @@
return;
}

-/*
-   Computes rp such that rp^2 * yp = 1 (mod 2^{b+1}).
-   Returns non-zero if such an integer rp exists.
+/* Compute r such that r^2 * y = 1 (mod 2^{b+1}).
+   Return non-zero if such an integer r exists.

-   Uses the iteration
+   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.

-     r' <-- (3r - r^3 y) / 2
+   FIXME:
+     (1) Simplify to do precision book-keeping in limbs rather than bits.

-   FIXME: Rewrite as
+     (2) Rewrite iteration as
+	   r' <-- r - r (r^2 y - 1) / 2
+	 and take advantage of zero low part of r^2 y - 1.

-     r' <-- r - r (r^2 y - 1) / 2
+     (3) Use wrap-around trick.

-   and take advantage of zero low part of r^2 y - 1.
+     (4) Use a small table to get starting value.
+
*/
static int
binv_sqroot (mp_ptr rp, mp_srcptr yp,
```