Fermat's factorization helper

Question

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\$.

The task

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
  );
});

Answer 1

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.

Answer 2

Python 2, 127 121 117 111 107 104 101 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.

Answer 3

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.

Answer 4

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