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: [math]\displaystyle{ \large U_0 = 0\,,\, U_1 = 1\,,\, V_0 = 2\,,\, V_1 = u }[/math] [math]\displaystyle{ \large U_n = uU_{n-1} - U_{n-2}\,,\, V_n = uV_{n-1} - V_{n-2} }[/math] for [math]\displaystyle{ n\geq 2 }[/math] where we can choose freely the integer [math]\displaystyle{ u }[/math].

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

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

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: [math]\displaystyle{ \large V_{2n} = V_n^2 - 2 }[/math] (duplication formula) [math]\displaystyle{ \large V_{m+n} = V_m V_n - V_{m-n} }[/math] (addition formula)

In order to perform index multiplication by a natural number [math]\displaystyle{ M }[/math], 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 [math]\displaystyle{ M }[/math] 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 [math]\displaystyle{ V_{21n} }[/math] using the recipe above we will get the sequence: [math]\displaystyle{ V_{2n} }[/math], [math]\displaystyle{ V_{4n} }[/math], [math]\displaystyle{ V_{5n} }[/math], [math]\displaystyle{ V_{10n} }[/math], [math]\displaystyle{ V_{20n} }[/math], [math]\displaystyle{ V_{21n} }[/math].

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

Then if we start with the pair ([math]\displaystyle{ V_{kn} }[/math], [math]\displaystyle{ V_{(k+1)n} }[/math]), when the bit equals zero we compute ([math]\displaystyle{ V_{2kn} }[/math], [math]\displaystyle{ V_{(2k+1)n} }[/math]) and when the bit equals one we compute ([math]\displaystyle{ V_{(2k+1)n} }[/math], [math]\displaystyle{ V_{2(k+1)n} }[/math]). 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 [math]\displaystyle{ q }[/math] between the bounds B₁ and B₂ and all the other prime factors less than B₁.

Using the formula [math]\displaystyle{ V_n = uV_{n-1} - V_{n-2} }[/math] we find that exchanging [math]\displaystyle{ V_n }[/math] and [math]\displaystyle{ V_{n+2} }[/math] and adding two to [math]\displaystyle{ n }[/math] we get: [math]\displaystyle{ V_n = uV_{n+1} - V_{n+2} }[/math]. From this we get: [math]\displaystyle{ V_{-1} = u = V_1 }[/math] and then [math]\displaystyle{ V_{-n} = V_n }[/math].

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

If [math]\displaystyle{ q=C+1+6k }[/math] the number computed will be equal to [math]\displaystyle{ V_{-r} }[/math] and if [math]\displaystyle{ q=C-1+6k }[/math] the number computed will be equal to [math]\displaystyle{ V_r }[/math]. Since [math]\displaystyle{ V_r = V_{-r} }[/math], we just multiply all [math]\displaystyle{ V_{(C+6k)r} - V_r\,\pmod N }[/math] 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: [math]\displaystyle{ 6^2 - 2\equiv 34 }[/math]
• Duplicate index: [math]\displaystyle{ 34^2 - 2\equiv 1154 }[/math]
• Duplicate index: [math]\displaystyle{ 1154^2 - 2\equiv 427936 }[/math]
• Multiply index by 3: [math]\displaystyle{ 427936^2 - 2\equiv 299066 }[/math], so the first pair is [math]\displaystyle{ (V_n , V_{2n})\equiv (427936,\, 299066) }[/math]. The second pair is [math]\displaystyle{ (V_{3n}, V_{4n})\equiv (V_n * V_{2n} - V_n,\, V_{2n}^2 - 2)\equiv (292372,\, 342029) }[/math]
• Multiply index by 3: [math]\displaystyle{ 292372^2 - 2\equiv 255586 }[/math], so the first pair is [math]\displaystyle{ (V_n , V_{2n})\equiv (292372,\, 255586) }[/math]. The second pair is [math]\displaystyle{ (V_{3n}, V_{4n})\equiv (V_n * V_{2n} - V_n,\, V_{2n}^2 - 2)\equiv (176913,\, 33332) }[/math]
• Multiply index by 5: [math]\displaystyle{ 176913^2 - 2\equiv 377427 }[/math], so the first pair is [math]\displaystyle{ (V_n , V_{2n})\equiv (176913,\, 377427) }[/math]. The second pair is [math]\displaystyle{ (V_{2n}, V_{3n})\equiv (V_n^2 - 2,\, V_n * V_{2n} - V_n)\equiv (377427,\, 447298) }[/math]. The third pair is [math]\displaystyle{ (V_{5n}, V_{6n})\equiv (V_{2n} * V_{3n} - V_n,\, V_{3n}^2 - 2)\equiv (50045,\, 290385) }[/math].
• Multiply index by 7: [math]\displaystyle{ 50045^2 - 2\equiv 133185 }[/math], so the first pair is [math]\displaystyle{ (V_n , V_{2n})\equiv (50045,\, 133185) }[/math]. The second pair is [math]\displaystyle{ (V_{3n}, V_{4n})\equiv (V_n * V_{2n} - V_n,\, V_{2n}^2 - 2)\equiv (282419,\, 245306) }[/math]. The third pair is [math]\displaystyle{ (V_{7n}, V_{8n})\equiv (V_{3n} * V_{4n} - V_n,\, V_{4n}^2 - 2)\equiv (374468,\, 138727) }[/math].

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 [math]\displaystyle{ V_r \equiv 252303 }[/math].

Initializing:

[math]\displaystyle{ C = 18 }[/math]

[math]\displaystyle{ V_{18r} \equiv 392380 }[/math]

[math]\displaystyle{ V_{24r} \equiv 224225 }[/math]

[math]\displaystyle{ V_{6r} \equiv 446313 }[/math]

Computing all other [math]\displaystyle{ V_{6kr} }[/math]:

[math]\displaystyle{ V_{30r} = V_{24r} V_{6r} - V_{18r} \equiv 157772 }[/math]

[math]\displaystyle{ V_{36r} = V_{30r} V_{6r} - V_{24r} \equiv 318875 }[/math]

[math]\displaystyle{ V_{42r} = V_{36r} V_{6r} - V_{30r} \equiv 430332 }[/math]

[math]\displaystyle{ V_{48r} = V_{42r} V_{6r} - V_{36r} \equiv 132372 }[/math]

[math]\displaystyle{ (V_{18r} - V_r)(V_{24r} - V_r)...(V_{48r} - V_r) \equiv 421170 }[/math]

[math]\displaystyle{ \gcd(421170,\, 451889) = 139 }[/math]

External links

• Williams's p+1 algorithm