# Fermat's factorization helper

We'd like to factorize a semiprime $N$. The goal of this challenge is to find two small integers $u$ and $v$ such that $uvN$ can be trivially factorized with Fermat's method, thus allowing to easily deduct the factors of $N$.

Given a semiprime $N$ and a positive integer $k$, we define $x$ and $y$ as:

$$x=\lceil\sqrt{kN}\rceil$$ $$y=x^2-kN$$

Step #1 - Find $k$

You first need to find the smallest possible value of $k$ such that $y$ is a square number (aka perfect square).

This allows to factorize $kN$ with a single iteration of Fermat's factorization method. More concretely, this immediately leads to:

$$kN=(x+\sqrt{y})\times(x-\sqrt{y})$$

(Update: this sequence is now published as A316780)

Step #2 - Factorize $k$

You then have to find the two positive integers $u$ and $v$ such that:

$$uv=k$$ $$cu=x+\sqrt{y}$$ $$dv=x-\sqrt{y}$$

where $c$ and $d$ are the prime factors of $N$.

Summary

Your task is to write a program or function that takes $N$ as input and prints or outputs $u$ and $v$ in any order and any reasonable format.

## Example

Let's consider $N = 199163$

Step #1

The smallest possible value of $k$ is $40$, which gives:

$$x = \lceil(\sqrt{40 \times 199163})\rceil = 2823$$ $$y = 2823^2 - 40 \times 199163 = 7969329 - 7966520 = 2809 = 53^2$$ $$kN = (2823 + 53) \times (2823 - 53)$$ $$kN = 2876 \times 2770$$

Step #2

The correct factorization of $k$ is $k = 4 \times 10$, because:

$$kN = 2876 \times 2770$$ $$kN = (719 \times 4) \times (277 \times 10)$$ $$N = 719 \times 277$$

So, the correct answer would be either $[ 4, 10 ]$ or $[ 10, 4 ]$.

## Rules

• It is not required to strictly apply the two steps described above. You're free to use any other method, as long as it finds the correct values of $u$ and $v$.
• You must support all values of $uvN$ up to the native maximum size of an unsigned integer in your language.
• The input is guaranteed to be a semiprime.
• This is code-golf, so the shortest answer in bytes wins.
• Standard loopholes are forbidden.

## Test cases

N          | k    | Output
-----------+------+------------
143        | 1    | [   1,  1 ]
2519       | 19   | [   1, 19 ]
199163     | 40   | [   4, 10 ]
660713     | 1    | [   1,  1 ]
4690243    | 45   | [   9,  5 ]
11755703   | 80   | [  40,  2 ]
35021027   | 287  | [   7, 41 ]
75450611   | 429  | [ 143,  3 ]
806373439  | 176  | [   8, 22 ]
1355814601 | 561  | [  17, 33 ]
3626291857 | 77   | [   7, 11 ]
6149223463 | 255  | [  17, 15 ]
6330897721 | 3256 | [  74, 44 ]


## Example implementation

In the snippet below, the $f$ function is an ungolfed implementation which takes $N$ as input and returns $u$ and $v$.

For illustrative purposes only, the snippet also includes the $g$ function which takes $N$, $u$ and $v$ as input and computes the factors of $N$ in $O(1)$.

f = N => {
for(k = 1;; k++) {
x = Math.ceil(Math.sqrt(k * N));
y = x * x - k * N;
ySqrt = Math.round(Math.sqrt(y));

if(ySqrt * ySqrt == y) {
p = x + ySqrt;

for(u = 1;; u++) {
if(!(p % u) && !(N % (p / u))) {
v = k / u;
return [ u, v ];
}
}
}
}
}

g = (N, u, v) => {
x = Math.ceil(Math.sqrt(u * v * N));
y = x * x - u * v * N;
ySqrt = Math.round(Math.sqrt(y));
p = x + ySqrt;
q = x - ySqrt;
return [ p / u, q / v ];
}

[
143, 2519, 199163, 660713, 4690243, 11755703,
35021027, 75450611, 806373439, 1355814601,
3626291857, 6149223463, 6330897721
]
.map(N => {
[u, v] = f(N);
[c, d] = g(N, u, v);
console.log(
'N = ' + N + ', ' +
'u = ' + u + ', ' +
'v = ' + v + ', ' +
'N = ' + c + ' * ' + d
);
});

• Are we guaranteed that the input N will in fact be a semiprime? – Greg Martin Apr 16 '17 at 18:56
• @GregMartin Yes you are. – Arnauld Apr 16 '17 at 19:00

# Mathematica, 81 79 bytes

Thanks to Martin Ender for saving 2 bytes!

(c=Ceiling;For[j=0;z=E,c@z>z,p=(x=c@Sqrt[j+=#])+{z=Sqrt[x^2-j],-z}];p/#~GCD~p)&


Pure function taking a semiprime as input and returning an ordered pair of positive integers. The For loop implements the exact procedure described in the question (using # for the input in place of n), with x as defined there, although we store j = k*n instead of k itself and z=Sqrt[y] instead of y itself. We also compute p={x+z,x-z} inside the For loop, which ends up saving one byte (on like the seventh try). Then the two desired factors are (x+z)/GCD[#,x+z] and (x-z)/GCD[#,x-z], which the concise expression p/#~GCD~p computes directly as an ordered pair.

Curiosities: we want to loop until z is an integer; but since we're going to use Ceiling already in the code, it saves two bytes over !IntegerQ@z to define c=Ceiling (which costs four bytes, as Mathematica golfers know) and then test whether c@z>z. We have to initialize z to something, and that something had better not be an integer so that the loop can start; fortunately, E is a concise choice.

• NB: I do not accept the shortest answer anymore as I used to do in my first code-golf challenges. I usually just let accepted answers as-is unless they are outgolfed, for the sake of consistency. This is what happened here. – Arnauld May 5 '20 at 8:27
• no worries, thanks for the note :) – Greg Martin May 5 '20 at 17:03

## JavaScript (ES7), 86 81 bytes

n=>(g=k=>(y=(n*k)**.5+1|0,y+=(y*y-n*k)**.5)%1?g(k+1):n*u++%y?g(k):[--u,k/u])(u=1)


Edit: Saved 4 bytes thanks to @Arnauld.

# APL (Dyalog Unicode), 6049 45 bytes

{×1|w←⍵∨(⌈+.5*⍨×⍨-⍨⌈×⌈).5*⍨⍵×⍺:⍺∇1+⍵⋄w,⍵÷w}∘1


Try it online!

-11b thanks to Adám, -4b courtesy of ngn along with some excellent tips from both on The APL Orchard

## Explanation:

First ever APL answer so might be able to be golfed further! Explanation below is for the 49b solution. Approach pretty much just follows the algorithm in the original post, but uses the GCD to compute the factors once k, x and y are known.

{×1|y←.5*⍨(×⍨x←⌈.5*⍨⍵×⍺)-⍵×⍺:⍺∇1+⍵⋄w,⍵÷w←⍵∨x+y}∘1 ⍝
{                                             }∘1 ⍝ dfn; ⍺ is left argument (N), ⍵ is right argument acting as a counter for k incremented each recursion from 1
(  x←⌈.5*⍨⍵×⍺)-⍵×⍺                      ⍝ Computes x as per question for current k
y←.5*⍨(×⍨x         )                          ⍝ Computes sqrt(y) from x as per question (×⍨x = x^2). Note that y here actually represents sqrt(y) from the question as we do not need to store actual y
×1|y                                             ⍝ Use 1 modulo y and 'direction' to return 0 if y is a whole number, else 1
:     ⋄               ⍝ 'guard' - if the evaluation of the left side of : is true, performs the statement immediately following it, otherwise the next statement (separated by ⋄)
⍺∇1+⍵                ⍝ if guard is true, recurse with the same left argument of N, and a right argument of the new k (k+1)
w,⍵÷w←⍵∨x+y    ⍝ if guard is false, computes gcd (∨) of k and (x+y) as one factor (w), and k÷w as the second factor


# Python 2, 127121117111107104101 99 bytes

-1 byte thanks to Neil & -3 bytes thanks to ovs

N=input()
k=u=1;p=m=.5
while p%1:p=1+(k*N)**m//1;p+=(p*p-k*N)**m;k+=1
while N*u%p:u+=1
print~-k/u,u


Try it Online!

Curiosities:

p is initialized to .5 so that the loop condition will be true on the first iteration. Note that it is shorter to store p (as x + sqrt(y)) than it is to store each of x and y separately.

• x*x instead of x**2? – Neil Apr 16 '17 at 21:49
• @Neil Yes, of course. Thanks – math junkie Apr 16 '17 at 22:15

# Axiom, 131 115 bytes

v(x)==floor(x^.5)::INT;r(n)==(k:=0;repeat(k:=k+1;x:=1+v(k*n);y:=v(x*x-k*n);x^2-y^2=k*n=>break);[w:=gcd(k,x+y),k/w])


The function that would resolve question is r(n) above. ungolf and test

vv(x)==floor(x^.5)::INT

--(x-y)*(x+y)=k*n
rr(n)==
k:=0
repeat
k:=k+1
x:=1+vv(k*n)
y:=vv(x*x-k*n)
x^2-y^2=k*n=>break
[w:=gcd(k,x+y),k/w]

(4) -> [[i,r(i)] for i in [143,2519,199163,660713,4690243,11755703]]
(4)
[[143,[1,1]], [2519,[1,19]], [199163,[4,10]], [660713,[1,1]],
[4690243,[9,5]], [11755703,[40,2]]]
Type: List List Any