I will use a typical sorting algo. something like a quicksort maybe. something like:

void quicksort(int start, int end, int prime_idx, int cum) {
  if(start == end) return;
     int i = (start - 1);
    //to handle case where arr[j]==1
    for(int j = start; j<=end; ++j) {
      if(arr[j]==cum) {
         swap(start,j);
         start++;
      }
    }
    for (int j = start; j <= end- 1; j++)
    {
        // 
        if (arr[j]%(cum*prime[prime_idx])==0)
        {
            i++;    // increment index of smaller element
            swap(arr[i],arr[j]);
        }
    }
    swap(arr[i+1], arr[end]);
    quicksort(start, i, prime_idx, cum*prime[prime[idx]]);
    quicksort(i+1, end, prime_idx+1, cum);
}

Beware this code is not tested and i havent written a quicksort in years. Also I am not sure about the complexity, but I believe it should be fast enough. EDIT: fix the arr[i]==1 part.

If your comparison function works in you don't need to waste memory on storing the factors, you can factor the numbers while comparing them.

Here's a nice solution:

First I'll denote the elements limit by lim. Find all the prime numbers from 1 to lim using sieve.

Now, we will make a vector keys of size lim + 1. every number X from 1 to lim will have in keys[X] a value, which is its location if we sorted all numbers from 1 to lim by the lexicographical order of their prime factorization (that is, 1 will have key 1, 2 will have key 2, 4 (2 * 2) will have key 3...).

To do this, we will create all the prime factorizations by lexicographical order. That is, we will first create all those that are just with 2's, then we will invlove 3 in the prime factorization, and so on. For instance, the first numbers we will go through are:

1, 2, 4, 8, 16, ...., 524288 (2^19), 786432 (2^18 * 3^1)...

The way to do this is with a very simple recursion: we will maintain the current value we're on, and the last prime factor we took. then, we go through all the primes from the last prime we took, until the end of the vector, or until (the current prime we're on) * (the current value in the recursion) is bigger than lim.

We maintain the last prime factor since we don't want things like this to happen: 12 = 2 * 3 * 2. (This is just by the definition of the prime factorization — we go by the increasing factors.)

The above can be proved to run in O(lim), and gives us a comparison function that works in O(1) (we just compare keys).

Here is a sample code:

const int lim = 1e6;
int n, nextkey = 1;
vector<int> a, p, keys;

void sieve() {
	vector<bool> b(1e6 + 1, true);
	int i = 2;
	for (; i * i < b.size(); i++) {
		if (b[i]) {
			p.push_back(i);
			for (int j = i * i; j < b.size(); j += i) b[j] = false;
		}
	}
	for (; i < b.size(); i++) if (b[i]) p.push_back(i);
}

void calc(int val = 1, int lastind = 0) {
	keys[val] = nextkey;
	nextkey++;

	while (lastind < p.size() && (long long)val * p[lastind] <= lim) {
		calc(val * p[lastind], lastind);
		lastind++;
	}
}

bool cmp(int &a, int &b) {
	return keys[a] < keys[b];
}

int main() {
	keys.resize(lim + 1);
	sieve();
	calc();

	cin >> n;
	a.resize(n);
	for (int i = 0; i < n; i++) cin >> a[i];
	
	sort(a.begin(), a.end(), cmp);

	for (auto &i : a) cout << i << " "; cout << endl;
}

This works in . (Note that you can also make the sorting in O(n) using countsort, because the keys are between 1 and lim).