ベクトル演算は強力だなぁ

「ピタゴラスの定理の整数解（原始ピタゴラス数）」ですが...

元のプログラムは，私のマシンだと1秒ほどかかる。
> system.time({
+     ## ピタゴラスの定理
+     library(gmp)
+     n = 1000
+     i = 1
+     ans = matrix(0, floor(n) * 3, ncol = 3)
+
+     for (a in 1:n) {
+         for (b in 1:n) {
+             if (a <= b) {
+                 if (floor(sqrt(a * a + b * b)) == sqrt(a * a + b * b)) {
+                     temp = c(a, b, sqrt(a^2 + b^2))
+                     if (Reduce(gcd, temp) == 1) {
+                         ans[i, ] = c(a, b, sqrt(a^2 + b^2))
+                         i = i + 1
+                     }
+                 }
+             }
+         }
+     }
+
+     ans = ans[apply(ans, 1, sum) != 0, ]
+ })
   ユーザ   システム       経過  
     1.004      0.007      1.003

あまりよいプログラムではないが，ベクトル演算を使うと3.7倍ほど速くなる。
> system.time({
+     library(gmp)
+     n = 1000
+     d = expand.grid(1:n, 1:n)
+     s = sqrt(rowSums(d^2))
+     d = cbind(d, s)
+     d = d[floor(s) == s & d[, 1] <= d[, 2], ]
+     GCD = apply(d, 1, function(x) Reduce(gcd, x))
+     d = unique(d/GCD)
+     d = d[order(d[, 1]), ]
+ })
   ユーザ   システム       経過  
     0.280      0.010      0.301