Prime tests / BPSW tests in various number theory packages
 

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 [code]4711930799906184953162487834760260422020574773409675520188634839616415335845034221205289256705544681972439104097777157991804380284218315038719444943990492579030720635990538452312528339864352999310398481791730017201031090[/code]

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

Regards,
Matt

[b][color=red]MODERATOR NOTE:[/b] Please put [noparse][code][/noparse] tags around long numbers.[/color]

[QUOTE=MattcAnderson;600118]in the Maple isprime() command
[/QUOTE]

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

[QUOTE=paulunderwood;600121]According to [url]https://www.maplesoft.com/support/help/Maple/view.aspx?path=isprime[/url] (sic) is a probablistic test. I wonder if it is like Pari/GP's 1+2 Selfridges BPSW variation, called ispseudoprime()?[/QUOTE]

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

[QUOTE=MattcAnderson;600118]in the Maple isprime() command, they use 100#, that is the product of the first 100 primes, as a precheck.
<snip>[/QUOTE]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[sub]n[/sub]# 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.

[b][u]Meaningless curiosity[/u]:[/b] 541 is the 100[sup]th[/sup] 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! :rolleyes:

[QUOTE=retina;600116]What about the poor forgotten real numbers? Are they of lesser value? :sad:

[size=1]And ℚ? ℂ?[/size][/QUOTE]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 [i]all[/i] called "primes." The Archimedean valuations are called "infinite primes."

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

[QUOTE=paulunderwood;600121]According to [url]https://www.maplesoft.com/support/help/Maple/view.aspx?path=isprime[/url] (sic) is a probablistic test. I wonder if it is like Pari/GP's 1+2 Selfridges BPSW variation, called ispseudoprime()?[/QUOTE]

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 ([url]https://en.wikipedia.org/wiki/Baillie%E2%80%93PSW_primality_test[/url]), and there is no known counter-example.

[QUOTE=paulunderwood;600123]It would be interesting to compare timings for a 10k digit prime.[/QUOTE]

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

[QUOTE=MattcAnderson;600117]
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.
[/QUOTE]

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

[QUOTE=R. Gerbicz;600127]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 ([url]https://en.wikipedia.org/wiki/Baillie%E2%80%93PSW_primality_test[/url]), and there is no known counter-example.


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

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.

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. [url]https://reference.wolfram.com/language/tutorial/SomeNotesOnInternalImplementation.html#6849[/url]

isprime() guts
 

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



[CODE]> 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[/CODE]
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.

Thanks for that Matt. So it seems:
[LIST][*]Maple; 1+2 Selfridges[*]Mathematica: 1+1+2 Selfridges[*]Pari/GP: 1+2 Selfridges[*]GMP: n*1+3 Selfridges (AFAICT)[/LIST]
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.