Prime tests / BPSW tests in various number theory packages

[MattcAnderson]

in the Maple isprime() command, they use 100#, that is the product of the first 100 primes, as a precheck. Then
they come back and use 1000#.

Just for our edification,
100# is calculated, by me, in Maple, as

4711930799906184953162487834760260422020574773409675520188634839616415335845034221205289256705544681972439104097777157991804380284218315038719444943990492579030720635990538452312528339864352999310398481791730017201031090

Now that number is big. How many digits does it have?

Regards,
Matt

MODERATOR NOTE: Please put [code] tags around long numbers.

[paulunderwood]

Originally Posted by MattcAnderson

in the Maple isprime() command

According to https://www.maplesoft.com/support/he...x?path=isprime (sic) is a probablistic test. I wonder if it is like Pari/GP's 1+2 Selfridges BPSW variation, called ispseudoprime()?

[paulunderwood]

Originally Posted by paulunderwood

According to https://www.maplesoft.com/support/he...x?path=isprime (sic) is a probablistic test. I wonder if it is like Pari/GP's 1+2 Selfridges BPSW variation, called ispseudoprime()?

It would be interesting to compare timings for a 10k digit prime.

[Dr Sardonicus]

Originally Posted by MattcAnderson

in the Maple isprime() command, they use 100#, that is the product of the first 100 primes, as a precheck.
<snip>

Unfortunately, there are two conflicting notations about such products. If p is a prime number, p# ("p primorial" or "primorial p") often means the product of the primes up to p. Sometimes p_n# is used to mean the product of the first n primes.

BTW the product of the primes up to 541 (the first hundred primes) has 220 decimal digits. Generally, the product of the primes up to p has roughly p/log(10) decimal digits, where log() is the natural log.

Meaningless curiosity: 541 is the 100^th prime, and its index, 100, is one of the square summands in the expression of 541 as the sum of two squares, 541 = 100 + 441.

Stop the presses!

[Dr Sardonicus]

Originally Posted by retina

What about the poor forgotten real numbers? Are they of lesser value?

And ℚ? ℂ?

The real numbers are the completion of the rational numbers WRT the usual (Archimedean) valuation. The p-adic valuations are non-Archimedean.

No, the real and complex numbers are definitely not forgotten. In algebraic number theory, the valuations, be they Archimedean or not, are all called "primes." The Archimedean valuations are called "infinite primes."

See "Ostrowki's Theorem" and "product formula."

[R. Gerbicz]

Originally Posted by paulunderwood

According to https://www.maplesoft.com/support/he...x?path=isprime (sic) is a probablistic test. I wonder if it is like Pari/GP's 1+2 Selfridges BPSW variation, called ispseudoprime()?

On the given link: "It returns false if n is shown to be composite within one strong pseudo-primality test and one Lucas test. It returns true otherwise."

Likely they are doing some trial divisions also, and in that case this is a single Baillie–PSW primality test (https://en.wikipedia.org/wiki/Bailli...primality_test), and there is no known counter-example.

Originally Posted by paulunderwood

It would be interesting to compare timings for a 10k digit prime.

Do you know that Maple, and Mathematica is using gmp?
https://www.maplesoft.com/support/he....aspx?path=GMP
https://www.wolfram.com/legal/third-...enses/gmp.html
so for big numbers you won't get a faster test using these costly softwares.

[VBCurtis]

Originally Posted by MattcAnderson

I simply stated that small primes are not
unimportant.
......

This data is rare because if you want to see it,
you have to calculate it yourself. It is not
yet shared.

Importance and "rare"ness concern how difficult a finding is to reproduce. If I can calculate something faster than I can download a file and check it for viruses, the file is useless (but the code to make the file may be useful).

Lists of small primes are exactly this problem- they can be reproduced nearly instantly, so there is no point in publishing lists of them. If you want such lists, by all means collect them at home- but the files take up server space while providing use to exactly nobody (because anyone can generate them themselves rather than trusting someone else's code and / or virus-free-ness).

[paulunderwood]

Originally Posted by R. Gerbicz

On the given link: "It returns false if n is shown to be composite within one strong pseudo-primality test and one Lucas test. It returns true otherwise."

Likely they are doing some trial divisions also, and in that case this is a single Baillie–PSW primality test (https://en.wikipedia.org/wiki/Bailli...primality_test), and there is no known counter-example.


Do you know that Maple, and Mathematica is using gmp?
https://www.maplesoft.com/support/he....aspx?path=GMP
https://www.wolfram.com/legal/third-...enses/gmp.html
so for big numbers you won't get a faster test using these costly softwares.

I think GMP does a 1+3 Selfridges BPSW with \(Q\ne \pm 1\) for its Lucas part, whereas Pari/GP uses \(Q=1\) which makes BPSW 1+2 Selfridges. I wonder what Maple (or Mathematica) does.

I find GMP's documentation of BPSW very confusing and I am trying to find its source code to see the implementation. Perhaps n-rounds of Miller-Rabin followed by the (2.X Selfridges) BPSW variant of mine: Mod(Mod(2*x),x^2-P*x+1)^((n+1)/2) ==2*kronecker(2*(P+2),n) with kronecker(P^2-4,n)==-1, P minimal, checked to 2^64 would be better than what GMP has done -- or even just a simple Lucas chain with Q=1.

I suppose a few extra milliseconds is no big deal when you are trying to eliminate pseudoprimes.

[R. Gerbicz]

For Mathematica:
"PrimeQ first tests for divisibility using small primes, then uses the Miller–Rabin strong pseudoprime test base 2 and base 3, and then uses a Lucas test."
ref. https://reference.wolfram.com/langua...tion.html#6849

[MattcAnderson]

Thank you everybody for all the interest and inputs.

For those who are interested, get a load of this.

Here is internal code for Maple's isprime() command. I use showstat().

> showstat(isprime);
%;

isprime := proc(n)
local btor, nr, p, r;
   1   if not type(n,'integer') then
   2     if type(n,('complex')('numeric')) then
   3       error "argument must be an integer"
         else
   4       return 'isprime(n)'
         end if
       end if;
   5   if n < 2 then
   6     return false
       elif member(n,isprime:-special_primes) then
   7     return true
       elif igcd(2305567963945518424753102147331756070,n) <> 1 then
   8     return false
       elif n < 10201 then
   9     return true
       elif igcd(8496969489233418110532339909187349965926062586648932736611545426342203893270769390909069477309509137509786917118668028861499333825097682386722983737962963066757674131126736578936440788157186969893730633113066478620448624949257324022627395437363639038752608166758661255956834630697220447512298848222228550062683786342519960225996301315945644470064720696621750477244528915927867113,n) <> 1 then
  10     return false
       elif n < 1018081 then
  11     return true
       else
  12     r := gmp_isprime(n);
  13     if not r or n <= 5000000000 then
  14       return r
         end if;
  15     nr := igcd(408410100000,n-1);
  16     nr := igcd(nr^5,n-1);
  17     r := iquo(n-1,nr);
  18     btor := modp(('power')(2,r),n);
  19     if cyclotest(n,btor,2,r) = false or irem(nr,3) = 0 and cyclotest(n,btor,3,r) = false or irem(nr,5) = 0 and cyclotest(n,btor,5,r) = false or irem(nr,7) = 0 and cyclotest(n,btor,7,r) = false then
  20       return false
         end if;
  21     if isqrt(n)^2 = n then
  22       return false
         end if;
  23     for p from 3 while numtheory:-jacobi(p^2-4,n) <> -1 do
  24       NULL
         end do;
  25     return evalb(TraceModQF(p,n+1,n) = [2, p])
       end if
end proc

So 25 lines of Maple internal code. I don't understand all of it. But, keep reading.

Now show isprime(n) in action.
> isprime(7);

true
> isprime(8);

false

If you have read this far, then you get a grade of A+.
Here is a summary of what isprime() does.

1. check for positive integer
2. If the number is less than 100, compare to a list.
3. check for divisors less than 100.
(this uses the product of prime numbers less than 100 and integer greatest common divisor command (igcd)
4. check for divisors between 101 and 1,000.
(again a hard coded, large, minimal, 1000-smooth number, and igcd().)
oops, I know I stated this differently in a previous post, but I'm pretty sure I got it right this time.

5. More that one probable prime test.
These separate tests probably don't overlap.
It has been shown that if there are any counterexamples,
that is, cases when the isprime function fails,
they are > 10^100.
You can take my word on that.
There are academic papers where people have tested isprime(n) to some limit.
They probably tested specially chosen n values, looking for a failure.
If a number is found that causes isprime to fail, the software team at Maple will probably change their code, if only to patch the specific values that might, someday, be found.

The GMP library is apparently used by Maple.

[paulunderwood]

Thanks for that Matt. So it seems:

• Maple; 1+2 Selfridges
• Mathematica: 1+1+2 Selfridges
• Pari/GP: 1+2 Selfridges
• GMP: n*1+3 Selfridges (AFAICT)

What I'd like to see is n*1+m*2 Selfridges, although the chance of failing any of the above implementations is highly remote.