diff --git a/Makefile.dep b/Makefile.dep index 1f26262fd6..667dad0729 100644 --- a/Makefile.dep +++ b/Makefile.dep @@ -2144,6 +2144,8 @@ gecode/int/count$(OBJSUFFIX) gecode/int/count$(SBJSUFFIX): \ ./gecode/support/sort.hpp ./gecode/support/static-stack.hpp ./gecode/support/thread.hpp \ ./gecode/support/thread/thread.hpp ./gecode/support/timer.hpp gecode/int/arithmetic$(OBJSUFFIX) gecode/int/arithmetic$(SBJSUFFIX): \ + ./gecode/int/arithmetic/gcd.hpp ./gecode/int/arithmetic/divides.hpp \ + ./gecode/int/arithmetic/product.hpp ./gecode/int/arithmetic/product-mod.hpp \ ./gecode/int.hh ./gecode/int/arithmetic.hh ./gecode/int/arithmetic/abs.hpp \ ./gecode/int/arithmetic/argmax.hpp ./gecode/int/arithmetic/divmod.hpp ./gecode/int/arithmetic/max.hpp \ ./gecode/int/arithmetic/mult.hpp ./gecode/int/arithmetic/nroot.hpp ./gecode/int/arithmetic/pow-ops.hpp \ @@ -4441,6 +4443,8 @@ gecode/int/branch/chb$(OBJSUFFIX) gecode/int/branch/chb$(SBJSUFFIX): \ ./gecode/support/sort.hpp ./gecode/support/static-stack.hpp ./gecode/support/thread.hpp \ ./gecode/support/thread/thread.hpp ./gecode/support/timer.hpp gecode/int/arithmetic/mult$(OBJSUFFIX) gecode/int/arithmetic/mult$(SBJSUFFIX): \ + ./gecode/int/arithmetic/gcd.hpp ./gecode/int/arithmetic/divides.hpp \ + ./gecode/int/arithmetic/product.hpp ./gecode/int/arithmetic/product-mod.hpp \ ./gecode/int.hh ./gecode/int/arithmetic.hh ./gecode/int/arithmetic/abs.hpp \ ./gecode/int/arithmetic/argmax.hpp ./gecode/int/arithmetic/divmod.hpp ./gecode/int/arithmetic/max.hpp \ ./gecode/int/arithmetic/mult.hpp ./gecode/int/arithmetic/nroot.hpp ./gecode/int/arithmetic/pow-ops.hpp \ diff --git a/Makefile.in b/Makefile.in index b446405125..2c67889a1f 100755 --- a/Makefile.in +++ b/Makefile.in @@ -352,7 +352,8 @@ INTHDR0 = \ idx-view.hh idx-view.hpp div.hh div.hpp \ exec.hh exec/when.hpp \ arithmetic/abs.hpp arithmetic/max.hpp arithmetic/argmax.hpp \ - arithmetic/mult.hpp arithmetic/divmod.hpp \ + arithmetic/mult.hpp arithmetic/gcd.hpp arithmetic/divides.hpp \ + arithmetic/product.hpp arithmetic/product-mod.hpp arithmetic/divmod.hpp \ arithmetic/pow-ops.hpp arithmetic/pow.hpp arithmetic/nroot.hpp \ bool/or.hpp bool/eq.hpp bool/lq.hpp bool/eqv.hpp bool/base.hpp \ bool/clause.hpp bool/ite.hpp \ diff --git a/changelog.in b/changelog.in index 5356b24303..26f0c092b1 100755 --- a/changelog.in +++ b/changelog.in @@ -73,6 +73,19 @@ Date: unreleased [DESCRIPTION] This is the development changelog for the next Gecode release. +[ENTRY] +Module: int +What: new +Rank: major +[DESCRIPTION] +Add ordinary and reified gcd, n-ary product, and n-ary product_mod +constraints, and reified divides. GCD uses nonnegative results and +gcd(0,0)=0; zero divides zero. Modular products accept a fixed positive +modulus or a variable modulus and use nonnegative Euclidean residues, +including for negative products. All reification modes are supported. +Propagation uses bounds and algebraic reasoning; the propagation-level +argument does not select different consistency strengths. + [RELEASE] Version: 6.4.0 Date: 2026-07-15 diff --git a/docs/integer-number-theory.md b/docs/integer-number-theory.md new file mode 100644 index 0000000000..323a35f1c1 --- /dev/null +++ b/docs/integer-number-theory.md @@ -0,0 +1,42 @@ +# Integer number theory constraints + +These APIs require ``. They use bounds and algebraic reasoning; +the `IntPropLevel` argument currently does not select different strengths. +They do not promise bounds or domain consistency. + +Inside a `Space` constructor, this model has `g=6`, `p=-216`, and `r=1`: + +```cpp +IntVar x(*this,-12,-12), y(*this,18,18); +IntVar g(*this,0,18), p(*this,-300,300), r(*this,0,6); +gcd(*this,x,y,g); +product(*this,IntVarArgs({x,y}),p); +product_mod(*this,IntVarArgs({x,y}),7,r); +``` + +Call `status()` on the space to run propagation. `product_mod` uses the +Euclidean residue: `-216 = -31*7 + 1`. In contrast, Gecode's `mod` uses a +dividend-signed remainder and would return `-6` for `-216 mod 7`. + +To require divisibility, pass a true reification variable: + +```cpp +BoolVar yes(*this,1,1); +IntVar divisor(*this,6,6), dividend(*this,-30,30); +divides(*this,divisor,dividend,Reify(yes)); +``` + +This relation means that some integer multiplier exists. Zero divides zero, +but does not divide a nonzero value. GCD is always nonnegative, with +`gcd(0,0)=0`. The exact product of an empty array is one; its modular product +is `1 mod m`, which is zero when `m=1`. + +Every reified overload supports equivalence (`RM_EQV`), `b` implying the +relation (`RM_IMP`), and the relation implying `b` (`RM_PMI`). For example, +`product_mod(*this,factors,m,result,Reify(enabled,RM_IMP))` requires the +modular relation only when `enabled=1`. + +A fixed integer modulus must lie in `1..Int::Limits::max`; an invalid constant +throws `Int::OutOfLimits` even when an implication is inactive. With an +`IntVar` modulus, positivity is part of the reified proposition. Consequently, +an inactive implication does not constrain that modulus or the result. diff --git a/gecode/int.hh b/gecode/int.hh index 2b6c2ea796..0c4d31282d 100755 --- a/gecode/int.hh +++ b/gecode/int.hh @@ -3031,6 +3031,118 @@ namespace Gecode { mult(Home home, IntVar x0, IntVar x1, IntVar x2, IntPropLevel ipl=IPL_DEF); + /** \brief Post propagator for \f$\gcd(x_0,x_1)=x_2\f$ + * + * The greatest common divisor is nonnegative, with + * \f$\gcd(0,0)=0\f$. Negative operands are interpreted by absolute + * value. Uses sound bounds and algebraic propagation; bounds consistency + * is not guaranteed. The propagation level \a ipl is currently ignored. + */ + GECODE_INT_EXPORT void + gcd(Home home, IntVar x0, IntVar x1, IntVar x2, + IntPropLevel ipl=IPL_DEF); + + /** \brief Post propagator for + * \f$(\gcd(x_0,x_1)=x_2)\leftrightarrow r\f$ + * + * Supports all reification modes. The greatest common divisor is + * nonnegative, with \f$\gcd(0,0)=0\f$. Uses conservative algebraic + * entailment and disentailment tests. The propagation level \a ipl is + * currently ignored. + */ + GECODE_INT_EXPORT void + gcd(Home home, IntVar x0, IntVar x1, IntVar x2, Reify r, + IntPropLevel ipl=IPL_DEF); + + /** \brief Reify whether \a divisor divides \a dividend + * + * Divisibility means that an integer \f$k\f$ exists such that + * \f$dividend=divisor\cdot k\f$. Consequently, zero divides zero, but + * zero does not divide a nonzero integer. Supports all reification modes. + * Uses bounds and conservative algebraic propagation. The propagation + * level \a ipl is currently ignored. + */ + GECODE_INT_EXPORT void + divides(Home home, IntVar divisor, IntVar dividend, Reify r, + IntPropLevel ipl=IPL_DEF); + + /** \brief Constrain \a y to the exact product of the variables in \a x + * + * The product of an empty array is one. + * Uses bounds and algebraic propagation, without guaranteeing bounds + * consistency. The propagation level \a ipl is currently ignored. + * \ingroup TaskModelInt + */ + GECODE_INT_EXPORT void + product(Home home, const IntVarArgs& x, IntVar y, + IntPropLevel ipl=IPL_DEF); + + /** \brief Reify whether \a y is the exact product of the variables in \a x + * + * The product of an empty array is one. + * Supports all reification modes using conservative algebraic tests. + * The propagation level \a ipl is currently ignored. + * \ingroup TaskModelInt + */ + GECODE_INT_EXPORT void + product(Home home, const IntVarArgs& x, IntVar y, Reify r, + IntPropLevel ipl=IPL_DEF); + + /** \brief Constrain \a y to the product of \a x modulo \a m + * + * The modulus \a m must be positive. The result uses the canonical + * Euclidean residue in the range zero through \a m minus one. The product + * of an empty array is one. + * Unlike mod(), a negative product still has a nonnegative residue. + * Uses bounds and algebraic propagation; \a ipl is currently ignored. + * Throws Int::OutOfLimits if \a m is nonpositive or exceeds Int::Limits::max. + * \ingroup TaskModelInt + */ + GECODE_INT_EXPORT void + product_mod(Home home, const IntVarArgs& x, int m, IntVar y, + IntPropLevel ipl=IPL_DEF); + + /** \brief Reify whether \a y is the product of \a x modulo \a m + * + * The modulus \a m must be positive. The result uses the canonical + * Euclidean residue in the range zero through \a m minus one. The product + * of an empty array is one. + * Supports all reification modes using conservative algebraic tests; + * \a ipl is currently ignored. + * Throws Int::OutOfLimits if \a m is nonpositive or exceeds Int::Limits::max, + * even for an inactive implication. + * \ingroup TaskModelInt + */ + GECODE_INT_EXPORT void + product_mod(Home home, const IntVarArgs& x, int m, IntVar y, Reify r, + IntPropLevel ipl=IPL_DEF); + + /** \brief Constrain \a y to the product of \a x modulo \a m + * + * The variable modulus \a m is constrained to be positive and \a y uses + * the canonical Euclidean residue, so that \f$0\leq y0\f$, the canonical range + * \f$0\leq y + ::post(home,x0,x1,x2,r.var()))); + break; + case RM_IMP: + GECODE_ES_FAIL((Arithmetic::ReGcd + ::post(home,x0,x1,x2,r.var()))); + break; + case RM_PMI: + GECODE_ES_FAIL((Arithmetic::ReGcd + ::post(home,x0,x1,x2,r.var()))); + break; + default: GECODE_NEVER; + } + } + + void + divides(Home home, IntVar divisor, IntVar dividend, Reify r, + IntPropLevel) { + using namespace Int; + GECODE_POST; + switch (r.mode()) { + case RM_EQV: + GECODE_ES_FAIL((Arithmetic::ReDivides + ::post(home,divisor,dividend,r.var()))); + break; + case RM_IMP: + GECODE_ES_FAIL((Arithmetic::ReDivides + ::post(home,divisor,dividend,r.var()))); + break; + case RM_PMI: + GECODE_ES_FAIL((Arithmetic::ReDivides + ::post(home,divisor,dividend,r.var()))); + break; + default: GECODE_NEVER; + } + } + + void + product(Home home, const IntVarArgs& x, IntVar y, IntPropLevel) { + using namespace Int; + GECODE_POST; + ViewArray xv(home,x); + GECODE_ES_FAIL(Arithmetic::Product::post(home,xv,y)); + } + + void + product(Home home, const IntVarArgs& x, IntVar y, Reify r, + IntPropLevel) { + using namespace Int; + GECODE_POST; + ViewArray xv(home,x); + switch (r.mode()) { + case RM_EQV: + GECODE_ES_FAIL((Arithmetic::ReProduct + ::post(home,xv,y,r.var()))); + break; + case RM_IMP: + GECODE_ES_FAIL((Arithmetic::ReProduct + ::post(home,xv,y,r.var()))); + break; + case RM_PMI: + GECODE_ES_FAIL((Arithmetic::ReProduct + ::post(home,xv,y,r.var()))); + break; + default: GECODE_NEVER; + } + } + + void + product_mod(Home home, const IntVarArgs& x, int m, IntVar y, + IntPropLevel) { + using namespace Int; + Limits::positive(m,"Int::product_mod"); + GECODE_POST; + IntView yv(y); + ViewArray xv(home,x); + GECODE_ES_FAIL(Arithmetic::ProductMod::post(home,xv,m,yv)); + } + + void + product_mod(Home home, const IntVarArgs& x, int m, IntVar y, Reify r, + IntPropLevel) { + using namespace Int; + Limits::positive(m,"Int::product_mod"); + GECODE_POST; + IntView yv(y); + ViewArray xv(home,x); + switch (r.mode()) { + case RM_EQV: + GECODE_ES_FAIL((Arithmetic::ReProductMod + ::post(home,xv,m,yv,r.var()))); + break; + case RM_IMP: + GECODE_ES_FAIL((Arithmetic::ReProductMod + ::post(home,xv,m,yv,r.var()))); + break; + case RM_PMI: + GECODE_ES_FAIL((Arithmetic::ReProductMod + ::post(home,xv,m,yv,r.var()))); + break; + default: GECODE_NEVER; + } + } + + void + product_mod(Home home, const IntVarArgs& x, IntVar m, IntVar y, + IntPropLevel) { + using namespace Int; + GECODE_POST; + ViewArray xv(home,x); + GECODE_ES_FAIL(Arithmetic::ProductModVar::post(home,xv,m,y)); + } + + void + product_mod(Home home, const IntVarArgs& x, IntVar m, IntVar y, Reify r, + IntPropLevel) { + using namespace Int; + GECODE_POST; + ViewArray xv(home,x); + switch (r.mode()) { + case RM_EQV: + GECODE_ES_FAIL((Arithmetic::ReProductModVar + ::post(home,xv,m,y,r.var()))); + break; + case RM_IMP: + GECODE_ES_FAIL((Arithmetic::ReProductModVar + ::post(home,xv,m,y,r.var()))); + break; + case RM_PMI: + GECODE_ES_FAIL((Arithmetic::ReProductModVar + ::post(home,xv,m,y,r.var()))); + break; + default: GECODE_NEVER; + } + } + void divmod(Home home, IntVar x0, IntVar x1, IntVar x2, IntVar x3, diff --git a/gecode/int/arithmetic.hh b/gecode/int/arithmetic.hh index 0f24c953bd..7638e7b949 100644 --- a/gecode/int/arithmetic.hh +++ b/gecode/int/arithmetic.hh @@ -766,6 +766,315 @@ namespace Gecode { namespace Int { namespace Arithmetic { #include +namespace Gecode { namespace Int { namespace Arithmetic { + + /** + * \brief Bounds propagator for \f$\gcd(x_0,x_1)=x_2\f$ + * + * Requires \code #include \endcode + * \ingroup FuncIntProp + */ + class Gcd : public TernaryPropagator { + protected: + using TernaryPropagator::x0; + using TernaryPropagator::x1; + using TernaryPropagator::x2; + /// Constructor for cloning \a p + Gcd(Space& home, Gcd& p); + /// Constructor for posting + Gcd(Home home, IntView x0, IntView x1, IntView x2); + public: + /// Post propagator + static ExecStatus post(Home home, IntView x0, IntView x1, IntView x2); + /// Copy propagator during cloning + virtual Actor* copy(Space& home); + /// Cost function + virtual PropCost cost(const Space& home, const ModEventDelta& med) const; + /// Perform propagation + virtual ExecStatus propagate(Space& home, const ModEventDelta& med); + }; + + /** \brief Reified bounds propagator for \f$\gcd(x_0,x_1)=x_2\f$ + * + * Requires \code #include \endcode + * \ingroup FuncIntProp + */ + template + class ReGcd : public Propagator { + protected: + /// Operands and greatest common divisor + IntView x0, x1, x2; + /// Reification control + BoolView b; + /// Constructor for posting + ReGcd(Home home, IntView x0, IntView x1, IntView x2, BoolView b); + /// Constructor for cloning \a p + ReGcd(Space& home, ReGcd& p); + public: + /// Post propagator + static ExecStatus post(Home home, IntView x0, IntView x1, IntView x2, + BoolView b); + /// Copy propagator during cloning + virtual Actor* copy(Space& home); + /// Cost function + virtual PropCost cost(const Space& home, const ModEventDelta& med) const; + /// Reschedule propagator + virtual void reschedule(Space& home); + /// Perform propagation + virtual ExecStatus propagate(Space& home, const ModEventDelta& med); + /// Delete propagator and return its size + virtual size_t dispose(Space& home); + }; + +}}} + +#include + +namespace Gecode { namespace Int { namespace Arithmetic { + + /** + * \brief Reified bounds propagator for divisibility + * + * The relation is \f$\exists k\in\mathbb Z:x_1=x_0k\f$. + * Requires \code #include \endcode + * \ingroup FuncIntProp + */ + template + class ReDivides : + public ReBinaryPropagator { + protected: + using ReBinaryPropagator::x0; + using ReBinaryPropagator::x1; + using ReBinaryPropagator::b; + /// Constructor for posting + ReDivides(Home home, IntView x0, IntView x1, BoolView b); + /// Constructor for cloning \a p + ReDivides(Space& home, ReDivides& p); + public: + /// Post propagator + static ExecStatus post(Home home, IntView x0, IntView x1, BoolView b); + /// Copy propagator during cloning + virtual Actor* copy(Space& home); + /// Cost function + virtual PropCost cost(const Space& home, const ModEventDelta& med) const; + /// Perform propagation + virtual ExecStatus propagate(Space& home, const ModEventDelta& med); + /// Delete propagator and return its size + virtual size_t dispose(Space& home); + }; + +}}} + +#include + +namespace Gecode { namespace Int { namespace Arithmetic { + + /** \brief Bounds propagator for an exact n-ary product + * + * Requires \code #include \endcode + * \ingroup FuncIntProp + */ + class Product : public NaryOnePropagator { + protected: + using NaryOnePropagator::x; + using NaryOnePropagator::y; + /// Whether the remaining product is negated + bool neg; + /// Constructor for posting + Product(Home home, ViewArray& x, IntView y, bool neg); + /// Constructor for cloning \a p + Product(Space& home, Product& p); + public: + /// Post propagator + static ExecStatus post(Home home, ViewArray& x, IntView y, + bool neg=false); + /// Copy propagator during cloning + virtual Actor* copy(Space& home); + /// Cost function + virtual PropCost cost(const Space& home, const ModEventDelta& med) const; + /// Perform propagation + virtual ExecStatus propagate(Space& home, const ModEventDelta& med); + /// Delete propagator and return its size + virtual size_t dispose(Space& home); + }; + + /** \brief Reified bounds propagator for an exact n-ary product + * + * Requires \code #include \endcode + * \ingroup FuncIntProp + */ + template + class ReProduct : public Propagator { + protected: + /// Factors + ViewArray x; + /// Exact product + IntView y; + /// Reification control + BoolView b; + /// Constructor for posting + ReProduct(Home home, ViewArray& x, IntView y, BoolView b); + /// Constructor for cloning \a p + ReProduct(Space& home, ReProduct& p); + public: + /// Post propagator + static ExecStatus post(Home home, ViewArray& x, IntView y, + BoolView b); + /// Copy propagator during cloning + virtual Actor* copy(Space& home); + /// Cost function + virtual PropCost cost(const Space& home, const ModEventDelta& med) const; + /// Reschedule propagator + virtual void reschedule(Space& home); + /// Perform propagation + virtual ExecStatus propagate(Space& home, const ModEventDelta& med); + /// Delete propagator and return its size + virtual size_t dispose(Space& home); + }; + +}}} + +#include + +namespace Gecode { namespace Int { namespace Arithmetic { + + /** \brief Bounds propagator for a fixed-modulus n-ary product + * + * Requires \code #include \endcode + * \ingroup FuncIntProp + */ + class ProductMod : public NaryOnePropagator { + protected: + using NaryOnePropagator::x; + using NaryOnePropagator::y; + /// Positive modulus + int m; + /// Constructor for posting + ProductMod(Home home, ViewArray& x, int m, IntView y); + /// Constructor for cloning \a p + ProductMod(Space& home, ProductMod& p); + public: + /// Post propagator + static ExecStatus post(Home home, ViewArray& x, int m, IntView y); + /// Copy propagator during cloning + virtual Actor* copy(Space& home); + /// Cost function + virtual PropCost cost(const Space& home, const ModEventDelta& med) const; + /// Perform propagation + virtual ExecStatus propagate(Space& home, const ModEventDelta& med); + }; + + /** \brief Bounds propagator for a variable-modulus n-ary product + * + * Requires \code #include \endcode + * \ingroup FuncIntProp + */ + class ProductModVar : public Propagator { + protected: + /// Factors + ViewArray x; + /// Positive modulus + IntView m; + /// Canonical residue + IntView y; + /// Constructor for posting + ProductModVar(Home home, ViewArray& x, IntView m, IntView y); + /// Constructor for cloning \a p + ProductModVar(Space& home, ProductModVar& p); + public: + /// Post propagator + static ExecStatus post(Home home, ViewArray& x, + IntView m, IntView y); + /// Copy propagator during cloning + virtual Actor* copy(Space& home); + /// Cost function + virtual PropCost cost(const Space& home, const ModEventDelta& med) const; + /// Reschedule propagator + virtual void reschedule(Space& home); + /// Perform propagation + virtual ExecStatus propagate(Space& home, const ModEventDelta& med); + /// Delete propagator and return its size + virtual size_t dispose(Space& home); + }; + + /** \brief Reified bounds propagator for a variable-modulus product + * + * Requires \code #include \endcode + * \ingroup FuncIntProp + */ + template + class ReProductModVar : public Propagator { + protected: + /// Factors + ViewArray x; + /// Modulus (positivity is part of the reified proposition) + IntView m; + /// Canonical residue + IntView y; + /// Reification control + BoolView b; + /// Constructor for posting + ReProductModVar(Home home, ViewArray& x, IntView m, IntView y, + BoolView b); + /// Constructor for cloning \a p + ReProductModVar(Space& home, ReProductModVar& p); + public: + /// Post propagator + static ExecStatus post(Home home, ViewArray& x, IntView m, + IntView y, BoolView b); + /// Copy propagator during cloning + virtual Actor* copy(Space& home); + /// Cost function + virtual PropCost cost(const Space& home, const ModEventDelta& med) const; + /// Reschedule propagator + virtual void reschedule(Space& home); + /// Perform propagation + virtual ExecStatus propagate(Space& home, const ModEventDelta& med); + /// Delete propagator and return its size + virtual size_t dispose(Space& home); + }; + + /** \brief Reified bounds propagator for a fixed-modulus n-ary product + * + * Requires \code #include \endcode + * \ingroup FuncIntProp + */ + template + class ReProductMod : public Propagator { + protected: + /// Factors + ViewArray x; + /// Canonical residue + IntView y; + /// Reification control + BoolView b; + /// Positive modulus + int m; + /// Constructor for posting + ReProductMod(Home home, ViewArray& x, int m, IntView y, + BoolView b); + /// Constructor for cloning \a p + ReProductMod(Space& home, ReProductMod& p); + public: + /// Post propagator + static ExecStatus post(Home home, ViewArray& x, int m, IntView y, + BoolView b); + /// Copy propagator during cloning + virtual Actor* copy(Space& home); + /// Cost function + virtual PropCost cost(const Space& home, const ModEventDelta& med) const; + /// Reschedule propagator + virtual void reschedule(Space& home); + /// Perform propagation + virtual ExecStatus propagate(Space& home, const ModEventDelta& med); + /// Delete propagator and return its size + virtual size_t dispose(Space& home); + }; + +}}} + +#include + namespace Gecode { namespace Int { namespace Arithmetic { /** @@ -857,4 +1166,3 @@ namespace Gecode { namespace Int { namespace Arithmetic { #endif // STATISTICS: int-prop - diff --git a/gecode/int/arithmetic/divides.hpp b/gecode/int/arithmetic/divides.hpp new file mode 100644 index 0000000000..fd44f0584f --- /dev/null +++ b/gecode/int/arithmetic/divides.hpp @@ -0,0 +1,230 @@ +/* -*- mode: C++; c-basic-offset: 2; indent-tabs-mode: nil -*- */ +/* + * Main authors: + * Mikael Zayenz Lagerkvist + * + * Copyright: + * Mikael Zayenz Lagerkvist, 2026 + * + * This file is part of Gecode, the generic constraint + * development environment: + * http://www.gecode.dev + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE + * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION + * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + * + */ + +#include + +namespace Gecode { namespace Int { namespace Arithmetic { + + forceinline bool + divides_value(int divisor, int dividend) { + return (divisor == 0) ? (dividend == 0) : (dividend % divisor == 0); + } + + /// Tighten a dividend to the extreme multiples of a nonzero divisor. + inline ExecStatus + divides_multiple_bounds(Home home, IntView dividend, int divisor) { + const int a=(divisor < 0) ? -divisor : divisor; + assert(a > 0); + const long long int l = + ceil_div_xx(static_cast(dividend.min()), + static_cast(a)) * a; + const long long int u = + floor_div_xx(static_cast(dividend.max()), + static_cast(a)) * a; + if (l > u) + return ES_FAILED; + GECODE_ME_CHECK(dividend.gq(home,static_cast(l))); + GECODE_ME_CHECK(dividend.lq(home,static_cast(u))); + return ES_OK; + } + + /// Enforce divisibility with bounds reasoning only. + inline ExecStatus + divides_bnd(Space& home, IntView divisor, IntView dividend) { + if (divisor == dividend) + return ES_OK; + if (dividend.assigned() && (dividend.val() == 0)) + return ES_OK; + if (divisor.assigned()) { + const int d=divisor.val(); + if (d == 0) { + GECODE_ME_CHECK(dividend.eq(home,0)); + return ES_OK; + } + if ((d == 1) || (d == -1)) + return ES_OK; + GECODE_ES_CHECK(divides_multiple_bounds(home,dividend,d)); + } + + // A nonzero dividend excludes a bounds-visible zero divisor. + if ((dividend.min() > 0) || (dividend.max() < 0)) { + if ((divisor.min() == 0) && (divisor.max() > 0)) + GECODE_ME_CHECK(divisor.gq(home,1)); + else if ((divisor.max() == 0) && (divisor.min() < 0)) + GECODE_ME_CHECK(divisor.lq(home,-1)); + } + + // A divisor of an assigned nonzero value has no greater magnitude. + if (dividend.assigned() && (dividend.val() != 0)) { + const int a=(dividend.val() < 0) ? -dividend.val() : dividend.val(); + GECODE_ME_CHECK(divisor.gq(home,-a)); + GECODE_ME_CHECK(divisor.lq(home,a)); + if ((divisor.min() == 0) && (divisor.max() > 0)) + GECODE_ME_CHECK(divisor.gq(home,1)); + else if ((divisor.max() == 0) && (divisor.min() < 0)) + GECODE_ME_CHECK(divisor.lq(home,-1)); + } + + if (divisor.assigned() && dividend.assigned()) + return divides_value(divisor.val(),dividend.val()) ? ES_OK : ES_FAILED; + return ES_FIX; + } + + inline RelTest + divides_status(const IntView& x0, const IntView& x1) { + if (x0 == x1) + return RT_TRUE; + if (x1.assigned() && (x1.val() == 0)) + return RT_TRUE; + if (x0.assigned() && ((x0.val() == 1) || (x0.val() == -1))) + return RT_TRUE; + if (x0.assigned() && (x0.val() == 0) && !x1.in(0)) + return RT_FALSE; + if (x0.assigned() && (x0.val() != 0)) { + const int a=(x0.val() < 0) ? -x0.val() : x0.val(); + const long long int l = + ceil_div_xx(static_cast(x1.min()), + static_cast(a)) * a; + const long long int u = + floor_div_xx(static_cast(x1.max()), + static_cast(a)) * a; + if (l > u) + return RT_FALSE; + } + if (x1.assigned() && (x1.val() != 0)) { + const int a=(x1.val() < 0) ? -x1.val() : x1.val(); + if ((x0.min() > a) || (x0.max() < -a)) + return RT_FALSE; + } + if (x0.assigned() && x1.assigned()) + return divides_value(x0.val(),x1.val()) ? RT_TRUE : RT_FALSE; + return RT_MAYBE; + } + + template + forceinline + ReDivides::ReDivides(Home home, IntView y0, IntView y1, BoolView b0) + : ReBinaryPropagator(home,y0,y1,b0) { + home.notice(*this,AP_WEAKLY); + } + + template + inline ExecStatus + ReDivides::post(Home home, IntView x0, IntView x1, BoolView b) { + if (b.one() && (rm == RM_PMI)) + return ES_OK; + if (b.zero() && (rm == RM_IMP)) + return ES_OK; + RelTest rt = divides_status(x0,x1); + if (rt == RT_TRUE) { + if (rm != RM_IMP) + GECODE_ME_CHECK(b.one(home)); + return ES_OK; + } + if (rt == RT_FALSE) { + if (rm != RM_PMI) + GECODE_ME_CHECK(b.zero(home)); + return ES_OK; + } + (void) new (home) ReDivides(home,x0,x1,b); + return ES_OK; + } + + template + forceinline + ReDivides::ReDivides(Space& home, ReDivides& p) + : ReBinaryPropagator(home,p) {} + + template + forceinline Actor* + ReDivides::copy(Space& home) { + return new (home) ReDivides(home,*this); + } + + template + forceinline PropCost + ReDivides::cost(const Space&, const ModEventDelta&) const { + return PropCost::binary(PropCost::LO); + } + + template + inline ExecStatus + ReDivides::propagate(Space& home, const ModEventDelta&) { + if (b.one()) { + if (rm == RM_PMI) + return home.ES_SUBSUMED(*this); + const unsigned int s0=x0.size(), s1=x1.size(); + GECODE_ES_CHECK(divides_bnd(home,x0,x1)); + const RelTest rt=divides_status(x0,x1); + if (rt == RT_FALSE) + return ES_FAILED; + if (rt == RT_TRUE) + return home.ES_SUBSUMED(*this); + // Revisit sparse endpoints and divisors assigned by this call. + return ((s0 != x0.size()) || (s1 != x1.size())) ? ES_NOFIX : ES_FIX; + } + if (b.zero()) { + if (rm == RM_IMP) + return home.ES_SUBSUMED(*this); + const RelTest rt=divides_status(x0,x1); + if (rt == RT_TRUE) return ES_FAILED; + if (rt == RT_FALSE) + return home.ES_SUBSUMED(*this); + return ES_FIX; + } + + switch (divides_status(x0,x1)) { + case RT_TRUE: + if (rm != RM_IMP) + GECODE_ME_CHECK(b.one_none(home)); + break; + case RT_FALSE: + if (rm != RM_PMI) + GECODE_ME_CHECK(b.zero_none(home)); + break; + case RT_MAYBE: + return ES_FIX; + default: GECODE_NEVER; + } + return home.ES_SUBSUMED(*this); + } + + template + forceinline size_t + ReDivides::dispose(Space& home) { + home.ignore(*this,AP_WEAKLY); + (void) ReBinaryPropagator::dispose(home); + return sizeof(*this); + } + +}}} diff --git a/gecode/int/arithmetic/gcd.hpp b/gecode/int/arithmetic/gcd.hpp new file mode 100644 index 0000000000..4bc028cd2e --- /dev/null +++ b/gecode/int/arithmetic/gcd.hpp @@ -0,0 +1,347 @@ +/* -*- mode: C++; c-basic-offset: 2; indent-tabs-mode: nil -*- */ +/* + * Main authors: + * Mikael Zayenz Lagerkvist + * + * Copyright: + * Mikael Zayenz Lagerkvist, 2026 + * + * This file is part of Gecode, the generic constraint + * development environment: + * http://www.gecode.dev + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE + * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION + * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + * + */ + +#include + +namespace Gecode { namespace Int { namespace Arithmetic { + + forceinline int + gcd_value(int a, int b) { + unsigned int x = static_cast((a < 0) ? -a : a); + unsigned int y = static_cast((b < 0) ? -b : b); + while (y != 0U) { + unsigned int t = x % y; + x = y; y = t; + } + return static_cast(x); + } + + forceinline int + gcd_max_abs(const IntView& x) { + return std::max((x.min() < 0) ? -x.min() : x.min(), + (x.max() < 0) ? -x.max() : x.max()); + } + + forceinline bool + gcd_excludes_zero_bnd(const IntView& x) { + return (x.min() > 0) || (x.max() < 0); + } + + /// Upper bound on a gcd from the operand bounds. + forceinline int + gcd_upper(const IntView& x0, const IntView& x1) { + const int a0=gcd_max_abs(x0); + const int a1=gcd_max_abs(x1); + int u=std::max(a0,a1); + if (gcd_excludes_zero_bnd(x0)) u=std::min(u,a0); + if (gcd_excludes_zero_bnd(x1)) u=std::min(u,a1); + return u; + } + + /// Tighten an operand to the extreme multiples of a positive gcd. + inline ExecStatus + gcd_multiple_bounds(Home home, IntView x, int g) { + assert(g > 0); + const long long int l = + ceil_div_xx(static_cast(x.min()), + static_cast(g)) * g; + const long long int u = + floor_div_xx(static_cast(x.max()), + static_cast(g)) * g; + if (l > u) + return ES_FAILED; + GECODE_ME_CHECK(x.gq(home,static_cast(l))); + GECODE_ME_CHECK(x.lq(home,static_cast(u))); + return ES_OK; + } + + /// Status of abs(x0)=x1 using assignments and interval bounds only. + inline RelTest + gcd_abs_status(const IntView& x0, const IntView& x1) { + if (x1.max() < 0) + return RT_FALSE; + const int l = ((x0.min() <= 0) && (x0.max() >= 0)) ? 0 : + std::min((x0.min() < 0) ? -x0.min() : x0.min(), + (x0.max() < 0) ? -x0.max() : x0.max()); + const int u=gcd_max_abs(x0); + if ((x1.max() < l) || (x1.min() > u)) + return RT_FALSE; + if (x0 == x1) { + if (x0.min() >= 0) + return RT_TRUE; + if (x0.max() < 0) + return RT_FALSE; + } + if (x0.assigned()) { + const int a=(x0.val() < 0) ? -x0.val() : x0.val(); + if (!x1.in(a)) + return RT_FALSE; + return x1.assigned() ? RT_TRUE : RT_MAYBE; + } + if (x1.assigned() && (x1.val() == 0) && !x0.in(0)) + return RT_FALSE; + return RT_MAYBE; + } + + inline RelTest + gcd_status(const IntView& x0, const IntView& x1, const IntView& x2) { + if (x2.max() < 0) + return RT_FALSE; + if (x2.min() > gcd_upper(x0,x1)) + return RT_FALSE; + if (x0 == x1) + return gcd_abs_status(x0,x2); + if (x0.assigned() && (x0.val() == 0)) + return gcd_abs_status(x1,x2); + if (x1.assigned() && (x1.val() == 0)) + return gcd_abs_status(x0,x2); + if ((x0.assigned() && ((x0.val() == 1) || (x0.val() == -1))) || + (x1.assigned() && ((x1.val() == 1) || (x1.val() == -1)))) { + if (!x2.in(1)) return RT_FALSE; + return x2.assigned() ? RT_TRUE : RT_MAYBE; + } + if (x2.assigned()) { + const int g=x2.val(); + if (g == 0) { + if (!x0.in(0) || !x1.in(0)) return RT_FALSE; + } else if (g > 0) { + if ((x0.assigned() && ((x0.val() % g) != 0)) || + (x1.assigned() && ((x1.val() % g) != 0))) + return RT_FALSE; + } + } + if (x0.assigned() && x1.assigned()) { + const int g=gcd_value(x0.val(),x1.val()); + if (!x2.in(g)) return RT_FALSE; + return x2.assigned() ? RT_TRUE : RT_MAYBE; + } + return RT_MAYBE; + } + + forceinline + Gcd::Gcd(Home home, IntView y0, IntView y1, IntView y2) + : TernaryPropagator(home,y0,y1,y2) {} + + inline ExecStatus + Gcd::post(Home home, IntView x0, IntView x1, IntView x2) { + if (x0 == x1) + return AbsBnd::post(home,x0,x2); + if (x0.assigned() && (x0.val() == 0)) + return AbsBnd::post(home,x1,x2); + if (x1.assigned() && (x1.val() == 0)) + return AbsBnd::post(home,x0,x2); + if ((x0.assigned() && ((x0.val() == 1) || (x0.val() == -1))) || + (x1.assigned() && ((x1.val() == 1) || (x1.val() == -1)))) { + GECODE_ME_CHECK(x2.eq(home,1)); + return ES_OK; + } + GECODE_ME_CHECK(x2.gq(home,0)); + if (gcd_excludes_zero_bnd(x0) || gcd_excludes_zero_bnd(x1)) + GECODE_ME_CHECK(x2.gq(home,1)); + GECODE_ME_CHECK(x2.lq(home,gcd_upper(x0,x1))); + if (x2.assigned()) { + const int g=x2.val(); + if (g == 0) { + GECODE_ME_CHECK(x0.eq(home,0)); + GECODE_ME_CHECK(x1.eq(home,0)); + return ES_OK; + } + GECODE_ES_CHECK(gcd_multiple_bounds(home,x0,g)); + GECODE_ES_CHECK(gcd_multiple_bounds(home,x1,g)); + if (x0.assigned() && (x0.val() == 0)) + return AbsBnd::post(home,x1,x2); + if (x1.assigned() && (x1.val() == 0)) + return AbsBnd::post(home,x0,x2); + } + if (x0.assigned() && x1.assigned()) { + GECODE_ME_CHECK(x2.eq(home,gcd_value(x0.val(),x1.val()))); + return ES_OK; + } + (void) new (home) Gcd(home,x0,x1,x2); + return ES_OK; + } + + forceinline + Gcd::Gcd(Space& home, Gcd& p) + : TernaryPropagator(home,p) {} + + forceinline Actor* + Gcd::copy(Space& home) { + return new (home) Gcd(home,*this); + } + + forceinline PropCost + Gcd::cost(const Space&, const ModEventDelta&) const { + return PropCost::ternary(PropCost::LO); + } + + inline ExecStatus + Gcd::propagate(Space& home, const ModEventDelta&) { + const unsigned int s0=x0.size(), s1=x1.size(), s2=x2.size(); + if (x0 == x1) + GECODE_REWRITE(*this,AbsBnd::post(home(*this),x0,x2)); + if (x0.assigned() && (x0.val() == 0)) + GECODE_REWRITE(*this,AbsBnd::post(home(*this),x1,x2)); + if (x1.assigned() && (x1.val() == 0)) + GECODE_REWRITE(*this,AbsBnd::post(home(*this),x0,x2)); + if ((x0.assigned() && ((x0.val() == 1) || (x0.val() == -1))) || + (x1.assigned() && ((x1.val() == 1) || (x1.val() == -1)))) { + GECODE_ME_CHECK(x2.eq(home,1)); + return home.ES_SUBSUMED(*this); + } + GECODE_ME_CHECK(x2.gq(home,0)); + if (gcd_excludes_zero_bnd(x0) || gcd_excludes_zero_bnd(x1)) + GECODE_ME_CHECK(x2.gq(home,1)); + GECODE_ME_CHECK(x2.lq(home,gcd_upper(x0,x1))); + if (x2.assigned()) { + const int g=x2.val(); + if (g == 0) { + GECODE_ME_CHECK(x0.eq(home,0)); + GECODE_ME_CHECK(x1.eq(home,0)); + return home.ES_SUBSUMED(*this); + } + GECODE_ES_CHECK(gcd_multiple_bounds(home,x0,g)); + GECODE_ES_CHECK(gcd_multiple_bounds(home,x1,g)); + if (x0.assigned() && (x0.val() == 0)) + GECODE_REWRITE(*this,AbsBnd::post(home(*this),x1,x2)); + if (x1.assigned() && (x1.val() == 0)) + GECODE_REWRITE(*this,AbsBnd::post(home(*this),x0,x2)); + } + if (x0.assigned() && x1.assigned()) { + GECODE_ME_CHECK(x2.eq(home,gcd_value(x0.val(),x1.val()))); + return home.ES_SUBSUMED(*this); + } + // A bounds update can jump over a hole to another non-multiple. + return ((s0 != x0.size()) || (s1 != x1.size()) || (s2 != x2.size())) + ? ES_NOFIX : ES_FIX; + } + + template + forceinline + ReGcd::ReGcd(Home home, IntView y0, IntView y1, IntView y2, BoolView b0) + : Propagator(home), x0(y0), x1(y1), x2(y2), b(b0) { + home.notice(*this,AP_WEAKLY); + x0.subscribe(home,*this,PC_INT_BND); + x1.subscribe(home,*this,PC_INT_BND); + x2.subscribe(home,*this,PC_INT_BND); + b.subscribe(home,*this,PC_BOOL_VAL); + } + + template + inline ExecStatus + ReGcd::post(Home home, IntView x0, IntView x1, IntView x2, BoolView b) { + if (b.one()) { + if (rm == RM_PMI) + return ES_OK; + return Gcd::post(home,x0,x1,x2); + } + if (b.zero() && (rm == RM_IMP)) + return ES_OK; + (void) new (home) ReGcd(home,x0,x1,x2,b); + return ES_OK; + } + + template + forceinline + ReGcd::ReGcd(Space& home, ReGcd& p) + : Propagator(home,p) { + x0.update(home,p.x0); x1.update(home,p.x1); x2.update(home,p.x2); + b.update(home,p.b); + } + + template + forceinline Actor* + ReGcd::copy(Space& home) { + return new (home) ReGcd(home,*this); + } + + template + forceinline PropCost + ReGcd::cost(const Space&, const ModEventDelta&) const { + return PropCost::ternary(PropCost::LO); + } + + template + forceinline void + ReGcd::reschedule(Space& home) { + x0.reschedule(home,*this,PC_INT_BND); + x1.reschedule(home,*this,PC_INT_BND); + x2.reschedule(home,*this,PC_INT_BND); + b.reschedule(home,*this,PC_BOOL_VAL); + } + + template + inline ExecStatus + ReGcd::propagate(Space& home, const ModEventDelta&) { + if (b.one()) { + if (rm == RM_PMI) + return home.ES_SUBSUMED(*this); + GECODE_REWRITE(*this,Gcd::post(home(*this),x0,x1,x2)); + } + if (b.zero()) { + if (rm == RM_IMP) + return home.ES_SUBSUMED(*this); + const RelTest rt=gcd_status(x0,x1,x2); + if (rt == RT_TRUE) return ES_FAILED; + if (rt == RT_FALSE) return home.ES_SUBSUMED(*this); + return ES_FIX; + } + + switch (gcd_status(x0,x1,x2)) { + case RT_FALSE: + if (rm != RM_PMI) + GECODE_ME_CHECK(b.zero(home)); + return home.ES_SUBSUMED(*this); + case RT_TRUE: + if (rm != RM_IMP) + GECODE_ME_CHECK(b.one(home)); + return home.ES_SUBSUMED(*this); + case RT_MAYBE: + return ES_FIX; + default: GECODE_NEVER; + } + } + + template + forceinline size_t + ReGcd::dispose(Space& home) { + x0.cancel(home,*this,PC_INT_BND); + x1.cancel(home,*this,PC_INT_BND); + x2.cancel(home,*this,PC_INT_BND); + b.cancel(home,*this,PC_BOOL_VAL); + home.ignore(*this,AP_WEAKLY); + (void) Propagator::dispose(home); + return sizeof(*this); + } + +}}} diff --git a/gecode/int/arithmetic/product-mod.hpp b/gecode/int/arithmetic/product-mod.hpp new file mode 100644 index 0000000000..5aa52c49e3 --- /dev/null +++ b/gecode/int/arithmetic/product-mod.hpp @@ -0,0 +1,874 @@ +/* -*- mode: C++; c-basic-offset: 2; indent-tabs-mode: nil -*- */ +/* + * Main authors: + * Mikael Zayenz Lagerkvist + * + * Copyright: + * Mikael Zayenz Lagerkvist, 2026 + * + * This file is part of Gecode, the generic constraint + * development environment: + * http://www.gecode.dev + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE + * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION + * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include + +namespace Gecode { namespace Int { namespace Arithmetic { + + /// Compute the canonical residue without signed overflow. + inline int + product_mod_value(const int* v, int n, int m) { + long long int p = 1 % m; + for (int i=0; i(v[i]) % m; + if (q < 0) + q += m; + // Both factors are smaller than Limits::max, so this fits in int64. + p = (p*q) % m; + } + return static_cast(p); + } + + /// Return exact representable product bounds, or false on saturation. + inline bool + product_mod_interval(const ViewArray& x, ProductInterval& p, + int omit=-1) { + p.min=1; p.max=1; + for (int i=0; i(p.min) * x[i].min(), + static_cast(p.min) * x[i].max(), + static_cast(p.max) * x[i].min(), + static_cast(p.max) * x[i].max() + }; + const long long int l = *std::min_element(q,q+4); + const long long int u = *std::max_element(q,q+4); + if (!Limits::valid(l) || !Limits::valid(u)) + return false; + p.min=static_cast(l); p.max=static_cast(u); + } + return true; + } + + /// Return whether an assigned factor makes the residue identically zero. + inline bool + product_mod_zero(const ViewArray& x, int m) { + for (int i=0; i& x, int m, int& residue) { + Region r; + int* v = r.alloc(x.size()); + for (int i=0; i& x, const IntView& y, int m) { + if ((y.max() < 0) || (y.min() >= m)) + return RT_FALSE; + const bool zero = (m == 1) || product_mod_zero(x,m); + if (zero) { + if (!y.in(0)) return RT_FALSE; + return y.assigned() ? RT_TRUE : RT_MAYBE; + } + int residue; + if (product_mod_assigned(x,m,residue)) { + if (!y.in(residue)) return RT_FALSE; + return y.assigned() ? RT_TRUE : RT_MAYBE; + } + ProductInterval p; + if (product_mod_interval(x,p)) { + const long long int kl = floor_div_xx + (static_cast(p.min),static_cast(m)); + const long long int ku = floor_div_xx + (static_cast(p.max),static_cast(m)); + if (kl == ku) { + const long long int shift=kl*m; + if ((static_cast(p.max)-shift < y.min()) || + (static_cast(p.min)-shift > y.max())) + return RT_FALSE; + } + } + return RT_MAYBE; + } + + /// Extended Euclid for a modular inverse (arguments are coprime). + inline long long int + product_mod_inverse(long long int a, long long int m) { + long long int old_r=a, r=m, old_s=1, s=0; + while (r != 0) { + const long long int q=old_r/r; + const long long int nr=old_r-q*r; old_r=r; r=nr; + const long long int ns=old_s-q*s; old_s=s; s=ns; + } + old_s %= m; + return old_s < 0 ? old_s+m : old_s; + } + + inline long long int + product_mod_gcd(long long int a, long long int b) { + while (b != 0) { + const long long int r=a%b; a=b; b=r; + } + return a < 0 ? -a : a; + } + + forceinline + ProductMod::ProductMod(Home home, ViewArray& z, int modulus, + IntView w) + : NaryOnePropagator(home,z,w), m(modulus) {} + + inline ExecStatus + ProductMod::post(Home home, ViewArray& x, int m, IntView y) { + GECODE_ME_CHECK(y.gq(home,0)); + GECODE_ME_CHECK(y.lq(home,m-1)); + if (m == 1) { + GECODE_ME_CHECK(y.eq(home,0)); + return ES_OK; + } + if (x.size() == 0) { + GECODE_ME_CHECK(y.eq(home,1 % m)); + return ES_OK; + } + if (product_mod_zero(x,m)) { + GECODE_ME_CHECK(y.eq(home,0)); + return ES_OK; + } + int residue; + if (product_mod_assigned(x,m,residue)) { + GECODE_ME_CHECK(y.eq(home,residue)); + return ES_OK; + } + (void) new (home) ProductMod(home,x,m,y); + return ES_OK; + } + + forceinline + ProductMod::ProductMod(Space& home, ProductMod& p) + : NaryOnePropagator(home,p), m(p.m) {} + + forceinline Actor* + ProductMod::copy(Space& home) { + return new (home) ProductMod(home,*this); + } + + forceinline PropCost + ProductMod::cost(const Space&, const ModEventDelta&) const { + return PropCost::quadratic(PropCost::LO,x.size()+1); + } + + inline ExecStatus + ProductMod::propagate(Space& home, const ModEventDelta&) { + GECODE_ME_CHECK(y.gq(home,0)); + GECODE_ME_CHECK(y.lq(home,m-1)); + + if (product_mod_zero(x,m)) { + GECODE_ME_CHECK(y.eq(home,0)); + return home.ES_SUBSUMED(*this); + } + + int residue; + if (product_mod_assigned(x,m,residue)) { + GECODE_ME_CHECK(y.eq(home,residue)); + return home.ES_SUBSUMED(*this); + } + + bool modified; + do { + modified=false; + ProductInterval p; + if (product_mod_interval(x,p)) { + const long long int kl = floor_div_xx + (static_cast(p.min),static_cast(m)); + const long long int ku = floor_div_xx + (static_cast(p.max),static_cast(m)); + if (kl == ku) { + const long long int shift=kl*m; + { + const long long int l=static_cast(p.min)-shift; + ModEvent me=y.gq(home,static_cast(l)); + if (me_failed(me)) return ES_FAILED; + modified |= me_modified(me); + } + { + const long long int u=static_cast(p.max)-shift; + ModEvent me=y.lq(home,static_cast(u)); + if (me_failed(me)) return ES_FAILED; + modified |= me_modified(me); + } + + // Within one quotient band, invert the exact-product interval. + const long long int tmin=static_cast(y.min())+shift; + const long long int tmax=static_cast(y.max())+shift; + for (int i=0; i 0) || (q.max < 0))) { + long long int c[4] = { + ceil_div_xx(tmin,static_cast(q.min)), + ceil_div_xx(tmin,static_cast(q.max)), + ceil_div_xx(tmax,static_cast(q.min)), + ceil_div_xx(tmax,static_cast(q.max)) + }; + long long int f[4] = { + floor_div_xx(tmin,static_cast(q.min)), + floor_div_xx(tmin,static_cast(q.max)), + floor_div_xx(tmax,static_cast(q.min)), + floor_div_xx(tmax,static_cast(q.max)) + }; + long long int l=*std::min_element(c,c+4); + long long int u=*std::max_element(f,f+4); + l=std::max(l,static_cast(Limits::min)); + u=std::min(u,static_cast(Limits::max)); + if (l > u) return ES_FAILED; + { + ModEvent me=x[i].gq(home,static_cast(l)); + if (me_failed(me)) return ES_FAILED; + modified |= me_modified(me); + } + { + ModEvent me=x[i].lq(home,static_cast(u)); + if (me_failed(me)) return ES_FAILED; + modified |= me_modified(me); + } + } + } + } + } + } while (modified); + + // If only one factor is free and the result is assigned, solve the + // resulting linear congruence and tighten to its extreme solutions. + bool congruence_modified=false; + if (y.assigned()) { + int free=-1; + long long int cofactor=1 % m; + for (int i=0; i(x[i].val()) % m; + if (r < 0) r += m; + cofactor=(cofactor*r) % m; + } else if (free < 0) { + free=i; + } else { + free=-2; break; + } + } + if (free >= 0) { + const long long int g=product_mod_gcd(cofactor,m); + if ((y.val() % g) != 0) + return ES_FAILED; + const long long int step=m/g; + long long int r=0; + if (step > 1) { + const long long int a=cofactor/g; + const long long int b=y.val()/g; + r=(product_mod_inverse(a % step,step) * b) % step; + if (r < 0) r += step; + } + const long long int l=r + ceil_div_xx + (static_cast(x[free].min())-r,step)*step; + const long long int u=r + floor_div_xx + (static_cast(x[free].max())-r,step)*step; + if (l > u) return ES_FAILED; + { + ModEvent me=x[free].gq(home,static_cast(l)); + if (me_failed(me)) return ES_FAILED; + congruence_modified |= me_modified(me); + } + { + ModEvent me=x[free].lq(home,static_cast(u)); + if (me_failed(me)) return ES_FAILED; + congruence_modified |= me_modified(me); + } + } + } + + if (product_mod_assigned(x,m,residue)) { + GECODE_ME_CHECK(y.eq(home,residue)); + return home.ES_SUBSUMED(*this); + } + return congruence_modified ? ES_NOFIX : ES_FIX; + } + + /// Whether the modulus occurs among the factors. + inline bool + product_mod_var_mod_factor(const ViewArray& x, const IntView& m) { + for (int i=0; i& x, long long int& p) { + long long int q=1; + for (int i=0; i(-(d+1))+1ULL + : static_cast(d); + assert(ad > 0); + + least=lower; + greatest=upper; + + // A positive divisor cannot exceed the absolute difference. + if (ad <= static_cast(Limits::max)) + greatest=std::min(greatest,static_cast(ad)); + if (least > greatest) + return false; + + // If floor(ad/m) is fixed throughout the remaining interval, there is + // at most one possible divisor. This is constant-time quotient reasoning; + // it does not inspect the values in the modulus domain. + const unsigned long long int q0 = + ad / static_cast(greatest); + const unsigned long long int q1 = + ad / static_cast(least); + if (q0 == q1) { + if ((q0 == 0) || ((ad % q0) != 0)) + return false; + const unsigned long long int candidate=ad/q0; + if ((candidate < static_cast(least)) || + (candidate > static_cast(greatest)) || + (candidate > static_cast(Limits::max))) + return false; + least=greatest=static_cast(candidate); + } + return true; + } + + /// Propagate a nonnegative representable product with a variable modulus. + inline ExecStatus + product_mod_var_ranges(Home home, ViewArray& x, IntView m, + IntView y, bool& modified) { + modified=false; + for (int i=0; i(p.min) / m.max(); + const long long int kmax = + static_cast(p.max) / m.min(); + if (kmin != kmax) + return ES_OK; + + const long long int k=kmin; + const long long int rmin = + static_cast(p.min)-k*m.max(); + const long long int rmax = + static_cast(p.max)-k*m.min(); + { + ModEvent me=y.gq(home,static_cast(rmin)); + if (me_failed(me)) return ES_FAILED; + modified |= me_modified(me); + } + { + ModEvent me=y.lq(home,static_cast(rmax)); + if (me_failed(me)) return ES_FAILED; + modified |= me_modified(me); + } + { + ModEvent me=m.gq(home,y.min()+1); + if (me_failed(me)) return ES_FAILED; + modified |= me_modified(me); + } + + const long long int tmin=k*m.min()+y.min(); + const long long int tmax=k*m.max()+y.max(); + for (int i=0; i 0) + return ES_FAILED; + continue; + } + const long long int candidate_min=ceil_div_xx(tmin, + static_cast(q.max)); + const long long int candidate_max=(q.min == 0) ? x[i].max() : + floor_div_xx(tmax,static_cast(q.min)); + const long long int l=std::max( + candidate_min,static_cast(x[i].min())); + const long long int u=std::min( + candidate_max,static_cast(x[i].max())); + if (l > u) + return ES_FAILED; + { + ModEvent me=x[i].gq(home,static_cast(l)); + if (me_failed(me)) return ES_FAILED; + modified |= me_modified(me); + } + { + ModEvent me=x[i].lq(home,static_cast(u)); + if (me_failed(me)) return ES_FAILED; + modified |= me_modified(me); + } + } + return ES_OK; + } + + /// Determine algebraic status; evaluate the residue only when fully assigned. + inline RelTest + product_mod_var_status(const ViewArray& x, const IntView& m, + const IntView& y) { + if ((m == y) || (m.max() <= 0) || (y.max() < 0) || + (y.min() >= m.max())) + return RT_FALSE; + bool zero=product_mod_var_mod_factor(x,m) || + (m.assigned() && (m.val() == 1)); + for (int i=0; !zero && (i 0) && y.assigned()) return RT_TRUE; + return RT_MAYBE; + } + if (x.size() == 0) { + const bool zero=m.in(1) && y.in(0); + const bool one=(m.max() >= 2) && y.in(1); + if (!zero && !one) return RT_FALSE; + if ((m.min() > 0) && + ((m.assigned() && (m.val() == 1) && y.assigned() && + (y.val() == 0)) || + ((m.min() >= 2) && y.assigned() && (y.val() == 1)))) + return RT_TRUE; + return RT_MAYBE; + } + if (y.assigned()) { + long long int p; + if (product_mod_var_exact(x,p) && + !Limits::overflow_sub(p,static_cast(y.val()))) { + const long long int d=p-y.val(); + if (d == 0) + return (y.val() >= 0) && (m.min() > y.val()) + ? RT_TRUE : RT_MAYBE; + int least, greatest; + if (!product_mod_var_divisor_bounds + (std::max(m.min(),y.val()+1),m.max(),d,least,greatest)) + return RT_FALSE; + if ((least == greatest) && !m.in(least)) + return RT_FALSE; + } + } + bool assigned=m.assigned() && y.assigned(); + for (int i=0; assigned && (i= m.val())) + return RT_FALSE; + int residue; + (void) product_mod_assigned(x,m.val(),residue); + return residue == y.val() ? RT_TRUE : RT_FALSE; + } + + forceinline + ProductModVar::ProductModVar(Home home, ViewArray& z, + IntView modulus, IntView w) + : Propagator(home), x(z), m(modulus), y(w) { + home.notice(*this,AP_WEAKLY); + x.subscribe(home,*this,PC_INT_BND); + m.subscribe(home,*this,PC_INT_BND); + y.subscribe(home,*this,PC_INT_BND); + } + + inline ExecStatus + ProductModVar::post(Home home, ViewArray& x, IntView m, IntView y) { + if (m == y) + return ES_FAILED; + GECODE_ME_CHECK(m.gq(home,1)); + GECODE_ME_CHECK(y.gq(home,0)); + GECODE_ME_CHECK(y.lq(home,m.max()-1)); + GECODE_ME_CHECK(m.gq(home,y.min()+1)); + if (product_mod_var_mod_factor(x,m)) { + GECODE_ME_CHECK(y.eq(home,0)); + return ES_OK; + } + if (x.size() == 0) { + GECODE_ME_CHECK(y.lq(home,1)); + if (!y.in(0)) { + GECODE_ME_CHECK(m.gq(home,2)); + GECODE_ME_CHECK(y.eq(home,1)); + } else if (!y.in(1)) { + GECODE_ME_CHECK(m.eq(home,1)); + GECODE_ME_CHECK(y.eq(home,0)); + } + } + if (m.assigned()) + return ProductMod::post(home,x,m.val(),y); + (void) new (home) ProductModVar(home,x,m,y); + return ES_OK; + } + + forceinline + ProductModVar::ProductModVar(Space& home, ProductModVar& p) + : Propagator(home,p) { + x.update(home,p.x); m.update(home,p.m); y.update(home,p.y); + } + + forceinline Actor* + ProductModVar::copy(Space& home) { + return new (home) ProductModVar(home,*this); + } + + forceinline PropCost + ProductModVar::cost(const Space&, const ModEventDelta&) const { + return PropCost::quadratic(PropCost::HI,x.size()+2); + } + + forceinline void + ProductModVar::reschedule(Space& home) { + x.reschedule(home,*this,PC_INT_BND); + m.reschedule(home,*this,PC_INT_BND); + y.reschedule(home,*this,PC_INT_BND); + } + + inline ExecStatus + ProductModVar::propagate(Space& home, const ModEventDelta&) { + GECODE_ME_CHECK(m.gq(home,1)); + GECODE_ME_CHECK(y.gq(home,0)); + GECODE_ME_CHECK(y.lq(home,m.max()-1)); + GECODE_ME_CHECK(m.gq(home,y.min()+1)); + if (product_mod_var_mod_factor(x,m)) { + GECODE_ME_CHECK(y.eq(home,0)); + return home.ES_SUBSUMED(*this); + } + for (int i=0; i= 2) { + GECODE_ME_CHECK(y.eq(home,1)); + return home.ES_SUBSUMED(*this); + } + if (!m.in(1)) { + GECODE_ME_CHECK(y.eq(home,1)); + return home.ES_SUBSUMED(*this); + } + } + if (m.assigned()) + GECODE_REWRITE(*this,ProductMod::post(home(*this),x,m.val(),y)); + + bool modified; + do { + GECODE_ES_CHECK(product_mod_var_ranges(home,x,m,y,modified)); + } while (modified); + if (m.assigned()) + GECODE_REWRITE(*this,ProductMod::post(home(*this),x,m.val(),y)); + + // With an assigned product and result, m must divide product-result. + long long int p; + if (y.assigned() && product_mod_var_exact(x,p)) { + if (!Limits::overflow_sub(p,static_cast(y.val()))) { + const long long int d=p-y.val(); + if (d == 0) + return home.ES_SUBSUMED(*this); + int least, greatest; + if (!product_mod_var_divisor_bounds + (m.min(),m.max(),d,least,greatest)) + return ES_FAILED; + GECODE_ME_CHECK(m.gq(home,least)); + GECODE_ME_CHECK(m.lq(home,greatest)); + } + } + if (m.assigned()) + GECODE_REWRITE(*this,ProductMod::post(home(*this),x,m.val(),y)); + bool assigned = m.assigned() && y.assigned(); + for (int i=0; assigned && (i + forceinline + ReProductModVar::ReProductModVar(Home home, ViewArray& z, + IntView modulus, IntView w, BoolView c) + : Propagator(home), x(z), m(modulus), y(w), b(c) { + home.notice(*this,AP_WEAKLY); + x.subscribe(home,*this,PC_INT_VAL); + m.subscribe(home,*this,PC_INT_BND); + y.subscribe(home,*this,PC_INT_BND); + b.subscribe(home,*this,PC_BOOL_VAL); + } + + template + inline ExecStatus + ReProductModVar::post(Home home, ViewArray& x, IntView m, + IntView y, BoolView b) { + if (b.one()) { + if (rm == RM_PMI) return ES_OK; + return ProductModVar::post(home,x,m,y); + } + if (b.zero() && (rm == RM_IMP)) + return ES_OK; + (void) new (home) ReProductModVar(home,x,m,y,b); + return ES_OK; + } + + template + forceinline + ReProductModVar::ReProductModVar(Space& home, + ReProductModVar& p) + : Propagator(home,p) { + x.update(home,p.x); m.update(home,p.m); y.update(home,p.y); + b.update(home,p.b); + } + + template + forceinline Actor* + ReProductModVar::copy(Space& home) { + return new (home) ReProductModVar(home,*this); + } + + template + forceinline PropCost + ReProductModVar::cost(const Space&, const ModEventDelta&) const { + return PropCost::linear(PropCost::HI,x.size()+3); + } + + template + forceinline void + ReProductModVar::reschedule(Space& home) { + x.reschedule(home,*this,PC_INT_VAL); + m.reschedule(home,*this,PC_INT_BND); + y.reschedule(home,*this,PC_INT_BND); + b.reschedule(home,*this,PC_BOOL_VAL); + } + + template + inline ExecStatus + ReProductModVar::propagate(Space& home, const ModEventDelta&) { + if (b.one()) { + if (rm == RM_PMI) + return home.ES_SUBSUMED(*this); + GECODE_REWRITE(*this,ProductModVar::post(home(*this),x,m,y)); + } + if (b.zero()) { + if (rm == RM_IMP) + return home.ES_SUBSUMED(*this); + const RelTest rt=product_mod_var_status(x,m,y); + if (rt == RT_TRUE) + return ES_FAILED; + if (rt == RT_FALSE) + return home.ES_SUBSUMED(*this); + return ES_FIX; + } + + const RelTest rt=product_mod_var_status(x,m,y); + switch (rt) { + case RT_TRUE: + if (rm != RM_IMP) GECODE_ME_CHECK(b.one_none(home)); + return home.ES_SUBSUMED(*this); + case RT_FALSE: + if (rm != RM_PMI) GECODE_ME_CHECK(b.zero_none(home)); + return home.ES_SUBSUMED(*this); + case RT_MAYBE: + return ES_FIX; + default: GECODE_NEVER; + } + } + + template + forceinline size_t + ReProductModVar::dispose(Space& home) { + x.cancel(home,*this,PC_INT_VAL); + m.cancel(home,*this,PC_INT_BND); + y.cancel(home,*this,PC_INT_BND); + b.cancel(home,*this,PC_BOOL_VAL); + home.ignore(*this,AP_WEAKLY); + (void) Propagator::dispose(home); + return sizeof(*this); + } + + template + forceinline + ReProductMod::ReProductMod(Home home, ViewArray& z, + int modulus, IntView w, BoolView c) + : Propagator(home), x(z), y(w), b(c), m(modulus) { + home.notice(*this,AP_WEAKLY); + x.subscribe(home,*this,PC_INT_BND); + y.subscribe(home,*this,PC_INT_BND); + b.subscribe(home,*this,PC_BOOL_VAL); + } + + template + inline ExecStatus + ReProductMod::post(Home home, ViewArray& x, int m, IntView y, + BoolView b) { + if (b.one() && (rm == RM_PMI)) + return ES_OK; + if (b.zero() && (rm == RM_IMP)) + return ES_OK; + if (x.size() == 0) + return Rel::ReEqDomInt::post(home,y,1 % m,b); + switch (product_mod_status(x,y,m)) { + case RT_TRUE: + if (rm != RM_IMP) GECODE_ME_CHECK(b.one(home)); + return ES_OK; + case RT_FALSE: + if (rm != RM_PMI) GECODE_ME_CHECK(b.zero(home)); + return ES_OK; + case RT_MAYBE: + break; + default: GECODE_NEVER; + } + (void) new (home) ReProductMod(home,x,m,y,b); + return ES_OK; + } + + template + forceinline + ReProductMod::ReProductMod(Space& home, ReProductMod& p) + : Propagator(home,p), m(p.m) { + x.update(home,p.x); y.update(home,p.y); b.update(home,p.b); + } + + template + forceinline Actor* + ReProductMod::copy(Space& home) { + return new (home) ReProductMod(home,*this); + } + + template + forceinline PropCost + ReProductMod::cost(const Space&, const ModEventDelta&) const { + return PropCost::linear(PropCost::HI,x.size()+2); + } + + template + forceinline void + ReProductMod::reschedule(Space& home) { + x.reschedule(home,*this,PC_INT_BND); + y.reschedule(home,*this,PC_INT_BND); + b.reschedule(home,*this,PC_BOOL_VAL); + } + + template + inline ExecStatus + ReProductMod::propagate(Space& home, const ModEventDelta&) { + if (b.one()) { + if (rm == RM_PMI) + return home.ES_SUBSUMED(*this); + GECODE_REWRITE(*this,ProductMod::post(home(*this),x,m,y)); + } + if (b.zero()) { + if (rm == RM_IMP) + return home.ES_SUBSUMED(*this); + switch (product_mod_status(x,y,m)) { + case RT_TRUE: return ES_FAILED; + case RT_FALSE: return home.ES_SUBSUMED(*this); + case RT_MAYBE: return ES_FIX; + default: GECODE_NEVER; + } + } + switch (product_mod_status(x,y,m)) { + case RT_TRUE: + if (rm != RM_IMP) GECODE_ME_CHECK(b.one_none(home)); + return home.ES_SUBSUMED(*this); + case RT_FALSE: + if (rm != RM_PMI) GECODE_ME_CHECK(b.zero_none(home)); + return home.ES_SUBSUMED(*this); + case RT_MAYBE: + return ES_FIX; + default: GECODE_NEVER; + } + } + + template + forceinline size_t + ReProductMod::dispose(Space& home) { + x.cancel(home,*this,PC_INT_BND); + y.cancel(home,*this,PC_INT_BND); + b.cancel(home,*this,PC_BOOL_VAL); + home.ignore(*this,AP_WEAKLY); + (void) Propagator::dispose(home); + return sizeof(*this); + } + +}}} diff --git a/gecode/int/arithmetic/product.hpp b/gecode/int/arithmetic/product.hpp new file mode 100644 index 0000000000..0d676409a2 --- /dev/null +++ b/gecode/int/arithmetic/product.hpp @@ -0,0 +1,627 @@ +/* -*- mode: C++; c-basic-offset: 2; indent-tabs-mode: nil -*- */ +/* + * Main authors: + * Mikael Zayenz Lagerkvist + * + * Copyright: + * Mikael Zayenz Lagerkvist, 2026 + * + * This file is part of Gecode, the generic constraint + * development environment: + * http://www.gecode.dev + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice shall be + * included in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE + * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION + * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION + * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include + +namespace Gecode { namespace Int { namespace Arithmetic { + + /// A representable, conservative interval for an exact product. + struct ProductInterval { + int min; + int max; + }; + + /// Multiply two representable intervals, clipping out-of-range endpoints. + forceinline ProductInterval + product_interval_mul(ProductInterval a, int bmin, int bmax) { + long long int p[4] = { + static_cast(a.min) * bmin, + static_cast(a.min) * bmax, + static_cast(a.max) * bmin, + static_cast(a.max) * bmax + }; + long long int l = *std::min_element(p,p+4); + long long int u = *std::max_element(p,p+4); + l = std::max(static_cast(Limits::min), + std::min(l,static_cast(Limits::max))); + u = std::max(static_cast(Limits::min), + std::min(u,static_cast(Limits::max))); + ProductInterval r = {static_cast(l),static_cast(u)}; + return r; + } + + /// Multiply two representable intervals. + forceinline ProductInterval + product_interval_mul(ProductInterval a, ProductInterval b) { + return product_interval_mul(a,b.min,b.max); + } + + /// Compute a nonnegative power, saturating at the integer limit. + forceinline int + product_power_abs(int x, int n) { + const long long int a = x < 0 ? -static_cast(x) : x; + long long int p=1; + for (int i=0; i static_cast(Limits::max)/a)) + return Limits::max; + p *= a; + } + return static_cast(p); + } + + /// Bounds for a repeated occurrence viewed as an integer power. + forceinline ProductInterval + product_power_interval(IntView x, int n) { + const int l=product_power_abs(x.min(),n); + const int u=product_power_abs(x.max(),n); + if ((n & 1) != 0) { + ProductInterval p = {x.min() < 0 ? -l : l, + x.max() < 0 ? -u : u}; + return p; + } + ProductInterval p = { + ((x.min() <= 0) && (x.max() >= 0)) ? 0 : std::min(l,u), + std::max(l,u) + }; + return p; + } + + forceinline bool + product_power_le(int x, int n, int limit) { + long long int p=1; + for (int i=0; i static_cast(limit)/x)) + return false; + p *= x; + } + return true; + } + + /// Floor and ceiling of a nonnegative integer root. + inline int + product_floor_root(int x, int n) { + int l=0, u=x; + while (l < u) { + const int m=l+(u-l+1)/2; + if (product_power_le(m,n,x)) l=m; else u=m-1; + } + return l; + } + + forceinline int + product_ceil_root(int x, int n) { + if (x <= 0) return 0; + return product_floor_root(x-1,n)+1; + } + + /// Compute the hull of quotients by a zero-free interval. + forceinline ProductInterval + product_quotient_interval(int ymin, int ymax, + int qmin, int qmax) { + long long int c[4] = { + ceil_div_xx(static_cast(ymin), + static_cast(qmin)), + ceil_div_xx(static_cast(ymin), + static_cast(qmax)), + ceil_div_xx(static_cast(ymax), + static_cast(qmin)), + ceil_div_xx(static_cast(ymax), + static_cast(qmax)) + }; + long long int f[4] = { + floor_div_xx(static_cast(ymin), + static_cast(qmin)), + floor_div_xx(static_cast(ymin), + static_cast(qmax)), + floor_div_xx(static_cast(ymax), + static_cast(qmin)), + floor_div_xx(static_cast(ymax), + static_cast(qmax)) + }; + const long long int l = std::max + (static_cast(Limits::min),std::min + (*std::min_element(c,c+4),static_cast(Limits::max))); + const long long int u = std::max + (static_cast(Limits::min),std::min + (*std::max_element(f,f+4),static_cast(Limits::max))); + ProductInterval r = {static_cast(l),static_cast(u)}; + return r; + } + + /// Compute product bounds, optionally omitting one factor. + inline ProductInterval + product_interval(const ViewArray& x, int omit=-1) { + ProductInterval r = {1,1}; + for (int i=0; i(v[i]))) + return false; + q *= static_cast(v[i]); + } + if (!Limits::valid(q)) + return false; + p = static_cast(q); + return true; + } + + /// Determine the relation when bounds or assignments prove it. + inline RelTest + product_status(const ViewArray& x, const IntView& y) { + if ((x.size() == 1) && (x[0] == y)) + return RT_TRUE; + for (int i=0; i y.max())) + return RT_FALSE; + bool assigned = true; + for (int i=0; assigned && (i(x.size()); + for (int i=0; i& z, IntView w, bool n) + : NaryOnePropagator(home,z,w), neg(n) { + home.notice(*this,AP_WEAKLY); + } + + inline ExecStatus + Product::post(Home home, ViewArray& x, IntView y, bool neg) { + for (int i=x.size(); i--;) { + if (!x[i].assigned()) + continue; + if (x[i].val() == 0) { + GECODE_ME_CHECK(y.eq(home,0)); + return ES_OK; + } + if ((x[i].val() == 1) || (x[i].val() == -1)) { + neg ^= x[i].val() == -1; + x.move_lst(i); + } + } + if (x.size() == 0) { + GECODE_ME_CHECK(y.eq(home,neg ? -1 : 1)); + return ES_OK; + } + if (x.size() == 1) { + if (neg) + return Rel::EqBnd::post(home,x[0],MinusView(y)); + return Rel::EqBnd::post(home,x[0],y); + } + // Keep zero-aware inverse bounds when both intervals contain zero; + // binary multiplication does less pruning in that case. + if ((x.size() == 2) && !neg && + ((x[0] == x[1]) || (x[0].min() > 0) || (x[0].max() < 0) || + (x[1].min() > 0) || (x[1].max() < 0))) + return MultBnd::post(home,x[0],x[1],y); + (void) new (home) Product(home,x,y,neg); + return ES_OK; + } + + forceinline + Product::Product(Space& home, Product& p) + : NaryOnePropagator(home,p), neg(p.neg) {} + + forceinline Actor* + Product::copy(Space& home) { + return new (home) Product(home,*this); + } + + forceinline PropCost + Product::cost(const Space&, const ModEventDelta&) const { + return PropCost::quadratic(PropCost::LO,x.size()+1); + } + + forceinline size_t + Product::dispose(Space& home) { + home.ignore(*this,AP_WEAKLY); + (void) NaryOnePropagator::dispose(home); + return sizeof(*this); + } + + inline ExecStatus + Product::propagate(Space& home, const ModEventDelta&) { + // Absorb a fixed zero and rewrite when new units have appeared. + int units=0; + bool next_neg=neg; + for (int i=x.size(); i--;) { + if (!x[i].assigned()) + continue; + if (x[i].val() == 0) { + GECODE_ME_CHECK(y.eq(home,0)); + return home.ES_SUBSUMED(*this); + } + if ((x[i].val() == 1) || (x[i].val() == -1)) { + next_neg ^= x[i].val() == -1; + units++; + } + } + if (units > 0) { + ViewArray z(home,x.size()-units); + int j=0; + for (int i=0; i= 0) { + if (y.assigned() && (y.val() == 0)) + return home.ES_SUBSUMED(*this); + const int unit=neg ? -1 : 1; + if (!y.in(0)) { + ViewArray z(home,x.size()-1); + for (int i=0, j=0; i q.max)) { + GECODE_ME_CHECK(y.eq(home,0)); + return home.ES_SUBSUMED(*this); + } + } + + // A nonzero result excludes a bounds-visible zero endpoint. + if (!y.in(0)) + for (int i=0; i 0)) + GECODE_ME_CHECK(x[i].gq(home,1)); + else if ((x[i].max() == 0) && (x[i].min() < 0)) + GECODE_ME_CHECK(x[i].lq(home,-1)); + } + + // Zero can only arise from a zero factor. + if (y.assigned() && (y.val() == 0)) { + int zero=-1; + for (int i=0; i= 0) { zero=-2; break; } + zero=i; + } + if (zero == -1) + return ES_FAILED; + if (zero >= 0) { + GECODE_ME_CHECK(x[zero].eq(home,0)); + return home.ES_SUBSUMED(*this); + } + } + + // Stable one-sided signs determine aggregate parity and zero possibility. + bool stable=true, strict=true, negative=neg; + for (int i=0; i= 0) { + strict &= x[i].min() > 0; + } else { + stable=false; strict=false; break; + } + } + if (stable) { + if (negative) + GECODE_ME_CHECK(y.lq(home,strict ? -1 : 0)); + else + GECODE_ME_CHECK(y.gq(home,strict ? 1 : 0)); + } + + // View identities do not change during propagation. Group once, then + // reuse the powers for forward and inverse bounds in each iteration. + Region r; + int* representative=r.alloc(x.size()); + int* exponent=r.alloc(x.size()); + int groups=0; + for (int i=0; i(groups); + ProductInterval* prefix=r.alloc(groups+1); + ProductInterval* suffix=r.alloc(groups+1); + bool modified; + do { + modified = false; + prefix[0] = ProductInterval{1,1}; + for (int g=0; g= 0) && y.in(0)) + continue; + if ((q.min == 0) && (q.max == 0)) + return ES_FAILED; + + bool have=false; + ProductInterval d = {Limits::max,Limits::min}; + if (q.min < 0) { + ProductInterval n = product_quotient_interval + (neg ? -y.max() : y.min(),neg ? -y.min() : y.max(), + q.min,std::min(q.max,-1)); + d=n; have=true; + } + if (q.max > 0) { + ProductInterval p = product_quotient_interval + (neg ? -y.max() : y.min(),neg ? -y.min() : y.max(), + std::max(q.min,1),q.max); + if (have) { + d.min=std::min(d.min,p.min); d.max=std::max(d.max,p.max); + } else { + d=p; have=true; + } + } + if (!have || (d.min > d.max)) + return ES_FAILED; + const int n=exponent[g]; + if (n == 1) { + ModEvent me=xi.gq(home,d.min); + if (me_failed(me)) return ES_FAILED; + modified |= me_modified(me); + me=xi.lq(home,d.max); + if (me_failed(me)) return ES_FAILED; + modified |= me_modified(me); + } else if ((n & 1) != 0) { + const int l = d.min < 0 ? + -product_floor_root(-d.min,n) : product_ceil_root(d.min,n); + const int u = d.max < 0 ? + -product_ceil_root(-d.max,n) : product_floor_root(d.max,n); + if (l > u) return ES_FAILED; + ModEvent me=xi.gq(home,l); + if (me_failed(me)) return ES_FAILED; + modified |= me_modified(me); + me=xi.lq(home,u); + if (me_failed(me)) return ES_FAILED; + modified |= me_modified(me); + } else { + if (d.max < 0) return ES_FAILED; + const int lo=product_ceil_root(std::max(d.min,0),n); + const int hi=product_floor_root(d.max,n); + if (lo > hi) return ES_FAILED; + bool have_negative=xi.min() <= -lo; + bool have_positive=xi.max() >= lo; + int l=Limits::max, u=Limits::min; + if (have_negative) { + l=std::max(xi.min(),-hi); u=std::min(xi.max(),-lo); + have_negative=l <= u; + } + if (have_positive) { + const int pl=std::max(xi.min(),lo); + const int pu=std::min(xi.max(),hi); + have_positive=pl <= pu; + if (have_positive) { + if (have_negative) { l=std::min(l,pl); u=std::max(u,pu); } + else { l=pl; u=pu; } + } + } + if (!have_negative && !have_positive) return ES_FAILED; + ModEvent me=xi.gq(home,l); + if (me_failed(me)) return ES_FAILED; + modified |= me_modified(me); + me=xi.lq(home,u); + if (me_failed(me)) return ES_FAILED; + modified |= me_modified(me); + } + } + } while (modified); + + // Result propagation above can have excluded zero during this call. + bool endpoint_modified=false; + if (!y.in(0)) + for (int i=0; i 0)) + me=x[i].gq(home,1); + else if ((x[i].max() == 0) && (x[i].min() < 0)) + me=x[i].lq(home,-1); + if (me_failed(me)) return ES_FAILED; + endpoint_modified |= me_modified(me); + } + if (endpoint_modified) + return ES_NOFIX; + + bool factors_assigned = true; + for (int i=0; factors_assigned && (i(x.size()); + for (int i=0; i + forceinline + ReProduct::ReProduct(Home home, ViewArray& z, IntView w, + BoolView c) + : Propagator(home), x(z), y(w), b(c) { + home.notice(*this,AP_WEAKLY); + x.subscribe(home,*this,PC_INT_BND); + y.subscribe(home,*this,PC_INT_BND); + b.subscribe(home,*this,PC_BOOL_VAL); + } + + template + inline ExecStatus + ReProduct::post(Home home, ViewArray& x, IntView y, BoolView b) { + if (b.one() && (rm == RM_PMI)) return ES_OK; + if (b.zero() && (rm == RM_IMP)) return ES_OK; + if (x.size() == 0) + return Rel::ReEqDomInt::post(home,y,1,b); + if (x.size() == 1) + return Rel::ReEqBnd::post(home,x[0],y,b); + (void) new (home) ReProduct(home,x,y,b); + return ES_OK; + } + + template + forceinline + ReProduct::ReProduct(Space& home, ReProduct& p) + : Propagator(home,p) { + x.update(home,p.x); y.update(home,p.y); b.update(home,p.b); + } + + template + forceinline Actor* + ReProduct::copy(Space& home) { + return new (home) ReProduct(home,*this); + } + + template + forceinline PropCost + ReProduct::cost(const Space&, const ModEventDelta&) const { + return PropCost::quadratic(PropCost::HI,x.size()+2); + } + + template + forceinline void + ReProduct::reschedule(Space& home) { + x.reschedule(home,*this,PC_INT_BND); + y.reschedule(home,*this,PC_INT_BND); + b.reschedule(home,*this,PC_BOOL_VAL); + } + + template + inline ExecStatus + ReProduct::propagate(Space& home, const ModEventDelta&) { + if (b.one()) { + if (rm == RM_PMI) return home.ES_SUBSUMED(*this); + GECODE_REWRITE(*this,Product::post(home(*this),x,y)); + } + if (b.zero()) { + if (rm == RM_IMP) return home.ES_SUBSUMED(*this); + switch (product_status(x,y)) { + case RT_TRUE: return ES_FAILED; + case RT_FALSE: return home.ES_SUBSUMED(*this); + case RT_MAYBE: return ES_FIX; + default: GECODE_NEVER; + } + } + switch (product_status(x,y)) { + case RT_TRUE: + if (rm != RM_IMP) GECODE_ME_CHECK(b.one_none(home)); + return home.ES_SUBSUMED(*this); + case RT_FALSE: + if (rm != RM_PMI) GECODE_ME_CHECK(b.zero_none(home)); + return home.ES_SUBSUMED(*this); + case RT_MAYBE: + return ES_FIX; + default: GECODE_NEVER; + } + } + + template + forceinline size_t + ReProduct::dispose(Space& home) { + x.cancel(home,*this,PC_INT_BND); + y.cancel(home,*this,PC_INT_BND); + b.cancel(home,*this,PC_BOOL_VAL); + home.ignore(*this,AP_WEAKLY); + (void) Propagator::dispose(home); + return sizeof(*this); + } + +}}} diff --git a/test/int/arithmetic.cpp b/test/int/arithmetic.cpp index 59ca93388e..6a9bb35c78 100644 --- a/test/int/arithmetic.cpp +++ b/test/int/arithmetic.cpp @@ -35,6 +35,7 @@ #include #include +#include #include @@ -48,6 +49,1521 @@ namespace Test { namespace Int { * \ingroup TaskTestInt */ //@{ + /// Compute the mathematical greatest common divisor for testing. + int gcd_value(int a, int b) { + a = (a < 0) ? -a : a; + b = (b < 0) ? -b : b; + while (b != 0) { + int t = a % b; + a = b; b = t; + } + return a; + } + + /// %Test for the ternary greatest-common-divisor constraint + // Reposting comparisons below are disabled where bounds subscriptions + // deliberately miss interior removals used by algebraic status tests. + // NumberTheoryLifecycle checks activation, cloning, and rescheduling + // separately; NumberTheorySparseBounds checks own-update fixpoints. + class GcdXYZ : public Test { + public: + /// Create and register test + GcdXYZ(const std::string& s, const Gecode::IntSet& d, + Gecode::IntPropLevel ipl, bool r=true) + : Test("Arithmetic::Gcd::XYZ::"+str(ipl)+"::"+s+(r ? "" : "::Plain"),3,d,r,ipl) { + contest=CTL_NONE; testfix=!r; + } + /// %Test whether \a x is solution + virtual bool solution(const Assignment& x) const { + return gcd_value(x[0],x[1]) == x[2]; + } + /// Post constraint on \a x + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { + Gecode::gcd(home,x[0],x[1],x[2],ipl); + } + /// Post reified constraint on \a x + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x, + Gecode::Reify r) { + Gecode::gcd(home,x[0],x[1],x[2],r,ipl); + } + }; + + /// %Test for gcd with identical operands + class GcdXXY : public Test { + public: + /// Create and register test + GcdXXY(const std::string& s, const Gecode::IntSet& d, + Gecode::IntPropLevel ipl, bool r=true) + : Test("Arithmetic::Gcd::XXY::"+str(ipl)+"::"+s+(r ? "" : "::Plain"),2,d,r,ipl) { + contest=CTL_NONE; testfix=!r; + } + /// %Test whether \a x is solution + virtual bool solution(const Assignment& x) const { + return gcd_value(x[0],x[0]) == x[1]; + } + /// Post constraint on \a x + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { + Gecode::gcd(home,x[0],x[0],x[1],ipl); + } + /// Post reified constraint on \a x + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x, + Gecode::Reify r) { + Gecode::gcd(home,x[0],x[0],x[1],r,ipl); + } + }; + + /// %Test for gcd with result aliased to the first operand + class GcdXYX : public Test { + public: + /// Create and register test + GcdXYX(const std::string& s, const Gecode::IntSet& d, + Gecode::IntPropLevel ipl, bool r=true) + : Test("Arithmetic::Gcd::XYX::"+str(ipl)+"::"+s+(r ? "" : "::Plain"),2,d,r,ipl) { + contest=CTL_NONE; testfix=!r; + } + /// %Test whether \a x is solution + virtual bool solution(const Assignment& x) const { + return gcd_value(x[0],x[1]) == x[0]; + } + /// Post constraint on \a x + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { + Gecode::gcd(home,x[0],x[1],x[0],ipl); + } + /// Post reified constraint on \a x + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x, + Gecode::Reify r) { + Gecode::gcd(home,x[0],x[1],x[0],r,ipl); + } + }; + + /// %Test for gcd with all variables aliased + class GcdXXX : public Test { + public: + /// Create and register test + GcdXXX(const std::string& s, const Gecode::IntSet& d, + Gecode::IntPropLevel ipl) + : Test("Arithmetic::Gcd::XXX::"+str(ipl)+"::"+s,1,d,true,ipl) { + contest=CTL_NONE; + } + /// %Test whether \a x is solution + virtual bool solution(const Assignment& x) const { + return gcd_value(x[0],x[0]) == x[0]; + } + /// Post constraint on \a x + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { + Gecode::gcd(home,x[0],x[0],x[0],ipl); + } + /// Post reified constraint on \a x + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x, + Gecode::Reify r) { + Gecode::gcd(home,x[0],x[0],x[0],r,ipl); + } + }; + + /// %Test for the reified divisibility constraint + // With divisor zero, removing interior zero from the dividend can + // strengthen a reposted status test without a bounds wakeup. + class DividesXY : public Test { + public: + /// Create and register test + DividesXY(const std::string& s, const Gecode::IntSet& d, + Gecode::IntPropLevel ipl) + : Test("Arithmetic::Divides::XY::"+str(ipl)+"::"+s, + 2,d,true,ipl) { contest=CTL_NONE; testfix=false; } + /// %Test whether \a x is solution + virtual bool solution(const Assignment& x) const { + return (x[0] == 0) ? (x[1] == 0) : (x[1] % x[0] == 0); + } + /// The constraint is deliberately exposed only as reified + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { + Gecode::BoolVar b(home,1,1); + Gecode::divides(home,x[0],x[1],Gecode::Reify(b),ipl); + } + /// Post reified constraint on \a x + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x, + Gecode::Reify r) { + Gecode::divides(home,x[0],x[1],r,ipl); + } + }; + + /// %Test divisibility with aliased operands + class DividesXX : public Test { + public: + /// Create and register test + DividesXX(const std::string& s, const Gecode::IntSet& d, + Gecode::IntPropLevel ipl) + : Test("Arithmetic::Divides::XX::"+str(ipl)+"::"+s, + 1,d,true,ipl) { contest=CTL_NONE; } + /// Every integer divides itself, including zero + virtual bool solution(const Assignment&) const { + return true; + } + /// The constraint is deliberately exposed only as reified + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { + Gecode::BoolVar b(home,1,1); + Gecode::divides(home,x[0],x[0],Gecode::Reify(b),ipl); + } + /// Post reified constraint on \a x + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x, + Gecode::Reify r) { + Gecode::divides(home,x[0],x[0],r,ipl); + } + }; + + /// Report a targeted arithmetic failure through the standard test log. + bool arithmetic_failed(const std::string& what, + const Gecode::VarArgArray& x) { + if (opt.log) + olog << what << ": " << x << std::endl; + return false; + } + + /// Sparse endpoints must reach a fixpoint after a divisor is assigned. + class NumberTheorySparseBounds : public ::Test::Base { + class TestSpace : public Gecode::Space { + public: + virtual Gecode::Space* copy(void) { return nullptr; } + }; + public: + NumberTheorySparseBounds(void) + : ::Test::Base("Int::Arithmetic::NumberTheorySparseBounds") {} + virtual bool run(void) { + using namespace Gecode; + for (int division=0; division<2; division++) { + TestSpace home; + IntVar x(home,IntSet({1,3,5,6})), y(home,10,10), g(home,1,2); + BoolVar b(home,1,1); + if (division) + divides(home,g,x,Reify(b)); + else + gcd(home,x,y,g); + if (home.status() == SS_FAILED) + return arithmetic_failed + ("initial propagation: x, y, g", + Gecode::IntVarArgs() << x << y << g); + rel(home,g,IRT_EQ,2); + if ((home.status() == SS_FAILED) || !x.assigned() || (x.val()!=6)) + return arithmetic_failed + ("assigned divisor: x, y, g", + Gecode::IntVarArgs() << x << y << g); + } + return true; + } + }; + + /// Activation and rescheduling survive cloning and arithmetic rewrites. + class NumberTheoryLifecycle : public ::Test::Base { + class TestSpace : public Gecode::Space { + public: + Gecode::IntVar x, c, y, m; + Gecode::BoolVar b; + TestSpace(int kind) + : x(*this,Gecode::IntSet({1,3,5,6})), + c(*this,kind == 0 ? 10 : 2,kind == 0 ? 10 : 2), + y(*this,kind == 0 ? 2 : (kind == 4 ? 1 : 12), + kind == 0 ? 2 : (kind == 4 ? 1 : 12)), + m(*this,kind == 4 ? 5 : 13,kind == 4 ? 7 : 17), b(*this,0,1) { + using namespace Gecode; + switch (kind) { + case 0: gcd(*this,x,c,y,Reify(b)); break; + case 1: divides(*this,c,x,Reify(b)); break; + case 2: product(*this,IntVarArgs({x,c}),y,Reify(b)); break; + case 3: product_mod(*this,IntVarArgs({x,c}),13,y,Reify(b)); break; + case 4: product_mod(*this,IntVarArgs({x,c}),m,y,Reify(b)); break; + default: GECODE_NEVER; + } + } + TestSpace(TestSpace& s) : Gecode::Space(s) { + x.update(*this,s.x); c.update(*this,s.c); y.update(*this,s.y); + m.update(*this,s.m); b.update(*this,s.b); + } + virtual Gecode::Space* copy(void) { return new TestSpace(*this); } + }; + public: + NumberTheoryLifecycle(void) + : ::Test::Base("Int::Arithmetic::NumberTheoryLifecycle") {} + virtual bool run(void) { + using namespace Gecode; + for (int kind=0; kind<5; kind++) { + TestSpace home(kind); + if (home.status() == SS_FAILED) + return arithmetic_failed + ("initial propagation: home.x, home.c, home.y, home.m", + Gecode::IntVarArgs() << home.x << home.c << home.y << home.m); + std::unique_ptr clone + (static_cast(home.clone())); + PropagatorGroup::all.disable(*clone); + rel(*clone,clone->b,IRT_EQ,1); + PropagatorGroup::all.enable(*clone); + if ((clone->status() == SS_FAILED) || + ((kind != 4) && (!clone->x.assigned() || (clone->x.val() != 6)))) + return arithmetic_failed + ("activation after cloning, kind "+str(kind)+": x, c, y, m", + Gecode::IntVarArgs() << clone->x << clone->c + << clone->y << clone->m); + // Clone the activated actor before fixing a variable modulus. + std::unique_ptr child + (static_cast(clone->clone())); + rel(*child,child->m,IRT_EQ,kind == 4 ? 5 : 13); + if ((child->status() == SS_FAILED) || !child->x.assigned() || + (child->x.val() != (kind == 4 ? 3 : 6)) || + home.b.assigned() || (home.x.size() != 4) || clone->m.assigned()) + return arithmetic_failed + ("fixed modulus in child, kind "+str(kind)+": child x, m; parent x, m", + Gecode::IntVarArgs() << child->x << child->m + << clone->x << clone->m); + } + return true; + } + }; + + /// Reified identities do not require the remaining factors to be fixed. + class NumberTheoryIdentities : public ::Test::Base { + class TestSpace : public Gecode::Space { + public: + virtual Gecode::Space* copy(void) { return nullptr; } + }; + public: + NumberTheoryIdentities(void) + : ::Test::Base("Int::Arithmetic::NumberTheoryIdentities") {} + virtual bool run(void) { + using namespace Gecode; + for (ReifyMode rm : {RM_EQV,RM_IMP,RM_PMI}) + for (int kind=0; kind<4; kind++) { + TestSpace home; + IntVar x(home,-7,7), zero(home,0,0), m(home,0,7), one(home,1,1); + BoolVar b(home,0,1); + if (kind == 0) product(home,IntVarArgs({zero,x}),zero,Reify(b,rm)); + if (kind == 1) product(home,IntVarArgs({x}),x,Reify(b,rm)); + if (kind == 2) product_mod(home,IntVarArgs({zero,x}),m,zero,Reify(b,rm)); + if (kind == 3) product_mod(home,IntVarArgs({x}),one,zero,Reify(b,rm)); + if (home.status() == SS_FAILED) + return arithmetic_failed + ("initial identity: x, zero, m, one", + Gecode::IntVarArgs() << x << zero << m << one); + if (kind == 2) { + // Zero residue alone does not establish a positive modulus. + if (b.assigned() || (m.min() != 0)) + return arithmetic_failed + ("modulus positivity unknown: x, zero, m, one", + Gecode::IntVarArgs() << x << zero << m << one); + rel(home,m,IRT_GQ,1); + if (home.status() == SS_FAILED) + return arithmetic_failed + ("positive modulus: x, zero, m, one", + Gecode::IntVarArgs() << x << zero << m << one); + } + if ((rm == RM_IMP) ? b.assigned() : + (!b.assigned() || (b.val() != 1))) + return arithmetic_failed + ("reified truth: x, zero, m, one", + Gecode::IntVarArgs() << x << zero << m << one); + if ((x.min() != -7) || (x.max() != 7)) + return arithmetic_failed + ("unconstrained factor: x, zero, m, one", + Gecode::IntVarArgs() << x << zero << m << one); + } + return true; + } + }; + + /// Known identities at the integer limits, independent of the oracles. + class NumberTheoryLimits : public ::Test::Base { + class TestSpace : public Gecode::Space { + public: + virtual Gecode::Space* copy(void) { return nullptr; } + }; + public: + NumberTheoryLimits(void) + : ::Test::Base("Int::Arithmetic::NumberTheoryLimits") {} + virtual bool run(void) { + using namespace Gecode; + const int hi=Gecode::Int::Limits::max; + const int lo=Gecode::Int::Limits::min; + const int cases[][4] = { + {lo,hi,hi,1}, {hi,lo,hi,1}, {lo,0,hi,1}, + {0,lo,hi,0}, {0,0,0,1}, {lo,1,1,0} + }; + for (const auto& c : cases) { + TestSpace home; + IntVar x(home,c[0],c[0]), y(home,c[1],c[1]), g(home,0,hi); + BoolVar b(home,0,1), correct(home,0,1), wrong(home,0,1); + gcd(home,x,y,g); + divides(home,x,y,Reify(b)); + IntVar expected(home,c[2],c[2]); + IntVar other(home,c[2] == 0 ? 1 : 0,c[2] == 0 ? 1 : 0); + gcd(home,x,y,expected,Reify(correct)); + gcd(home,x,y,other,Reify(wrong)); + if ((home.status() == SS_FAILED) || !g.assigned() || + (g.val() != c[2]) || !b.assigned() || (b.val() != c[3]) || + !correct.assigned() || (correct.val() != 1) || + !wrong.assigned() || (wrong.val() != 0)) + return arithmetic_failed + ("signed GCD and divisibility: x, y, g, expected, other", + Gecode::IntVarArgs() << x << y << g << expected << other); + } + // hi-1 is -1 modulo hi. Its square is 1 and its cube is -1; + // negating one factor reverses these residues. The cube exceeds + // signed 64-bit range, but its residue remains exactly known. + for (int n=2; n<=3; n++) + for (int negative=0; negative<=1; negative++) + for (int variable=0; variable<=1; variable++) { + const int expected=((n+negative)%2 == 0) ? 1 : hi-1; + { + TestSpace home; + IntVarArgs x(home,n,hi-1,hi-1); + if (negative) x[0]=IntVar(home,1-hi,1-hi); + IntVar m(home,hi-1,hi), y(home,0,hi-1); + if (variable) product_mod(home,x,m,y); + else product_mod(home,x,hi,y); + if (home.status() == SS_FAILED) + return arithmetic_failed + ("large product before rewrite: x, m, y", + Gecode::IntVarArgs() << x << m << y); + // Force variable-modulus rewriting after initial propagation. + rel(home,m,IRT_EQ,hi); + if ((home.status() == SS_FAILED) || !y.assigned() || + (y.val() != expected)) + return arithmetic_failed + ("large product residue: x, m, y", + Gecode::IntVarArgs() << x << m << y); + } + for (ReifyMode rm : {RM_EQV,RM_IMP,RM_PMI}) + for (int truth=0; truth<=1; truth++) + for (int control=0; control<=1; control++) { + TestSpace home; + IntVarArgs x(home,n,hi-1,hi-1); + if (negative) x[0]=IntVar(home,1-hi,1-hi); + IntVar m(home,hi,hi); + IntVar y(home,truth ? expected : 0,truth ? expected : 0); + BoolVar b(home,control,control); + if (variable) product_mod(home,x,m,y,Reify(b,rm)); + else product_mod(home,x,hi,y,Reify(b,rm)); + const bool holds=(rm == RM_EQV) ? (control == truth) : + ((rm == RM_IMP) ? (!control || truth) : (!truth || control)); + if ((home.status() != SS_FAILED) != holds) + return arithmetic_failed + ("large product reification: x, m, y", + Gecode::IntVarArgs() << x << m << y); + } + } + return true; + } + }; + + /// Evaluate an exact product for testing without overflowing. + bool product_value(const Assignment& x, int n, int& product) { + for (int i=0; i(x[i]))) + return false; + p *= static_cast(x[i]); + } + if (!Gecode::Int::Limits::valid(p)) + return false; + product = static_cast(p); + return true; + } + + /// %Test for an ordinary and reified three-factor product + // Zero membership and assigned-product membership can change without + // a bounds event, including when factors or the result are aliased. + class ProductXYZR : public Test { + public: + ProductXYZR(const std::string& s, const Gecode::IntSet& d, + Gecode::IntPropLevel ipl) + : Test("Arithmetic::Product::XYZR::"+str(ipl)+"::"+s, + 4,d,true,ipl) { contest = CTL_NONE; testfix=false; } + virtual bool solution(const Assignment& x) const { + int p; + return product_value(x,3,p) && (p == x[3]); + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { + Gecode::IntVarArgs f({x[0],x[1],x[2]}); + Gecode::product(home,f,x[3],ipl); + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x, + Gecode::Reify r) { + Gecode::IntVarArgs f({x[0],x[1],x[2]}); + Gecode::product(home,f,x[3],r,ipl); + } + }; + + /// %Test for the empty product + class ProductEmpty : public Test { + public: + ProductEmpty(const std::string& s, const Gecode::IntSet& d, + Gecode::IntPropLevel ipl) + : Test("Arithmetic::Product::Empty::"+str(ipl)+"::"+s, + 1,d,true,ipl) { + contest = CTL_NONE; + } + virtual bool solution(const Assignment& x) const { + return x[0] == 1; + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { + Gecode::product(home,Gecode::IntVarArgs(),x[0],ipl); + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x, + Gecode::Reify r) { + Gecode::product(home,Gecode::IntVarArgs(),x[0],r,ipl); + } + }; + + /// %Test for the singleton product + class ProductSingleton : public Test { + public: + ProductSingleton(const std::string& s, const Gecode::IntSet& d, + Gecode::IntPropLevel ipl) + : Test("Arithmetic::Product::Singleton::"+str(ipl)+"::"+s, + 2,d,true,ipl) { contest = CTL_NONE; } + virtual bool solution(const Assignment& x) const { + return x[0] == x[1]; + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { + Gecode::product(home,Gecode::IntVarArgs({x[0]}),x[1],ipl); + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x, + Gecode::Reify r) { + Gecode::product(home,Gecode::IntVarArgs({x[0]}),x[1],r,ipl); + } + }; + + /// %Test repeated factors and result aliasing + class ProductXXYAlias : public Test { + public: + ProductXXYAlias(const std::string& s, const Gecode::IntSet& d, + Gecode::IntPropLevel ipl) + : Test("Arithmetic::Product::XXYAlias::"+str(ipl)+"::"+s, + 2,d,true,ipl) { contest = CTL_NONE; testfix=false; } + virtual bool solution(const Assignment& x) const { + long long int p = static_cast(x[0]) * x[0] * x[1]; + return (p >= Gecode::Int::Limits::min) && + (p <= Gecode::Int::Limits::max) && (p == x[1]); + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { + Gecode::product(home,Gecode::IntVarArgs({x[0],x[0],x[1]}),x[1],ipl); + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x, + Gecode::Reify r) { + Gecode::product(home,Gecode::IntVarArgs({x[0],x[0],x[1]}),x[1],r,ipl); + } + }; + + /// %Test product bounds propagation beyond the support-enumeration cap + class ProductBoundsLarge : public ::Test::Base { + protected: + class TestSpace : public Gecode::Space { + public: + virtual Gecode::Space* copy(void) { return nullptr; } + }; + public: + ProductBoundsLarge(void) + : ::Test::Base("Int::Arithmetic::Product::BoundsLarge") {} + virtual bool run(void) { + using namespace Gecode; + { + TestSpace home; + IntVarArgs x(home,3,2,100); + IntVar y(home,0,Gecode::Int::Limits::max); + BoolVar b(home,1,1); + product(home,x,y,Reify(b)); + if ((home.status() == SS_FAILED) || + (y.min() != 8) || (y.max() != 1000000)) + return arithmetic_failed + ("forward bounds: x, y", + Gecode::IntVarArgs() << x << y); + } + { + TestSpace home; + IntVarArgs x(home,3,10,100); + IntVar y(home,1000,1000); + product(home,x,y); + if (home.status() == SS_FAILED) + return arithmetic_failed + ("inverse bounds failure: x, y", + Gecode::IntVarArgs() << x << y); + for (int i=0; i(home.clone()); + PropagatorGroup::all.disable(*clone); + rel(home,home.x,IRT_NQ,0); + rel(*clone,clone->x,IRT_NQ,0); + (void) home.status(); + (void) clone->status(); + rel(home,home.y,IRT_GQ,2); + rel(*clone,clone->y,IRT_GQ,2); + PropagatorGroup::all.enable(*clone); + const bool ok=(home.status()!=SS_FAILED) && + (clone->status()!=SS_FAILED) && (home.x.min()==1) && + (clone->x.min()==1) && (home.y.min()==2) && + (clone->y.min()==2); + delete clone; + if (!ok) + return arithmetic_failed + ("rescheduling: home.x, home.q, home.y", + Gecode::IntVarArgs() << home.x << home.q << home.y); + } + return true; + } + }; + + /// Test repeated powers and result-alias cancellation + class ProductPowerAlias : public ::Test::Base { + protected: + class TestSpace : public Gecode::Space { + public: + virtual Gecode::Space* copy(void) { return nullptr; } + }; + public: + ProductPowerAlias(void) + : ::Test::Base("Int::Arithmetic::Product::PowerAlias") {} + virtual bool run(void) { + using namespace Gecode; + { + TestSpace home; + IntVar x(home,-10,10), y(home,-10,100); + product(home,IntVarArgs({x,x}),y); + if ((home.status() == SS_FAILED) || (y.min() != 0)) + return arithmetic_failed + ("square sign: x, y", + Gecode::IntVarArgs() << x << y); + } + { + TestSpace home; + IntVar x(home,-10,10), y(home,20,30); + product(home,IntVarArgs({x,x}),y); + if ((home.status() == SS_FAILED) || + (x.min() != -5) || (x.max() != 5)) + return arithmetic_failed + ("square inverse: x, y", + Gecode::IntVarArgs() << x << y); + } + { + TestSpace home; + IntVar x(home,-10,10), y(home,20,30); + product(home,IntVarArgs({x,x,x}),y); + if ((home.status() == SS_FAILED) || !x.assigned() || + (x.val() != 3)) + return arithmetic_failed + ("cube inverse: x, y", + Gecode::IntVarArgs() << x << y); + } + { + TestSpace home; + IntVar x(home,-10,10), y(home,-30,-20); + product(home,IntVarArgs({x,x,x}),y); + if ((home.status() == SS_FAILED) || !x.assigned() || + (x.val() != -3)) + return arithmetic_failed + ("negative cube inverse: x, y", + Gecode::IntVarArgs() << x << y); + } + { + TestSpace home; + IntVar x(home,-10,10), two(home,2,2), y(home,50,72); + product(home,IntVarArgs({x,x,two}),y); + if ((home.status() == SS_FAILED) || + (x.min() != -6) || (x.max() != 6)) + return arithmetic_failed + ("square cofactor: x, two, y", + Gecode::IntVarArgs() << x << two << y); + } + { + // x*y=y has only the zero branch when x cannot be one. + TestSpace home; + IntVar x(home,2,4), y(home,-10,10); + product(home,IntVarArgs({x,y}),y); + if ((home.status() == SS_FAILED) || !y.assigned() || + (y.val() != 0)) + return arithmetic_failed + ("zero result alias: x, y", + Gecode::IntVarArgs() << x << y); + } + { + // A nonzero result permits cancellation of one result occurrence. + TestSpace home; + IntVar x(home,0,2), y(home,2,10); + product(home,IntVarArgs({x,y}),y); + if ((home.status() == SS_FAILED) || !x.assigned() || + (x.val() != 1)) + return arithmetic_failed + ("nonzero cancellation: x, y", + Gecode::IntVarArgs() << x << y); + } + { + // Cancelling one of two result occurrences leaves y*x=1. + TestSpace home; + IntVar x(home,-1,0), y(home,-2,-1); + product(home,IntVarArgs({y,y,x}),y); + if ((home.status() == SS_FAILED) || !x.assigned() || + !y.assigned() || (x.val() != -1) || (y.val() != -1)) + return arithmetic_failed + ("repeated result cancellation: x, y", + Gecode::IntVarArgs() << x << y); + } + { + // Direct n-ary evaluation retains zero after an overflowing prefix. + TestSpace home; + const int hi=Gecode::Int::Limits::max; + IntVar a(home,hi,hi), z(home,0,0), y(home,0,0); + product(home,IntVarArgs({a,a,z}),y); + if (home.status() == SS_FAILED) + return arithmetic_failed + ("overflow followed by zero: a, z, y", + Gecode::IntVarArgs() << a << z << y); + } + { + TestSpace home; + IntVar x(home,-10,10), y(home,20,30); + BoolVar b(home,1,1); + product(home,IntVarArgs({x,x}),y,Reify(b,RM_EQV)); + if ((home.status() == SS_FAILED) || + (x.min() != -5) || (x.max() != 5)) + return arithmetic_failed + ("reified square: x, y", + Gecode::IntVarArgs() << x << y); + } + return true; + } + }; + + /// Test zero-aware inverse bounds for every cofactor sign class + class ProductInverseBounds : public ::Test::Base { + protected: + class TestSpace : public Gecode::Space { + public: + virtual Gecode::Space* copy(void) { return nullptr; } + }; + static bool bounds(const Gecode::IntVar& x, int min, int max) { + return (x.min() == min) && (x.max() == max); + } + public: + ProductInverseBounds(void) + : ::Test::Base("Int::Arithmetic::Product::InverseBounds") {} + virtual bool run(void) { + using namespace Gecode; + struct Case { int qmin; int qmax; int xmin; int xmax; }; + const Case cases[] = { + { 2, 4, 5, 12}, // Positive + { 0, 4, 5, 24}, // Nonnegative + {-4,-2, -12, -5}, // Negative + {-4, 0, -24, -5}, // Nonpositive + {-4, 3, -24, 24} // Mixed across zero + }; + for (unsigned int i=0; i(x[i]) % m; + if (q < 0) + q += m; + p = (p*q) % m; + } + return static_cast(p); + } + + /// %Test for an ordinary and reified three-factor modular product + // Reified status tests observe interior residue membership, including + // the empty-product residue, but subscribe only to bounds changes. + class ProductModXYZR : public Test { + protected: + int m; + public: + ProductModXYZR(const std::string& s, const Gecode::IntSet& d, + int m0, Gecode::IntPropLevel ipl, bool r=true) + : Test("Arithmetic::ProductMod::XYZR::"+str(ipl)+"::"+s+(r ? "" : "::Plain"), + 4,d,r,ipl), m(m0) { contest=CTL_NONE; testfix=!r; } + virtual bool solution(const Assignment& x) const { + return product_mod_value(x,3,m) == x[3]; + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { + Gecode::product_mod(home,Gecode::IntVarArgs({x[0],x[1],x[2]}), + m,x[3],ipl); + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x, + Gecode::Reify r) { + Gecode::product_mod(home,Gecode::IntVarArgs({x[0],x[1],x[2]}), + m,x[3],r,ipl); + } + }; + + /// %Test for the empty modular product, including modulus one + class ProductModEmpty : public Test { + protected: + int m; + public: + ProductModEmpty(const std::string& s, const Gecode::IntSet& d, + int m0, Gecode::IntPropLevel ipl) + : Test("Arithmetic::ProductMod::Empty::"+str(ipl)+"::"+s, + 1,d,true,ipl), m(m0) { contest=CTL_NONE; } + virtual bool solution(const Assignment& x) const { + return x[0] == (1 % m); + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { + Gecode::product_mod(home,Gecode::IntVarArgs(),m,x[0],ipl); + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x, + Gecode::Reify r) { + Gecode::product_mod(home,Gecode::IntVarArgs(),m,x[0],r,ipl); + } + }; + + /// %Test for a singleton modular product + class ProductModSingleton : public Test { + protected: + int m; + public: + ProductModSingleton(const std::string& s, const Gecode::IntSet& d, + int m0, Gecode::IntPropLevel ipl, bool r=true) + : Test("Arithmetic::ProductMod::Singleton::"+str(ipl)+"::"+s+(r ? "" : "::Plain"), + 2,d,r,ipl), m(m0) { contest=CTL_NONE; testfix=!r; } + virtual bool solution(const Assignment& x) const { + int r = x[0] % m; + if (r < 0) + r += m; + return r == x[1]; + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { + Gecode::product_mod(home,Gecode::IntVarArgs({x[0]}),m,x[1],ipl); + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x, + Gecode::Reify r) { + Gecode::product_mod(home,Gecode::IntVarArgs({x[0]}),m,x[1],r,ipl); + } + }; + + /// %Test repeated factors with the result aliased to a factor + class ProductModXXYAlias : public Test { + protected: + int m; + public: + ProductModXXYAlias(const std::string& s, const Gecode::IntSet& d, + int m0, Gecode::IntPropLevel ipl, bool r=true) + : Test("Arithmetic::ProductMod::XXYAlias::"+str(ipl)+"::"+s+(r ? "" : "::Plain"), + 2,d,r,ipl), m(m0) { contest=CTL_NONE; testfix=!r; } + virtual bool solution(const Assignment& x) const { + long long int a = x[0] % m; + long long int b = x[1] % m; + if (a < 0) a += m; + if (b < 0) b += m; + return (((a*a) % m)*b) % m == x[1]; + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { + Gecode::product_mod(home,Gecode::IntVarArgs({x[0],x[0],x[1]}), + m,x[1],ipl); + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x, + Gecode::Reify r) { + Gecode::product_mod(home,Gecode::IntVarArgs({x[0],x[0],x[1]}), + m,x[1],r,ipl); + } + }; + + /// Targeted algebraic and quotient-band propagation for fixed modulus + class ProductModAlgebraic : public ::Test::Base { + protected: + class TestSpace : public Gecode::Space { + public: + virtual Gecode::Space* copy(void) { return nullptr; } + }; + public: + ProductModAlgebraic(void) + : ::Test::Base("Int::Arithmetic::ProductMod::Algebraic") {} + virtual bool run(void) { + using namespace Gecode; + { + TestSpace home; + IntVarArgs x(home,3,-100,100); + IntVar y(home,-10,10); + product_mod(home,x,1,y); + if ((home.status() == SS_FAILED) || !y.assigned() || + (y.val() != 0)) + return arithmetic_failed + ("modulus one: x, y", + Gecode::IntVarArgs() << x << y); + } + { + TestSpace home; + IntVar z(home,14,14), x(home,-1000,1000), y(home,0,6); + product_mod(home,IntVarArgs({z,x}),7,y); + if ((home.status() == SS_FAILED) || !y.assigned() || + (y.val() != 0)) + return arithmetic_failed + ("zero residue factor: z, x, y", + Gecode::IntVarArgs() << z << x << y); + } + { + TestSpace home; + IntVar one(home,1,1), x(home,20,30), y(home,0,99); + product_mod(home,IntVarArgs({one,x}),100,y); + if ((home.status() == SS_FAILED) || + (y.min() != 20) || (y.max() != 30)) + return arithmetic_failed + ("nonwrapping product: one, x, y", + Gecode::IntVarArgs() << one << x << y); + } + { + TestSpace home; + IntVar x(home,15,19), y(home,0,6); + product_mod(home,IntVarArgs({x}),7,y); + if ((home.status() == SS_FAILED) || + (y.min() != 1) || (y.max() != 5)) + return arithmetic_failed + ("positive quotient band: x, y", + Gecode::IntVarArgs() << x << y); + } + { + TestSpace home; + IntVar x(home,-20,-16), y(home,0,6); + product_mod(home,IntVarArgs({x}),7,y); + if ((home.status() == SS_FAILED) || + (y.min() != 1) || (y.max() != 5)) + return arithmetic_failed + ("negative quotient band: x, y", + Gecode::IntVarArgs() << x << y); + } + { + TestSpace home; + IntVar c(home,6,6), x(home,-100,100), y(home,9,9); + product_mod(home,IntVarArgs({c,x}),15,y); + if ((home.status() == SS_FAILED) || + (x.min() != -96) || (x.max() != 99)) + return arithmetic_failed + ("linear congruence: c, x, y", + Gecode::IntVarArgs() << c << x << y); + } + { + TestSpace home; + IntVar c(home,6,6), x(home,-100,100), y(home,8,8); + product_mod(home,IntVarArgs({c,x}),15,y); + if (home.status() != SS_FAILED) + return arithmetic_failed + ("inconsistent congruence: c, x, y", + Gecode::IntVarArgs() << c << x << y); + } + { + TestSpace home; + IntVar z(home,0,0), x(home,-100,100), y(home,0,1); + BoolVar b(home,0,1); + product_mod(home,IntVarArgs({z,x}),7,y,Reify(b)); + rel(home,y,IRT_EQ,0); + if ((home.status() == SS_FAILED) || !b.assigned() || + (b.val() != 1)) + return arithmetic_failed + ("reified zero residue: z, x, y", + Gecode::IntVarArgs() << z << x << y); + } + return true; + } + }; + + /// %Test for an ordinary two-factor product with a variable modulus + // Empty-product status observes membership of 0/1 in the result, and + // divisor status observes interior modulus membership. Reposting can + // therefore strengthen propagation without a subscribed bounds event. + class ProductModVarXYMR : public Test { + public: + ProductModVarXYMR(const std::string& s, const Gecode::IntSet& d, + Gecode::IntPropLevel ipl) + : Test("Arithmetic::ProductModVar::XYMR::"+str(ipl)+"::"+s, + 4,d,true,ipl) { contest=CTL_NONE; testfix=false; } + virtual bool solution(const Assignment& x) const { + return (x[2] > 0) && + (product_mod_value(x,2,x[2]) == x[3]); + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { + Gecode::product_mod(home,Gecode::IntVarArgs({x[0],x[1]}), + x[2],x[3],ipl); + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x, + Gecode::Reify r) { + Gecode::product_mod(home,Gecode::IntVarArgs({x[0],x[1]}), + x[2],x[3],r,ipl); + } + }; + + /// %Test for the empty product with a variable modulus + class ProductModVarEmpty : public Test { + public: + ProductModVarEmpty(const std::string& s, const Gecode::IntSet& d, + Gecode::IntPropLevel ipl) + : Test("Arithmetic::ProductModVar::Empty::"+str(ipl)+"::"+s, + 2,d,true,ipl) { contest=CTL_NONE; testfix=false; } + virtual bool solution(const Assignment& x) const { + return (x[0] > 0) && (x[1] == (1 % x[0])); + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { + Gecode::product_mod(home,Gecode::IntVarArgs(),x[0],x[1],ipl); + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x, + Gecode::Reify r) { + Gecode::product_mod(home,Gecode::IntVarArgs(),x[0],x[1],r,ipl); + } + }; + + /// %Test for a singleton product with a variable modulus + class ProductModVarSingleton : public Test { + public: + ProductModVarSingleton(const std::string& s, const Gecode::IntSet& d, + Gecode::IntPropLevel ipl) + : Test("Arithmetic::ProductModVar::Singleton::"+str(ipl)+"::"+s, + 3,d,true,ipl) { contest=CTL_NONE; testfix=false; } + virtual bool solution(const Assignment& x) const { + if (x[1] <= 0) + return false; + int r = x[0] % x[1]; + if (r < 0) + r += x[1]; + return r == x[2]; + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { + Gecode::product_mod(home,Gecode::IntVarArgs({x[0]}), + x[1],x[2],ipl); + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x, + Gecode::Reify r) { + Gecode::product_mod(home,Gecode::IntVarArgs({x[0]}), + x[1],x[2],r,ipl); + } + }; + + /// %Test repeated factors with a variable modulus + class ProductModVarRepeated : public Test { + public: + ProductModVarRepeated(const std::string& s, const Gecode::IntSet& d, + Gecode::IntPropLevel ipl) + : Test("Arithmetic::ProductModVar::Repeated::"+str(ipl)+"::"+s, + 3,d,true,ipl) { contest=CTL_NONE; testfix=false; } + virtual bool solution(const Assignment& x) const { + if (x[1] <= 0) + return false; + long long int a = x[0] % x[1]; + if (a < 0) a += x[1]; + return (a*a) % x[1] == x[2]; + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { + Gecode::product_mod(home,Gecode::IntVarArgs({x[0],x[0]}), + x[1],x[2],ipl); + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x, + Gecode::Reify r) { + Gecode::product_mod(home,Gecode::IntVarArgs({x[0],x[0]}), + x[1],x[2],r,ipl); + } + }; + + /// %Test modulus/factor aliasing + class ProductModVarModFactorAlias : public Test { + public: + ProductModVarModFactorAlias(const std::string& s, + const Gecode::IntSet& d, + Gecode::IntPropLevel ipl) + : Test("Arithmetic::ProductModVar::ModFactorAlias::"+ + str(ipl)+"::"+s,2,d,true,ipl) { contest=CTL_NONE; testfix=false; } + virtual bool solution(const Assignment& x) const { + return (x[0] > 0) && (x[1] == 0); + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { + Gecode::product_mod(home,Gecode::IntVarArgs({x[0]}), + x[0],x[1],ipl); + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x, + Gecode::Reify r) { + Gecode::product_mod(home,Gecode::IntVarArgs({x[0]}), + x[0],x[1],r,ipl); + } + }; + + /// %Test factor/result aliasing + class ProductModVarFactorResultAlias : public Test { + public: + ProductModVarFactorResultAlias(const std::string& s, + const Gecode::IntSet& d, + Gecode::IntPropLevel ipl) + : Test("Arithmetic::ProductModVar::FactorResultAlias::"+ + str(ipl)+"::"+s,2,d,true,ipl) { contest=CTL_NONE; testfix=false; } + virtual bool solution(const Assignment& x) const { + if (x[1] <= 0) + return false; + int r = x[0] % x[1]; + if (r < 0) r += x[1]; + return r == x[0]; + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { + Gecode::product_mod(home,Gecode::IntVarArgs({x[0]}), + x[1],x[0],ipl); + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x, + Gecode::Reify r) { + Gecode::product_mod(home,Gecode::IntVarArgs({x[0]}), + x[1],x[0],r,ipl); + } + }; + + /// %Test modulus/result aliasing, which is necessarily inconsistent + class ProductModVarModResultAlias : public Test { + public: + ProductModVarModResultAlias(const std::string& s, + const Gecode::IntSet& d, + Gecode::IntPropLevel ipl) + : Test("Arithmetic::ProductModVar::ModResultAlias::"+ + str(ipl)+"::"+s,2,d,true,ipl) { contest=CTL_NONE; } + virtual bool solution(const Assignment&) const { return false; } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { + Gecode::product_mod(home,Gecode::IntVarArgs({x[0]}), + x[1],x[1],ipl); + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x, + Gecode::Reify r) { + Gecode::product_mod(home,Gecode::IntVarArgs({x[0]}), + x[1],x[1],r,ipl); + } + }; + + /// Targeted bounds and arithmetic checks for a variable modulus + class ProductModVarBounds : public ::Test::Base { + protected: + class TestSpace : public Gecode::Space { + public: + virtual Gecode::Space* copy(void) { return nullptr; } + }; + public: + ProductModVarBounds(void) + : ::Test::Base("Int::Arithmetic::ProductModVar::Bounds") {} + virtual bool run(void) { + using namespace Gecode; + { + TestSpace home; + const int vm[5] = {-7,0,1,4,9}; + IntVar m(home,IntSet(vm,5)); + IntVar y(home,-5,12); + product_mod(home,IntVarArgs(),m,y); + if ((home.status() == SS_FAILED) || (m.min() != 1) || + (y.min() != 0) || (y.max() >= m.max())) + return arithmetic_failed + ("empty product range: m, y", + Gecode::IntVarArgs() << m << y); + } + { + TestSpace home; + IntVar x(home,Gecode::Int::Limits::max,Gecode::Int::Limits::max); + IntVar m(home,Gecode::Int::Limits::max,Gecode::Int::Limits::max); + IntVar y(home,0,0); + product_mod(home,IntVarArgs({x}),m,y); + if (home.status() == SS_FAILED) + return arithmetic_failed + ("integer limit modulus: x, m, y", + Gecode::IntVarArgs() << x << m << y); + } + { + TestSpace home; + IntVar x(home,3,5), c(home,2,2), m(home,5,6), y(home,2,2); + product_mod(home,IntVarArgs({x,c}),m,y); + if (home.status() == SS_FAILED) + return arithmetic_failed + ("delayed factor bounds: x, c, m, y", + Gecode::IntVarArgs() << x << c << m << y); + rel(home,x,IRT_LQ,4); + if ((home.status() == SS_FAILED) || !x.assigned() || + (x.val() != 4)) + return arithmetic_failed + ("nonwrapping variable modulus: x, c, m, y", + Gecode::IntVarArgs() << x << c << m << y); + } + { + TestSpace home; + IntVar x(home,10,20), c(home,10,20), m(home,401,509); + IntVar y(home,0,600); + product_mod(home,IntVarArgs({x,c}),m,y); + if ((home.status() == SS_FAILED) || (y.min() != 100) || + (y.max() != 400)) + return arithmetic_failed + ("check 5: x, c, m, y", + Gecode::IntVarArgs() << x << c << m << y); + } + return true; + } + }; + + /// Targeted algebraic propagation for a variable modulus + class ProductModVarAlgebraic : public ::Test::Base { + protected: + class TestSpace : public Gecode::Space { + public: + virtual Gecode::Space* copy(void) { return nullptr; } + }; + static bool domain(const Gecode::IntVar& x, const int* v, int n) { + if (x.size() != static_cast(n)) + return false; + for (int i=0; i(n)) + return false; + for (int i=0; i=-1; m--) { + try { + Gecode::product_mod(home,Gecode::IntVarArgs(),m,y); + return arithmetic_failed + ("accepted invalid modulus: y", + Gecode::IntVarArgs() << y); + } catch (const Gecode::Int::OutOfLimits&) { + } + } + return true; + } + }; + /// %Test for multiplication constraint class MultXYZ : public Test { public: @@ -408,7 +1924,7 @@ namespace Test { namespace Int { : Test("Arithmetic::DivMod::"+s,4,d) {} /// %Test whether \a x is solution virtual bool solution(const Assignment& x) const { - return x[0] == x[1]*x[2]+x[3] && + return x[0] == static_cast(x[1])*x[2]+x[3] && abs(x[3]) < abs(x[1]) && (x[3] == 0 || sgn(x[3]) == sgn(x[0])); } @@ -456,7 +1972,7 @@ namespace Test { namespace Int { divsign * static_cast(floor(static_cast(std::abs(x[0]))/ static_cast(std::abs(x[1])))); - return x[0] == x[1]*divresult+x[2]; + return x[0] == static_cast(x[1])*divresult+x[2]; } /// Post constraint on \a x virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { @@ -1085,6 +2601,13 @@ namespace Test { namespace Int { Gecode::IntSet b(vb,9); Gecode::IntSet c(-8,8); Gecode::IntSet d(-70,70); + const int vg[7] = {-12,-6,-1,0,4,9,12}; + Gecode::IntSet g(vg,7); + const int vp[5] = { + Gecode::Int::Limits::min,-1,0,1,Gecode::Int::Limits::max + }; + Gecode::IntSet p(vp,5); + Gecode::IntSet q(-2,2); (void) new DivMod("A",a); (void) new DivMod("B",b); @@ -1100,6 +2623,47 @@ namespace Test { namespace Int { for (IntPropLevels ipls; ipls(); ++ipls) { + (void) new GcdXYZ("C",c,ipls.ipl()); + (void) new GcdXYZ("C",c,ipls.ipl(),false); + (void) new GcdXYZ("Sparse",g,ipls.ipl()); + (void) new GcdXYZ("Sparse",g,ipls.ipl(),false); + (void) new GcdXXY("C",c,ipls.ipl()); + (void) new GcdXXY("C",c,ipls.ipl(),false); + (void) new GcdXYX("C",c,ipls.ipl()); + (void) new GcdXYX("C",c,ipls.ipl(),false); + (void) new GcdXXX("C",c,ipls.ipl()); + + (void) new DividesXY("C",c,ipls.ipl()); + (void) new DividesXY("Sparse",g,ipls.ipl()); + (void) new DividesXX("C",c,ipls.ipl()); + + (void) new ProductXYZR("C",q,ipls.ipl()); + (void) new ProductXYZR("Sparse",g,ipls.ipl()); + (void) new ProductXYZR("Limits",p,ipls.ipl()); + (void) new ProductEmpty("C",q,ipls.ipl()); + (void) new ProductSingleton("C",q,ipls.ipl()); + (void) new ProductXXYAlias("C",q,ipls.ipl()); + + (void) new ProductModXYZR("C",q,5,ipls.ipl()); + (void) new ProductModXYZR("C",q,5,ipls.ipl(),false); + (void) new ProductModXYZR("Sparse",g,7,ipls.ipl()); + (void) new ProductModXYZR("Sparse",g,7,ipls.ipl(),false); + (void) new ProductModEmpty("C",q,5,ipls.ipl()); + (void) new ProductModEmpty("ModulusOne",q,1,ipls.ipl()); + (void) new ProductModSingleton("C",g,5,ipls.ipl()); + (void) new ProductModSingleton("C",g,5,ipls.ipl(),false); + (void) new ProductModXXYAlias("C",q,5,ipls.ipl()); + (void) new ProductModXXYAlias("C",q,5,ipls.ipl(),false); + + (void) new ProductModVarXYMR("C",q,ipls.ipl()); + (void) new ProductModVarXYMR("Sparse",g,ipls.ipl()); + (void) new ProductModVarEmpty("C",q,ipls.ipl()); + (void) new ProductModVarSingleton("C",q,ipls.ipl()); + (void) new ProductModVarRepeated("C",q,ipls.ipl()); + (void) new ProductModVarModFactorAlias("C",q,ipls.ipl()); + (void) new ProductModVarFactorResultAlias("C",q,ipls.ipl()); + (void) new ProductModVarModResultAlias("C",q,ipls.ipl()); + (void) new AbsXY("A",a,ipls.ipl()); (void) new AbsXY("B",b,ipls.ipl()); (void) new AbsXY("C",c,ipls.ipl()); @@ -1249,6 +2813,19 @@ namespace Test { namespace Int { (void) new ArgMinBool(i,1,false); (void) new ArgMinBoolShared(i,false); } + (void) new NumberTheorySparseBounds; + (void) new NumberTheoryLifecycle; + (void) new NumberTheoryLimits; + (void) new NumberTheoryIdentities; + (void) new ProductModInvalidModulus; + (void) new ProductModAlgebraic; + (void) new ProductModVarBounds; + (void) new ProductModVarAlgebraic; + (void) new ProductModVarInactive; + (void) new ProductBoundsLarge; + (void) new ProductSimplifySign; + (void) new ProductPowerAlias; + (void) new ProductInverseBounds; } };