CORDICアルゴリズムに必要なarctan(2^{-i})のテーブルを整数演算で作成する

前回の記事ではCORDICアルゴリズムで三角関数を求めたのですが，定数として$\arctan 2^{-i}$の表が必要になりますので，本記事では$\arctan 2^{-i}$の表の計算をElixirで記述しました．

ソースコードは下記です．

実装にあたり，下記の公式集を参考にしました．

実装の方針

下記のTaylor展開式を用います．

$$\arctan\frac{1}{n} = \sum_{k=0}^{\infty}\frac{(-1)^k}{2k+1}\frac{1}{n^{2k+1}}=\frac{1}{n}-\frac{1}{3}\frac{1}{n^3}+\frac{1}{5}\frac{1}{n^5}-\frac{1}{7}\frac{1}{n^7}+\cdots$$

$n$が小さい時には収束が遅いので，下記の公式を使います．

$$\arctan 1 = 12\arctan\frac{1}{49} + 32\arctan\frac{1}{57} - 5\arctan\frac{1}{239} + 12\arctan\frac{1}{110443}$$

$$\arctan\frac{1}{2} = \arctan\frac{1}{3} + \arctan\frac{1}{7}$$

$$\arctan\frac{1}{3} = \arctan\frac{1}{5} + \arctan\frac{1}{8}$$

$$\arctan\frac{1}{7} = \arctan\frac{1}{8} + \arctan\frac{1}{57}$$

$$\arctan\frac{1}{57} = \arctan\frac{1}{32} - \arctan\frac{1}{73}$$

また次の式も用います．

$$\lim_{i \rightarrow 0}\arctan \frac{1}{i} = \frac{\pi}{2} = 2\arctan 1$$

$$ \lim_{i \rightarrow \infty}\arctan \frac{1}{i} = \arctan 0 = 0 $$

実装

$0 \leq x \leq 1$ の実数をnビットの固定小数点数表記として，$2^{n-2}$を乗じて表現します．

まず全体をGenServerで構成してキャッシュを形成します．

defmodule NxAtanTable do
  use GenServer

  @impl true
  def init(initial_state) do
    {:ok, initial_state}
  end

  def start_link(initial_state \\ %{}) do
    GenServer.start_link(__MODULE__, initial_state, name: __MODULE__)
  end

  @impl true
  def handle_call({:atan_of_reciprocal, n, b}, _from, state) do
    cache_atan_of_reciprocal({n, b}, state, Map.get(state, {n, b}))
  end

  def atan_of_reciprocal(n, b), do: GenServer.call(__MODULE__, {:atan_of_reciprocal, n, b})

  defp cache_atan_of_reciprocal({n, b}, state, nil) do
    r = atan_of_reciprocal_s(n, state, b)
    {:reply, r, Map.put(state, {n, b}, r)}
  end

  defp cache_atan_of_reciprocal(_, state, r) do
    {:reply, r, state}
  end

表の生成は次のようにStream.unfoldとStream.map，Enum.reverseを用います．

  def table(n, b) do
    Stream.unfold(n, fn
      0 -> nil
      n -> {n - 1, n - 1}
    end)
    |> Stream.map(fn n -> Bitwise.bsl(1, n) end)
    |> Stream.map(&atan_of_reciprocal(&1, b))
    |> Enum.reverse()
    |> Nx.tensor(type: {:s, b})
  end

前述のTaylor展開式は次のようにStream.unfoldとEnum.reduceで実装します．

  defp atan_of_reciprocal_s(n, _state, bit) when n > 0 and n < Bitwise.bsl(1, bit - 1) do
    n2 = n * n

    Stream.unfold({0, 0, n, Bitwise.bsl(1, bit - 2), bit, n2}, &atan_of_reciprocal_body/1)
    |> Enum.reduce(fn
      {_, a, _, _, _, _}, _acc -> a
    end)
  end

  defp atan_of_reciprocal_body({k, a, b, c, bit, n2}) do
    if c > 0 do
      c = div(Bitwise.bsl(1, bit - 2), b * (Bitwise.bsl(k, 1) + 1))
      b = b * n2

      a =
        case Bitwise.band(k, 1) do
          0 -> a + c
          1 -> a - c
        end

      {
        {k + 1, a, b, c, bit, n2},
        {k + 1, a, b, c, bit, n2}
      }
    end
  end

あとは前述の公式に沿って，関数パターンマッチでatan_of_reciprocal_sを実装していきます．

  defp atan_of_reciprocal_s(0, state, b) do
    {:reply, r1, _state} = cache_atan_of_reciprocal({1, b}, state, Map.get(state, {1, b}))
    Bitwise.bsl(r1, 1)
  end

  defp atan_of_reciprocal_s(1, state, b) do
    {:reply, r1, state} = cache_atan_of_reciprocal({49, b}, state, Map.get(state, {49, b}))
    {:reply, r2, state} = cache_atan_of_reciprocal({57, b}, state, Map.get(state, {57, b}))
    {:reply, r3, state} = cache_atan_of_reciprocal({239, b}, state, Map.get(state, {239, b}))

    {:reply, r4, _state} =
      cache_atan_of_reciprocal({110_443, b}, state, Map.get(state, {110_443, b}))

    12 * r1 + 32 * r2 - 5 * r3 + 12 * r4
  end

  defp atan_of_reciprocal_s(2, state, b) do
    {:reply, r1, state} = cache_atan_of_reciprocal({3, b}, state, Map.get(state, {3, b}))
    {:reply, r2, _state} = cache_atan_of_reciprocal({7, b}, state, Map.get(state, {7, b}))
    r1 + r2
  end

  defp atan_of_reciprocal_s(3, state, b) do
    {:reply, r1, state} = cache_atan_of_reciprocal({5, b}, state, Map.get(state, {5, b}))
    {:reply, r2, _state} = cache_atan_of_reciprocal({8, b}, state, Map.get(state, {8, b}))
    r1 + r2
  end

  defp atan_of_reciprocal_s(7, state, b) do
    {:reply, r1, state} = cache_atan_of_reciprocal({8, b}, state, Map.get(state, {8, b}))
    {:reply, r2, _state} = cache_atan_of_reciprocal({57, b}, state, Map.get(state, {57, b}))
    r1 + r2
  end

  defp atan_of_reciprocal_s(57, state, b) do
    {:reply, r1, state} = cache_atan_of_reciprocal({32, b}, state, Map.get(state, {32, b}))
    {:reply, r2, _state} = cache_atan_of_reciprocal({73, b}, state, Map.get(state, {73, b}))
    r1 - r2
  end

  defp atan_of_reciprocal_s(n, _state, bit) when n > 0 and n < Bitwise.bsl(1, bit - 1) do
    n2 = n * n

    Stream.unfold({0, 0, n, Bitwise.bsl(1, bit - 2), bit, n2}, &atan_of_reciprocal_body/1)
    |> Enum.reduce(fn
      {_, a, _, _, _, _}, _acc -> a
    end)
  end

  defp atan_of_reciprocal_s(n, _state, bit) when n >= Bitwise.bsl(1, bit - 1) do
    0
  end

最初に実装する時には浮動小数点数表記で動くものを作ってから固定小数点数表記に少しずつ改めていきました．