SCIP in C++11 ― 2.5.3節その4

代数操作その7：有理関数（完全版）と問題2.96～2.97

汎用の冪乗手続きには、1.2節の高速ルーチンを使用。

reduceの結果で、分母、分子の各第1項の計数が両方共負だと見づらいので、

その時には分母・分子に-1を掛けるようにした。

----

namespace Generic{

（略）

const bool isReal(const List& x)

    {return(typeID(typeTag(x))<3);}

    const bool isNegative(const List& x)

    {

        if(!isReal(x)){return(false);}

        return(makeLeaf(0)!=applyGeneric("negative?",x));

    }

（略）

const List greatestCommonDivisor(const List& x,const List& y)

    {return(applyGeneric("gcd",x,y));}

    const List greatestCommonDivisorList(const List& numberList)

    {

        function<List(List,List)> gcdProc(greatestCommonDivisor);

        return(accumulate(gcdProc,makeNumber(0),numberList));

    }

    template<typename...ListTypes>

    const List greatestCommonDivisor

    (const List& x1, const List& x2, const ListTypes ... args)

    {

        return(greatestCommonDivisor(makeList(x1,x2,args...)));

    }

   

    const List reduce(const List& n,const List& d)

    {return(applyGeneric("reduce",n,d));}

    const List power(const List base,const int exponent)

    {

        if(isEqu(drop(base),makeNumber(1))||isZero(base)){return(base);}

        if(exponent<=0){

            return(makeNumber(1));

        }else if(0==exponent%2){

        return(square(power(base, exponent/2)));

        }

        return(mul(base,power(base,exponent-1)));

    }

（略）

}

class NumberArithmetic{

public:

NumberArithmetic(const string tagIn="number"):tagString(tagIn){

（略）

put(makeList("reduce",makeList(this->getTag(),this->getTag())),

                makeLeaf(function<List(List,List)>

                         ([this](const List& x,const List& y)

                          {

                              const auto reduceResult(this->reduce(x,y));

                              return(makeList(this->tag(car(reduceResult)),

                                              this->tag(cadr(reduceResult))));})));

put(makeList("negative?",makeList(this->getTag())),

            makeLeaf(function<List(List)>

                     ([this](const List& x)

                      {return(makeLeaf(this->isNegative(x)));})));

（略）

    }

const bool isNegative(const List& x)const

    {return(x<makeLeaf(0));}

    const List negate(const List& x)const

    {return(x*makeLeaf(-1));}

    const List reduce(const List& n,const List& d)const

    {

        const auto g(this->gcd(this->abs(n),this->abs(d)));

        return(makeList(this->div(n,g),this->div(d,g)));

    }

（略）

};

class RationalArithmetic{

public:

（略）

   const List makeRational

    (const List& numerator,const List& denominator)const

    {

        if(typeTag(numerator)!="polynomial"

           || typeTag(denominator)!="polynomial"){

            if(typeTag(numerator)=="number"

               && typeTag(denominator)=="number"){

                if(contents(denominator)<makeLeaf(0)){

                    return(Generic::reduce(Generic::negate(numerator),

                                           Generic::negate(denominator)));

                }else{

                    return(Generic::reduce(numerator,denominator));

                }

            }else{

                return(makeList(numerator,denominator));

            }

        }

        //polynomial

        const auto reduceResult(Generic::reduce(numerator,denominator));

        if(Generic::isNegative

           (Generic::leadingCoefficient(car(reduceResult)))

           && Generic::isNegative

           (Generic::leadingCoefficient(cadr(reduceResult))))

            {return(Generic::reduce

                    (Generic::negate(numerator),Generic::negate(denominator)));}

        return(Generic::reduce(numerator,denominator));

    }

（略）

};

class SparseTermList{

public:

    SparseTermList(const TagType tagIn="sparse"):tagString(tagIn)

    {

    （略）

        put(makeList("reduce",makeList(this->getTag(),this->getTag())),

            makeLeaf(function<List(List,List)>

                     ([this](const List& x,const List& y)

                      {

                          const List reduceResult(this->reduceTerms(x,y));

                          return(makeList(this->tag(car(reduceResult)),

                                          this->tag(cadr(reduceResult))));})));

（略）

}

    const List pseudoRemainderTerms(const List& L1, const List& L2)

    {

     const int factorPower(1+value<int>(this->termListOrder(L1))

                           -value<int>(this->termListOrder(L2)));

     if(factorPower<=0){

         return(this->remainderTerms(L1,L2));

     }

     return(this->mulTermsByAllTerms

            (this->makeTerm

             (makeLeaf(0),

              (Generic::power(this->firstCoefficient(L2),factorPower))),

             this->remainderTerms(L1,L2)));

    }

    const List gcdTerms(const List& L1, const List& L2)

    {

     if(this->isZero(L2)){return(L1);}

     return(this->gcdTerms

            (L2,this->pseudoRemainderTerms(L1,L2)));

    }

    const List reduceTerms(const List& n, const List& d)

    {

     const auto gcdPseudo(this->gcdTerms(n,d));

     const int factorPower

         (1+max(value<int>(this->termListOrder(n)),

                value<int>(this->termListOrder(d)))

          -value<int>(this->termListOrder(gcdPseudo)));

     const auto numberizeFactorTerm

         (this->makeTerm

          (makeLeaf(0),

           Generic::power

           (this->firstCoefficient(gcdPseudo),factorPower)));

     const auto nQuotient

         (car(this->divTerms

              (this->mulTermsByAllTerms(numberizeFactorTerm,n),

               gcdPseudo)));

     const auto dQuotient

         (car(this->divTerms

              (this->mulTermsByAllTerms(numberizeFactorTerm,d),

               gcdPseudo)));

    

     const auto commonTermFactorProc

         =[this](const List& termList){

         const auto commonTermFactor

         (Generic::makeRational

          (Generic::makeNumber(1),

           Generic::greatestCommonDivisorList

           (mapping

            (function<List(List)>

             ([](const List& term){

                 return(Generic::abs(cadr(term)));}),

             termList))));

         return(this->mulTermsByAllTerms

                (this->makeTerm(makeLeaf(0),commonTermFactor),termList));

     };

     return(makeList(commonTermFactorProc(nQuotient),

                     commonTermFactorProc(dQuotient)));

    }

(

（略）

};

class DenseTermList{

public:

    DenseTermList(const TagType tagIn="dense"):tagString(tagIn)

    {

（略）

put(makeList("reduce",makeList(this->getTag(),this->getTag())),

         makeLeaf(function<List(List,List)>

                  ([this](const List& x,const List& y)

                   {

                       const List reduceResult(this->reduceTerms(x,y));

                       return(makeList(this->tag(car(reduceResult)),

                                       this->tag(cadr(reduceResult))));})));

（略）

}

const List pseudoRemainderTerms(const List& L1, const List& L2)

    {

     const int factorPower(1+value<int>(this->termListOrder(L1))

                           -value<int>(this->termListOrder(L2)));

     if(factorPower<=0){

         return(this->remainderTerms(L1,L2));

     }

     return(this->mulTermsByAllTerms

            (Generic::power(this->firstCoefficient(L2),factorPower),

             makeLeaf(0),

             this->remainderTerms(L1,L2)));

    }

    const List gcdTerms(const List& L1, const List& L2)

    {

     if(this->isZero(L2)){

         return(L1);

     }else if(this->isConstant(L1)||this->isConstant(L2)){

         return(makeList(Generic::makeNumber(1)));

     }

     return(this->gcdTerms(L2,this->pseudoRemainderTerms(L1,L2)));

    }

   

    const List reduceTerms(const List& n, const List& d)

    {

     const auto gcdPseudo(this->gcdTerms(n,d));

     const int factorPower

         (1+max(value<int>(this->termListOrder(n)),

                value<int>(this->termListOrder(d)))

          -value<int>(this->termListOrder(gcdPseudo)));

     const auto numberizeFactorTerm

         (Generic::power

          (this->firstCoefficient(gcdPseudo),factorPower));

     const auto nQuotient

         (car(this->divTerms

              (this->mulTermsByAllTerms

               (numberizeFactorTerm,makeLeaf(0),n),

               gcdPseudo)));

     const auto dQuotient

         (car(this->divTerms

              (this->mulTermsByAllTerms

               (numberizeFactorTerm,makeLeaf(0),d),

               gcdPseudo)));

    

     const auto commonTermFactorProc

         =[this](const List& termList){

         const auto commonTermFactor

         (Generic::makeRational

          (Generic::makeNumber(1),

           Generic::greatestCommonDivisorList(termList)));

         return(this->mulTermsByAllTerms

                (commonTermFactor,makeLeaf(0),termList));

     };

     return(makeList(commonTermFactorProc(nQuotient),

                     commonTermFactorProc(dQuotient)));

    }

（略）

};

class PolynomialArithmetic{

public:

    PolynomialArithmetic(const TagType tagIn="polynomial")

     :tagString(tagIn),

      _sparseTermList(new SparseTermList()),

      _denseTermList(new DenseTermList())

    {

（略）

put(makeList("reduce",makeList(this->getTag(),this->getTag())),

         makeLeaf(function<List(List,List)>

                  ([this](const List& x,const List& y)

                   {

                       const auto reduceResult(this->reducePolynomial(x,y));

                       return(makeList(drop(this->tag(car(reduceResult))),

                                       drop(this->tag(cadr(reduceResult)))));})));

（略）

}

const List reducePolynomial(const List& p1,const List& p2)

    {

     if(this->isSameVariable

        (this->variable(p1),this->variable(p2))){

         const auto reduceResult

             (Generic::reduce(this->termList(p1),this->termList(p2)));

         return(makeList

                (this->makePolynomial

                 (this->variable(p1),car(reduceResult)),

                 this->makePolynomial

                 (this->variable(p2),cadr(reduceResult))));

     }

     cerr<<"Polynomials not in same variable -- GCD-POLY "

         <<listString(p1)<<listString(p2)<<endl;

     exit(1);

     return(makeList());

    }

（略）

};

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

int main(int argc, char** argv)

{

    installNumberPackage();

    installRationalPackage();

    installRealPackage();

    installComplexPackage();

    installPolynomialPackage();

    installCoercion();

    using namespace Generic;

    cout<<"Excersize 2.96 & 2.97:"<<endl;

    const auto p1(makePolynomial

               ("x",

                makeSparseTermList

                (makeList

                 (makeList(1,makeNumber(1)),

                  makeList(0,makeNumber(1))))));

    const auto p2(makePolynomial

               ("x",

                makeSparseTermList

                (makeList

                 (makeList(3,makeNumber(1)),

                  makeList(0,makeNumber(-1))))));

    const auto p3(makePolynomial

               ("x",

                makeSparseTermList

                (makeList

                 (makeList(1,makeNumber(1))))));

    const auto p4(makePolynomial

               ("x",

                makeSparseTermList

                (makeList

                 (makeList(2,makeNumber(1)),

                  makeList(0,makeNumber(-1))))));

    const auto rf1(makeRational(p1,p2));

    const auto rf2(makeRational(p3,p4));

    cout<<"p1 = "<<expressionString(p1)<<endl;

    cout<<"p2 = "<<expressionString(p2)<<endl;

    cout<<"p3 = "<<expressionString(p3)<<endl;

    cout<<"p4 = "<<expressionString(p4)<<endl;

    cout<<"rf1 = "<<expressionString(rf1)<<endl;

    cout<<"rf2 = "<<expressionString(rf2)<<endl;

    cout<<"(add rf1 rf2) = "<<expressionString(add(rf1,rf2))<<endl;

    uninstallNumberPackage();

    uninstallRationalPackage();

    uninstallRealPackage();

    uninstallComplexPackage();

    uninstallPolynomialPackage();

    return(0);

}

----

出力

----

Excersize 2.96 & 2.97:

p1 = x+1

p2 = x^3-1

p3 = x

p4 = x^2-1

rf1 = (x+1)/(x^3-1)

rf2 = (x)/(x^2-1)

(add rf1 rf2) = (x^3+2x^2+3x+1)/(x^4+x^3-x-1)