gmp prime testing

Jason Moxham
Sat, 2 Nov 2002 02:00:12 +0000

On Sunday 27 Oct 2002 9:52 pm, Kevin Ryde wrote:
> Jason Moxham <> writes:
> > You claim gmp is not really a number theory library , well looking at=
> > number of number theory functions implemented , yes I would have to a=
> > , but except for the lack of functions , gmp is perfect for number th=
> > , 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=
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 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(=
> > 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=
> > product can be done mod N , and this would be best if we expect divis=
> > 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=
> > rqft tests + 1 mulmod includes the sprp test , and it can be arranged=
> > that this part is done first.
> That still takes 3x longer than an sprp done alone though?

eg consider case N=3D1 mod 4 then the rqf test consists of a couple of ch=
, the order that you do the checks in doesn't matter , so you can arrange=
the "sprp looking" part is done first, therefore getting the earlier abor=
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=
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,=
> > 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=
testing , there are slower and less accurate than rqft , but they are fai=
useful . But if everyone uses slightly different definitions then they lo=
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=
in prime testing , also the rqft test includes all the underlying code=20

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


> > 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.


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