SCIP in C++11 ― 2.4節その2

代数操作その3：型による微分演算の振り分けと問題2.73

複素数で苦労して、データ主導プログラミングをそれなりにしっかり実装した分、

問題2.73では本質的な苦労はあまりなかったのだが、

C++の型を整形するインターフェース部分を書くのがだるい・・・。

これはかなり後を引く問題になるかも。

問題2.73dはgetに対応してputの順番を入れ替えるだけだから略。

----

typedef List Algebraic;

typedef TagType Operator;

typedef double Field;

//variable?

const bool isVariable(const Algebraic& x)

{return(isSymbol(x));}

//same-variable?

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

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

const function<Operator(Algebraic)> algebraicOperator=typeTag;

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

class AlgebraicOperation{

public:

    AlgebraicOperation(void)=delete;

    AlgebraicOperation(const Operator operatorIn)

        :operatorString(operatorIn){}

    virtual ~AlgebraicOperation(void){};

    const Operator getOperator(void)const

    {return(this->operatorString);}

    virtual const Algebraic operate(const Algebraic& expression)const

    {return(attachTag(this->getOperator(),expression));}

private:

    const Operator operatorString;

};

template <typename FieldType>

class AlgebraicSum :public AlgebraicOperation{

public:

    AlgebraicSum(void):AlgebraicOperation("+"){

        put("derivative",this->getOperator(),

            makeLeaf(function<Algebraic(Algebraic,Algebraic)>

                     ([this](const Algebraic& expression,

                             const Algebraic& variable)

                      {return(this->derivativeSum

                              (expression,variable));})));

       

        put("make-sum",this->getOperator(),

            makeLeaf(function<Algebraic(Algebraic,Algebraic)>

                     ([this](const Algebraic& addend,

                             const Algebraic& augend)

                      {return(this->makeSum(addend,augend));})));

    }

   

    const Algebraic derivativeSum

    (const Algebraic& sumExpression, const Algebraic& variable)

        const{

        return(this->makeSum

               (derivative(this->addend(sumExpression),variable),

                derivative(this->augend(sumExpression),variable)));

    }

   

    const Algebraic makeSum

    (const Algebraic& addend, const Algebraic& augend)const

    {

        if(isEqNumber(addend,0)){return(augend);}

        else if(isEqNumber(augend,0)){return(addend);}

        else if(isNumber(addend)&&isNumber(augend))

            {return(addend+augend);}

        return(this->operate(makeList(addend,augend)));

    }

    const Algebraic addend(const Algebraic& sumExpression)const

    {return(car(sumExpression));}

    const Algebraic augend(const Algebraic& sumExpression)const

    {return(cadr(sumExpression));}

};

template <typename FieldType>

class AlgebraicProduct :public AlgebraicOperation{

public:

    AlgebraicProduct(void):AlgebraicOperation("*"){

        put("derivative",this->getOperator(),

            makeLeaf(function<Algebraic(Algebraic,Algebraic)>

                     ([this](const Algebraic& expression,

                             const Algebraic& variable)

                      {return(this->derivativeProduct

                              (expression,variable));})));

       

        put("make-product",this->getOperator(),

            makeLeaf(function<Algebraic(Algebraic,Algebraic)>

                     ([this](const Algebraic& multiplier,

                             const Algebraic& multiplicand)

                      {return(this->makeProduct

                              (multiplier,multiplicand));})));

    }

   

    const Algebraic derivativeProduct

    (const Algebraic& productExpression, const Algebraic& variable)

        const{

        return(makeSum

               (this->makeProduct

                (this->multiplier(productExpression),

                 derivative

                 (this->multiplicand(productExpression),variable)),

                this->makeProduct

                (derivative

                 (this->multiplier(productExpression),variable),

                 this->multiplicand(productExpression))));

    }

    const Algebraic makeProduct

    (const Algebraic& multiplier, const Algebraic& multiplicand)

        const{

        if(isEqNumber(multiplier,0)||isEqNumber(multiplicand,0))

            {return(makeNumber(0));}

        else if(isEqNumber(multiplier,1)){return(multiplicand);}

        else if(isEqNumber(multiplicand,1)){return(multiplier);}

        else if(isNumber(multiplier)&&isNumber(multiplicand)){

            return(multiplier*multiplicand);

        }else if(isNumber(multiplicand)){

            return(this->operate(makeList(multiplicand,multiplier)));

        }

        return(this->operate(makeList(multiplier,multiplicand)));

    }

    const Algebraic multiplier

    (const Algebraic& productExpression)const

    {return(car(productExpression));}

    const Algebraic multiplicand

    (const Algebraic& productExpression)const

    {return(cadr(productExpression));}

};

template <typename FieldType>

class AlgebraicExponentiation :public AlgebraicOperation{

public:

    AlgebraicExponentiation(void):AlgebraicOperation("**"){

        put("derivative",this->getOperator(),

            makeLeaf(function<Algebraic(Algebraic,Algebraic)>

                     ([this](const Algebraic& expression,

                             const Algebraic& variable)

                      {return(this->derivativeExponentiation

                              (expression,variable));})));

       

        put("make-exponentiation",this->getOperator(),

            makeLeaf(function<Algebraic(Algebraic,Algebraic)>

                     ([this](const Algebraic& base,

                             const Algebraic& exponent)

                      {return(this->makeExponentiation

                              (base,exponent));})));

    }

   

    const Algebraic derivativeExponentiation

    (const Algebraic& exponentiationExpression,

     const Algebraic& variable)

        const{

        return(makeProduct

               (this->exponent(exponentiationExpression),

               

                makeProduct

                (this->makeExponentiation

                 (this->base(exponentiationExpression),

                  this->exponent

                  (exponentiationExpression)+makeNumber(-1)),

                 derivative(this->base

                            (exponentiationExpression),variable))));

    }

   

    const Algebraic makeExponentiation

    (const Algebraic& base, const Algebraic& exponent)const

    {

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

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

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

        else if(isNumber(exponent)&&isNumber(base)){

            return(power(base,exponent));

        }

        return(this->operate(makeList(base,exponent)));

    }

    const Algebraic base

    (const Algebraic& exponentiationExpression)const

    {return(car(exponentiationExpression));}

   

    const Algebraic exponent

    (const Algebraic& exponentiationExpression)const

    {return(cadr(exponentiationExpression));}

};

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

AlgebraicOperation* _algebraicPackage1(nullptr);

AlgebraicOperation* _algebraicPackage2(nullptr);

AlgebraicOperation* _algebraicPackage3(nullptr);

const void installSumPackage(void){

    _algebraicPackage1=new AlgebraicSum<Field>();

}

const void installProductPackage(void){

    _algebraicPackage2=new AlgebraicProduct<Field>();

}

const void installExponentiationPackage(void){

    _algebraicPackage3=new AlgebraicExponentiation<Field>();

}

const void uninstallAlgebraicPackages(void){

    if(nullptr!=_algebraicPackage1) delete _algebraicPackage1;

    if(nullptr!=_algebraicPackage2) delete _algebraicPackage2;

    if(nullptr!=_algebraicPackage3) delete _algebraicPackage3;

}

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

// make-sum

const Algebraic makeSum

(const Algebraic& addend,const Algebraic& augend)

{

    const Algebraic sumLeaf(get("make-sum","+"));

    if(isFunction(sumLeaf)){

        const auto sumProcedure

            (executable<Algebraic,Algebraic,Algebraic>(sumLeaf));

        return(sumProcedure(addend,augend));

    }

    return(makeList());

}

template<typename AType, typename BType>

const Algebraic makeSum(const AType& addend, const BType& augend)

{return(makeSum(makeLeaf(addend),makeLeaf(augend)));}

//make-product

const Algebraic makeProduct

(const Algebraic& multiplier,const Algebraic& multiplicand)

{

    const Algebraic productLeaf(get("make-product","*"));

    if(isFunction(productLeaf)){

        const auto productProcedure

            (executable<Algebraic,Algebraic,Algebraic>(productLeaf));

        return(productProcedure(multiplier,multiplicand));

    }

    return(makeList());

}

template<typename AType, typename BType>

const Algebraic makeProduct

(const AType& multiplier, const BType& multiplicand)

{return(makeProduct(makeLeaf(multiplier),makeLeaf(multiplicand)));}

//make-exponentiation

const Algebraic makeExponentiation

(const Algebraic& base,const Algebraic& exponent)

{

    const Algebraic exponentiationLeaf(get("make-exponentiation","**"));

    if(isFunction(exponentiationLeaf)){

        const auto exponentiationProcedure

            (executable<Algebraic,Algebraic,Algebraic>

             (exponentiationLeaf));

        return(exponentiationProcedure(base,exponent));

    }

    return(makeList());

}

template<typename AType, typename BType>

const Algebraic makeExponentiation

(const AType& base, const BType& exponent)

{return(makeExponentiation(makeLeaf(base),makeLeaf(exponent)));}

//deriv

const Algebraic getOperator(const Algebraic& expression)

{return(car(expression));}

const Algebraic operands(const Algebraic& expression)

{return(cdr(expression));}

const Algebraic derivative(const Algebraic& expression,

                           const Algebraic& variable)

{

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

    else if(isVariable(expression)){

        if(isSameVariable(expression,variable)){

            return(makeNumber(1));

        }else{

            return(makeNumber(0));

        }

    }

    const Algebraic derivativeLeaf

        (get("derivative",value<Operator>(getOperator(expression))));

    if(isFunction(derivativeLeaf)){

        const auto derivativeProcedure

            (executable<Algebraic,Algebraic,Algebraic>(derivativeLeaf));

        return(derivativeProcedure(operands(expression),variable));

    }

    return(makeList());

}

template<typename LeafType>

const Algebraic derivative(const Algebraic& expression,

                      const LeafType& variable)

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

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

int main(int argc, char** argv)

{

    installSumPackage();

    installProductPackage();

    installExponentiationPackage();

    const string variable("x");

    const auto polynomial1(makeSum("x",3));

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

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

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

   

    const auto polynomial2(makeProduct("x","y"));

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

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

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

   

    const auto polynomial3

        (makeProduct

         (makeProduct("x","y"),

          makeSum("x",3)));

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

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

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

    const auto polynomial4

        (makeExponentiation("x","4"));

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

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

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

    const auto polynomial5

        (makeSum(makeExponentiation("x","4"),

                 makeProduct(3,

                             makeExponentiation("x",3))));

    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)))

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

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