Lucas-Lehmer test for Mersenne numbers (J)

From LiteratePrograms
Jump to: navigation, search
Other implementations: Erlang | Haskell | J | Java | Lisp | Python | Ruby | Scheme

[edit] theory

According to the wikipedia article, first we define a sequence as follows:


  s_i=
   \begin{cases}
    4 & \mbox{if }i=0;
   \\
    s_{i-1}^2-2 & \mbox{otherwise.}
   \end{cases}

Then a Mersenne number (Mp = 2p− 1) is prime iff

s_{p-2}\equiv0\pmod{M_p};

[edit] practice

The sequence recurrence is straightfoward to define: NB. we define as a dyad so that we may do the modulo arithmetic on every iteration.

seq =: dyad : 'x|_2+*:y'

and, sure enough, the first few terms of the sequence agree with OEIS A00301. NB. 4x is 4 represented as an extended (arbitrary-precision) integer.

        0 seq^:(i.8) 4x
4 14 194 37634 1416317954 2005956546822746114 4023861667741036022825635656102100994 16191462721115671781777559070120513664958590125499158514329308740975788034

Now we define a function, using the iteration operator (^:) to find the residue sp − 2(mod Mp):

res =: monad : '(_1+2x^y) seq^:(y-2) 4x'

[edit] wrapping up

We can double check our implementation against the examples in the article:

         res each 3 11
+-+----+
|0|1736|
+-+----+
        (0=res) each 3 11
+-+-+
|1|0|
+-+-+

Or we can test primality in a different manner by counting (#) factors (q:) of the first few Mersenne numbers (Mp = 2p− 1), then checking that two sets of indications (zero residue or single factor) agree:

        mpf =: monad : '#q:_1+2x^y'
        >((0=res)=(1=mpf)) each p:1+i.15
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

Alternatively, we could wrap everything up into the following single-line definition:

llmp =: monad : '0=(_1+2x^y) ([|_2+[:*:])^:(y-2) 4x'
Download code
hijacker
hijacker
hijacker
hijacker