this repo has no description
dekker.one
1/* -*- mode: C++; c-basic-offset: 2; indent-tabs-mode: nil -*- */
2/*
3 * Main authors:
4 * Vincent Barichard <Vincent.Barichard@univ-angers.fr>
5 *
6 * Copyright:
7 * Vincent Barichard, 2012
8 *
9 * This file is part of Gecode, the generic constraint
10 * development environment:
11 * http://www.gecode.org
12 *
13 * Permission is hereby granted, free of charge, to any person obtaining
14 * a copy of this software and associated documentation files (the
15 * "Software"), to deal in the Software without restriction, including
16 * without limitation the rights to use, copy, modify, merge, publish,
17 * distribute, sublicense, and/or sell copies of the Software, and to
18 * permit persons to whom the Software is furnished to do so, subject to
19 * the following conditions:
20 *
21 * The above copyright notice and this permission notice shall be
22 * included in all copies or substantial portions of the Software.
23 *
24 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
25 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
26 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
27 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
28 * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
29 * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
30 * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
31 *
32 */
33
34namespace Gecode { namespace Float { namespace Trigonometric {
35
36
37 /*
38 * ASin projection function
39 *
40 */
41template<class V>
42void aSinProject(Rounding& r, const V& aSinIv, FloatNum& iv_min, FloatNum& iv_max, int& n_min, int& n_max) {
43 #define I0_PI_2I FloatVal(0,pi_half_upper())
44 #define IPI_2_PII FloatVal(pi_half_lower(),pi_upper())
45 #define IPI_3PI_2I FloatVal(pi_lower(),3*pi_half_upper())
46 #define I3PI_2_2PII FloatVal(3*pi_half_lower(),pi_twice_upper())
47 #define POS(X) ((I0_PI_2I.in(X))?0: (IPI_2_PII.in(X))?1: (IPI_3PI_2I.in(X))?2: 3 )
48 #define ASININF_DOWN r.asin_down(aSinIv.min())
49 #define ASINSUP_UP r.asin_up(aSinIv.max())
50
51 // 0 <=> in [0;PI/2]
52 // 1 <=> in [PI/2;PI]
53 // 2 <=> in [PI;3*PI/2]
54 // 3 <=> in [3*PI/2;2*PI]
55 switch ( POS(iv_min) )
56 {
57 case 0:
58 if (r.sin_down(iv_min) > aSinIv.max()) { n_min++; iv_min = -ASINSUP_UP; }
59 else if (r.sin_up(iv_min) < aSinIv.min()) { iv_min = ASININF_DOWN; }
60 break;
61 case 1:
62 if (r.sin_down(iv_min) > aSinIv.max()) { n_min++; iv_min = -ASINSUP_UP; }
63 else if (r.sin_up(iv_min) < aSinIv.min()) { n_min+=2; iv_min = ASININF_DOWN; }
64 break;
65 case 2:
66 if (r.sin_down(iv_min) > aSinIv.max()) { n_min++; iv_min = -ASINSUP_UP; }
67 else if (r.sin_up(iv_min) < aSinIv.min()) { n_min+=2; iv_min = ASININF_DOWN; }
68 break;
69 case 3:
70 if (r.sin_down(iv_min) > aSinIv.max()) { n_min+=3; iv_min = -ASINSUP_UP; }
71 else if (r.sin_up(iv_min) < aSinIv.min()) { n_min+=2; iv_min = ASININF_DOWN; }
72 break;
73 default:
74 GECODE_NEVER;
75 break;
76 }
77
78 // 0 <=> in [0;PI/2]
79 // 1 <=> in [PI/2;PI]
80 // 2 <=> in [PI;3*PI/2]
81 // 3 <=> in [3*PI/2;2*PI]
82 switch ( POS(iv_max) )
83 {
84 case 0:
85 if (r.sin_down(iv_max) > aSinIv.max()) { iv_max = ASINSUP_UP; }
86 else if (r.sin_up(iv_max) < aSinIv.min()) { n_max--; iv_max = -ASININF_DOWN; }
87 break;
88 case 1:
89 if (r.sin_down(iv_max) > aSinIv.max()) { iv_max = ASINSUP_UP; }
90 else if (r.sin_up(iv_max) < aSinIv.min()) { n_max++; iv_max = -ASININF_DOWN; }
91 break;
92 case 2:
93 if (r.sin_down(iv_max) > aSinIv.max()) { iv_max = ASINSUP_UP; }
94 else if (r.sin_up(iv_max) < aSinIv.min()) { n_max++; iv_max = -ASININF_DOWN; }
95 break;
96 case 3:
97 if (r.sin_down(iv_max) > aSinIv.max()) { n_max+=2; iv_max = ASINSUP_UP; }
98 else if (r.sin_up(iv_max) < aSinIv.min()) { n_max++; iv_max = -ASININF_DOWN; }
99 break;
100 default:
101 GECODE_NEVER;
102 break;
103 }
104 #undef ASININF_DOWN
105 #undef ASINSUP_UP
106 #undef POS
107 #undef I0_PI_2I
108 #undef IPI_2_PII
109 #undef IPI_3PI_2I
110 #undef I3PI_2_2PII
111}
112
113/*
114 * Bounds consistent sinus operator
115 *
116 */
117
118 template<class A, class B>
119 ExecStatus
120 Sin<A,B>::dopropagate(Space& home, A x0, B x1) {
121 GECODE_ME_CHECK(x1.eq(home,sin(x0.val())));
122 Rounding r;
123 int n_min = 2*static_cast<int>(r.div_up(x0.min(), pi_twice_upper()));
124 int n_max = 2*static_cast<int>(r.div_up(x0.max(), pi_twice_upper()));
125 if (x0.min() < 0) n_min-=2;
126 if (x0.max() < 0) n_max-=2;
127 FloatNum iv_min = r.sub_down(x0.min(),r.mul_down(n_min, pi_upper()));
128 FloatNum iv_max = r.sub_up (x0.max(),r.mul_down(n_max, pi_upper()));
129 aSinProject(r,x1,iv_min,iv_max,n_min,n_max);
130 FloatNum n_iv_min = r.add_down(iv_min,r.mul_down(n_min, pi_upper()));
131 FloatNum n_iv_max = r.add_up (iv_max,r.mul_down(n_max, pi_upper()));
132 if (n_iv_min > n_iv_max) return ES_FAILED;
133 GECODE_ME_CHECK(x0.eq(home,FloatVal(n_iv_min,n_iv_max)));
134 GECODE_ME_CHECK(x1.eq(home,sin(x0.val()))); // Redo sin because with x0 reduction, sin may be more accurate
135 return ES_OK;
136 }
137
138 template<class A, class B>
139 forceinline
140 Sin<A,B>::Sin(Home home, A x0, B x1)
141 : MixBinaryPropagator<A,PC_FLOAT_BND,B,PC_FLOAT_BND>(home,x0,x1) {}
142
143 template<class A, class B>
144 ExecStatus
145 Sin<A,B>::post(Home home, A x0, B x1) {
146 if (x0 == x1) {
147 GECODE_ME_CHECK(x0.eq(home,0.0));
148 } else {
149 GECODE_ME_CHECK(x1.gq(home,-1.0));
150 GECODE_ME_CHECK(x1.lq(home,1.0));
151 GECODE_ES_CHECK(dopropagate(home,x0,x1));
152 (void) new (home) Sin<A,B>(home,x0,x1);
153 }
154
155 return ES_OK;
156 }
157
158
159 template<class A, class B>
160 forceinline
161 Sin<A,B>::Sin(Space& home, Sin<A,B>& p)
162 : MixBinaryPropagator<A,PC_FLOAT_BND,B,PC_FLOAT_BND>(home,p) {}
163
164 template<class A, class B>
165 Actor*
166 Sin<A,B>::copy(Space& home) {
167 return new (home) Sin<A,B>(home,*this);
168 }
169
170 template<class A, class B>
171 ExecStatus
172 Sin<A,B>::propagate(Space& home, const ModEventDelta&) {
173 GECODE_ES_CHECK(dopropagate(home,x0,x1));
174 return (x0.assigned()) ? home.ES_SUBSUMED(*this) : ES_FIX;
175 }
176
177 /*
178 * Bounds consistent cosinus operator
179 *
180 */
181
182 template<class A, class B>
183 ExecStatus
184 Cos<A,B>::dopropagate(Space& home, A x0, B x1) {
185 GECODE_ME_CHECK(x1.eq(home,cos(x0.val())));
186 Rounding r;
187 FloatVal x0Trans = x0.val() + FloatVal::pi_half();
188 int n_min = 2*static_cast<int>(r.div_up(x0Trans.min(), pi_twice_upper()));
189 int n_max = 2*static_cast<int>(r.div_up(x0Trans.max(), pi_twice_upper()));
190 if (x0Trans.min() < 0) n_min-=2;
191 if (x0Trans.max() < 0) n_max-=2;
192 FloatNum iv_min = r.sub_down(x0Trans.min(),r.mul_down(n_min, pi_upper()));
193 FloatNum iv_max = r.sub_up (x0Trans.max(),r.mul_down(n_max, pi_upper()));
194 aSinProject(r,x1,iv_min,iv_max,n_min,n_max);
195 FloatNum n_iv_min = r.add_down(iv_min,r.mul_down(n_min, pi_upper()));
196 FloatNum n_iv_max = r.add_up (iv_max,r.mul_down(n_max, pi_upper()));
197 if (n_iv_min > n_iv_max) return ES_FAILED;
198 GECODE_ME_CHECK(x0.eq(home,FloatVal(n_iv_min,n_iv_max) - FloatVal::pi_half()));
199 GECODE_ME_CHECK(x1.eq(home,cos(x0.val()))); // Redo sin because with x0 reduction, sin may be more accurate
200 return ES_OK;
201 }
202
203 template<class A, class B>
204 forceinline
205 Cos<A,B>::Cos(Home home, A x0, B x1)
206 : MixBinaryPropagator<A,PC_FLOAT_BND,B,PC_FLOAT_BND>(home,x0,x1) {}
207
208 template<class A, class B>
209 ExecStatus
210 Cos<A,B>::post(Home home, A x0, B x1) {
211 if (x0 == x1) {
212 GECODE_ME_CHECK(x0.gq(home,0.7390851332151));
213 GECODE_ME_CHECK(x0.lq(home,0.7390851332152));
214 bool mod;
215 do {
216 mod = false;
217 GECODE_ME_CHECK_MODIFIED(mod,x0.eq(home,cos(x0.val())));
218 } while (mod);
219 } else {
220 GECODE_ME_CHECK(x1.gq(home,-1.0));
221 GECODE_ME_CHECK(x1.lq(home,1.0));
222 GECODE_ES_CHECK(dopropagate(home,x0,x1));
223 (void) new (home) Cos<A,B>(home,x0,x1);
224 }
225 return ES_OK;
226 }
227
228
229 template<class A, class B>
230 forceinline
231 Cos<A,B>::Cos(Space& home, Cos<A,B>& p)
232 : MixBinaryPropagator<A,PC_FLOAT_BND,B,PC_FLOAT_BND>(home,p) {}
233
234 template<class A, class B>
235 Actor*
236 Cos<A,B>::copy(Space& home) {
237 return new (home) Cos<A,B>(home,*this);
238 }
239
240 template<class A, class B>
241 ExecStatus
242 Cos<A,B>::propagate(Space& home, const ModEventDelta&) {
243 GECODE_ES_CHECK(dopropagate(home,x0,x1));
244 return (x0.assigned()) ? home.ES_SUBSUMED(*this) : ES_FIX;
245 }
246
247}}}
248
249// STATISTICS: float-prop
250