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

[mercurial at gmplib.org]

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:
Further edit comments.

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,