全ての記事をそのまま表示している。
とても重くなってしまっているかもしれない。
Ubuntu 16.04 で emacs 25.3 から octave 4.0 を使おうとするとハングする (org-babel から実行しても実行が終わらない)。
octave を最新版 (4.4) にしたら直った。最新版は ppa にないので手元でビルドした。
octave の最新版をビルドするときにライブラリが足りないという警告がいくつか出るのでパッケージをインストールして対処したが、 以下の警告は最後まで残った。
configure: WARNING: SUNDIALS IDA library not configured with IDAKLU, ode15i and ode15s will not support the sparse Jacobian feature今は ode15 コマンドを使わないのでほっておいたが、後で必要になるかもしれない。
タイトルの通りで、 std::is_arithmetic_v を含むコードを clang++-5.0 でコンパイルしたらエラーが出た。
clang はバージョン 5.0 でC++17 に対応しているはずなのにおかしい。
原因は C++ の標準ライブラリが C++17 に対応していないことにあった。 clang++ はデフォルトだと標準ライブラリとして libstdc++ を用いるが、 このバージョンが自分の環境だと 5.4 で C++17 に対応していなかった。
というわけで新しい libstdc++ を入れたら直った。
sudo add-apt-repository ppa:ubuntu-toolchain-r/test
sudo apt update
sudo apt install libstdc++-8-dev標準ライブラリとして llvm の libc++ を使う手もあるが、 ppa が用意されていないためインストールがめんどくさいのでやめた。
emacsの M-C-% (query-replace-regexp) で
a[index_i] = 1
のような文章を
a(index_i) = 1
文章に変換したいことがあった。
このとき置換のクエリを \[\([^[]]*\)\] -> (\1) としても上手く行かない。
なぜなら [^[]] という正規表現が表すのは (”[“以外の1文字 + “]”) であって、
こちらの意図する (”[”, “]” 以外の1文字) とは異なっているからだ。
じゃあエスケープすればいいだろうと思って [^\[\]] としてもこれも上手く行かない。
[] の中にある文字はエスケープされるので、
この正規表現が表すのは (”\”, “[”, “]” 以外の1文字) だからだ。
解決法はマニュアルにあった。
emacs の正規表現では [...] 中の最初の文字が “]” だった場合はエスケープされるとのことなので、
(”[”, “]” 以外の1文字)を表す正規表現は [^][] となる。
よって、冒頭の変換を行うクエリは \[\([^][*\)\] -> (\1) となる。
IEEE 754 方式の浮動小数点数ではビット列がある有理数をちょうど表す。
例えば倍精度の場合、
\begin{align*} \text{s} \ \text{e}_{11} \text{e}_{10} \ldots \text{e}_{1} \ \text{s}_{52} \text{s}_{51} \ldots \text{s}_{1} \end{align*}
というビット列が
\begin{align*} (-1)^{\text{s}} (1.\text{s}_{52} \text{s}_{51} \ldots \text{s}_{1})_{2} \times 2^{(\text{e}_{11} \text{e}_{10} \ldots \text{e}_{1})_{2} - 1023} \end{align*}
という有理数を表す。
前の記事で浮動小数点数を扱うクラスを作り、 別の記事で任意精度の有理数クラスを作ったので、 この２つを組み合わせれば浮動小数点数が表す有理数を計算することができる。
浮動小数点数を扱うクラスを拡張してこの処理を実装した。
#include <bitset>
#include <cassert>
#include <type_traits>
#include "big_int.hpp"
#include "fraction.hpp"
template <typename T>
class FloatingPoint;
using Single = FloatingPoint<float>;
using Double = FloatingPoint<double>;
template <typename...>
constexpr bool false_v = false;
// FloatingPoint assumes that sign bit, exponent bits, and significand bits are
// in this order.
template <typename T>
class FloatingPoint {
public:
constexpr static unsigned width_exp = [] {
if (std::is_same_v<T, float>) {
return 8;
} else if (std::is_same_v<T, double>) {
return 11;
} else {
assert(false_v<>);
}
}();
constexpr static unsigned width_sig = [] {
if (std::is_same_v<T, float>) {
return 23;
} else if (std::is_same_v<T, double>) {
return 52;
} else {
assert(false_v<>);
}
}();
constexpr static unsigned bias = [] {
if (std::is_same_v<T, float>) {
return 127;
} else if (std::is_same_v<T, double>) {
return 1023;
} else {
assert(false_v<>);
}
}();
constexpr static unsigned width = 1 + width_exp + width_sig;
private:
std::bitset<width> m_bits_all;
bool m_bit_sign;
std::bitset<width_exp> m_bits_exp;
std::bitset<width_sig> m_bits_sig;
T m_val;
static std::bitset<width> bits_of_val(T val);
static T val_of_bits(std::bitset<width> bits);
void set_bits() {
m_bit_sign = m_bits_all[width - 1];
for (unsigned i = 0; i < width_exp; ++i) {
m_bits_exp[i] = m_bits_all[width_sig + i];
}
for (unsigned i = 0; i < width_sig; ++i) {
m_bits_sig[i] = m_bits_all[i];
}
}
public:
FloatingPoint(T val) : m_bits_all(bits_of_val(val)), m_val(val) {
set_bits();
}
FloatingPoint(bool sign, std::bitset<width_exp> exp,
std::bitset<width_sig> sig)
: m_bits_all(std::to_string(sign) + exp.to_string() + sig.to_string()),
m_bit_sign(sign),
m_bits_exp(exp),
m_bits_sig(sig),
m_val(val_of_bits(m_bits_all)) {}
explicit FloatingPoint(const std::string &val)
: m_bits_all(val), m_val(val_of_bits(std::bitset<width>(val))) {
set_bits();
}
const T &value() { return m_val; }
std::bitset<width> bits_all() const { return m_bits_all; }
std::bitset<width_exp> bits_exp() const { return m_bits_exp; }
std::bitset<width_sig> bits_sig() const { return m_bits_sig; }
bool bit_sign() const { return m_bit_sign; }
// assumes that width_sig < 64 && width_exp < 64.
Fraction<BigInt> to_fraction() {
if (m_bits_all == std::bitset<width>(0)) {
return Fraction<BigInt>(0);
}
uint64_t exp = m_bits_exp.to_ullong();
uint64_t sig = m_bits_sig.to_ullong();
uint64_t total_bias = bias + width_sig;
if (exp > total_bias) {
BigInt a(sig + (1ul << width_sig));
BigInt b = BigInt(1).shift_left(exp - total_bias);
return m_bit_sign ? -Fraction<BigInt>(a * b) : Fraction<BigInt>(a * b);
} else {
BigInt numer(sig + (1ul << width_sig));
BigInt denom = BigInt(1).shift_left(total_bias - exp);
return m_bit_sign ? -Fraction<BigInt>(numer, denom)
: Fraction<BigInt>(numer, denom);
}
}
friend bool operator==(const FloatingPoint &lhs, const FloatingPoint &rhs) {
return rhs.bits_all() == lhs.bits_all();
}
};
template <>
float Single::val_of_bits(std::bitset<Single::width> val) {
auto temp = val.to_ulong();
char *c = reinterpret_cast<char *>(&temp);
return *reinterpret_cast<float *>(c);
}
template <>
std::bitset<Single::width> Single::bits_of_val(float val) {
char *c = reinterpret_cast<char *>(&val);
return *reinterpret_cast<uint32_t *>(c);
}
template <>
double Double::val_of_bits(std::bitset<Double::width> val) {
auto temp = val.to_ullong();
char *c = reinterpret_cast<char *>(&temp);
return *reinterpret_cast<double *>(c);
}
template <>
std::bitset<Double::width> Double::bits_of_val(double val) {
char *c = reinterpret_cast<char *>(&val);
return *reinterpret_cast<uint64_t *>(c);
}
以下のように使える。
#include "floating_point.hpp"
#include <iostream>
int main() {
std::cout << "1234.567f = " << Single(1234.567f).to_fraction() << "\n";
std::cout << "12345678.90123456 = " << Double(12345678.90123456).to_fraction() << "\n";
std::cout << "1.0e20f = " << Double(1.0e20f).to_fraction() << "\n";
std::cout << "1.0e50 = " << Double(1.0e50).to_fraction() << "\n";
return 0;
}
1234.567f = 10113573 / 8192
12345678.90123456 = 1657008972741239 / 134217728
1.0e20f = 100000002004087734272
1.0e50 = 100000000000000007629769841091887003294964970946560この答えがあっていることは出力された文字列をググれば分かる。
Python の Decimal や Haskell の toRational の使用例が出てくるので多分あっているだろう。
応用として、以下のようにすると減算で生じる桁落ち誤差を計算することができる。
#include "floating_point.hpp"
#include <iostream>
// return the error that a - b will have.
template <typename T>
Fraction<BigInt> error_by_sub(T a, T b) {
auto exact = FloatingPoint(a).to_fraction() - FloatingPoint(b).to_fraction();
return exact - FloatingPoint(a - b).to_fraction();
}
int main() {
double x = 1.3;
std::cout << "x = 1.3 = " << Double(x).to_fraction() << "\n";
std::cout << "error in 1.3 - x : " << error_by_sub(1.3, 1.3) << "\n";
std::cout << "error in 12.3 - x : " << error_by_sub(12.3, 1.3) << "\n";
std::cout << "error in 123.3 - x : " << error_by_sub(123.3, 1.3) << "\n";
std::cout << "error in 123456.3 - x : " << error_by_sub(123456.3, 1.3) << "\n";
std::cout << "error in 1.0e25 - x : " << error_by_sub(1.0e25, 1.3) << "\n";
return 0;
}
x = 1.3 = 5854679515581645 / 4503599627370496
error in 1.3 - x : 0
error in 12.3 - x : 3 / 4503599627370496
error in 123.3 - x : -13 / 4503599627370496
error in 123456.3 - x : 13107 / 4503599627370496
error in 1.0e25 - x : -5854679515581645 / 4503599627370496 前の記事で作った BigInt クラスが上手く動くことを確かめるために有理数クラスを作って動かした。
有理数クラスをテンプレート化して作る。
#include <iostream>
#include <ostream>
#include <stdexcept>
template <typename T>
T gcd_calc(T a, T b) {
T c;
while (a != 0) {
c = a;
a = b % a;
b = c;
}
return b;
}
template <typename T = int64_t>
class Fraction {
private:
// 0 is always represented with "+0 / 1"
bool m_sign;
T m_numer;
T m_denom;
Fraction(bool sign, T numerator, T denominator)
: m_sign(sign), m_numer(numerator), m_denom(denominator) {}
void reduce() {
auto gcd = gcd_calc(m_numer, m_denom);
m_numer /= gcd;
m_denom /= gcd;
}
public:
T numerator() const { return m_numer; }
T denominator() const { return m_denom; }
T sign() const { return m_sign ? -1 : 1; }
bool is_zero() const { return m_numer == 0; }
bool is_positive() const { return !is_zero() && !is_negative(); }
bool is_negative() const { return m_sign; }
Fraction() : Fraction(static_cast<T>(0)) {}
Fraction(T numerator, T denominator = 1)
: m_sign(numerator != 0l && ((numerator > 0l) != (denominator > 0l))),
m_numer(numerator < 0l ? -numerator : numerator),
m_denom(denominator < 0l ? -denominator : denominator) {
if (denominator == 0l) {
throw std::runtime_error("denominator must not be 0.");
}
if (m_numer == 0l) {
m_denom = 1;
}
reduce();
}
Fraction inv() const { return Fraction(m_sign, m_denom, m_numer); }
friend Fraction operator-(const Fraction &frac) {
if (frac.numerator() == static_cast<T>(0)) {
return frac;
}
return Fraction(!frac.m_sign, frac.numerator(), frac.denominator());
}
friend Fraction operator+(const Fraction &lhs, const Fraction &rhs) {
if (lhs == static_cast<T>(0)) {
return rhs;
}
if (rhs == static_cast<T>(0)) {
return lhs;
}
auto gcd = gcd_calc(lhs.denominator(), rhs.denominator());
auto denom = lhs.denominator() / gcd * rhs.denominator();
auto numer = lhs.sign() * lhs.numerator() * (denom / lhs.denominator()) +
rhs.sign() * rhs.numerator() * (denom / rhs.denominator());
return Fraction(numer, denom);
}
friend Fraction operator-(const Fraction &lhs, const Fraction &rhs) {
if (lhs == static_cast<T>(0)) {
return -rhs;
}
if (rhs == static_cast<T>(0)) {
return lhs;
}
auto gcd = gcd_calc(lhs.denominator(), rhs.denominator());
auto denom = lhs.denominator() / gcd * rhs.denominator();
auto numer = lhs.sign() * lhs.numerator() * (denom / lhs.denominator()) -
rhs.sign() * rhs.numerator() * (denom / rhs.denominator());
return Fraction(numer, denom);
}
friend Fraction operator*(const Fraction &lhs, const Fraction &rhs) {
if (lhs == static_cast<T>(0) || rhs == static_cast<T>(0)) {
return static_cast<T>(0);
}
auto gcd1 = gcd_calc(lhs.numerator(), rhs.denominator());
auto gcd2 = gcd_calc(lhs.denominator(), rhs.numerator());
auto numer = (lhs.numerator() / gcd1) * (rhs.numerator() / gcd2);
auto denom = (lhs.denominator() / gcd2) * (rhs.denominator() / gcd1);
return Fraction(lhs.m_sign != rhs.m_sign, numer, denom);
}
friend Fraction operator/(const Fraction &lhs, const Fraction &rhs) {
if (rhs == static_cast<T>(0)) {
throw std::runtime_error("Fraction::operator/ : devision by zero.");
}
if (lhs == static_cast<T>(0)) {
return static_cast<T>(0);
}
return lhs * rhs.inv();
}
Fraction &operator+=(const Fraction &other) {
*this = *this + other;
return *this;
}
Fraction &operator-=(const Fraction &other) {
*this = *this - other;
return *this;
}
Fraction &operator*=(const Fraction &other) {
*this = *this * other;
return *this;
}
Fraction &operator/=(const Fraction &other) {
*this = *this / other;
return *this;
}
friend bool operator==(const Fraction &lhs, const Fraction &rhs) {
return lhs.numerator() == rhs.numerator() &&
lhs.denominator() == rhs.denominator() && lhs.sign() == rhs.sign();
}
friend bool operator!=(const Fraction &lhs, const Fraction &rhs) {
return !(lhs == rhs);
}
friend bool operator<(const Fraction &lhs, const Fraction &rhs) {
if (lhs.is_negative() && rhs.is_positive()) {
return true;
}
if (lhs.is_positive() && rhs.is_negative()) {
return false;
}
return lhs.sign() * lhs.numerator() * rhs.denominator() <
lhs.sign() * lhs.denominator() * rhs.numerator();
}
friend bool operator>(const Fraction &lhs, const Fraction &rhs) {
if (lhs.is_negative() && rhs.is_positive()) {
return false;
}
if (lhs.is_positive() && rhs.is_negative()) {
return true;
}
return lhs.sign() * lhs.numerator() * rhs.denominator() >
lhs.sign() * lhs.denominator() * rhs.numerator();
}
friend bool operator<=(const Fraction &lhs, const Fraction &rhs) {
return lhs == rhs || lhs < rhs;
}
friend bool operator>=(const Fraction &lhs, const Fraction &rhs) {
return lhs == rhs || lhs > rhs;
}
friend std::ostream &operator<<(std::ostream &os, const Fraction &rhs) {
if (rhs.denominator() == 1 || rhs.numerator() == 0) {
os << rhs.numerator();
} else {
if (rhs.is_negative()) {
os << "-";
}
os << rhs.numerator() << " / " << rhs.denominator();
}
return os;
}
};
この有理クラスを使って、級数
\begin{align*} \frac{\pi}{4} = 1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} \cdots \end{align*}
を任意精度で計算する。
BigInt の実装が上手く行っていれば、以下のコードで正しく計算ができるはずだ。
#include <iostream>
#include "fraction.hpp"
#include "big_int.hpp" // 前の記事を参照
int main() {
using Frac = Fraction<BigInt>;
Frac sum(0);
for (int i = 0; i< 100; ++i) {
if (i % 2 == 0) {
sum += Frac(BigInt(4), BigInt(2 * i + 1));
} else {
sum -= Frac(BigInt(4), BigInt(2 * i + 1));
}
}
std::cout << sum;
return 0;
}
このプログラムのコンパイルは無事通り、結果は
8252079759413970386664454687621174435983101115012912631997769614579677862845786070667088 / 2635106162757236442495826303084698495565581115509040892412867358728390766099042109898375となる。
この答えがあっているかどうかは出力結果をググるとなんとなく分かる。 同じ計算をしたサイトが出てくるので、多分あっているだろう。