Как вычислить все простые множители lcm из n целых чисел? - PullRequest
Новая
       17

Как вычислить все простые множители lcm из n целых чисел?

4 голосов
/ 17 марта 2012

У меня есть n целых чисел, хранящихся в массиве a, скажем, a [0], a [1], ....., a [n-1], где каждый a[i] <= 10^12 и n <100.Теперь мне нужно найти все простые множители LCM для этих n целых чисел, т. Е. LCM для {a [0], a [1], ....., a [n-1]}

У меня есть метод, но мне нужен более эффективный.

Мой метод: 

 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

Есть ли лучший подход к этой проблеме?

Я пишуссылка на проблему:

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

Ссылка на мой код: http://pastebin.com/R8TMYxNz

Решение:

По предложению Даниэля Фишера мой кодпотребовались некоторые оптимизации, такие как более быстрое сито и некоторые незначительные модификации.После всех этих модификаций я смог решить проблему.Это мой принятый код на SPOJ, который занял 1,05 секунды: 

#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;
}

В случае, если кто-то сможет сделать это лучше или более быстрым способом, пожалуйста, дайте мне знать.

Ответы [ 4 ]

2 голосов
/ 19 марта 2012

FWIW, в ответ на первоначальный запрос о более быстром способе получить все простые числа до одного миллиона, вот намного более быстрый способ сделать это.

При этом используется оконное сито эратосфена на колесе с размером колеса 30 и размером окна, равным квадратному корню из верхнего предела поиска (1000 для поиска до 1 000 000).

Поскольку я не обладаю достаточными знаниями в C ++, я кодировал его в C #, исходя из предположения, что он должен быть легко преобразован в C ++.Однако даже в C # он может перечислять все простые числа до 1 000 000 за 10 миллисекунд.Даже генерация всех простых чисел до миллиарда занимает всего 5,3 секунды, и я думаю, что это будет еще быстрее в 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 голосов
/ 17 марта 2012

С этими ограничениями, несколькими не слишком большими числами, лучший способ найти первичную факторизацию их наименьшего общего множителя - это факторизация каждого числа. Поскольку ниже 10 6 есть только 78498 простых чисел, пробное деление будет достаточно быстрым (если только вы действительно не отчаянно нуждаетесь в последнем падении производительности) и просеивает простые числа до 10 6 . это также вопрос нескольких миллисекунд.

Если скорость имеет первостепенное значение, комбинированный подход деления проб и детерминистического простого теста типа Миллера-Рабина с методом факторизации, таким как алгоритм Ро Полларда или метод факторизации эллиптической кривой, возможно, немного быстрее (но цифры настолько малы, что разница не будет большой, и вам понадобится числовой тип с более чем 64 битами, чтобы тест на простоту и факторизацию выполнялся быстро).

При факторизации вы, конечно, должны удалить основные факторы, так как они найдены 

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
}

чтобы уменьшить предел, до которого нужно проверять простые числа.

Насколько я понимаю, ваша главная проблема заключается в том, что ваше сито слишком медленное и использует слишком много места (что частично отвечает за его медлительность). Используйте немного сита, чтобы получить лучшую локальность кэша, удалите четные числа из сита и прекратите проверять, нужно ли вычеркивать кратные числа в квадратном корне из max. И вы выделяете слишком много места для простого массива. 

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

Вы должны проверить j < size, прежде чем получить доступ к prime[j]. 

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

Не устанавливать m[prime[j]] несколько раз. Даже если std::map с довольно быстро, это медленнее, чем установка только один раз.

1 голос
/ 17 марта 2012

Кажется, есть несколько полезных алгоритмов на http://en.wikipedia.org/wiki/Least_common_multiple

В частности http://en.wikipedia.org/wiki/Least_common_multiple#A_method_using_a_table
кажется уместным.

Он выворачивает вложенные циклы "наизнанку" и работает со всеми числами одновременно, по одному простому за раз.

Поскольку он использует одно простое число за раз, вы можете найти следующее простое число по мере необходимости, избегая генерации 10 ^ 6 простых чисел перед началом. Поскольку каждое число уменьшается на свои основные факторы, максимальное простое число, необходимое для проверки чисел, может быть уменьшено, поэтому ему требуется еще меньше работы.

Изменить: Это также делает завершение однозначным и легко проверить, потому что число уменьшается до единицы, когда все его факторы были найдены. Фактически, когда все числа сводятся к одному, оно может завершаться, хотя я не использовал это свойство в своем коде.

Редактировать: Я прочитал проблему, алгоритм на http://en.wikipedia.org/wiki/Least_common_multiple#A_method_using_a_table решает ее напрямую.

SWAG: Есть 78 498 простых чисел ниже 10 ^ 6 (http://primes.utm.edu/howmany.shtml#table) В худшем случае, есть 100 чисел, которые необходимо проверить на 78 498 простых чисел, = 7,849,800 «мод» операций

Ни одно число не может быть успешно учтено простым числом (один мод и одно деление) больше, чем log2 (10 ^ 12) = 43 мода и деления, поэтому 4300 делений и 4300 мод, в которых преобладают простые тесты. Для простоты назовем это 8 000 000 целочисленных делителей и модов. Он должен генерировать простые числа, но, как уже сказал Дэниел Фишер, это быстро. Остальное - бухгалтерский учет.

Итак, на современном процессоре я бы выиграл около 1 000 000 000 делений или модов в секунду, поэтому время работы около 10 мс х 2?

Edit: Я использовал алгоритм на http://en.wikipedia.org/wiki/Least_common_multiple#A_method_using_a_table

Смарта нет, именно так, как там объяснено.

Я сильно ошибся в оценке, примерно в 10 раз, но все еще только 20% от максимально допустимого времени выполнения.

Производительность (с печатью для подтверждения результатов) 

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

на 100 номеров: 

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

10 раз, чтобы убедиться, что почти все простые числа должны быть проверены, так как это кажется основным вычислением.

, а также с тем же объемом печати, но значением почти 10 ^ 12 

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

на 100 999962000357L, // ((long)999979L)*((long)999983L)

gcc - версия i686-apple-darwin10-gcc-4.2.1 (GCC) 4.2.1 (Apple Inc., сборка 5666) (точка 3) Copyright (C) 2007 Free Software Foundation, Inc. Это бесплатное программное обеспечение; см. источник для условий копирования. Здесь нет гарантия; даже не для ИЗДЕЛИИ или ФИТНЕСА ДЛЯ ОСОБЕННОЙ ЦЕЛИ.

Название модели: MacBook Pro Название процессора: Intel Core 2 Duo Скорость процессора: 2,16 ГГц

Резюме: он четко работает, и время выполнения составляет около 20% от допустимого максимума на относительно старом процессоре, что сопоставимо с реализацией Даниэля Фишера.

Вопрос: Я новый участник, поэтому мне кажется немного резким -2 мой ответ, когда:
а. кажется точным, полным и удовлетворяющим всем критериям, и
б. Я написал код, проверил его и предоставил результаты
Что я сделал не так? Как получить обратную связь, чтобы я мог улучшить?

0 голосов
/ 18 марта 2012 
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

Примечание: здесь приводится только список основных факторов LCM, а не фактическое значение LCM (т. Е. Не рассчитывается экспонента), что, как я полагаю, является всем необходимым вопросом.

...