juliaで可変ワード長の演算を試す

juliaが結構便利なので、勉強がてら有限ワード長のモジュールを作っています。 その過程をメモ代わりにそのまま吐き出してみようと思います。

  • 概要を考える

ワード長は可変にしたいので、値とワード長の二つのフィールドが必要になりそうです。値の方は、簡単のため UInt で済ませます。 演算をするのであれば、コンストラクタと加減算、あとは、連結はほしいところです。

というざっくりしたイメージでスタートしました。

  • コンストラクタ

こちらは、特に悩む必要もなくtypeを使うだけで終わります。 任意ワード長の RegVarのほかに、8bit固定の Reg8 も用意しておこうかと思い、これらをまとめる抽象型 AbstractReg も作りました。

abstract AbstractReg
type RegVar <: AbstractReg
  val :: UInt
  wl :: Int
end
type Reg8 <: AbstractReg
  val :: Uint
end
  • 加算

Baseの+関数をimportして、ワード長を考慮しながら計算します。

import Base: +
export +
+(a::RegVar, b::RegVar) = addreg(a, b)
function addreg%(a::RegVar, b::RegVar)
  n = max(a.wl, b.wl)
  ub = 2^(n-1)-1
  lb = -2^(n-1)
  y = a.val + b.val
  if (y < lb)
    y = y + 2^n
  else (y > ub)
    y = y - 2^n
  end
  RegVar(y, n)
end

RegVar としてはこれでいいのですが、Reg8に適用しようとしてもReg8はwlをメンバに持っていないのでうまくいきません。かといって、Reg8にwlを追加するのも嫌だな……。 と考えたところで、wl を総称関数にすればよいことに気づきました。

wl(::Reg8) = 8
wl(x::RegVar) = x.wl

これをつかって先ほどのaddregを書き換えると、こうなります。

import Base: +
export +
+(a::AbstractReg, b::AbstractReg) = addreg(a, b)
function addreg(a::AbstractReg, b::AbstractReg)
  n = max(wl(a), wl(b))
  ub = 2^(n-1)-1
  lb = -2^(n-1)
  y = a.val + b.val
  if (y < lb)
    y = y + 2^n
  else (y > ub)
    y = y - 2^n
  end
  RegVar(y, n)
end

a.wl などとしていたところを wl(a) にしてあげただけですが、これで Reg8 にも RegVar にも対応できるようになりました。

Juliaの場合、型のフィールドに対して自動的にアクセサが生成されるのですが、総称関数を定義するとこういうときに便利なのですね。納得。

  • 減算

減算は、単項演算子のマイナスを定義して、先ほどの加算と組み合わせることにします。

import Base: -
export -
-(x::AbstractReg) = RegVar((1 << wl(x)) - x.val, wl(x))
-(a::AbstractReg, b::AbstractReg) = a + (-b)

ワード長が64の時にどうするんだという話はありますが、とくに問題ないようです。

julia> convert(UInt, 1) << 64 - 1
0xffffffffffffffff
  • 連接

このあたりから自信がなくなってきます。 演算子は ++ にすることにして、二項演算子としてはすぐに書けました。ただ、これを連続して使おうとした場合にどうするかで手が止まります。 よくわからないので、julia本体から似たようなコードを探してみることに。sum なら可変長引数をとるかとおもって確認してみたところ、tuple.jlに記載がありそうです。

julia> methods(sum)
# 13 methods for generic function "sum":
sum(x::Tuple{Any,Vararg{Any,N<:Any}}) at tuple.jl:205
...

そこを確認したところ、ドット3つでタプルを展開してくれるような感じです。いわれてみれば、MATLABにそんなのがあったような気もします。

export ++
++(args...) = assoc(args)
assoc(x::Tuple{Any, Vararg{Any}}) = assoc(x...)
assoc(a::RegVar, args...) = assoc(a, assoc(args))
...

こんな感じで展開できました。

  • 補足:表示、等号

書き始めた当初に、完全に忘れていたのが表示と等号。これらがないと REPL やテストで困ることに、後で気づきました。 表示については、juliaの本体を検索すると、showという関数を定義すればよさそうなことがわかりました。

> ag "print\(" *.jl
...
version.jl
48:function print(io::IO, v::VersionNumber)
49:    v == typemax(VersionNumber) && return print(io, "∞")
50:    print(io, v.major)
51:    print(io, '.')
52:    print(io, v.minor)
53:    print(io, '.')
54:    print(io, v.patch)
...

というわけで、定義します。

import Base: show
function show(io::IO, x::AbstractReg)
  mroundup(x, n) = Int(ceil(x / n)) * n;
  print(io, wl(x));
  print(io, "'h");
  octmask = 1 << 4 -1;
  for i=mroundup(wl(x), 4)-4:-4:0
    z = octmask & (x.val >> i);
    if z < 10
      print(io, z);
    else
      print(io, 'a'+ (z-10));
    end
  end
end

これだけ定義すれば、REPL や print 関数などから自動的に呼んでもらえます。

julia> Reg8(10)
8'h0a

等号はこんな感じ。

==(a::AbstractReg, b::AbstractReg) = (a.val == b.val) && (wl(a) == wl(b))

できた。

配列の配列から行列への変換

配列の配列を、行列に変換する方法がわからなかったのでメモ。
結果としては以下のようにすれば良いようです。

import Base.convert
convert{T}(::Array{T,2}, x :: Array{Array{T,1},1}) = hcat(x...)'

今回学んだのは2点。

一つは、型自体で多重ディスパッチしたい時は ::型 でいけるということ。
もう一つは、... の挙動が理解できていないこと。
後者は、現段階では「配列やタプルを分解してくれるもの」というぼんやりした認識です。そういえば MATLAB にも同じ文法がありますね。
詳細について、いずれドキュメントを確認しようと思います。

Julia で scanl を使いたい

Julia で scanl に相当するものが見当たらなかったので、勉強がてら、考えてみました。

foldl が参考になるのではないかと思ってソースを見てみたところ、以下のようになっています。

foldl(op, v0, itr) = mapfoldl(identity, op, v0, itr)
mapfoldl(f, op, v0, itr) = mapfoldl_impl(f, op, v0, itr, start(itr))
function mapfoldl_impl(f, op, v0, itr, i)
  if done(itr, i)
    return r_promote(op, v0)
  else
    (x, i) = next(itr, i)
    v = op(r_promote(op, v0), f(x))
    while !done(itr, i)
      @inbounds (x, i) = next(itr, i)
      v = op(v, f(x))
    end
  end
end

色々呼び出していますが、f は map 用の関数で、op は v に結果を次々と付け加えていく関数と考えれば良さそう。
Haskell と(雰囲気だけですが)対応させると、こんな感じでしょうか。

itr :: (Ix i) => (b, i)
f :: b -> a
v :: a
op :: a -> a -> a  -- mappend 相当
mapfoldl f op v0 itr = foldl op v0 $ map f itr

scanl を作るためには、2つの引数を取る関数とmappend が必要ですが、 mapreduce には前者がないので、これで scanl を実現するのはちょっと難しそうです。

以上を踏まえ、素直に書くことにしました。こんな感じです。

 function scanl{T}(op, v0 :: T, xs)
   function scanList(op, vs, xs)
     isempty(xs) ? vs : scanList(op, [vs; op(vs[end], xs[1])], xs[2:end])
   end
   scanList(op, [v0], xs)
 end

再帰を使っていますが、foldl の実装を見ると素直にループを回した方が良いのかもしれません。この辺りも、Julia なら簡単に検討できるはずなので、いずれやってみようと思います。

Julia で関数合成マクロをつくる

Julia は便利そうなのですが、Haskell から移ってくるとカッコの数が気になってきます。
せっかく Julia ではマクロが作れるので、関数合成を簡単に書けないかと考えてみました。一気に書くことができなかったので、関数列と引数を与えると合成した結果を返すマクロ @rc を用意し、さらに @compose で無名関数にするという2段構えです。
(もっとスマートな方法があれば教えてください)

macro rc(args ...)
  if length(args) == 1
    :($(args[1]))
  else
    :($(args[1])(@rc($(args[2:end]...))))
  end
end

macro compose(args ...)
  :(x -> @rc($(args ...), x))
end

rc は、recursive compose のつもり。使い方はこんな感じになります。

>>> @rc(asin, sin, 1)  # => 1
>>> @compose(asin, sin)  # =>(generic function with 1 method)
>>> @compose(asin, sin)(1)  # => 1

MATLABでカリー化

一度関数型プログラムの考え方に触れると、どの言語でも同じようなことがしたくなってしまいます。
MATLAB は便利ですが、当然関数型ではないので勉強がてらちょっとずつ自作しています。
その過程で、カリー化をしようとしてハマったのでメモ。

結論からいうと、セル配列が使えました。セル配列は型の違うものをひとまとめに扱うしくみで、
取り出すときには小括弧と中括弧の2つの方法があります。
中括弧を使って範囲指定をして取り出すと、多値で答えが帰ってくるので、これを使って curry がつくれました。

function f = curry(f, varargin)
%curry do currying
%   feval(curry(@plus, 1), 10) => 11
%   feval(curry(@strcat, 'a', 'b', 'c'), 'z') => 'abcz'

if nargin >= 2
    xs = varargin;
    f = curry(@(varargin) f(xs{1}, varargin{:}), xs{2:end});
end

第1引数の f をそのまま返値にしているので、 nargin == 1 の場合はそのまま f が返ります。
この辺り、お行儀が良いのかはよくわかりません。明示的に分岐を書いた方が読みやすい気もします。

# 追記
と思ったら、まんま stack overflow に書いてあった。。。

Stackの味見

他人事だと思っていた cabal hell についにはまったので、cabal から stack に乗り換えて見ることにしました。

brew uninstall ghc
brew install ghc
brew haskell-stack
stack new test-proj
cd ./test-proj
vim ./test-proj.cabal
# library の build-depends にカンマ区切りで使用するパッケージを追加
stack build
stack ghci

ソースファイルに import を書くのに加え、.cabal ファイルにも記載が必要なので二度手間のような感じです。
多分、ソースから .cabal に反映させる方法はあるのでしょうが、まだ見つけられていません。

haskell での単位付き計算を味見

昔から計算ミスの多かった私は、高校時代に塾の先生に「次元の確認をするとよい」とアドバイスされて以来、必ず次元チェックをするようにしています。
市販ツールでは Mathcad が便利だと思うのですが、個人で購入するには少々高価です。どうにかならないかと常々思っていたのですが、haskell に dimensional というモジュールがあったので、味見してみました。

School of Haskell で使い方を確認しつつ、試してみます。ghci上でこんな感じ。

import Numeric.Units.Dimensional.Prelude
import qualified Prelude

l = 1.0 *~ kilo meter
speeds = [10,20,30] *~~ (kilo meter / hour)
putStrLn . show $ l / speeds

できた。でも、記法が長ったらしいので、もうちょっとなんとかならないかと以下のようにしてみました。

(#) val unit = val *~ unit
mV = milli volt
uA = micro ampere
kOhm = kilo ohm

演算子はシャープにしていますが、頻出するのでもうすこし存在感が薄い文字が使いたいところです。
で、たとえば

1#mV / 0.5#uA / 1#kOhm

とすれば概ね 2 がかえってきます。