C言語で円周率の値を使う
円周率とは3.14のことです。
本当はもっと桁数があるので、前回の投稿では次のようにしていました。
#define PI 3.141592653589793238462643383276
小数桁数が30の例です。
FreeBSDで次のように実行すれば出てきます。
$ echo 'scale=30; a(1)*4' | bc -l 3.141592653589793238462643383276
では、どうして30桁なんでしょうか。
理由は16桁以上あればいい、くらいの予備知識があったからです。
本当にそれだけで、適当に30にしました。
それなら20桁で十分ですよね、とか、
100桁や1000桁にすればいいのではという話になるわけですが、
100桁はいくら何でも多すぎるのではないでしょか。
なら何桁にすればいいんだ？
何を根拠に言ってるんだ？
ということをちゃんと考えましょう。
円周率を使いたいときはどうしたらいいのか調べて行きます。
1. 標準ではない定数に頼る
標準にこだわらないのであれば、M_PIという定数があるそうです。
これは初めて知りました。
すでにあるものを使えるなら、わざわざ桁数なんかに頭を悩ませる必要はありませんね。
使っていいならぜひぜひ使ってください。
ただ、私は標準にこだわる人なので多分使うことはないと思います。
移植性みたいなことに苦労した経験があるので、やめておこうかなと。
先にネタバレしておきますが、C言語の仕様書では浮動小数の厳密な内容を定義していませんので、
円周率の値は絶対にCPU依存の方法でしか定義できません。
標準にこだわったところでC言語汎用の円周率なんて定義できないのです。
C言語を使うってそういうことなんだろうとあきらめるしかありません。
だから、どうせ何かしら標準以外の方法に頼らなければならないなら
この方法は少しだけおススメできます。
それでは、M_PIを見ていきましょう。
いったい何が書かれているのか確認してみましょう。
まずはFreeBSDのclang環境から。
/usr/include/math.h #define M_E 2.7182818284590452354 /* e */ #define M_LOG2E 1.4426950408889634074 /* log 2e */ #define M_LOG10E 0.43429448190325182765 /* log 10e */ #define M_LN2 0.69314718055994530942 /* log e2 */ #define M_LN10 2.30258509299404568402 /* log e10 */ #define M_PI 3.14159265358979323846 /* pi */ #define M_PI_2 1.57079632679489661923 /* pi/2 */ #define M_PI_4 0.78539816339744830962 /* pi/4 */ #define M_1_PI 0.31830988618379067154 /* 1/pi */ #define M_2_PI 0.63661977236758134308 /* 2/pi */ #define M_2_SQRTPI 1.12837916709551257390 /* 2/sqrt(pi) */ #define M_SQRT2 1.41421356237309504880 /* sqrt(2) */ #define M_SQRT1_2 0.70710678118654752440 /* 1/sqrt(2) */
すげえ、いっぱいある！
続けて、Gentoo Linuxのgcc環境も見てみます。
/usr/include/math.h # define M_E 2.7182818284590452354 /* e */ # define M_LOG2E 1.4426950408889634074 /* log_2 e */ # define M_LOG10E 0.43429448190325182765 /* log_10 e */ # define M_LN2 0.69314718055994530942 /* log_e 2 */ # define M_LN10 2.30258509299404568402 /* log_e 10 */ # define M_PI 3.14159265358979323846 /* pi */ # define M_PI_2 1.57079632679489661923 /* pi/2 */ # define M_PI_4 0.78539816339744830962 /* pi/4 */ # define M_1_PI 0.31830988618379067154 /* 1/pi */ # define M_2_PI 0.63661977236758134308 /* 2/pi */ # define M_2_SQRTPI 1.12837916709551257390 /* 2/sqrt(pi) */ # define M_SQRT2 1.41421356237309504880 /* sqrt(2) */ # define M_SQRT1_2 0.70710678118654752440 /* 1/sqrt(2) */
同じものが書かれていますが、gccは上記以外にも色々と定義されています。
全部書き出しても意味がないので、M_PIを全通り抜き出してみます。
# define M_PI 3.14159265358979323846 /* pi */ # define M_PIl 3.141592653589793238462643383279502884L /* pi */ # define M_PIf16 __f16 (3.141592653589793238462643383279502884) /* pi */ # define M_PIf32 __f32 (3.141592653589793238462643383279502884) /* pi */ # define M_PIf64 __f64 (3.141592653589793238462643383279502884) /* pi */ # define M_PIf128 __f128 (3.141592653589793238462643383279502884) /* pi */ # define M_PIf32x __f32x (3.141592653589793238462643383279502884) /* pi */ # define M_PIf64x __f64x (3.141592653589793238462643383279502884) /* pi */
こんなに使えるんですね。
ちなみに、M_PIlは128bitの浮動小数で使えますよと書かれています。
どういうことなのか、これから説明していきますので覚えておきましょう。
これらの定数を使った場合の利点は、CPUのことなどまるで考慮しなくても良いことです。
たぶん、コンパイラは十分な精度を提供してくれるはず、と思い込むことができます。
非常に良い利点なのですけど、
FreeBSDではlong doubleの円周率が用意されていないので注意してくださいね。
つまり、十分な精度が用意されていません。
困ったものだ。
2. 値を自分で用意する
標準だけを使いたい場合は、自分で定義するしかありません。
ではどうやって定義していくかを見ていきます。
まずは浮動小数型に対して、10進数で表すことができるぎりぎり限界を調べます。
floatは、8桁以上。
doubleは、16桁以上。
long double (Intel 80bit)は、20桁以上。
long double (IEEE754 binary128)は、35桁以上。
ただし、上記の内容は、私が使っているCPUに合った内容です。
C言語の仕様書に載っているわけではありません。
そこはもうあきらめてください。
私のCPUがあたなのCPUと一致しているとは限りませんので、
もしあなたのCPUがIEEE 754に従ってないとか、
あるいは基数10とか16ですよとか、そういう場合は自分自身で計算して桁数を求めてください。
話を戻しましょう。
求めた桁数にすこし余裕を持たせることにします。
結果は次のようになります。
#define PI_FLOAT 3.141592654f #define PI_DOUBLE 3.14159265358979324 #define PI_LONG 3.141592653589793238462643383279502884L
注意点としては、PI_LONGの最後にLを忘れないようにしてください。
忘れるとdoubleに切り捨てられます。
3. 桁数の求め方
どうしてこの値になったのかを詳細に書いていきます。
floatはIEEE754 binary32形式、仮数23bit、hidden bitあり、合計24bit。
doubleはIEEE754 binary64形式、仮数52bit、hidden bitあり、合計53bit。
long doubleのIntel形式は、仮数64bitの、hidden bitなし、合計64bit。
long doubleのIEEE 754 binary128形式は、仮数112bit、hidden bitあり、合計113bit。
基数はすべて2です。
long doubleは、IntelのCPUだと拡張倍精度浮動小数点数と呼ばれる
80bit長の浮動小数を実装しています。
それとは別に、IEEE754ではbinary128形式という別の形式も存在します。
どちらも考慮に入れて考えてみましょう。
まとめると、
floatは、24bit
doubleは、53bit
long double(Intel)は、64bit
long double(IEEE754)は、113bit
では順番に、浮動小数の最大値が10進数で何桁かを求めます。
24bitということは、2の24乗が最大なわけで、出してみればいいんです。
* (expt 2 24) 16777216
これ、何桁かわかりますか？
わからないならこうしましょう。
* (length (princ-to-string (expt 2 24))) 8
8桁だそうです。
もっと頭良さそうに出す方法もあって、
logを使って10の何乗かを出してみます。
こうやって出します。
* (* 24 (log 2 10)) 7.22472
おおよそ7.22桁だと思ってもらえばいいと思います。
中途半端で意味不明なので、繰り上げて8桁にしましょう。
同じようにして全部出してみます。
floatは、24bit -> 7.22 (8桁)
doubleは、53bit -> 15.95 (16桁)
long double(Intel)は、64bit -> 19.27 (20桁)
long double(IEEE754)は、113bit -> 34.02 (35桁)
最低限、この桁数だけ確保すればいいわけです。
実は10進数で表すことができる最低限の桁数は
float.hに定数として設定されています。
ちょっと見てみましょう。
#include <stdio.h> #include <float.h> int main() { printf("%d\n", FLT_DIG); printf("%d\n", DBL_DIG); printf("%d\n", LDBL_DIG); return 0; }
実行結果は下記の通り。
6 15 18
これより大きな桁数であればいいわけなので、 まあまあ合っていると言えるのではないでしょうか。
それではPIに戻ります。
floatの値を出してみます。
多い分には構わないので、小数の部分を8桁+1桁の9桁にします。
実数部も1桁あるので、たぶん十分だと思います。
まず小数10桁だして、最後を四捨五入して9桁にします。
ちなみに、bcコマンドで出すと、最後の桁周辺が怪しい値になるので、
$ echo 'scale=10; a(1)*4' | bc -l 3.1415926532 ★おかしかも
みたいにはしないでください。
10桁と20桁の出力を比べてみます。
3.1415926532 3.14159265358979323844
最後の桁が2ではなく5になっているのが分かります。
ということで、いろいろやって9桁の円周率
#define PI_FLOAT 3.141592654f
が得られました。
同様にdoubleとlong doubleもやっていきます。
#define PI_FLOAT 3.141592654f #define PI_DOUBLE 3.14159265358979324 #define PI_INTEL 3.141592653589793238463L #define PI_LONG 3.141592653589793238462643383279502884L
これらの値が正しく格納されているかどうかは、 浮動小数を16進数で出力すると分かりやすいです。
#include <stdio.h> #define PI_FLOAT 3.141592654f #define PI_DOUBLE 3.14159265358979324 #define PI_INTEL 3.141592653589793238463L #define PI_LONG 3.141592653589793238462643383279502884L int main() { float x; double y; long double z, w; x = PI_FLOAT; y = PI_DOUBLE; z = PI_INTEL; w = PI_LONG; printf("%23.20A\n", x); printf("%23.20A\n", y); printf("%23.20LA\n", z); printf("%23.20LA\n", w); return 0; }
たしかprintfの%Aはc99の機能だったと思います。
実行結果は下記の通り。
0X1.921FB600000000000000P+1 0X1.921FB54442D180000000P+1 0X1.921FB54442D1846A0000P+1 0X1.921FB54442D1846A0000P+1
よさそうですね。
一応補足しますが、実行したPCのCPUはIntel製のため、
IEE754 binary128は確認できません。
あとちょっと話題を広げますが、PI_INTELとPI_LONGって二つも必要でしょうか。
IEEE754のbinary 128というは、
CPUもプログラミング言語もほとんど対応してないようで、
めったにお目にかかれないそうです。
（とはいえgccでは非標準で使えます。mpfrの機能なんでしょうね。）
あまり存在しないものに合わせてもしょうがないので、
PI_INTELの方をlong doubleにしても誰も文句言わないと思います。
しかし桁数は多いに越したことはないので、
わざわざ桁数の少ないPI_INTELを採用する必要もありません。
PI_LONGだけにしましょう。
ということで削除すると次の結果になります。
#define PI_FLOAT 3.141592654f #define PI_DOUBLE 3.14159265358979324 #define PI_LONG 3.141592653589793238462643383279502884L
これで完了です！
ちなみに、私のお気に入りはこれ。
/* $ echo 'scale=40; a(1)*4' | bc -l */ #define PI_FLOAT 3.1415926535897932384626433832795028841968f #define PI_DOUBLE 3.1415926535897932384626433832795028841968 #define PI_LONG 3.1415926535897932384626433832795028841968L
適当に40桁がいちばん楽ですね。
だったら何のためにここまで詳しく語ってきたんだ。
4. 定数の型
適当にしていても何とかなるのが「型」です。
定数の型とは、例えばこんな感じ
3.14はdouble型3.14fはfloat型3.14Lはlong double型
標準ではないM_PIという定数はdouble型です。
C言語には暗黙のキャスト機能がありますので、floatでも使えます。
しかし、long doubleでは使えません。
正確には、使えるのですけど精度が十分ではありません。
それなら初めからlong double型で定数を宣言しておけばいいと思いませんか？
つまり、M_PIなんていらずにM_PIlだけにすればいいのでは？
正しいとも間違ってるとも言えます。
例ですがCommon Lispの定数piは、そういう考えで最大の長さを持つ型の値が格納されています。
見てみましょう。
* pi 3.14159265358979324L0
値がlong-floatになっているのが分かります。
sbclとcclのpiはdouble-float型で格納されていますが、
どちらもそもそもlong-floatという型そのものが存在せず、
double-floatとして扱われます。
C言語もlong double型のPIを一つだけ用意しておけばいいのでは？
それをしない理由は、定数の型によって出力されるコードが変わる可能性があるためです。
ここでは話題にしませんが、当然Common Lispも同じです。
浮動小数の演算は、型が大きい方に自動的にキャストされます。
例えばdoubleとfloatの足し算は、doubleです。
floatとlong doubleの掛け算は、long doubleです。
この暗黙のキャストを、次の例で見てみます。
double multi_long(double x) { return x * PI_LONG; }
関数multi_longは、引数にPIを乗算する関数です。
定数はPI_LONGを使っているので、long doubleで演算されます。
その様子を見てみましょう。
実行例はFreeBSD, amd64, clang, gdbです。
$ cc -g main.c $ gdb a.out (gdb) disassemble multi_long Dump of assembler code for function multi_long: 0x00000000002019e0 <+0>: push %rbp 0x00000000002019e1 <+1>: mov %rsp,%rbp 0x00000000002019e4 <+4>: movsd %xmm0,-0x8(%rbp) 0x00000000002019e9 <+9>: fldl -0x8(%rbp) ★ここを見てほしい★ 0x00000000002019ec <+12>: fldt -0x14d2(%rip) # 0x200520 0x00000000002019f2 <+18>: fmulp %st,%st(1) 0x00000000002019f4 <+20>: fstpl -0x10(%rbp) 0x00000000002019f7 <+23>: movsd -0x10(%rbp),%xmm0 0x00000000002019fc <+28>: pop %rbp 0x00000000002019fd <+29>: ret End of assembler dump.
「★ここを見てほしい★」では、
fldl命令からx87 FPUのSTレジスタが使われているのが分かります。
つまりlong doubleで演算されています。
それでは、doubleの定数を乗算する方も見ていきます。
double multi_double(double x) { return x * PI_DOUBLE; }
逆アセンブルしてみます。
(gdb) disassemble multi_double Dump of assembler code for function multi_double: 0x0000000000201a00 <+0>: push %rbp 0x0000000000201a01 <+1>: mov %rsp,%rbp 0x0000000000201a04 <+4>: movsd -0x1504(%rip),%xmm1 # 0x200508 0x0000000000201a0c <+12>: movsd %xmm0,-0x8(%rbp) 0x0000000000201a11 <+17>: mulsd -0x8(%rbp),%xmm1 0x0000000000201a16 <+22>: movaps %xmm1,%xmm0 0x0000000000201a19 <+25>: pop %rbp 0x0000000000201a1a <+26>: ret End of assembler dump.
multi_longとはちがい、STレジスタではなくxmmレジスタが使われており、
doubleで乗算されているのが分かります。
定数の型をちゃんと意識しないと、
doubleで済むところをlong doubleで演算してしまうということです。
でもそれの何が問題なの？ という話ですが、難しいですね。
何も問題ないかもしれません。
はたしてfloatとdoubleとlong doubleは、
演算の速度に差があるのでしょうか？
その答えは、当たり前なのですがCPU依存です。
大昔にこんなことを聞いたことがあります。
「C言語はdouble使ってりゃいいんだよ。」
たぶんIntelのCPUの話ではないと思うのですが、
CPUによってはdoubleしか演算できず、
floatの場合はわざわざdoubleにキャストしてから演算し、
結果をfloatにキャストして戻すんだそうです。
手間がかかるため、doubleの方が圧倒的に早かったそうです。
amd64の逆アセンブルの例でも、doubleとlong doubleで出力するコードが違っていました。
doubleはxmmレジスタ。
long doubleはSTレジスタです。
演算する場所すら違うので、速度に差が生じる可能性はあります。
確認はしてみたのですが、私の環境でテストした限りでは
multi_doubleとmulti_longの演算時間に差は全くありませんでした。
でも、そんなに簡単な話でもないとは思います。
せっかくなのでfloatの例も見てみましょう。
double multi_float(double x) { return x * PI_FLOAT; }
(gdb) disassemble multi_float Dump of assembler code for function multi_float: 0x00000000002019c0 <+0>: push %rbp 0x00000000002019c1 <+1>: mov %rsp,%rbp 0x00000000002019c4 <+4>: movsd -0x14bc(%rip),%xmm1 # 0x200510 0x00000000002019cc <+12>: movsd %xmm0,-0x8(%rbp) 0x00000000002019d1 <+17>: mulsd -0x8(%rbp),%xmm1 0x00000000002019d6 <+22>: movaps %xmm1,%xmm0 0x00000000002019d9 <+25>: pop %rbp 0x00000000002019da <+26>: ret End of assembler dump.
面白いことに、multi_floatとmulti_doubleは全く同一です。
PI_FLOATをdoubleの定数として格納しているんでしょうね。
定数の型がfloatとdoubleのどちらでもdoubleとして扱われるため、
速度差に影響することはなさそうです。
まあdoubleの演算にfloatの定数を入れたって、
いいこと一つもありませんけどね。
念のため注意書きしておきますが、あくまでも私の実験での話です。
一般的な話ではございません。
適当にまとめます。
何も考えたくないなら、円周率の定数はdouble型を使えばいいと思います。
でも速度が要求されるような場合は、
浮動小数の型を全体的に考えた方がいいかもしれません。
自分が使おうとしているCPUの仕様を少しだけでも調べておいて、
必要に応じて逆アセンブルしてみて下さい。
本当に適当にまとめましたが以上です。
