sincos.hpp

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