P+1 factorization method

The p+1 factorization method of factoring integers was discovered by Hugh Williams in 1982 and it is based in the p-1 method.

Instead of using multiplications modulo p, this method uses Lucas sequences. They have very interesting properties that make them useful in the Lucas-Lehmer test.

Let's define the following Lucas sequence: ${\displaystyle U_0=0,U_1=1,V_0=2,V_1=u}$ ${\displaystyle U_n=u{U}_{n-1}-U_{n-2},V_n=u{V}_{n-1}-V_{n-2}}$ for ${\displaystyle n\ge 2}$ where we can choose freely the integer ${\displaystyle u}$.

If we define ${\displaystyle D=u^2-4}$, then for any odd prime ${\displaystyle p}$, ${\displaystyle p}$ divides both gcd(${\displaystyle N}$, ${\displaystyle U_M}$) and gcd(${\displaystyle N}$, ${\displaystyle V_M-2}$) whenever ${\displaystyle M}$ is a multiple of ${\displaystyle p-(D/p)}$ where ${\displaystyle (D/p)}$ is the Legendre symbol, which is equal to 1 if ${\displaystyle D}$ is a quadratic residue mod ${\displaystyle p}$ and -1 otherwise (we suppose that ${\displaystyle D}$ is not multiple of ${\displaystyle p}$).

Since we do not know in advance the value of the Legendre symbol (because in order to compute it we need the prime factor ${\displaystyle p}$ of ${\displaystyle N}$ which is not known), it is possible that ${\displaystyle (D/p)=1}$ making this method a slow variant of p-1. So when the gcd does not reveal the prime factor of ${\displaystyle N}$ we retry with a different value of ${\displaystyle u}$ at least three times in order to have a high probability of selecting ${\displaystyle (D/p)=-1}$.

From the formulas above, the values U_n are not needed to compute the value of V_n, so they are not used in this algorithm.

Step 1

The step 1 works in the same way as explained in the P-1 factorization method except that instead of computing aⁿ (mod N) we compute V_n (mod N) and instead of computing gcd(aⁿ - 1, N) we compute gcd(V_n - 2, N).

The main computation in step 1 is to find V_mn from V_n (index multiplication). In order to compute the Lucas sequence for high values of the index, the formulas given at the beginning of this article are not useful. Instead, we will use the two formulas shown below: ${\displaystyle V_{2n}=V_n^2-2}$ (duplication formula) ${\displaystyle V_{m+n}=V_m{V}_n-V_{m-n}}$ (addition formula)

In order to perform index multiplication by a natural number ${\displaystyle M}$, there are two approaches using the formulas shown above.

Montgomery's PRAC is an almost optimal algorithm, but this is too complicated to explain it here.

Another algorithm is the same used in exponentiation. We represent the value ${\displaystyle M}$ in binary form. Then erase the first "1" (all numbers when represented in binary start with the digit "1"). Then for every "1" write the string "DA" and for every "0" write "D".

For example for the number 21 which is 10101 when written in binary, you will have the string DDADDA. The letter D means duplication and A means addition.

So if we have to compute ${\displaystyle V_{21n}}$ using the recipe above we will get the sequence: ${\displaystyle V_{2n}}$, ${\displaystyle V_{4n}}$, ${\displaystyle V_{5n}}$, ${\displaystyle V_{10n}}$, ${\displaystyle V_{20n}}$, ${\displaystyle V_{21n}}$.

A problem with this approach is that the formula for addition shown above for ${\displaystyle V_{q+r}}$ needs the value of ${\displaystyle V_{q-r}}$. This can be overcomed by using another variable. When we compute ${\displaystyle V_{kn}}$ we will be also computing ${\displaystyle V_{(k+1)n}}$.

Then if we start with the pair (${\displaystyle V_{kn}}$, ${\displaystyle V_{(k+1)n}}$), when the bit equals zero we compute (${\displaystyle V_{2kn}}$, ${\displaystyle V_{(2k+1)n}}$) and when the bit equals one we compute (${\displaystyle V_{(2k+1)n}}$, ${\displaystyle V_{2(k+1)n}}$). In both cases we need one duplication and one addition requiring only two modular multiplications.

The index multiplication explained above can be translated to pseudocode as follows (let B = V_n):

 x=B
 y=(B^2-2) mod N
  for each bit of M to the right of the most significant bit
    if the bit is 1
      x=(x*y-B) mod N
      y=(y^2-2) mod N
    else
      y=(x*y-B) mod N
      x=(x^2-2) mod N
  V=x

The PRAC algorithm is 30% faster in average than the algorithm shown above to multiply a point by a natural number.

Step 2

The step 2 is useful when p+1 has a prime factor ${\displaystyle q}$ between the bounds B₁ and B₂ and all the other prime factors less than B₁.

Using the formula ${\displaystyle V_n=u{V}_{n-1}-V_{n-2}}$ we find that exchanging ${\displaystyle V_n}$ and ${\displaystyle V_{n+2}}$ and adding two to ${\displaystyle n}$ we get: ${\displaystyle V_n=u{V}_{n+1}-V_{n+2}}$. From this we get: ${\displaystyle V_{-1}=u=V_1}$ and then ${\displaystyle V_{-n}=V_n}$.

Let ${\displaystyle V_r}$ be the number computed at the end of step 1. Now let ${\displaystyle C}$ be the maximum multiple of 6 less than or equal to B₁. The idea is to compute ${\displaystyle V_{Cr}}$ and ${\displaystyle V_{6r}}$. Then using the addition formula we can compute ${\displaystyle V_{(C+6)r},V_{(C+12)r},V_{(C+18)r},...,V_{B_{2}r}}$.

If ${\displaystyle q=C+1+6k}$ the number computed will be equal to ${\displaystyle V_{-r}}$ and if ${\displaystyle q=C-1+6k}$ the number computed will be equal to ${\displaystyle V_r}$. Since ${\displaystyle V_r=V_{-r}}$, we just multiply all ${\displaystyle V_{(C+6k)r}-V_r(\mathrm{mod}N)}$ and then take the GCD of the product in order to reveal the factor.

Example

Let's factor N = 451889 with u = 6 and B1 = 10. All congruences below will be modulo N.

As in the example of the P-1 factorization method we get:

E = 2³ × 3² × 5 × 7

V₁ = u = 6.

• Duplicate index: ${\displaystyle 6^2-2\equiv 34}$
• Duplicate index: ${\displaystyle {34}^2-2\equiv 1154}$
• Duplicate index: ${\displaystyle {1154}^2-2\equiv 427936}$
• Multiply index by 3: ${\displaystyle {427936}^2-2\equiv 299066}$, so the first pair is ${\displaystyle (V_n,V_{2n})\equiv (427936,299066)}$. The second pair is ${\displaystyle (V_{3n},V_{4n})\equiv (V_n\ast V_{2n}-V_n,V_{2n}^2-2)\equiv (292372,342029)}$
• Multiply index by 3: ${\displaystyle {292372}^2-2\equiv 255586}$, so the first pair is ${\displaystyle (V_n,V_{2n})\equiv (292372,255586)}$. The second pair is ${\displaystyle (V_{3n},V_{4n})\equiv (V_n\ast V_{2n}-V_n,V_{2n}^2-2)\equiv (176913,33332)}$
• Multiply index by 5: ${\displaystyle {176913}^2-2\equiv 377427}$, so the first pair is ${\displaystyle (V_n,V_{2n})\equiv (176913,377427)}$. The second pair is ${\displaystyle (V_{2n},V_{3n})\equiv (V_n^2-2,V_n\ast V_{2n}-V_n)\equiv (377427,447298)}$. The third pair is ${\displaystyle (V_{5n},V_{6n})\equiv (V_{2n}\ast V_{3n}-V_n,V_{3n}^2-2)\equiv (50045,290385)}$.
• Multiply index by 7: ${\displaystyle {50045}^2-2\equiv 133185}$, so the first pair is ${\displaystyle (V_n,V_{2n})\equiv (50045,133185)}$. The second pair is ${\displaystyle (V_{3n},V_{4n})\equiv (V_n\ast V_{2n}-V_n,V_{2n}^2-2)\equiv (282419,245306)}$. The third pair is ${\displaystyle (V_{7n},V_{8n})\equiv (V_{3n}\ast V_{4n}-V_n,V_{4n}^2-2)\equiv (374468,138727)}$.

So gcd(374468 - 2, 451889) = 139 which is a proper factor of the original number.

Notice that some values computed above are not really useful if you had to perform the computation by hand, for example the second element of the pair in the last step of the index multiplication. When using a computer, these additional multiplication can be avoided by using conditional sentences.

If instead of u=6 we would have selected u=7, the final gcd would have been: gcd(252303-2, 451889) = 1, so this particular value u=7 is not useful in the p+1 method.

Since u=7 does not factor 451889 in step 1, we have to continue with step 2. Let the upper bound B₂ be 50. From the previous paragraph ${\displaystyle V_r\equiv 252303}$.

Initializing:

${\displaystyle C=18}$

${\displaystyle V_{18r}\equiv 392380}$

${\displaystyle V_{24r}\equiv 224225}$

${\displaystyle V_{6r}\equiv 446313}$

Computing all other ${\displaystyle V_{6kr}}$:

${\displaystyle V_{30r}=V_{24r}{V}_{6r}-V_{18r}\equiv 157772}$

${\displaystyle V_{36r}=V_{30r}{V}_{6r}-V_{24r}\equiv 318875}$

${\displaystyle V_{42r}=V_{36r}{V}_{6r}-V_{30r}\equiv 430332}$

${\displaystyle V_{48r}=V_{42r}{V}_{6r}-V_{36r}\equiv 132372}$

${\displaystyle (V_{18r}-V_r)(V_{24r}-V_r)...(V_{48r}-V_r)\equiv 421170}$

${\displaystyle gcd(421170,451889)=139}$

External links

• Williams's p+1 algorithm