cats-collection の PairingHeap を使った素数の無限リスト

Cats 関連モジュールの cats-collection に入ってる、コレクションクラス `PairingHeap` を使って、素数の無限リストを書いてみる。
(サンプルコードを Scala 3 に直した。2022-07-27)
## cats-colection について
[Cats の ドキュメント](https://typelevel.org/cats/)から "Data structures which facilitate pure functional programming"として参照されている [cats-collections](https://github.com/typelevel/cats-collections) というライブラリがある。

関数型なデータ構造の Scala 実装を集めたもので、内容は以下のものに大別される。

- Scala 標準ライブラリのデータ構造クラスを、FP 的に再実装したもの。Set(AvlSet)、Map(AvlMap)、Range など。
- Scala 標準ライブラリに対応物がないが、関数型データ構造として知られているもの。Discrete、TreeList など[^to-be-implemented]。

この記事では、後者の `PairingHeap` を試してみる。

## PairingHeap について
Chris Okasaki の 『[Purely Functional Data Structures](https://www.amazon.co.jp/Purely-Functional-Data-Structures-Okasaki/dp/0521663504/)』で、「5.5 Pairing Heap」として紹介されている[ヒープ](https://ja.wikipedia.org/wiki/%E3%83%92%E3%83%BC%E3%83%97)の一種で、シンプルな実装ながらなぜかパフォーマンスが良いので理論家を長年悩ませていたというもの。

cats-collection の中では比較的新しく、[公式ドキュメント](https://typelevel.org/cats-collections//index.html)にはまだ記述がない。Cats との関連では、`CommutativeMonoid` のインスタンスと `Order` のインスタンスが提供されている。また`PairingHeap[A]` を操作するために、`Order[A]`を前提とする。

ヒープなので、コレクション内のその時点の最小値が、常に$O(1)$で取り出せる形で配置されているが、これを利用すると素数のストリームを生成するプログラムの中で、除去すべき合成数を小さい順に取り出すことができるデータ構造として利用できる。

以下、具体的に見てみる。

## Primes サンプル
素数生成の基本といえばエラトステネスの篩になるが、典型的な実装は $N<10,000$ のように事前に定めた範囲の書き換え可能なメモリ領域を用意して、小さい素数から順にその倍数を消していって、残った数を素数の集まりとするものだった。

ここでは `PairingHeap` を使って、範囲を限定しない遅延評価リスト、ないしは無限リストを生成するプログラムを書いてみる。実装スタイルの面でも、自然とイミュータブルで関数型らしいコードになる。

アルゴリズムとしては、Wheel Sieve と PairingHeap を併用する。

Wheel Seive は、素数の候補を得るために使う。具体的には、既知の少数個の素数列から得た数値の間隔の列[^number_interval]＝*Wheel（車輪）* を、いわば「回転」させるような要領で、繰り返し足し込んで素数候補の無限リスト(以下、Wheelストリーム)を生成し、この予め絞り込んだ素数候補ストリームを篩にかける。つまり合成数を判別して除外していく。[^about-wheel]


PairingHeap は、合成数の置き場所として使う。ヒープ内の個々の要素は、その時点までに得られた素数を Wheel ストリーム中の各素数候補にかけた合成数の列で、それぞれが小さい順に並んだ`LazyList[Long]` 型の無限リストになる。この無限リストの大小比較を先頭要素同士の比較と定義しておけば、最も小さい先頭要素を持つ LazyList が PairingHeap 内で待機している状態となり、従って最小の合成数が順々に得られるようになる。

こうして Wheel ストリームの先頭要素と、PairingHeap から取り出す最小の合成数を順次比較することで、素数の無限リストを得ることができるようになる。

```scala
type Longs = LazyList[Long]

def spin(start: Long, wheel: Longs): Longs =
  def cycle(ns: Longs): Longs = ns.lazyAppendedAll(cycle(ns))
  cycle(wheel).scan(start)(_ + _)

def toComposites(ns: Longs): Longs =
  def loop(m: Long, ms: Longs): Longs = (m * ms.head) #:: loop(m, ms.tail)
  loop(ns.head, ns)

def sieve(candidates: Longs, composites: PairingHeap[Longs]): Longs =
  val n #:: ns = candidates 
  composites.minimumOption.fold {
    n #:: sieve(ns, PairingHeap(toComposites(candidates)))
  } { case m #:: ms =>
      if   n < m
      then n #:: sieve(ns, composites + toComposites(candidates))
      else sieve(if m == n then ns else candidates, composites.remove + ms)
  }
def primes: Longs =
  val startingPrimes = LazyList(2L, 3L, 5L)
  val wheel          = LazyList(4L, 2L, 4L, 2L, 4L, 6L, 2L, 6L)
  sieve(startingPrimes.lazyAppendedAll(spin(7, wheel)), empty)
```
本当は、サイズ 6くらいの初期素数列が速いらしいが、これでも 100万番目の素数 15485863 が、手元のマシンでだいたい16秒くらいで得られた。

- バージョン等: Scala 3.1.x, cats 2.8, cats-collection 0.8.0, [その他](https://github.com/yasuabe/qiita-cats3-sample/blob/main/build.sbt)
- [ソース](https://github.com/yasuabe/qiita-cats3-sample/blob/main/src/main/scala/prime_number_pairing_heap/prime_number_pairing_heap.worksheet.sc)

## リンク
- [cats-collection ドキュメント](https://typelevel.org/cats-collections//index.html)
- [cats-collection GitHub](https://github.com/typelevel/cats-collections)

[^to-be-implemented]: [TODO.md](https://github.com/typelevel/cats-collections/blob/master/TODO.md) に、Zipper、FingerTree もあるので、いつか実装されるかもしれない。

[^number_interval]: 例えば [2, 3, 5] なら 5 〜 5+2*3=11 の範囲内の 2, 3 の倍数でない数[5, 7, 11]の間隔[2, 4]を得る。これを5 に繰り返し足していくと 5,7,11,13,17,19,23,25,29〜といった数列が得られる。5以上の素数からなる合成数は含まれるが、候補はかなり絞られる。
[^about-wheel]: このサンプル実装では、Wheel を [2, 3, 5, 7] に固定したが、指定サイズでプログラム的に生成することももちろんできる。(参照: [Haskell の Data.Numbers.Primes](http://hackage.haskell.org/package/primes-0.2.1.0/docs/src/Data-Numbers-Primes.html#wheel))