|
いろいろあってBigDecimalを実装したので公開します。
バグってたら教えて下さい。
verify
参考文献
// 固定小数点で実装されている // 内部はlong longの配列で下31bitのみ使用 // ゼロ除算・オーバーフロー等は例外を投げる // 比較はデフォルトでEPSILON付き // EPSILONはライブラリ内部用でかなり小さいため、呼び出す側でEPSを指定したほうが良い // doubleから変換すると結局doubleの精度になってしまうためstringから変換することを推奨 #include<iostream> #include<cmath> using namespace std; #define rep(i, n) for (int i = 0; i < (int)(n); ++i) // ID*2^31,DD*2^31で実装されている // 桁数は必要最低限の倍以上を推奨 // 整数部桁数 const int ID = 6; // 小数部桁数 const int DD = 9; const int BS = 31; const bool PLUS = false; const bool MINUS = true; struct BD { bool sign; long long d[ID + DD]; BD(); ~BD(); BD(int); BD(long long); BD(string); BD(double); BD normal(); BD operator-() const; BD operator<<(int) const; BD operator>>(int) const; BD operator+(const BD&) const; BD operator-(const BD&) const; BD operator*(unsigned int) const; BD operator*(const BD&) const; BD operator/(const BD&) const; BD operator%(const BD&) const; BD operator<<=(int); BD operator>>=(int); BD operator+=(const BD&); BD operator-=(const BD&); BD operator*=(unsigned int); BD operator*=(const BD&); BD operator/=(const BD&); BD operator%=(const BD&); BD operator++(); BD operator++(int); BD operator--(); BD operator--(int); // EPSILONを考慮した比較 bool operator==(const BD&) const; bool operator!=(const BD&) const; bool operator<(const BD&) const; bool operator>(const BD&) const; bool operator<=(const BD&) const; bool operator>=(const BD&) const; // EPSILONなしの厳密な比較 // 0になるまで繰り返すなどの処理に用いる bool equals(const BD&) const; int toInt() const; long long toLongLong() const; string toString(int, string) const; double toDouble() const; }; ostream &operator<<(ostream&, BD); istream &operator>>(istream&, BD&); // cmathから基本的な関数を抽出 BD sin(BD); BD cos(BD); BD tan(const BD&); BD asin(BD); BD acos(const BD&); BD atan(BD); BD atan2(const BD&, const BD&); BD sinh(BD); BD cosh(BD); BD tanh(const BD&); BD exp(const BD&); BD log(const BD&); BD log10(const BD&); BD pow(const BD&, const BD&); BD sqrt(const BD&); BD abs(const BD&); BD ceil(const BD&); BD floor(BD); // 何度もコンストラクタを呼ぶと無駄なので定数にしてしまう // ライブラリ外では考慮する必要はほぼない const BD BD0 = BD(0); const BD BD1 = BD(1); const BD BD2 = BD(2); const BD BD3 = BD(3); const BD REV10 = BD1 / BD(10); const BD PI = (atan(BD1 / BD(5)) << 4) - (atan(BD1 / BD(239)) << 2); // EPSILONの値はかなり怪しいので適宜調整が必要 const BD EPSILON = BD1 >> (BS - 4) * DD; // complexでデフォルトコンストラクタは0を返すことが仮定されている BD::BD() {*this = BD0;} BD::~BD() {} BD::BD(int n) { sign = PLUS; rep (i, ID + DD) d[i] = 0; d[DD] = n; normal(); } BD::BD(long long n) { sign = PLUS; rep (i, ID + DD) d[i] = 0; d[DD] = n; normal(); } // TODO どう見てもオーバーフローするので後で修正 BD::BD(string str) { *this = BD0; bool minus = false; if (str[0] == '-') { minus = true; str = str.substr(1); } int rp = 0; rep (i, str.size()) { if (str[i] == '.') rp = str.size() - i - 1; else *this = *this * 10 + (str[i] - '0'); } while (rp--) *this *= REV10; if (minus) sign = MINUS; } BD::BD(double r) { *this = BD0; int n; BD bd = r >= 0 ? BD1 : -BD1; r = 2 * abs(frexp(abs(r), &n)); bd <<= n - 1; rep (i, 64) { if (r >= 1) { *this += bd; r -= 1; } r *= 2; bd >>= 1; } } BD BD::normal() { rep (i, ID + DD - 1){ d[i + 1] += d[i] >> BS; d[i] &= (1LL << BS) - 1; } if (d[ID + DD - 1] < 0) { sign = !sign; rep (i, ID + DD) d[i] = -d[i]; normal(); } if (d[ID + DD - 1] >= (1LL << BS)) throw "overflow"; rep (i, ID + DD) if (d[i] != 0) return *this; sign = PLUS; return *this; } BD BD::operator-() const { BD bd(*this); bd.sign = !bd.sign; return bd; } BD BD::operator<<(int a) const {return BD(*this) <<= a;} BD BD::operator>>(int a) const {return BD(*this) >>= a;} BD BD::operator+(const BD &a) const {return BD(*this) += a;} BD BD::operator-(const BD &a) const {return BD(*this) -= a;} BD BD::operator*(unsigned int a) const {return BD(*this) *= a;} BD BD::operator*(const BD &a) const {return BD(*this) *= a;} BD BD::operator/(const BD &a) const {return BD(*this) /= a;} BD BD::operator%(const BD &a) const {return BD(*this) %= a;} BD BD::operator<<=(int a) { if (a < 0) return *this >>= -a; while (a >= BS) { if (d[ID + DD - 1]) throw "overflow"; for (int i = ID + DD; --i > 0; ) d[i] = d[i - 1]; d[0] = 0; a -= BS; } if (d[ID + DD - 1] >= (1LL << (BS - a))) throw "overflow"; rep (i, ID + DD) d[i] <<= a; return normal(); } BD BD::operator>>=(int a) { if (a < 0) return *this <<= -a; while (a >= BS) { rep (i, ID + DD - 1) d[i] = d[i + 1]; d[ID + DD - 1] >>= BS; a -= BS; } rep (i, ID + DD - 1) { d[i] |= d[i + 1] << BS; d[i + 1] = 0; } rep (i, ID + DD) d[i] >>= a; return normal(); } BD BD::operator+=(const BD &a) { if (sign == a.sign) rep (i, ID + DD) d[i] += a.d[i]; else rep (i, ID + DD) d[i] -= a.d[i]; return normal(); } BD BD::operator-=(const BD &a) {return *this += -a;} BD BD::operator*=(unsigned int a) { rep (i, ID + DD) d[i] *= a; return normal(); } BD BD::operator*=(const BD &a) { BD res = BD0; rep (i, ID + DD) { if (i < DD) res = (res + *this * a.d[i]) >> BS; else res += *this * a.d[i] << (i - DD) * BS; } res.sign = (sign == a.sign ? PLUS : MINUS); return *this = res.normal(); } BD BD::operator/=(const BD &a) { if (a == BD0) throw "divide by zero"; BD rev = (double)1 / a.toDouble(); rep (i, 7) rev = (rev << 1) - a * rev * rev; rev.sign = a.sign; return *this *= rev; } BD BD::operator%=(const BD &a) { if (a == BD0) throw "modulo by zero"; return *this -= floor(*this / a) * a; } BD BD::operator++() { return *this += 1; } BD BD::operator++(int) { BD bd = *this; *this += 1; return bd; } BD BD::operator--() { return *this -= 1; } BD BD::operator--(int) { BD bd = *this; *this -= 1; return bd; } // TODO 比較周りの整合性がとれていない可能性があるので後で調べる bool BD::operator==(const BD &a) const {return *this <= a && a <= *this;} bool BD::operator!=(const BD &a) const {return !(*this == a);} bool BD::operator<(const BD &a) const { if (*this == a) return false; BD aa = a - EPSILON; if (sign == MINUS) { if (aa.sign == MINUS) return -a < -*this; else return true; } if (aa.sign == MINUS) return false; for (int i = ID + DD; i-- > 0; ) if (d[i] != aa.d[i]) return d[i] < aa.d[i]; return false; } bool BD::operator>(const BD &a) const {return a < *this;} bool BD::operator<=(const BD &a) const { BD aa = a + EPSILON; if (sign == MINUS) { if (aa.sign == MINUS) return -aa <= -*this; else return true; } if (aa.sign == MINUS) return false; for (int i = ID + DD; i-- > 0; ) if (d[i] != aa.d[i]) return d[i] < aa.d[i]; return true; } bool BD::operator>=(const BD &a) const {return a <= *this;} bool BD::equals(const BD &a) const { if (sign != a.sign) return false; rep (i, ID + DD) if (d[i] != a.d[i]) return false; return true; } int BD::toInt() const { int res = 0; rep (i, ID) res += d[DD + i] << BS * i; if (sign == MINUS) return -res; return res; } long long BD::toLongLong() const { long long res = 0; rep (i, ID) res += d[DD + i] << BS * i; if (sign == MINUS) return -res; return res; } string BD::toString(int digit = 100, string mode = "near") const { string str; BD a = *this, bd = BD1; if (a.sign == MINUS) { str += "-"; a.sign = PLUS; } if (mode == "near") { BD round = BD1 >> 1; rep (i, digit) round *= REV10; a += round + EPSILON; } if (mode == "ceil") { BD round = BD1; rep (i, digit) round *= REV10; a += round - EPSILON; } for (; bd <= a; bd *= 10) ++digit; if (bd > BD1) { bd *= REV10; --digit; } rep (i, digit + 1) { if (bd == BD0) { str += "0"; continue; } if (bd == REV10) str += "."; int d = 0; while (bd < a) { a -= bd; ++d; } if (d > 9) { d -= 10; string::iterator itr = str.end(); while (true) { if (itr == str.begin()) { str = "1" + str; break; } --itr; if (*itr == '.') continue; ++*itr; if (*itr > '9') *itr = '0'; else break; } } str += '0' + d; bd *= REV10; } return str; } double BD::toDouble() const { double res = 0; rep (i, ID + DD) res += d[i] * pow(2, (i - DD) * BS); if (sign == MINUS) return -res; return res; } ostream &operator<<(ostream &os, BD a) { os << a.toString(); return os; } istream &operator>>(istream &is, BD &a) { string str; is >> str; a = BD(str); return is; } BD sin(BD x) { BD res = BD0, xx = - x * x; for (BD i = BD1; ; i += BD2) { x /= max(i * (i - 1), BD1); if (x.equals(BD0)) break; res += x; x *= xx; } return res; } BD cos(BD x) { BD res = BD0, xx = - x * x; x = BD1; for (BD i = BD0; ; i += BD2) { x /= max(i * (i - BD1), BD1); if (x.equals(BD0)) break; res += x; x *= xx; } return res; } BD tan(const BD &x) {return sin(x) / cos(x);} BD asin(BD x) { if (abs(x) > BD1) throw "out of domain"; if (x > BD1 / sqrt(BD2)) return (PI >> 1) - asin(sqrt(BD1 - x * x)); if (x < -BD1 / sqrt(BD2)) return -(PI >> 1) + asin(sqrt(BD1 - x * x)); BD res = BD0, xx = x * x >> 2; for (BD i = BD0; ; ++i) { x *= max(i * BD2 * (i * BD2 - BD1), BD1); x /= max(i * i, BD1); BD add = x / (i * BD2 + BD1); if (add.equals(BD0)) break; res += add; x *= xx; } return res; } BD acos(const BD &x) {return (PI >> 1) - asin(x);} BD atan(BD x) { if (x.sign == MINUS) return -atan(-x); if (abs(x) > sqrt(BD2) + BD1) return (PI >> 1) - atan(BD1 / x); if (abs(x) > sqrt(BD2) - BD1) return (PI >> 2) + atan((x - BD1) / (x + BD1)); BD res = BD0, xx = - x * x; for (BD i = BD1; ; i += BD2) { BD add = x / i; if (add.equals(BD0)) break; res += add; x *= xx; } return res; } BD atan2(const BD &y, const BD &x) { if (x == BD0) { if (y > BD0) return PI / BD2; if (y < BD0) return -PI / BD2; throw "origin can't define argument"; } if (x.sign == PLUS) return atan(y / x); if (y.sign == PLUS) return atan(y / x) + PI; return atan(y / x) - PI; } BD sinh(BD x) { BD res = BD0, xx = x * x; for (BD i = BD1; ; i += BD2) { x /= max(i * (i - BD1), BD1); if (x.equals(BD0)) break; res += x; x *= xx; } return res; } BD cosh(BD x) { BD res = BD0, xx = x * x; x = BD1; for (BD i = BD0; ; i += BD2) { x /= max(i * (i - BD1), BD1); if (x.equals(BD0)) break; res += x; x *= xx; } return res; } BD tanh(const BD &x) {return sinh(x) / cosh(x);} BD exp(const BD &x) { BD res = BD0, xx = BD1; for (int i = 0; ; ++i) { xx /= max(i, 1); if (xx.equals(BD0)) break; res += xx; xx *= x; } return res; } BD log(const BD &x) { BD y = log(x.toDouble()); BD z = x / exp(y); BD a = (z - 1) / (z + 1); BD res = BD0, b = a, aa = a * a; for (int i = 1; ; i += 2) { if (b.equals(BD0)) break; res += b / i; b *= aa; } return y + res * BD2; } BD log10(const BD &x) {return log(x) / log(BD(10));} BD pow(const BD &x, const BD &y) { if (x.sign == MINUS) { if (floor(y) == y) return floor(y).d[DD] % 2 ? -pow(-x, floor(y)) : pow(-x, floor(y)); throw "power of negative number"; } return exp(y * log(x)); } BD sqrt(const BD &x) { BD r = 1 / sqrt(x.toDouble()); rep (i, 7) r *= (BD3 - x * r * r) >> 1; return BD1 / r; } BD abs(const BD &x) {return x.sign == PLUS ? x : -x;} BD ceil(const BD &x) { if (x.sign == MINUS) return -floor(-x); BD f = floor(x); return f == x ? f : x + 1; } BD floor(BD x) { if (x.sign == MINUS) return -ceil(-x); x += EPSILON; rep (i, DD) x.d[i] = 0; return x; }
