# Speeding up FactorInteger for product of two primes

I have a large integer $$N$$ of size about 10^150. I know that $$N$$ is a product of two primes $$p$$ and $$q$$. I also know that both $$p$$ and $$q$$ are of roughly equal size, so one of them is not, for example, 2 or 3.

Can I use this information to somehow speed up FactorInteger ? Or is factoring such a large number already beyond the scope of a average desktop computer. I also noticed that FactorInteger uses one of my cpu cores only, is it possible to parallelize this?

• Maybe start around Sqrt[n] (avoiding N because of its built in definition), which, since the prime counting function is about $\text{li}(x)$, is approximately equal to Prime[Floor[LogIntegral[Sqrt[n]]]]? Then work your way down Prime[k] manually (as long as you start above Sqrt[n])? No idea if this is remotely a good idea, though, since I don't know how mathematica deals with Prime. But you could probably parallelize it easily, at least! :P – thorimur Mar 3 at 8:28
• How about implementing the general number field sieve as it's the most efficient: en.wikipedia.org/wiki/General_number_field_sieve – Dominic Mar 3 at 10:27
• @Dominic thanks, reading this led me to Fermat's optimization method which seems exactly what I need – spore234 Mar 3 at 13:26
• I'm sure you know RSA numbers take a long time to factor even with efficient algorithms and yours seems in this group. Note the list here: en.wikipedia.org/wiki/RSA_numbers are all composed of equal-sized factors. Even RSA-100 I think is practically beyond desktop computers. – Dominic Mar 3 at 13:37
• @Dominic : Until polynomial selection stops being a black art, this is not a productive suggestion. See, e.g., [math.ttu.edu/~cmonico/software/ggnfs/] for GGNFS, and other implementations have similar weaknesses in polynomial selection. – Eric Towers Mar 3 at 20:50

It all depends on "p and q are of roughly equal size". If they differ no more than say 10^6, the following works rather fast. E.g.:

fac[big_] := Module[{ne = NextPrime[Sqrt[big] - 1]},
While[Mod[big, ne] > 0, ne = NextPrime[ne]];
{ne,big/ne}
]
n = 80;
fac[NextPrime[10^n] NextPrime[10^n + 10^6] ]

(* {100000000000000000000000000000000000000000000000000000000000000000000\
000001000507, \
1000000000000000000000000000000000000000000000000000000000000000000000\
00000000129} )*

• ! Divisible[big, ne] is slightly faster than Mod[big, ne] > 0 because it does not need to actually compute the modulus. – Roman Mar 3 at 10:29

If x is around 10^150, let y = floor(sqrt(x)), and if y^2≠x then solve (y+k+1)(y-k)=x, (y+k+2)(y-k) = x, (y + k + 3)(y - k) = x etc. until you find a k that is an integer. This will cover a huge range, much larger than 10^6.

We know 0 <= d = x - y^2 <= 2y ≈ 2 x 10^75.

Take for example (y+k+2)(y-k) = y^2 + yk + 2y - ky - k^2 - 2k. We solve 2y - k^2 - 2k = d or k^2 + 2k + 1 = 2y +1 - d. k will be between sqrt(y) and sqrt(2y), that is between 4.5 x 10^37 and 6.3 x 10^37. You can halve the number of calculations by not calculating k if the solution would have to consist of an even and an odd number.

(y+k+c)(y-k) = x will have a solution around sqrt (cy) = sqrt(2c) * x^(1/4). j calculations let you cover even c up to 2c, with factors up to 4 sqrt(j) x^(1/4) apart.

• Have you tried to program your idea? – yarchik Mar 3 at 19:54
• @gnasher729, I assume this a known algorithm. If it is, what is it called? – Ted Ersek Mar 4 at 21:38

Translating some of @gnasher729's suggestions to Mathematica gives indeed a superfast recipe:

$MaxExtraPrecision = 10^3; f[x_Integer] := Module[{y, i, z, s}, y = Floor[Sqrt[x]]; Catch[ i = 0; While[True, z = (2 y + i)^2 - 4 x; If[z >= 0 && IntegerQ[Sqrt[z]], s = i/2 + y + {1, -1} Sqrt[z]/2; If[IntegerQ[s[[1]]] && IntegerQ[s[[2]]], Throw[s]]]; i++]]] f[15] (* {5, 3} *) f[16] (* {4, 4} *) f[17] (* {17, 1} *) n = 80; f[NextPrime[10^n] NextPrime[10^n + 10^6]] // AbsoluteTiming (* {0.010925, {100000000000000000000000000000000000000000000000000000000000000000000000001000507, 100000000000000000000000000000000000000000000000000000000000000000000000000000129}} *) f[NextPrime[10^n] NextPrime[10^n + 10^12]] // AbsoluteTiming (* {0.007086, {100000000000000000000000000000000000000000000000000000000000000000001000000000191, 100000000000000000000000000000000000000000000000000000000000000000000000000000129}} *)  Things get a lot slower once the difference between the two numbers becomes larger than their square root: f[NextPrime[10^n] NextPrime[10^n + 10^43]] // AbsoluteTiming (* {132.132, {100000000000000000000000000000000000010000000000000000000000000000000000000000357, 100000000000000000000000000000000000000000000000000000000000000000000000000000129}} *)  The remaining point is a fast method to determine IntegerQ[Sqrt[z]] for very large values of z. Using @MichaelE2's fast perfect-square test gives a 150-fold speedup, which is however no match for the exponentially increasing difficulty with factor-spacing: $MaxExtraPrecision = 10^3;
f[x_Integer] := Module[{y, i, z, s},
y = Floor[Sqrt[x]];
Catch[
i = 0;
While[True,
z = (2 y + i)^2 - 4 x;
If[z >= 0 && FractionalPart[Sqrt[z + 01]] == 0,
s = i/2 + y + {1, -1} Sqrt[z]/2;
If[IntegerQ[s[[1]]], Throw[s]]];
i++]]]

n = 80;
f[NextPrime[10^n] NextPrime[10^n + 10^43]] // AbsoluteTiming
(*    {0.894308,
{100000000000000000000000000000000000010000000000000000000000000000000000000000357,
100000000000000000000000000000000000000000000000000000000000000000000000000000129}}    *)

f[NextPrime[10^n] NextPrime[10^n + 10^44]] // AbsoluteTiming
(*    {92.1891,
{100000000000000000000000000000000000100000000000000000000000000000000000000000179,
100000000000000000000000000000000000000000000000000000000000000000000000000000129}}

• For IntegerQ[Sqrt[z]] try (Count[JacobiSymbol[#, Prime /@ Range[1, Ceiling[1 + Log[2, #]]]], -1] == 0 &) to get strong probably square-ness. Test with, e.g., the non-square 853973422267356706546355086954657449503488853576566997801404568608559958358411214284271742350323575809430922219 and the square 853973422267356706546355086954657449503488853576566997800752173762579354251844397804033888338676797897369495524. – Eric Towers Mar 3 at 21:16
• It's easy to prove that a False result is correct. I don't know if there are (named analogously) Euler-Jacobi pseudo-squares. Also the upper limit, Ceiling[...], could be lowered substantially since most primes are greater than $2$. – Eric Towers Mar 3 at 21:19
• Thanks @EricTowers for the suggestion. Your method as presented is not much faster than the plain IntegerQ[Sqrt[z]] and I've decided to go with @MichaelE2's simple-yet-superfast method for now. – Roman Mar 4 at 11:53
• It was >30-times faster in RepeatedTiming[]s on my hardware, but of course, I'm not using your hardware. – Eric Towers Mar 4 at 20:58