How to calculate all lcm prime factors of n integers?

I have n integers stored in array a, for example a [0], a [1], ....., a [n-1], where each is a[i] <= 10^12

and n <100

. Now I need to find all the LCM prime factors of these n integers, that is, 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 a better approach to this problem?

I am posting a link to the problem:

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

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

Decision:

As Daniel Fischer suggested, my code needed some optimizations like a faster sieve and some minor modifications. After making all these changes, I can fix 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 can do it better or with a faster method, please let me know.

+3


source to share


4 answers


With these constraints, a few not too large numbers, the best way to find the primary factorization of their least common multiple is, in fact, factorization of each number. Since there are only 78498 primes below 10 6, trial division will be fast enough (unless you really are desperate for the last performance hit) and sifting primes down to 10 6 is also a matter of a few milliseconds.

If speed is of the utmost importance, then a combined trial division approach and a deterministic simple Miller-Rabin type test using a factorization method like Rollall Pollard's algorithm or elliptic curve factorization method is probably slightly faster (but the numbers are so small that the difference won't be big and you need a number type with more than 64 bits to make the strength test and factorization fast).

In factorization, you must, of course, remove the main 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 on which prime numbers need to be checked.



Your main problem, as far as I can see, is that your sieve is too slow and uses up too much space (which is partly due to its slowness). Use the sieve bit to improve cache locality, remove even numbers from the sieve, and stop checking to see if multiples of the square root should be flushed from max

. And you are allocating too much space for a simple 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 install m[prime[j]]

multiple times. Even if it is std::map

pretty fast, it is slower than installing it just once.

+2


source


FWIW, in response to an initial request for a faster way to get all prime numbers up to a million, that's a quick way to much .

For this, a window sieve from Eratosthenes on wheels is used with a wheel size of 30 and a window size set to the square root of the upper search limit (1000 for searches up to 1,000,000).



Since I am not fluent in C ++, I have coded it in C # assuming it should be easily convertible to C ++. However, even in C #, it can enumerate all primes up to 1,000,000 in 10 milliseconds. Even generating all primes up to a billion takes only 5.5 seconds, and I would imagine it 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))));
    }
}

      

+2


source


There seems to be some 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 nested loops "from the inside" and works on all numbers at the same time, one at a time.

Since it uses one stroke at a time, you can find the next prime as needed, avoiding generating 10 ^ 6 primes before starting. Since each number is reduced by its simple coefficients, the maximum number needed to test the numbers can be reduced, so it requires even less work.

Edit: It also makes the ending unambiguous and easy to test because the number converges to one when all of its factors have been found. In fact, when all numbers come down to one, it might terminate, although I haven't used this 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 resolves it directly.

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

No number can be successfully handled by a simple (one mod and one division) more than log2 (10 ^ 12) = 43 mods and divides, so 4300 divisions and 4300 mods, so simple tests prevail. To keep it simple, call it 8,000,000 whole divisions and mods. It should generate prime numbers, but as already stated by Daniel Fischer, it's fast. The rest is accounting.

So, on a modern processor, I would WAG around 1,000,000,000 divisions or mods per second, so the runtime is around 10ms x 2?

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

No cleverness, exactly as explained there.

I poorly rated my estimate at about 10x, but still only 20% of the maximum allowed execution time.

Performance (with some stamp 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 that almost all primes should be tested as this appears to be the main computation.

and also with the same print volume, but the value is 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) (point 3) Copyright (C) 2007 Free Software Foundation, Inc. This is free software; see source for copy conditions. There is no guarantee here; even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

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

Summary. It works well, and the startup time is about 20% of the allowed maximum on a relatively old processor, which is comparable to Daniel Fischer's implementation.

Q: I am a new member here, so it seems a little harsh on my answer when:
a. it appears to be accurate, complete, and satisfies all criteria, and
b. I wrote the code, tested it and presented the results.
What did I do wrong? How can I get feedback to improve the situation?

+1


source


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 a list of simple LCM ratios, not the actual LCM value (i.e. it does not calculate the scores) and I believe all of this is necessary.

0


source







All Articles