galois

defmodule Polynom do
  defstruct [:coeffs, :modulo]

  defp zip(lhs, rhs) when length(lhs) == length(rhs), do: Enum.zip([lhs, rhs])

  defp zip(lhs, rhs) do
    pad_len = abs(length(lhs) - length(rhs))
    pad = for _ <- 1..pad_len, do: 0
    [lhs, rhs] = if length(lhs) < length(rhs), do: [lhs ++ pad, rhs], else: [lhs, rhs ++ pad]
    Enum.zip([lhs, rhs])
  end

  def add(lhs, rhs),
    # component wise addition
    do:
      (for {l, r} <- zip(lhs, rhs), do: l + r)
      |> reduce()

  def sub(lhs, rhs),
    # component wise subtraction
    do:
      (for {l, r} <- zip(lhs, rhs), do: l - r)
      |> reduce()

  def mul(lhs, rhs) do
    # cartesian product of all coeffs with their exponent
    for {l, le} <- Enum.with_index(lhs),
        {r, re} <- Enum.with_index(rhs) do
      {{l, le}, {r, re}}
    end
    # build map containing "exponent => [coefficients]"
    |> Enum.reduce(%{}, fn {{l, le}, {r, re}}, acc ->
      addends = [l * r | Map.get(acc, le + re, [])]
      Map.put(acc, le + re, addends)
    end)
    |> Map.to_list()
    # sort terms ascending, so we can insert the coeffs in the next step
    # --> polynom is built from high exp to low exp
    |> Enum.sort(fn {le, _l}, {re, _r} -> le > re end)
    |> Enum.reduce([], fn {_exp, addends}, coeffs ->
      [Enum.sum(addends) | coeffs]
    end)
  end

  # division by zero
  def div(poly, [0]) when is_list(poly), do: raise(ArgumentError, "divisor is zero")

  # cannot divide, rhs is too big
  def div(lhs, rhs) when length(lhs) < length(rhs), do: {[0], lhs}
  def div(lhs, rhs) when length(lhs) == length(rhs) and hd(lhs) < hd(rhs), do: {[0], lhs}

  # same polynom
  def div(poly, poly), do: {[1], [0]}

  # regular main entry
  #   polynoms are sorted from x^0 to x^n, we need to reverse,
  #   so we can easily get the coefficient with the highest exponent with pattern matching
  def div(lhs, rhs) do
    {q, rest} = div(lhs |> Enum.reverse(), rhs |> Enum.reverse(), [])
    q = q |> Enum.reverse() |> Enum.drop_while(&(&1 == 0)) |> Enum.reverse()

    rest = rest |> Enum.reverse() |> Enum.drop_while(&(&1 == 0)) |> Enum.reverse()
    rest = if rest == [], do: [0], else: rest

    {q, rest}
  end

  # no remainder -> division without rest
  defp div([], _rhs, result), do: {result, [0]}
  # remainder left and irr poly has higher exp -> division with rest
  defp div(lhs, rhs, result) when length(lhs) < length(rhs),
    do: {result, lhs |> Enum.reverse()}
  # division by single zero-exp coefficient
  defp div([l], [r], result),
    do: {[Kernel.div(l, r) | result], [Kernel.rem(l, r)]}
  # division in progress
  defp div([l | ltl], [r | rtl] = rhs, result) do
    # divide highest power coeffs
    quotient = Kernel.div(l, r)

    # multiply rhs with quotient and subtract element-wise from lhs
    #   we do not want to pad this zip, only pair up as many coeffs, as the divisor has
    head =
      Enum.zip(ltl, rtl)
      |> Enum.reduce([], fn {l, r}, acc ->
        [l - quotient * r | acc]
      end)
      |> Enum.reverse()

    # build remaining lhs
    #   drop outdated coeffs
    tail = Enum.drop(ltl, length(head))

    div(head ++ tail, rhs, [trunc(quotient) | result])
  end

  def poly_mod([], _irr_pol), do: []
  def poly_mod(_pol, []), do: []

  def poly_mod(pol, irr_pol) do

    [ll, rl] = [pol, irr_pol] |> Enum.map(&length(&1))

    cond do
      ll < rl -> pol

      ll >= rl ->
        {_p, remainder} = Polynom.div(pol, irr_pol)
        if length(remainder) < length(irr_pol) do
          remainder
        else
          [e | _] = remainder |> Enum.reverse()
          [ie | _] = irr_pol |> Enum.reverse()
          difference = abs(Kernel.div(e, ie))
          sub(remainder, mul(irr_pol, [difference]))
        end
    end
  end

  def int_mod(pol, modulus) do
    for c <- pol do
      if c < 0,
        do: rem(c + modulus, modulus),
        else: rem(c, modulus)
    end
  end

  def print(poly) do
    term =
      poly
      |> Enum.with_index()
      # drop zeros
      |> Enum.filter(fn {c, _e} -> c > 0 end)
      |> Enum.map(fn {c, e} ->
        cond do
          e == 0 -> "#{c}"
          e == 1 and c == 1 -> "x"
          e == 1 -> "#{c}x"
          c == 1 -> "x" <> exponents(e)
        end
      end)
      |> Enum.join(" + ")

    if term == "", do: "0", else: term
  end

  defp exponents(e), do: Enum.at(~w"⁰ ¹ ² ³ ⁴ ⁵ ⁶ ⁷ ⁸ ⁹", e)

  defp reduce(polynom) do
    dropped = polynom
      |> Enum.reverse()
      |> Enum.drop_while(&(&1 == 0))
      |> Enum.reverse()
    if dropped == [], do: [0], else: dropped
  end
end