📌

printfの再実装をしたたかにやってみる(1)

に公開

C学習の定番

C言語の学習におけるありがちなテーマの1つにprintfの再実装があります。
Cを触る者であれば、テキスト処理のみならず、デバッグ用途などでもprintfをガシガシ入れ込むことは日常茶飯事ですし、その実装を自ら書いてみるというのはなかなか楽しく勉強になるものです。

実際、ネット上には様々な学習者たちによるprintfの実装体験が見受けられます。
それらの多くはprintfの軽量版の実装例を示して、こんな感じに書けるよ!と読み手のハードルをやさしく下げてくれています。

私自身もこれまで何回も軽めのprintfみたいな関数を書いてきました。簡易なものでは%s, %c, %d/i, %u, %x/X,%p/Pあたりに対応するもので、フィールド幅/精度に未対応であればちょっとした空き時間に書ける程度のコード量かと思います。フィールド幅/精度と各種プリフィックスはそれなりにややこしい仕様となっており、この変換指定子にこのプリフィックスは無効だとか、このプリフィックスがあるとこちらは無効になるとか、気にすべき項目が増えてきます。

さらに浮動小数点数に対応するとなるとコード量は倍以上になります。3年ほど前にIEEE 754の文字列変換コードを書いたときには結構手間取りました。無知だったからというのが大きいですが。

今回、mallocの再実装を行う過程でデバッグ出力にprintfを使うとその中でmallocを呼び出してしまってよろしくない、という事象が発生したため、mallocを一切使わないprintfが必要になりました。
必要な仕様としては%pが出せたり、size_tあたりに素直に対応できる程度でよかったのですが、実装しているうちに欲が出てしまい、あれもこれもと手を伸ばすことになりました。いやー楽しい。現在のところ以下の変換指定子/フラグ/長さ修飾子に対応しています。

フラグ文字
- #
- 0
- -
- ' ' (スペース)
- +

対応変換指定子
- d,i,o,u,x,X
- e,E
- f,F
- g,G
- a,A
- c
- s
- p,P
- %
- b (独自拡張: 2進数表示)

対応長さ修飾子
- hh
- h
- l
- ll
- L
- z

再実装のルール

再実装にあたって使用する関数はwriteなどの環境上どうしても必須となるシステムコールのみとしています。なのでstrlenを始めstring系の関数も自身で実装します。なんでそんな面倒なと思うかもしれませんが、そもそも楽をしたいならprintfなんてもの書く必要ありません。自作したいから自作するわけです。「Dr. STONE」が好きなのです。
とはいえ、このシリーズでは掲載コード量を少なくするためstrncpyあたりの全ての関数をいちいち再実装してみせることはしません。

浮動小数点数

なにはともあれ一番おもしろいのはdoubleやlong doubleの文字列化でしょう。整数の文字化は実はそれほど難しいことはないですし、文字列出力はポインタもらってstrlenしてそのまま出力するだけ。一方、浮動小数点数は値を取得したところでdouble型はそう簡単に文字列化できません。

現代の一般的なPCで利用される浮動小数点数はIEEE 754が実質的な標準となっています。C規格では浮動小数点数の形式に指定はないものの(オプションとしてはある)、まだIEEE 754がドラフトだった段階でハードウェアとして実装されたx87コプロセッサのIBM PCなどへの搭載とその商業的な大成功によって実現したデファクトスタンダードがIEEE 754(および国際標準規格ISO/IEC 60559)です。というわけで、現代の一般的な環境を前提とするならばprintfで出力すべき浮動小数点数はIEEE 754形式となります。なお、現在のIEEE 754には10進浮動小数点数形式も存在しますが、こちらはまだまだ一般的ではありませんので対応しません。以降、IEEE 754形式といったら2進浮動小数点数形式と読み替えてください。

IEEE 754形式のデコード実装についてはそれこそQiitaあたりに記事がたくさんあるだろうからリンク貼って省略、と思って検索したらあまりない。

以下に(簡単のため)float型の変換についてコード付きで説明されているものがあるくらい?
https://qiita.com/daitai-daidai/items/44632da01d4b0a18c33e

printfの浮動小数点数の引数はfloatではなくdoubleであり、floatを渡してもdoubleに昇格します。ですのでprintfの再実装という目的から今回はfloatガン無視でdoubleに対応することになりますが、IEEE 754の構造としてはfloatでもdoubleでも同じなので上記の記事は参考になるはずです(自分が取った処理方法とは少し違いましたがとっつきやすそうです)。

倍精度浮動小数点数

floatは32ビット幅のところ、doubleでは64ビット幅となります。倍精度と言われる所以ですね。

符号: 1ビット
指数部: 11ビット
仮数部: 52ビット

図を用意するのが面倒なのでwikipediaを貼っておきます。
https://ja.wikipedia.org/wiki/倍精度浮動小数点数

IEEE 754の浮動小数点数には、正規化数、非正規化数(subnormal numbers)、ゼロ、無限大(±∞)、NaNが存在します。指数部により非常に広いダイナミックレンジを持ち、その主要な範囲を正規化数がカバーします。倍精度であれば以下の範囲です。

正規化数の範囲(正)

2^{-1022}〜2^{1023}

正規化数

doubleの正規化数の場合、

(-1)^{\text{符号部}} \times 2^{\text{指数部} - 1023} \times (1.\text{仮数部})

となります。
ここでの仮数部52ビットは当然2進数の小数部。指数部も2進数なので、double型の文字列化に際してまず行うべきことは、それぞれの2進数→10進数化です。

指数部と仮数部を取り出すのはシフトするなり、ビットフィールドで取り出すなり(余計なパディングされないように注意が必要ですが)すればシンプルに実現します。次に、取り出した仮数部を10進数にしましょう。

正規化数では52ビットの最上位にさらに1を加えて53ビットとし、新たに加えた1のビットを1 * 2^0、その下のビット * 2^-1、さらにその下のビット * 2^-2、さらにさらにその下のビット * 2^-3としてすべてのビットについて計算して合計します。

たとえば仮数部(に暗黙の1ビットを先頭に足して)11010000....(残りは0)なら

1 * 2^0 = 1 * 1 = 1
1 * 2^-1 = 1 * 0.5 = 0.5
0 * 2^-2 = 0 * 0.25 = 0
1 * 2^-3 = 1 * 0.125 = 0.125

ですから、合計は1.625になります。とてもシンプルです。ですが、この計算をするための容れ物がありません。doubleを文字列に変換するためにdoubleを使うことはできないのです(そもそもそれが可能ならビット列を眺める必要はありませんね)。それなら整数として扱えばいいのでは?とまずは誰しも考えるでしょう。しかしこれに見合う容れ物=型もみつかりません。なぜなら2の負の累乗は1進むごとに1桁増えるからです。52ビットあると52桁です。そんな大きな桁を扱える型はCには存在しません。さらに指数部があります。指数部は「11ビットしかない」ように見えますが、いえいえ「11ビットも」あります。

11ビットで表現可能な最大値は2047ですが、IEEE 754の倍精度では1023のバイアスを差し引き、正規化数では−1022から1023までの範囲を取り得ます。負の指数が1小さくなるごとに1桁増加ですから、つまり2^-1022!というおそろしく桁の長い数もあるのです。

たとえばDBL_MINなら以下のようなビット列になります。

DBL_MIN: 0b0000000000010000000000000000000000000000000000000000000000000000

符号が0で正、指数部が1、仮数部が0です。バイアスを引くと指数は−1022であり、暗黙の1(ケチ表現)を先頭に付加した仮数部は1.0を意味しますので、十分に大きな精度を指定した場合のprintf("%f", DBL_MIN)の出力結果は以下のようになります(末尾の0群は除去。以下同)。

0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000002225073858507201383090232717332404064219215980462331830553327416887204434813918195854283159012511020564067339731035811005152434161553460108856012385377718821130777993532002330479610147442583636071921565046942503734208375250806650616658158948720491179968591639648500635908770118304874799780887753749949451580451605050915399856582470818645113537935804992115981085766051992433352114352390148795699609591288891602992641511063466313393663477586513029371762047325631781485664350872122828637642044846811407613911477062801689853244110024161447421618567166150540154285084716752901903161322778896729707373123334086988983175067838846926092773977972858659654941091369095406136467568702398678315290680984617210924625396728515625

整数部を入れると1023桁。文字という意味ではドットをいれて1024文字あります。

非正規化数

しかもですよ、なんとこれが最大長ではありません! DBL_MIN / 2やDBL_MIN / 4のようにDBL_MINを下回る値もあります。この領域は正規化数ではなく非正規化数と呼ばれます。

もし非正規化数が存在しなければ、正規化数の正の最小値2^{-1022}を下回ると0に丸められてしまいます。これを避けるために作られたのが非正規化数です。 非正規化数は、正規化数の最小絶対値とゼロの間の範囲を扱います。

正規化数の正の最小値 

2^{-1022}

非正規化数の範囲(正) 
2^{-1074}〜2^{-1022}

仮数部に暗黙の1を加えることはしません。非正規化数の指数部のビット表現は0。よって、double型における非正規化数の正の最小値DBL_TRUE_MINは仮数部のLSBが1になります。

DBL_TRUE_MIN: 0b0000000000000000000000000000000000000000000000000000000000000001

非正規化数は以下の式によって表します。

(-1)^{\text{符号部}} \times 2^{1-1023} \times (0.{\text{仮数部}})

指数の項が1-1023になっており、実際の指数部が0であることに対して1が当てはめられているのが気になりますが、これは

(-1)^{\text{符号部}} \times 2^{0-1023} \times ({\text{仮数部52ビット}})

と同じです。この2つの式の違いは、仮数部の最上位に0のビットを加えた53ビットとして正規化数との処理の共通化を図るために1ビット分の調整をしているかしていないかです。1つ目の式はソフトウェア処理よりもむしろFPUとしてハードウェア実装することを念頭に置いて書かれた式と言えるでしょう。

話を正の最小の値であるDBL_TRUE_MINに戻します。
2^{(1-1023)} \times 2^{\text{-52}} = 2^{-1074}となるこの値のprintfの出力は以下の通りです。桁数は1022にさらに52を足した1074に整数部を含めて1075桁。文字としては小数点を含めた1076文字になります。

0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000004940656458412465441765687928682213723650598026143247644255856825006755072702087518652998363616359923797965646954457177309266567103559397963987747960107818781263007131903114045278458171678489821036887186360569987307230500063874091535649843873124733972731696151400317153853980741262385655911710266585566867681870395603106249319452715914924553293054565444011274801297099995419319894090804165633245247571478690147267801593552386115501348035264934720193790268107107491703332226844753335720832431936092382893458368060106011506169809753078342277318329247904982524730776375927247874656084778203734469699533647017972677717585125660551199131504891101451037862738167250955837389733598993664809941164205702637090279242767544565229087538682506419718265533447265625

当然ながらDBL_MINにDBL_TRUE_MINを加算した正規化数もまた同じ桁数を持ちます。

なお、DBL_TRUE_MIN(2^{-1074})のprintf("%.1074f")による出力桁数と有効桁数は一致しません。DBL_TRUE_MINの仮数部は最下位ビットのみが1で、2進数の有効桁数は1ビットです。このため、10進数変換後(約4.94 × 10^{-324})の最初の非ゼロ桁「4」は有効桁ですが、後続の桁(9, 4, 0, …)については信頼できません。言い換えると、4.94 × 10^{-324}4.95 × 10^{-324}の間の値はIEEE 754倍精度では表現できず、次の値は9.88 × 10^{-324}です。一方、printfの出力はフォーマット指定に依存していて有効桁数とは独立です。printfによって有効桁数を超える桁が表示されるのは、10進数変換の完全性を保つためです。printfを再実装するなら有効桁数に関係なくこれらの数値文字列を正しく出力する必要があります。つまり極めて長い数値の演算が必要になります。というわけで多倍長整数演算の登場です。

長くなりそうなので続きは別記事で。

Discussion