1.4.5.2 test-interval.cc を読む

   1 #include "interval.hpp"
   2 
   3 int main()
   4 {
   5         interval x;
   6         interval y = 1.;
   7         interval z(1.);
   8 
   9         x = 1.;
  10         y = 10.;
  11         z = x / y;
  12 
  13         std::cout << z << "\n";
  14         std::cout.precision(17);
  15         std::cout << z << "\n";
  16 
  17         // copy
  18         x = interval(1., 2.);
  19         y = interval(3., 4.);
  20         // assign
  21         y.assign(3., 4.);
  22 
  23         // basic four operations
  24         std::cout << x + y << "\n";
  25         std::cout << x - y << "\n";
  26         std::cout << x * y << "\n";
  27         std::cout << x / y << "\n";
  28 
  29         // operation with constant
  30         std::cout << x + 1 << "\n";
  31         std::cout << x + 1. << "\n";
  32 
  33         // compound assignment operator
  34         z += x;
  35         z += 1.;
  36 
  37         z = interval(3., 4.);
  38         // access to endpoints
  39         std::cout << z.lower() << "\n";
  40         std::cout << z.upper() << "\n";
  41         z.lower() = 3.5;
  42         std::cout << z << "\n";
  43 
  44         // static functions
  45         // whole and hull
  46         std::cout << interval::whole() << "\n";
  47         std::cout << interval::hull(1., 2.) << "\n";
  48         std::cout << interval::hull(2., 1.) << "\n";
  49         std::cout << interval::hull(1., z) << "\n";
  50         std::cout << interval::hull(z, 1.) << "\n";
  51         std::cout << interval::hull(x, y) << "\n";
  52 
  53         // friend functions
  54         std::cout << width(z) << "\n";
  55         std::cout << rad(z) << "\n";
  56         std::cout << median(z) << "\n";
  57         std::cout << mid(z) << "\n";
  58         std::cout << norm(z) << "\n";
  59         std::cout << mag(z) << "\n";
  60         std::cout << std::boolalpha;
  61         std::cout << in(3.9, z) << "\n";
  62         std::cout << in(4.1, z) << "\n";
  63         std::cout << subset(x, y) << "\n";
  64         std::cout << subset(y, z) << "\n";
  65         std::cout << proper_subset(y, z) << "\n";
  66         std::cout << overlap(x, y) << "\n";
  67         std::cout << overlap(y, z) << "\n";
  68         std::cout << intersect(y, z) << "\n";
  69         std::cout << abs(interval(2., 3.)) << "\n";
  70         std::cout << abs(interval(-2., 3.)) << "\n";
  71         std::cout << abs(interval(-2., -1.)) << "\n";
  72         std::cout << mig(interval(-2., 1.)) << "\n";
  73         std::cout << mig(interval(-2., -1.)) << "\n";
  74         std::cout << mig(interval(2., 3.)) << "\n";
  75 
  76         std::cout << max(interval(2., 4.), interval(3., 5.)) << "\n";
  77         std::cout << min(interval(2., 4.), interval(3., 5.)) << "\n";
  78 
  79         // test for midrad
  80         std::cout << "test for midrad\n";
  81         double m, r;
  82         z = interval(1., 1 + std::numeric_limits<double>::epsilon() * 3);
  83         std::cout << mid(z) << "\n";
  84         std::cout << rad(z) << "\n";
  85         midrad(z, m, r);
  86         std::cout << m << "\n";
  87         std::cout << r << "\n";
  88 
  89         // comparison operators
  90         std::cout << (x < y) << "\n";
  91         std::cout << (x < 1.) << "\n";
  92         std::cout << (1. < x) << "\n";
  93         std::cout << (interval(1.) == 1.) << "\n";
  94         std::cout << (interval(1., 2.) != 3.) << "\n";
  95 
  96         // division_part1, division_part2
  97         // calculate X / (Y \ 0). the result may be divided into two part.
  98         bool parted;
  99         std::cout << division_part1(interval(1., 2.), interval(-3., 4.), parted) << "\n";
 100         // if the result has division_part2, parted is set to true.
 101         if (parted) std::cout << division_part2(interval(1., 2.), interval(-3., 4.)) << "\n";
 102 
 103         // integer power
 104         std::cout << pow(interval(2., 3.), 2) << "\n";
 105         std::cout << pow(interval(2., 3.), -2) << "\n";
 106         std::cout << pow(interval(-2., 3.), 2) << "\n";
 107         std::cout << pow(interval(-2., 3.), 3) << "\n";
 108         // general power
 109         std::cout << pow(interval(2., 3.), interval(2., 3)) << "\n";
 110         // mathematical functions
 111         std::cout << sqrt(interval(2.5, 3.5)) << "\n";
 112         std::cout << exp(interval(2.5, 3.5)) << "\n";
 113         std::cout << exp(interval(-std::numeric_limits<double>::infinity(), 0.)) << "\n";
 114         std::cout << exp(interval(0., std::numeric_limits<double>::infinity())) << "\n";
 115         std::cout << expm1(interval(-0.25, 0.25)) << "\n";
 116         std::cout << log(interval(0.75, 1.25)) << "\n";
 117         std::cout << log(interval(1., std::numeric_limits<double>::infinity())) << "\n";
 118         std::cout << log(interval(0., 1.)) << "\n";
 119         std::cout << log1p(interval(-0.25, 0.25)) << "\n";
 120         std::cout << sin(interval(-0.25, 0.25)) << "\n";
 121         std::cout << cos(interval(-0.25, 0.25)) << "\n";
 122         std::cout << tan(interval(-0.25, 0.25)) << "\n";
 123         std::cout << atan(interval(-0.25, 0.25)) << "\n";
 124         std::cout << asin(interval(-0.25, 0.25)) << "\n";
 125         std::cout << acos(interval(-0.25, 0.25)) << "\n";
 126         std::cout << atan2(interval(1.), interval(1.)) << "\n";
 127         std::cout << sinh(interval(-0.25, 0.25)) << "\n";
 128         std::cout << cosh(interval(-0.25, 0.25)) << "\n";
 129         std::cout << tanh(interval(-0.25, 0.25)) << "\n";
 130         std::cout << asinh(interval(-0.25, 0.25)) << "\n";
 131         std::cout << acosh(interval(1.5, 2.)) << "\n";
 132         std::cout << atanh(interval(-0.25, 0.25)) << "\n";
 133 
 134         // string is converted to interval which includes the number represented
 135         // by the string
 136 
 137         // initialize by string. 
 138         x = "0.1";
 139         std::cout << x << "\n";
 140         // operation with string
 141         std::cout << "0.1" * x << "\n";
 142         x += "0.1";
 143         // comparison with string
 144         std::cout << ("0.20001" > x) << "\n";
 145 
 146         // numeric constants
 147         std::cout << interval::pi() << "\n";
 148         std::cout << interval::e() << "\n";
 149         std::cout << interval::ln2() << "\n";
 150 }

以下解読してみる (かなりの急ぎ仕事なので、間違えている可能性がある)。

• 5〜7行目 interval 名前; で区間の下限、 上限が double 型である区間の変数が定義できる。 初期化も = で double 型データを代入したり、 () でコンストラクトしたりできる。
• 6, 7, 9, 10 行目 1 や 10 のように double で正確に表せる数値ならば、 double の代入で区間 (いわゆる点区間 $a=[a,a]=\{a\}$になる)　を設定できる。
• 13〜15行目 C++ の普通の出力の仕方 (std::cout << ) で出力できる。
• 18〜19行目 interval($a$, $b$) で区間の下限 $a$、 上限 $b$ を指定して区間 $[a,b]$ が作れる (区間の左端、右端とか言いたくなるが、 区間演算業界では、下限、上限というのが普通)。
  以下 (色々計算しているが一貫して) x は $[1,2]$, y は $[3,4]$ である。
• 24〜27行目 四則演算は分かるでしょう。
• 30〜31行目 あたり前だけど、定数と演算できる (ように作られている)。
• 34〜35行目 これは C 言語ファミリーならできて自然。
• 39〜41行目 .lower(), .upper() として、区間の下限、 上限が得られる。これらは代入もできる。
  以下 (色々計算しているが一貫して) z は $[3.5,4]$ である。
• 46行目 whole() で全体 $[-\infty,\infty]$ を表す。 何に使うんだろう。
• 47〜51行目 hull(a,b) は、 a とb の区間包と訳すのかな、 2つのものを含む最小の区間を求める関数。 $[1,2]$ と $[3,4]$ の区間包は $[1,4]$ である。
• 54行目 width() は区間の幅。$[3.5,4]$ の幅は $0.5$.
• 55行目 rad() は区間の半径。$[3.5,4]$ の半径は $0.25$.
• 56〜57行目 median() も mid() も区間の中心。 $[3.5,4]$ の中心は $3.75$.
• 58〜59行目 norm() も mag() も区間に含まれる数の絶対 値の最大値。$[3.5,4]$ の場合は $4$.
• 60行目 ？
• 61〜62行目 in($a$,$X$) で $a\in X$ であるかどうか。
• 63〜64行目 subset($X$,$Y$) で $X\subset Y$ であるかどうか。
• 65行目 proper_subset($X$,$Y$) で $X\Subset Y$ ($X=[a,b]$ が $Y=[c,d]$ の内部 $(c,d)$ に含まれるか、すなわち $c<a\land b<d$ である) かどうか。
• 66〜67行目 overlap($X$,$Y$) で $X\cap Y\ne\emptyset$ であるかどうか。
• 64行目 intersect($X$,$Y$) で $X\cap Y$ を返す。 (intersect だと動詞のようで紛らわしい！ intersection() にして欲しかった。)
• 69〜71行目 abs($X$) で $\left\{\vert x\vert\relmiddle\vert x\in X\right\}$. $X=[2,3]$ ならば $[2,3]$, $X=[-2,3]$ ならば $[0,3]$, $X=[-2,-1]$ ならば $[1,2]$.
• 69〜71行目 mig($X$) で $\min\left\{\vert x\vert\relmiddle\vert x\in X\right\}$ かなあ？ $X=[-2,1]$ ならば 0, $X=[-2,-1]$ ならば $1$, $X=[2,3]$ ならば $2$.
• 76〜77行目 max($X$, $Y$) で $\left\{\max\{x,y\}\relmiddle\vert x\in X, y\in Y\right\}$, min($X$, $Y$) で $\left\{\min\{x,y\}\relmiddle\vert x\in X, y\in Y\right\}$ かなあ。
• 79〜87行目 z は $[1,1+3\eps_M]$. ($\eps_M$ は machine epsilon というもので、IEEE 754 の binary 64 の場 合は、
  $$\eps_M=2^{-52}=2.220446049250313080847263336181640625\times 10^{-16}.$$
  である。) z の中心は数学的には $1+1.5\eps$ であるが、 コンピューターで mid(z) と計算すると、 $1$ と比べて $\eps_M$ が小さすぎて、$1$ になってしまう。 z の半径 rad(z) は $1.5\eps_M$ であるが、 [mid(z)-rad(z),mid(z)+rad(z)] が z を含まないことになる。 midrad(z,m,r) で m と r を計算すると、 [m-r,m+r] が z を含むようになる。
• 89〜94行目は(大小)比較演算子。 $X<Y$ は ( $\forall x\in X$) ( $\forall y\in Y$) $x<y$ という意味である。 $[1,2]<[3,4]$ は真、 $[1,2]<[1,1]$ は偽、$1<[1,2]$ も偽、$[1,1]=1$ は真 ($1$ は区間と比較するときは $[1,1]$ という区間にされるからだろう)、 $[1,2]\ne3$ は真である。
• 96〜101行目は 0 を含む区間での割り算。 $X$, $Y$が区間で $0\in Y$とするとき、普通 $X/Y$ は定義されないとするが (区間演算の常識？)、
  $$\left\{x/y\relmiddle\vert x\in X, y\in Y\setminus\{0\}\right\} ... ...right\} \cup \left\{x/y\relmiddle\vert x\in X, y\in Y\cap(0,\infty)\right\}$$
  は意味がある。 $X=[1,2]$, $Y=[-3,4]$とするとき、

  $$\left\{x/y\relmiddle\vert x\in X, y\in Y\setminus\{0\}\right\} =\left\{x/y\relmiddle\vert 1\le x\le 2, -3\le y<0\right\}\cup \left\{x/y\relmiddle\vert 1\le x\le 2, 0<y\le 4\right\}$$

  $$=(-\infty,-1/3]\cup[1/4,\infty).$$

  
  

• 103〜107行目は整数冪乗、つまり$X$は区間, $n\in\mathbb{Z}$とするとき $X^n= \left\{x^n\relmiddle\vert x\in X\right\}$. $X$ が負の数を含んでいても大丈夫。

  $$[2,3]^2=\left\{x^2\relmiddle\vert 2\le x\le 3\right\}=[2^2,3^2]=[4,9],$$

  $$[2,3]^{-2}=\left\{x^{-2}\relmiddle\vert 2\le x\le 3\right\}=[3^{-2},2^{-2}] =[1/9,1/4]=[1.11111\cdots,0.25],$$

  $$[-2,3]^2=\left\{x^{-2}\relmiddle\vert -2\le x\le 3\right\}=[0^2,3^2]=[4,9],$$

  $$[-2,3]^3=\left\{x^{-2}\relmiddle\vert -2\le x\le 3\right\}=[(-2)^3,3^3]=[-8,27].$$

  
  

• 108〜109行目は一般の冪乗、$X$ が正の数しか含まない区間の場合、 区間 $Y$ に対して $X^Y=\left\{x^y\relmiddle\vert x\in X, y\in Y\right\}$ は 意味を持つ。
  $$[2,3]^{[2,3]}=[2^2,3^3]=[4,27].$$

• 110〜132行目 数学関数。 ほとんどは説明が不要と思われる (作る人は大変だろうが)。
  $$\exp [-\infty,0]=[0,1],\quad \exp[0,\infty]=[1,\infty],$$
  $$\log [1,\infty]=[0,\infty],\quad \log[0,1]=[-\infty,0]$$
  は少しぎょっとするかもしれないが。
• 137〜144行目 文字列 (string) は区間と演算するときは、その文字列の表す 数を含む区間となる、ということである。 文字列でない 0.1 は、$0.1$ という実数に近い double 型の値 であるが、"0.1" は $0.1$ を含む機械区間 (下限と上限が 浮動小数点数である区間) である。 この機能がないと interval(1)/interval(10) のような、 ちょっと見にくいコーディングが必要になる。
• 146〜149行目 よく使う $\pi$, $e$, $\log 2$ を含む区間。 何でもないようだが、これがないと不便である。ありがたい。

数学で区間 $[a,b]$ と言うと、 $a,b\in\mathbb{R}$, $a<b$ であることが仮定されているが、 区間演算では、$a=b$ の場合も含める。 特にコンピューター上の区間演算では、 $a$, $b$ は (広い意味の)浮動小数点数とする、いわゆる機械区間が使われる。 浮動小数点数の全体 $\mathbb{F}$ には、 普通は、実数のうちで有限桁の2進数で表現できるものしか含まれない。 例えば $0.1\not\in\mathbb{F}$. 一方で、 $\mathbb{F}$ に $+\infty$ と $-\infty$ が含まれる (これは柏木先生の趣味というべきか？)。