BigQueryã§ç´ æ•°ã‚„é«˜åº¦åˆæˆæ•°ã®åˆ—æŒ™ã‚’SQLで計算(実用性はない)

素数列挙

以前書いた記事では\(10^4\)ぐらいのオーダーまでしか素数を列挙できませんでしたが、クエリを工夫して\(10^7\)ぐらいまでは行けるようになりました
sucrose.hatenablog.com

以前のクエリ

上の記事に書いた試し割りによる\(10^4\)個ぐらいまでが限界のクエリ

#standardSQL

#数列を作る
WITH numbers AS (
  SELECT
    i
  FROM UNNEST(GENERATE_ARRAY(3, 20000, 2)) AS i
)
#割り切れるか試す
, divisible AS (
  SELECT
    a.i
    , MOD(a.i, b.i) = 0 AS is_divisible
  FROM numbers AS a
  JOIN numbers AS b
  ON SQRT(a.i) + 1 > b.i
)

#1回も割り切れていないものが素数
SELECT * FROM UNNEST([2, 3]) AS i
UNION ALL
SELECT
  i
FROM divisible
GROUP BY i
HAVING NOT LOGICAL_OR(is_divisible)
ORDER BY i

JOINをやめて平方根までの割る数を生成する

割る数を用意するためにJOINしていたのですがこれをやめてUNNEST(GENERATE_ARRAY(3, CAST(SQRT(i) + 1 AS INT64), 2)にすると\(10^7\)までの素数列挙もできるようになりました。ただし10分ぐらいかかります
計算量は\(O(n \sqrt n)\)ぐらいになっているっぽいです

#standardSQL

#数列を作る
WITH numbers AS (
  SELECT
    i + 100000 * j AS i
  FROM UNNEST(GENERATE_ARRAY(1, 99999, 2)) AS i, UNNEST(GENERATE_ARRAY(0, 99)) AS j
  WHERE i + 100000 * j > 2
)
#割り切れるか試す
, primes AS (
  SELECT
    i
    , LOGICAL_OR(MOD(i, j) = 0) AS is_divisible
  FROM numbers, UNNEST(GENERATE_ARRAY(3, CAST(SQRT(i) + 1 AS INT64), 2)) AS j
  GROUP BY i
  HAVING NOT is_divisible
)

SELECT * FROM UNNEST([2, 3]) AS i
UNION ALL
SELECT
  i
FROM primes
ORDER BY i

エラトステネスのふるいに似た方法を使う

エラトステネスのふるいは素数の倍数は素数でないものとして弾いていく方法で素数かどうかの判定が他のループの結果に依存しているのでSQLで書くのは難しいです
なので、素数だけに絞るのは諦めてすべての整数の倍数について計算するSQLを書きます
このときの計算量は調和級数になるので以下のようになります
$$O(\sum_{i=1}^n \frac{n}{i} = \frac{n}{1} + \frac{n}{2} + \frac{n}{3} + \frac{n}{4} + \frac{n}{5} + \dots + \frac{n}{n}) = O(n \log n)$$
この方法だと\(10^7\)までの素数列挙が3分弱ぐらいでできます

#standardSQL

#数列を作る
WITH numbers AS (
  SELECT
    i + 100000 * j AS i
  FROM UNNEST(GENERATE_ARRAY(0, 99999)) AS i, UNNEST(GENERATE_ARRAY(0, 99)) AS j
  WHERE i + 100000 * j > 1
)
, primes AS (
  SELECT
    i * (j1 + 100000 * j2) AS x
  FROM numbers,
    UNNEST(GENERATE_ARRAY(0, MOD(CAST(FLOOR(9999999 / i) AS INT64), 100000))) AS j1,
    UNNEST(GENERATE_ARRAY(0, CAST(FLOOR(CAST(FLOOR(9999999 / i) AS INT64) / 100000) AS INT64))) AS j2
  WHERE j1 + j2 > 0
  GROUP BY x
  HAVING COUNT(*) = 1
)

SELECT
  x
FROM primes
ORDER BY x

ちょっとした高速化

ある数が素数で割り切れるときの最大の素因数は割られる数の平方根より小さいです
なのでエラトステネスのふるい風に数を列挙するときにも平方根までのかける数を考えればよいことになります
この高速化を入れると\(10^7\)までの素数列挙が1.5分ぐらいでできます

#standardSQL

#数列を作る
WITH numbers AS (
  SELECT
    i + 100000 * j AS i
  FROM UNNEST(GENERATE_ARRAY(0, 99999)) AS i, UNNEST(GENERATE_ARRAY(0, 99)) AS j
  WHERE i + 100000 * j > 1
)
, divisible AS (
  SELECT
    i * j AS x
  FROM numbers,
    UNNEST(GENERATE_ARRAY(1, CAST(FLOOR(IF(9999999 / i < SQRT(9999999) + 1, 9999999 / i, SQRT(9999999) + 1)) AS INT64))) AS j
  GROUP BY x
  HAVING COUNT(*) = 1
)

SELECT
  x
FROM divisible
ORDER BY x

今回の話のオチ

クエリの速度的には\(10^8\)個以上も行けるぞ！と思っていたのですが、これ以上数列の数を多くするとメモリが足りなくなってエラーになってしまいました、残念
出力が正しいかはあまりちゃんと検証していないので間違っていたら教えてください

おまけ: ついでに高度合成数を計算するクエリをもだしたので載せておきます

以下の記事で書きましたが高度合成数は約数の個数がそれより小さないずれかの正の整数よりも多い整数のことです
sucrose.hatenablog.com
エラトステネスのふるいに似た方法でやると\(10^7\)までの高度合成数の列挙が3分弱ぐらいでできます

#standardSQL

#数列を作る
WITH numbers AS (
  SELECT
    i + 100000 * j AS i
  FROM UNNEST(GENERATE_ARRAY(0, 99999)) AS i, UNNEST(GENERATE_ARRAY(0, 99)) AS j
  WHERE i + 100000 * j > 1
)
#約数の数を調べる
, divisible AS (
  SELECT
    i * (j1 + 100000 * j2) AS x,
    COUNT(*) AS count
  FROM numbers,
    UNNEST(GENERATE_ARRAY(0, MOD(CAST(FLOOR(9999999 / i) AS INT64), 100000))) AS j1,
    UNNEST(GENERATE_ARRAY(0, CAST(FLOOR(CAST(FLOOR(9999999 / i) AS INT64) / 100000) AS INT64))) AS j2
  WHERE j1 + j2 > 0
  GROUP BY x
)

SELECT
  1 AS x,
  1 AS count
UNION ALL
SELECT
  x,
  count + 1 AS count
FROM
(
  SELECT
    x,
    count,
    MAX(count) OVER less_than_x AS prev_max
  FROM divisible
  WINDOW less_than_x AS (ORDER BY x ROWS BETWEEN UNBOUNDED PRECEDING and 1 PRECEDING)
)
WHERE prev_max IS NULL OR count > prev_max
ORDER BY x