30のプログラミング言語でFast inverse square rootを実装してみました！

あなたの好きな言語は何ですか。そして、あなたの好きなアルゴリズムは何ですか。

好きな言語があると、その言語でどんな問題でも解決しようとなりがちになります。その言語を極めるのは素晴らしいことですが、その言語や似たような言語でしかコードが書けなくなったり、他の言語に対して見向きもしなくなってしまう可能性があります。

勇気を出して新しい言語にチャレンジしてみませんか？色々な言語に挑戦してみませんか？

何から始めればいいか分からない。次にどの言語を学べばいいか分からない。いま特に何も困っていない。何でも得意な言語で書けてしまう。そういう人が多いのではないでしょうか。

新しい言語にチャレンジするきっかけを作る一つの方法は、ある特定のアルゴリズムを一つ決めて、あらゆる言語で実装してみることです。解く問題が大きすぎると力尽きてしまうので、せいぜい20〜30行程度で書ける簡単なものが良いでしょう。大事なことは、そのコードが実際に実行できることを手元で確認することです。実行処理系を実際にインストールし、コードをコンパイルして実行してみる、それだけでも多くのことを学ぶことができます。

実際、私もこういうチャレンジをするのは今回が初めてです。このエントリーでは、Fast inverse squre rootを30の言語で書いてみたので、それぞれのコードを紹介したいと思います。

なお、以下のコードは全てGitHubに公開しています。

Fast inverse square rootアルゴリズム

まずC言語による実装を見てみましょう。

#include <stdio.h>
#include <stdlib.h>

float fastInvSqrt(float x) {
  int i = *(int*)&x;
  i = 0x5f3759df - (i >> 1);
  float y = *(float*)&i;
  return y * (1.5F - 0.5F * x * y * y);
}

int main(void) {
  char *line = NULL, *endptr = NULL;
  size_t size;
  float x;
  while (getline(&line, &size, stdin) != -1) {
    x = strtof(line, &endptr);
    if (*endptr == '\n') {
      printf("%f\n", fastInvSqrt(x));
    }
  }
  return 0;
}

32bit浮動小数点数の平方根の逆数の近似を求めるアルゴリズムです。

あるベクトルからその方向の単位ベクトルを求める (各要素の二乗和の平方根の逆数が必要になる) ために、かつて実際に製品で使われたことがあるそうです。 あくまで近似値ではありますが、およそ0.3%未満の誤差で求まります。

アルゴリズムのミソは、float⇔intのpointer castingにあります。

  int i = *(int*)&x; // float x

  float y = *(float*)&i; // int i

この変換によって x と i (y と i) はどういう関係にあるのでしょうか。

32bitの浮動小数点数は、下のbitから23bitの仮数部と、127のゲタを履かせた8bitの指数部と1bitの符号部から成ります。 従って、正の32bit浮動小数点数 と、そのビット列を32bitの整数値として解釈した値 は、次の近似式が成り立ちます。

特別なケースとしてが2の冪である場合 (仮数部 = 0) は、近似なしに次のようになります。

例えば、floatの16.0ã‚’intにpointer castingすると 、0.125ã‚’intにすると のようになります。

#include <assert.h>
#include <stdio.h>

int main(void) {
  int i, j = 30, passed = 0;
  float x = (float)(1 << j);
  while (j >= -30) {
    i = *(int*)&x;
    assert(i == ((127 + j) << 23));
    x /= 2; j--; passed++;
  }
  printf("%d test passed\n", passed);
  return 0;
}

求めたい値を とすると、 となりますので、次の関係式が得られます。

これが i = 0x5f3759df - (i >> 1); の行に対応します。 この値をまたpointer castingでfloatに戻すことで、最初の近似値が得られます。 さらに、ニュートン法のiterationを一回行うことで精度を高めています。

更に詳しい説明は、いくらでも文献はでてきますので、そちらをご参照ください。 https://www.google.com/search?q=Fast+inverse+square+root

さて、あらゆる言語で一つのアルゴリズムを書こうとなった時に、私はなぜこのアルゴリズムを選んだのでしょうか。以下の理由があげられます。

• コードが小さすぎず大きすぎない。ただHello worldするだけならチュートリアルの最初を見ただけで終わってしまい、面白くありません。一方で、問題が難しいと様々な言語で書いている途中で挫折して投げ出してしまいそうです。
• float⇔intのpointer castingがアルゴリズムのポイントです。これをどう実装するか、初めて触る言語においていかに速く見つけるかというところが重要になってきます。pointerがない言語ではbinary表現を介して変換するわけですが、これもドキュメントをわりと探さないと見つかりません。この変換がちゃんとできないと、めちゃくちゃ答を吐くことになるでしょう。
• なによりも、このアルゴリズムが好きだから。浮動小数点がどう表わされるかを実にうまく使っています。最初に気がついた人は頭がいいなぁと思います。そして、今やもう何の役にも立たないという儚さにも心惹かれるところです。
• ほとんどの言語において、このアルゴリズムを書いた人間はおそらく自分が初めてで (マイナーなアルゴリズムだしスクリプト言語で書いても無意味)、前人未踏の地を制覇していく感覚が楽しい。無意味なことはなぜか楽しい。

改めて、問題のレギュレーションをプログラミングコンテスト風に並べてみます。

• 正の浮動小数点の数が改行区切りで与えられるので、その平方根の逆数を改行区切りで出力せよ。
• float⇔int32のFast inverse square rootアルゴリズムを用いること。
• seq 10000000 を入力としても動作すること。
• 数値として不正な入力は無視して良いが、例外などで落ちないようになっていると好ましい。
• 許される誤差は、真値に対して0.5%未満とする。

サンプル入力

100
10
1
0.1
0.01
0.0025

サンプル出力

0.099845
0.315686
0.998307
3.157232
9.982522
19.965044

さあ、問題設定ができました。 他の言語でこの問題を解くとどうなるのか、違いを楽しみながらご覧ください。

C++

#include <iostream>
#include <string>

float fastInvSqrt(float x) {
  int i = *reinterpret_cast<int*>(&x);
  i = 0x5f3759df - (i >> 1);
  float y = *reinterpret_cast<float*>(&i);
  return y * (1.5F - 0.5F * x * y * y);
}

int main() {
  std::string line;
  while (std::getline(std::cin, line)) {
    try {
      std::cout << fastInvSqrt(std::stof(line, NULL)) << std::endl;
    } catch (...) {}
  }
  return 0;
}

C++には4つのcast演算子がありますが、今回は最も危険なcastであるreinterpret_castを使わざるを得ない場面です。 また、文字列からfloatへの変換にはstr::stofを使ってみました。 この関数は、入力が正しくないと例外を吐きます。

C#

using System;

class FastInvSqrt
{

    public static void Main()
    {
        string line;
        float number;
        while ((line = Console.ReadLine()) != null)
        {
            if (float.TryParse(line, out number))
            {
                Console.WriteLine(fastInvSqrt(number));
            }
        }
    }

    unsafe public static float fastInvSqrt(float x)
    {
        int i = *(int*)&x;
        i = 0x5f3759df - (i >> 1);
        float y = *(float*)&i;
        return y * (1.5F - 0.5F * x * y * y);
    }

}

Microsoftによって開発された.NET系言語です。 Cの関数にunsafeキーワードをつけるとそのまま動きました。 float.TryParseが変換できたかを例外ではなく返り値で扱えるのは素敵です。 参照渡しを表すrefキーワードがあることも、言語として面白いところです。

Visual Basic

Imports System
Imports System.IO

Public Class FastInvSqrt

    Shared Sub Main(args As String())
        Dim s As String = Console.ReadLine()
        Dim x As Single
        While s IsNot Nothing
            If Single.TryParse(s, x)
                Console.WriteLine(fastInvSqrt(x))
            End If
            s = Console.ReadLine()
        End While
    End Sub

    Shared Function fastInvSqrt(x As Single) As Single
        Dim i As Integer = BitConverter.ToInt32(BitConverter.GetBytes(x), 0)
        i = &H5f3759df - (i >> 1)
        Dim y As Single = BitConverter.ToSingle(BitConverter.GetBytes(i), 0)
        Return y * (1.5 - 0.5 * x * y * y)
    End Function

End Class

.NET系言語です。 BitConverterを使って数値の変換を行ってみました。 Main関数はC#とほぼ同じです。 DimキーワードやEnd Subなどが懐かしい感じがします。 16進数表現リテラルの &Hという表記はユニークです。

F#

let fastinvsqrt(x: float32): float32 =
    let i = System.BitConverter.ToInt32(System.BitConverter.GetBytes(x), 0)
    let j = 0x5f3759df - (i >>> 1)
    let y = System.BitConverter.ToSingle(System.BitConverter.GetBytes(j), 0)
    y * (1.5F - 0.5F * x * y * y)

seq { while true do yield stdin.ReadLine() }
    |> Seq.takeWhile ((<>) null)
    |> Seq.iter (fun x -> try float32 x |> fastinvsqrt |> printfn "%f"
                          with :? System.FormatException -> ())

.NETの言語が続きます。 BitConverterを使っていることなど、fastinvsqrtの実装はVisual Basicとほぼ同じです。 floatが64bitの浮動小数点 (=double) を表すというのが少し違和感があります。 1.5と書くとfloat (=double) になりコンパイルエラーになるので1.5F (float32 = single) と書く必要があります。 例外のtry ... with :?や、pipe operator (|>) が特徴的です。

Perl

use strict;
use warnings;
use feature qw(say);

while (<>) {
  say fastInvSqrt($_);
}

sub fastInvSqrt {
  my ($x) = @_;
  my $i = unpack("l", pack("f", $x));
  $i = 0x5f3759df - ($i >> 1);
  my $y = unpack("f", pack("l", $i));
  $y * (1.5 - 0.5 * $x * $y * $y);
}

歴史のある言語ですが、テキスト処理が得意で、今なおwebアプリケーションに広く使われています。 pointerを扱えない言語では、pack/unpackによりbinary表現を介して実装することができます。 厳密にはpointer castingしていないコードはFast inverse square rootと言うべきではないかもしれませんが、それでは書ける言語がかなり限られてしまいますのでお許しください。

PHP

<?php
while ($line = fgets(STDIN)) {
  echo fastInvSqrt(floatval($line)), "\n";
}

function fastInvSqrt($x) {
  $i = unpack('l', pack('f', $x))[1];
  $i = 0x5f3759df - ($i >> 1);
  $y = unpack('f', pack('l', $i))[1];
  return $y * (1.5 - 0.5 * $x * $y * $y);
}
?>

webアプリケーションに特化した言語。 Perlのコードを直しただけです。 フォーマット指定も同じです。

Ruby

def fastInvSqrt(x)
  i = [x].pack("e").unpack("i")[0]
  i = 0x5f3759df - (i >> 1)
  y = [i].pack("i").unpack("e")[0]
  y * (1.5 - 0.5 * x * y * y)
end

if __FILE__ == $0
  STDIN.each_line do |line|
    p fastInvSqrt(line.to_f)
  end
end

オブジェクト指向の言語で、webアプリケーションフレームワークRuby on Railsを中心として広く使われています。 Perl譲りのpack/unpackが使えるので、少し書き直すだけで動きました。 一行ずつ処理できるSTDIN.each_lineも素敵です。

Python

import struct
import sys

def fastInvSqrt(x):
    i = struct.unpack('>i', struct.pack('>f', x))[0]
    i = 0x5f3759df - (i >> 1);
    y = struct.unpack('>f', struct.pack('>i', i))[0]
    return y * (1.5 - 0.5 * x * y * y)

if __name__ == '__main__':
    for line in iter(sys.stdin.readline, ""):
        print(fastInvSqrt(float(line)))

webアプリケーションから数式処理や機械学習まで幅広く使用されている言語。 やはりpack/unpackがあるので、binary表現を介して変換しました。 フォーマット文字に着目すると、一文字目の>ã‚„< (他に @, =, !) などによりendianを指定する方法は、上の3つの言語は異なっています。

R

fastInvSqrt <- function(x) {
  i <- readBin(writeBin(as.numeric(x), raw(), size=4), integer(), size=4)
  i <- 0x5f3759df - bitwShiftR(i, 1)
  y <- readBin(writeBin(as.integer(i), raw(), size=4), double(), size=4)
  y * (1.5 - 0.5 * x * y * y)
}

f <- file("stdin")
open(f)
while (length(line <- readLines(f, n=1)) > 0) {
  write(fastInvSqrt(as.numeric(line)), stdout())
}

統計処理に広く用いられている言語です。 ポインタもありませんし、浮動小数点型に32bitと64bitの区別がないので苦労します。 readBinにdouble, size=4を指定すると32bitの数値として読み込んでくれます。 他のシステムが吐いたデータファイルを読み込むという需要はあるためでしょうか。 主にファイルとのやり取りに使う関数ですが、raw vectorを引数に取ることもできます。

JavaScript

require('readline').createInterface({
    input: process.stdin,
    output: process.null
}).on('line', function(line) {
    console.log(fastInvSqrt(parseFloat(line)));
});

function fastInvSqrt(x) {
    var i = floatToUInt32(x);
    i = 0x5f3759df - (i >> 1);
    var y = uint32ToFloat(i);
    return y * (1.5 - 0.5 * x * y * y);
}

function floatToUInt32(x) {
    var buf = new ArrayBuffer(4);
    new Float32Array(buf)[0] = x;
    return new Uint32Array(buf)[0];
}

function uint32ToFloat(i) {
    var buf = new ArrayBuffer(4);
    new Uint32Array(buf)[0] = i;
    return new Float32Array(buf)[0];
}

ブラウザーで動くことはもちろん、JavaScriptã‚„CSSへのトランスパイラがJavaScriptで書かれていることもあり、web開発には欠かせない言語。 pointerなどもちろんなく、また数値にfloatとdoubleの区別はおろか、整数型と浮動小数点型の区別もなく、Rよりもさらに不利な状況です。 TypedArrayをうまく使うとなんとか実装することができました。 ArrayBufferを引数にしてUint32Arrayã‚„Float32Arrayを作れるのがポイントです。

CoffeeScript

require 'readline'
  .createInterface
    input: process.stdin
    output: process.null
  .on 'line', (line) =>
    console.log fastInvSqrt parseFloat line

fastInvSqrt = (x) ->
  i = floatToUInt32 x
  i = 0x5f3759df - (i >> 1)
  y = uint32ToFloat i
  y * (1.5 - 0.5 * x * y * y)

floatToUInt32 = (x) ->
  buf = new ArrayBuffer 4
  (new Float32Array(buf))[0] = x
  new Uint32Array(buf)[0]

uint32ToFloat = (i) ->
  buf = new ArrayBuffer 4
  (new Uint32Array(buf))[0] = i
  new Float32Array(buf)[0]

JavaScriptに変換されるAltJSの一つで、JavaScript自体の構文にも大きな影響を与えました。 上のJavaScriptのコードが吐かれるように書くだけです。 (new Float32Array(buf))[0] = x ã‚’ new Float32Array(buf)[0] = x と書けないのは残念なところです。

LiveScript

require \readline
  .createInterface do
    input: process.stdin
    output: process.null
  .on \line, (line) ->
    line |> parseFloat |> fastInvSqrt |> console.log

fastInvSqrt = (x) ->
  i = floatToUInt32 x
  i = 0x5f3759df - (i .>>. 1)
  y = uint32ToFloat i
  y * (1.5 - 0.5 * x * y * y)

floatToUInt32 = (x) ->
  buf = new ArrayBuffer 4
  (new Float32Array buf).0 = x
  (new Uint32Array buf).0

uint32ToFloat = (i) ->
  buf = new ArrayBuffer 4
  (new Uint32Array buf).0 = i
  (new Float32Array buf).0

CoffeeScriptとよく似ているAltJSですが、さらに先進的な機能が入っています。 バックスラッシュから始まる単語が文字列になったりする (\line == 'line') のは、他の言語では見たことがないsyntaxです (しいて言えばRubyのsymbolに似ています)。 pipe operator |> ã‚„ composition operator >> など、F#に大きな影響を受けていることが分かります。 >> が関数合成なので、シフト演算子が.>>.になっています。

PureScript

module Main where

import Prelude
import Control.Monad.Eff (Eff)
import Control.Monad.Eff.Console (CONSOLE, logShow)
import Control.Monad.Eff.Exception (EXCEPTION)
import Data.ArrayBuffer.ArrayBuffer (create)
import Data.ArrayBuffer.DataView (READER, WRITER, whole, getInt32, getFloat32, setInt32, setFloat32)
import Data.Int.Bits (shr)
import Data.Maybe (fromJust)
import Data.Monoid (mempty)
import Global (readFloat)
import Node.Process (stdin)
import Node.ReadLine (READLINE, setLineHandler, createInterface)
import Partial.Unsafe (unsafePartial)

main :: forall e. Eff (console :: CONSOLE, readline :: READLINE, err :: EXCEPTION, reader :: READER, writer :: WRITER | e) Unit
main = do
  interface <- createInterface stdin mempty
  setLineHandler interface $ \line ->
    logShow =<< fastInvSqrt (readFloat line)

fastInvSqrt :: forall e. Number -> Eff (reader :: READER, writer :: WRITER | e) Number
fastInvSqrt x = do
  i <- floatToUInt32 x
  let j = 0x5f3759df - unsafePartial i `shr` 1
  y <- uint32ToFloat j
  pure $ y * (1.5 - 0.5 * x * y * y)

floatToUInt32 :: forall e. Number -> Eff (reader :: READER, writer :: WRITER | e) Int
floatToUInt32 x = do
  let dataView = whole $ create 4
  setFloat32 dataView x 0
  unsafePartial $ map fromJust $ getInt32 dataView 0

uint32ToFloat :: forall e. Int -> Eff (reader :: READER, writer :: WRITER | e) Number
uint32ToFloat i = do
  let dataView = whole $ create 4
  setInt32 dataView i 0
  unsafePartial $ map fromJust $ getFloat32 dataView 0

Haskellに似たsyntaxを持つAltJSです。 JavaScriptにコンパイルされる言語ですので、基本的にはJavaScriptのコードを直すだけなのですが、Haskellの経験がないと難しいでしょう。 Eff monadの型定義やPartial type classなど、Haskellと異なるところも多くあります。

J

fastInvSqrt =: 3 : 0
  i =. _2 ic 1 fc y
  i =. (dfh'5f3759df') - _1 (33 b.) i
  z =. _1 fc 2 ic i
  z * (1.5 - 0.5 * y * z * z)
)

echo & fastInvSqrt & (0 & ".) & (-. & CRLF) ;. 2 &. stdin ''
exit ''

書くのも読むのも難しい言語ですが、数日ほど家に籠もってドキュメントを読みこめばそれなりに分かるようになります (Vocabularyがとても便利)。 関数ではなく動詞と表現したり、動詞の挙動を変える関数を副詞と表現したりする言い回しはかなりユニーク。 演算子が常に右結合だったり、ic =: 3!:4ã‚„fc =: 3!:5などによる数字とbinary列の変換、dfhによる16進数表現の変換、0 ". yによる文字列から数字への変換、3 : 0による動詞の定義など、上のたった数行のコードが書けるようになるまでに学ばなければならないことは多い。

Go

package main

import (
    "bufio"
    "fmt"
    "os"
    "strconv"
    "unsafe"
)

func fastInvSqrt(x float32) float32 {
    i := *(*int32)(unsafe.Pointer(&x))
    i = 0x5f3759df - i>>1
    y := *(*float32)(unsafe.Pointer(&i))
    return y * (1.5 - 0.5*x*y*y)
}

func main() {
    scanner := bufio.NewScanner(os.Stdin)
    for scanner.Scan() {
        if x, err := strconv.ParseFloat(scanner.Text(), 32); err == nil {
            fmt.Println(fastInvSqrt(float32(x)))
        }
    }
}

Googleが開発したコンパイル言語です。 ここまでしばらくスクリプト言語が続いていましたが、ここからコンパイル言語が続きます。 普通にpointer castingできるだけで安心感が出てきますね。 unsafe.Pointerという字面は、いかにも危険なことをやっているんだということを喚起してくれます。 strconv.ParseFloatの返り値の扱いにも見られるように、エラーは例外ではなく複数返り値で扱うのが慣習となっています。 標準のコードフォーマッターgofmtでコードを整形する慣習も、いい文化です。

Swift

func fastInvSqrt(x: Float) -> Float {
    var z: Float = x
    var i: Int = UnsafeMutablePointer<Int>(withUnsafeMutablePointer(&z, { $0 })).memory
    i = 0x5f3759df - (i >> 1)
    let y: Float = UnsafeMutablePointer<Float>(withUnsafeMutablePointer(&i, { $0 })).memory
    return y * (1.5 - 0.5 * x * y * y)
}

while let line = readLine() {
    if let x = Float(line) {
        print(fastInvSqrt(x))
    }
}

Appleが開発し、WWDC 2014で発表された言語です。 普通はこんな危なっかしいコード書かないでしょうけど、低レイヤーとのインターフェースも用意されており、ポインタにアクセスすることができます (たすかりました)。 if letã‚„while letなど、optional bindingは素敵だと思います。 変数 var と定数 let は区別されます。

Rust

use std::io::BufRead;

fn fast_inv_sqrt(x: f32) -> f32 {
    let i = unsafe { *std::mem::transmute::<&f32, &i32>(&x) };
    let j = 0x5f3759df - (i >> 1);
    let y = unsafe { *std::mem::transmute::<&i32, &f32>(&j) };
    y * (1.5 - 0.5 * x * y * y)
}

fn main() {
    let stdin = std::io::stdin();
    for line in stdin.lock().lines().filter_map(|x| x.ok()) {
        match line.parse::<f32>() {
            Ok(x) => println!("{}", fast_inv_sqrt(x)),
            Err(_) => {},
        }
    }
}

Mozillaにより開発されている言語です。 std::mem::transmuteという関数により、pointer castingを行うことができます。 unsafeキーワードはC#に少し似ています。 危険なことをしている感じがよくでています。

D

import std.conv;
import std.stdio;
import std.string;

void main()
{
    foreach (line; stdin.byLine())
    {
        try
        {
            writeln(fastInvSqrt(to!float(chomp(line))));
        }
        catch (ConvException e) {}
    }
}

float fastInvSqrt(float x)
{
    int i = *cast(int*)(&x);
    i = 0x5f3759df - (i >> 1);
    float y = *cast(float*)(&i);
    return y * (1.5F - 0.5F * x * y * y);
}

赤いDの文字の形をしたキャラクターが特徴的なコンパイル言語です。 綺麗なC/C++という印象です。 上のコードでは取り立てて特徴なところはありません。

Crystal

def fastinvsqrt(x : Float32) : Float32
  i = pointerof(x).as(Int32*).value
  i = 0x5f3759df - (i >> 1)
  y = pointerof(i).as(Float32*).value
  y * (1.5 - 0.5 * x * y * y)
end

while line = gets
  puts fastinvsqrt(line.to_f32)
end

Rubyに似た構文を持つ、静的型付けのコンパイル言語です。 C言語との運用も考慮されていて、pointerを扱うインターフェースが用意されています。 このsyntaxは他の言語ではあまり見かけませんが、私はかなり好きですね。 C言語風の&やら*よりも読みやすい気がするのですが、どうでしょうか。 コード全体も引き締まっていてシンプルです。 コンパイル言語で、実行時のパフォーマンスも優秀です。

Nim

import strutils

proc fastInvSqrt(x: float32): float32 =
  let i = cast[int32](x)
  let j = 0x5f3759df - i shr 1
  let y = cast[float32](j)
  return y * (1.5 - 0.5 * x * y * y)

while true:
  try:
    let input = stdin.readLine
    try:
      echo fastInvSqrt(parseFloat(input))
    except ValueError:
      continue
  except IOError:
    break

こちらはPythonに似た構文を持つ、静的型付けの言語です。 この言語もCrystal同様、わりと低レイヤーが触れる言語です。 bit patternを保存するtype castはcast[T](v)のようにとてもすっきりしています。 あえてpointerでのcastであることを強調するために、let i = cast[ptr int32](addr(z))[] と書くこともできます。 addrが取れるのは変更可能な変数だけというのはSwiftに似ています。 pointer dereferenceがptr[]なのは他では見たことがありません (Cでは *ptr) が、配列と似ていてこちらのほうが一貫性があるように思えてきます。 Crystalと同様、実行時パフォーマンスは優秀です。

Java

import java.util.Scanner;

public class FastInvSqrt {

    public static void main(String[] args) {
        Scanner scanner = new Scanner(System.in);
        while (scanner.hasNext()) {
            try {
                float x = Float.parseFloat(scanner.nextLine());
                System.out.printf("%f\n", fastInvSqrt(x));
            } catch (NumberFormatException e) {}
        }
    }

    public static float fastInvSqrt(float x) {
        int i = Float.floatToRawIntBits(x);
        float y = Float.intBitsToFloat(0x5f3759df - (i >> 1));
        return y * (1.5F - 0.5F * x * y * y);
    }

}

業務系、webアプリケーション、スマホアプリ開発など、広く使われているオブジェクト指向の言語です。 生のpointerはありませんが、標準ライブラリの関数でfloatとintを変換する関数があるのでこれを用いることができます。 文字列からfloatへの変換関数Float.parseFloatが例外NumberFormatExceptionを投げるので、これをcatchする必要があります。

Groovy

Scanner scanner = new Scanner(System.in)
while (scanner.hasNext()) {
  try {
    println fastInvSqrt(Float.parseFloat(scanner.nextLine()))
  } catch (NumberFormatException e) {}
}

def fastInvSqrt(float x) {
  int i = Float.floatToRawIntBits(x)
  float y = Float.intBitsToFloat(0x5f3759df - (i >> 1))
  y * (1.5F - 0.5F * x * y * y)
}

JVM上の言語でJavaと親和性が高い言語。 上記のコードもJavaのコードをそっくりそのまま持ってきたものです。 セミコロンやreturnを書かなくてよくてすっきりしています。

Kotlin

fun main(args: Array<String>) {
  while (true) {
    val line = readLine()
    if (line == null) break
    try {
      println(fastInvSqrt(line.toFloat()))
    } catch (e: NumberFormatException) {}
  }
}

fun fastInvSqrt(x: Float): Float {
  val i = java.lang.Float.floatToRawIntBits(x)
  val y = java.lang.Float.intBitsToFloat(0x5f3759df - (i shr 1))
  return y * (1.5F - 0.5F * x * y * y)
}

これまたJVM上で実装されていて、Javaとの相互運用が可能な言語です。 右ビットシフトがshr中置演算子なのがユニークです。 nullの扱いがSwiftと似ています。

Scala

object FastInvSqrt {

  def main(args: Array[String]) {
    Iterator.continually(scala.io.StdIn.readLine()).takeWhile(_ != null).foreach { line =>
      util.control.Exception.catching(classOf[NumberFormatException]) opt line.toFloat map { x =>
        println(fastInvSqrt(x))
      }
    }
  }

  def fastInvSqrt(x: Float): Float = {
    val i = java.lang.Float.floatToRawIntBits(x)
    val y = java.lang.Float.intBitsToFloat(0x5f3759df - (i >> 1))
    y * (1.5F - 0.5F * x * y * y)
  }

}

同じくJVM上の言語です。 上記のJVM系の言語と同様ですが、String#toFloatが例外を投げるので、これをcatchする必要があります。 セミコロンもreturnも不要です。 valとvarがあります (他にKotlin, Swift, Rust, Nimなどがこの特徴を持ちます)。

Clojure

(defn fast-inv-sqrt [x]
  (let [i (Float/floatToRawIntBits x)
        y (Float/intBitsToFloat (- 0x5f3759df (bit-shift-right i 1)))]
    (* y (- 1.5 (* 0.5 x y y)))))

(doseq [line (line-seq (java.io.BufferedReader. *in*))]
  (try (println (fast-inv-sqrt (read-string line)))
      (catch java.lang.RuntimeException e ())))

JVM言語繋がりで、Clojureです。 Javaで書ければClojureで書くのも簡単です。 LISP系ですがletやdefnあたりの括弧の数が少なく、コードはすっきりしています。

OCaml

let fast_inv_sqrt x =
  let i = Int32.bits_of_float x in
  let j = Int32.sub (Int32.of_int 0x5f3759df) (Int32.shift_right i 1) in
  let y = Int32.float_of_bits j in
  y *. (1.5 -. 0.5 *. x *. y *. y)

let () =
  try
    while true do
      let line = input_line stdin in
      try
        let x = float_of_string line in
        Printf.printf "%f\n%!" (fast_inv_sqrt x)
      with
        Failure _ -> ()
    done;
  with
    End_of_file -> ()

floatは64bitの数字ですが、Int32.bits_of_floatが32bitでのbinary表現にしたものをint32で返してくれます。 float上の演算子は-. *.となっており、int上の演算子とは区別されています。 printfの%!でfflushするのは、他の言語ではあまり見ない特徴です。

Haskell

import Control.Monad ((<=<))
import Data.Bits (shiftR)
import Data.Traversable (mapM)
import Data.Word (Word32)
import Foreign.Marshal.Alloc (alloca)
import Foreign.Ptr (castPtr)
import Foreign.Storable (peek, poke)
import Prelude hiding (mapM)
import Text.Read (readMaybe)

main :: IO ()
main = mapM_ (mapM (print <=< fastInvSqrt) . readMaybe) . lines =<< getContents

fastInvSqrt :: Float -> IO Float
fastInvSqrt x =
  alloca $ \ptr -> do
    poke ptr x
    i <- peek (castPtr ptr)
    poke (castPtr ptr) $ 0x5f3759df - (i :: Word32) `shiftR` 1
    y <- peek ptr
    return $ y * (1.5 - 0.5 * x * y * y)

遅延評価が特徴的な純粋関数型言語です。 メモリー領域を用意して、そこに書き込んで(poke)、別の型で読み取る(peek)という発想が必要です。 Float -> Word32とその逆で二回の変換が必要ですが、サイズが同じなので同じメモリー領域を使い回すことができます。 getContentsしてlinesで分割してそれぞれの行を… と書いたコードが巨大な入力にもクラッシュしないのは遅延評価のおかげです。

Erlang

-module(fastinvsqrt).
-export([main/0, fast_inv_sqrt/1]).

main() ->
  case io:get_line("") of
    eof -> ok;
    {error, _} -> ok;
    Line -> catch io:format("~f~n", [fast_inv_sqrt(parse_float(Line))]), main()
  end.

fast_inv_sqrt(X) ->
  <<I:32/integer>> = <<X:32/float>>,
  J = 16#5f3759df - (I bsr 1),
  <<Y:32/float>> = <<J:32/integer>>,
  Y * (1.5 - 0.5 * X * Y * Y).

parse_float(Str) ->
  case string:to_float(Str) of
    {error, _} ->
      case string:to_integer(Str) of
        {error, _} -> error;
        {Int, _} -> float(Int)
      end;
    {Float, _} -> Float
  end.

並行処理を重視した動的型付けの関数型言語です。 このエントリーのコードの中で、動いた時の喜びが最も大きかったコードです。 << ... >>はバイナリー列を表します。 バイナリー列のパターンマッチによりfloatとintergerとの変換を行っています。 16進数表現リテラルの16#という表記、io:formatのフォーマット指定~f~n、ビットシフトのbsrなどはユニークです。 string:to_float("100")がエラーを返すのが残念です (いい方法があれば教えてください)。 カンマやセミコロンやドットなどが最初は難しい。

Elixir

defmodule Main do

  def main() do
    case IO.gets "" do
      :eof -> :ok
      {:error, _} -> :ok
      line ->
        case Float.parse(line) do
          :error -> :ok
          {float, _} -> IO.puts FastInvSqrt.fast_inv_sqrt(float)
        end
        main()
    end
  end

end

defmodule FastInvSqrt do

  use Bitwise, only_operators: true

  def fast_inv_sqrt(x) do
    <<i::size(32)-integer>> = <<x::float-32>>
    j = 0x5f3759df - (i >>> 1)
    <<y::float-32>> = <<j::size(32)-integer>>
    y * (1.5 - 0.5 * x * y * y)
  end

end

Erlangの仮想環境上で動作する言語、Elixir。 Erlangのコードとほぼ同じ。 ビットシフトは>>>なのは、>>がバイナリー列の閉じ括弧で既に使われていたからでしょうか。 syntaxはRubyに似ていますが、Erlangから学ばないと書きにくい。

Scheme

(use binary.pack)

(define (main args)
  (define (main-loop)
    (let ((line (read-line)))
      (cond ((not (eof-object? line))
             (let ((x (string->number line)))
               (cond ((not (eqv? x #f))
                      (display (fast-inv-sqrt (* 1.0 x)))
                      (newline))))
             (main-loop)))))
  (main-loop)
  0)

(define (fast-inv-sqrt x)
  (let* ((i (car (unpack "l" :from-string (pack "f" (list x) :to-string? #t))))
         (i (- #x5f3759df (ash i -1)))
         (y (car (unpack "f" :from-string (pack "l" (list i) :to-string? #t)))))
    (* y (- 1.5 (* 0.5 x y y)))))

LISP系の言語の一つです。 歴史が長く登場したのは1975年ですが、今でも教育の現場で使われています。 私も大学時代に触ったのを思い出しました。 binary.packライブラリーがあって助かりました。

Assembly

さて、せっかくここまで様々な言語で書いてきたわけですが、以上のコードは平方根の逆数を高速に求めるという目的において、もはや意味がありません。 X86の命令にrsqrtssという、まさに実装したいことそのままの命令があり、上記のあらゆるコードよりも速い (はずだ) からです。

fmt0: .string   "%f"
fmt1:    .string   "%f\n"

.globl    main
main:
    pushl    %ebp
    movl %esp, %ebp
    andl $-16, %esp
    subl $32, %esp
    jmp  .L4
.L5:
    movl 28(%esp), %eax
    movl %eax, (%esp)
    call fastInvSqrt
    fstpl    4(%esp)
    movl $fmt1, (%esp)
    call printf
.L4:
    leal 28(%esp), %eax
    movl %eax, 4(%esp)
    movl $fmt0, (%esp)
    call __isoc99_scanf
    cmpl $-1, %eax
    jne  .L5
    movl $0, %eax
    leave
    ret

fastInvSqrt:
    pushl    %ebp
    movl %esp, %ebp
    subl $8, %esp
    movd 8(%ebp), %xmm0
    rsqrtps  %xmm0, %xmm0
    movd %xmm0, -4(%ebp)
    flds -4(%ebp)
    leave
    ret

私はアセンブリを書くのは慣れておらず、C言語のコードをgcc -Sしたコードを元に書いてみたので変なところがあるかもしれません。 ただ注目して欲しいのはこの一行です。

    rsqrtps   %xmm0, %xmm0

この一命令によって平方根の逆数の近似値を求めています。 憶測に過ぎませんが、速さが勝負のゲームエンジンの実装では、この命令への最適化が使われていることでしょう。

ちなみにCでの実装をgcc -Sしたものは次のような感じでした。

f0:   .float    0.5
f1:  .float    1.5

fastInvSqrt:
    pushl    %ebp
    movl %esp, %ebp
    subl $16, %esp
    leal 8(%ebp), %eax
    movl (%eax), %eax
    sarl %eax
    movl %eax, %edx
    movl $1597463007, %eax
    subl %edx, %eax
    movl %eax, -8(%ebp)
    leal -8(%ebp), %eax
    movl (%eax), %eax
    movl %eax, -4(%ebp)
    flds 8(%ebp)
    flds f0
    fmulp    %st, %st(1)
    fmuls    -4(%ebp)
    fmuls    -4(%ebp)
    flds f1
    fsubp    %st, %st(1)
    fmuls    -4(%ebp)
    leave
    ret

lealしてmovlしているところが面白いところです。

まとめ

30個のプログラミング言語でFast inverse square rootというアルゴリズムを実装してみました。 コードは全て以下のリポジトリからcloneできます。 この言語がないじゃないか！とかいうツッコミがありそうですが、許してください。 SmalltalkやLua, Cleanなどは書こうとしましたが、float⇔int32の変換をどうすればよいか分かりませんでした。 AdaやDelphi, Hexe, Io, Miranda, Tcl, Racket, Pony, Forthなど、まだまだ書きたい言語は残っているので、時間があるときに書いていこうと思います。

自分が書いたことがない言語の処理系をインストールし、目的のアルゴリズムをいかに速く実装するかというのは、さながらスポーツのようです。 このエントリーのために色々調べていたら、最初のCの実装を書いてから二か月くらいかかってしまいました。 単純に30で割ると2日に1言語ですが、一時間ほどで書けるものもあれば、数日間ドキュメントを読まないと書けないものもありました。 長く孤独な戦いでした。 それでも、自分が初めて触る言語で動いた瞬間はとても嬉しく、それが何度も何度も起きるわけですから、本当にこのチャレンジは楽しかったです。

どんな言語にも、得意不得意ががあります。 どんな言語にも、良いところと悪いところがあります。 素晴らしい型システムの代償がコンパイル時間だったり、実行時パフォーマンスの代償がメモリー領域に関する脆弱性だったりします。 歴史のある言語もあれば、歴史の浅い言語もあります。 様々な言語を研究した上で作られた言語もあれば、様々な言語に多大な影響を及ぼした言語もあります。

大事なことは、目的に沿って適切な道具を選ぶことです。 仕事で新しい挑戦ができなかったとしても、プログラマーとして生きていく上で、色々な言語の特性を知っておくことは大事だと思います。 少なくとも処理系を自分のラップトップに入れてコードを実行したり、ドキュメントを探しまわったという経験があるのとないのでは、言語に対する心理的障壁がだいぶ違ってきます。

ある問題を様々な言語で書いてみるという課題は、技術を深める第一歩に過ぎません。 次のステップは、これをきっかけに興味を持った言語を更に調べて身に付けるということになるでしょう。 私はこれまでHaskell, Scala, Go, JavaScriptといった言語を書いてきましたが、今回のチャレンジで新たにRust, OCaml, Erlang, Clojureなどにも興味を持つようになりました。 これまで書いてみようと思うことすらなかった言語です。 それぞれ活躍する場面が異なっており、これらを身につけることで、より多くの場面で適切な言語選択ができるようになれたらいいなぁと思っています。

あなたも色々な言語を触ってみませんか？九九表でもいいし、Project Eulerでもいいし、数式処理系でもいいし、JSONパーサーでも、数独ソルバーでも、HTTP echo serverでもいいと思います。大きすぎない問題を選ぶことと、一つの言語に時間をかけ過ぎないこと、30程度の言語 (あるいはそれ以上) で書いてみることが大事です。そして、公式のドキュメントを見ること、きちんと処理系をインストールして実行すること、できればテストを通すこと。実装言語が増えていき、全ての言語でテストが通ることはとても楽しいことです。

さあ、問題は設定できましたか？次はあなたの番ですよ。

素敵な言語との出会いがあらんこと

そして、あなたのブログ記事をお待ちしております。

あとがき ― 比較プログラミング言語学を夢見て

以前、moznionさんという方がYAPC Asia Hachioji 2016でのトークの感想エントリーとして次の記事を書かれていました。

やられた。 そう思いました。 このエントリーを見た時点で、fastinvsqrtも15言語くらいまで進んでいて、このブログ記事の用意を進めているところでした。

しかし、冷静になってみれば、一つのアルゴリズムをいろいろな言語で書きましょうというのは定期的に誰かが言っていることですし、慌てて出すよりも、もっと言語を増やしてからしっかり言語比較しつつ出したいという気持ちが高まり、公開するまでに時間が少しかかってしまいました。 20言語を超えたあたりから進捗が悪くなり、たった10行程度のコードを書くのに数日かかったりしていました。 実装しようと思って必死になったけど結局書けなかった言語もありました。 でも、あれもこれも自分にとっていい経験になったと思います。

自然言語の研究分野の一つに、比較言語学という分野があります。 民族の移動により言語は分裂し、また文化の交流により影響しあってきました。 関連している言語の文法や音韻を比較研究することで言語のルーツを辿るのが比較言語学の主たる目的です。 私は大学時代に一般教養の科目でこの分野の講義を受けて、謎解きみたいでとてもおもしろいと感じました。

あらゆる言語で一つのアルゴリズムを書いてみるというのは、スキルの向上や言語理解よりも、もっと深い意義があるように思います。 私はこのエントリーを用意する中で、「比較プログラミング言語学」というものを妄想しました。 言語は他の言語の影響を受けて生まれ、他の言語に影響を与えながら成長し、そして廃れていきます。 構文やキーワードから、演算子や関数名など細部を比較していくと、言語計者が影響を受けた言語や、言語の設計思想があらわになってきます。

なんか偉そうなこと書いてきましたが、似たようなこと言っている人は他にいないのかなと調べてみたら、以下の記事を見つけてしまいました。

流石です、西尾先生。 当時この記事を見たかどうか覚えていませんが、こういう風に言語を観察する事の面白さに、私はようやく気が付きましたよ！ どんなに人に言われてもやろうとしなかったことでも、ふと始めたら面白くて夢中になってしまいました。 でも、人が成長するってそういうことでしょう？

エンジニア募集

はてなでは、自ら研鑽し技術を深め、あらゆる場面で適切な道具を使い分けられるようになりたい、そんな意欲的なエンジニアを募集しています。