はじめに

少し前に seq コマンドを awk を使って実装してみました。表面上はシェルスクリプトになっていますが、シェル関数として使えるようにラッパーにしているだけでオプション解析も含めて全て awk で実装しています。

実装した理由は seq コマンドを POSIX で標準化するのは可能だろうか？という研究のためです。元々 seq コマンドは 1985 年の Version 8 Unix で初めて実装されたものであり、1995 年に GNU で再実装されました。BSD Unix では実装が遅れましたが現在ではすべての BSD 系 Unix で実装されています。Solaris でも GNU 版が標準で組み込まれており、今なら seq コマンドは POSIX で標準化が可能になったのではないか？という検証のためです。

この記事ではこれ以上 POSIX での標準化の話はしませんが、この検証作業の一環として seq コマンドを実験的に実装してみて、小数点の誤差の扱いはめんどくせぇと思った点、誤差をなくすにはどうするか、ちゃんと動くものを作るのは簡単なことではないと思った話をしたいと思います。

seq コマンドの最低限の実装なら簡単

seq コマンドを知らない人もいるだろうということで簡単に説明すると、seq コマンドは連続する数値を出力するコマンドで、開始値から終了値まで指定したステップ毎に出力するコマンドです。引数は省略できますが、この記事では3つの引数をすべて指定しています。

$ seq 1 2 10 # seq 開始値 ステップ値 終了値 （1から10まで2つ置きに出力）
1
3
5
7
9

さて、小数に対応していない seq コマンドを実装するのは簡単です。補足ですが %g は macOS 版の seq がデフォルト（-f で書式を指定したなかったとき）で使用する書式です。ちなみに GNU 版の seq は少し複雑な方法でデフォルトの書式を決めています。

seq.awk

BEGIN {
  first = ARGV[1]; incr = ARGV[2]; last = ARGV[3];
  for (i = first; i <= last; i += incr) {
    printf "%g\n", i
  }
}

$ awk -f seq.awk 1 2 10
1
3
5
7
9

シェルスクリプトだけで作ることも簡単です。

seq.sh

#!/bin/sh
set -eu

seq() {
  first=$1 incr=$2 last=$3
  i=$first
  while [ "$i" -le "$last" ]; do
    printf '%g\n' "$i"
    i=$((i + incr))
  done
}

seq 1 2 10

もっと短く書くとこうなります。

seq.sh

# 簡略化してローカル変数をなくしたバージョン
seq() {
  while [ "$1" -le "$3" ]; do
    printf '%g\n' "$1"
    set -- $(($1 + $2)) "$2" "$3"
  done
}


# 補足: もちろん bash 4.0 以降ならこれで良い
printf "%g\n" {1..10..2}

ただし小数の計算は（ksh、yash、zsh を除き）POSIX シェルの機能の範囲では実現できません。seq を awk で実装した理由はシェルが小数に対応していないからです。

小数を扱う場合はこれでは上手く行かない

先程の awk の実装は、awk が小数にも対応しているため、小数を指定しても動きます。

$ awk -f seq.awk 0 0.1 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1

コードのメイン部分の実装はこれでよし、さて seq コマンドのさまざまなオプションを実装してみるかーと考えていましたが、いくつかテストをしているとすぐにこれでは動かないことに気づきました。

$ awk -f seq.awk 1 0.1 2
1
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
　　　　　👈 2 が出力されてない！

まあね、考えてみれば当然の話ですね。awk は10進数の小数を2進数の浮動小数点数として扱っています。

なぜ最後の2は出力されなかったのか？

その理由は printf のフォーマット書式を %g から %.30g に変更して小数点以下の数値の桁を増やすとすぐに分かります。

BEGIN {
  first = ARGV[1]; incr = ARGV[2]; last = ARGV[3];
  for (i = first; i <= last; i += incr) {
    printf "%g (%.30g)\n", i, i
  }
}

$ awk -f seq.awk 1 0.1 2
1 (1)
1.1 (1.10000000000000008881784197001)
1.2 (1.20000000000000017763568394003)
1.3 (1.30000000000000026645352591004)
1.4 (1.40000000000000035527136788005)
1.5 (1.50000000000000044408920985006)
1.6 (1.60000000000000053290705182008)
1.7 (1.70000000000000062172489379009)
1.8 (1.8000000000000007105427357601)
1.9 (1.90000000000000079936057773011)
          👈 最後の出力されない 2 は 2 よりも少し大きな数になっている
             (2.00000000000000088817841970013)

左側の数値が %g で出力された数値で、右側のカッコ内の数値が小数点以下含めて30桁まで出力した結果です。ここから目で見えていた数値（%g で出力された数値）は正確な値ではなく、内部的にはわずかに大きい数値として扱われていた事がわかります。そして最後 2 は 2 よりも少し大きな数になってしまいます。ループの終了条件は i <= last (last = 2) なので、最後の 2よりも少し大きな数はループ条件を満たさないので出力されません。

しかしここで一つ不思議なことがあります。0 から 0.1 ステップで 1 まで出力した場合は、最後の 1 が出力されたということです。その理由も 30 桁まで出力するとわかります。

$ awk -f seq.awk 0 0.1 1
0 (0)
0.1 (0.100000000000000005551115123126)
0.2 (0.200000000000000011102230246252)
0.3 (0.300000000000000044408920985006)
0.4 (0.400000000000000022204460492503)
0.5 (0.5)
0.6 (0.599999999999999977795539507497)
0.7 (0.699999999999999955591079014994)
0.8 (0.799999999999999933386618522491)
0.9 (0.899999999999999911182158029987)
1 (0.999999999999999888977697537484)     👈 最後の 1 は 1 よりも少し小さな数

この結果から小数の値は必ずしも大きい数値に丸められるとは限らないということがわかります。

$ awk 'BEGIN { printf "%.30g\n", 0.4 }' # 0.4 は 0.4より少し大きい値に丸められる
0.400000000000000022204460492503

$ awk 'BEGIN { printf "%.30g\n", 0.5 }' # 0.5 は 0.5 と完全に同じ
0.5

$ awk 'BEGIN { printf "%.30g\n", 0.6 }' # 0.6 は 0.6 より少し小さい値に丸められる
0.599999999999999977795539507497

ここで 0.5 の出力に誤差がないのは 0.5 が2進数の浮動小数点数でも正確に扱える値だからです。2 進数を 10 進数で表す場合次のような表のイメージで変換します（きちんとした計算方法は他の記事を参照してください）。

2³	2²	2¹	2⁰	2⁻¹	2⁻²	2⁻³
8	4	2	1	0.5	0.25	0.125
<---------------- ２倍にする ------------------ ------------------ 半分にする ---------------->

たとえば 2 進数で 1010 の場合、2³ + 2¹ = 8 + 2 で 10 です。同じように 2進数で 0.101 は 2⁻¹ + 2⁻³ = 0.5 + 0.125 で 0.625 です。合わせると 1010.101 は 10.625 になります。

2³	2²	2¹	2⁰	2⁻¹	2⁻²	2⁻³
8	4	2	1	0.5	0.25	0.125
1	0	1	0 .	1	0	1	→ 1010.101 = 10.625

0.5 は 2⁻¹ なので、2 進数で 0.1 となります。ではこのルールにしたがって 0.6 は 2 進数でどう表現すればよいでしょうか？答えは「表現できません」です。正確には 0.100110011... と小数点以下が無限に続くため有限の桁数しか扱えないコンピュータで浮動小数点数を使って表現しようとする限り途中で打ち切らないといけないため表現できません。もし 10 進数の小数値を正確に表現したい場合は、内部的にも 10 進数で表現する必要があります。

さて、0.5 は 2進数で表現できる、0.6 は 2進数で表現できない、と説明しましたが、1.0 や 2.0 はどうでしょうか？答えは「表現できます」です。小数点以下がないのですから当然です。

$ awk 'BEGIN { printf "%.30g\n", 1.0 }'
1

$ awk 'BEGIN { printf "%.30g\n", 2.0 }'
2

しかしここまでの話では 1.0 や 2.0 にも誤差が含まれていました。つまりこれは単純な小数の誤差とは別の問題が含まれているということを意味しており一種のバグです。

小数値の誤差を積み重ねてはいけない

さて、それでは awk 実装版の seq コマンドをより正確に出力するにはどうしたら良いでしょうか？ Deciaml 型などの 10 進数を使う？たしかにそれも一つの手なのですが、awk に Deciaml 型や 10進数ライブラリがない（探せば見つかるかもしれませんが）ので簡単にはできません。これは awk を使って消費税の計算などを行う場合には誤差に注意が必要だということです。たいていは 1 円程度の誤差しか出なかったりするのですが、銀行レベルでは 1 円の誤差も許されないでしょう。

seq コマンドの場合、小数値の加算しか行わず掛け算や割り算がないので、10進数として計算するのはさほど難しくはありませんが、ただでさえ awk 実装で遅いのがさらに遅くなってしまいます。なによりその他の seq コマンドも 2 進数の浮動小数点数で処理を行っています。以下のように seq コマンドでフォーマットを指定すると、seq コマンドも内部的に 2 進数で扱っているのだろうということがわかります。それなら awk で実装しても同等の出力を得ることは可能なはずです。

macOS 版 seq と自作 seq の出力の違い

$ seq -f '%.30g' 1 0.1 2           $ awk -f seq.awk 1 0.1 2
1                                  1
1.10000000000000008881784197001    1.10000000000000008881784197001
1.19999999999999995559107901499    1.20000000000000017763568394003
1.30000000000000004440892098501    1.30000000000000026645352591004
1.39999999999999991118215802999    1.40000000000000035527136788005
1.5                                1.50000000000000044408920985006
1.60000000000000008881784197001    1.60000000000000053290705182008
1.70000000000000017763568394003    1.70000000000000062172489379009
1.80000000000000004440892098501    1.8000000000000007105427357601 
1.89999999999999991118215802999    1.90000000000000079936057773011
2                                  なし

どちらも小数値の誤差がありますがその誤差の値が異なります。1.1 は誤差が同じですが、よく見ると誤差がどんどん大きくなっています。そして 2 進数で正しく扱えるはずの、1.0、1.5、2.0 でさえ自作 seq では誤差が発生しています。

最初に書いた、大きな誤差が発生してしまう seq の awk 実装コードを再掲します。

大きな誤差が発生してしまうseqの実装

BEGIN {
  first = ARGV[1]; incr = ARGV[2]; last = ARGV[3];
  for (i = first; i <= last; i += incr) {
    printf "%.30g\n", i, i
  }
}

感のいい人ならすでに気づいていると思いますが、このコードの問題点は変数 iに誤差が含まれる小数値 (0.1) を足し込んでいることです。0.1 は 2 進数で表現すると 0.100000000000000005551115123126 と 0.1 よりも少し大きい値になってしまうため、ループを繰り返すたびに誤差を積み重ねてしまい、本来誤差が発生しないはずの 1.5 でさえ誤差が含まれてしまいます。つまり、誤差を大きくさせないためには i に小数値を足し込むようなコードを書いてはいけないということになります。

ではどうすればいいか？ここで「ループ回数」に注目します。1 から 0.1 を加えながら 2 になるまで出力するということは、合計で 11 回のループを行うということです。つまり 1 から 2 になるまで繰り返す代わりに 11 回ループするという形になるように変換し、その都度計算を行えば良いということです。開始値と終了値とステップ値が分かればループ回数は計算で求めることができます。ループ回数に小数値は登場しないため原理的に誤差の積み重ねは発生しません。

誤差の発生を小さくしたseqの実装

BEGIN {
  first = ARGV[1]; incr = ARGV[2]; last = ARGV[3];
  n = int((last - first) / incr) # ループ回数 = int((2 - 1) / 0.1)
  for (i = 0; i <= n; i++) {
    printf "%.30g\n", first + (i * incr)
  }
}

ちなみにこの話は「浮動小数点変数をループカウンタに使用しない」という話で知られています。

誤差を積み重ねることをやめたため出力結果が改善されました。

macOS 版 seq（左）と改善された自作 seq の出力の違い

$ seq -f '%.30g' 1 0.1 2           $ awk -f seq.awk 1 0.1 2
1                                  1
1.10000000000000008881784197001    1.10000000000000008881784197001
1.19999999999999995559107901499    1.19999999999999995559107901499
1.30000000000000004440892098501    1.30000000000000004440892098501
1.39999999999999991118215802999    1.39999999999999991118215802999
1.5                                1.5
1.60000000000000008881784197001    1.60000000000000008881784197001
1.70000000000000017763568394003    1.70000000000000017763568394003
1.80000000000000004440892098501    1.80000000000000004440892098501
1.89999999999999991118215802999    1.89999999999999991118215802999
2                                  2

素晴らしい。完全一致です。と言いたいところなのですが、実は GNU 版 seq は違った結果になったりします。

macOS 版 seq（左）と GNU 版 seq（右）との比較

$ seq -f '%.30g' 1 0.1 2
1                                  1
1.10000000000000008881784197001    1.10000000000000000002168404345 
1.19999999999999995559107901499    1.2000000000000000000433680869
1.30000000000000004440892098501    1.2999999999999999999566319131
1.39999999999999991118215802999    1.39999999999999999997831595655
1.5                                1.5
1.60000000000000008881784197001    1.60000000000000000002168404345
1.70000000000000017763568394003    1.7000000000000000000433680869
1.80000000000000004440892098501    1.7999999999999999999566319131
1.89999999999999991118215802999    1.9000000000000000000867361738
2                                  2

今回の私の検証（seq コマンドは POSIX で標準化が可能になったのではないか？）とはテーマが異なるので深追いはしていませんが、どうやら以下の検証から GNU seq は倍精度 (double) よりも高い精度を持つ拡張精度 (long double) を使っているようです。よく見ると GNU 版 seq の方が本来の値に対して誤差が小さく（つまりより正確に）なっているように見えます。

参考として Linux と macOS で雑に調べた結果です。

Linux での C 言語の出力結果

#include <stdio.h>

int main(void) {
  printf("%.30f\n", 1.7);
  // => 1.699999999999999955591079014994
  printf("%.30g\n", 1.7);
  // => 1.69999999999999995559107901499
  printf("%.30Lg\n", 1.7l);
  // => 1.7000000000000000000433680869
  return 0;
}

Linux での出力

$ gawk 'BEGIN { printf "%.30g\n" ,1.7 }'
=> 1.69999999999999995559107901499
$ gawk 'BEGIN { printf "%.30Lg\n" ,1.7l }' # エラーにならないが結果も変わらない
=> 1.69999999999999995559107901499

$ mawk 'BEGIN { printf "%.30g\n", 1.7 }'
=> 1.69999999999999995559107901499
$ original-awk 'BEGIN { printf "%.30g\n", 1.7 }'
=> 1.69999999999999995559107901499

$ printf "%.30g\n" 1.7
=> 1.7000000000000000000433680869
$ seq -f "%.30g" 1.7 1.7
=> 1.7000000000000000000433680869

Ubuntu 18.04.6 on WSL1 での謎現象

$ seq -f "%.30g" 1.6 0.1 1.6
1.60000000000000000002168404345
$ seq -f "%.30g" 1.7 0.1 1.7
1.7000000000000000000433680869

$ seq -f "%.30g" 1.6 0.1 1.7
1.60000000000000000002168404345
1.69999999999999995559107901499
👆 なんでループしたら異なる結果になるのさ？
Ubuntu 18.04.6 on docker では発生しない

macOS での出力

$ awk 'BEGIN { printf "%.30g\n", 1.7 }'
=> 1.69999999999999995559107901499
$ printf "%.30g\n" 1.7
=> 1.69999999999999995559107901499
$ seq -f "%.30g" 1.7 1.7
=> 1.69999999999999995559107901499

$ gawk 'BEGIN { printf "%.30g\n", 1.7 }'
=> 1.69999999999999995559107901499
$ gprintf "%.30g\n" 1.7
=> 1.7000000000000000000433680869
$ gseq -f "%.30g" 1.7 1.7
=> 1.7000000000000000000433680869

さいごに

10 進数の小数を 2 進数で表すと大抵の場合、誤差が発生します。誤差がでないようにするには 10 進数で扱う方法がありますが、必ずしもその方法が正解というわけではなく 2 進数の浮動小数点数として扱わなければいけない場合もあります。その場合のテクニックとして他にも数値が同じであるかを比較する際にイコールではなくある程度の誤差を許容して比較する方法もありますが、そのテクニックが最善とは限りません。一番の方法は誤差があることを前提に、その誤差が大きくならないようなアルゴリズムを使うことではないでしょうか。誤差があることを前提とすると、整数だけを考慮すればよい場合の普通のやり方とは変わる場合があります。

小数の誤差の話はプログラミングの基礎なのですが、わかっていても普段から使っていなければ失念するものです（言い訳）。この記事では 0.1 刻みで 0 から 1 の数値を出力する場合には問題が発生しませんでしたが、1 から 2 の数値を出力する場合には問題が発生しました。テストが甘いとこのようなバグを見逃してしまいます。しかし小数の誤差の話を知らなければ十分なテストはできないでしょう。簡単に思える seq コマンドの実装でも油断すると簡単にバグを入れてしまいます。seq コマンド程度簡単に作れると甘く見ないことです。バグが有る実装を長年使っていて後で気づいたときは悲惨です。プログラミングの基礎をしっかり学んだ上で十分なテストを行うことが重要です。