Bill, I don't see your post to me at google group, I think the address
was bad
to mpir list. So I post your last message also with my 3 answers.
And uploaded your new code perfpow.c at google group.
I've attached a first attempt at an mpn version of Robert Gerbicz'
fast perfect power test code, minus some functions I need to make it
work. It doesn't compile yet, and I have probably introduced oodles of
bugs to the code. However I'm going to stop working on it for now,
until we can get some of the functions described below, implemented at
the mpn level in MPIR.
My impression of the code is that this is a *really* clever algorithm
and should beat pretty much anything out there when coded efficiently.
Therefore it is really important that we get this into MPIR.
I did correct what I think was a minor bug. On what is now line 328 I
replaced smallprimes[j]*smallprimes[j] < i + bound_for_p with
smallprimes[j]*smallprimes[j] < i + sq, which I think is correct.
Robert can you verify this!
RE:
Yes, you are correct, (we are sieving in the [i,i+sq-1] interval).
I also replaced the primes q with just unsigned longs. I can imagine
one day someone having hardware to test numbers n of N = 10*2^40
*bits* for whether they are perfect powers. As the space requirement
for this function is approximately N log N, this seems to be a
sensible bound for the forseeable future. To perform the test in this
case we need primes up to 2^40. At one point we look for primes of the
form q = p*2k+1, starting with q = 2*p + 1. But we only generate
primes q = p*2k + 1 for as long as p^k <= bound_for_p < B (i.e. fits
in a limb). In general q will be much smaller than p^k, and thus if
the latter fits in a limb, so does q. Therefore I changed quite a bit
of the code to use single limbs instead of mpz_t's. I think this will
be OK, but if Robert could also verify this, I would be more
comfortable.
RE:
OK, that raises my limits.
There is some misunderstanding in your post:
If p is prime then I'm using exactly c primes of the form q=2*k*p+1,
where c is the largest integer for that p^c<=bound_for_p,
so for example if p>sqrt(bound_for_p),
then I'm using only one such prime, c=1 and in fact q can be bigger
than bound_for_p, the good news is that I think: for the first q prime
in this form:
q<2*p*log(p)^2 if p>2. It is similar to a conjecture in prime gaps.
Anyhow, we now need the following functions:
* b = powm(a, e, m) which sets b = a^e mod m, where a, b, e and m are
all mp_limb_t's
(we should pull a function out of FLINT long_extras module for this).
* mpn_root(r, ap, an, p) returns 1 is {ap, an} is a p-th power (else
0), setting r to the root if it is
(actually we don't need the root itself)
Currently we have mpz_root which basically uses mpn_rootrem which
computes the root and
remainder and checks if the remainder is 0. This is probably
sufficient for this code which does
O(1) root computations. (In general one can do much better than
mpz_root to test if something
is a root. One only needs to work with an/p bits and compute a root.
Then check the p-th power
of the bottom limb of this root gives the correct bottom limb of
{ap, an}. If not it is not a root. Of
course then we need to perform a complete root with remainder check).
* mpn_remove_1(bp, ap, an, p) remove the highest possible power of p
(which is 1 limb) from
{ap, an} and return the exponent.
Currently we have mpz_remove, but nothing optimised for the case
when p is 1 limb.
* multimod_reduce_ui(mp_limb_t * res, ap, an, mp_limb_t * primes,
count)
Given an array of count primes, each of which fits in a limb, reduce
{ap, an} by each of the primes
in the array and return the residues in the array of limbs res.
Robert did this with a tree, but it can be done slightly more
efficiently. It turns out that you can do
the bottom layer or two of the tree using umul_ppmm when multiplying
up, and then udiv_qrnnd for
the last layer when doing the reductions. This gives a substantial
speedup in practice.
It is important that this function allows aliasing between primes
and res to save space.
Some optimisations for later might include:
* Factor out the prime sieve. MPIR could have a globally defined array
of primes, computed up to some
limit, which could be extended as needed. It seems silly to compute
such a list of primes every time
mpn_perfect_power_p is called. We need to be careful about doing
this however, as global objects are
not threadsafe.
* Have an optimised version of mpn_perfect_power_p which takes a
precomputed "tree", then make
mpn_perfect_power_p just a wrapper for this function which keeps
track of some global precomputed tree
so that if someone calls mpn_perfect_power_p many times, the tree
does not need to be computed
multiple times. Again this needs to be done carefully for thread
safety reasons. This was handled in
FLINT with thread local storage, but that is not available on Sun or
Cygwin and it also means the
precomputation needs to be done per thread, which is unsatisfactory.
* It's not clear to me that we need to take the "tree" out as far as
we do. We should probably check at
each step that we haven't already exceeded the size of the original
n when multiplying up. The reason
is that reduction is trivial once our modulus is bigger than n
itself. But this is an optimisation that can be
used in the multimodular reduction code.
RE:
After the multiplication I think in the root of the tree: tree[1]
is smaller than n for all inputs. If you just mulitple all primes
up to log(n)/10 then that's about n^(1/10), we are using more than
one primes for small p, but in general only one, and for avarage case
q=2kp+1 is about log(p)*p. I've tried some values, and for all of
them
prod[1]<n was true.
* An additional copy of the input is made (n --> u) if no shifting is
required, to avoid destroying the input
to this function. This copy might be avoided somehow. It's probably
not a big deal, as it only happens
one time in 64 (or 32 on a 32 bit machine), and it is certainly not
where the complexity in the algorithm
lies, so I left it for now.
Bill.
--~--~---------~--~----~------------~-------~--~----~
You received this message because you are subscribed to the Google Groups
"mpir-devel" group.
To post to this group, send email to mpir-devel@googlegroups.com
To unsubscribe from this group, send email to
mpir-devel+unsubscr...@googlegroups.com
For more options, visit this group at
http://groups.google.com/group/mpir-devel?hl=en
-~----------~----~----~----~------~----~------~--~---
