kashiの日記

はじめに

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

#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公式を用いた数値積分のための関数群を追加。

• dd (double-double)用のfrexpが、入力が2のベキ乗よりも僅かに小さい値のときに正しく動作していなかったバグの修正。
• interval.hppで、内部型としてfrexpが誤差を含む可能性があるような型を用いていたとき (例えばdd)、logがごく稀に正しい値を返していなかったバグの修正。
• ddの(精度保証付きでない)sinとcosで、引き戻し処理が雑だったのを修正。

ddのfrexpのバグ

    double x, y;
    int i;
    y = frexp(x, &i);

    kv::dd x, y;
    int i;
    y = frexp(x, &i);

    friend dd frexp(const dd& x, int* m) {
        double z1, z2;
        z1 = std::frexp(x.a1, m);
        z2 = std::ldexp(x.a2, -(*m));
        return dd(z1, z2);
    }

    friend dd frexp(const dd& x, int* m) {
        double z1, z2;
        z1 = std::frexp(x.a1, m);
        z2 = std::ldexp(x.a2, -(*m));
        if ((z1 == 0.5 && z2 < 0) || (z1 == -0.5 && z2 > 0)) {
            z1 *= 2;
            z2 *= 2;
            (*m)--;
        }
        return dd(z1, z2);
    }

x: 0.999999999999999944488848768742172978818416595458984375
y: 0.4999999999999999722444243843710864894092082977294921875
i: 1

x: 0.999999999999999944488848768742172978818416595458984375
y: 0.999999999999999944488848768742172978818416595458984375
i: 0

interval<dd>のlogのバグ

• y = 0.5+2^-1075
• i = 1

[0,0]

[0,9.881312916824930883531375857364427447301196052286495288511713650013510145404
17503730599672723271984759593129390891435461853313420711879592797549592021563756
25260142638062280905569163433569796420737743727211399746144610001277481830712996
87746249467945463392302800634307707961482524771311823420533171133735363740791206
21249863890543182984910658610913088802254960259419999083863978818160833126649049
51429573802945356031871047722310026960705298694403875805362142149834066644536895
06671441664863872184765786916736120212023012339619506156684554636658495809965049
46155275185449574931216955640746893939906729403594535543517025132110239826300978
22029020757254763345019116747794671979873296198823284114052741805584855350891304
5817507736501283943653106689453125e-324]

ddのsin,cosの引き戻し問題

おわりに

• 辺に特異性を持つ2変数関数の二重積分
• 一点に特異性を持つ2変数関数の二重積分