[Gmpcommit] /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:
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 bookkeeping in
 limbs rather than bits.

 FIXME: Rewrite iteration as
+ FIXME:
+ (1) Simplify to do precision bookkeeping 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 wraparound 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^{k1}  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 nonzero if such an integer rp exists.
+/* Compute r such that r^2 * y = 1 (mod 2^{b+1}).
+ Return nonzero 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 bookkeeping 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 wraparound 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,
More information about the gmpcommit
mailing list