Integer Factorization via Continued Fractions and Quadratic Forms
Nadir Murru Giulia Salvatori



§ ABSTRACT 
We propose a novel factorization algorithm that leverages the theory underlying the SQUFOF method, including reduced quadratic forms, infrastructural distance, and Gauss composition. We also present an analysis of our method, which has a computational complexity of O ( exp( 3/√(8)√(ln N lnln N)) ), making it more efficient than the classical SQUFOF and CFRAC algorithms. Additionally, our algorithm is polynomial-time, provided knowledge of a (not too large) multiple of the regulator of ℚ(N).


Quantum features of the transport through ion channels in the soft knock-on model
    @
    September 9, 2024
=================================================================================





§ INTRODUCTION


The integer factorization problem is a fascinating challenge in number theory, with many important theoretical aspects and practical applications (e.g., in cryptography, where the most important public key cryptosystems are based on the hardness of solving this problem for large composite numbers). Indeed, currently, there does not exist a polynomial algorithm for factorizing integers and thus the research in this field is fundamental and active. 

To date, the most efficient algorithm known for factoring integers larger than 10^100 is the General Number Field Sieve designed by Pomerance <cit.>, with a heuristic running time of exp (  (√(64/9) + o(1)  )(ln N)^1/3(lnln N)^2/3 ). 
However, for smaller numbers other algorithms perform better, such as the Quadratic Sieve and SQUare FOrm Factorization (SQUFOF).

SQUFOF algorithm (the best method for numbers between 10^10 and 10^18) was proposed by Shanks in <cit.> and it is based on the properties of square forms and continued fractions. 
The algorithm is discussed, for example, in <cit.> and <cit.>. A rigorous and complete description of the method and its complexity is provided in the well-regarded paper by Gower and Wagstaff <cit.>, where the details are meticulously presented and the algorithm is examined in depth.
Recently, the SQUFOF algorithm has been revisited by Elia <cit.> who proposed an improvement whose complexity is based on the computation of the regulator of a quadratic field. Another improvement, which exploits a sieve inspired by the Quadratic Sieve, can be found in <cit.>.

The other main algorithm, which exploits the theory of continued fractions, is CFRAC <cit.> which was implemented and used for factorizing large numbers (such as the 7th Fermat number) by Morrison and Brillhart <cit.>.

In this paper, we focus on the underlying theory of SQUFOF and, starting from the work by Elia <cit.>, we improve it, obtaining a novel factorization algorithm whose complexity for factoring the integer N is O  ( exp (3/√(8)√(ln N lnln N) )  ), making it more efficient than the classical CFRAC and SQUFOF algorithms.
The time complexity is similar to that obtained in <cit.>, although their method uses a different approach.


The paper is structured as follows. Section <ref> is devoted to introducing the notation and developing the auxiliary results utilized in the design of the factorization algorithm. Specifically, in Section <ref>, we deal with the theory of continued fractions, focusing on the expansion of square roots and the properties of particular sequences arising from these expansions. Section <ref> discusses when a nontrivial factor of N can be found in the case that the period of the continued fraction expansion of √(N) is even. Finally, in Section <ref>, we introduce the tools regarding quadratic forms, including the notion of distance between forms, the reduction operator, and the Gauss composition, focusing on the properties of particular sequences of quadratic forms used in the algorithm. In Section <ref>, we present and discuss all the details of our new algorithm and analyze the time complexity, highlighting also the fundamental role played by the computation of the regulator of ℚ(√(N)).



§ PRELIMINARIES AND AUXILIARY RESULTS
 


 §.§ Continued fractions
 
It is well-known that the continued fraction expansion of quadratic irrationals is periodic and in this case the Lagrange algorithm can be used for obtaining such expansion. Let us consider, without loss of generality, quadratic irrationals

    α_0 = P_0 + √(N)Q_0,

with  P_0, Q_0, N ∈ℤ, N > 0 not square, Q_0 ≠0 and Q_0 | N - P_0^2. The continued fractions expansion [a_0, a_1, …] of α_0 can be obtained computing

    a_m = ⌊α_m ⌋
    
    P_m+1 = a_m Q_m - P_m 
    
    Q_m+1 = (N-P_m+1^2)/Q_m,   where  α_m = P_m + √(N)/Q_m,    m ≥ 0.

We recall that the continued fraction expansion of √(N) is periodic and has the following particular form

    √(N) = [a_0, a_1, a_2, a_3, …, a_τ -1, 2a_0],

where the sequence (a_1, …, a_τ-1) is a palindrome. 



Kraitchik <cit.> showed that the period τ of the continued fraction expansion of √(N) is upper bounded by
0.72 √(N)ln N, for N > 7.
However, the period length has irregular behavior as a function of N: it may assume any value from 1, when N = M^2 + 1, to values greater than √(N)lnln N (see <cit.> and <cit.>, respectively).

From now on, we always consider the continued fraction expansion of √(N) as in (<ref>) (i.e., quadratic irrationals with Q_0 = 1 and P_0 = 0). Let {p_n}_n ≥ -1 and {q_n}_n ≥ -1 be the sequences of numerators and denominators of convergents of √(N), defined by

    p_-1 = 1,    p_0 = a_0,    q_-1 = 0,    q_0 = 1

and

    p_m = a_m p_m-1 + p_m-2,    q_m = a_m q_m-1 + q_m-2,   ∀  m ≥ 1.

We also recall the following two properties

    [  a_0 1  1    0 ]⋯[  a_m 1  1    0 ] = [       p_m p_m-1 q_m     q_m-1 ]  ∀  m ≥ 0,

and

    p_τ -2 = -a_0 p_τ -1 + Nq_τ -1
    
    q_τ -2 = p_τ -1 -a_0 q_τ -1,

see, e.g., <cit.>.

Next, we examine the sequence {𝔠_n  }_n ≥ -1, defined by

    𝔠_m := p_m + q_m √(N).

The result in the next proposition can also be found in <cit.>. Here, we provide a slightly different proof for completeness.


The sequence {𝔠_n  }_n ≥ -1 satisfies the relation


    𝔠_m + k τ = 𝔠_m𝔠_τ -1^k   for all  k ∈ℕ and m≥-1.

 



The claimed equality is trivial for k = 0.
First, we prove by induction on m the equality for k=1, and then we generalize for k>1. 

The case m=-1 is trivial, since 𝔠_-1 = 1.
Now, we proceed by induction. Using the inductive 
hypothesis, we consider the following chain of equalities

    𝔠_τ +m+1   = a_m+1𝔠_τ +m + 𝔠_τ +m-1 = a_m+1𝔠_m𝔠_τ -1 + 𝔠_m-𝔠_τ -1
       = (a_m+1𝔠_m + 𝔠_m-1 )𝔠_τ -1 = 𝔠_m+1𝔠_τ -1

which concludes the proof in the case k=1.

For the case k>1 we iterate as follows

    𝔠_m + kτ = 𝔠_m + (k-1)τ𝔠_τ - 1 = ⋯ = 𝔠_m𝔠_τ - 1^k.



We recall that the minimal positive solution of the Pell Equation X^2 - N Y^2 = 1 is (p_τ-1, q_τ -1) if τ is even, and (p_2τ-1, q_2τ -1) if τ is odd. 
We denote by R^+(N) the logarithm of the minimal positive solution of the Pell equation, by R(N) the regulator of ℚ(√(N)), and by 𝒩 the field norm of ℚ(√(N)). Moreover, given x + y √(N)∈ℚ(√(N)), we denote by x + y √(N) the element x - y √(N).



We have

    (-1)^m P_m+1 = p_mp_m-1 - Nq_mq_m-1  ∀  m ≥ 0

and

    (-1)^m+1 Q_m+1 =p_m^2 -Nq_m^2= 𝒩(𝔠_m)   ∀  m ≥ -1.



The proof is straightforward by induction.




Since we will exploit these sequences to factor the integer N, it is computationally important to bound their elements. 


We have

    0< Q_m < 2/a_m√(N) ,    0 ≤ P_m < √(N)  ∀  m ≥ 0.



See <cit.>.


The following lemma proves two equalities useful in Theorem <ref>.

It holds that

    p_τ -m-2 = (-1)^m-1p_τ -1p_m + (-1)^mNq_τ -1q_m
    
        q_τ -m-2 = (-1)^mp_τ -1q_m + (-1)^m-1q_τ -1p_m
    ∀  -1 ≤ m ≤τ -1.



Using Equation (<ref>) and the fact that (a_1, …, a_τ -1) is palindrome, we obtain, for all 0 ≤ m ≤τ - 2,

    [        p_τ -1 p_τ -2 q_τ -1        q_τ -2 ]   = [         p_τ-m-2 p_τ-m-3 q_τ-m-2         q_τ-m-3 ][ a_τ -m-1     1  1        0 ]⋯[ a_τ -1   1  1      0 ]
       = [         p_τ-m-2 p_τ-m-3 q_τ-m-2         q_τ-m-3 ][ a_m+1  1  1     0 ]⋯[  a_1 1  1    0 ]
       = [         p_τ-m-2 p_τ-m-3 q_τ-m-2         q_τ-m-3 ] [ [  a_0 1  1    0 ]^-1[     p_m+1 p_m q_m+1       q_m ] ]^T
       = [         p_τ-m-2 p_τ-m-3 q_τ-m-2         q_τ-m-3 ][     p_m+1 q_m+1 p_m       q_m ][    0 1  1 -a_0 ].

Multiplying by the inverse of the matrices [     p_m+1 q_m+1 p_m       q_m ] and [    0 1  1 -a_0 ], and using Equation (<ref>), we obtain

    [         p_τ-m-2 p_τ-m-3 q_τ-m-2         q_τ-m-3 ]   = (-1)^m[        p_τ -1 p_τ -2 q_τ -1        q_τ -2 ][  a_0 1  1    0 ][         q_m -q_m+1 -p_m       p_m+1 ]
       =(-1)^m [             Nq_τ-1q_m - p_mp_τ-1 -Nq_τ-1q_m+1 p_τ-1q_m - p_mq_τ-1           -p_τ-1q_m+1+q_τ-1p_m+1 ].

The cases m=-1 and m=τ -1 are straightforward to verify.

The transformation defined by (<ref>) is identified by the matrix


    M_τ -1 =
    [ -p_τ -1 Nq_τ -1; -q_τ -1  p_τ -1 ]
    .


The results that follow in this subsection are those found by Elia in <cit.>, but they are further extended, approaching also the case of odd periods, and we include more detailed proofs.
The sequences {Q_n}_n ≥ 0 and {P_n}_n ≥ 1 are periodic of period τ, where τ is the period of the sequence of partial quotients { a_n}_n ≥ 0 of the continued fraction expansion of √(N). Further, within a period, there exist interesting symmetries.


The sequence {Q_n}_n ≥ 0 is periodic with period τ. The elements of the first block { Q_n }_n=0^τ satisfy the
symmetry relation 
    Q_m = Q_τ -m,   ∀  0 ≤ m ≤τ.




Using Equation (<ref>) and the fact that 𝒩(p_τ -1 + q_τ -1√(N)) = (-1)^τ, the following chain of equalities holds

    Q_m+τ   =  |𝒩(𝔠_m-1+ τ)  | =  | 𝒩(𝔠_m-1𝔠_τ -1)  | =  | 𝒩(𝔠_m-1) 𝒩( 𝔠_τ -1)  | = Q_m,

from which we deduce that the period of {Q_n}_n ≥ 0 is τ.

The symmetry of the sequence {Q_n}_n ≥ 0 within the τ elements of the first period follows from Equation (<ref>).
We have

    p_τ -m-2^2 -Nq_τ -m-2^2    = (p_τ -1p_m -Nq_τ -1q_m)^2 -N(-p_τ -1q_m + q_τ -1p_m)^2 
       =  (p_m^2 - Nq_m^2)(p_τ -1^2 - Nq_τ -1^2) 
       =  (-1)^τ(p_m^2 - Nq_m^2),

implying that (-1)^τ -m-1Q_τ -m-1 = (-1)^τ(-1)^m+1 Q_m+1 for all -1 ≤ m ≤τ -1. 




The sequence {P_n}_n ≥ 1 is periodic with period τ. The elements of the first block { P_m }_m=1^τ -1 satisfy the symmetry relation

    P_τ -m-1 = P_m,   ∀  1 ≤ m ≤τ -2.




The periodicity of the sequence {P_n}_n ≥ 1 follows from the property expressed by Equation (<ref>) and Equation (<ref>), noting that

    (-1)^m P_m+1   = 1/2(𝔠_m𝔠_m-1 + 𝔠_m𝔠_m-1) 
       =1/2 (𝔠_m+τ/𝔠_τ - 1𝔠_m-1+τ/𝔠_τ - 1 + 𝔠_m + τ/𝔠_τ - 1𝔠_m-1+τ/𝔠_τ - 1 ) 
       = (-1)^τ (-1)^m+τ P_m+1+τ.

The next chain of equalities proves the symmetry property

    (-1)^τ -m-1P_τ -m   = p_τ -m-1p_(τ -1)-m-1 - Nq_τ-1-mq_(τ -1) -m-1
       = -(p_τ -1p_m - Nq_τ -1q_m)(p_τ -1p_m-1 - Nq_τ -1q_m-1) 
       + N(p_τ -1q_m - q_τ -1p_m)(p_τ -1q_m-1 - p_m -1q_τ-1) 
       = -(p_τ -1^2 -Nq_τ -1^2)(p_mp_m-1 -Nq_mq_m-1) 
       = -(-1)^τ(p_mp_m-1 -Nq_mq_m-1) 
       = (-1)^τ +1 (-1)^m P_m+1,

where, in the second-to-last equality, we used Equation (<ref>).




Note that M_τ -1^2 = (-1)^τI_2, with I_2 the identity matrix, and, if τ is even, the eigenvalues of M_τ -1 are λ_0 =1 and λ_1 =-1, with eigenvectors

    𝐗^(h) =  [ p_τ -1 + λ_h/d , q_τ -1/d ]^T  for  λ_h = (-1)^h,

where d = (p_τ -1 + λ_h, q_τ -1) for h ∈{ 0,1 }.




If the period τ of the continued fraction expansion of √(N) is even, a factor of 2N is located at positions τ/2 + jτ with j = 0, 1, …, in the sequence {Q_n}_n ≥ 0.



It is sufficient to consider j = 0, due to the periodicity of {Q_n}_n ≥ 0. Since τ is even, M_τ -1 is involutory and has eigenvalues λ_0 =1 and λ_1 =-1 with corresponding eigenvectors shown in (<ref>).
Considering Equation (<ref>) written as

    [ p_τ - j-2; q_τ - j-2 ] = (-1)^j-1M_τ -1[ p_j; q_j ],

we see that 𝐘^(j) =  [ p_j, q_j  ]^T is an eigenvector of M_τ -1, of eigenvalue (-1)^j-1, if and only if j satisfies the condition τ -j-2 = j, that is j=τ -2/2= τ_0. From the comparison of 𝐘^(j) and 𝐗^(h), we have

    p_τ_0 = p_τ -1 + (-1)^τ_0/d    q_τ_0 = q_τ -1/d,

where the equalities are fully motivated because gcd ( p_τ_0, q_τ_0 )=1, recalling that d = (p_τ -1 + (-1)^τ_0,q_τ -1). Direct computation yield

    (-1)^τ_0+1 Q_τ_0+1 = (p_τ -1 + (-1)^τ_0 -1)^2 - Nq_τ -1^2/d^2 = 2(-1)^τ_0p_τ -1 +1/d^2,

which can be written as p_τ_0^2 - Nq_τ_0^2 = 2(-1)^τ_0p_τ_0/d. Dividing this equality by 2p_τ_0/d we have

    d p_τ_0/2- N1/2p_τ_0/dq_τ_0^2 = (-1)^τ_0.

Noting that gcd ( p_τ_0, q_τ_0 )=1, it follows that 2p_τ_0/d is a divisor of 2N, i.e. Q_τ_0+1 = Q_τ/2| 2N.



In the case where τ is odd, we can state the following two results.


Let N be a positive integer such that the continued fraction expansion of √(N) has odd period τ. The representation of N = x^2 + y^2 is given by x = Q_(τ +1)/2 and y = P_(τ +1)/2.


Since τ is odd, by the anti-symmetry in the sequence { Q_n}_n=0^τ -1, we have Q_(τ +1)/2 = Q_(τ -1)/2, so that the quadratic form Q_(τ -1)/2 X^2 + 2P_(τ +1)/2XY - Q_(τ +1)/2 Y^2 has discriminant 4P_(τ +1)/2^2 - 4Q_(τ +1)/2 Q_(τ -1)/2 = 4N, which shows the assertion.

From this, we can deduce a result similar to that in Theorem <ref> for the case of an odd period.


Let N>0 be a composite non-square integer such that the continued fraction expansion of √(N) has odd period τ. If -1 is a quadratic nonresidue modulo N, then Q_(τ +1)/2 contains a nontrivial factor of N.


Using the previous theorem, N = Q_(τ +1)/2^2 + P_(τ +1)/2^2, and so P_(τ +1)/2^2 ≡ - Q_(τ +1)/2^2 N. If (N, Q_(τ +1)/2) = 1, then Q_(τ +1)/2^-1N exists. Therefore, (Q^-1_(τ +1)/2 P_(τ +1)/2 )^2 ≡ -1 N, and so -1 is a quadratic residue modulo N.



The following identity holds

    √(N) + P_m+1/Q_m+1 = - p_m-1 - q_m-1√(N)/p_m - q_m√(N)  ∀  m ≥ 0.



The proof is straightforward.


The following result will be useful in the next section.


If τ is even, defining γ as

    γ = ∏_m=0^τ -1 (√(N) + P_m+1),

we have γ/γ =(p_τ -1 + q_τ -1√(N))^2 = 𝔠_τ -1^2.
If τ is odd, defining ω as

    ω = ∏_m=0^2τ -1 (√(N) + P_m+1),

we have ω/ω =(p_2τ -1 + q_2τ -1√(N))^2 = 𝔠_2τ -1^2.



We provide a proof for the case τ even; the proof for the odd case follows the same procedure. We have

    γ/γ   = ∏_m=0^τ -1√(N) + P_m+1/-√(N) + P_m+1
       = ∏_m=0^τ -1(√(N) + P_m+1)^2/ P_m+1^2 -N
       = ∏_m=0^τ -1(√(N) + P_m+1)^2/-Q_m+1 Q_m.

Noting that ∏_m=0^τ -1 - Q_m Q_m+1 = ∏_m=0^τ -1 - Q_m+1^2 = ∏_m=0^τ -1 Q_m+1^2 due to the periodicity of the sequence {Q_m}_m ≥ 0 and the parity of τ, we deduce that γ/γ is a perfect square. Considering the identity of the Lemma <ref>, we have that the base of the square giving γ/γ is

    ∏_m=0^τ -1√(N) + P_m+1/Q_m+1   = ∏_m=0^τ -1 - p_m-1 - q_m-1√(N)/p_m - q_m√(N)
       = p_-1 - q_-1√(N)/p_τ -1 - q_τ -1√(N)
       = p_τ -1 + q_τ -1√(N).

Therefore, 

    ∏_m=0^τ -1√(N) + P_m+1/Q_m+1 = p_τ -1 + q_τ -1√(N) = 𝔠_τ -1,

and in conclusion γ/γ = 𝔠_τ -1^2. 

Similarly, if τ is odd, the following equation holds

    ∏_m=0^2τ -1√(N) + P_m+1/Q_m+1 = p_2τ -1 + q_2τ -1√(N) = 𝔠_2τ -1.







 §.§ Even period and nontrivial factor
 
In this subsection, we provide sufficient conditions on the integer N to ensure that the period τ is even and that Q_τ/2 2.

Theorem <ref> provides a sufficient condition on N for an even period: if N cannot be written as sum of two squares, then the period of the continued fraction expansion of √(N) is even. Using Fermat's theorem on sums of two squares, we derive the following sufficient condition on the factors of N.

Let N>0 be a non-square integer. If there exists a prime p ≡ 3 4 that divides N with odd exponent, then τ≡ 0 2.


The following theorem, proven by Mollin in <cit.>, provides necessary and sufficient conditions, involving Diophantine equations, to have τ even and Q_τ/2=2.


The following statements are equivalent for N > 2.

    
  * X^2 - NY^2 = ± 2 is solvable, indicating that at least one of the equations X^2 - NY^2 = 2 and X^2 - NY^2 =-2 has a solution.
    
  * τ is even and Q_τ/2=2.


This result implies that, if τ≡ 0 2 and the two Diophantine equations 

    X^2 - N Y^2 = 2   and   X^2 - N Y^2 = -2

are insoluble, then the central term Q_τ/2 2, and so it contains a proper factor of N. The following result, due to Yokoi, and presented in <cit.>, gives us sufficient conditions for the insolubility of the two equations in (<ref>). 


For any positive non-square integer N, if the Diophantine equation X^2 - N Y^2 = ± 2 has an integral solution, then 

    N ≡ 2 4  or   N ≡ 3 4.



Hence, if N has the following form 

    N=p^k_1q^k_2  with   p≡ q ≡ 3 4  and   k_1 ≡ k_2 ≡ 1 2,

then the period τ of the continued fraction expansion of √(N) is even and Q_τ/2 contains a proper factor of N. Integers of that form, also known as Blum integers, have been broadly employed in cryptographic applications. They were first introduced by L. Blum, M. Blum, and Shub in <cit.> in the context of pseudo-random number generators, and later applied in zero-knowledge proofs (see, e.g., <cit.>, <cit.>, <cit.>, <cit.>).



 §.§ Quadratic forms


An overview about binary quadratic forms can be found in <cit.>. 

A binary quadratic form is a polynomial F(x, y) = ax^2+bxy+cy^2, with a, b, c ∈ℤ.
The matrix associated with F is

    M_F =[   a b/2; b/2   c ].


We abbreviate a binary quadratic form with coefficients a,b,c as (a,b,c).


Two quadratic forms F and F' are equivalent if there exists a matrix C ∈ℤ^2 × 2 such that 

    M_F' = C^T M_F C

and (C) = ± 1. If (C) = 1, the forms are properly equivalent and we write F ∼ F'.



The discriminant Δ of a quadratic form (a,b,c) is Δ = b^2 -4ac. We define 𝔽_Δ as the set of all quadratic forms of discriminant Δ.

The discriminant is an invariant for the equivalence relation of quadratic forms ∼: if F ∈𝔽_Δ, and F ∼ F', then F' ∈𝔽_Δ.


A quadratic form (a,b,c), with positive discriminant Δ= b^2 -4ac is reduced if

    | √(Δ) -2   | a  |  | < b < √(Δ).



Given a quadratic form F, it is always possible to find a reduced quadratic form equivalent to F. In the following we are going to prove it giving a reduction algorithm on quadratic forms of positive non-square discriminant. To do so, we first need the next definition.


For any form F = (a, b, c) with ac  0 of discriminant Δ, a non-square positive integer, we define the standard reduction operator ρ by

    ρ(a, b, c) =  (c, r(-b, c), r(-b, c)^2 - Δ/4c ),

where r(-b, c) is defined to be the unique integer r such that r +b ≡ 0 2c and

    -|c| < r ≤ |c|   if  √(Δ)< |c|,


    √(Δ) - 2|c| < r < √(Δ)  if   |c| < √(Δ).

ρ(F) is called the reduction of F. The inverse reduction operator is defined by

    ρ^-1(a,b,c) =  ( r(-b, c)^2 - Δ/4c, r(-b, a), a  ).


We denote ρ^n(F) the result of n applications of ρ on F. The identities ρ(ρ^-1(F)) = ρ^-1(ρ(F)) = F hold when F is reduced. We point out the fact that (a,b,c)∼ρ(a,b,c) through the transformation given by the matrix 

    [  0 -1;  1  t ],

where r(-b,c)=-b+2ct.
The proof of the following, fundamental, proposition can be found in <cit.>.

 

    
  * The number of iterations of ρ which are necessary to reduce a form (a, b, c) is at most 2 + ⌈log_2 ( |c|/ √(Δ) ) ⌉.
    
  * If F = (a, b, c) is a reduced form, then ρ(a, b, c) is again a reduced form.





If (a, b, c), of discriminant Δ >0, is reduced, then | a  |, b and | c  | are less than √(Δ), and a and c are of opposite signs (<cit.>). This implies that the number of reduced quadratic forms of discriminant Δ is finite.



Two forms F(x, y) = ax^2+bxy+cy^2 and F'(x, y) = a'x^2+b'xy+c'y^2 are adjacent if c = a' and b + b' ≡ 0 2c.
 

Given a reduced quadratic form F, there exists a unique reduced quadratic form equivalent to F and adjacent to F. This form is ρ(F). As we have seen in Remark <ref>, there exists a finite number of reduced quadratic forms of positive discriminant Δ, and so this process eventually repeats, forming a cycle. The important aspect of this is that the cycle is actually all of the reduced forms equivalent to the first form, as proven in <cit.>.


We call the principal form the reduced form of discriminant Δ having as first coefficient 1. It is denoted with 1 and the cycle in which it lies is called the principal cycle. 




Let Υ= {_n }_n ≥ 0 be the sequence of binary quadratic forms defined as

    _m(x,y) = (-1)^mQ_m x^2 + 2 P_m+1 xy + (-1)^m+1Q_m+1y^2,   for m ≥0.





If τ is even, then the sequence Υ is periodic of period τ, and if τ is odd Υ is periodic of period 2τ. This is due to Theorem <ref> and Theorem <ref>.


 
Let F=(a_1, b_1, c_1) and G=(a_2, b_2, c_2) two quadratic forms having same discriminant Δ. The Gauss composition of F and G is

    F ∘ G = (a_3, b_3, c_3) =  ( d_0 a_1 a_2/n^2, b_1 + 2 a_1/n ( s(b_2 - b_1)/2 - c_1 v ), b_3^2 - Δ/4 a_3 ),

where β = (b_1 + b_2)/2, n = (a_1, a_2, β), s, u, v such that a_1 s + a_2 u + β v = n, and d_0 = (a_1,a_2,β, c_1,c_2,(b_1-b_2)/2). Although the composition is not unique, all compositions of given forms F and G are equivalent.

We remark that all quadratic forms in Υ have the same discriminant Δ =4N, where N>0 is the non-square integer we want to factorize. This implies that for all (a,b,c) ∈Υ, we have

    Δ≡ b^2 ≡ b 2,

and so b ≡ 0 2. Therefore, the value β in the Definition <ref> is an integer.
A quadratic form F=(a,b,c) is primitive if gcd( a, b, c ) = 1. As proven in the next proposition, the forms in Υ are primitive.


The forms _n are primitive for all n≥0.


We prove the statement by induction, for n ≥0. 

Base step (n=0): The base step is proven by the following chain of equalities

    (Q_0, 2 P_1, Q_1) = (1, 2a_0, a_0^2 - N) =1.

Inductive step (n ⇒ n+1): The thesis follows from the next equalities

    (Q_n, 2 P_n+1, Q_n+1)    = (Q_n, 2(a_nQ_n -P_n), (N-P_n+1^2)/Q_n) 
       = (Q_n,-2P_n, Q_n-1 -a_n^2Q_n +2a_nP_n) 
       = (Q_n,  -2P_n, Q_n-1) =1.


This implies that in the Gauss composition of two elements of Υ, the coefficient d_0 is always equal to 1.


Using the definition of ρ and Equation (<ref>), it is straightforward to prove that

    ρ^n(_0)=_n.

The form _0 is reduced and equal to 1, so the quadratic forms in Υ are reduced, and Υ is the principal cycle. Moreover, for any pair of forms _n, _m ∈Υ, their Gauss composition _n ∘_m is equivalent to _0. This follows from the property proven by Gauss in Article 237-239 of <cit.>: if F ∼ G, then H ∘ F ∼ H ∘ G, for all quadratic forms F, G, H having same discriminant. In particular, we obtain

    _n ∘_m ∼_n ∘_0 ∼_n ∼_0,
 
using also the fact that _0 ∘_n ∼_n for all n ≥ 0. Consequently, applying the Gauss composition to any couple of elements of Υ, followed by a sufficient number of applications of ρ to obtain a reduced form, results in an element of Υ.


As mentioned previously, we are interested in quickly finding the coefficient Q_τ/2 when τ is even, or Q_(τ + 1)/2 when τ is odd. Consequently, we aim to determine the quadratic form _τ/2 or _(τ - 1)/2 in an efficient manner (i.e. with time complexity O(ln(N)^α), with α constant). The value of τ could be too large (see Remark <ref>), so we need a way to make longer jumps within the principal cycle. 
As we will see, Gauss composition, followed by the reduction, will allow us to make long jumps in Υ. To estimate the length of these jumps we use the (well-known) infrastructural distance δ. A comprehensive definition and detailed description of distance is provided in <cit.>.  


Given a quadratic form F=(a,b,c) of discriminant Δ>0, the infrastructural distance δ of F and ρ(F) is

    δ(F, ρ(F)) = 1/2ln | b + √(Δ)/b - √(Δ) |.

Given n>0, the  distance δ of F and ρ^n(F) is

    δ(F,ρ^n(F)) = ∑_i=1^n δ(ρ^i-1(F),ρ^i(F)).




We now restrict ourselves to forms in the principal cycle. We then have the following proposition.


Let _n and _m be two reduced forms in the principal cycle, and let _0 be the unit form. Then if we define G = _n ∘_m, G may not be reduced, but let _r be a (non unique) form obtained from G by the reduction algorithm, i.e. by successive
applications of ρ. Then we have

    δ(_0,_r) = δ(_0,_n) + δ(_0,_m) + δ(G,_r),

and furthermore,

    | δ(G,_r)  | < 2 lnΔ,

where Δ is the discriminant of these forms.

The above proposition follows from the property that δ is exactly additive under composition (before any reductions are made) (see <cit.>) and the estimation of the bound for | δ(G,_r) | discussed in Section 12 of <cit.>.

The next proposition is fundamental for a computational point of view. Indeed, given _i, _j, with j > i, and their distance δ(_i, _j), it gives an estimation of j-i. In particular, if δ(_i, _j) = D, then 2D/ln (4N)<j-i< 2D/ln 2 + 1.


Let F ∈𝔽_Δ reduced. The following two bounds hold (see <cit.> for the proofs):

    
  * δ(F, ρ(F)) < 1/2lnΔ,
    
  * δ( F , ρ^2( F )) > ln 2, and the same holds for ρ^-1.





In our case, the discriminant of the forms in Υ is Δ=4N, where N is an odd non-square integer. Theorem <ref> (proven also in <cit.>) and Theorem <ref> show that the distance between quadratic forms is considered modulo 

    R^+(N) = ln(p_τ -1 + q_τ -1√(N))= ln(𝔠_τ -1)    if τ≡ 0 2
    ln(p_2τ -1 + q_2τ -1√(N))= ln(𝔠_2τ -1)    if τ≡ 1 2
     .



If τ is even, the distance δ(_0, _τ) (the distance of a period) is exactly equal to ln(𝔠_τ -1) and the distance δ(_0,_τ/2) is exactly equal to 1/2δ(_0, _τ).


The distance between _τ and _0 is the summation

    d(_0,_τ) = ∑_i=0^τ -1 d(_i, _i+1) = ∑_i=1^τ1/2ln( √(N) + P_i/√(N) - P_i) = 1/2ln ( ∏_i=1^τ√(N) + P_i/√(N) - P_i ).

Recalling that N- P_i^2 = Q_i Q_i-1 >0, and taking into account the periodicity of the sequence {Q_n}_n ≥ 0, the last expression can be written with rational denominator as

    1/2ln ( ∏_i=1^τ(√(N) + P_i)^2/Q_i Q_i-1 ) = 1/2ln ( ∏_i=1^τ(√(N) +P_i)^2/ Q_i^2 ) = ln ( ∏_i=1^τ√(N) + P_i/ Q_i ).

The conclusion follows from Equation (<ref>). The equality d(_0, _τ/2) = 1/2d(_0,_τ) is an immediate consequence of the symmetry of the sequence {P_n}_n ≥ 1 within a period.


We now give the same result for the case of odd period.


If τ is odd, the distance δ(_0, _2τ) (the distance of a period) is exactly equal to ln(𝔠_2τ -1) and the distance δ(_0, _τ) is equal to ln(𝔠_2τ -1)/2.


The distance between _2τ and _0 is the summation

    d(_0,_2τ) = ∑_i=0^2τ -1 d(_i, _i+1) = ∑_i=1^2τ1/2ln( √(N) + P_i/√(N) - P_i) = 1/2ln ( ∏_i=1^2τ√(N) + P_i/√(N) - P_i ).

Recalling that N- P_i^2 = Q_i Q_i-1 >0, and taking into account the periodicity of the sequence {Q_n}_n ≥ 0, the last expression can be written with rational denominator as

    1/2ln ( ∏_i=1^2τ(√(N) + P_i)^2/Q_i Q_i-1 ) = 1/2ln ( ∏_i=1^2τ(√(N) +P_i)^2/ Q_i^2 ) = ln ( ∏_i=1^2τ√(N) + P_i/ Q_i ).

The conclusion follows from Equation (<ref>) and the periodicity of { P_n }_n ≥ 1.



If τ is odd, the distance δ(_0, _(τ -1)/2) (the distance of a target form) is equal to ln(𝔠_2τ -1)/4 + O(ln (N)).


From the previous theorem, we know that δ(_0, _τ)=ln(𝔠_2τ -1)/2. Moreover, using the symmetry (<ref>), we obtain the following equality

    δ(_0,_τ)    = 2 δ(_0,_(τ - 5)/2) + 1/2ln( √(N) + P_(τ - 1)/2/√(N) - P_(τ - 1)/2) 
       + 1/2ln( √(N) + P_τ - 1/√(N) - P_τ - 1) + 1/2ln( √(N) + P_τ/√(N) - P_τ).

Therefore,

    δ(_0,_(τ - 5)/2)    = 1/4ln(𝔠_2τ -1) - 1/4ln( √(N) + P_(τ - 1)/2/√(N) - P_(τ - 1)/2) 
       - 1/4ln( √(N) + P_τ - 1/√(N) - P_τ - 1) - 1/4ln( √(N) + P_τ/√(N) - P_τ),

and so

    δ(_0,_(τ - 1)/2)    = 1/4ln(𝔠_2τ -1) - 1/4ln( √(N) + P_(τ - 1)/2/√(N) - P_(τ - 1)/2) 
       - 1/4ln( √(N) + P_τ - 1/√(N) - P_τ - 1) - 1/4ln( √(N) + P_τ/√(N) - P_τ) 
       + 1/2ln( √(N) + P_(τ - 1)/2/√(N) - P_(τ - 1)/2) + 1/2ln( √(N) + P_(τ + 1)/2/√(N) - P_(τ + 1)/2).

Using Proposition <ref>, the upper bound holds


    1/4 | ln( √(N) + P_(τ - 1)/2/√(N) - P_(τ - 1)/2) - ln( √(N) + P_τ - 1/√(N) - P_τ - 1) -  ln( √(N) + P_τ/√(N) - P_τ)
    + 2 ln( √(N) + P_(τ + 1)/2/√(N) - P_(τ + 1)/2)  | 
       < 5/4ln (4N),

from which, we deduce that

    δ(_0,_(τ - 1)/2) = 1/4ln(𝔠_2τ -1) + O(ln (N)).


The following remark is fundamental from a computational point of view, because it provides an upper bound on the distance of a cycle in Υ.



We have that R^+(N) = n R(N), with n ≤ 6. Hua, in <cit.>, proves that 
    R(N) ≤√(N) (1/2ln N + 1  )    if N ≡ 1 4
    
    2√(N) (1/2ln (4N) + 1  )    if N ≡ 3 4

and so, R^+(N) = O(√(N)ln N). However, we do not know which is the largest value that R(N) can attain as a function of N. It is conjectured that there exists an infinite set of values of N such that R(N) ≫√(N)lnln N (see <cit.> for large-scale numerical experiments and more details).






§ THE FACTORIZATION ALGORITHM


In this section, we present our factorization method.
The integer N>0 to be factorized is odd, non-square and composite.
In the first part of this section, we describe the method and provide the pseudocodes (Algorithms <ref> and <ref>). We then prove the correctness of our approach and analyze its computational cost.
This method is a modification of the one presented by Elia <cit.>. 
We assume that R^+(N) has been preliminarily computed. In the final part of this section, we mention a method for computing an integer multiple of R^+(N).

To simplify the notation, we introduce the following definition.


Given two forms _n, _m∈Υ, the giant step of _n and _m is the composition 
    _n∙_m = ρ^r(_n∘_m),
 realized through the Gauss composition _n∘_m, followed by the minimum number r of reduction operations ρ to obtain a reduced form. The notation _n^t represents t successive applications of the giant step of _n with itself, i.e., _n∙⋯∙_n (repeated t times).


We define our method for both the even-period and odd-period cases. To enhance readability, we define the following quantity

    (N) =  R^+(N)/2    if τ≡ 0 2
    
    R^+(N)/4    if τ≡ 1 2
    ,

which represents the distance of the quadratic form we want to reach (or an approximation of it). Indeed, if τ≡ 0 2, then δ(_0, _τ/2) = R^+(N)/2 = (N) and if τ≡ 1 2, then δ(_0, _(τ + 1)/2) = R^+(N)/4 + O(ln(N)) = (N) + O(ln(N)), using Theorem <ref> and Corollary <ref>.

We distinguish two cases: R^+(N) ≤ (ln N)^2 and R^+(N) > (ln N)^2. In the first case, we compute _i = ρ^i(_0) until a non-trivial factor of N is found among their coefficients, if such a factor exists. By Proposition <ref>, the number of reduction steps ρ is at most 2 δ(_0, _τ/2)/ln 2 + 1 = R^+(N)/ln 2 + 1 = O((ln N)^2) when the period is even, and 2 δ(_0, _(τ - 1)/2)/ln 2 + 1 ≤R^+(N)/2 ln 2 + 5/2 ln 2ln (4N) + 1 = O((ln N)^2), when the period is odd. If the number of iterations exceeds R^+(N)/ln 2 + 5/2 ln 2ln (4N) + 1 the procedure is stopped: our algorithm cannot find a factor of N. The pseudocode for this method is given in Algorithm <ref>.

If R^+(N) > (ln N)^2 we proceed in the following way.


    
  * First phase: In this phase we compute an approximation of _τ/2, if τ is even, or _(τ + 1)/2, if τ is odd. Starting from _0, we compute the forms _i in the principal cycle, for i=0, …, ℓ, until δ(_0, _ℓ) ≥ 2ln (4N) +1 and δ(_0, _ℓ) ≤ 4ln N (this is possible using Proposition <ref> and the definition of distance). Then, we compute the quadratic forms _ℓ^2^i, using giant steps, and their exact distance d_i = δ(_0, _ℓ^2^i), using Proposition <ref>, for i=1, …, t, with t such that d_t-1≤(N) < d_t. We point out the fact that d_i+1>d_i for all i ≥ 0. Then, using the forms _ℓ, …, _ℓ^2^t-1, we compute F, an approximation of _τ/2.
    To do so, first we set F= _ℓ^2^t -1 and d = d_t-1. Then we start by computing d + d_t-2 and we check if it is greater or smaller than (N); in the second case we update F with F∙_ℓ^2^t -2 and d with d + d_t-2. We iterate this procedure for i=t -3, …, 0 by computing d+ d_i, comparing it with (N), and, if it is smaller, updating F with F∙_ℓ^2^i. 

    
  * Second phase: Starting from F, we iterate the operators ρ and ρ^-1 until a factor of N is found. An upper bound on the number of iterations of ρ and ρ^-1 needed to find a factor (both in the case of even and odd period) is given by: 
    
    Ψ(R^+(N), N) := 2/ln 2 ( 4 ln (4N) log_2 ( R^+(N)/2 ) + 13/4ln (4N)  ) + 1.

    This bound derives from the results demonstrated later in this section. 

A priori, we do not know the parity of τ, so we proceed as follows (as outlined in Algorithm <ref>). First, we run the procedure assuming τ≡ 0 2, which implies (N) = R^+(N)/2. If, at the end of the second phase, after Ψ(R^+(N), N) steps, a factor is not found, we then try again assuming τ≡ 1 2, which implies (N) = R^+(N)/4. If, even in this case, no factor is found after Ψ(R^+(N), N) steps during the second phase, the output is -1: our method cannot factor N.





1.5em

1.5em



1.5em

1.5em

 

In the following we present two propositions, which are fundamental to the analysis of our method for the case R^+(N)> (ln N)^2. The first shows that t (the number of powers of G_0) is always “small", the second proves that our approximation of _τ/2, if τ even, or _(τ + 1)/2 if τ is odd, is good.


The value of t in Algorithm <ref> is at most ⌈log_2(N) ⌉.


Using Proposition <ref> and the above notation, we have that

    δ(_0, G_i) > 2δ(_0, G_i-1) - 
    2ln (4N)   ∀  i >0,

and so

    δ(_0, G_i)    > 2^iδ(_0, G_0) - 
    2 ∑_k=0^i-1 2^k ln(4N) 
       = 2^iδ(_0, G_0) - 2(2^i -1)ln(4N) 
       ≥ 2^i(2 ln(4N) +1) - 2(2^i -1)ln(4N) 
       = 2^i + 2ln(4N).

Therefore, for i ≥⌈log_2 (N) ⌉, we have δ(_0, G_i) > (N).


This proposition implies that t = O(ln N), thanks to Remark <ref>.


Let F be the quadratic form obtained at the end of the second phase of the method (using the notation of Algorithm <ref>). Then, the following holds

    | (N) - δ(_0, F)  | = O((ln N)^2).




In the for loop at line 20 of Algorithm <ref>, we have at most t-1 giant steps. Therefore, using the previous proposition and Proposition <ref>, we have that

    | d - δ(_0, F)  | = O((ln N)^2).

Now, we prove that | d- (N)  | = O((ln N)^2).
We define I ⊆{0, …, t-1 } the set of indexes of the distances d_0, …, d_t-1 that appear in the computation of d, i.e.

    d = ∑_i ∈ I d_i.

We distinguish two cases:

    
  * Case 0 ∉ I: Then we have d + d_0 > (N), and so
    
    d≤(N) < d + d_0,

    from which we deduce that
    
    0 ≤(N) - d < d_0 ≤ 4 ln N.

    
  * Case 0 ∈ I: Let j = min{ i ∈ℕ| i ∉ I }, then 1 ≤ j ≤ t-2. We have that
    
    d + d_0    = ∑_i ∈ I∖{ 0 } d_i + 2 d_0 
       = ∑_i ∈ I∖{ 0 } d_i + d_1 + O(ln N) 
       = ∑_i ∈ I∖{ 0,1 } d_i + 2d_1 + O(ln N) 
       = ∑_i ∈ I∖{ 0,1 } d_i + d_2 + O(ln N) 
         ⋮
       = ∑_i ∈ I∖{ 0, …, j-1} d_i + d_j + γ(N)

    where γ(N)=O((ln N)^2).
    By construction, we have that 
    ∑_i ∈ I∖{ 0, …, j-1} d_i + d_j > (N),
 from which
    
    d≤(N) ≤d + d_0 - γ(N),

    and so 
    
    0 ≤(N) - d≤ d_0 - γ(N) = O((ln N)^2).


Therefore,

    |(N) - δ(_0, F)  | ≤ |(N) - d | +  | d - δ(_0, F)  | = O((ln N)^2).


Therefore, if τ is even, then | δ(_0, _τ /2) - δ(_0, F)  | = O((ln N)^2). If τ is odd, we have that 
    | δ(_0, _(τ + 1 )/2) - δ(_0, F)  |    =  | (N) -δ(_0, F)  | +  O(ln N) = O((ln N)^2).

since δ(_0, _(τ + 1 )/2) = (N) + O(ln N) by Corollary <ref>.
We analyze the computational cost of this algorithm, assuming that R^+(N)> (ln N)^2 is preliminarily computed. The cost of the computation of the form G_0 such that δ(_0, G_0) ≥ 2 ln(4N)+1 is at most 2(2 ln(4N)+1)/ln 2 + 1, using Proposition <ref>.

We then analyze the cost of the while loop at line 9. The number t of the giant steps is at most ⌈log_2(N) ⌉ = O(ln N), proven in Proposition <ref>. Each giant step requires: O(ln N) elementary operations for the extended Euclidean algorithm, used to compute s, u and v in the Gauss composition, and at most O(ln N) applications of ρ. 
Indeed, if we apply the Gauss composition of two forms in Υ, using the extended Euclidean algorithm, we obtain (a,b,c) such that | c  | = O(N^4). This follows from Proposition <ref> and the classical bounds on the solution of Bézout's identity via the extended Euclidean algorithm. Hence, by Proposition <ref>, the number of applications of ρ is O(ln N). Therefore, the cost of the while loop at line 9 is O((ln N)^2). 

The computation of F requires at most O((ln N)^2) steps: at most O(ln N) giant steps, and at most O(ln N) application of ρ for each giant step.

Finally, as proven in Proposition <ref>, the distance between the approximation F and _τ/2, or _(τ +1)/2, is at most O((ln N)^2), and so, using again the fact that, for each reduced form F, δ(F, ρ^2(F))> ln 2, are needed only O((ln N)^2) applications of ρ to reach _τ/2 or _(τ +1)/2.

In conclusion, the method has a computational complexity of O((ln N)^2).
It is remarked that the cost of elementary arithmetic operations (i.e. additions, subtractions, multiplications and divisions of big integers) and logarithm valuations are not counted.





Using Proposition <ref>, Proposition <ref>, Corollary <ref>, and Proposition <ref>, it is possible to derive the following upper bound on the number of iterations of ρ and ρ^-1 in the second phase of the method

    Ψ(R^+(N), N) = 2/ln 2 ( 4 ln (4N) log_2 ( R^+(N)/2 ) + 13/4ln (4N)  ) + 1.

This bound holds for both the cases when τ is even and τ is odd.




We point out that it is not necessary to have R^+(N) precomputed; it is sufficient to have an integer multiple of it: R'(N) = k R^+(N), with k ∈ℕ. For simplicity, we describe how the method is modified in the case of an even period. In this case, running Algorithm <ref> with R'(N) instead of R^+(N) is equivalent to considering the principal cycle with multiplicity k (i.e. k times the principal cycle). Our target form is located in the middle of some period of distance R^+(N)/2, so:

    
  * if k is odd, a factor of 2N is found (as coefficient of a form) at the position at distance k R^+(N)/2, from the beginning;
    
  * if k is even, the quadratic form _τ-1 is found in a position at distance k R^+(N)/2 (which reveals a posteriori that k is even); in this case, the procedure can be repeated, targeting the form at position at distance k R^+(N)/4 from _0. Again, either a factor of 2N is found, or k is found to be a multiple of 4. Clearly the process can be iterated h times until k R^+(N)/2^h is an odd multiple of R^+(N), and a factor of 2N is found.

If k, as a function of N, is O(N^α) with α constant, then the computational cost of the algorithm does not change. Indeed, in this case, the value of t in Algorithm <ref> is at most ⌈log_2 ( k(N) ) ⌉ = O(ln N) (see Proposition <ref>), and the number of iterations of ρ and ρ^-1 is at most Ψ(kR^+(N), N) = O((ln N)^2). The odd-period case follows similarly.


As outlined in Remark <ref>, we are focused on studying and researching methods to efficiently calculate or approximate R^+(N), or one of its integer multiples that is “not too large". In particular, we are seeking a method that efficiently computes kR^+(N), where k, as a function of N, is O(N^α), with α constant. Since R^+(N) = n R(N), with n ≤ 6, our problem is equivalent to finding an efficient algorithm that computes (a multiple of) the regulator of ℚ(√(N)). The method due to Vollmer, described in <cit.>, currently has the best known complexity. 
It is a Monte Carlo algorithm that computes R(N) in time O  ( exp (  3/√(8)√(ln N lnln N) )  ), assuming the GRH. 
Other methods for computing the regulator are detailed in <cit.>.



In conclusion, using this approach for the precomputation of R(N) and the algorithm previously described, we have obtained a factorization method of (conjectured) time complexity O  ( exp (3/√(8)√(ln N lnln N) )  ), which is more efficient than CFRAC and SQUFOF.



§ ACKNOWLEDGMENTS

The first author is member of GNSAGA of INdAM and acknowledges support from Ripple’s University Blockchain Research Initiative. The second author is partially supported by project SERICS (PE00000014 - https://serics.eu) under the MUR National Recovery and Resilience Plan funded by European Union - NextGenerationEu.










