対数関数( log )の数値的解法
今回は対数の数値的解法を整理しておこうと思います。まずは自然対数です。
テーラー展開での解法となりますが、ただ \(log(x)\) では０近傍でのテーラー展開ができませんので、工夫のためとりあえず\(log(1+x)\)及び、\(log(1-x)\)を展開します。
\begin{align*} f(x)=log(1+x)\end{align*}
\begin{align*}g(x)=log(1-x)\end{align*}
\begin{align*} f(x)= x - \frac{x^2}{2} + \frac{x^3}{3}-\frac{x^4}{4}+\dots = \sum_{n=1}^{\infty}(-1)^{n-1}\frac{x^n}{n} \tag{1}\end{align*}
\begin{align*} g(x)= -x - \frac{x^2}{2} - \frac{x^3}{3}-\frac{x^4}{4}-\dots = \sum_{n=1}^{\infty}(-1)\frac{x^n}{n} \tag{2}\end{align*}
(1)から(2)を引きます。
\begin{align*} log\frac{1+x}{1-x}= 2x + \frac{2x^3}{3} + \frac{2x^5}{5}+\dots = \sum_{n=0}^{\infty}\frac{2x^{(2n+1)}}{(2n+1)} \tag{3}\end{align*}
ここで
\begin{align*} y=\frac{1+x}{1-x}\end{align*}
とおくと
\begin{align*} x=\frac{y -1}{y+1}\end{align*}
となることに留意して、(３)を書き換えると
\begin{align*} log(y)= 2\left(\frac{y-1}{y+1}\right) + \frac{2}{3}\left(\frac{y-1}{y+1}\right)^3 + \frac{2}{5}\left(\frac{y-1}{y+1}\right)^5+\dots = \sum_{n=0}^{\infty}\frac{2}{2n+1}\left(\frac{y-1}{y+1}\right)^{(2n+1)} \tag{4}\end{align*}
といかにもこのままプログラムを書けばよさそうな形に変化させる事が出来ました。しかし\(1-x\)の０近傍というのは、\(x<1\)であることを前提としており(一般二項定理の係数と同じ)\(y\)に直したところで1の近傍である必要があります。ここで平方根の時に出くわした式変形の壁に出くわすのですが、対数の場合は掛け算を足し算に変形できます、具体には
\begin{align*} log(y) = log(\alpha \cdot 2^{n} ) = n\cdot log(2)+log(\alpha)\end{align*}
として計算します。\( log(2) \)は常に変わりませんので定数として扱います。では、どうやって\(y\)を\(\alpha \cdot 2^{n} \)の形にするかですが、実数型はそもそもこの形をしているので今回はdoubleを掘り下げるのエントリを応用してプログラムを見てみましょう。
#include <iostream>
typedef union UnDouble_{
int iNum;
struct StAna{
unsigned int lfraction:20;
unsigned int exponent:11;
unsigned int sign:1;
} stAna;
} UnDouble;
typedef struct StDouble_{
UnDouble uLsb;
UnDouble uMsb;
} StDouble;
double log(double x );
#define LN2 0.6931471805599453094172321
int _tmain(int argc, _TCHAR* argv[])
{
std::cout << log(2.0) << std::endl;
return 0;
}
double log(double x)
{
int i = 1, n;
double a, r, r2;
StDouble *stAc = (StDouble*)&x;
if( stAc->uMsb.stAna.sign ) return 0;
if( x == 0 ) return 0;
n = stAc->uMsb.stAna.exponent - 1023;
r = n * LN2;
stAc->uMsb.stAna.exponent = 1023;
x = ( x - 1.0 ) / ( x + 1.0 );
a = 2.0 * x;
x = x * x;
do{
r2 = r;
r = r + a / i;
a = a * x;
i = i + 2;
}while( r != r2 );
return r;
}
ホーナー法を使ってしまっているのでちょっと見通しが悪いですが(４)式をそのままプログラムにしています。バイトオーダが変わると結果が変わってしまうと思うので注意ください。実際は標準関数では割り算のiterationを使用して\(\alpha \cdot 2^{n}\)を作っているようです。なので標準関数より高速なプログラムになっていると思います。また常用対数にしたい場合は底の変換公式 \( log_{10}x = log_{e} x / log_{e}10 \) を使用し\(log10\)を定数として扱えば簡単に実現できます。
平方根( sqrt )の数値的解法③
既に平方根の計算方法として最も有用なニュートン法を紹介していますが、その他にも興味深い方法がありますので紹介しておきます。
＜テーラー展開＞
平方根の計算方法としては難がありますが、計算機科学はテーラー展開の発見なしには語れません。\( sin(x), cos(x), exp(x), log(x)\) など、平方根以外の初等関数はそのほとんどがテーラー展開によるものですし、力学シミュレーションの漸化式を与えているのもテーラー展開です。
平方根の計算においては次の式の形に直します。(０近傍で無限回微分が可能な形)
\begin{align*} \sqrt{\alpha} = f(x) = (x+1)^{\frac{1}{2}} \tag{1}\end{align*}
(１)を0近傍で展開(マクローリン展開)してみましょう。
\begin{align*} f(x) = \sum_{n=0}^{\infty}\frac{f^{(n)}(0)}{n!}x^n = 1 + \frac{1}{2}x-\frac{1}{8}x^2+\frac{1}{16}x^3 - \frac{5}{128}x^4+\dots \end{align*}
係数の一般化をしておくと
\begin{align*} \binom{z}{m}:=\frac{z(z-1)(z-2)\dots(z-m+1)}{m!} \end{align*}
と定義して
\begin{align*} f(x) = 1 + \binom{1/2}{1}x+\binom{1/2}{2}x^2+\binom{1/2}{3}x^3 + \binom{1/2}{4}x^4+\dots \tag{2}\end{align*}
と書けます。この形はテーラー展開が発見される以前からニュートンの一般二項定理として知られた形でした。\(x\)が０に近いほど収束が早く\(x < 1\)である必要があります。ですので残念なことに例えば \(\sqrt{1+1}\)として２の平方根を取ろうとすると収束が悪いです。この制約からテーラー展開による平方根は入力値によって形を変えてやる必要があります。例えば、３、５、６、７の平方根であれば、
\begin{align*} \sqrt{3} = 2(1-\frac{1}{4})^{\frac{1}{2}} \end{align*}
\begin{align*} \sqrt{5} = 2(1+\frac{1}{4})^{\frac{1}{2}} \end{align*}
\begin{align*} \sqrt{6} = 2(1+\frac{1}{2})^{\frac{1}{2}} \end{align*}
\begin{align*} \sqrt{7} = 2(1+\frac{3}{4})^{\frac{1}{2}} \end{align*}
といったように式を直します。(実際はもっと分数部分を小さくする工夫ができます) これらの式の手直しを数値計算で対応するのは難しく、テーラー展開は平方根の計算には向いていません。では代表として\(\sqrt{2}\)のプログラムを見てみましょう。係数の性質からホーナー法を用いています。
#include "stdafx.h"
#include <iostream>
double sqrt2(void);
int _tmain(int argc, _TCHAR* argv[ ])
{
std::cout << "sqrt(2) =" << sqrt2() << std::endl;
return 0;
}
double sqrt2(void)
{
double x1=0 ,x2, a;
int n;
a = 1.0;
x2 = a ;
for(n=1; ( (a>=0)? a: -a) > 0.000001; n++){
x1 = x2;
a = a * ( (1.5 -(double)n )/(double)n);
x2 = x1 + a ;
}
std::cout << "n =" << n-1 << std::endl;
return x2;
}
このプログラムの結果は下記になります。
n =4302 sqrt(2) =1.41421
収束が悪いですが、テーラー展開の問題ではなく初期値として\(x=1\)としてしまっているからです。テーラー展開は収束半径に気を付けて使う必要があります。
＜連分数＞
連分数は通常の学習指導要領には含まれないかと思いますので学校では習はなかったかと思いますが、分数では表せないはずの無理数を分数の形で表現する事が出来ます。黄金比を表現する方法として知られています。数値計算では割り算が続くので計算量が多く好まれないのかと思われます。
さてここでは例として\(\sqrt{2}\)と\(\sqrt{5}\)を連分数の形で表して一般化してみます。
\begin{align*} \sqrt{2} = 1 + \sqrt{2} -1 \tag{1}\end{align*}
と表せます。外に出した1 は\( a^2 \le \alpha: \alpha = 2 \) である最大の\(a\)を選びます。ここで\(\sqrt{2}-1\)の項に対して乗法の公式を使って有理化っぽい事(本来は分母を有理数にする意味)をしてみます(分母に無理数を作り出すので無理化かもしれません)。
\begin{align*} \sqrt{2}-1 = \frac{\sqrt{2}+1}{\sqrt{2}+1} \times (\sqrt{2} -1) = \frac{1}{\sqrt{2} + 1} \end{align*}
となります。式①はこれにより次のように表せます。
\begin{align*} \sqrt{2} = 1 + \frac{1}{\sqrt{2} + 1} \tag{2}\end{align*}
ここで右辺の\(\sqrt{2}\) にも式②を代入し繰り返してみましょう。(再起呼び出しと同じこと)
\begin{align*} \sqrt{2} = 1 + \frac{1}{ 2 + \frac{1}{ 2 + \frac{1}{2 + \frac{1}{2 + \dots}}}} \tag{3}\end{align*}
このように連分数として表現する事が出来ました。無限に分母が続いてしまいますが、それでは計算できませんので適当なところで\(\sqrt{2}\)には１を入れます。 同様に\(\sqrt{5}\)でもやってみましょう。
\begin{align*} \sqrt{5} = 2 + \sqrt{5} -2 = 2 + \frac{3}{ 4 + \frac{3}{ 4 + \frac{3}{ 4 + \frac{3}{ 4 + \dots}}}} \tag{4}\end{align*}
こちらの場合の最終的な分母の\(\sqrt{5}\)には\( a^2 \le \alpha: \alpha = 5 \) である最大の\(a\)である２を入れます。 通常連分数では分子に１を持ってくる正則連分数が好まれるようですが、ここではこの形でOKとさせて下さい。さて、それでは代数的に一般化した形で表現してみたいと思います。
\begin{align*} \sqrt{\alpha} = a + \frac{\alpha - a^2 }{ 2a + \frac{\alpha - a^2}{ 2a + \frac{\alpha - a^2}{2a + \frac{\alpha - a^2}{2a + \frac{\alpha - a^2}{a+a}}}}} \tag{5}\end{align*}
\(\alpha\)が入力値、\(a\)は\( a^2 \le \alpha\)である最大の\(a\)を選びます。数値計算においてはこのaを選ぶのに直観というわけには行かないので、iterationで探す手間があります。それではプログラムを見てみましょう。
#include <iostream>
double sqrt(double alpha);
int _tmain(int argc, _TCHAR* argv[])
{
std::cout << "sqrt(2)=" << sqrt(2.0) << std::endl;
return 0;
}
double sqrt( double alpha ){
double numerator, x1=0, x2;
int a=0, n;
while(alpha > (a+1)*(a+1) ) a++;
numerator = alpha - (double)(a*a);
x2 = numerator / (double)(a+a);
for(n=0;x1!=x2;++n){
x1 = x2;
x2 = numerator / ( (double)2.0 * a + x1 );
}
x2 = x2 + (double)a;
std::cout<<"n="<< n << std::endl;
return(x2);
}
このプログラムの結果は下記になります。
n=21 sqrt(2)=1.41421
doubleで計算できるマックスまでiterationさせて２１回で収束しているので、二分法より早く収束しています。決して悪くない方法に思えます。
＜開平法＞
ルートを取る事を開平するといったり、ルートを開くと言ったりします。昔は学校でも必ず習う計算手法だったみたいでルートの計算といえば開平法という印象があります。幾何的な特徴からルートを求めます。基本的には筆算用の方法なので、プログラム化すると複雑化してしまいわざわざ開平法を再現するプログラムというのは見かけません(2進だとすっきりしたプログラムになりますがここでは10進でのプログラムを紹介します)。まず前提として10進の数は分解すると次のようになります。
\begin{align*} r = r_n 10^n + r_{n-1} 10^{n-1} + \dots + r_0 10^0 + r_{-1} 10^{-1}\dots = \sum_{k}^{n}r_k 10^k\tag{1}\end{align*}
\(n\)は整数、\(r_k\)は\(0 \le r_k \le 9 \)である整数です。何進数であろうとこの考え方は同じです。もう一つ前提としてルートを取りたい数\(\alpha\)が\(10^{2n} \le \alpha \le 10^ {2(n-1)}\) とすると、\( (10^{i})^2 \le \alpha \)である最大の\(i\)が\(n\)となります。
さて、では例として144の開平を見ていきます。144は\( 10^{2 \times 1} < 144 < 10^{2 \times (1+1)} \)なので開平後の最大の\(n\)は 1 となります。開平後を①式のように表してみると、
\begin{align*} \sqrt{144} = r_1 10^1 + r_0 10^0 \tag{2}\end{align*}
ルートの解として知りたいのが\(r_1\)と\(r_0\)という問題にすり替わります。②式を二乗すると、
\begin{align*} 144 = r_1^2 10^2 + 20 r_1 r_0 + r_0^2 \tag{3}\end{align*}
です。二乗というのは面積計算なので\(144＝S_k\)と置くと③式は\(S_k = S_1 + S_0 \)
\begin{align*} S_1 = r_1^2 10^2 \tag{4}\end{align*}
\begin{align*} S_0 = 20 r_1 r_0 + r_0^2 \tag{5}\end{align*}
となります。この後は段階的に\(r\)を検算して求めます。結局人間の手の終えるレベルにして検算で求めるという方法なのですが、試しに\(r_1\) に2を入れると、S1=200になって144より大きくなってしまうので\(r_1\)は１になります。とすると、\(S_0=44\)となりますので、⑤式を書き直すと、
\begin{align*} 44 = 20 r_0 + r_0^2 \tag{6} \end{align*}
です。\(r_0\)に合いそうな数を類推して当てはめます。３だと大きくなりすぎるので２です。よって答えは１２となります。開平法というと割り算みたいな計算を思い出しますが、実際やっている事はこの手続きになります。
それではプログラムを見てみましょう。
#include <iostream>
#include <deque>
#include <cmath>
double kaihei(double alpha);
int _tmain(int argc, _TCHAR* argv[ ])
{
std::cout << kaihei(2) << std::endl;
return 0;
}
double kaihei(double alpha)
{
int i, j, k, t1, t2;
std::deque<unsigned char> v_p;
std::deque<unsigned char> v_a;
double ans = 0;
i=0;
while( pow(10, 2*(i+1)) <= alpha ) i++;
j=0;
while(alpha != 0 ){
v_p.push_back( (unsigned char)(alpha / pow(10, 2*(i-j))));
alpha = alpha - v_p[j]*pow(10, 2*(i-j));
j++;
}
k=0;
while( (k+1)*(k+1) <= v_p[0]) k++;
v_a.push_back( k );
t2 = ( v_p[0] - v_a[0] * v_a[0] ) * 100;
if( j-1 > 0 ) t2 = t2 + v_p[1];
t1 = v_a[0] + v_a[0];
for(int l=1; l<10; l++){
k=0;
while( ( ( t1 * 10 ) + (k+1) ) * (k+1) <= t2 ) k++;
t2 = ( t2 - ( ( t1 * 10 ) + k ) * k) * 100;
if( j-1 > l ) t2 = t2 + v_p[l];
v_a.push_back( k );
t1 = ( ( t1 * 10 ) + k ) + k;
}
for(int m=0; m<10; m++) ans = ans + (double)( v_a[m] * pow(10, i-m) );
return (ans);
}
このプログラムの結果は下記になります。
sqrt(2)=1.41421
cmathをincludeしてしまっていますし、アロケートして配列を増やしています。さらに今までとは毛色が違いiterationの数を10と決めています。iterationするだけ下位の桁まで計算できます。基本的には標準関数には向かないと思いますが頭の体操と思ってプログラム内容を考えてみてください。
まとめ
５種類の方法を見てきたが平方の計算にはニュートン法が向いている。
平方根( sqrt )の数値的解法②
引き続き平方根の数値的解法を見ていきたいと思います。ここではニュートン法での解法をまとめておきます。ほかの方法よりも収束が早く平方根の解法としてmath.hで定義されているのがニュートン法です。
＜ニュートン法＞
ニュートン法も二分法と同じく平方根の値を方程式の解として導き出します。
\begin{align*} f(x) = x^2 - \alpha \tag{1}\end{align*}
\(\alpha\)を入力値として(1)の関数の\(x\)軸との交点が平方根の値となります。ここで任意の初期値 \(x_{0}\)を与え、\(x_{0}\)での式(1)の接線を求めます。
接線は一次関数\(g(x)\)として次のようにかけます。
\begin{align*} g(x) = f'(x_{0}) x + \chi \tag{2}\end{align*}
切片\(\chi\)を求めておきましょう。\(g(x)\)は\(x=x_{0}\)の時、(1)より\(x_{0}^{2}-\alpha\)なので(2)に代入して
\begin{align*} x_{0}^{2} -\alpha = f'(x_{0})x_{0} + \chi \end{align*}
ここで\( f'(x_{0}) = 2x_{0} \) なので、\(\chi = -(x_{0}^{2} + \alpha ) \) となります。(2)を書き直しておきます。
\begin{align*} g(x) = 2x_{0} x - x_{0}^{2} - \alpha \tag{3} \end{align*}
この接線とx軸との交点を\(x_1\)として\(g(x_1) = 0 \)から、
\begin{align*} 0 = 2x_{0} x_1 - x_{0}^{2} - \alpha \end{align*}
\begin{align*} x_1 = \frac{1}{2} (x_0 + \frac{\alpha} {x_{0}}) \tag{4}\end{align*}
(4)のように漸化式の形で\(x_1\)が求まります。同様に\(x_2, x_3, \dots \) と求めていくことで解\(x_n = \sqrt{\alpha }\) が求まります。初期値\(x_0\)は二分法と同じく\(\frac{\alpha + 1}{2}\)にしておきます。ではプログラムで確かめてみましょう。
#include <iostream>
double sqrt(double alpha);
int _tmain(int argc, _TCHAR* argv[])
{
std::cout << sqrt(2.0) << std::endl;
return 0;
}
double sqrt(double alpha){
double x0=0, x1=0.5*(1+alpha);
int n;
for( n=0; x0!=x1; ++n ){
x0 = x1;
x1 = 0.5*( x0 + alpha / x0 );
}
std::cout << n << std::endl;
return(x1);
}
このプログラムの結果は下記になります。
n = 5 sqrt(2) = 1.41421
doubleで計算できるMAXまでiterationさせていますが5回しかかかっていません。二分法が54回かかりましたから、圧倒的に収束スピードが速いことが分かります。\(\sqrt{2}\)の最終桁の値が違いますが、二分法の時と同じくstd::coutの制約です。実際はもっと多くの桁まで正しい値になっています。補足ですが二次方程式の解でニュートン法を使う場合は問題になりませんが、接線を使いますので関数の形、初期値の選び方によっては収束しないので注意してください。
平方根( sqrt )の数値的解法①
math.h で定義されている初等関数の数値計算に関してまとめておこうかと思います。(必ずしもmath.hで定義されているものと同じとは限りませんので注意ください)まずは一番初歩的な平方根( sqrt )の解法を幾つか紹介します。
<二分法>
平方根を方程式の問題として扱い\(f(x) = 0 \) となる\(x\)を挟み込んで段々と探す方法です。計算量が多いのですが、数値計算手法の中で最も基本的(力技で値を探す方法として基本的)なものとなります。それでは二分法で平方根を求めてみましょう。
\begin{align*} x = \sqrt{\alpha} \tag{1} \end{align*}
\(x\)が知りたい値です。両辺二乗して右辺を左辺に持ってきて\( x^2 - \alpha = 0 \)を得ます。この二次方程式の解が平方根の値になります。関数の形に直しておきましょう。
\begin{align*} f(x) = x^2 - \alpha \tag{2}\end{align*}
(2) は下記のような形をしています。x軸と曲線との交点を数値的に求めます。
ここからは手続きの話しになります。交点 となる\(x_{\alpha}\)に対して\(a\lt x_{\alpha}\lt b\)となる\(a, b\)を初期値として与えます。
手続き①
\(a\)と\(b\)の和の２分の１となる数値を式(2)で評価してみます。
\begin{align*} x_{1}=\frac{a+b}{2}\end{align*}
\begin{align*} f(x_{1}) = x_{1}^{2} -\alpha \end{align*}
手続き②
ここで\(f(x_{1}) \lt 0\)であるならば、解 \(x_{\alpha}\)は \(x_{1} \lt x_{\alpha} \lt b \)の領域にいるはずなので
\begin{align*} x_{2}=\frac{x_{1}+b}{2}\end{align*}
として再び式(2)で評価します。
ここで\( f(x_{1}) \gt 0\)であるならば、解 \(x_{\alpha}\)は \(a \lt x_{\alpha} \lt x_{1} \)の領域にいるはずなので
\begin{align*} x_{2}=\frac{a+x_{1}}{2}\end{align*}
として再び式(2)で評価します。
手続き③
\begin{align*} f(x_{2}) = x_{2}^{2} -\alpha \end{align*}
上記のように式(2)で評価して手続き②に戻ります。
このように、手続き②と③を繰り返す事で解がどんどん近似されていきます。許容誤差 \(\beta\) を最初に決めておき \( x_{n} - x_{n-1} \lt \beta \) となった時点でiterationをやめにします。さてここで問題になるのが初期値\(a, b\)の決め方です。平方根においては次の関係を使用します。
\begin{align*} 0 \le \sqrt{x} \le \frac{x+1}{2} \end{align*}
平方根においてよく知られた関係です。高校数学でよく出る証明問題かと思います。では実際にプログラムを見てみましょう。
#include <iostream>
double sqrt(double alpha);
int _tmain(int argc, _TCHAR* argv[])
{
std::cout << "sqrt(2) = " << sqrt(2) << std::endl;
return 0;
}
double sqrt(double alpha)
{
int n;
double a = 0, b = (alpha+1)/2;
double x1, x2, y;
x2 = b;
x1 = (a + b) / 2 ;
for(n=0;x2-x1 != 0;++n){
x2 = x1;
y = x1*x1-alpha;
if( y < 0 ){
a = x1;
x1 = ( x1 + b )/2;
}else if( y > 0){
b = x1;
x1 = ( x1 + a )/2;
}else{
break;
}
}
std::cout << "n = " << n << std::endl;
return(x1);
}
このプログラムの結果は下記になります。
n = 54 sqrt(2) = 1.41421
doubleで計算できるMAXまでiterationさせていますが54回かかっています。\(\sqrt{2}\)の最終桁の値が違いますが、これは単なるstd::coutの制約です。実際はもっと多くの桁まで正しい値になっています。
まとめ
・二分法は力技で解を探すので計算量が多いが計算機での解法の基礎が凝縮されている。
doubleを掘り下げる
doubleの中身を確認し整理しておこうと思います。浮動小数点演算標準としてIEEE754で定義されていますが、少し難解な印象があります。まず前提としてすべての実数\(x\)は
\begin{align*} x=(-1)^{(2-\delta)} \cdot \alpha \cdot 2 ^n \tag{1}\end{align*}
として表せます。ここで\(\delta\)は0か1(符号)、nは整数(指数)、\(\alpha\)は\(1\le\alpha\lt 2\)である実数(端数)です。\(1\le\alpha\lt 2\)は
\begin{align*} S_n = 1 + \sum_{k=1}^n \frac{1}{2^k} \tag{2}\end{align*}
が\( \lim\limits_{n->\infty}S_n=2 \)ですので、(2)の\( \frac{1}{2^k} \)の和の組み合わせで表現されます。
doubleはこの符号(sign) と指数(exponent)と端数(fraction)を決められたbitに配置する型です。IEEE754ではこれらを下記のように配置する事を決めています。
| 63 | 62 | 61 | 60 | 59 | 58 | 57 | 56 | 55 | 54 | 53 | 52 | 51 | 50 | 49 | 48 | ~ | 3 | 2 | 1 | 0 |
| sign | exponent | fraction | ||||||||||||||||||
1. sign
符号は1をマイナス、0をプラスと決めています。
2. exponent
指数部は2進数11桁なので2048個の数値を持ちます、値の持ち方が独特で実際にはunsignedで値を保持しながら、オフセット(バイアス)で-1023する事になっています。それにより -1023～1024 までを表現する事が出来ます。例えば、
| 62 | 61 | 60 | 59 | 58 | 57 | 56 | 55 | 54 | 53 | 52 |
| 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
はunsingedでは\(2^{10}\)で1024ですが、オフセットで-1023するので実際は１となります。
3. fraction
端数部は52bitの中に小数点以下の情報を保持します。オフセットとして＋１を持ちます。例えば
| 51 | 50 | 49 | 48 | 47 | 46 | 45 | ~ | 2 | 1 | 0 |
| \(2^{-1}\) | \(2^{-2}\) | \(2^{-3}\) | \(2^{-4}\) | \(2^{-5}\) | \(2^{-6}\) | \(2^{-7}\) | ~ | \(2^{-50}\) | \(2^{-51}\) | \(2^{-52}\) |
| 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
であれば、\(1 + 2^{-1} + 2^{-3} = 1.625 \)となります。
ではテストプログラムで確かめてみましょう。
#include <iostream>
#include <bitset>
typedef union UnDouble_{
int iNum;
struct StAna{
unsigned int lfraction:20;
unsigned int exponent:11;
unsigned int sign:1;
} stAna;
} UnDouble;
typedef struct StDouble_{
UnDouble uLsb;
UnDouble uMsb;
} StDouble;
int _tmain(int argc, _TCHAR* argv[])
{
double foo = -11.0;
StDouble *stAc = (StDouble*)&foo;
std::bitset<32> bit[2];
bit[0] = stAc->uMsb.iNum;
bit[1] = stAc->uLsb.iNum;
std::cout << bit[0] << bit[1] << std::endl;
std::cout << "sign=" << stAc->uMsb.stAna.sign << std::endl;
std::cout << "exponent=" << stAc->uMsb.stAna.exponent-1023 << std::endl;
stAc->uMsb.stAna.sign = 0;
stAc->uMsb.stAna.exponent = 1023;
std::cout << "fraction=" << foo << std::endl;
return 0;
}
構造体変数のポインタでキャストする事で直接signとexponentを取り出しています。fractionは52桁もありますのでdoubleじゃないと表現できませんsignとexponentの部分を意図的に消し込むことで抽出しています。なおこのプログラムはバイトオーダに依存するかと思いますので注意ください。このプログラムの入力はfoo = -11.0 です。実行結果は
1100000000100110000000000000000000000000000000000000000000000000 sign=1 exponent=3 fraction=1.375
となります。この結果を(1)の式に入れてみましょう。
\begin{align*} x=(-1)^{(2-1)} \cdot 1.375 \cdot 2 ^3 \end{align*}
ちゃんと -11.0 になっています。
まとめ
・ doubleは実数を(１)の式のsign, exponent, fractionに分けて、それぞれに決まったbitを割り当てたもの。
ビットフィールドを掘り下げる②
ビットフィールドを掘り下げる①では主にビットフィールドで宣言した構造体のサイズに注目しました。今回は符号修飾 signed, unsigned, 符号なしに注目したいと思います。
struct StBit{
unsigned char foo: 4;
unsigned char bar: 4;
};
通常、上記のようにビットフィールドのメンバ変数はunsigned で宣言します。ではこれを、
struct StBit{
unsigned char foo: 4;
signed char bar: 4;
};
としたらどうなるでしょうか？結論から言うと、メンバ変数barは4bitの中でちゃんとsingedとして動作します。符号bitありという事ですからsinged char bar: 4; で表現できる数は -8 ～ 7です。一方unsigned char foo :4; で表現できる数は 0 ～ 15 になります。
ではテストプログラムを書いて確認してみましょう。
#include <stdio.h>
struct StBit{
unsigned char foo: 4;
signed char bar: 4;
};
int _tmain(int argc, _TCHAR* argv[])
{
StBit stBit;
stBit.foo = 15;
stBit.bar = 15;
printf("sum = %d\n", stBit.foo + stBit.bar );
return 0;
}
このテストプログラムのsumの値は何になりますか？ 15は0b1111ですからfooにもbarにも0b1111が代入されます。しかしbarはsignedですから最上位bitが立っている場合、マイナスをつけた2の補数で表されます。この場合0b1111 の2の補数なので0b0001(ビット反転して1足す) なのでシステム的には-1という認識をします。よって15＋(-1)となりコンソールには sum = 14 と表示されます。なお、signed char bar : 1; とした場合bar = 1は-1と認識されます(1の2の補数は1)基本的にビットフィールドでsignedにする事なんてないかと思うのでunsignedをつけるのを標準化しているわけです。そうでなければ演算をした場合に思わぬバグが出る可能性があります。最後に
struct StBit{
char foo: 4;
char bar: 4;
};
これだと、fooとbarの符号はどうなりますか？これは何度か説明した通り、システム依存になります。Visual Studioではcharが標準でsignedになるので、やっぱりunsignedをつけておかないと演算で失敗する可能性があります。
まとめ
・ビットフィールドにはunsignedをつける。
ビットフィールドを掘り下げる①
ビットフィールドを指定した場合、一般には下のような書き方をします。
struct StBit0 {
unsigned btFirst: 1;
unsigned btSecond: 1;
unsigned btNull: 5;
unsigned btLast: 1;
};
今回注目したいのはビットフィールドで宣言した構造体変数のサイズについてです。上の例では構造体変数一個につき何バイトの領域が確保されるのでしょうか？ビットフィールドとしては 1+1+5+1 で8bitしか使っていませんから1byteだとうれしいですが、この場合はint型の大きさ(4byte)分取られます。おそらくこれは処理系に依存しません。注目すべきはunsigned です。これは unsigned int 表記の省略形でした。という事はunsigned intで領域が確保されているという事です。さてでは次の場合はどうなるでしょうか？
struct StBit0 {
unsigned char btFirst: 1;
unsigned char btSecond: 1;
unsigned char btNull: 5;
unsigned char btLast: 1;
};
unsigned charとなっていますから、1byteの領域確保のような気がします。もしくは良く勉強している人だと4byteと答えるかもしれません。しかし答えは処理系依存です。書籍によってはビットフィールドで確保する領域というのは、どの型で指定したところでint型の大きさの倍数になるという書き方をしているものがありますが、実際は処理系に依存するようです。なおVisualStudioではこの宣言は1byteの領域確保になりました。 デュアルターゲット開発をする場合の落とし穴になることがありますから認識しておく必要のある項目です。
#include <stdio.h>
struct StBit1 {
unsigned char a: 8;
};
struct StBit2 {
unsigned short a: 8;
};
struct StBit3 {
unsigned long a: 8;
};
int _tmain(int argc, _TCHAR* argv[])
{
/*sizeof で調べる*/
printf("StBit1 = %d\n",sizeof(StBit1));
printf("StBit2 = %d\n",sizeof(StBit2));
printf("StBit3 = %d\n",sizeof(StBit3));
/*素直に差分を求める*/
StBit1 stBit1;
StBit2 stBit2;
StBit3 stBit3;
printf("real StBit1 = %d\n", ( (int)( (&stBit1)+1)-(int)(&stBit1) ) );
printf("real StBit2 = %d\n", ( (int)( (&stBit2)+1)-(int)(&stBit2) ) );
printf("real StBit3 = %d\n", ( (int)( (&stBit3)+1)-(int)(&stBit3) ) );
return 0;
}
このプログラムの結果は
StBit1 = 1 StBit2 = 2 StBit3 = 4 real StBit1 = 1 real StBit2 = 2 real StBit3 = 4
となります。
まとめ
・ビットフィールドにビット演算をするようなプログラム、たとえばstBit2 &= 0x0100;のようなプログラムがあると環境によって結果が変わってしまう事があります。こういった処理はしないようにしましょう。
