Elixir (とJavaScriptとRuby) でMöbius関数を実装した

Elixirでメビウス函數 (Mobius function) を実装する積りだったのだけれど、余りにもElixirを書けな過ぎて、JavaScriptとRubyでも実装した。

ElixirでMöbius函數

素數列 (素数列) は、エラトステネスの篩 (Sieve of Eratosthenes) で求めた。遅延listが無いので、其れの実装に一番手間を掛けた。普通にlistを使った遅延list (defmodule Sequence) はまずまず出来たが、EVM (Erlang Virtual Machine) processを使った実装 (defmodule Stream) は、ほんとに手間であった。未だ迂遠な気がしてならない。
詰まりHaskellならもっと素早く出来たと云ふ事だ。Elixir (Erlang) に不得手な事をやらせたのだから仕方が無い。
Mobius.mobius_with_sequence, Prime.PrimeSequence, Sequenceが、遅延listを使った実装である。Mobius.mobius_with_stream, Prime.PrimeStream, Streamが、EVM processを使った実装である。
遅延listの実装で知ったが、Elixir (Erlang) には型が無いので、[head | tail]とした時、tailが別にlistでなくても好い。単に何らかのheadと何らかのtailとがconsで結び附いてゐれば好いのだ。其処で、tailを、無引數で評価するとlistを返す函數を入れる事にした。Elixir (Erlang) は正格評価で副作用も有る為、無引數の函數は充分意味を持つ。函數と云ふより全てactionなのだ。
勿論tailがlistでなければ、Elem.mapや:lists.filter等のlist函數には食わせられない。よって同等の函數を自分で実装しなければならなかったが、pattern matchに支障は起らない。
と云ふ事で、Elixirで無限listが必要な時は、遅延listを使おうと思った……。外部環境と遣り取りする種の無限streamは、EVM processの方が好い気がするけど。

cf. Haskell (とElixir) でもMöbius関数を実装した http://c4se.hatenablog.com/entry/2013/07/14/185745 短くして再掲した。

# license: Public Domain

defmodule Mobius do
  @moduledoc """
  """

  @doc ""
  @spec mobius_with_sequence(integer) :: integer
  def mobius_with_sequence n do
    factors = Prime.PrimeSequence.prime_factors n
    mobius_from_factors factors
  end

  @doc ""
  @spec mobius_with_stream(integer) :: integer
  def mobius_with_stream n do
    factors = Prime.PrimeStream.prime_factors n
    mobius_from_factors factors
  end

  @spec mobius_from_factors([{integer, integer}]) :: integer
  defp mobius_from_factors factors do
    cond do
      :lists.any(fn {_, n} -> n > 1 end, factors) -> 0
      rem(length(factors), 2) == 1                -> -1
      true                                        -> 1
    end
  end
end

# {{{ Prime
defmodule Prime do
  @moduledoc """
  """

  @doc ""
  @spec prime_index(integer, integer) :: [index: integer, number: integer]
  def prime_index n, prime do
    prime_index n, prime, 0
  end
  defp prime_index(n, prime, acc) when
    rem(n, prime) != 0 do
    [index: acc, number: n]
  end
  defp prime_index n, prime, acc do
    prime_index div(n, prime), prime, acc + 1
  end

  # {{{ Prime.PrimeSequence
  defmodule PrimeSequence do
    @moduledoc """
    """

    @doc "Integer factorization."
    @spec prime_factors(integer) :: [{integer, integer}]
    def prime_factors n do
      primes = sieve
      prime_factors n, [], primes
    end
    @spec prime_factors(integer, [{integer, integer}], [integer | function() :: List]) :: [{integer, integer}]
    defp prime_factors 1, factors, _ do
      :lists.reverse factors
    end
    defp prime_factors n, factors, [prime|xs] do
      [index: index, number: n] = Prime.prime_index n, prime
      if index > 0, do: factors = [{prime, index} | factors]
      prime_factors n, factors, xs.()
    end

    @doc "Primes generator using Sieve of Eratosthenes."
    @spec sieve() :: [integer | function() :: List]
    def sieve do
      sieve Sequence.succ(2)
    end
    defp sieve [prime|xs] do
      pred = fn n -> rem(n, prime) != 0 end
      xs = Sequence.filter pred, [prime | xs]
      [prime | fn -> sieve xs end]
    end
  end
  # }}} Prime.PrimeSequence

  # {{{ Prime.PrimeStream
  defmodule PrimeStream do
    @moduledoc """
    """

    @doc "Integer factorization."
    @spec prime_factors(integer) :: [{integer, integer}]
    def prime_factors n do
      primes = sieve
      prime_factors n, [], primes
    end
    @spec prime_factors(integer, [{integer, integer}], StreamRecord) :: [{integer, integer}]
    defp prime_factors 1, factors, _ do
      :lists.reverse factors
    end
    defp prime_factors n, factors, primes do
      prime = Stream.current primes
      [index: index, number: n] = Prime.prime_index n, prime
      if index > 0, do: factors = [{prime, index} | factors]
      Stream.next primes
      primes = sieve primes
      prime_factors n, factors, primes
    end

    @doc "Primes generator using Sieve of Eratosthenes."
    @spec sieve() :: StreamRecord
    def sieve do
      sieve Stream.succ(2)
    end
    @spec sieve(StreamRecord) :: StreamRecord
    def sieve stream do
      prime = trunc Stream.current(stream)
      pred = fn
        n ->
          n = trunc n
          n == prime or rem(n, prime) != 0
      end
      stream = Stream.filter pred, stream
      func = fn acc -> [acc: sieve(acc), item: Stream.next acc] end
      Stream.start func, stream, prime
    end
  end
  # }}} Prime.PrimeStream
end
# }}} Prime

# {{{ Sequence
defmodule Sequence do
  @moduledoc """
  Lazy evaluated sequence.
  [x | xs/0]
  """

  @doc "Lazy map."
  @spec map(function(T) :: any, [T | List]) :: [T | List]
  def map func, [], do: []
  def map func, [x], do: [func.(x)]
  def map func, [x|xs] do
    [func.(x) | fn -> map xs.(), func end]
  end

  @doc "Lazy filter."
  @spec filter(function(T) :: boolean, [T | List]) :: [T | List]
  def filter pred, [], do: []
  def filter pred, [x] do
    if pred.(x), do: [x], else: []
  end
  def filter pred, [x|xs] do
    if pred.(x) do
      [x | fn -> filter pred, xs.() end]
    else
      filter pred, xs.()
    end
  end

  @doc "Lazy list of integers."
  @spec succ() :: [integer | function() :: List]
  def succ do
    succ 0
  end
  @doc "Lazy list of integers that start from n."
  @spec succ(integer) :: [integer | function() :: List]
  def succ n do
    [n | fn -> succ n + 1 end]
  end
end
# }}} Sequence

# {{{ Stream
defmodule Stream do
  @moduledoc """
  Lazy evaluated stream with an independent process.
  cf. M.Hiroi's Home Page / お気楽 Erlang プログラミング入門
      http://www.geocities.jp/m_hiroi/func/abcerl06.html
  """

  @doc ""
  defrecord StreamRecord,
    pid: nil :: pid,
    func: nil :: function(T) :: [acc: T, item: S],
    acc: nil :: T,
    item: nil :: S

  @doc ""
  @spec start(function(T) :: [acc: T, item: S], T, S) :: StreamRecord
  def start func, acc0, item0 do
    stream = StreamRecord.new func: func, acc: acc0, item: item0
    # stream.pid Process.spawn_link(Stream, function(loop/1), [stream])
    stream.pid Process.spawn_link(fn -> loop(stream) end)
  end

  @doc ""
  @spec next(StreamRecord) :: any
  def next stream do
    stream = send_message stream, :next
    receive_item stream
  end

  @doc ""
  @spec current(StreamRecord) :: any
  def current stream do
    stream = send_message stream, :current
    receive_item stream
  end

  @doc ""
  @spec stop(StreamRecord) :: any
  def stop stream do
    send_message stream, :stop
    receive_item stream
  end

  @doc ""
  @spec dropwhile(function(any) :: boolean, StreamRecord) :: StreamRecord
  def dropwhile pred, stream do
    if pred.(current stream) do
      next stream
      dropwhile pred, stream
    else
      stream
    end
  end

  @doc ""
  @spec map(function(any) :: [acc: any, item: any], StreamRecord) :: StreamRecord
  def map func, stream do
    start fn acc -> [acc: acc, item: func.(next acc)] end, stream, func.(current stream)
  end

  @doc ""
  @spec filter(function(any) :: boolean, StreamRecord) :: StreamRecord
  def filter pred, stream do
    func = fn
      acc ->
        next acc
        acc = dropwhile fn n -> not pred.(n) end, acc
        [acc: acc, item: current acc]
    end
    start func, stream, current dropwhile(fn n -> not pred.(n) end, stream)
  end

  @doc "Lazy list of integers."
  @spec succ() :: StreamRecord
  def succ do
    succ 0
  end
  @doc "Lazy list of integers that start from n."
  @spec succ(integer) :: StreamRecord
  def succ n do
    start fn n -> [acc: n + 1, item: n + 1] end, n, n
  end

  @spec loop(StreamRecord) :: any
  defp loop stream do
    StreamRecord[func: func, acc: acc, item: item] = stream
    receive do
      {:next, pid} ->
        [acc: acc, item: item] = func.(acc)
        stream = stream.update acc: acc, item: item
        pid <- {:item, item}
        loop stream
      {:current, pid} ->
        pid <- {:item, item}
        loop stream
      {:stop, pid} ->
        pid <- :ok
      _ ->
        exit :unknown_message
    end
  end

  @spec send_message(StreamRecord, atom) :: StreamRecord | {:error, String.t}
  defp send_message stream, sym do
    StreamRecord[pid: pid] = stream
    if :erlang.is_process_alive pid do
      pid <- {sym, Kernel.self}
      stream
    else
      {:error, "The process is dead."}
    end
  end

  @spec receive_item(StreamRecord | {:error, String.t}) :: any | {:error, String.t}
  defp receive_item {:error, message} do
    {:error, message}
  end
  defp receive_item _ do
    receive do
      {:item, value}     -> value
      :ok                -> :ok
      {:EXIT, _, reason} -> {:error, "EXIT: #{reason}"}
      _                  -> {:error, "Unknown message."}
    after
      1000 -> {:error, "Timeout."}
    end
  end
end
# }}} Stream

IO.puts Mobius.mobius_with_sequence(3 * 19 * 19 * 23) == 0
IO.puts Mobius.mobius_with_sequence(3 * 19 * 23) == -1
IO.puts Mobius.mobius_with_sequence(3 * 19 * 23 * 151) == 1

IO.puts Mobius.mobius_with_stream(3 * 19 * 19 * 23) == 0
IO.puts Mobius.mobius_with_stream(3 * 19 * 23) == -1
# IO.puts Mobius.mobius_with_stream(3 * 19 * 23 * 151) == 1
# vim:set et sw=2 sts=2 ff=unix foldmethod=marker:

JavaScriptでMöbius函數

JavaScriptは状態持ちなので、素數列をcacheして割り算していった。屑である。

/**
 * @license Public Domain
 */

(function(scope) {
/** @type {function(number):number} */
scope.mobius = mobius;

Object.defineProperty(Array.prototype, 'last',
  {
    get: function() { return this[this.length - 1]; }
  });

/**
 * Mobius function.
 * cf. メビウス関数 - Wikipedia http://ja.wikipedia.org/wiki/%E3%83%A1%E3%83%93%E3%82%A6%E3%82%B9%E9%96%A2%E6%95%B0
 *
 * @param {number} n
 * @return {number}
 */
function mobius(n) {
  var doseContainSquare = false,
      numberOfPrimes = 0;

  eachPrimeFactors(n,
    function(prime, index) {
      if (index > 1) {
        doseContainSquare = true;
        return false;
      }
      numberOfPrimes += 1;
    });
  if (doseContainSquare)
    return 0;
  else if (numberOfPrimes % 2 === 1)
    return -1
  else
    return 1;
}

/**
 * Call a function on each prime factors.
 *
 * @param {number} n
 * @param {function(number,number):boolean} callback
 *   Break iterator when return false.
 */
function eachPrimeFactors(n, callback) {
  /**
   * @param {number} prime
   * @return {number}
   */
  function indexOfPrime(prime) {
    var index = 0;

    while (n % prime === 0) {
      n /= prime;
      index += 1;
    }
    return index;
  }

  eachPrimes(
    function(prime) {
      var index = 0,
          doseContinue = true;

      index = indexOfPrime(prime);
      if (index > 0)
        doseContinue = callback(prime, index);
      if (doseContinue === false || n === 1)
        return false;
    });
}

/**
 * Call a function on each primes.
 *
 * @param {function(number):boolean} callback
 *   Break iterator when return false.
 */
function eachPrimes(callback) {
  var primes = new Primes,
      doseContinue = true;

  while (doseContinue !== false)
    doseContinue = callback(primes.next());
}

// {{{ Primes
/**
 * Primes iterator.
 *
 * @constructor
 */
function Primes() {
  /** @type {number} */
  this.currentIndex = -1;
}

/**
 * Known primes.
 *
 * @type {Array.<number>}
 */
Primes.primes = [2, 3];

Primes.prototype = {
  /**
   * Next prime.
   *
   * @return {number}
   */
  next: function() {
    var n, i, iz,
        _primes = Primes.primes,
        is_prime = false;

    if (_primes[this.currentIndex + 1]) {
      this.currentIndex += 1;
      return _primes[this.currentIndex];
    }
    n = _primes.last;
    while (! is_prime) {
      n += 2;
      for (i = 0, iz = _primes.length; i < iz; ++i) {
        if (n % _primes[i] === 0)
          break;
        if (_primes[i] * _primes[i] > n)
          is_prime = true;
      }
    }
    _primes.push(n);
    this.currentIndex += 1;
    return n;
  },

  /**
   * Reset the iterator pointer to the first.
   */
  first: function() { this.currentIndex = -1; }
};
// }}} Primes
}());

assert = require('assert');
assert.equal(mobius(3 * 19 * 19 * 23), 0);
assert.equal(mobius(3 * 19 * 23), -1);
assert.equal(mobius(3 * 19 * 23 * 151), 1);
// vim:set et sw=2 sts=2 ff=unix foldmethod=marker:

RubyでMöbius函數

JavaScript版とほぼ同じ実装にしたが、素數列をEnumerableにした。

# encoding=utf-8
# license: Public Domain

# Primes iterator.
#
# == Usage
# Primes.new.each{|prime| p prime }
# Primes.new.take(100).each{|prime| p prime }
# Primes.new.take_while{|prime| prime < 1000 }.each{|prime| p prime }
# Primes.new.lazy.select{|prime| prime > 1000 }.take(10).each{|prime| p prime }
class Primes
  include Enumerable
  @@primes = [2, 3]

  # Iterate each primes.
  def each
    current_index = 0
    loop do
      @@primes << next_prime unless @@primes[current_index]
      yield @@primes[current_index]
      current_index += 1
    end
  end

private
  # Find a next prime.
  #
  # return [Number]
  def next_prime
    n = @@primes.last + 2
    n += 2 while @@primes.
      take_while{|prime| prime ** 2 <= n }.
      any?{|prime| n % prime == 0 }
    n
  end
end

# Prime factors iterator.
class PrimeFactors
  include Enumerable

  # param [Number]
  def initialize n
    @number = n
  end

  # Iterate each prime factors.
  def each
    n = @number
    Primes.new.each do |prime|
      index, n = prime_index prime, n
      yield prime, index if index > 0
      break if n == 1
    end
  end

  # Divide number with the prime.
  #
  # == Usage
  # prime_index 3, 162 == [4, 2]
  # 3 ** 4 * 2 == 162
  #
  # param [Number] prime
  # param [Number] n
  # return [Array] [index, n]
  def prime_index prime, n=@number
    index = 0
    while n % prime == 0
      n /= prime
      index += 1
    end
    [index, n]
  end
end

class Integer
  # Generate an iterator on each prime factors.
  #
  # == Usage
  # 360.prime_factors.each{|prime, index| p [prime, index] }
  #
  # return [PrimeFactors]
  def prime_factors
    PrimeFactors.new self
  end
end

# Mobius function.
#
# param [Number] n
# return [Number]
def mobius n
  if n.prime_factors.any?{|prime, index| index > 1 }
    0
  elsif n.prime_factors.to_a.length % 2 == 1
    -1
  else
    1
  end
end

p mobius(3 * 19 * 19 * 23) == 0
p mobius(3 * 19 * 23) == -1
p mobius(3 * 19 * 23 * 151) == 1
# vim:set et sw=2 sts=2 ff=unix foldmethod=marker: