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

Question

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

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

Мой метод:

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

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

Answer 1

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

Если скорость имеет первостепенное значение, то комбинированный подход пробного деления и детерминированный первичный тест типа Миллера-Рабина с использованием метода факторизации, такого как ро-алгоритм Полларда или метод факторизации эллиптической кривой, вероятно, немного быстрее (но числа настолько малы, что разница не будет большой, и вам нужен тип номера с более чем 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 s довольно быстр, это медленнее, чем установка его только один раз.

Answer 2

Там, как представляется, несколько полезных алгоритмов на 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 модов и водоразделов, так 4,300 водоразделы и 4,300 модов, так не доминирует основные тесты. Чтобы сделать его простым, назовем его 8 000 000 целых делений и модов. Ему нужно генерировать простые числа, но, как уже было сказано Даниэлем Фишером, это быстро. Остальное - бухгалтерский учет.

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

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

Нет СМАРТС, в точности, как описано здесь.

Я был плохо на моей оценке, около 10x, но все же только 20% от максимально допустимого времени выполнения.

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

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 --version i686-яблочно-darwin10-GCC-4.2 .1 (GCC) 4.2.1 (Apple Inc. build 5666) (точка 3) Copyright (C) 2007 Free Software Foundation, Inc. Это бесплатное программное обеспечение; см. источник условий копирования. Существует NO гарантия; даже для КОММЕРЧЕСКОЙ ЦЕННОСТИ или ПРИГОДНОСТИ ДЛЯ ОПРЕДЕЛЕННОЙ ЦЕЛИ.

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

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

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

Answer 3

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, а не фактическое значение МДК (т.е. не рассчитывает показатели), которые я считаю это все вопрос требуется.

Answer 4

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

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

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