Generated on Tue May 22 09:39:50 2018 for Gecode by doxygen 1.6.3

sincos.hpp

Go to the documentation of this file.
00001 /* -*- mode: C++; c-basic-offset: 2; indent-tabs-mode: nil -*- */
00002 /*
00003  *  Main authors:
00004  *     Vincent Barichard <Vincent.Barichard@univ-angers.fr>
00005  *
00006  *  Copyright:
00007  *     Vincent Barichard, 2012
00008  *
00009  *  This file is part of Gecode, the generic constraint
00010  *  development environment:
00011  *     http://www.gecode.org
00012  *
00013  *  Permission is hereby granted, free of charge, to any person obtaining
00014  *  a copy of this software and associated documentation files (the
00015  *  "Software"), to deal in the Software without restriction, including
00016  *  without limitation the rights to use, copy, modify, merge, publish,
00017  *  distribute, sublicense, and/or sell copies of the Software, and to
00018  *  permit persons to whom the Software is furnished to do so, subject to
00019  *  the following conditions:
00020  *
00021  *  The above copyright notice and this permission notice shall be
00022  *  included in all copies or substantial portions of the Software.
00023  *
00024  *  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
00025  *  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
00026  *  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
00027  *  NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
00028  *  LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
00029  *  OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
00030  *  WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
00031  *
00032  */
00033 
00034 namespace Gecode { namespace Float { namespace Trigonometric {
00035 
00036 
00037   /*
00038    * ASin projection function
00039    *
00040    */
00041 template<class V>
00042 void aSinProject(Rounding& r, const V& aSinIv, FloatNum& iv_min, FloatNum& iv_max, int& n_min, int& n_max) {
00043   #define I0__PI_2I    FloatVal(0,pi_half_upper())
00044   #define IPI_2__PII   FloatVal(pi_half_lower(),pi_upper())
00045   #define IPI__3PI_2I  FloatVal(pi_lower(),3*pi_half_upper())
00046   #define I3PI_2__2PII FloatVal(3*pi_half_lower(),pi_twice_upper())
00047   #define POS(X) ((I0__PI_2I.in(X))?0: (IPI_2__PII.in(X))?1: (IPI__3PI_2I.in(X))?2: 3 )
00048   #define ASININF_DOWN r.asin_down(aSinIv.min())
00049   #define ASINSUP_UP r.asin_up(aSinIv.max())
00050 
00051   // 0 <=> in [0;PI/2]
00052   // 1 <=> in [PI/2;PI]
00053   // 2 <=> in [PI;3*PI/2]
00054   // 3 <=> in [3*PI/2;2*PI]
00055   switch ( POS(iv_min) )
00056   {
00057     case 0:
00058       if (r.sin_down(iv_min) > aSinIv.max())    { n_min++; iv_min = -ASINSUP_UP;   }
00059       else if (r.sin_up(iv_min) < aSinIv.min()) {          iv_min = ASININF_DOWN;  }
00060     break;
00061     case 1:
00062       if (r.sin_down(iv_min) > aSinIv.max())    { n_min++;  iv_min = -ASINSUP_UP;   }
00063       else if (r.sin_up(iv_min) < aSinIv.min()) { n_min+=2; iv_min = ASININF_DOWN;  }
00064     break;
00065     case 2:
00066       if (r.sin_down(iv_min) > aSinIv.max())    { n_min++;  iv_min = -ASINSUP_UP;   }
00067       else if (r.sin_up(iv_min) < aSinIv.min()) { n_min+=2; iv_min = ASININF_DOWN; }
00068     break;
00069     case 3:
00070       if (r.sin_down(iv_min) > aSinIv.max())    { n_min+=3; iv_min = -ASINSUP_UP;    }
00071       else if (r.sin_up(iv_min) < aSinIv.min()) { n_min+=2; iv_min = ASININF_DOWN; }
00072     break;
00073     default:
00074       GECODE_NEVER;
00075     break;
00076   }
00077 
00078   // 0 <=> in [0;PI/2]
00079   // 1 <=> in [PI/2;PI]
00080   // 2 <=> in [PI;3*PI/2]
00081   // 3 <=> in [3*PI/2;2*PI]
00082   switch ( POS(iv_max) )
00083   {
00084     case 0:
00085       if (r.sin_down(iv_max) > aSinIv.max())    {           iv_max = ASINSUP_UP;    }
00086       else if (r.sin_up(iv_max) < aSinIv.min()) { n_max--;  iv_max = -ASININF_DOWN; }
00087     break;
00088     case 1:
00089       if (r.sin_down(iv_max) > aSinIv.max())    {          iv_max = ASINSUP_UP;    }
00090       else if (r.sin_up(iv_max) < aSinIv.min()) { n_max++; iv_max = -ASININF_DOWN; }
00091     break;
00092     case 2:
00093       if (r.sin_down(iv_max) > aSinIv.max())    {          iv_max = ASINSUP_UP;    }
00094       else if (r.sin_up(iv_max) < aSinIv.min()) { n_max++; iv_max = -ASININF_DOWN; }
00095     break;
00096     case 3:
00097       if (r.sin_down(iv_max) > aSinIv.max())    { n_max+=2; iv_max = ASINSUP_UP;    }
00098       else if (r.sin_up(iv_max) < aSinIv.min()) { n_max++;  iv_max = -ASININF_DOWN; }
00099     break;
00100     default:
00101       GECODE_NEVER;
00102     break;
00103   }
00104   #undef ASININF_DOWN
00105   #undef ASINSUP_UP
00106   #undef POS
00107   #undef I0__PI_2I
00108   #undef IPI_2__PII
00109   #undef IPI__3PI_2I
00110   #undef I3PI_2__2PII
00111 }
00112 
00113 /*
00114    * Bounds consistent sinus operator
00115    *
00116    */
00117 
00118   template<class A, class B>
00119   ExecStatus
00120   Sin<A,B>::dopropagate(Space& home, A x0, B x1) {
00121     GECODE_ME_CHECK(x1.eq(home,sin(x0.val())));
00122     Rounding r;
00123     int n_min = 2*static_cast<int>(r.div_up(x0.min(), pi_twice_upper()));
00124     int n_max = 2*static_cast<int>(r.div_up(x0.max(), pi_twice_upper()));
00125     if (x0.min() < 0) n_min-=2;
00126     if (x0.max() < 0) n_max-=2;
00127     FloatNum iv_min = r.sub_down(x0.min(),r.mul_down(n_min, pi_upper()));
00128     FloatNum iv_max = r.sub_up  (x0.max(),r.mul_down(n_max, pi_upper()));
00129     aSinProject(r,x1,iv_min,iv_max,n_min,n_max);
00130     FloatNum n_iv_min = r.add_down(iv_min,r.mul_down(n_min, pi_upper()));
00131     FloatNum n_iv_max = r.add_up  (iv_max,r.mul_down(n_max, pi_upper()));
00132     if (n_iv_min > n_iv_max) return ES_FAILED;
00133     GECODE_ME_CHECK(x0.eq(home,FloatVal(n_iv_min,n_iv_max)));
00134     GECODE_ME_CHECK(x1.eq(home,sin(x0.val()))); // Redo sin because with x0 reduction, sin may be more accurate
00135     return ES_OK;
00136   }
00137 
00138   template<class A, class B>
00139   forceinline
00140   Sin<A,B>::Sin(Home home, A x0, B x1)
00141     : MixBinaryPropagator<A,PC_FLOAT_BND,B,PC_FLOAT_BND>(home,x0,x1) {}
00142 
00143   template<class A, class B>
00144   ExecStatus
00145   Sin<A,B>::post(Home home, A x0, B x1) {
00146     if (same(x0,x1)) {
00147       GECODE_ME_CHECK(x0.eq(home,0.0));
00148     } else {
00149       GECODE_ME_CHECK(x1.gq(home,-1.0));
00150       GECODE_ME_CHECK(x1.lq(home,1.0));
00151       GECODE_ES_CHECK(dopropagate(home,x0,x1));
00152       (void) new (home) Sin<A,B>(home,x0,x1);
00153     }
00154 
00155     return ES_OK;
00156   }
00157 
00158 
00159   template<class A, class B>
00160   forceinline
00161   Sin<A,B>::Sin(Space& home, Sin<A,B>& p)
00162     : MixBinaryPropagator<A,PC_FLOAT_BND,B,PC_FLOAT_BND>(home,p) {}
00163 
00164   template<class A, class B>
00165   Actor*
00166   Sin<A,B>::copy(Space& home) {
00167     return new (home) Sin<A,B>(home,*this);
00168   }
00169 
00170   template<class A, class B>
00171   ExecStatus
00172   Sin<A,B>::propagate(Space& home, const ModEventDelta&) {
00173     GECODE_ES_CHECK(dopropagate(home,x0,x1));
00174     return (x0.assigned()) ? home.ES_SUBSUMED(*this) : ES_FIX;
00175   }
00176 
00177   /*
00178    * Bounds consistent cosinus operator
00179    *
00180    */
00181 
00182   template<class A, class B>
00183   ExecStatus
00184   Cos<A,B>::dopropagate(Space& home, A x0, B x1) {
00185     GECODE_ME_CHECK(x1.eq(home,cos(x0.val())));
00186     Rounding r;
00187     FloatVal x0Trans = x0.val() + FloatVal::pi_half();
00188     int n_min = 2*static_cast<int>(r.div_up(x0Trans.min(), pi_twice_upper()));
00189     int n_max = 2*static_cast<int>(r.div_up(x0Trans.max(), pi_twice_upper()));
00190     if (x0Trans.min() < 0) n_min-=2;
00191     if (x0Trans.max() < 0) n_max-=2;
00192     FloatNum iv_min = r.sub_down(x0Trans.min(),r.mul_down(n_min, pi_upper()));
00193     FloatNum iv_max = r.sub_up  (x0Trans.max(),r.mul_down(n_max, pi_upper()));
00194     aSinProject(r,x1,iv_min,iv_max,n_min,n_max);
00195     FloatNum n_iv_min = r.add_down(iv_min,r.mul_down(n_min, pi_upper()));
00196     FloatNum n_iv_max = r.add_up  (iv_max,r.mul_down(n_max, pi_upper()));
00197     if (n_iv_min > n_iv_max) return ES_FAILED;
00198     GECODE_ME_CHECK(x0.eq(home,FloatVal(n_iv_min,n_iv_max) - FloatVal::pi_half()));
00199     GECODE_ME_CHECK(x1.eq(home,cos(x0.val()))); // Redo sin because with x0 reduction, sin may be more accurate
00200     return ES_OK;
00201   }
00202 
00203   template<class A, class B>
00204   forceinline
00205   Cos<A,B>::Cos(Home home, A x0, B x1)
00206     : MixBinaryPropagator<A,PC_FLOAT_BND,B,PC_FLOAT_BND>(home,x0,x1) {}
00207 
00208   template<class A, class B>
00209   ExecStatus
00210   Cos<A,B>::post(Home home, A x0, B x1) {
00211     if (same(x0,x1)) {
00212       GECODE_ME_CHECK(x0.gq(home,0.7390851332151));
00213       GECODE_ME_CHECK(x0.lq(home,0.7390851332152));
00214       bool mod;
00215       do {
00216         mod = false;
00217         GECODE_ME_CHECK_MODIFIED(mod,x0.eq(home,cos(x0.val())));
00218       } while (mod);
00219     } else {
00220       GECODE_ME_CHECK(x1.gq(home,-1.0));
00221       GECODE_ME_CHECK(x1.lq(home,1.0));
00222       GECODE_ES_CHECK(dopropagate(home,x0,x1));
00223       (void) new (home) Cos<A,B>(home,x0,x1);
00224     }
00225     return ES_OK;
00226   }
00227 
00228 
00229   template<class A, class B>
00230   forceinline
00231   Cos<A,B>::Cos(Space& home, Cos<A,B>& p)
00232     : MixBinaryPropagator<A,PC_FLOAT_BND,B,PC_FLOAT_BND>(home,p) {}
00233 
00234   template<class A, class B>
00235   Actor*
00236   Cos<A,B>::copy(Space& home) {
00237     return new (home) Cos<A,B>(home,*this);
00238   }
00239 
00240   template<class A, class B>
00241   ExecStatus
00242   Cos<A,B>::propagate(Space& home, const ModEventDelta&) {
00243     GECODE_ES_CHECK(dopropagate(home,x0,x1));
00244     return (x0.assigned()) ? home.ES_SUBSUMED(*this) : ES_FIX;
00245   }
00246 
00247 }}}
00248 
00249 // STATISTICS: float-prop
00250