Go to Qiita Advent Calendar 2024 Top
LoginSignup
1
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@NokonoKotlin(おこ ていゆ)

C++ で Fraction (有理数) ライブラリを作ってみた

Posted at

まえがき

この記事は投稿者(NokonoKotlin)の個人サイトの記事から Qiita 用に移植 & 加筆したものです。 (この文言は記事が剽窃でないことを周知するためのものであり、個人サイトの宣伝ではないことをあらかじめご了承ください。)

はじめに

ABC385 の F 問題 で、数値計算の誤差に注意する場面がありました。この問題は __float128 などでゴリゴリ計算すれば良いのですが、この機会にと思って自分の Fraction ライブラリ (有理数ライブラリ) を放出です。

有理数ライブラリでは、分子と分母をわけて計算することで、計算結果を double などの数値にキャストする場面以外で誤差を出さずに計算を重なうことができます。また、特殊な場面ではありますが、$0$ 徐算を行なってもエラーを出さず、かつ計算結果から $0$ 徐算の存在がわかるような設計のため、条件分岐の手間などがなくなります。また $\infty$ などの特殊なインスタンスも導入してみました。

概要

整数型 (__int128_t などのオーバーフローしづらいもの) の $A,B$ を用いて有理数 $\frac{A}{B}$ を表現します (__int128_t だとまだ普通にオーバーフローしやすいから注意)。

  • 値を $\mathrm{T}$ 型で取り出したい場合は to_real<T>() を呼び出すか、$\mathrm{T}$ 型に明示的にキャストする ($\mathrm{T}$ が整数ならあまりは切り捨て)。
  • $A,B$ を取り出したいときは raw_pair() でアクセス。
  • 逆数を返す inv 関数や、絶対値を返す abs 関数、特殊なインスタンス $\infty,\varepsilon$ などを返す関数や、特殊なインスタンスかどうかをチェックする関数も搭載。

Fraction がとる値

冒頭で説明した通り、Fraction は整数 $A,B$ で $\frac{A}{B}$ の形で定義されます。Fraction は通常の数値に加え、無限大/微少量 などの値をとりうる特殊な型です。特殊な値 (インスタンス) は以下の通りです。

  • $\mathrm{amb} := \frac{0}{0}$ (ゼロ徐算や無限大同士の減算など)
  • $\mathrm{inf} := \infty$ (無限大)
  • $\varepsilon := \frac{1}{\infty}$ (微少量)
$\mathrm{amb},\mathrm{inf}$ は特殊インスタンスとして明示的に宣言/定義しますが、$\varepsilon$ は $\mathrm{inf}$ から導出されるので明示的に宣言/定義しません。

表現と計算

実数有理数の計算であれば、ただやるだけです。Fraction クラスのメインは、先に定義した $\mathrm{amb},\mathrm{inf},\varepsilon$ に対して計算規則をうまく決めて、矛盾がないような 有理数+特殊インスタンス の計算を実現することです。

Fraction クラスによる有理数 + $\alpha$ の表現は、メンバ $A,B$ および補助的なメンバ $\mathrm{eps\_sign}$ によって以下で定められます。

  1. $(A,B) = (0,0)$ の場合
    • $\frac{A}{B} = $ 不定値 ($\mathrm{amb}$) として定義する。
  2. $(A,B) = (0,1)$ の場合、$\mathrm{eps\_sign}$ というメンバの値で場合分けする。
    • $\mathrm{eps\_sign} = 0$ なら $\frac{A}{B} = 0$
    • $\mathrm{eps\_sign} = 1$ なら $\frac{A}{B} = \varepsilon$
    • $\mathrm{eps\_sign} = -1$ なら $\frac{A}{B} = -\varepsilon$
  3. $(A,B) = (\pm 1 ,0)$ の場合、$A$ の符号で場合分けする。
    • $A = 1$ なら $\frac{A}{B} = \infty$
    • $A = -1$ なら $\frac{A}{B} = -\infty$
    $|A|$ が $0$ 以外なら、符号だけ保存して $|A| = 1$ の形にフォーマットする。
  4. 上記の $1,2,3$ のいずれでもない場合、$\frac{A}{B}$ は実数有理数を表す。ただし符号は $A$ につける方針で統一する。

計算規則

Fraction オブジェクト $x,y$ に対して、数学的に矛盾がないように四則演算を定義します。結果が $\mathrm{amb}$ になる場合の計算さえ分かれば、あとは数学の規則に則って計算するだけです。

  1. $\mathrm{amb}$ に対してどの演算を行っても結果は $\mathrm{amb}$
  2. $x + y$ の計算
    • $\infty$ と $-\infty$ が混在するなら $x + y = \mathrm{amb}$
    • $\varepsilon$ と $-\varepsilon$ が混在するなら $x + y = \mathrm{amb}$
  3. $x \times y$ の計算
    • 片方が $\pm \varepsilon$ で他方が $\pm \infty$ なら $x \times y = \mathrm{amb}$
  4. $x-y$ の計算は $x+(-y)$ の計算と等しいです。$-y$ は $y\times (-1)$ です。
  5. $\frac{x}{y}$ は $x \times y^{-1}$ なので $y$ の逆数が分かればよく、$y$ の逆数は以下の計算規則で計算可能です。
    • $y = 0$ または $y = \mathrm{amb}$ なら $ y^{-1} =\mathrm{amb}$
    • $y = \pm\varepsilon $ なら符号に合わせた $y^{-1} =\infty$
    • $y = \pm\infty $ なら符号に合わせた $y^{-1} =\varepsilon$
    • 上記のいずれでもないなら $y$ のメンバ $A,B$ をスワップし、符号を $A$ につけたものが $y^{-1}$。

設計

キャストについて

  • 自作クラスの機能を優先させたいので、C++ 組み込み型 $\rightarrow$ Fraction のキャストは暗黙的。
  • Fraction $\rightarrow$ 別の何か のキャストは明示的。

メンバの持ち方について

以下の制約を満たすように $A,B$ を整形する format 関数を用意し、コンストラクタなどでの初期化時に呼び出す。演算などを行うたびに Fraction オブジェクトを生成し、コンストラクタを通すことで常に format されるようにする。

  • $\frac{A}{B}$ が有理数の場合、$\frac{A}{B}$ は常に既約であるように設計 (gcd で割る)
  • 符号は $A$ につける。分母 $B$ は常に非負整数。

比較演算

三方演算子 std::partial_ordering <=>(T a , T b) を使います (C++ 20 からの機能)。

  • <=> の実装が返す以下の結果から、< , <= , > , >= 演算の結果を作るすごいやつ。
    • std::partial_ordering::equivalent : a と b が等しい
    • std::partial_ordering::less : a が b 未満
    • std::partial_ordering::greater : a が b より大きい
    • std::partial_ordering::unordered : a と b の順序が未定義
  • 例えば、0 <=> ∞ が less を返すように実装する。すると、0 < ∞ の結果が true になる。
  • unordered について
    • a,b の順序比較ができない時、unordered を返し、全ての順序演算が false を返す(無限大同士の比較など)。
      • unordered な順序があるバブルソートは永遠に終わらない
      • unordered な順序がある列をソートすると、unordered な順序以外の要素も正しくソートされない
    • できるだけ unordered な順序がある要素を比較しない

その他

  • 各演算の両辺は必ず整数 or Fraction
  • 演算の規則は一般の数学の規則に準拠
  • コンストラクタで明示的に分子が非ゼロ、分母がゼロを与えた場合、それは無限大/無限小 になる。

ソースコード

使い方は後述します。あとコメントの英語が拙いのは気にしないでください。雰囲気がわかればいいのです。

ライセンスは守ってくださいね (このライブラリの冒頭に以下を書くだけで良いです)。

ライセンスに従わない利用を確認した場合、ライセンスに従わない利用をやめるよう指示するなどの処置を取る場合があります。

#include<cassert>
#include<iostream>
#include<compare>

using std::make_pair,std::pair;
using std::partial_ordering;
using std::istream,std::ostream;

/*   
    このコメントは消さないでください。
    Don't remove this comment!!!

    Copyright ©️ (c) NokonoKotlin (okoteiyu) 2024. 
    Released under the MIT license(https://opensource.org/licenses/mit-license.php) 
*/
class Fraction{
    protected:
    using Integer = __int128_t;
    Integer A,B;
    Integer m_eps_sign = 0;
    Integer gcd(Integer x, Integer y) { return (y!=0) ? gcd(y, x % y) : x; }

    void format(){
        if(A == 0 && B == 0){}
        else if(A == 0)B = 1; 
        else if(B == 0){
            Integer sign_ = 1;
            if(this->A < 0)sign_ = -1;
            this-> A = sign_;
        }else{ 
            Integer g = gcd(A,B);
            A/=g;B/=g;
            if(B < 0){A*=-1;B*=-1;} 
        }
    }

    public:
    inline static const Fraction inf(){return Fraction(1,0);}
    inline static const Fraction ninf(){return Fraction(-1,0);}
    inline static const Fraction amb(){return Fraction(0,0);}
    inline static const Fraction zero(){return Fraction(0,1);}

    Fraction(){}
    Fraction(Integer a):A(a),B(1),m_eps_sign(0){this->format();}
    Fraction(Integer a , Integer b):A(a),B(b),m_eps_sign(0){this->format();}
    Fraction(const Fraction& x):A(x.A),B(x.B),m_eps_sign(x.m_eps_sign){this->format();}
    Fraction& operator =(const Fraction& x){
        this->m_eps_sign = x.m_eps_sign;
        this->A = x.A;
        this->B = x.B;
        this->format();
        return (*this);
    } 

    Fraction abs() const& {return (*this)*this->sign();}
    bool is_ambiguous() const& {return bool(this->A == 0 && this->B == 0);}
    bool is_infinite() const& {return bool(this->B==0);}
    Integer sign() const& {
        assert(!this->is_ambiguous());
        if(this->A == 0)return this->m_eps_sign;
        if(this->A < 0)return -1;
        return 1;
    }

    template<typename real>constexpr real to_real() const& {return real(A)/real(B);}
    pair<Integer,Integer> raw_pair() const& {return make_pair(A,B);}
    template<typename real>explicit operator real() const& {return real(A)/real(B);}
    Fraction inv() const& {
        if(this->is_ambiguous())return amb();
        if(this->A == 0){
            if(this->sign() == 0)return amb();
            return Fraction(this->sign(),0);
        }
        if(this->B == 0){
            Fraction res(0,1);
            res.m_eps_sign = this->sign();
            return res;
        }
        return Fraction(this->B,this->A);
    }

    
    Fraction& operator +=(Fraction x){
        if(this->is_ambiguous() || x.is_ambiguous()) [[unlikely]] (*this) = amb();
        else if(this->is_infinite() && x.is_infinite()) [[unlikely]] {
            if(this->sign() != x.sign())(*this) = amb();
            else if(this->sign() == -1)(*this) = ninf();
            else (*this) = inf();
        }else if(this->is_infinite()) [[unlikely]] {
        }else if(x.is_infinite()) [[unlikely]] (*this) = x;
        else if(this->A == 0 && x.A == 0) [[unlikely]] {
            if(this->sign() == x.sign()){}
            else if(this->sign() == 0)(*this) = x;
            else if(x.sign() == 0){}
            else (*this) = amb();
        }else [[likely]]{
            (*this) = Fraction(this->A*x.B + this->B*x.A , this->B*x.B);
        }
        return (*this);
    }
    Fraction& operator *=(Fraction x){
        if(this->is_ambiguous() || x.is_ambiguous()) [[unlikely]] (*this) = amb();
        else if(this->is_infinite() || x.is_infinite()) [[unlikely]] {
            if(this->A != 0 && x.A != 0){
                if(this->sign() == x.sign())(*this) = inf();
                else (*this) = ninf();
            }else{
                if(this->sign() == 0 || x.sign() == 0)(*this) = 0;
                else (*this) = amb();
            }
        }else if(this->A == 0 || x.A == 0) [[unlikely]] {
            if(this->sign() == 0 || x.sign() == 0)(*this) = 0;
            else {
                Fraction new_eps = 0;
                new_eps.m_eps_sign = this->sign()*x.sign();
                (*this) = new_eps;
            }
        }else [[likely]] {(*this) = Fraction( this->A*x.A, this->B*x.B );}
        return (*this);
    }
    Fraction operator -(){Fraction res(*this);return res*Fraction(-1);}
    Fraction& operator -=(Fraction x){return ((*this) += (-x));}
    Fraction& operator /=(Fraction x){return ((*this) *= x.inv());}
    Fraction& operator %=(Fraction x){
        assert(!this->is_ambiguous() && !this->is_infinite() && !x.is_ambiguous() && !x.is_infinite());
        assert(x.A!= 0);
        Fraction tmp = (*this);
        tmp/=x;
        Integer m = tmp.raw_pair().first/tmp.raw_pair().second;
        (*this) -= m*x;
        return (*this);
    }

    friend Fraction operator +(Fraction a , Fraction b){return a+=b;}
    friend Fraction operator -(Fraction a , Fraction b){return a-=b;}
    friend Fraction operator *(Fraction a , Fraction b){return a*=b;}
    friend Fraction operator /(Fraction a , Fraction b){return a/=b;}
    friend Fraction operator %(Fraction a , Fraction b){return a%=b;}

    friend partial_ordering operator <=>(const Fraction& x , const Fraction& y){
        if(x.is_ambiguous() || y.is_ambiguous()) [[unlikely]] return partial_ordering::unordered;
        else if(x.sign() != y.sign()) [[unlikely]] return (x.sign() <=> y.sign());
        else if(x.A == 0 || y.A == 0) [[unlikely]] {
            if(x.sign() == 0 && y.sign() == 0)return partial_ordering::equivalent;
            if(x.A == 0 && y.A == 0)return partial_ordering::equivalent; 
            else if(x.A == 0){
                if(y.sign() == -1)return partial_ordering::greater;
                return partial_ordering::less;
            }else{
                if(x.sign() == -1)return partial_ordering::less;
                return partial_ordering::greater;
            }
        }
        else if(x.is_infinite() && y.is_infinite()) [[unlikely]] return partial_ordering::equivalent;
        else if(x.is_infinite()) [[unlikely]] {
            if(x.sign() == -1)return partial_ordering::less;
            return partial_ordering::greater;
        }else if(y.is_infinite()) [[unlikely]] {
            if(y.sign() == -1)return partial_ordering::greater;
            return partial_ordering::less;
        }
        else [[likely]] return ((x.A)*(y.B) <=> (y.A)*(x.B));
    }

    friend bool operator ==(const Fraction& x,const Fraction& y){
        if(x.is_ambiguous() || x.is_infinite() || y.is_ambiguous() || y.is_infinite())return false;
        if(x.A == 0 && y.A == 0){
            if(x.sign() == 0 && y.sign() == 0)return true;
            return false;
        }
        return bool(x.A == y.A && x.B == y.B);
    }
    friend bool operator !=(const Fraction& x,const Fraction& y){return !(x==y);}

    friend ostream& operator << (ostream& os, const Fraction& x){
        os.precision(40);
        if(x.is_ambiguous())os << "?";
        else if(x.A == 0){
            if(x.sign()<0)os << "-ε";
            else if(x.sign()>0)os << "ε";
            else os << 0;
        }else if(x.is_infinite()){
            if(x.sign() < 0)os << "-∞";
            else os << "∞";
        }else os << double(x.A) << "/" << (long long)(x.B);
        return os;
    }  
    friend istream& operator >> (istream& is, Fraction& x){
        long long tmp;
        is >> tmp;
        x = Fraction(tmp);
        return is;
    }
    /*
        このコメントは消さないでください。
        Don't remove this comment!!!

        Copyright ©️ (c) NokonoKotlin (okoteiyu) 2024. 
        Released under the MIT license(https://opensource.org/licenses/mit-license.php) 
    */
};
1
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?