Algorithm for checking if number is factorable into a set of prime numbers

Question

I was wondering if there is an algorithm that checks wether a given number is factorable into a set of prime numbers and if no give out the nearest number. The problem comes always up when I use the FFT.

Thanks a lot for your help guys.

Answer 1

In general this looks like a hard problem, particularly finding the next largest integer that factors into your set of primes. However, if your set of primes isn't too big, one approach would be to turn this into an integer optimization problem by taking the logs. Here is how to find the smallest number > n that factors into a set of primes p_1...p_k

choose integers x_1,...,x_k to minimize (x_1 log p_1 + ... + x_k log p_k - log n)
Subject to:
  x_1 log p_1 + ... + x_k log p_k >= log n
  x_i >= 0 for all i

The x_i will give you the exponents for the primes. Here is an implementation in R using lpSolve:

minfact<-function(x,p){
  sol<-lp("min",log(p),t(log(p)),">=",log(x),all.int=T)
  prod(p^sol$solution)
}

> p<-c(2,3,13,31)
> x<-124363183
> y<-minfact(x,p)
> y
[1] 124730112
> factorize(y)
Big Integer ('bigz') object of length 13:
 [1] 2  2  2  2  2  2  2  2  3  13 13 31 31
> y-x
[1] 366929
> 

Using big integers, this works pretty well even for large numbers:

> p<-c(2,3,13,31,53,79)
> x<-as.bigz("1243631831278461278641361")
> y<-minfact(x,p)
y
> 
Big Integer ('bigz') :
[1] 1243634072805560436129792
> factorize(y)
Big Integer ('bigz') object of length 45:
 [1] 2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2 
[26] 2  2  2  2  2  2  2  2  3  3  3  3  13 31 31 31 31 53 53 53
> 

Answer 2

Your question is about well-known factorization problem - which could not be resolved with 'fast' (polynomial) time. Lenstra's elliptic algorithm is the most efficient (known) way in common case, but it requires strong knowledge of numbers theory - and it's also sub-exponential (but not polynomial).

Other algorithms are listed in the page by first link in my post, but such things as direct try (brute force) are much more slower, of cause.

Please, note, that under "could not be resolved with polynomial time" - I mean that there's no way to do this now - but not that such way does not exist (at least now, number theory can not provide such solution for this problem)

Answer 3

Here is a brute force method in C++. It returns the factorization of the nearest factorable number. If N has two equidistant factorable neighbours, it returns the smallest one.

GCC 4.7.3: g++ -Wall -Wextra -std=c++0x factorable-neighbour.cpp

#include <iostream>
#include <vector>

using ints = std::vector<int>;

ints factor(int n, const ints& primes) {
  ints f(primes.size(), 0);
  for (int i = 0; i < primes.size(); ++i) {
    while (0< n && !(n % primes[i])) {
      n /= primes[i];
      ++f[i]; } }

  // append the "remainder"
  f.push_back(n);
  return f;
}

ints closest_factorable(int n, const ints& primes) {
  int d = 0;
  ints r;
  while (true) {
    r = factor(n + d, primes);
    if (r[r.size() - 1] == 1) { break; }
    ++d;
    r = factor(n - d, primes);
    if (r[r.size() - 1] == 1) { break; }
  }
  r.pop_back();
  return r; }

int main() {
  for (int i = 0; i < 30; ++i) {
    for (const auto& f : closest_factorable(i, {2, 3, 5, 7, 11})) {
      std::cout << f << " "; }
    std::cout << "\n"; }
}

Answer 4

I suppose that you have a (small) set of prime numbers S and an integer n and you want to know is n factors only using number in S. The easiest way seems to be the following:

P <- product of s in S
while P != 1 do
    P <- GCD(P, n)
    n <- n/P
return n == 1

You compute the GCD using Euclid's algorithm.

The idea is the following: Suppose that S = {p1, p2, ... ,pk}. You can write n uniquely as

n = p1^n1 p2^n2 ... pk^nk * R 

where R is coprime wrt the pi. You want to know whether R=1.

Then

GCD(n, P) = prod ( pi such that ni <> 0 ). 

Therefore n/p decrease all those non zeros ni by 1 so that they eventually become 0. At the end only R remains.

For example: S = {2,3,5}, n = 5600 = 2^5*5^2*7. Then P = 2*3*5 = 30. One gets GCD(n, p)=10=2*5. And therefore n/GCD(n,p) = 560 = 2^4*5*7.

You are now back to the same problem: You want to know if 560 can be factored using S = {2,5} hence the loop. So the next steps are

• GCD(560, 10) = 10. 560/GCD = 56 = 2^3 * 7.
• GCD(56, 10) = 2. 56/2 = 28 = 2^2 * 7
• GCD(28, 2) = 2. 28/2 = 14 = 2 * 7
• GCD(14, 2) = 2. 14/2 = 7
• GCD(7, 2) = 1 so that R = 7 ! Your answer if FALSE.

Answer 5

kissfft has a function

int kiss_fft_next_fast_size(int n)

that returns the next largest N that is an aggregate of 2,3,5.

Also related is a kf_factor function that factorizes a number n, pulling out the "nice" FFT primes first (e.g. 4's are pulled out before 2's)