雖然最終要求的是：$ x^2+y^2 ＝N $ 的整數解，過程中也回顧了 Eratosthenes's Sieve 求質數表，和用質數表做因式分解 的方法。所用的語言方面，一個目標是自己想要熟習，而且資源又常見的Java；另一個是如果自己在中學時期想這樣做的話，應該會用到的VBA。
Java的例子不少，但試過才知道即使實現的是同一套算法，效率也可以有很大差別。從前只知道應用數學中會看一套算法的Big O去比較時間複雜度，但原來實踐起來用什麼語言、什麼變數物件，結果用什麼什麼形式去表現⋯⋯都會大大影響速度。就是為了找出實際運行得更快的寫法，就一頭栽到這個質數的做法上。雖然一開始時用LinkedList、Map等寫法真的很方便理解，但之後還是儘量換成基本的陣列。
下面這樣的Java寫的 Eratosthenes Sieve 求 一億（$10^8$）之內的質數，大約是10秒。如果不計最後一個for loop用來製作傳回值的LinkedList<Integer>的話，主流程大約4秒。不過要求更大的質數，應該還要考慮該程式對使用整數的上限（int: [-2,147,483,648, 2,147,483,647]），大數字在記憶體的儲存方法，整數表的儲存方法等。很多事情都是這樣吧，開始時只要求做到是容易的，但要追求深究下去就會困難⋯⋯
最後還有兩個方法是Modified Eratosthenes Sieve 和 Euler Sieve，都是理論上應該更有效率的算法，不過實踐起來還是達不到應有的分別，要怪我對如何更好地寫代碼還是不熟識吧。 Eratosthenes Sieve的Modification的想法是一開始就不考慮2和3的倍數，只看6k+1和6k+5的情況。至於 Euler Sieve 是改善Eratosthenes Sieve中一個合成數會重覆被質數篩去而浪費的時間，例如，合成數6 會在篩選質數 2和3 的倍數時重覆考慮。理想中Eratosthenes Sieve的時間複雜度是 O( n*log(log(n)) )，而Euler Sieve的時間複雜度是O(n)。實測中，同樣求一億（$10^8$）之內的質數，大約是2-3秒的時間。
JAVA：
Prime.javapackage JGT; import java.util.*; public class Prime { // 用 Eratosthenes Sieve 求質數表 public static List<Integer> Eratosthenes(int max) { boolean[] a = new boolean[max+1]; LinkedList<Integer> primes = new LinkedList<>(); for (int i=2; i<=max; i++) a[i] = true; for (int i=2; i<=Math.sqrt(max); i++) { if (a[i]) for (int j=i; j*i<=max; j++) a[j*i] = false; } for (int i = 0; i <= max; i++) { if (a[i]) primes.add(i); } return primes; }// 選擇求質數表的方法，因為之後還有其他方法public static List<Integer> create(int max) { return Eratosthenes(max);// return Euler(max);} // 因式分解（一） public static List<Integer> factorize(int n) { final List<Integer> Primes; Primes = Prime.create(n / 2); return factorize(n, Primes); } //因式分解（二），提供public static List<Integer> factorize(int n, List<Integer> primes) { List<Integer> factors = new ArrayList<>(); if (primes.size() <= Math.sqrt(n)) { primes = Prime.create(n / 2); } for(int i = 0; primes.get(i) <= Math.sqrt(n);) { if(n % primes.get(i) == 0) { factors.add(primes.get(i)); n /= primes.get(i); } else { i++; } } factors.add(n); return factors; } }的質數表預先準備
測試用的 PrimeFactorize.java，將結果顯示到Console和檔案：
package JGT; import java.util.*; import java.io.*; public class PrimeFactorize { public static void main(String[] args) { final int N = 1000000; List<Integer> primes, factors = new ArrayList<>(); long startTime, endTime; //用 Eratosthenes Sieve 求質數表startTime = System.currentTimeMillis(); primes = Prime.Eratosthenes(N); endTime = System.currentTimeMillis(); System.out.println("Execution time: " + (endTime - startTime) + " ms." ); try { FileWriter fw = new FileWriter("primes.txt"); for (Object p : primes) fw.write(p +"\n"); fw.close(); } catch (IOException e) {} // 求因式分解（一） startTime = System.currentTimeMillis(); factors = Prime.factorize(1000545); endTime = System.currentTimeMillis(); System.out.println("Execution time: " + (endTime - startTime) + " ms." ); for(Integer f : factors) { System.out.printf("%d ", f); } System.out.println(); //因式分解（二），提供startTime = System.currentTimeMillis(); factors = Prime.factorize(1000545, primes); endTime = System.currentTimeMillis(); System.out.println("Execution time: " + (endTime - startTime) + " ms." ); for(Integer f : factors) { System.out.printf("%d ", f); } System.out.println(); } }的質數表預先準備
VBA：
'Public Function Eratosthenes(max As Integer) As Variant Dim a() As Boolean Dim primes() As Integer For i = 2 To max ReDim Preserve a(1 To i) a(i) = True Next i For i = 2 To Sqr(max) For j = i * i To max Step i a(j) = False Next j Next i For i = 2 To max If a(i) Then primeCount = primeCount + 1 ReDim Preserve primes(1 To primeCount) primes(primeCount) = i End If Next i Eratosthenes = primes End Function '用 Eratosthenes Sieve 求質數表求因式分解Public Function Factorize(n As Integer) As Variant Dim factors() As Integer primes = Eratosthenes(n / 2) factorCount = 0 For Each i In primes If (n Mod i = 0) Then Do While (n Mod i = 0) factorCount = factorCount + 1 ReDim Preserve factors(1 To factorCount) factors(factorCount) = i n = n / i Loop End If Next i Factorize = factors End Function
測試用的 Sub-procedure，將結果顯示到Column A, B：
Sub Testing()
primes = Eratosthenes(100)
For i = 1 To (UBound(primes) - LBound(primes) + 1)
Range("A" & i).Value = primes(i)
Next i
factors = Factorize(100)
For i = 1 To (UBound(factors) - LBound(factors) + 1)
Range("B" & i).Value = factors(i)
Next i
End Sub
其他：
另外兩套求質數表的做法。當中Modified Eratosthenes Sieve 是按 這篇 寫的，而Euler Sieve是按 這篇 的pseudo-code寫的。
Euler Sieve中留意 prime[j] < i 的關係、和這句：if (i % primeList[j] == 0) break; 讓合成數只會被它最小的一個質因數篩走。
當 i 可以整除 質數 primeList[j] 的時候，會跳出最內層的 for loop 而不必浪費時間去篩走 i 的倍數，因為它們應該在篩走 質數 primeList[j] 的倍數時已經處理了。如果用藍色表示會被篩走的合成數當時的(i * prime[j])，刪除線表示因為跳出 for-loop 而去略去的組合，大概是這樣？
i : prime[j]
1:
2: 2
3: 2,
4:
5: 2, 3,
6:
7: 2, 3, 5,
8:
9: 2,
10:
11: 2, 3, 5, 7,
...
15: 2,
...
//public static List<Integer> Eratosthenes2(int max) { LinkedList<Integer> primes = new LinkedList<>(); boolean[] a = new boolean[max+1]; a[2]=true; a[3]=true; a[5]=true; for(int i=1; i<=(max-5)/6; i++) { a[6*i+1] = true; a[6*i+5] = true; } if( (max-5)/6+2 <= max) a[(max-5)/6+2]=true; for (int i=2; i<=Math.sqrt(max); i++) { if (a[i]) for (int j=i; j*i<=max; j++) a[j*i] = false; } for (int i=0; i<=max; i++) { if (a[i]) primes.add(i); } return primes; } //用 修改的 Eratosthenes Sieve 求質數表public static List<Integer> Euler(int max) { boolean[] a = new boolean[max+1]; int[] primeList = new int[(int)( max/(Math.log(max)-4) )]; int primeCount = 0; LinkedList<Integer> primes = new LinkedList<>(); for (int i=2; i<=max; i++) a[i] = true; for (int i=1; i<=max; i++) { if (a[i]) primeList[primeCount++] = i; for (int j = 0; j < primeCount; j++) { if (i * primeList[j] > max) break; a[i * primeList[j] ] = false; if (i % primeList[j] == 0) break; } } for (int i=0; i<=max; i++) { if (a[i]) primes.add(i); } return primes; }用 Euler Sieve 求質數表
沒有留言:
張貼留言