2011年7月7日木曜日

[C/C++] マシンイプシロン

マシンイプシロンを求めるプログラム例

#include <stdio.h>

int main(void)
{
  double deps=1.0//Machine epsilon, type: double
  double dtmp;      //Working variable, type: double
  float feps=1.0;   //Machine epsilon, type: float
  float ftmp;       //Working variable, type: float

  //Compute machine epsilon of type double
  //During dtmp = deps + 1.0, divide dtmp by 2

  for(dtmp = deps + 1.0; dtmp > 1; deps /= 2.0, dtmp = deps + 1.0);
  printf("Machine Epsilon of type double is %-16g\n", 2.0*deps);
  printf("Unit roundoff of type double is %-16g\n", deps);

  //Compute machine epsilon of type float
  //During ftmp + 1 > 1, divide ftmp by 2

  for(ftmp = feps + 1.0; ftmp > 1; feps /= 2.0, ftmp = feps + 1.0);
  printf("Machine Epsilon of type float is %-16g\n", 2.0*feps);
  printf("Unit roundoff of type float is %-16g\n", feps);

  return 0;

}

上記のプログラムをコンパイル,実行すると以下のようになります.
Machine Epsilon of type double is 2.22045e-16     
Unit roundoff of type double is 1.11022e-16     
Machine Epsilon of type float is 1.19209e-07     
Unit roundoff of type float is 5.96046e-08     

2011年7月2日土曜日

[C/C++] 倍精度型浮動小数点数フォーマットを出力

キーボードから入力された数の倍精度型浮動小数点フォーマットを出力するプログラム例を以下に示します.

#include <stdio.h>
#include <string.h>

int main(void)
{
  double a;
  char c[sizeof(a)];  // The data size of char type is 1 byte.
  int i, j, k;

  printf("Please input a number ---> ");
  scanf("%lf", &a);

  memcpy(c, &a, sizeof(a)); // Copy the contents of a to array c.

  printf("sEEEEEEE EEEEdddd dddddddd... \n");

  for( i = sizeof(a) - 1 ; i >= 0; i--)
  {
    for ( j = sizeof(c[i]) * 8 - 1; j >= 0; j--)
    {
      k = c[i] & (1 << j); //Logically multiply j by one by one.
      printf("%d", k ? 1 : 0); //If k = 0, 1 is displayed. 0 is displayed if k = 0.
    }
    printf(" ");
  }
  printf("\n");

  return 0;

}

実行結果は以下のようになります.
Please input a number ---> 777
sEEEEEEE EEEEdddd dddddddd... 
01000000 10001000 01001000 00000000 00000000 00000000 00000000 00000000 

上記のプログラムでは,1バイトごとにビット演算を行うため,char型の配列を用意しています(char型のデータサイズは1バイト,int型のデータサイズは4バイト).
文字配列 c の各要素 c[i] のサイズは1バイトなので,64ビットの倍精度型変数を格納するには c の配列サイズを8とする必要があります.
ここでは,sizeof関数を使って配列サイズをsizeof(a)=8としています.

memcpy関数は,メモリに格納されている内容をコピーするための関数で,例えばmemcpy(a, b, n)とすると,バッファbのnバイトの内容をバッファaにコピーします.
aとbには領域の先頭アドレスが入ります.
プログラムでは,データサイズが8バイトであるdouble型変数aの内容を文字配列cにコピーするために,aの先頭アドレス&a,配列cの先頭アドレスcおよびデータサイズsizeof(a)=8をmemcpy関数の実引数としています.

2011年3月22日火曜日

[C/C++] void型とvoid型ポインタ

void型とは,型を持たない型,存在しない変数や関数,あるいは何もしない変数や関数を意味し,汎用型と呼ばれることもあります.例えば,void型関数といえば,値を返さない関数のことです.また,引数が全てvoid型の関数とは,引数がない関数のことです.

void型ポインタ(void*型)とは,どんな型へのポインタ変数にでも代入できるポインタ型です.どの型にも一致するポインタとも言えます.従って,ポインタ演算において型の不一致が起こらないので,void*型で定義された関数や変数をキャストすれば,色々な型のポインタ変数として利用することができます.ただし,void*型ポインタを格納できるだけなので,例えば

void *vp
double x = 3.14 ;
vp = &x ;

printf("The content pointed to by vp is %f. \n, *vp) ;

というような間接参照はできません.間接参照する場合は,以下のように目的の型へ変換する必要があります.

#include <stdio.h>

int main(void)
{
  double x = 3.14, *fp ;
  void *vp ;

  vp = &x ; //The void * type can assign a pointer variable to any type.
  fp = vp ;

  printf("The address of pointer variable fp is %p. \n", fp) ;
  printf("The address of pointer variable vp is %p. \n", vp) ;

  printf("The address of pointer variable fp is %f. \n", *fp) ;

  //Cast is necessary for indirect reference.
  printf("The address of pointer variable vp is %f. \n", *(double *)vp) ;

  return 0 ;

}





2011年3月21日月曜日

[C/C++] キャスト

配列でデータ領域を確保すると,それが不要になったときも,プログラムが終了するまで,その配列がメモリを占有してしまうので,配列に占有された部分が無駄になります.
そこで,記憶領域を効率よく利用するために,必要なときに記憶領域を確保し,不要になったらそれを開放する機能を動的メモリ確保(Dynamic memory allocation)としいます.
そして,記憶領域を確保することをアロケートといい,それを実現するために使用される関数をアロケート関数といいます.また,動的に確保した領域を解放するための関数がfree関数です.アロケート関数を利用するには"stdlib.h"が必要です.

malloc関数の返す型は,voidへのポインタ型,つまり,void*となっています.
この意味を理解するためにはvoid型ポインタの考えを知っておく必要があります.
また,malloc関数を使うにはキャストという考え方も知る必要があります.

キャストとは,変数や式といったオペラントの型を一時的に変更することです.
キャストを行うにはキャスト演算子を使用します.
これは,データの型変換を行う演算子で,明示的に型変換を行いたい場合に使用します.
例えば,int型の変数 i, j をdouble型へ変換したい場合には,

int i, j ;
(double)i ;
(double) i+j ;

とします.
この(double0 の部分がキャスト演算子です.
逆にdouble型の変数aをint型へ変換したいときには(int)aとします.

ただし,int型へのポインタ変数iをdouble型へのポインタ変数に変換しようとして

#include <stdio.h>

int main(void)
{
  int *i ;
  (double *)i = 1.0 / 2.0 ;

  return 0 ;
}

とするとエラーになります.
なお,文法上は以下のようにして変数を用意してポインタ変数のキャストが可能です.

#include <stdio.h>

int main(void)
{
  int *i ;
  double a ;

  i = (int *)&a ; //

  return 0 ;
}

しかし,このようなポインタのキャストはプログラムが読みにくくなるので,止めた方が良いとされています.

2011年3月14日月曜日

[C/C++] プログラムの実行時間の計測

プログラムの実行時間を測るにはclock関数を使います.このclock関数を使うには time.hが必要です.この関数を使ってプログラムの実行時間を計測するには,計測したい部分の始めとと終わりにclock関数を記述してその差を取ります.

clock関数は返り値としてclock_t型を返します.多くの場合は,clock.t型はlong型と同じです.以下に実行時間を測るためのプログラム例を示します.

#include <stdio.h>
#include <time.h>
#include <math.h>

#define N 10000000

int main(void)
{
  clock_t t ;
  double a ;
  long i ;

  t = clock() ;                     //Start measurement

  for( i = 1 ; i <= N ; i++)
  {
    sin(30.0) * cos(20.0) * tan(15.0) + sin(30.0) + cos(20.0) + tan(15.0) ;
  }

  t = clock() - t ;                 //Acquire execution time
  a = (double)t / CLOCKS_PER_SEC ;  //Convert execution time to seconds

  printf("Execution time is %f seconds (%ld μ s)\n", a, t) ;

  return 0 ;
}

上記のプログラムで"(double)"は計算結果をdouble型にするための命令です.また,clock関数で得られた値を秒単位に変換するには値をCLOCKS_PER_SECで割ります.上記のプログラムをコンパイル,実行すると以下のようになります.
Execution time is 0.034697 seconds (34697 μ s)

2011年3月1日火曜日

[C/C++] ホイン法

ホイン法は常微分方程式
        y' = f(x, y)        (1)
        y(a) = y0           (2)
を満たすy(x)の近似値を求めるための方法です.ここで,a, y0は与えられた定数,f(x, y)は与えられた関数です.今,上記の2式を閉区間[a, b]上で与え,区間[a, b]をn等分して,各点xiを
        h = (b-a)/n,    x0 = a,     xi = a + ih    (i = 0, 1, 2, ..., n)
とします.なお,nを分割数,hを刻み幅と呼びます.

式(1)より
  (3)
となります.ここで,2変数のテイラー展開より
       
(4)
が成り立つので,yi = f(xi) とし,式(4)において,α=1, β=f(xi, yi)として式(3)を適用すれば

(5)
となります.ただし,O(h^n)はランダウの記号で,今の場合は最低次数がnであるhの多項式です.

一方,yi = f(xi)とするとテイラー展開より
となるので,これと式(5)より
となります.この式において,hは十分に小さいとしてO(h^3)を無視して,Yi = yi,つまり,xi上のy(x)の近似値をYiとすれば,ホイン法
が得られます.

ホイン法のアルゴリズムは以下のようになります.

以下に,ホイン法で常微分方程式
        y' = x + y,        0<= x <= 1,        y(0) = 1
を解くプログラム例を示します.

#include <stdio.h>

double f(double x, double y) ;

int main(void)
{
  double x = 0.0, y = 1.0, b = 1.0, h, k1, k2 ;
  int i, n ;

  printf("Input the number of partitions. ---> ") ;
  scanf("%d", &n) ;
  h = (b - x) / n ;

  for (i = 0 ; i < n ; i++)
  {
    k1 = f(x, y) ; k2 = f(x + h, y + h * k1) ;
    y = y + h / 2.0 * (k1 + k2) ;
    x += h ;
    printf("x = %f \t y = %f \n", x, y) ;
  }
  return 0 ;
}

double f(double x, double y)
{
  return(x + y) ;
}

Input the number of partitions. ---> 10
x = 0.100000  y = 1.110000 
x = 0.200000  y = 1.242050 
x = 0.300000  y = 1.398465 
x = 0.400000  y = 1.581804 
x = 0.500000  y = 1.794894 
x = 0.600000  y = 2.040857 
x = 0.700000  y = 2.323147 
x = 0.800000  y = 2.645578 
x = 0.900000  y = 3.012364 
x = 1.000000  y = 3.428162 

2011年2月28日月曜日

[C/C++] ニュートン法

ニュートン法は方程式
f(x) = 0        (1)
の解を求めるための方法です.原理を以下に説明します.

関数f(x)は3回微分可能だとします.方程式
f(x)= 0
の解をαとして,αの近くにある点x0の周りでテイラー展開すると
  (2)
となります.これよりΔx = α-x0が十分に小さければ
なので,
f(x0) + f'(x0)Δx = 0
は,式(2)の近似式と考えることができます.式(3)より,Δxの近似としてΔx0をΔx0=-(f(x0) / f'(x0))と選ぶと,x1 = x0 + Δx0 ≒ x0 + Δx = αは,x0よりアルファに対する良い近似となっているはずです.これを繰り返して得られる次の反復式が,式(1)に対するニュートン法です.

ニュートン法はΔ = α - x0が十分に小さい,つまり,初期値x0がαの近くにあるということが前提となっているので,初期値の選び方によっては収束しない可能性があります.したがって,あらかじめ反復回数の上限を設定しておかなければ,プログラムが終了しない可能性があります.

以上より,ニュートン法のアルゴリズムは以下のようになります.
  1. 初期値x0,小さい数 ε > 0,最大反復回数kmaxを与える.
  2. n = 0, 1, 2, …kmaxに対してxn+1 = xn - (f(xn) / f'(xn))を計算する.ここで,f'(x)はxの導関数.
  3. 補正量 - (f(xn) / f'(xn))の絶対値がεより小さくなれば xn + 1を解とし,そうでなければ解は見つからなかったこととする.また,反復回数がkmaxに達した場合も解が見つからなかったこととする.
このアルゴリズムをC/C++で書くと以下のようになります.

  x,ε,kmaxの入力
  k <- 0
  do
  {
          d <- -f(x)/f'(x)
          x <- x+d
          k <- k+1
  }while (|d| > ε and k < kmax)
  if
          k = kmax
          解が見つからない
  else
          xを出力
  endif

上記のアルゴリズムに基づいて,x - cosx = 0を解くためのプログラム例を以下に示します.

#include <stdio.h>
#include <math.h>

#define EPS pow(10.0, -8) //Epsilon setting
#define KMAX 10           //Maximum iteration function

double f(double x) ;      //Calculation of f(x)
double df(double x) ;     //Calculation of f'(x)

int main(void)
{
  int k = 0 ;
  double x, d ;
  printf("Input the initial value. \n") ;
  scanf("%lf", &x) ;
  do
  {
    d = -f(x)/df(x) ;
    x = x + d ;
    k ++ ;
  } while(fabs(d) > EPS && k < KMAX) ;

  if ( k == KMAX )
  {
    printf("An answer was not found \n") ;
  }
  else
  {
    printf("The answer is x = %f. \n", x) ;
  }

  return 0 ;
}

double f(double x)
{
  return(x - cos(x)) ;
}

double df(double x)
{
  return(1.0 + sin(x)) ;

}

上記のプログラムをコンパイル,実行すると以下のようになります(プログラムを3回実行した結果).
Input the initial value. 
1
The answer is x = 0.739085. 
Input the initial value. 
3
The answer is x = 0.739085. 
Input the initial value. 
5
An answer was not found 

2011年2月21日月曜日

[C/C++] エラーのチェック

fscanf関数は,ファイルのエラーや終わりを検出するとEOFをその戻り値として返します.しかし,EOFが返されたとしても,それだけではファイルの終端なのか,エラーが発生したのか分かりません.これを調べる関数としてfeof関数とferror関数があります.
  • feof関数:ファイルの終わりを検出すると0でない数を返し,そうでなければ0を返します.
  • ferror関数:ファイルにエラーがあると0でない数を返し,エラーがなければ0を返します.
データを読み込む際にエラーの発生を調べたい場合は,以下のようにします.

#include <stdlib.h>    //It is necessary to use "foes" and "ferror".
...
  if(foes(fp)) //If "foes" is not 0, it reaches the end of the file.
  {
    printf("It is the end of the file. \n") ;
  }
  if (ferror(fp))
  {
    printf("A reading error occurred. \n") ;

  }

データの書き込み時にも同様のチェックは必要となりますが,これにはfeof関数やferror関数を使う必要はありません.例えば,fprintf関数を使った場合は,fprintf関数が負の値を出力したかどうかを調べればよく,fwrite関数を使った場合はfwrite関数が出力した値とプログラム中で指定した要素数が等しいか否かを調べれば良いからです.

2011年2月20日日曜日

[C/C++] fopen関数のモード

fopen関数のモードとしては,以下の表のようなものがあります.

モードファイル機能ファイルが
ないとき
"r"テキスト読み取りエラー
"w"テキスト書き込み新規作成
"a"テキスト追加書き込み新規作成
"rb"バイナリ読み取りエラー
"wb"バイナリ書き込み新規作成
"ab"バイナリ追加書き込み新規作成
"r+"テキスト更新(読み取り及び書き込み)エラー
"w+"テキスト更新(読み取り及び書き込み)新規作成
"a+"テキスト更新(追加書き込み)新規作成
"r+b" or "rb+"バイナリ更新(読み取り及び書き込み)エラー
"w+b" or "wb+"バイナリ更新(読み取り及び書き込み)新規作成
"a+b" or "ab+"バイナリ更新(追加書き込み)新規作成

ファイル"test1.dat","test2.dat","test3.dat"の内容が,以下のように同じものだとします.
9876 ABCDEFGHIJK

その上で,以下のプログラムを実行します.

#include <stdio.h>

int main(void)
{
  FILE *frPlus, *fwPlus, *faPlus ;

  frPlus = fopen("test1.dat", "r+") ;
  fwPlus = fopen("test2.dat", "w+") ;
  faPlus = fopen("test3.dat", "a+") ;

  fprintf(frPlus, "%d", 1234) ;
  fprintf(fwPlus, "%d", 1234) ;
  fprintf(faPlus, "%d", 1234) ;

  fclose(frPlus) ; fclose(fwPlus) ; fclose(faPlus) ;

  return 0 ;

}

すると,以下のような結果が得られます.

test1.dat("r+"モード)
1234 ABCDEFGHIJK

test2.dat("w+"モード)
1234

test3.dat("a+"モード)
9876 ABCDEFGHIJK1234

この結果より,

  • "r+"モード ファイルを読み込んだ上で,データを先頭から上書きされる(読み込みも可能)
  • "w+"モード データを先頭から上書きし,読み込みも可能
  • "a+"モード データを最後尾に追記し,読み込みも可能

ということがわかります.