Find divisors $n_i$ of $a_i$, mutualy coprime, with $\prod n_i=\operatorname{lcm}(a_i)$

arithmeticcomputational mathematicsgcd-and-lcmleast-common-multiple

We are given $k$ positive integers $a_i$, and want as many positive integers $n_i$ each dividing the respective $a_i$, with the $n_i$ mutually coprime, and the product of the $n_i$ equal to to the Least Common Multiple of the $a_i$.

Note: $1$ is coprime with every positive integer, so $n_i=1$ is sometime fine.

There's a simple method when we can put the $a_i$ as $\prod p_j^{m_{j,i}}$ for primes $p_j$. We initially set all the $n_i$ to $a_i$; then for each prime $p_j$ appearing in the factorizations of the $a_i$, we choose one of possibly several $i$ with maximal multiplicity $m_{j,i}$, and for each other $i$ divide $n_i$ by $p_j^{m_{j,i}}$.

Can we proceed with no factorization whatsoever, or if impossible with as little factorization work as possible (perhaps using a Greatest Common Divisor primitive and/or other computationally efficient methods)?

Given $a_0=554400=2^5\cdot 3^2\cdot 5^2\cdot 7\cdot 11$, $a_1=4422600=2^3\cdot 3^5\cdot 5^2\cdot 7\cdot 13$
$k=2$, the only 4 solutions are:

  • $n_0=2^5\cdot 5^2\cdot 7\cdot 11=61600$, $n_1=3^5\cdot 13=3159$
  • $n_0=2^5\cdot 5^2\cdot 11=8800$, $n_1=3^5\cdot 7\cdot 13=22113$
  • $n_0=2^5\cdot 7\cdot 11=2464$, $n_1=3^5\cdot 5^2\cdot 13=78975$
  • $n_0=2^5\cdot 11=352$, $n_1=3^5\cdot 5^2\cdot 7\cdot 13=552825$

Solving the problem for $k=2$ leads to a solution for larger $k$, by iteratively replacing $a_i$ and $a_{i'}$ with the $n_i$ and $n_{i'}$ obtained for the $k(k-1)/2$ pairs of $i,i'$ with $i<i'$. However that seems less than optimal.

The problem is of general interest in large CRT problems. The initial motivating context was applying the Chinese Remainder Theorem in the final step of the Pohlig-Hellman algorithm with composite moduli, there.

Best Answer

There is currently no known polynomial time algorithm for finding prime factorizations of integers, so here is an algorithm which, while far from perfect, still runs in polynomial time.

Let $A=\{a_1,\ldots,a_k\}$, let $M$ be the maximum element of $A$, and let $n=\log_2{M}$. We will determine the computational complexity in terms of $n$ and $k$.

Every prime $p$ has a largest power $m_p$ amongst elements of $A$, which is to say $m_p$ is the power of $p$ in $\text{lcm}\hspace{1pt}A$. Our first goal will be to find values $b_i$, which will be a multiple of only those primes $p$ for which the power of $p$ in $a_i$ is $m_p$. We will then use these to find values $c_i$, which will be the product of values $p^{m_p}$ such that $p^{m_p}$ divides $a_i$. Finally, we will consider these against each other to determine $n_i$. This will be completed in $\mathcal{O}(k\cdot n^5+k^2\cdot n)$ time, which is indeed a polynomial.

For each $a_i$, we define

$$b_i=\frac{a_i^{n+1}}{\gcd(a_i^{n+1},\text{lcm}\hspace{1pt}\{a_j^n\text{ s.t. }j\neq i\}\})}$$

Let us figure out what is happening here. What the denominator does is find the the maximum power of every prime which divides an element of $A\backslash\{a_i\}$, and multiplies that maximum by $n$. It then takes these powers and compares them with those found in $a_i^{n+1}$, and takes the smaller one. If the power of $p$ found in $a_i$ is $p^{m_p}$, where $m_p$ must be at most $n$, then $p^{m_p}\cdot(n+1)>p^{m_p}\cdot n$, so $p$ will divide $b_i$, and otherwise, since $(p^{m_p}-1)\cdot(n+1)<p^{m_p}\cdot n$, $p$ will not divide $b_i$. Thus $p$ divides $b_i$ if and only if $p^{m_p}$ divides $a_i$. This process will take $\mathcal{O}(k\cdot n^4)$ time, resulting from taking $k$ $\gcd$ functions on values with at most about $n^2$ bits.

Now we wish to find $c_i$, for which the power of $p$ is $m_p$ if $p^{m_p}$ divides $a_i$ and 0 otherwise. To do this, we will do the following (note the subscripts have been removed here for clarity):

x = gcd(a,b)
while (x != 1)
    c *= x
    a /= x
    x = gcd(a,b)

Every iteration, this will multiply the power of $p$ in $c_i$ by at least 1 so long as it is not yet $p^{m_p}$. Thus it will take at most $m_p\leq n$ iterations, each taking $\mathcal{O}(n^5)$ time due to the multiplication, so the process of finding values of $c_i$ takes at most $\mathcal{O}(k\cdot n^5)$ time.

Finally, we need to find our values of $n_i$. We find these in order. We simply let $n_1=c_1$, then take the $\gcd$ of this with the other values of $c_i$ and divide these old values of $c_i$ by that $\gcd$ to get new values. We then do this with every other value from 2 to $k$, taking $\mathcal{O}(k^2\cdot n+n^2)$ time, giving us working values of $n_i$ in polynomial time.

Related Question