@@ -115,6 +115,26 @@ void aSinProject(Rounding& r, const V& aSinIv, FloatNum& iv_min, FloatNum& iv_ma
115115 *
116116 */
117117
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+
118138 template <class A , class B >
119139 forceinline
120140 Sin<A,B>::Sin(Home home, A x0, B x1)
@@ -128,6 +148,7 @@ void aSinProject(Rounding& r, const V& aSinIv, FloatNum& iv_min, FloatNum& iv_ma
128148 } else {
129149 GECODE_ME_CHECK (x1.gq (home,-1.0 ));
130150 GECODE_ME_CHECK (x1.lq (home,1.0 ));
151+ GECODE_ES_CHECK (dopropagate (home,x0,x1));
131152 (void ) new (home) Sin<A,B>(home,x0,x1);
132153 }
133154
@@ -149,20 +170,7 @@ void aSinProject(Rounding& r, const V& aSinIv, FloatNum& iv_min, FloatNum& iv_ma
149170 template <class A , class B >
150171 ExecStatus
151172 Sin<A,B>::propagate(Space& home, const ModEventDelta&) {
152- GECODE_ME_CHECK (x1.eq (home,sin (x0.val ())));
153- Rounding r;
154- int n_min = 2 *static_cast <int >(r.div_up (x0.min (), pi_twice_upper ()));
155- int n_max = 2 *static_cast <int >(r.div_up (x0.max (), pi_twice_upper ()));
156- if (x0.min () < 0 ) n_min-=2 ;
157- if (x0.max () < 0 ) n_max-=2 ;
158- FloatNum iv_min = r.sub_down (x0.min (),r.mul_down (n_min, pi_upper ()));
159- FloatNum iv_max = r.sub_up (x0.max (),r.mul_down (n_max, pi_upper ()));
160- aSinProject (r,x1,iv_min,iv_max,n_min,n_max);
161- FloatNum n_iv_min = r.add_down (iv_min,r.mul_down (n_min, pi_upper ()));
162- FloatNum n_iv_max = r.add_up (iv_max,r.mul_down (n_max, pi_upper ()));
163- if (n_iv_min > n_iv_max) return ES_FAILED;
164- GECODE_ME_CHECK (x0.eq (home,FloatVal (n_iv_min,n_iv_max)));
165- GECODE_ME_CHECK (x1.eq (home,sin (x0.val ()))); // Redo sin because with x0 reduction, sin may be more accurate
173+ GECODE_ES_CHECK (dopropagate (home,x0,x1));
166174 return (x0.assigned ()) ? home.ES_SUBSUMED (*this ) : ES_FIX;
167175 }
168176
@@ -171,6 +179,27 @@ void aSinProject(Rounding& r, const V& aSinIv, FloatNum& iv_min, FloatNum& iv_ma
171179 *
172180 */
173181
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+
174203 template <class A , class B >
175204 forceinline
176205 Cos<A,B>::Cos(Home home, A x0, B x1)
@@ -190,6 +219,7 @@ void aSinProject(Rounding& r, const V& aSinIv, FloatNum& iv_min, FloatNum& iv_ma
190219 } else {
191220 GECODE_ME_CHECK (x1.gq (home,-1.0 ));
192221 GECODE_ME_CHECK (x1.lq (home,1.0 ));
222+ GECODE_ES_CHECK (dopropagate (home,x0,x1));
193223 (void ) new (home) Cos<A,B>(home,x0,x1);
194224 }
195225 return ES_OK;
@@ -210,21 +240,7 @@ void aSinProject(Rounding& r, const V& aSinIv, FloatNum& iv_min, FloatNum& iv_ma
210240 template <class A , class B >
211241 ExecStatus
212242 Cos<A,B>::propagate(Space& home, const ModEventDelta&) {
213- GECODE_ME_CHECK (x1.eq (home,cos (x0.val ())));
214- Rounding r;
215- FloatVal x0Trans = x0.val () + FloatVal::pi_half ();
216- int n_min = 2 *static_cast <int >(r.div_up (x0Trans.min (), pi_twice_upper ()));
217- int n_max = 2 *static_cast <int >(r.div_up (x0Trans.max (), pi_twice_upper ()));
218- if (x0Trans.min () < 0 ) n_min-=2 ;
219- if (x0Trans.max () < 0 ) n_max-=2 ;
220- FloatNum iv_min = r.sub_down (x0Trans.min (),r.mul_down (n_min, pi_upper ()));
221- FloatNum iv_max = r.sub_up (x0Trans.max (),r.mul_down (n_max, pi_upper ()));
222- aSinProject (r,x1,iv_min,iv_max,n_min,n_max);
223- FloatNum n_iv_min = r.add_down (iv_min,r.mul_down (n_min, pi_upper ()));
224- FloatNum n_iv_max = r.add_up (iv_max,r.mul_down (n_max, pi_upper ()));
225- if (n_iv_min > n_iv_max) return ES_FAILED;
226- GECODE_ME_CHECK (x0.eq (home,FloatVal (n_iv_min,n_iv_max) - FloatVal::pi_half ()));
227- GECODE_ME_CHECK (x1.eq (home,cos (x0.val ()))); // Redo sin because with x0 reduction, sin may be more accurate
243+ GECODE_ES_CHECK (dopropagate (home,x0,x1));
228244 return (x0.assigned ()) ? home.ES_SUBSUMED (*this ) : ES_FIX;
229245 }
230246
0 commit comments