SCIP in C++11 ― 2.3.2節その1

代数操作その1：記号微分と問題2.56～2.57

今後の構文解析の手始めになる、多項式の処理の実装。

結局ネタは構文木の構築なので、3章のデジタル回路が出てから4章以降の、

というかこの本の主要部となる、計算機のシミュレータ・Schemeインタープリタの実装、

といった事項と大きく関係してくる。

数式処理の本格的なのを自分で作るくらいならMaximaとか使う。

が、最終的には2.5.3節をちょっと応用すれば体上の有理関数までいけるわけで、

問題2.6を応用すればできそうな群の抽象表現とも組み合わせれば、

Galois理論の計算とかできちゃいそう。

まあ実際にやりはしないだろうけど、なんかワクワクするｗ

数と変数の識別はC++11の新機能regexを本当は使いたいが、

gcc-4.7ではまだ実装が不完全な模様。boostという手もなくはないが、

手っ取り早い苦肉の策でC++11のstodを使い、

invalid_argument例外をキャッチして判断する。

これはListTreeクラス内部に入れているので略。

その他、isNumberやpower,

またListTreeクラスのフレンド函数ととしてoperator+,*,==などを

実装しているが、これらも略。

問題2.57の出力では(* 3 (* 3...))なんてのが出てて、

これは　(* 9...)なんて風になって欲しくはあるが、

凝り出したら切りが無いのでそのまま。

----

//number?

const bool isNumber(const List& _listTree)

{return(_listTree->isNumber());}

//---------abstraction barrier---------

// make a polynomial

template<typename ... PolynomialExps>

const List makeExpression(const PolynomialExps ...elements)

{return(makeList(elements...));}

// print an expression

const string expString(const List& expression)

{return(listString(expression));}

//make number leaf

template <typename NumberType>

const List makeNumber(const NumberType num)

{return(makeLeaf(num));}

//make symbol leaf

template <typename SymbolType>

const List makeSymbol(const SymbolType symbol)

{return(makeLeaf(symbol));}

//equal?

template<typename LeafType>

const bool isEqual(const List& list1,

                   const List& list2)

{

    if(isNull(list1) && isNull(list2)){return(true);}

    else if(isNull(list1) || isNull(list2)){return(false);}

    else if(value<LeafType>(car(list1))

            !=value<LeafType>(car(list2)))

        {return(false);}

    return(isEqual<LeafType>(cdr(list1),cdr(list2)));

}

//symbol?

const bool isSymbol(const List& x)

{return(!isPair(x) && !isNumber(x));}

//=number?

const bool isEqNumber(const List& expression,

                      const List& number)

{return(isNumber(expression) && expression==number);}

template<typename LeafType>

const bool isEqNumber(const List& expression,

                      const LeafType& number)

{return(isEqNumber(expression,makeNumber(number)));}

//---------abstraction barrier---------

//variable?

const bool isVariable(const List& x)

{return(isSymbol(x));}

//same-variable?

const bool isSameVariable(const List& v1,const List& v2)

{return(isVariable(v1) && isVariable(v2) && isEq(v1,v2));}

//sum?

const bool isSum(const List& x)

{return(isPair(x) && isEq("+",car(x)));}

//addend

const List addend(const List& x)

{return(cadr(x));}

//augend

const List augend(const List& x)

{

    if(isNull(cdddr(x))){return(caddr(x));}

    return(cons("+",cddr(x)));

}

//make-sum

const List makeSum(const List& a1, const List& a2)

{

    if(isEqNumber(a1,0)){return(a2);}

    else if(isEqNumber(a2,0)){return(a1);}

    else if(isNumber(a1)&&isNumber(a2)){return(a1+a2);}

    /*else if(isNumber(a1)){

        return(makeExpression("+",a2,a1));

        }*/

    return(makeExpression("+",a1,a2));

}

//product?

const bool isProduct(const List& x)

{return(isPair(x) && isEq("*",car(x)));}

//multiplier

const List multiplier(const List& x)

{return(cadr(x));}

//multoplicand

const List multiplicand(const List& x)

{

    if(isNull(cdddr(x))){return(caddr(x));}

    return(cons("*",cddr(x)));

}

//make-product

const List makeProduct(const List& m1, const List& m2)

{

    if(isEqNumber(m1,0)||isEqNumber(m2,0)){return(makeNumber(0));}

    else if(isEqNumber(m1,1)){return(m2);}

    else if(isEqNumber(m2,1)){return(m1);}

    else if(isNumber(m1)&&isNumber(m2)){

        return(m1*m2);

    }else if(isNumber(m2)){

        return(makeExpression("*",m2,m1));

    }

    return(makeExpression("*",m1,m2));

}

//exponentiation?

const bool isExponentiation(const List& x)

{return(isPair(x) && isEq("**",car(x)));}

//base

const List base(const List& x)

{return(cadr(x));}

//exponent

const List exponent(const List& x)

{return(caddr(x));}

//make-sum

const List makeExponentiation(const List& b, const List& e)

{

    if(isEqNumber(b,0)){return(makeNumber(0));}

    else if(isEqNumber(e,0)){return(makeNumber(1));}

    else if(isEqNumber(e,1)){return(e);}

    else if(isNumber(e)&&isNumber(b)){

        return(power(b,e));

    }

    return(makeExpression("**",b,e));

}

//---------abstraction barrier---------

//deriv

const List derivative(const List& expression,

                      const List& variable)

{

    if(isNumber(expression)){return(makeNumber(0));}

    else if(isVariable(expression)){

        if(isSameVariable(expression,variable)){

            return(makeNumber(1));

        }else{

            return(makeNumber(0));

        }

    }else if(isSum(expression)){

        return(makeSum

               (derivative(addend(expression),variable),

                derivative(augend(expression),variable)));

    }else if(isProduct(expression)){

        return(makeSum

               (makeProduct

                (multiplier(expression),

                 derivative(multiplicand(expression),variable)),

                makeProduct

                (derivative(multiplier(expression),variable),

                 multiplicand(expression))));

    }else if(isExponentiation(expression)){

        return(makeProduct

               (exponent(expression),

                makeProduct

                (makeExponentiation

                 (base(expression),

                  exponent(expression)+makeNumber(-1)),

                 derivative(base(expression),variable))));

    }

   

    cerr<<"unknown exoression type --- DERIV "<<expression<<endl;

    return(makeExpression());

}

template<typename LeafType>

const List derivative(const List& expression,

                      const LeafType& variable)

{return(derivative(expression,makeSymbol(variable)));}

int main(int argc, char** argv)

{

    const string variable("x");

    const auto polynomial1(makeExpression("+","x",3));

    cout<<"(deriv "<<expString(polynomial1)

        <<" "<<variable<<")="

        <<expString(derivative(polynomial1,variable))<<endl;

   

    const auto polynomial2(makeExpression("*","x","y"));

    cout<<"(deriv "<<expString(polynomial2)

        <<" "<<variable<<")="

        <<expString(derivative(polynomial2,variable))<<endl;

   

    const auto polynomial3

        (makeExpression

         ("*",makeExpression("*","x","y"),

          makeExpression("+","x",3)));

    cout<<"(deriv "<<expString(polynomial3)

        <<" "<<variable<<")="

        <<expString(derivative(polynomial3,variable))<<endl;

    cout<<endl<<"Excersize 2.56:"<<endl;

    const auto polynomial4

        (makeExpression("**","x","4"));

    cout<<"(deriv "<<expString(polynomial4)

        <<" "<<variable<<")="

        <<expString(derivative(polynomial4,variable))<<endl;

    cout<<endl<<"Excersize 2.57:"<<endl;

    const auto polynomial5

        (makeExpression("+",makeExpression("**","x","4"),

                        makeExpression

                        ("*",3,makeExpression("**","x",3)),

                        makeExpression("*",3,"x")));

    cout<<"(deriv "<<expString(polynomial5)

        <<" "<<variable<<")="

        <<expString(derivative(polynomial5,variable))<<endl;

    return(0);

}

----

出力

----

(deriv (+ x 3) x)=1

(deriv (* x y) x)=y

(deriv (* (* x y) (+ x 3)) x)=(+ (* x y) (* y (+ x 3)))

Excersize 2.56:

(deriv (** x 4) x)=(* 4 (** x 3))

Excersize 2.57:

(deriv (+ (** x 4) (* 3 (** x 3)) (* 3 x)) x)=(+ (* 4 (** x 3)) (+ (* 3 (* 3 (** x 2))) 3))