Euler's factorization method

Theoretical basis

The Brahmagupta–Fibonacci identity states that the product of two sums of two squares is a sum of two squares. Euler's method relies on this theorem but it can be viewed as the converse, given ${\displaystyle n=a^2+b^2=c^2+d^2}$ we find ${\displaystyle n}$ as a product of sums of two squares.

First deduce that

${\displaystyle a^2-c^2=d^2-b^2}$

and factor both sides to get

${\displaystyle (a-c)(a+c)=(d-b)(d+b)}$ (1)

Now let ${\displaystyle k=\gcd (a-c,d-b)}$ and ${\displaystyle h=\gcd (a+c,d+b)}$ so that there exists some constants ${\displaystyle l,m,l^{\prime },m^{\prime }}$ satisfying

• ${\displaystyle (a-c)=kl}$,
• ${\displaystyle (d-b)=km}$,

${\displaystyle \gcd (l,m)=1}$

• ${\displaystyle (a+c)=h{m}^{\prime }}$,
• ${\displaystyle (d+b)=h{l}^{\prime }}$,

${\displaystyle \gcd (l^{\prime },m^{\prime })=1}$

Substituting these into equation (1) gives

${\displaystyle klh{m}^{\prime }=kmh{l}^{\prime }}$

Canceling common factors yields

${\displaystyle l{m}^{\prime }=l^{\prime }m}$

Now using the fact that ${\displaystyle (l,m)}$ and ${\displaystyle \left(l^{\prime },m^{\prime }\right)}$ are pairs of relatively prime numbers, we find that

• ${\displaystyle l=l^{\prime }}$
• ${\displaystyle m=m^{\prime }}$

So

• ${\displaystyle (a-c)=kl}$
• ${\displaystyle (d-b)=km}$
• ${\displaystyle (a+c)=hm}$
• ${\displaystyle (d+b)=hl}$

We now see that ${\displaystyle m=\gcd (a+c,d-b)}$ and ${\displaystyle l=\gcd (a-c,d+b)}$

Applying the Brahmagupta–Fibonacci identity we get

${\displaystyle \left(k^2+h^2\right)\left(l^2+m^2\right)=(kl+hm{)}^2+(km-hl{)}^2=((a-c)+(a+c){)}^2+((d-b)-(d+b){)}^2=(2a{)}^2+(2b{)}^2=4n.}$

As each factor is a sum of two squares, one of these must contain both even numbers: either ${\displaystyle (k,h)}$ or ${\displaystyle (l,m)}$. Without loss of generality, assume that pair ${\displaystyle (k,h)}$ is even. The factorization then becomes

${\displaystyle n=\left({\left(\frac{k}{2}\right)}^2+{\left(\frac{h}{2}\right)}^2\right)\left(l^2+m^2\right).}$

Worked example

Since: ${\displaystyle 1000009={1000}^2+3^2={972}^2+{235}^2}$

we have from the formula above:

Thus,

${\displaystyle 1000009=\left[{\left(\frac{4}{2}\right)}^2+{\left(\frac{34}{2}\right)}^2\right]\cdot \left[{\left(\frac{14}{2}\right)}^2+{\left(\frac{116}{2}\right)}^2\right]}$

${\displaystyle =\left(2^2+{17}^2\right)\cdot \left(7^2+{58}^2\right)}$

${\displaystyle =(4+289)\cdot (49+3364)}$

${\displaystyle =293\cdot 3413}$

Pseudocode

function Euler_factorize(int n) -> list[int] if is_prime(n) then print("Number is not factorable") exit function for-loop from a=1 to a=ceiling(sqrt(n)) b2 = n - a*a b = floor(sqrt(b2)) if b*b==b2 break loop preserving a,b if a*a+b*b!=n then print("Failed to find any expression for n as sum of squares") exit function for-loop from c=a+1 to c=ceiling(sqrt(n)) d2 = n - c*c d = floor(sqrt(d2)) if d*d==d2 then break loop preserving c,d if c*c+d*d!=n then print("Failed to find a second expression for n as sum of squares") exit function A = c-a, B = c+a C = b-d, D = b+d k = GCD(A,C)//2, h = GCD(B,D)//2 l = GCD(A,D)//2, m = GCD(B,C)//2 factor1 = k*k + h*h factor2 = l*l + m*m return list[ factor1, factor2 ]