SCIP in C++11 ― 2.1.4節

区間算術演算：問題2.7～2.16

誤差の伝搬とかの統計処理に使えそうな一節。

まあ書いていてあるままではちょっと難ありだけど、応用は効きそうな。

なおエラーでも処理は止めないように、エラーの場合はNaNを入れることにする。

問題2.13は線型近似すりゃいいだけだろうし、

問題2.14～2.16は誤差の伝播とかちゃんと考えなかったらそりゃそうだろう。

なので問題2.13～2.16はあまり興味を惹かない、てか特に実装することもないのでパス

問題2.7～12のコードのみ示す。

----

typedef double DataType;

typedef List Interval;

const Interval makeInterval

(const DataType& lower,const DataType& upper){

    return(cons(min(lower,upper),max(lower,upper)));

}

const DataType lowerBound(const Interval& interval)

{return car(interval)->getDouble();}

const DataType upperBound(const Interval& interval)

{return cadr(interval)->getDouble();}

const Interval unDefinedInterval(void){

    return(makeInterval(numeric_limits<DataType>::quiet_NaN(),

                        numeric_limits<DataType>::quiet_NaN()));

}

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

const DataType width(const Interval& x){

    return((upperBound(x)-lowerBound(x))/2);

}

const Interval addInterval(const Interval& x, const Interval& y)

{

    return(makeInterval(lowerBound(x)+lowerBound(y),

                        upperBound(x)+upperBound(y)));

}

const Interval subInterval(const Interval& x, const Interval& y)

{

    return(makeInterval(lowerBound(x)-lowerBound(y),

                        upperBound(x)-upperBound(y)));

}

const Interval mulInterval(const Interval& x, const Interval& y)

{

    const DataType lowerX(lowerBound(x));

    const DataType upperX(upperBound(x));

    const DataType lowerY(lowerBound(y));

    const DataType upperY(upperBound(y));

   

    if(lowerX>0){

        if(lowerY>0){

            return(makeInterval(lowerX*lowerY,upperX*upperY));

        }else if(upperY<0){

            return(makeInterval(upperX*lowerY,lowerX*upperY));

        }else{

            return(makeInterval(upperX*lowerY,lowerX*upperY));

        }

    }else if(upperX<0){

        if(lowerY>0){

            return(makeInterval(lowerX*upperY,upperX*lowerY));

        }else if(upperY<0){

            return(makeInterval(upperX*upperY,lowerX*lowerY));

        }else{

            return(makeInterval(lowerX*upperY,lowerX*lowerY));

        }

    }else{

        if(lowerY>0){

            return(makeInterval(lowerX*upperY,upperX*upperY));

        }else if(upperY<0){

            return(makeInterval(upperX*lowerY,lowerX*lowerY));

        }else{

            return(makeInterval(min(lowerX*upperY,upperX*lowerY),

                                max(lowerX*lowerY,upperX*upperY)));

        }

    }

}

const Interval divInterval(const Interval& x, const Interval& y)

{

    if(0==upperBound(y) || 0==lowerBound(y)

       ||(upperBound(y)>0&&lowerBound(y)<0)){

        return(unDefinedInterval());

    }

    return(makeInterval(1.0/upperBound(y),1.0/lowerBound(y)));

}

const DataType centre(const Interval& i)

{

    return((lowerBound(i)+(upperBound(i))/2.0));

}

const Interval makeCenterPercent

(const DataType& centerValue, const double percentage)

{

    return(makeInterval(centerValue*(1.0-percentage/100.0),

                        centerValue*(1.0+percentage/100.0)));

}

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

const double percent(const Interval& i)

{

    return(width(i)/centre(i)*100.0);

}

void printInterval(const Interval& x)

{

    cout<<"["<<lowerBound(x)<<","<<upperBound(x)<<"]";

}

int main(int argc, char** argv)

{

    const Interval a(makeInterval(1,2));

    const Interval b(makeInterval(3,4));

    const Interval c(makeInterval(-2,3));

    const Interval d(makeInterval(0,1));

    cout<<"Excersize 2.7-2.11:"<<endl;

    printInterval(a);cout<<"+";printInterval(b);cout<<"=";

    const Interval a_Plus_b(addInterval(a,b));

    printInterval(a_Plus_b);

    cout<<"  (Width: "<<width(a)<<"+"<<width(b)

        <<"="<<width(a_Plus_b)<<")"<<endl;

    cout<<endl;

    printInterval(a);cout<<"-";printInterval(b);cout<<"=";

    const Interval a_Minus_b(subInterval(a,b));

    printInterval(a_Minus_b);

    cout<<"  (Width: "<<width(a)<<"-"<<width(b)

        <<"="<<width(a_Minus_b)<<")"<<endl;

    cout<<endl;

    printInterval(a);cout<<"*";printInterval(b);cout<<"=";

    const Interval a_Times_b(mulInterval(a,b));

    printInterval(a_Times_b);

    cout<<"  (Width: "<<width(a)<<"*"<<width(b)

        <<"!="<<width(a_Times_b)<<")"<<endl;

    cout<<endl;

   

    printInterval(a);cout<<"/";printInterval(b);cout<<"=";

    const Interval a_DividedBy_b(divInterval(a,b));

    printInterval(a_DividedBy_b);

    cout<<"  (Width: "<<width(a)<<"/"<<width(b)

        <<"!="<<width(a_DividedBy_b)<<")"<<endl;

    cout<<endl;

    printInterval(a);cout<<"/";printInterval(c);cout<<"=";

    printInterval(divInterval(a,c));

    cout<<endl;

    printInterval(a);cout<<"/";printInterval(d);cout<<"=";

    printInterval(divInterval(a,d));

    cout<<endl;

    printInterval(a);cout<<"*";printInterval(c);cout<<"=";

    printInterval(mulInterval(a,c));

    cout<<endl;

    return(0);

}

----

出力

----

Excersize 2.7-2.11:

[1,2]+[3,4]=[4,6]  (Width: 0.5+0.5=1)

[1,2]-[3,4]=[-2,-2]  (Width: 0.5-0.5=0)

[1,2]*[3,4]=[3,8]  (Width: 0.5*0.5!=2.5)

[1,2]/[3,4]=[0.25,0.333333]  (Width: 0.5/0.5!=0.0416667)

[1,2]/[-2,3]=[nan,nan]

[1,2]/[0,1]=[nan,nan]

[1,2]*[-2,3]=[-4,3]