kashiの日記

• dd-interval.pdf
• inflate.pdf

#include <kv/interval.hpp>
#include <kv/dd.hpp>
#include <kv/rdd.hpp>

int main(){
        kv::interval<kv::dd> x;
        std::cout.precision(34);

        x = std::numeric_limits<kv::dd>::max();

        std::cout << x << "\n";
        std::cout << sqrt(x) << "\n";
}

[1.340780792994259561100831764080293e+154,1.340780792994259672743259815885529e+154]

[1.340780792994259672743259815885478e+154,1.340780792994259672743259815885529e+154]

はじめに

精度保証せずに普通に計算

#include <iostream>
#include <cmath>

int main()
{
	double s, tmp;
	long long int i;
	std::cout.precision(17);

	s = 0;
	i = 1;
	while (true) {
		tmp = s + 1.0 / i;
		if (floor(tmp) != floor(s)) {
			std::cout << i-1 << ":" << s << "\n";
			std::cout << i << ":" << tmp << "\n";
		}
		s = tmp;
		i++;
	}
}

区間演算してみる

#include <kv/interval.hpp>
#include <kv/rdouble.hpp>
#include <kv/dd.hpp>
#include <kv/rdd.hpp>

typedef kv::interval<double> itv;
// typedef kv::interval< kv::dd > itv;

int main()
{
	itv s, tmp;
	long long int i;
	std::cout.precision(17);

	s = 0;
	i = 1;
	while (true) {
		tmp = s + 1.0 / itv(i);
		if (floor(tmp.lower()) != floor(tmp.upper())) {
			std::cout << "warning\n";
			std::cout << i-1 << ":" << s << "\n";
			std::cout << i << ":" << tmp << "\n";
			break;
		} else if (floor(tmp.lower()) != floor(s.lower())) {
			std::cout << i-1 << ":" << s << "\n";
			std::cout << i << ":" << tmp << "\n";
		}
		s = tmp;
		i++;
	}
}

double-double区間演算で頑張る

もうちょっと数学っぽく

#include <iostream>
#include <kv/interval.hpp>
#ifdef DD
#include <kv/dd.hpp>
#include <kv/rdd.hpp>
#endif
#include <kv/gamma.hpp>

#ifdef DD
typedef kv::interval<kv::dd> itv;
#define MAX 66
#else
typedef kv::interval<double> itv;
#define MAX 32
#endif

int main() {
	#ifdef DD
	std::cout.precision(34);
	#else
	std::cout.precision(17);
	#endif

	itv gamma, x, y;
	int i;

	gamma = - kv::digamma(itv(1.));

	for (i=2; i<=MAX; i++) {
		std::cout << i << "\n";

		x = i;
		y = floor(mid(exp(x - gamma) + 0.5));
		std::cout << y.lower() << "\n";
		std::cout <<  kv::digamma(y) + gamma << "\n";;
		std::cout <<  kv::digamma(y + 1) + gamma << "\n";;
	}
}

• test-rounding.ccで、非正規化数を強制的に0にする"flush to zero"モードになっていないかどうかをチェックするようにした。
• iccでコンパイルできるように修正。
• optimize.hppで、「端にない」区間で傾きが単調ならばそこに最小値は存在しない というチェックを入れる改良を行った。
• defint-newtoncotes.hpp, double-newtoncotes.hppで、最適な分割数の 推定がうまくできていなかったバグを修正。
• defint.hpp, defint-singular.hppで、急減少関数などの特定の被積分関数 に対してstep sizeがうまく計算できていなかった問題を修正。
• psa.hppの中のinv関数を、係数のオーバーフローが発生しづらいように改良。
• gamma.hppの中のdigamma関数の計算アルゴリズムをより高精度なものに変更。
• その他細かい修正。

• 区間演算を行うinterval.hppにfloorとceilを追加。
• 区間演算と領域の再帰的分割による最適化を行うoptimize.hppが場合によって正しく終了しなかったバグの修正。
• 高階微分を計算するhighderiv.hppが定数関数を渡されたとき動作しなかったバグの修正。
• 昨年11月のNVR2020の発表で作成した、複合Newton-Cotes公式を用いた数値積分のための関数群を追加。