2015/12/23(水)kv-0.4.28

精度保証
kvライブラリを0.4.28にアップデートしました。

今回はライブラリ本体の変更はほとんどなく、highderiv.hpp(高階微分の計算)、test-rounding.cc(丸めモードの変更が正常に動いているか簡易チェック)を追加したくらいです。webのドキュメントは、psaやaffineにいくらか加筆しています。明日、応用数理セミナーの講演でこのライブラリの話をするのでその前に区切りとして全ての変更を放出しておきたかったという理由もあります。

test-rounding.ccの追加は、前回日記で触れたcygwinのsqrt問題があったので、とりあえず軽くチェックしようと思ったためです。いろいろな環境でテストしてみたところ、概ね以下のような感じです。
  • Linux, Macは64bit OSで64bitコンパイルしていれば問題なし
  • x86の32bitモードはFPUが使われてしまうためIEEE754に従っておらず、お勧め出来ない。
  • Windowsでは、cygwin, msys2で最適化をしないと、sqrtに丸めの向きが変わらない問題が発生。VC++の64bitなら問題なし。
sqrtはVC++の32bitモードでも問題が発生することが分かっており、VC++由来のライブラリを使っている環境ではsqrtの方向付き丸めをエミュレーションで何とかすることも考えてみます。
コメント（0件）

2015/12/16(水)cygwinのsqrtの丸めモード変更

精度保証
普段は自分は使っていないのですが、学生の一人がcygwinを使っていて、奇妙な現象を見つけたのでメモ代わりに記録しておきます。

問題は、cygwinのsqrtで丸めの向きが最適化をかけないときに変わらないというものです。環境は64bitのwindows7で、cygwinも64bitの最新のものです。gccは4.9.3でした。
#include <stdio.h>
#include <math.h>
#include <fenv.h>

int main()
{
        volatile double x, y, z;

        x = 2.;
        fesetround(FE_DOWNWARD);
        y = sqrt(x);
        fesetround(FE_UPWARD);
        z = sqrt(x);
        fesetround(FE_TONEAREST);

        if (y == z) {
                printf("sqrt error\n");
        }
}
このようなプログラムで、
cc cygwin-sqrt.c
のように最適化オプションをつけないでコンパイルしたとき、丸めの向きが上向きでも下向きでも結果が同じになってしまいます。
cc -O3 cygwin-sqrt.c
のように最適化をかければ問題ありません。いろいろ試すと、-O0ではアウト、-O1, -O2, -O3はセーフでした。

volatileを使わない場合最適化をかけることで結果がおかしくなるのはさほど珍しくありませんが、このように最適化をかけない状態で結果が異常になるのはとても珍しいです。

cc -Sでアセンブリソースを表示させてみると、最適化無しでは単に外部のsqrt関数を
        call    sqrt
と呼んでいるのに対して、-O3だと 
        sqrtsd  %xmm0, %xmm1
のようにSSE2で直接計算するコードが生成されているので、どうやら外部(libm?)のsqrtに問題がありそうです。

以前にVisual C++で似たようなことがあったので、それとも関係あるのでしょうか。どなたか詳しい方に教えて欲しいものです。

これで丸めの向きが変わらないのはIEEE754的にまずいので、そんな腐った環境は無視するというのも一つの考え方でしょうが、こういう事例が多くあるようならsqrtに関しては丸めモード変更に期待しないで何とかすることを検討するべきかも。kvライブラリでいえば-DKV_NOHWROUNDを付ければ問題なくなるので、それをデフォルトにするとか。
コメント（1件）

2015/11/27(金)区間演算の実装について(2)

精度保証
区間演算の実装について(1)の記事の続きです。前の記事でC++による区間演算の簡単な実装を示しましたが、それの問題点をいくつか挙げていきます。

まず、入出力の問題です。区間型の内部で使われているdoubleは2進数で表現されているので、区間型に数値を代入するとき10進数→2進数変換が、区間型の数値を表示するとき2進数→10進数変換が行われます。精度保証付き数値計算を行うためには、この変換が適切に行われる必要があります。例えば、
#include "interval.hpp"

int main()
{
	interval x;
	std::cout.precision(17);

	x = 0.1;
	std::cout << x << "\n";
}
のようなプログラムを実行すると、
[0.10000000000000001,0.10000000000000001]
となってしまい、なぜかxは0.1を含んでいません。これは精度保証付き数値計算のビギナーが陥りやすい罠で、ソーステキストに書いた「0.1」はコンパイラが2進数に変換しますが、0.1は2進数だと無限小数になるのでどこかで打ち切る必要があり、正確に0.1になりません。実際、2進数53桁で0.1に最も近い数は、
1.1001100110011001100110011001100110011001100110011010 x 2-4
になります。これは10進数だと、
0.1000000000000000055511151231257827021181583404541015625
という数です。頑張って工夫して、例えば
        fesetround(FE_DOWNWARD);
        x.lower() = 0.1;
        fesetround(FE_UPWARD);
        x.upper() = 0.1;
と書いたとしても、この0.1を2進数に変換するのはコンパイル時なので、どちらの0.1も上で示した近似値になってしまい、意図通りの0.1を含む区間にはなりません。(あまり大きくない)整数は2進数で正しく表現できるので、
        x = 1;
        x /= 10;
のように書けば一応0.1を含む区間をxに入れることは出来ますが、毎回このように書くのは大変です。根本的に解決するには、ソーステキストの「0.1」をコンパイラに変換させないこと、すなわちソーステキストのレベルでは文字列で「"0.1"」のように保持し、それの2進数への変換はこちらでやってあげることになります。

出力でも同じ問題があって、例えば
#include "interval.hpp"

int main()
{
	interval x, y, z;

	x = 1;
	y = 10;
	z = x / y;
	std::cout << z << "\n";
	std::cout.precision(17);
	std::cout << z << "\n";
}
のようなプログラムを実行すると
[0.1,0.1]
[0.099999999999999992,0.10000000000000001]
のようになり、前者はまるで幅の無い区間のように見えてしまいます。これは、表示をcout(というか標準ライブラリ)に任せているせいで、本当は10進文字列への変換を自力で行い、下端は下向き丸めで、上端は上向き丸めで行わなければなりません。

2つ目に、最適化の問題があります。以下に、
g++ (Ubuntu 4.8.4-2ubuntu1~14.04) 4.8.4
で実行した結果を示します。他のシステムやコンパイラでは少し違うかも知れません。
まず、
#include "interval.hpp"

int main()
{
	interval x, y;
	std::cout.precision(17);

	x = 1.;
	y = 10.;
	std::cout << x / y << "\n";
}
のようなプログラムを、オプション無しでコンパイルして実行すると、
[0.099999999999999992,0.10000000000000001]
となりますが、これを-O3をつけて最適化して実行すると、
[0.10000000000000001,0.10000000000000001]
のように誤った結果が得られてしまいます。驚いたことに、最終行を2回実行するだけの
#include "interval.hpp"

int main()
{
	interval x, y;
	std::cout.precision(17);

	x = 1.;
	y = 10.;
	std::cout << x / y << "\n";
	std::cout << x / y << "\n";
}
だと、最適化を行っても行わなくても
[0.099999999999999992,0.10000000000000001]
[0.099999999999999992,0.10000000000000001]
のように正しい結果が得られました。これは推測なのですが、これはコンパイラがプログラムを解析して、実行結果が分かりきっている部分を予めコンパイル時に計算しておく、いわゆる定数最適化のせいだと考えられます。コンパイル時には丸めの方向は変更されませんので、結果がおかしくなってしまいます。プログラムが単純だと除算の結果が常に1/10だと気づくが、ある程度以上複雑になるとそれに気づかないのでしょうか。最適化の影響を抑えるのはかなり難しい問題なのですが、この場合はその計算に関連する変数にvolatile修飾子を付けることがよく行われます。volatile修飾子は、その変数がいつどんなタイミングで変更されるか分からない特殊な変数であることをコンパイラに伝えるもので、コンパイラはその変数に関係する部分の最適化を抑制します。

第3の問題は、加減乗除と平方根以外の、例えばexpやsinなどの関数の問題です。sinやcosなどの非単調な関数の区間演算は元々難しいのですが、単調増加関数であるexpの区間演算を行うとき、sqrtと同様に
        friend interval exp(const interval& x) {
                interval r;

                fesetround(FE_DOWNWARD);
                r.inf = exp(x.inf);
                fesetround(FE_UPWARD);
                r.sup = exp(x.sup);
                fesetround(FE_TONEAREST);

                return r;
        }
と書いて良いのか? 結論から言うと、このように書いてはいけません。IEEE754の方向付き丸めは、加減乗除と平方根のみに有効で、それ以外の演算に関してはどんな結果になるか保証されていません。そもそも加減乗除と平方根以外の演算では計算結果の精度の保証が何もありませんので、丸めの方向の指定どころの話ではありません。実際調べた結果がここにあります。すなわち、加減乗除と平方根以外の演算に関して区間演算を行うには、その関数を自力で実装する必要があります。

以上により、区間演算をきちんと実装するにはかなり手間がかかることが分かります。以下に、これら3つの問題点を解決した実装例を挙げておきます。
interval-simple-0.4.28.tar.gz
([2015/12/23]より、kvライブラリの一部に含めました。)
これは、上で挙げた3つ問題点を全て解決したものです。
  • 文字列とdoubleの相互変換に関しては、dtostring, stringtodという関数で丸め方向の指定が出来る相互変換を実装しました。これにより、例えば x = "0.1"; と書けばxは0.1を含む区間になりますし、表示は必ず外側に丸めて行われます。
  • 最適化対策は、例えば上向き丸めの加算は
    	static double add_up(const double& x, const double& y) {
		volatile double r, x1 = x, y1 = y;
		fesetround(FE_UPWARD);
		r = x1 + y1;
		fesetround(FE_TONEAREST);
		return r;
	}
    のように加算に関係する全ての変数をvolatileとしています。
  • 数学関数は全て自力で実装しています。実装方法についてはここを参照して下さい。
kvライブラリの中のinterval.hppからtemplateとboost依存を取り除いたもの、でもあります。単一ファイルで使いやすいと思いますので、大掛かりなkvライブラリを使うほどのことでもない場合にご利用下さい。
コメント（0件）

2015/11/27(金)区間演算の実装について(1)

精度保証
区間演算は、数値を[下端, 上端]のペアで表して一種の集合値演算として数値計算を行う手法で、精度保証付き数値計算と基本となるものです。下端を計算するときは下向き丸め、上端を計算をするときは上向き丸めで計算することによって、丸め誤差の影響を把握します(真の値を区間内に包含できます)。両端点を倍精度浮動小数点数で持つ場合は、IEEE754に備わっている方向付き丸めの機能を利用できます。

区間演算は原理は簡単ですが、実際に実装となると案外ハマリポイントが多く、難しいものです。簡単な区間演算の実装例を通して、少し解説してみようと思います。

使いやすい区間演算を実装するには、区間型を定義してその区間型を組み込み型と同等に扱える、いわゆる演算子多重定義(operator overloading)のできるプログラミング言語が望ましいです。自分はC++を愛用していますので、ここでもC++で実装します。両端点の型はdouble、方向付き丸めはC99のfenv.hとfesetroundを使います。

最も単純で最低限の機能のみの実装を示します (interval.hpp)。
#ifndef INTERVAL_HPP
#define INTERVAL_HPP

#include <iostream>
#include <cmath>
#include <stdexcept>
#include <fenv.h>

class interval {
	double inf;
	double sup;

	public:

	interval() {
		inf = 0.;
		sup = 0.;
	}

	interval(const double& x) {
		inf = x;
		sup = x;
	}

	interval(const double& x, const double& y) {
		inf = x;
		sup = y;
	}

	friend interval operator+(const interval& x, const interval& y) {
		interval r;

		fesetround(FE_DOWNWARD);
		r.inf = x.inf + y.inf;
		fesetround(FE_UPWARD);
		r.sup = x.sup + y.sup;
		fesetround(FE_TONEAREST);

		return r;
	}

	friend interval operator+(const interval& x, const double& y) {
		interval r;

		fesetround(FE_DOWNWARD);
		r.inf = x.inf + y;
		fesetround(FE_UPWARD);
		r.sup = x.sup + y;
		fesetround(FE_TONEAREST);
		return r;
	}

	friend interval operator+(const double& x, const interval& y) {
		interval r;

		fesetround(FE_DOWNWARD);
		r.inf = x + y.inf;
		fesetround(FE_UPWARD);
		r.sup = x + y.sup;
		fesetround(FE_TONEAREST);

		return r;
	}

	friend interval& operator+=(interval& x, const interval& y) {
		x = x + y;
		return x;
	}

	friend interval& operator+=(interval& x, const double& y) {

		fesetround(FE_DOWNWARD);
		x.inf = x.inf + y;
		fesetround(FE_UPWARD);
		x.sup = x.sup + y;
		fesetround(FE_TONEAREST);

		return x;
	}

	friend interval operator-(const interval& x, const interval& y) {
		interval r;

		fesetround(FE_DOWNWARD);
		r.inf = x.inf - y.sup;
		fesetround(FE_UPWARD);
		r.sup = x.sup - y.inf;
		fesetround(FE_TONEAREST);

		return r;
	}

	friend interval operator-(const interval& x, const double& y) {
		interval r;

		fesetround(FE_DOWNWARD);
		r.inf = x.inf - y;
		fesetround(FE_UPWARD);
		r.sup = x.sup - y;
		fesetround(FE_TONEAREST);

		return r;
	}

	friend interval operator-(const double& x, const interval& y) {
		interval r;

		fesetround(FE_DOWNWARD);
		r.inf = x - y.sup;
		fesetround(FE_UPWARD);
		r.sup = x - y.inf;
		fesetround(FE_TONEAREST);

		return r;
	}

	friend interval& operator-=(interval& x, const interval& y) {
		x = x - y;
		return x;
	}

	friend interval& operator-=(interval& x, const double& y) {

		fesetround(FE_DOWNWARD);
		x.inf = x.inf - y;
		fesetround(FE_UPWARD);
		x.sup = x.sup - y;
		fesetround(FE_TONEAREST);

		return x;
	}

	friend interval operator-(const interval& x) {
		interval r;

		r.sup = - x.inf;
		r.inf = - x.sup;

		return r;
	}

	friend interval operator*(const interval& x, const interval& y) {
		interval r;
		double tmp;

		if (x.inf >= 0.) {
			if (y.inf >= 0.) {
				fesetround(FE_DOWNWARD);
				r.inf = x.inf * y.inf;
				fesetround(FE_UPWARD);
				r.sup = x.sup * y.sup;
			} else if (y.sup <= 0.) {
				fesetround(FE_DOWNWARD);
				r.inf = x.sup * y.inf;
				fesetround(FE_UPWARD);
				r.sup = x.inf * y.sup;
			} else {
				fesetround(FE_DOWNWARD);
				r.inf = x.sup * y.inf;
				fesetround(FE_UPWARD);
				r.sup = x.sup * y.sup;
			}
		} else if (x.sup <= 0.) {
			if (y.inf >= 0.) {
				fesetround(FE_DOWNWARD);
				r.inf = x.inf * y.sup;
				fesetround(FE_UPWARD);
				r.sup = x.sup * y.inf;
			} else if (y.sup <= 0.) {
				fesetround(FE_DOWNWARD);
				r.inf = x.sup * y.sup;
				fesetround(FE_UPWARD);
				r.sup = x.inf * y.inf;
			} else {
				fesetround(FE_DOWNWARD);
				r.inf = x.inf * y.sup;
				fesetround(FE_UPWARD);
				r.sup = x.inf * y.inf;
			}
		} else {
			if (y.inf >= 0.) {
				fesetround(FE_DOWNWARD);
				r.inf = x.inf * y.sup;
				fesetround(FE_UPWARD);
				r.sup = x.sup * y.sup;
			} else if (y.sup <= 0.) {
				fesetround(FE_DOWNWARD);
				r.inf = x.sup * y.inf;
				fesetround(FE_UPWARD);
				r.sup = x.inf * y.inf;
			} else {
				fesetround(FE_DOWNWARD);
				r.inf = x.inf * y.sup;
				tmp = x.sup * y.inf;
				if (tmp < r.inf) r.inf = tmp;
				fesetround(FE_UPWARD);
				r.sup = x.inf * y.inf;
				tmp = x.sup * y.sup;
				if (tmp > r.sup) r.sup = tmp;
			}
		}
		fesetround(FE_TONEAREST);

		return r;
	}

	friend interval operator*(const interval& x, const double& y) {
		interval r;

		if (y >= 0.) {
			fesetround(FE_DOWNWARD);
			r.inf = x.inf * y;
			fesetround(FE_UPWARD);
			r.sup = x.sup * y;
		} else {
			fesetround(FE_DOWNWARD);
			r.inf = x.sup * y;
			fesetround(FE_UPWARD);
			r.sup = x.inf * y;
		}
		fesetround(FE_TONEAREST);

		return r;
	}

	friend interval operator*(const double& x, const interval& y) {
		interval r;

		if (x >= 0.) {
			fesetround(FE_DOWNWARD);
			r.inf = x * y.inf;
			fesetround(FE_UPWARD);
			r.sup = x * y.sup;
		} else {
			fesetround(FE_DOWNWARD);
			r.inf = x * y.sup;
			fesetround(FE_UPWARD);
			r.sup = x * y.inf;
		}
		fesetround(FE_TONEAREST);

		return r;
	}

	friend interval& operator*=(interval& x, const interval& y) {
		x = x * y;
		return x;
	}

	friend interval& operator*=(interval& x, const double& y) {
		x = x * y;
		return x;
	}

	friend interval operator/(const interval& x, const interval& y) {
		interval r;

		if (y.inf > 0.) {
			if (x.inf >= 0.) {
				fesetround(FE_DOWNWARD);
				r.inf = x.inf / y.sup;
				fesetround(FE_UPWARD);
				r.sup = x.sup /  y.inf;
			} else if (x.sup <= 0.) {
				fesetround(FE_DOWNWARD);
				r.inf = x.inf / y.inf;
				fesetround(FE_UPWARD);
				r.sup = x.sup / y.sup;
			} else {
				fesetround(FE_DOWNWARD);
				r.inf = x.inf / y.inf;
				fesetround(FE_UPWARD);
				r.sup = x.sup / y.inf;
			}
		} else if (y.sup < 0.) {
			if (x.inf >= 0.) {
				fesetround(FE_DOWNWARD);
				r.inf = x.sup / y.sup;
				fesetround(FE_UPWARD);
				r.sup = x.inf / y.inf;
			} else if (x.sup <= 0.) {
				fesetround(FE_DOWNWARD);
				r.inf = x.sup / y.inf;
				fesetround(FE_UPWARD);
				r.sup = x.inf / y.sup;
			} else {
				fesetround(FE_DOWNWARD);
				r.inf = x.sup / y.sup;
				fesetround(FE_UPWARD);
				r.sup = x.inf / y.sup;
			}
		} else {
			fesetround(FE_TONEAREST);
			throw std::domain_error("interval: division by 0");
		}
		fesetround(FE_TONEAREST);

		return r;
	}

	friend interval operator/(const interval& x, const double& y) {
		interval r;

		if (y > 0.) {
			fesetround(FE_DOWNWARD);
			r.inf = x.inf / y;
			fesetround(FE_UPWARD);
			r.sup = x.sup / y;
		} else if (y < 0.) {
			fesetround(FE_DOWNWARD);
			r.inf = x.sup / y;
			fesetround(FE_UPWARD);
			r.sup = x.inf / y;
		} else {
			fesetround(FE_TONEAREST);
			throw std::domain_error("interval: division by 0");
		}
		fesetround(FE_TONEAREST);

		return r;
	}

	friend interval operator/(const double& x, const interval& y) {
		interval r;

		if (y.inf > 0. || y.sup < 0.) {
			if (x >= 0.) {
				fesetround(FE_DOWNWARD);
				r.inf = x / y.sup;
				fesetround(FE_UPWARD);
				r.sup = x / y.inf;
			} else {
				fesetround(FE_DOWNWARD);
				r.inf = x / y.inf;
				fesetround(FE_UPWARD);
				r.sup = x / y.sup;
			}
		} else {
			fesetround(FE_TONEAREST);
			throw std::domain_error("interval: division by 0");
		}
		fesetround(FE_TONEAREST);

		return r;
	}

	friend interval& operator/=(interval& x, const interval& y) {
		x = x / y;
		return x;
	}

	friend interval& operator/=(interval& x, const double& y) {
		x = x / y;
		return x;
	}

	friend std::ostream& operator<<(std::ostream& s, const interval& x) {
		s << '[';
		s << x.inf;
		s << ',';
		s << x.sup;
		s << ']';
		return s;
	}

	friend interval sqrt(const interval& x) {
		interval r;

		if (x.inf < 0.) {
			throw std::domain_error("interval: sqrt of negative value");
		}

		fesetround(FE_DOWNWARD);
		r.inf = sqrt(x.inf);
		fesetround(FE_UPWARD);
		r.sup = sqrt(x.sup);
		fesetround(FE_TONEAREST);

		return r;
	}

	const double& lower() const {
		return inf;
	}
	const double& upper() const {
		return sup;
	}
	double& lower() {
		return inf;
	}
	double& upper() {
		return sup;
	}
};

#endif // INTERVAL_HPP
使える演算は加減乗除と平方根で、coutでの表示にも対応してます。lowerとupperは下限と上限へのアクセサ。乗除算は手抜きせずに一応ちゃんと符号による場合分けを行っています。C++の初心者でも入門書を読みながら実装できるレベルでしょう。

次にこのライブラリの使用例を示します (test-interval.cc)。
#include "interval.hpp"

int main()
{
	interval x;
	interval y = 1.;
	interval z(1.);

	x = 1.;
	y = 10.;
	z = x / y;

	std::cout << z << "\n";
	std::cout.precision(17);
	std::cout << z << "\n";

	// copy
	x = interval(1., 2.);
	y = interval(3., 4.);

	// basic four operations
	std::cout << x + y << "\n";
	std::cout << x - y << "\n";
	std::cout << x * y << "\n";
	std::cout << x / y << "\n";

	// operation with constant
	std::cout << x + 1 << "\n";
	std::cout << x + 1. << "\n";

	// compound assignment operator
	z += x;
	z += 1.;

	z = interval(3., 4.);
	// access to endpoints
	std::cout << z.lower() << "\n";
	std::cout << z.upper() << "\n";
	z.lower() = 3.5;
	std::cout << z << "\n";
}
実行すると、
[0.1,0.1]
[0.099999999999999992,0.10000000000000001]
[4,6]
[-3,-1]
[3,8]
[0.25,0.66666666666666674]
[2,3]
[2,3]
3
4
[3.5,4]
のようになります。この2つのファイルを一応アップロードしておきます。

interval-supersimple.tar.gz

これを使って、大きく誤差が出る計算をしてみます。2次方程式ax2+bx+c=0でa=1, b=1015, c=1014として、2つの解のうちの大きい方を、
  • 普通に解の公式を使う
  • 解の公式の分子を有理化する
の2通りで解いてみます。

まず普通にdoubleで。
#include <iostream>
#include <cmath>

int main()
{
	double a, b, c, x, y;

	std::cout.precision(17);

	a = 1.;
	b = 1e15;
	c = 1e14;

	x = (-b + sqrt(b * b - 4. * a * c)) / (2. * a);
	y = 2 * c / (-b - sqrt(b * b - 4. * a * c));

	std::cout << x << "\n";
	std::cout << y << "\n";
}
のようなプログラムを実行すると、
-0.125
-0.10000000000000002
のようになります。この計算は後者がほぼ正しく、前者は大きな誤差が入っています。

これを区間演算で計算してみます。プログラムは、
#include <iostream>
#include <cmath>
#include "interval.hpp"

int main()
{
        interval a, b, c, x, y;

        std::cout.precision(17);

        a = 1.;
        b = 1e15;
        c = 1e14;

        x = (-b + sqrt(b * b - 4. * a * c)) / (2. * a);
        y = 2 * c / (-b - sqrt(b * b - 4. * a * c));

        std::cout << x << "\n";
        std::cout << y << "\n";
}
のような感じ。演算子多重定義のおかげで最小限の変更で区間演算が行えます。計算結果は、
[-0.1875,-0.0625]
[-0.10000000000000003,-0.099999999999999992]
となり、前者の区間幅が極端に広いことで、計算結果が怪しいと気づくことが出来ます。

さて、学習用としてはこれでいいのですが、実際に使おうとするとこのレベルの実装ではいろいろと問題点があります。少し記事が長くなったので、問題点は次の記事で。
コメント（0件）

2015/11/26(木)kv-0.4.27

精度保証
kvライブラリを0.4.27にアップデートしました。

今回の変更点は、
  • 公比が区間な等比級数の和の公式を修正し、それに伴いhypergeom.hppを修正。
  • optimize.hppで、求まった区間のうち隣接しているものを結合する機能を追加。
  • allsol.hppで、TRIMのアルゴリズムを少し変更。
  • constants<affine<T>>::pi()などが機能していなかった問題を修正。
  • webデモの問題記述でpi, e, ln2を使えるようにした。
  • webデモにaffineとintervalによる関数の値域評価を追加。
などです。久しぶりにwebデモに少しだけ機能追加しましたが、全く中身を覚えていなくて大変でした。
コメント（0件）
OK キャンセル 確認 その他