How to calculate all the prime factors of lcm of n integers?

Question

I have n integers stored in an array a, say a[0],a[1],.....,a[n-1] where each a[i] <= 10^12 and n <100. Now, I need to find all the prime factors of the LCM of these n integers i.e., LCM of {a[0],a[1],.....,a[n-1]}

I have a method but I need a more efficient one.

My method:

 First calculate all the prime numbers up to 10^6 using sieve of Eratosthenes.

 For each a[i]
      bool check_if_prime=1;
      For all prime <= sqrt(a[i])
             if a[i] % prime[i] == 0 {
                store prime[i]
                check_if_prime=0
             }
      if check_if_prime
             store a[i]     // a[i] is prime since it has no prime factor <= sqrt(n) 
  Print all the stored prime[i]'s

Is there any better approach to this problem?

I'm posting the link to the problem:

http://www.spoj.pl/problems/MAIN12B/

Link to my code: http://pastebin.com/R8TMYxNz

Solution:

As suggested by Daniel Fischer my code needed some optimizations like a faster sieve and some minor modifications. After doing all those modification, I'm able to solve the problem. This is my accepted code on SPOJ which took 1.05 seconds:

#include<iostream>
#include<cstdio>
#include<map>
#include<bitset>

using namespace std;

#define max 1000000

bitset <max+1> p;
int size;
int prime[79000];

void sieve(){
    size=0;
    long long i,j;
    p.set(0,1);
    p.set(1,1);
    prime[size++]=2;
    for(i=3;i<max+1;i=i+2){
        if(!p.test(i)){
            prime[size++]=i;
            for(j=i;j*i<max+1;j++){
                p.set(j*i,1);
            }
        }
    }
}

int main()
{
    sieve();

    int t;
    scanf("%d", &t);

    for (int w = 0; w < t; w++){
        int n;
        scanf("%d", &n);
        long long a[n];

        for (int i = 0; i < n; i++)
            scanf("%lld", &a[i]);

        map < long long, int > m;
        map < long long, int > ::iterator it;
        for (int i = 0; i < n; i++){
            long long num = a[i];
            long long pp;
            for (int j = 0; (j < size) && ((pp = prime[j]) * pp <= num); j++){
                int c = 0;
                for ( ; !(num % pp); num /= pp)
                    c = 1;
                if (c)
                    m[pp] = 1;
            }

            if ((num > 0) && (num != 1)){
                m[num] = 1;
            }
        }
        printf("Case #%d: %d\n", w + 1, m.size());
        for (it = m.begin(); it != m.end(); it++){
            printf("%lld\n", (*it).first);
        }
    }

    return 0;
}

In case, anyone is able to do it in a better way or by some faster method, please let me know.

Answer 1

There seems to be several useful algorithms at http://en.wikipedia.org/wiki/Least_common_multiple

In particular http://en.wikipedia.org/wiki/Least_common_multiple#A_method_using_a_table
seems to be appropriate.

It turns the nested loops 'inside out' and works on all the numbers simultaneously, one prime at a time.

Because it uses one prime at a time, you could find the next prime as it is needed, avoiding generating 10^6 primes before starting. As each number is reduced by its prime factors, the maximum prime needed to test the numbers might be reduced, so it needs even less work.

Edit: It also makes termination unambiguous and easy to check because the number is reduced to one when all of its factors have been found. In fact, when all numbers are reduced to one it can terminate, though I didn't use that property in my code.

Edit: I read the problem, the algorithm at http://en.wikipedia.org/wiki/Least_common_multiple#A_method_using_a_table solves it directly.

SWAG: There are 78,498 primes below 10^6 (http://primes.utm.edu/howmany.shtml#table) So worst case, there are 100 numbers which need to be tested for 78,498 primes, = 7,849,800 'mod' operations

No number can factored successfully by a prime (one mod and one divide) more than log2(10^12) = 43 mods and divides, so 4,300 divides and 4,300 mods, so dominated by the prime tests. To keep it simple, lets call it 8,000,000 integer divides and mods. It needs to generate the primes, but as already stated by Daniel Fischer that is quick. The rest is book keeping.

So, on a modern processor, I'd WAG about 1,000,000,000 divides or mods/second, so run time around 10ms x 2?

Edit: I used the algorithm at http://en.wikipedia.org/wiki/Least_common_multiple#A_method_using_a_table

No smarts, just exactly as explained there.

I was badly out on my estimate, about 10x, but still only 20% of maximum allowed run time.

Performance (with some printing to confirm results)

real 0m0.074s
user 0m0.062s
sys  0m0.004s

for 100 numbers:

999979, 999983, 999979, 999983, 999979, 999983, 999979, 999983, 999979, 999983,

10 times, to ensure almost all primes had to be tested as that seems to be the major computation.

and also with the same amount of printing, but a value of almost 10^12

real 0m0.207s
user 0m0.196s
sys  0m0.005s

for 100 of 999962000357L, // ((long)999979L)*((long)999983L)

gcc --version i686-apple-darwin10-gcc-4.2.1 (GCC) 4.2.1 (Apple Inc. build 5666) (dot 3) Copyright (C) 2007 Free Software Foundation, Inc. This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

Model Name: MacBook Pro Processor Name: Intel Core 2 Duo Processor Speed: 2.16 GHz

Summary: It clearly works, and run time is about 20% of allowable maximum, on a relatively old processor, which is comparable to Daniel Fischer's implementation.

Question: I am a new contributor here, so it seems a bit harsh to -2 my answer when:
a. it seems to be accurate, complete, and satisfies all of the criteria, and
b. I wrote the code, tested it, and provided the results
What did I do wrong? How do I get feedback so that I can improve?

Answer 2

With those constraints, a few not too large numbers, the best way to find the prime factorisation of their least common multiple is indeed the factorisation of each number. Since there are only 78498 primes below 10⁶, trial division will be fast enough (unless you're really desperate for the last drop of performance), and sieving the primes to 10⁶ is also a matter of few milliseconds.

If speed is of the utmost importance, a combined approach of trial division and a deterministic Miller-Rabin type prime test with a factorisation method like Pollard's rho-algorithm or an elliptic curve factorisation method is probably a bit faster (but the numbers are so small that the difference won't be large, and you need a number type with more than 64 bits for the primality test and factorisation to be fast).

In the factorisation, you should of course remove the prime factors as they are found

if (a[i] % prime[k] == 0) {
    int exponent = 0;
    do {
        a[i] /= prime[k];
        ++exponent;
    }while(a[i] % prime[k] == 0);
    // store prime[k] and exponent
    // recalculate bound for factorisation
}

to reduce the limit to which you need to check the primes.

Your main problem, as far as I can see, is that your sieve is too slow and uses too much space (which is in part responsible for its slowness). Use a bit sieve to get better cache locality, remove even numbers from the sieve and stop checking whether to cross off multiples at the square root of max. And you're allocating way too much space for the prime array.

for(int j=0;(prime[j]*prime[j] <= num) && (j<size);j++){

You must check j < size before accessing prime[j].

    while(num%prime[j]==0){
        c=1;
        num /= prime[j];
        m[prime[j]]=1;
    }

Don't set m[prime[j]] multiple times. Even if std::maps are pretty fast, that's slower than setting it only once.

Answer 3

FWIW, in response to the original request for a faster way to get all of the primes up to one million, here is a much faster way to do that.

This uses a windowing Wheel-based Sieve of Eratosthenes, with a Wheel size of 30 and the window size set to the square root of the search's upper limit (1000 for a search up to 1,000,000).

As I am not proficient in C++, I have coded it in C#, on the assumption that this should be easily convertible to C++. However, even in C# it can enumerate all of the primes up to 1,000,000 in about 10 milliseconds. Even generating all of the primes up to a billion only takes 5.3 seconds, and I'd imagine that this would be even faster in C++.

public class EnumeratePrimes
{
    /// <summary>
    /// Determines all of the Primes less than or equal to MaxPrime and
    /// returns then, in order, in a 32bit integer array.
    /// </summary>
    /// <param name="MaxPrime">The hishest prime value allowed in the list</param>
    /// <returns>An array of 32bit integer primes, from least(2) to greatest.</returns>
    public static int[] Array32(int MaxPrime)
    {
        /*  First, check for minimal/degenerate cases  */
        if (MaxPrime <= 30) return LT30_32_(MaxPrime);

        //Make a copy of MaxPrime as a double, for convenience
        double dMax = (double)MaxPrime;

        /*  Get the first number not less than SQRT(MaxPrime)  */
        int root = (int)Math.Sqrt(dMax);
        //Make sure that its really not less than the Square Root
        if ((root * root) < MaxPrime)  root++;

        /*  Get all of the primes <= SQRT(MaxPrime) */
        int[] rootPrimes = Array32(root);
        int rootPrimeCount = rootPrimes.Length;
        int[] primesNext = new int[rootPrimeCount];

        /*  Make our working list of primes, pre-allocated with more than enough space  */
        List<int> primes = new List<int>((int)Primes.MaxCount(MaxPrime));
        //use our root primes as our starting list
        primes.AddRange(rootPrimes);

        /*  Get the wheel  */
        int[] wheel = Wheel30_Spokes32();

        /*  Setup our Window frames, starting at root+1 */
        bool[] IsComposite; // = new bool[root];
        int frameBase = root + 1;
        int frameMax = frameBase + root; 
        //Pre-set the next value for all root primes
        for (int i = WheelPrimesCount; i < rootPrimeCount; i++)
        {
            int p = rootPrimes[i];
            int q = frameBase / p;
            if ((p * q) == frameBase) { primesNext[i] = frameBase; }
            else { primesNext[i] = (p * (q + 1)); }
        }

        /*  sieve each window-frame up to MaxPrime */
        while (frameBase < MaxPrime)
        {
            //Reset the Composite marks for this frame
            IsComposite = new bool[root];

            /*  Sieve each non-wheel prime against it  */
            for (int i = WheelPrimesCount; i < rootPrimeCount; i++)
            {
                // get the next root-prime
                int p = rootPrimes[i];
                int k = primesNext[i] - frameBase;
                // step through all of its multiples in the current window
                while (k < root) // was (k < frameBase) ?? //
                {
                    IsComposite[k] = true;  // mark its multiple as composite
                    k += p;                 // step to the next multiple
                }
                // save its next multiple for the next window
                primesNext[i] = k + frameBase;
            }

            /*  Now roll the wheel across this window checking the spokes for primality */
            int wheelBase = (int)(frameBase / 30) * 30;
            while (wheelBase < frameMax)
            {
                // for each spoke in the wheel
                for (int i = 0; i < wheel.Length; i++)
                {
                    if (((wheelBase + wheel[i] - frameBase) >= 0)
                        && (wheelBase + wheel[i] < frameMax))
                    {
                        // if its not composite
                        if (!IsComposite[wheelBase + wheel[i] - frameBase])
                        {
                            // then its a prime, so add it to the list
                            primes.Add(wheelBase + wheel[i]);
                        }
                        // // either way, clear the flag
                        // IsComposite[wheelBase + wheel[i] - frameBase] = false;
                    }
                }
                // roll the wheel forward
                wheelBase += 30;
            }

            // set the next frame
            frameBase = frameMax;
            frameMax += root;
        }

        /*  truncate and return the primes list as an array  */
        primes.TrimExcess();
        return primes.ToArray();
    }

    // return list of primes <= 30
    internal static int[] LT30_32_(int MaxPrime)
    {
        // As it happens, for Wheel-30, the non-Wheel primes are also
        //the spoke indexes, except for "1":
        const int maxCount = 10;
        int[] primes = new int[maxCount] {2, 3, 5, 7, 11, 13, 17, 19, 23, 29 };

        // figure out how long the actual array must be
        int count = 0;
        while ((count <= maxCount) && (primes[count] < MaxPrime)) { count++; }

        // truncte the array down to that size
        primes = (new List<int>(primes)).GetRange(0, count).ToArray();
        return primes;
    }
    //(IE: primes < 30, excluding {2,3,5}.)

    /// <summary>
    /// Builds and returns an array of the spokes(indexes) of our "Wheel".
    /// </summary>
    /// <remarks>
    /// A Wheel is a concept/structure to make prime sieving faster.  A Wheel
    /// is sized as some multiple of the first three primes (2*3*5=30), and
    /// then exploits the fact that any subsequent primes MOD the wheel size
    /// must always fall on the "Spokes", the modulo remainders that are not
    /// divisible by 2, 3 or 5.  As there are only 8 spokes in a Wheel-30, this
    /// reduces the candidate numbers to check to 8/30 (4/15) or ~27%.
    /// </remarks>
    internal static int[] Wheel30_Spokes32() {return new int[8] {1,7,11,13,17,19,23,29}; }
    // Return the primes used to build a Wheel-30
    internal static int[] Wheel30_Primes32() { return new int[3] { 2, 3, 5 }; }
    // the number of primes already incorporated into the wheel
    internal const int WheelPrimesCount = 3;
}

/// <summary>
/// provides useful methods and values for working with primes and factoring
/// </summary>
public class Primes
{
    /// <summary>
    /// Estimates PI(X), the number of primes less than or equal to X,
    /// in a way that is never less than the actual number (P. Dusart, 1999)
    /// </summary>
    /// <param name="X">the upper limit of primes to count in the estimate</param>
    /// <returns>an estimate of the number of primes between 1  and X.</returns>
    public static long MaxCount(long X)
    {
        double xd = (double)X;
        return (long)((xd / Math.Log(xd)) * (1.0 + (1.2762 / Math.Log(xd))));
    }
}

Answer 4

result := []

for f in <primes >= 2>:
   if (any a[i] % f == 0):
      result = f:result
   for i in [0..n-1]:
      while (a[i] % f == 0):
         a[i] /= f
   if (all a[i] == 1):
      break

Note: this only gives the list of prime factors of the LCM, not the actual value of the LCM (i.e. does not calculate exponents), which I believe is all the question required.