[squeak-dev] Integer>>isPrime why reverting Knuth's algorithm P (probabilistic primality test) with iteration over division with even numbers?

David T. Lewis lewis at mail.msen.com
Sat Jan 23 20:51:54 UTC 2010

```To follow up on this discussion from last month, I updated Squeak trunk
such that LargePositiveInteger uses the probabilistic algorithm for
#isPrime, and added method comments to explain. A couple of folks suggested
this change, and Enrico concurred.

It turns out that SmallInteger maxVal is a reasonable point at which
to switch from use of #isPrime to #isProbablyPrime. On my system before
the change:

count := 1000.
largeMin := SmallInteger maxVal + 1.

"SmallInteger values up to maxVal"
Time millisecondsToRun: [(SmallInteger maxVal - count to: SmallInteger maxVal) do: [:e | e isPrime]]. ==> 120
Time millisecondsToRun: [(SmallInteger maxVal - count to: SmallInteger maxVal) do: [:e | e isProbablyPrime]]. ==> 984

"LargePositiveInteger values just above SmallInteger maxVal"
Time millisecondsToRun: [(largeMin to: largeMin + count) do: [:e | e isPrime]]. ==> 6599
Time millisecondsToRun: [(largeMin to: largeMin + count) do: [:e | e isProbablyPrime]]. ==> 714

After changing LargePositiveInteger>>isPrime, we have:

"LargePositiveInteger values just above SmallInteger maxVal"
Time millisecondsToRun: [(largeMin to: largeMin + count) do: [:e | e isPrime]]. ==> 719

So the implementation of LargePositiveInteger>>isPrime is now this:

isPrime
"Answer true if the receiver is a prime number. Use a probabilistic implementation	 that
is much faster for large integers, and that is correct to an extremely high statistical
level of confidence (effectively deterministic)."

^ self isProbablyPrime

Thanks to Enrico for his patience ;-)

Dave

On Sat Dec 19 13:58:16 UTC 2009, Enrico Spinielli wrote:
>
> The original implementation is was has been renamed by ul to isProbablyPrime
> Note that the probability is according to Knuth's words
>
> ...less than (1/4)^25 that such a 25-time-in-a-row procedure gives the wrong
> information about its input. This is less than one chance in a quadrillion;
> [...]
> It's much more likely that our computer has dropped a bit in its
> calculations,
> due to hardware malfunctions or cosmic radiation, than that Algorithm P has
> repeatedly guessed wrong!
>
> So 'probabilistic' in this case is very much deterministic!
> For accessible rationale about the algoritm (and a non probabilistic better
> one as well)
> you can also see "1.2.6 Example: Testing for
> Primality<http://mitpress.mit.edu/sicp/full-text/book/book-Z-H-11.html#%_sec_1.2.6>"
> in SICP.
> A possible improvement could have been to keep a list of the first N
> prime numbers (i.e. N=1000 or whatever integrer where gain in performance
> and or memory) and resort to the probabilistic test if self
> is greater than the bigger in the list.
>
> The value of original algorithm is all research and reasoning that went into
> it from
> Knuth et al. (note the order of growth is log(n), where n is the integer
> under scrutiny)
> The problem with the new implementation is that it goes thru testing numbers
> which are
> clearly not prime, 5, 15, 25, 35...just to mention multiples of 5.
> It can possibly be faster for small integers hence the above possible
> improvement suggestion...but for the rest it should be thrown away IMHO.
>
> Primality testing is very important for modern electronic communications,
> i.e. encryption
> and as such it has to be reliable and performant for big integers.
>
> Hope this clarifies
> Bye
> Enrico
> PS: the comment in the code was explicit enough to allow to research for
>       the rationale about the algorithm...we should not fix what works
> (well)
>       there are so many other fixes waiting...
> On Sat, Dec 19, 2009 at 12:56 AM, David T. Lewis <lewis at mail.msen.com>wrote:
>
> > On Fri, Dec 18, 2009 at 05:25:45PM +0100, Enrico Spinielli wrote:
> > > Hi all,
> > > I am back to checking Squeak after quite a while and got latest trunk.
> > > I looked after one of my contributions Integer>>isPrime
> > > and I found my implementation of Algorithm P from Knuth's AOCP vol 2
> > > substituted by an iteration of dividing self by all even numbers starting
> > > from 3
> > > and (correctly) stopping at self sqrtFloor.
> > > This is IMHO a questionable/useless "improvement", not even looking to
> > try
> > > to implement the Sieve of Eratostene...!
> > >
> > > Again IMHO isPrime should be reverted back to what has been renamed
> > > isProbablyPrime
> > >
> > > Not being anymore used to contribute I just signal it here...
> > >
> > > Hope it helps
> > > Bye
> > > --
> > > Enrico Spinielli
> >
> > Enrico,
> >
> > Is this your original implementation?
> >
> >
> >   isPrime
> >        "See isProbablyPrimeWithK:andQ: for the algoritm description."
> >        | k q |
> >        self <= 1 ifTrue: [^self error: 'operation undefined'].
> >        self even ifTrue: [^self = 2].
> >        k := 1.
> >
> >        q := self - 1 bitShift: -1.
> >        [q odd] whileFalse:
> >                        [q := q bitShift: -1.
> >                        k := k + 1].
> >
> >        25 timesRepeat: [(self isProbablyPrimeWithK: k andQ: q) ifFalse:
> > [^false]].
> >        ^true
> >
> > It was recently changed as follows:
> >
> > > Name: Kernel-ul.305
> > > Author: ul
> > > Time: 25 November 2009, 2:55:43.339 am
> > > Ancestors: Kernel-nice.304
> > >
> > > - added Integer >> #sqrtFloor, which returns the floor of the square root
> > > - renamed Integer >> #isPrime to #isProbablyPrime.
> > > - added Integer >> #isPrime which is implemented as a deterministic
> > primality test
> > > - both #isPrime and #isProbablyPrime return false for receivers <= 1
> > instead of raising an error
> >
> > Is this a reasonable change?
> >
> > Dave
> >
> >

```