gmp prime testing

[Jason Moxham]

On Sunday 27 Oct 2002 9:52 pm, Kevin Ryde wrote:
> Jason Moxham <J.L.Moxham@maths.soton.ac.uk> writes:
> > You claim gmp is not really a number theory library , well looking at=
 the
> > number of number theory functions implemented , yes I would have to a=
gree
> > , but except for the lack of functions , gmp is perfect for number th=
eory
> > , fast,portable,accurate,free and in C.
>
> Remember other libraries like LiDIA and NTL exist on top of gmp, doing
> lots of number theory stuff, very likely more and better than we can.
>

they are C++ though , yuck

> > " but basic things well-implemented would suit" , to me solving
> > x^2+dy^2=3Dp in the integers is a pretty basic function ,
>
> Does that fall out of units in a field of quadratic integers?  I'm not
> well up with such stuff, but it might be the sort of thing I'd imagine
> a computer algebra system could do more flexibly than we can.
>

just an example of a simple , faily common equations , that are best solv=
ed=20
NOT using an algebra system=20


> Apropos of that sort of equation, I'd thought at one time of a demo
> program doing the cattle problem of archimedes, showing how a PC with
> gmp can do in less than a second what was reported in 1965 in Math
> Comp as taking nearly 8 hours :-).  A demo program for gmpxx.h is
> wanted too, a cattle.cc could kill two birds with one stone.
>
> > if we do 3 reps of random sprp then 2 reps of rqft do we get
> > 1/64 * 1/7710^2 or just the highest ?
>
> No idea.  You'd like to hope yes, but there's probably some heavy
> theoretical questions there.
>

I will see if I can find out

> > for the trial division before a sprp test on N, it seems that k*n*lg(=
n)^5
> > is a good bound(to minimize total time for random inputs) where n=3Dl=
g(N) ,
> > calculated with a bit of theory and some simulated runs.
>
> Right now I guess we're on (log2(N)/30)^2, following Knuth maybe, but
> with some possibly premature rounding.
>
> > A better trial div for large N might be to multiply all the divisors
> > together (using binary splitting) and take the gcd with N , obviously=
 the
> > product can be done mod N , and this would be best if we expect divis=
ors
> > to be rare.
>
> If you mean an Nx1 gcd, yes that might be an option since it can use
> the faster modexact_1, though the 1x1 gcd at the end will be slower
> than 1x1 remainders (by the primes in the product).
>
> If you mean an NxM gcd, remember gcd becomes very slow as M increases.
> Might still win, I don't know.
>

I was thinking of both cases , some experimental code may have to be writ=
en to=20
see if we can get a speedup

> > Also doing the single sprp before the rqft test is not neccesairy, th=
e
> > rqft tests + 1 mulmod includes the sprp test , and it can be arranged=
 so
> > that this part is done first.
>
> That still takes 3x longer than an sprp done alone though?
>

No=20
eg consider case N=3D1 mod 4 then the rqf test consists of a couple of ch=
ecks
, the order that you do the checks in doesn't matter , so you can arrange=
 that=20
the "sprp looking" part is done first, therefore getting the earlier abor=
t=20
that a separate sprp test would give on a composite but without the 3x=20
runtime. One minor drawback is that the parameter must chosen such that=20
jacobi(b,N)=3D1 , so time to first abort is the same as sprp test but wit=
h a=20
jacobi symbol computation as well . As jacobi if order n^2 ( power is ord=
er=20
n^3) this isn't to significant , and in practice you can hardly tell.



> > int=09mpz_probable_prime_p(mpz_srcptr N,int reps,gmp_randstate_t rnd,=
int
> > expected_type,unsigned long trial_div_lim)
>
> Maybe.  I had in mind all the pieces broken out and available to be
> run separately.  Right now that might merely mean an mpz_millerrabin
> Torbjorn was thinking of exposing, and your mpz_trial_div proposal.
>
> > I have temporaliry dropped doing the lucas
> > probable prime tests for the reason that I guessed they would not get
> > accepted and also there is no consistent definition in the books !!!!
>
> The definition doesn't matter, only the power of the test :-).
>

The lucas prob prime tests I have dropped , were not ment for prob. prime=
=20
testing , there are slower and less accurate than rqft , but they are fai=
rly=20
useful . But if everyone uses slightly different definitions then they lo=
se=20
there appeal , we would have to code for each definition. Would still be=20
usefull to have a fn to calculate lucas fn's mod n , as these are very co=
mmon=20
in prime testing , also the rqft test includes all the underlying code=20
allready.

> > Currently gmp has prime tables to 1000 , which takes 168*2 bytes , pr=
ime
> > table to 65536 would take 13K , prime difference to 65536 would take =
6.4K
> > . if this is too much
>
> Don't know.  I'm hoping Torbjorn will make that decision.  :-)
>

yes

> > construct_and_use_prime_table(global_prime_table_pointer,limits);
>
> PARI does something like that, it's nice and flexible for
> applications.
>
> > you would call this at a "thread safe time"
>
> Yes, we don't want to get intimate with the thread manager.
>
> > ul=09mpz_trial_div(mpz_srcptr N , ul start , ul stop)
>
> Sounds fair to me, if Torbjorn likes it.
>
>

yes

> PS. Some or all of this discussion could probably move to gmp-devel@swo=
x.