n個の整数の1cmのすべての素因数を計算するには？

Question

私は、配列aにa [0]、a [1]、.....、a [n-1]のような整数が格納されていて、それぞれa[i] <= 10^12とn <100です。今、これらの整数のLCMのすべての素因数、すなわち{a [0]、a [1]、...、a [n-1]のLCMを見つける必要がある。

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

ソリューションより速いふるいといくつかの小さな修正のように。これらすべての変更を行った後、私は問題を解決することができます。これは、1.05秒かかったSPOJの私の受け入れコードです：ケースで

#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

これらの制約条件では、数があまりにも多くない場合、最小公倍数の素因数分解を見つける最良の方法は、実際には各数値の因数分解です。 1の下には78498個の素数しかないので、試行の分割は（実際に最後のパフォーマンス低下に必死でない限り）十分速くなり、素数を1にふるい分けることは数ミリ秒の問題です。

速度が最も重要である場合、試行分割と、Pollardのrhoアルゴリズムや楕円曲線因子分解法などの因子分解法を使用した決定論的Miller-Rabin型素数検定の組み合わせアプローチはおそらく少し速いです数は非常に小さく、差は大きくならないので、素数性テストと素因数分解を高速化するには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++){ 

あなたはprime[j]にアクセスする前にj < sizeをチェックする必要があります。

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

m[prime[j]]を複数回設定しないでください。 std::mapが非常に高速であっても、1回だけ設定するよりも時間がかかります。

Answer 2

特定http://en.wikipedia.org/wiki/Least_common_multiple#A_method_using_a_table
でhttp://en.wikipedia.org/wiki/Least_common_multiple

でいくつかの有用なアルゴリズムがあるように思わ適切であると思われます。

ネストされたループを「内側に」折り返し、すべての数字に対して同時に機能します。一度に1つずつプライムします。

一度に1つの素数を使用するため、必要に応じて次の素数を見つけることができます。開始する前に10^6素数を生成しないでください。それぞれの数がその素因数だけ減少するにつれて、数をテストするために必要な最大の素数は減少する可能性があるため、作業がさらに少なくてすみます。

編集：すべての要因が検出されたときに番号が1に減少するため、終了を明白に簡単に確認できます。実際には、すべての数値が1に減少すると終了することができますが、コード内でそのプロパティを使用しませんでした。

編集：私は問題を読んで、http://en.wikipedia.org/wiki/Least_common_multiple#A_method_using_a_tableのアルゴリズムはそれを直接解決します。

SWAG： 10^6（http://primes.utm.edu/howmany.shtml#table） だから、最悪の場合、以下の78498個の素数がありますが、78498個の素数のためにテストする必要がある100個の番号があり、 = 7849800「MOD」操作

LOG2（10^12）= 43台の改造及び分割以上（1つのMOD一分割）プライムによって正常因数分解することができません数は、そう4300の分割及び4300の改造、これによって支配ありませんプライムテスト。 それを簡単に保つために、8,000,000の整数除算および改造と呼ぶことができます。 それは素数を生成する必要がありますが、Daniel Fischerがすでに述べたように素早くです。 残りは本棚です。

最新のプロセッサでは、約10億分の1秒または1秒あたりの改造時間があるので、実行時間は10ms×2です。

編集： が説明したように、私はちょうどまさに、http://en.wikipedia.org/wiki/Least_common_multiple#A_method_using_a_table

ない知性でアルゴリズムを使用していました。

私の推定では約10倍でしたが、実行可能な最大実行時間の20％に過ぎませんでした。それがあるように思わとしてほぼすべての素数をテストする必要がありました

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

10回、 を確保するために：

パフォーマンス

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

100の番号の

（一部印刷して結果を確認するため）主要な計算。

と、印刷の同じ量が、100のためのほぼ10^12

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

の値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フリーソフトウェア財団、 これはフリーソフトウェアです。コピー条件についてはソースを参照してください。いいえ 保証はありません。商品性や特定の目的への適合性についてさえも、

モデル名：MacBook Proの プロセッサ名：インテルCore 2 Duoプロセッサ プロセッサ速度：2.16 GHzの

概要：それは明らかに動作し、実行時間が比較的古いプロセッサ上で、許容最大の約20％であり、これはDaniel Fischerの実装に匹敵するものです。

質問：私はここに新たに寄稿しているので、-2の回答には少し厳しいようです。
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、LCMないの実際の値の素因数（すなわち指数を計算しない）、のリストを与えます。

Answer 4

FWIWは、百万までの素数のすべてを取得するためのより高速な方法のため、元の要求に応じて、ここでそれを行うために多くのより高速な方法です。

これは、ホイールのサイズが30で、ウィンドウのサイズが検索の上限（1000,000までの検索では1000）の平方根に設定された、エラトステネスのホイールベースのふるいを使用しています。

私はC++に堪能ではないので、これをC++に簡単に変換できると仮定してC＃でコード化しました。しかし、C＃でも約10ミリ秒で100万までのすべての素数を列挙できます。 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)))); 
    } 
}