スキップしてメイン コンテンツに移動

多項式 in Rust

環 R 上の多項式 R[x] とは, 有限個の R の元の組の全体 R[x] に以下の代数構造を入れたもの.
  • 多項式 p(x) in R[x] と元 r in R の和 p(x) + r
  • 多項式 p(x), q(x) in R[x] の和 p(x) q(x)
  • 多項式 p(x) in R[x] と元 r in R の積 r p(x)
  • 多項式 p(x), q(x) in R[x] の和 p(x) + q(x)
  • 多項式 p(x) in R[x] に元 r in R を代入する: p(r)
  • 多項式 p(x) in R[x] に多項式 q in R[x] を代入する: f(q(x))
R[x] はこの和と積に関して環をなす. より詳細な多項式環の性質については wikipedia 参照.
これを Rust で実装する. ジェネリクスとトレイト境界のよい練習になった.
特に多項式への代入として環 R の元と多項式 R[x] のどちらをも許す部分で手間取った. 若干 clone を使い過ぎた気もするけど, 他にどうしようもなさそう.

/// 多項式を表す構造体.
#[derive(Clone, Debug, PartialEq, Eq)]
pub struct Polynomial<T> {
    coef: Vec<T>, 
}



extern crate num_traits;
use num_traits::Zero;

use std::ops::{Add, Sub, Neg, Mul};
impl<T> Polynomial<T> where T: Clone + Zero + PartialEq
{
    /// 係数の配列から新しく多項式を生成する. 係数はタプルでもベクトルでもよい.
    /// 係数 [ a_0, a_1, ... ] の $n$ 番目の要素が $n$ 次項の係数となる.
    ///
    /// 引数の配列が要素を持たない可能性を考慮し, `Polynomial::new` 関数の戻り値
    /// は `Option<Polynomial>` である. 使用する前に `.unwrap()` する.
    /// 
    /// ```
    /// use polynomial::Polynomial;
    /// 
    /// let coef: [f64;3] = [ 9., 3., 1. ];
    /// let polynomial = Polynomial::new(&coef).unwrap(); // x^2 + 3x + 9
    ///
    /// let n: usize = 16;
    /// let mut coef: Vec<f64> = vec![ 0.; n ];
    /// coef[n-1] = 4.;
    /// let polynomial = Polynomial::new(&coef).unwrap(); // 4 x^{15}
    /// ```
    /// 
    pub fn new( coef: &[T] ) -> Option<Self> {
        if coef.len() == 0 {
            None
        } else {
            let mut p = Polynomial { coef: coef.to_vec() };
            concat( &mut p );
            Some(p)
        }
    }

    /// 係数の数 (多項式の次数 + 1) を返すメソッド. ゼロ多項式の場合も 1 になる.
    ///
    /// ```
    /// use polynomial::Polynomial;
    /// 
    /// let p1 = Polynomial::new(&[ 0, 0, 1 ]).unwrap();
    /// assert_eq!( p1.len(), 3 );
    /// 
    /// let p2 = Polynomial::new(&[ 0, 0, 0 ]).unwrap();
    /// assert_eq!( p2.len(), 1 );
    /// ```
    pub fn len(&self) -> usize {
        self.coef.len()
    }
}

impl<T> Polynomial<T> where T: Clone + Zero + Add + Mul<Output=T> {
    /// 与えられた引数 `x` を代入したときの多項式の値を返す.
    /// 引数はスカラー値でも多項式でもよい.
    pub fn eval<X>(&self, x: X) -> X
        where X: Clone + Zero + Add<T, Output=X> + Mul<T, Output=X> + Mul<X, Output=X>
    {
        let mut value: X = X::zero();
        for c in self.coef.iter().rev() {
            value = value.clone() * x.clone() + c.clone();
        }
        value
    }
}


impl<T> Zero for Polynomial<T> where T: Clone + Zero + PartialEq {
    fn zero() -> Polynomial<T> {
        Polynomial { coef: vec![ T::zero() ] }
    }
    fn is_zero(&self) -> bool {
        self == &Self::zero()
    }
}

// Add トレイトの実装 (スカラー値を足す)
impl<T> Add<T> for Polynomial<T> where T: Clone + Add<Output=T>
{
    type Output = Polynomial<T>;
   
    fn add(self, other: T) -> Polynomial<T> {
        let mut poly = self.clone();
        poly.coef[0] = poly.coef[0].clone() + other;

        poly
    }
}
// Add トレイトの実装 (多項式を足す)
impl<T> Add<Polynomial<T>> for Polynomial<T> where T: Clone + Zero + Add + PartialEq
{
    type Output = Polynomial<T>;

    fn add(self, other: Self) -> Self {
        let mut poly = Polynomial { coef: Vec::new() };

        for i in 0..self.coef.len().max(other.coef.len()) {
            let a: T = match self.coef.get(i) {
                Some(a) => a.clone(),
                None    => T::zero()
            };
            let b: T = match other.coef.get(i) {
                Some(b) => b.clone(),
                None    => T::zero()
            };
            poly.coef.push( a + b );
        }

        concat( &mut poly );
       
        poly
    }
}


// Sub トレイトの実装
impl<T> Sub for Polynomial<T> where T: Clone + Zero + Sub<Output=T> + PartialEq
{
    type Output = Polynomial<T>;

    fn sub(self, other: Self) -> Self {
        let mut poly = Polynomial { coef: Vec::new() };

        for i in 0..self.coef.len().max(other.coef.len()) {
            let a: T = match self.coef.get(i) {
                Some(a) => a.clone(),
                None    => T::zero()
            };
            let b: T = match other.coef.get(i) {
                Some(b) => b.clone(),
                None    => T::zero()
            };
            poly.coef.push( a - b );
        }

        concat( &mut poly );
        poly
    }
}


// Neg トレイトの実装
impl<T> Neg for Polynomial<T> where T: Clone + Zero + Neg<Output=T> + PartialEq
{
    type Output = Polynomial<T>;

    fn neg(self) -> Polynomial<T> {
        Polynomial {
            coef: self.coef.into_iter().map(|c| c.neg() ).collect()
        }
    }
}



// Mul トレイトの実装 (スカラー倍)
impl<T> Mul<T> for Polynomial<T> where T: Clone + Zero + Mul<Output=T> + PartialEq
{
    type Output = Polynomial<T>;

    fn mul(self, scalar: T) -> Polynomial<T> {
        let mut p = Polynomial {
            coef: self.coef.into_iter().map(|c| c * scalar.clone()).collect()
        };
        concat( &mut p );
        p
    }
}
// Mul トレイトの実装 (多項式の積)
impl<T> Mul for Polynomial<T> where T: Clone + Zero + Add + Mul<Output=T> + PartialEq
{
    type Output = Polynomial<T>;

    fn mul(self, other: Polynomial<T>) -> Polynomial<T> {
        // 新しく多項式インスタンスを生成
        let mut poly = Polynomial { coef: Vec::new() };
        // 係数を順番に計算
        for k in 0..self.coef.len()+other.coef.len()-1 {
            let mut tot = T::zero();
            for n in (k+1).saturating_sub(other.coef.len())..self.coef.len().min(k+1) {
                tot = tot + self.coef[n].clone() * other.coef[k-n].clone();
            }
            poly.coef.push( tot );
        }
        // 最高次の係数がノンゼロになることを保証
        concat( &mut poly );
        // 多項式を返す
        poly
    }
}




// 最高次の係数がゼロである場合に, その係数を除去して最高次の係数がノンゼロになるように調整する.
// すべての係数がゼロであるならばゼロ多項式にする.
fn concat<T>( p: &mut Polynomial<T> ) where T: Zero + PartialEq {
    loop {
        match p.coef.last() {
            Some(c) => {
                if *c != T::zero() { break; } else {
                    p.coef.pop().unwrap();
                }
            },
            None    => { p.coef.push(T::zero()); break; }
        }
    }
}


#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_new() {
        let p = Polynomial::new(&[ 1, 0, 0 ]);
        assert_eq!( p, Polynomial::new(&[1]) );
    }

    #[test]
    fn test_eval() {
        // スカラー値の代入
        let p = Polynomial::new(&[ 0, 2, 1 ]).unwrap();
        assert_eq!( p.eval(0), 0 );
        assert_eq!( p.eval(1), 3 );
        assert_eq!( p.eval(2), 8 );

        // 多項式の代入
        let q = Polynomial::new(&[ 3, 1 ]).unwrap();
        assert_eq!( p.eval(q), Polynomial::new(&[ 15, 8, 1 ]).unwrap() );
    }
   
    #[test]
    fn test_add() {
        // スカラー値を足す
        let p = Polynomial::new(&[ 0, 1 ]).unwrap();
        assert_eq!( p + 1, Polynomial::new(&[ 1, 1 ]).unwrap() );

        // 多項式を足す
        let p1 = Polynomial::new(&[ 1., 2. ]).unwrap();
        let p2 = Polynomial::new(&[ 1., 2., 3. ]).unwrap();

        assert_eq!( p1 + p2, Polynomial::new(&[ 2., 4., 3. ]).unwrap() );
    }

    #[test]
    fn test_sub() {
        let p1 = Polynomial::new(&[ 1, 2 ]).unwrap();
        let p2 = Polynomial::new(&[ 1, 2, 3 ]).unwrap();

        assert_eq!( p1 - p2, Polynomial::new(&[ 0, 0, -3 ]).unwrap() );
    }

    #[test]
    fn test_neg() {
        let p1 = Polynomial::new(&[ 0, 0, 4 ]).unwrap();
        let p2 = - p1.clone();
        assert!( ( p1 + p2 ).is_zero() );
    }

    #[test]
    fn test_mul() {
        // スカラーの積
        let p: Polynomial<i32> = Polynomial::new(&[ 1, 2 ]).unwrap();
        let scalar: i32 = 3;
        assert_eq!( p * scalar, Polynomial::new(&[ 3, 6 ]).unwrap() );

        // 多項式の積
        let p1 = Polynomial::new(&[ 1, 2 ]).unwrap();
        let p2 = Polynomial::new(&[ 1, 2, 3 ]).unwrap();
        assert_eq!( p1 * p2, Polynomial::new(&[ 1, 4, 7, 6 ]).unwrap() );
    }
}
ラベル:

コメント

コメントを投稿

このブログの人気の投稿

matplotlib.histのnormedが変

以下の内容はPython 3.5.2 + matplotlib 1.5.1およびPython 3.6.2 + matplotlib 2.0.2で検証した. 作業日は2017年10月22日, 前者はUbuntu 16.04 on Win10 (WSL), 後者はDebian 8.9 (Anaconda). 規格化したいのにできない なにか数値の列 data があったとして, そのヒストグラムをmatplotlibでプロットしたいとする. 普通に plt.hist( data ) とすると, これは縦軸が各bin内に入るデータ点が何個あるかを表すことになる. これをデータ総数 len(data) で規格化したプロットにしようと思って plt.hist( data, normed=True ) または normed=1 とかやっちゃうと, 思った通りのアウトプットにならずに頭を傾げることになる. 例えば: import numpy as np import matplotlib.pyplot as plt data = np.random.normal(0,0.1,1000) weights = np.ones(len(data))/len(data) plt.hist( data, weights=weights ) plt.show() アウトプットは で, 縦軸が1を超えるとか, 意味がわからない. 原因 matplotlibのドキュメント を見ても何も書いてない. これは numpyのドキュメント に答えが書いてあるからで, 要するに normed オプションは事実上 density オプションと等しく, これは縦軸を確率分布関数と思って規格化するオプションである, と. 従って, normed=True オプションを指定すると, binの 面積 が1に規格化されることになる. いま欲しいものは値の 総和 が1に規格化されたアウトプットなのだから, binの幅が1でない限り, 欲しい結果は得られない. 対策 代わりに weights オプションを指定すればこの問題は解決できる. これは data の1つの値の重みを指定するパラメー
コメントを投稿
続きを読む

UbuntuのPDF編集ツールの使い方まとめ

PDFtkやpoppler-utilsに含まれるツールを使ってPDFを編集するコマンドのまとめです. 0. インストール sudo apt install pdftk sudo apt install poppler-utils UbuntuあるいはBash on Windowsならaptコマンドで一発. 1. PDFの文書情報の表示 pdfinfo (ファイル名) 2. PDFの分割と結合 pdftk (入力ファイル) cat (ページの指定) output (出力ファイル) ページ数の指定は1-12 14-endみたいな形で書けばよい. 入力ファイルを複数指定してページ数の指定を省略すればPDFの結合ができる. 3. PDFをJPEG/PNGに変換 pdftoppm [オプション] (入力ファイル名) (出力ファイル名の接頭辞) JPEGが欲しい場合-jpegを, PNGが欲しい場合-pngをオプションに指定する. デフォルトの解像度はDPI150で粗すぎるのでオプションで-r 300ないし-r 600を指定すべき. 白黒画像にしたい場合は-monoないし-grayを指定 (個人的には-grayのが好み). 複数ページのPDFを変換すると, 出力ファイルは (出力ファイル名の接頭辞)-1.jpg (出力ファイル名の接頭辞)-2.jpg みたいな形で生成される. 4. PDFからテキストを抽出 pdftotext (入力ファイル名) (出力ファイル名) 出力ファイルの文字コードはUTF-8, 改行コードはLF. 出力ファイル名は省略可能. 5. PDFから画像を抽出 pdfimages [オプション] (入力ファイル名) (出力ファイル名の接頭辞) オプションについて: -pngでPNG, -jでJPEG, -tiffでTIFF形式で出力される (オプションなしだとPPM形式) 6. PDFを回転 pdftk (入力ファイル名) cat 1-end(向き) output (出力ファイル名) 向きとしてはleft, right, downまたはnorth, east, west, southが使用可能. ※向きを指定する前にスペースを入れるとエラーになる. 関連ページ
コメントを投稿
続きを読む

源ノ明朝/源ノ角ゴシックをLuaLaTeXで使用する (Windows)

源ノ明朝 (Source Han Serif), 源ノ角ゴシック (Source Han Sans) はAdobeとGoogleが作成したCJKフォントで, オープンソースフォントとして公開されています (ライセンスはSIL Open Source License 1.1). 以下では, LuaTeXエンジンでこのフォントを (日本語で) 使用する方法を説明します. あるいは, otfファイルがダウンロードできれば, どのようなフォントにも適用可能です (PDFにフォント埋め込みすることに関してライセンスに注意が必要ですが). Windows10を前提としますが, 他の環境 (Linux/Mac) でもほぼ同じ方法で導入できるんじゃないでしょうか. 使用したのはTeX Live 2016 (LuaTeX-0.95) です.
コメントを投稿
続きを読む