diff --git a/Makefile.dep b/Makefile.dep index 1f26262fd6..41da67a38a 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 \ @@ -19260,6 +19264,20 @@ test/int/channel$(OBJSUFFIX) test/int/channel$(SBJSUFFIX): \ ./test/int.hh ./test/int.hpp ./test/test.hh \ ./test/test.hpp test/int/arithmetic$(OBJSUFFIX) test/int/arithmetic$(SBJSUFFIX): \ + ./gecode/int/arithmetic.hh ./gecode/int/arithmetic/abs.hpp ./gecode/int/arithmetic/argmax.hpp \ + ./gecode/int/arithmetic/divides.hpp ./gecode/int/arithmetic/divmod.hpp ./gecode/int/arithmetic/gcd.hpp \ + ./gecode/int/arithmetic/max.hpp ./gecode/int/arithmetic/mult.hpp ./gecode/int/arithmetic/nroot.hpp \ + ./gecode/int/arithmetic/pow-ops.hpp ./gecode/int/arithmetic/pow.hpp ./gecode/int/arithmetic/product-mod.hpp \ + ./gecode/int/arithmetic/product.hpp ./gecode/int/bool.hh ./gecode/int/bool/base.hpp \ + ./gecode/int/bool/clause.hpp ./gecode/int/bool/eq.hpp ./gecode/int/bool/eqv.hpp \ + ./gecode/int/bool/ite.hpp ./gecode/int/bool/lq.hpp ./gecode/int/bool/or.hpp \ + ./gecode/int/idx-view.hh ./gecode/int/idx-view.hpp ./gecode/int/linear.hh \ + ./gecode/int/linear/bool-int.hpp ./gecode/int/linear/bool-scale.hpp ./gecode/int/linear/bool-view.hpp \ + ./gecode/int/linear/int-bin.hpp ./gecode/int/linear/int-dom.hpp ./gecode/int/linear/int-nary.hpp \ + ./gecode/int/linear/int-noview.hpp ./gecode/int/linear/int-ter.hpp ./gecode/int/linear/post.hpp \ + ./gecode/int/rel.hh ./gecode/int/rel/eq.hpp ./gecode/int/rel/lex.hpp \ + ./gecode/int/rel/lq-le.hpp ./gecode/int/rel/nq.hpp ./gecode/int/support-values.hh \ + ./gecode/int/support-values.hpp \ ./gecode/float.hh ./gecode/float/array-traits.hpp ./gecode/float/array.hpp \ ./gecode/float/branch.hpp ./gecode/float/branch/action.hpp ./gecode/float/branch/afc.hpp \ ./gecode/float/branch/assign.hpp ./gecode/float/branch/chb.hpp ./gecode/float/branch/traits.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..555d4ddd8c --- /dev/null +++ b/docs/integer-number-theory.md @@ -0,0 +1,96 @@ +# 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. + +## Propagation and cost + +For a fixed positive modulus `m`, assigned factors are reduced to a coefficient +`c`. Every result is a multiple of `gcd(c,m)`, even with several unassigned +factors. A zero coefficient fixes the result to zero. Thus `6*x*z mod 12 = 1` +fails without branching, and assigned factors 2 and 3 force a zero residue +modulo 6. With one free factor and an assigned result, a linear congruence +also restricts the factor bounds. + +A nonzero dividend bounds the magnitude of every possible divisor, including +when the dividend is not assigned. Negating `divides(0,x)` removes zero from +`x`. For GCD, an assigned nonzero operand bounds the result to its magnitude; +a fixed quotient band can identify its only possible divisor. When both an +operand and the GCD are assigned, endpoint probes check that the remaining +quotient is coprime. For example, `gcd(6,x)=3` with `x=6..12` assigns `x=9`. +The probes inspect at most 64 arithmetic candidates per endpoint and retain +the original bound if that budget is exhausted. They use the interval hull, +so interior holes do not change a bounds-subscribed actor's fixpoint. + +Exact-product posting sorts equal views together. The grouping survives +cloning, and subsequent propagation groups adjacent occurrences in linear +time. Assigned units are compacted in place while preserving that order. +Fixed- and variable-modulus products use prefix/suffix intervals to +compute all cofactors in linear time per filtering pass, with linear temporary +storage. The temporary intervals are rebuilt after bounds change. Repeated +factors in exact products additionally require integer-root bounds; their +arithmetic work is at most proportional to the factor count times the number +of bits in the integer limit. Local closure can require several passes, +especially across sparse domain endpoints. + +Modular interval reasoning uses exact representable endpoints. If an endpoint +overflows, that interval deduction is skipped; saturation is not valid for +modular arithmetic. Assigned modular products remain safe because each +multiplication is reduced before the next factor. Cheap coefficient and +zero checks precede allocation and cofactor filtering. Fixed-modulus actors +fold assigned occurrences into a persistent residue coefficient and remove +them from the active view array. They also retain the exact coefficient while +it fits in the integer limits, for interval reasoning; an unavailable exact +coefficient never replaces the modular one. Later scans and clones visit only +the remaining occurrences. Variable-modulus actors remove assigned ones, but +retain other constants until the modulus is fixed. Both forms exclude zero +endpoints from factors when the result is strictly positive, including factors +with negative domains. + +These are deliberately incomplete rules. They do not enumerate all divisors +or all supported residues. Undecided or negated reified constraints generally +perform less filtering than asserted relations. The variable-modulus interval +stage handles nonnegative factors in a single quotient band; fixing the +modulus enables the fixed-modulus rules. Bounds subscriptions can miss +opportunities caused solely by interior domain removal in the propagators +that declare weak monotonicity. Fully assigned tuples are still checked +exactly. Faster propagation passes alone do not establish faster solving; +compare search and end-to-end time on representative models. 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..6e24584070 100644 --- a/gecode/int/arithmetic.hh +++ b/gecode/int/arithmetic.hh @@ -766,6 +766,322 @@ 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; + /// Canonical residue of removed factors + int c; + /// Exact product of removed factors, or zero when not representable + int a; + /// Constructor for posting + ProductMod(Home home, ViewArray& x, int m, IntView y, + int c, int a); + /// 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); + /// Delete propagator and return its size + virtual size_t dispose(Space& home); + }; + + /** \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 +1173,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..50749e3ed8 --- /dev/null +++ b/gecode/int/arithmetic/divides.hpp @@ -0,0 +1,223 @@ +/* -*- 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)); + } + + // Every divisor of a nonzero dividend has no greater magnitude. + if (!dividend.in(0)) { + const int a=gcd_max_abs(dividend); + GECODE_ME_CHECK(divisor.gq(home,-a)); + GECODE_ME_CHECK(divisor.lq(home,a)); + GECODE_ME_CHECK(divisor.nq(home,0)); + } + + 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.in(0)) { + const int a=gcd_max_abs(x1); + 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); + if (x0.assigned() && (x0.val() == 0)) { + GECODE_ME_CHECK(x1.nq(home,0)); + 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..169880d147 --- /dev/null +++ b/gecode/int/arithmetic/gcd.hpp @@ -0,0 +1,415 @@ +/* -*- 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; + } + + /// Probe at most 64 multiples for a coprime quotient at an endpoint. + template + inline bool + gcd_support_min(View x, int a, int g, int& bound) { + bound=x.min(); + long long int v=ceil_div_xx(static_cast(x.min()), + static_cast(g))*g; + // Use the interval hull: interior holes must not affect a BND actor's + // fixpoint. The candidates are isolated arithmetic supports, not a + // domain enumeration. On exhaustion retain the original endpoint. + for (int budget=64; (budget > 0) && (v <= x.max()); --budget, v+=g) + if (gcd_value(a,static_cast(v)) == g) { + bound=static_cast(v); + return true; + } + return v <= x.max(); + } + + /// With an assigned operand, exclude endpoints lacking a coprime quotient. + inline ExecStatus + gcd_coprime_bounds(Home home, IntView x, int a, int g) { + if ((a % g) != 0) + return ES_FAILED; + int l, negative_u; + if (!gcd_support_min(x,a,g,l) || + !gcd_support_min(MinusView(x),a,g,negative_u)) + return ES_FAILED; + // On budget exhaustion retain the original endpoint: never advance a + // partially checked bound, which could cause an unbounded closure walk. + GECODE_ME_CHECK(x.gq(home,l)); + GECODE_ME_CHECK(x.lq(home,-negative_u)); + return ES_OK; + } + + /// A fixed quotient of an assigned operand admits at most one GCD divisor. + inline ExecStatus + gcd_divisor_bounds(Home home, IntView g, int a) { + a=(a < 0) ? -a : a; + assert(a > 0); + GECODE_ME_CHECK(g.gq(home,1)); + GECODE_ME_CHECK(g.lq(home,a)); + const int q0=a/g.max(), q1=a/g.min(); + if (q0 == q1) { + if ((a % q0) != 0) + return ES_FAILED; + GECODE_ME_CHECK(g.eq(home,a/q0)); + } + 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))); + // Result narrowing can assign an aliased operand to zero. + if (x0.assigned() && (x0.val() != 0)) + GECODE_ES_CHECK(gcd_divisor_bounds(home,x2,x0.val())); + if (x1.assigned() && (x1.val() != 0)) + GECODE_ES_CHECK(gcd_divisor_bounds(home,x2,x1.val())); + 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()) + GECODE_ES_CHECK(gcd_coprime_bounds(home,x1,x0.val(),g)); + if (x1.assigned()) + GECODE_ES_CHECK(gcd_coprime_bounds(home,x0,x1.val(),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))); + // Result narrowing can assign an aliased operand to zero. + if (x0.assigned() && (x0.val() != 0)) + GECODE_ES_CHECK(gcd_divisor_bounds(home,x2,x0.val())); + if (x1.assigned() && (x1.val() != 0)) + GECODE_ES_CHECK(gcd_divisor_bounds(home,x2,x1.val())); + 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()) + GECODE_ES_CHECK(gcd_coprime_bounds(home,x1,x0.val(),g)); + if (x1.assigned()) + GECODE_ES_CHECK(gcd_coprime_bounds(home,x0,x1.val(),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..9269bdbfb7 --- /dev/null +++ b/gecode/int/arithmetic/product-mod.hpp @@ -0,0 +1,907 @@ +/* -*- 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 { + + /// Multiply intervals exactly, rejecting overflow rather than clipping. + inline bool + product_mod_interval_mul(ProductInterval a, ProductInterval b, + ProductInterval& p) { + long long int q[4] = { + static_cast(a.min)*b.min, + static_cast(a.min)*b.max, + static_cast(a.max)*b.min, + static_cast(a.max)*b.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 exact representable product bounds, or false on overflow. + inline bool + product_mod_interval(const ViewArray& x, ProductInterval& p, + int a=1) { + p.min=a; p.max=a; + for (int i=0; i& x, + ProductInterval* prefix, ProductInterval* suffix, + int a=1) { + prefix[0]=ProductInterval{a,a}; + suffix[x.size()]=ProductInterval{1,1}; + for (int i=0; i& x, int m, + int& free_count) { + long long int p=1 % m; + free_count=0; + for (int i=0; i(x[i].val()) % m; + if (q < 0) q+=m; + p=(p*q) % m; + } + } + return static_cast(p); + } + + /// Fold assigned occurrences into modular and representable exact coefficients. + inline void + product_mod_fold(ViewArray& x, int m, int& c, int& a) { + for (int i=x.size(); i--;) + if (x[i].assigned()) { + const int v=x[i].val(); + c=static_cast((static_cast(c)*v) % m); + if (c < 0) c+=m; + const long long int p=static_cast(a)*v; + // Zero marks an unavailable exact coefficient; c remains exact modulo m. + a=Limits::valid(p) ? static_cast(p) : 0; + x.move_lst(i); + } + } + + /// Drop assigned units without changing a product with an unknown modulus. + inline void + product_mod_units(ViewArray& x) { + for (int i=x.size(); i--;) + if (x[i].assigned() && (x[i].val() == 1)) + x.move_lst(i); + } + + /// A nonzero residue excludes bounds-visible zero factors. + inline ExecStatus + product_mod_nonzero(Space& home, ViewArray& x, bool& modified) { + for (int i=0; i& x, const IntView& y, int m) { + if ((y.max() < 0) || (y.min() >= m)) + return RT_FALSE; + int free_count; + const int c=product_mod_coefficient(x,m,free_count); + if ((c == 0) || (free_count == 0)) { + if (!y.in(c)) return RT_FALSE; + return y.assigned() ? RT_TRUE : RT_MAYBE; + } + const int g=gcd_value(c,m); + const long long int l=ceil_div_xx + (static_cast(std::max(0,y.min())), + static_cast(g))*g; + if (l > std::min(m-1,y.max())) + return RT_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; + 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; + } + + /// Extreme solutions to c*x = y modulo m, or an empty interval. + inline ProductInterval + product_mod_congruence(ProductInterval x, int c, int m, int y) { + const int g=gcd_value(c,m); + if ((y % g) != 0) + return ProductInterval{1,0}; + const long long int step=m/g; + long long int r=0; + if (step > 1) { + const long long int a=c/g, b=y/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.min)-r,step)*step; + const long long int u=r + floor_div_xx + (static_cast(x.max)-r,step)*step; + if (l > u) + return ProductInterval{1,0}; + // Nonempty bounds lie inside x, so narrowing to int is safe. + return ProductInterval{static_cast(l),static_cast(u)}; + } + + forceinline + ProductMod::ProductMod(Home home, ViewArray& z, int modulus, + IntView w, int c0, int a0) + : NaryOnePropagator(home,z,w), + m(modulus), c(c0), a(a0) {} + + 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; + } + int c=1, a=1; + product_mod_fold(x,m,c,a); + if ((c == 0) || (x.size() == 0)) { + GECODE_ME_CHECK(y.eq(home,c)); + return ES_OK; + } + (void) new (home) ProductMod(home,x,m,y,c,a); + return ES_OK; + } + + forceinline + ProductMod::ProductMod(Space& home, ProductMod& p) + : NaryOnePropagator(home,p), + m(p.m), c(p.c), a(p.a) {} + + forceinline Actor* + ProductMod::copy(Space& home) { + return new (home) ProductMod(home,*this); + } + + forceinline PropCost + ProductMod::cost(const Space&, const ModEventDelta&) const { + return PropCost::linear(PropCost::HI,x.size()+1); + } + + forceinline size_t + ProductMod::dispose(Space& home) { + (void) NaryOnePropagator::dispose(home); + return sizeof(*this); + } + + inline ExecStatus + ProductMod::propagate(Space& home, const ModEventDelta&) { + GECODE_ME_CHECK(y.gq(home,0)); + GECODE_ME_CHECK(y.lq(home,m-1)); + + Region region; + ProductInterval* prefix=NULL; + ProductInterval* suffix=NULL; + bool modified; + do { + modified=false; + product_mod_fold(x,m,c,a); + if ((c == 0) || (x.size() == 0)) { + GECODE_ME_CHECK(y.eq(home,c)); + return home.ES_SUBSUMED(*this); + } + const unsigned int old_size=y.size(); + GECODE_ES_CHECK(gcd_multiple_bounds(home,y,gcd_value(c,m))); + modified |= old_size != y.size(); + if (y.min() > 0) + GECODE_ES_CHECK(product_mod_nonzero(home,x,modified)); + + // One remaining occurrence gives a linear congruence. + if ((x.size() == 1) && y.assigned()) { + const ProductInterval bounds={x[0].min(),x[0].max()}; + const ProductInterval solutions= + product_mod_congruence(bounds,c,m,y.val()); + if (solutions.min > solutions.max) + return ES_FAILED; + ModEvent me=x[0].gq(home,solutions.min); + if (me_failed(me)) return ES_FAILED; + modified |= me_modified(me); + me=x[0].lq(home,solutions.max); + if (me_failed(me)) return ES_FAILED; + modified |= me_modified(me); + } + ProductInterval p; + if ((a != 0) && product_mod_interval(x,p,a)) { + 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); + } + + // Allocate cofactors only when interval inversion can contribute. + if (prefix == NULL) { + prefix=region.alloc(x.size()+1); + suffix=region.alloc(x.size()+1); + } + if (!product_mod_intervals(x,prefix,suffix,a)) + continue; + + // 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))) { + const ProductWideInterval d= + product_quotient_bounds(tmin,tmax,q.min,q.max); + const long long int l=std::max + (d.min,static_cast(Limits::min)); + const long long int u=std::min + (d.max,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); + + return 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) { + Region region; + ProductInterval* prefix=NULL; + ProductInterval* suffix=NULL; + bool modified; + do { + modified=false; + if (y.min() > 0) + GECODE_ES_CHECK(product_mod_nonzero(home,x,modified)); + for (int i=0; i(p.min) / m.max(); + const long long int kmax = + static_cast(p.max) / m.min(); + if (kmin != kmax) + continue; + + 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); + } + + if (prefix == NULL) { + prefix=region.alloc(x.size()+1); + suffix=region.alloc(x.size()+1); + } + if (!product_mod_intervals(x,prefix,suffix)) + continue; + + 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); + } + } + } while (modified); + 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; + if (m.assigned()) + return product_mod_status(x,y,m.val()); + bool zero=product_mod_var_mod_factor(x,m); + 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() >= 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; + } + } + return RT_MAYBE; + } + + 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)); + product_mod_units(x); + 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::linear(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)); + product_mod_units(x); + 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)); + + GECODE_ES_CHECK(product_mod_var_ranges(home,x,m,y)); + 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. + const unsigned int modulus_size=m.size(); + 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)); + // A sparse-domain update can expose a new fixed quotient band. + return (modulus_size != m.size()) ? ES_NOFIX : ES_FIX; + } + + forceinline size_t + ProductModVar::dispose(Space& home) { + x.cancel(home,*this,PC_INT_BND); + m.cancel(home,*this,PC_INT_BND); + y.cancel(home,*this,PC_INT_BND); + home.ignore(*this,AP_WEAKLY); + (void) Propagator::dispose(home); + return sizeof(*this); + } + + template + 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_BND); + 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_BND); + 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_BND); + 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..43e6954fdb --- /dev/null +++ b/gecode/int/arithmetic/product.hpp @@ -0,0 +1,628 @@ +/* -*- 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; + } + + /// An interval whose endpoints have not been restricted to integer limits. + struct ProductWideInterval { + long long int min; + long long int max; + }; + + /// Quotient bounds for a zero-free divisor, before clipping the endpoints. + forceinline ProductWideInterval + product_quotient_bounds(long long int ymin, long long int ymax, + int qmin, int qmax) { + long long int c[4] = { + ceil_div_xx(ymin,static_cast(qmin)), + ceil_div_xx(ymin,static_cast(qmax)), + ceil_div_xx(ymax,static_cast(qmin)), + ceil_div_xx(ymax,static_cast(qmax)) + }; + long long int f[4] = { + floor_div_xx(ymin,static_cast(qmin)), + floor_div_xx(ymin,static_cast(qmax)), + floor_div_xx(ymax,static_cast(qmin)), + floor_div_xx(ymax,static_cast(qmax)) + }; + ProductWideInterval r = { + *std::min_element(c,c+4), *std::max_element(f,f+4) + }; + return r; + } + + /// Compute a saturated hull of quotients by a zero-free interval. + forceinline ProductInterval + product_quotient_interval(int ymin, int ymax, int qmin, int qmax) { + const ProductWideInterval q=product_quotient_bounds(ymin,ymax,qmin,qmax); + const long long int l = std::max + (static_cast(Limits::min),std::min + (q.min,static_cast(Limits::max))); + const long long int u = std::max + (static_cast(Limits::min),std::min + (q.max,static_cast(Limits::max))); + ProductInterval r = {static_cast(l),static_cast(u)}; + return r; + } + + /// Hull of quotients by the nonzero parts of a cofactor interval. + forceinline ProductInterval + product_quotient_hull(ProductInterval y, ProductInterval q) { + ProductInterval d = {1,0}; + if (q.min < 0) + d=product_quotient_interval(y.min,y.max,q.min,std::min(q.max,-1)); + if (q.max > 0) { + const ProductInterval p=product_quotient_interval + (y.min,y.max,std::max(q.min,1),q.max); + if (q.min < 0) { + d.min=std::min(d.min,p.min); d.max=std::max(d.max,p.max); + } else { + d=p; + } + } + return d; + } + + /// Invert integer power bounds; min greater than max denotes no support. + inline ProductInterval + product_root_interval(ProductInterval x, int n, ProductInterval d) { + if (n == 1) + return d; + if ((n & 1) != 0) { + ProductInterval r = { + d.min < 0 ? -product_floor_root(-d.min,n) : + product_ceil_root(d.min,n), + d.max < 0 ? -product_ceil_root(-d.max,n) : + product_floor_root(d.max,n) + }; + return r; + } + if (d.max < 0) + return ProductInterval{1,0}; + 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 ProductInterval{1,0}; + const ProductInterval negative = { + std::max(x.min,-hi), std::min(x.max,-lo) + }; + const ProductInterval positive = { + std::max(x.min,lo), std::min(x.max,hi) + }; + if (negative.min > negative.max) + return positive; + if (positive.min > positive.max) + return negative; + return ProductInterval{negative.min,positive.max}; + } + + /// Compute product bounds, optionally omitting one factor. + inline ProductInterval + product_interval(const ViewArray& x, int omit=-1) { + ProductInterval r = {1,1}; + // Equal views are adjacent after posting and remain adjacent in clones. + for (int i=0; i= i) && (omit < j) ? 1 : 0); + if (n > 0) + r = product_interval_mul(r,product_power_interval(x[i],n)); + i=j; + } + return r; + } + + /// Evaluate a tuple without overflowing internal arithmetic. + inline bool + product_value(const int* v, int n, int& p) { + 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); + Support::quicksort(&x[0],x.size()); + (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::linear(PropCost::HI,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&) { + // Preserve adjacent equal views while removing assigned units in place. + int n=0; + bool zero=false; + 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; + + const ProductInterval result = { + neg ? -y.max() : y.min(), neg ? -y.min() : y.max() + }; + const ProductInterval d=product_quotient_hull(result,q); + if (d.min > d.max) + return ES_FAILED; + const ProductInterval bounds={xi.min(),xi.max()}; + const ProductInterval roots=product_root_interval(bounds,exponent[g],d); + if (roots.min > roots.max) + return ES_FAILED; + ModEvent me=xi.gq(home,roots.min); + if (me_failed(me)) return ES_FAILED; + modified |= me_modified(me); + me=xi.lq(home,roots.max); + 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); + Support::quicksort(&x[0],x.size()); + (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::linear(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..8ac25e965f 100644 --- a/test/int/arithmetic.cpp +++ b/test/int/arithmetic.cpp @@ -35,8 +35,10 @@ #include #include +#include #include +#include namespace Test { namespace Int { @@ -48,6 +50,1857 @@ 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 { + /// Whether to exchange the two operands + bool swap; + public: + /// Create and register test + GcdXYX(const std::string& s, const Gecode::IntSet& d, + Gecode::IntPropLevel ipl, bool r=true, bool sw=false) + : Test("Arithmetic::Gcd::XYX::"+str(ipl)+"::"+s+ + (r ? "" : "::Plain")+(sw ? "::Swap" : ""),2,d,r,ipl), + swap(sw) { + 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) { + if (swap) + Gecode::gcd(home,x[1],x[0],x[0],ipl); + else + 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) { + if (swap) + Gecode::gcd(home,x[1],x[0],x[0],r,ipl); + else + 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; + } + + /// Algebraic deductions must not wait for unrelated variables to be fixed. + class NumberTheoryAlgebraic : public ::Test::Base { + class TestSpace : public Gecode::Space { + public: + virtual Gecode::Space* copy(void) { return nullptr; } + }; + int kind; + public: + NumberTheoryAlgebraic(int k) + : ::Test::Base("Int::Arithmetic::NumberTheoryAlgebraic::"+str(k)), + kind(k) {} + virtual bool run(void) { + using namespace Gecode; + TestSpace home; + if (kind <= 2) { + IntVar a(home,6,6), g(home,3,3); + IntVar x(home,kind == 1 ? -12 : 6,kind == 1 ? -6 : 12); + if (kind == 2) dom(home,x,IntSet({6,9,12})); + gcd(home,a,x,g); + if ((home.status() == SS_FAILED) || !x.assigned() || + (x.val() != (kind == 1 ? -9 : 9))) + return arithmetic_failed("coprime quotient: a, x, g", + IntVarArgs() << a << x << g); + } else if (kind == 3) { + IntVar a(home,12,12), x(home,1,24), g(home,7,12); + gcd(home,a,x,g); + if ((home.status() == SS_FAILED) || !g.assigned() || (g.val()!=12)) + return arithmetic_failed("divisor quotient band: a, x, g", + IntVarArgs() << a << x << g); + } else if (kind <= 5) { + IntVar d(home,kind == 4 ? 6 : -10,kind == 4 ? 10 : -6); + IntVar x(home,1,5); BoolVar b(home,1,1); + divides(home,d,x,Reify(b)); + if (home.status() != SS_FAILED) + return arithmetic_failed("divisor magnitude: d, x", + IntVarArgs() << d << x); + } else if (kind == 6) { + IntVar d(home,0,0), x(home,-5,5); BoolVar b(home,0,0); + divides(home,d,x,Reify(b)); + if ((home.status() == SS_FAILED) || x.in(0)) + return arithmetic_failed("negated zero divisor: d, x", + IntVarArgs() << d << x); + } else if (kind == 7) { + IntVar c(home,6,6), x(home,1,100), z(home,1,100), y(home,1,1); + product_mod(home,IntVarArgs({c,x,z}),12,y); + if (home.status() != SS_FAILED) + return arithmetic_failed("coefficient gcd failure: c, x, z, y", + IntVarArgs() << c << x << z << y); + } else if (kind == 8) { + IntVar c(home,-6,-6), x(home,-100,100), z(home,-100,100); + IntVar y(home,4,9); + product_mod(home,IntVarArgs({c,x,z}),12,y); + if ((home.status() == SS_FAILED) || !y.assigned() || (y.val()!=6)) + return arithmetic_failed("coefficient result bounds: y", + IntVarArgs() << y); + } else if (kind == 9) { + IntVar a(home,2,2), b(home,3,3), x(home,-100,100), y(home,0,5); + product_mod(home,IntVarArgs({a,x,b}),6,y); + if ((home.status() == SS_FAILED) || !y.assigned() || (y.val()!=0)) + return arithmetic_failed("collective zero coefficient: y", + IntVarArgs() << y); + } else if (kind == 10) { + IntVar c(home,6,6), x(home,1,100), z(home,1,100); + IntVar y(home,IntSet({1,5,7,11})); + product_mod(home,IntVarArgs({c,x,z}),12,y); + if (home.status() != SS_FAILED) + return arithmetic_failed("sparse result multiples: y", + IntVarArgs() << y); + } else if (kind == 11) { + const int hi=Gecode::Int::Limits::max; + IntVar c(home,6,6), x(home,1,hi), z(home,1,hi), y(home,1,1); + product_mod(home,IntVarArgs({c,x,z}),hi,y); + if (home.status() != SS_FAILED) + return arithmetic_failed("large coefficient gcd failure: y", + IntVarArgs() << y); + } else { + const int hi=Gecode::Int::Limits::max; + IntVar a(home,-hi,-hi), x(home,hi-1,hi), g(home,1,1); + gcd(home,a,x,g); + if ((home.status() == SS_FAILED) || !x.assigned() || (x.val()!=hi-1)) + return arithmetic_failed("limit coprimality: a, x, g", + IntVarArgs() << a << x << g); + } + return true; + } + }; + + /// 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); + } + }; + + /// Interleaved powers exercise grouping, cofactors, and cloned rewrites. + class ProductInterleaved : public Test { + int kind; + Gecode::IntVarArgs factors(const Gecode::IntVarArray& x) const { + return Gecode::IntVarArgs({x[0],x[1],x[0],x[2],x[1],x[0]}); + } + public: + ProductInterleaved(int k) + : Test("Arithmetic::Product::Interleaved::"+str(k),k == 2 ? 5 : 4, + Gecode::IntSet(-2,3),true,Gecode::IPL_BND), kind(k) { + contest=CTL_NONE; testfix=false; + } + virtual bool solution(const Assignment& x) const { + const int p=x[0]*x[0]*x[0]*x[1]*x[1]*x[2]; + if (kind == 0) return p == x[3]; + const int m=kind == 1 ? 7 : x[4]; + return (m > 0) && (((p % m)+m) % m == x[3]); + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { + if (kind == 0) Gecode::product(home,factors(x),x[3],ipl); + else if (kind == 1) + Gecode::product_mod(home,factors(x),7,x[3],ipl); + else Gecode::product_mod(home,factors(x),x[4],x[3],ipl); + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x, + Gecode::Reify r) { + if (kind == 0) Gecode::product(home,factors(x),x[3],r,ipl); + else if (kind == 1) + Gecode::product_mod(home,factors(x),7,x[3],r,ipl); + else Gecode::product_mod(home,factors(x),x[4],x[3],r,ipl); + } + }; + + /// %Test forward and inverse product bounds on large factor domains + 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); + } + }; + + /// Assigned coefficients with several free factors and all reification modes. + class ProductModCoefficient : public Test { + bool variable; + public: + ProductModCoefficient(bool v) + : Test("Arithmetic::ProductMod::Coefficient::"+ + std::string(v ? "Variable" : "Fixed"),3, + Gecode::IntSet(-3,6),true,Gecode::IPL_BND), variable(v) { + contest=CTL_NONE; testfix=false; + } + virtual bool solution(const Assignment& x) const { + const int p=6*x[0]*x[1]; + return ((p%12)+12)%12 == x[2]; + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { + Gecode::IntVar c(home,6,6); + const Gecode::IntVarArgs factors({x[0],c,x[1]}); + if (variable) { + Gecode::IntVar m(home,12,12); + Gecode::product_mod(home,factors,m,x[2],ipl); + } else { + Gecode::product_mod(home,factors,12,x[2],ipl); + } + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x, + Gecode::Reify r) { + Gecode::IntVar c(home,6,6); + const Gecode::IntVarArgs factors({x[0],c,x[1]}); + if (variable) { + Gecode::IntVar m(home,12,12); + Gecode::product_mod(home,factors,m,x[2],r,ipl); + } else { + Gecode::product_mod(home,factors,12,x[2],r,ipl); + } + } + }; + + /// Folded coefficients, repeated factors, and products beyond machine range. + class ProductModFolded : public Test { + int c, m; + public: + ProductModFolded(int c0, int m0, bool r) + : Test("Arithmetic::ProductMod::Folded::"+str(c0)+"::"+str(m0)+ + (r ? "::Reified" : "::Plain"),3,Gecode::IntSet(-2,6), + r,Gecode::IPL_BND), c(c0), m(m0) { + contest=CTL_NONE; testfix=!r; + } + virtual bool solution(const Assignment& x) const { + long long int p=1; + const int v[]={c,x[0],c,x[0],x[1],c}; + for (int i=0; i<6; i++) + p=(p*((v[i]%m+m)%m))%m; + return p == x[2]; + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x) { + Gecode::IntVar a(home,c,c); + Gecode::product_mod(home, + Gecode::IntVarArgs({a,x[0],a,x[0],x[1],a}), + m,x[2],ipl); + } + virtual void post(Gecode::Space& home, Gecode::IntVarArray& x, + Gecode::Reify r) { + Gecode::IntVar a(home,c,c); + Gecode::product_mod(home, + Gecode::IntVarArgs({a,x[0],a,x[0],x[1],a}), + m,x[2],r,ipl); + } + }; + + /// A nonzero modular result excludes every zero factor. + class ProductModNonzero : public ::Test::Base { + class TestSpace : public Gecode::Space { + public: + virtual Gecode::Space* copy(void) { return nullptr; } + }; + public: + ProductModNonzero(void) + : ::Test::Base("Int::Arithmetic::ProductMod::Nonzero") {} + virtual bool run(void) { + using namespace Gecode; + for (int negative=0; negative<2; negative++) + for (int variable=0; variable<2; variable++) { + TestSpace home; + IntVarArgs x(home,4,negative ? -1 : 0,negative ? 0 : 1); + IntVar y(home,0,6), m(home,7,11); + if (variable) + product_mod(home,x,m,y); + else + product_mod(home,x,7,y); + if (home.status() == SS_FAILED) + return false; + rel(home,y,IRT_EQ,1); + if (home.status() == SS_FAILED) + return arithmetic_failed("nonzero result failed",x); + for (int i=0; i& x, + Gecode::Int::IntView y, size_t& d) + : ProductMod(home,x,7,y,1,1), disposed(d) {} + virtual size_t dispose(Gecode::Space& home) { + disposed=ProductMod::dispose(home); + // The probe has its own field, beyond the actor being checked. + return sizeof(*this); + } + }; + public: + ProductModDispose(void) + : ::Test::Base("Int::Arithmetic::ProductMod::Dispose") {} + virtual bool run(void) { + using namespace Gecode; + size_t disposed=0; + TestSpace home; + IntVarArgs factors(home,2,1,3); + IntVar y(home,0,6); + ViewArray x(home,factors); + (void) new (home) DisposeProbe(home,x,y,disposed); + PropagatorGroup::all.kill(home); + if (disposed != sizeof(Gecode::Int::Arithmetic::ProductMod)) { + if (opt.log) + olog << "disposed " << disposed << " bytes, expected " + << sizeof(Gecode::Int::Arithmetic::ProductMod) << std::endl; + return false; + } + return true; + } + }; + + /// %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; + } + }; + + /// Modulus bounds must reach a fixpoint after jumping over a domain hole. + class ProductModVarSparseModulus : public ::Test::Base { + class TestSpace : public Gecode::Space { + public: + virtual Gecode::Space* copy(void) { return nullptr; } + }; + /// Whether the remaining modulus domain contains a solution + bool supported; + public: + ProductModVarSparseModulus(bool s) + : ::Test::Base("Int::Arithmetic::ProductModVar::SparseModulus::"+ + std::string(s ? "Assignment" : "Failure")), + supported(s) {} + virtual bool run(void) { + using namespace Gecode; + TestSpace home; + IntVar p(home,60,60), y(home,4,4); + IntVar m(home,IntSet({20,supported ? 28 : 25,100})); + product_mod(home,IntVarArgs({p}),m,y); + branch(home,m,INT_VAL_MIN()); + const SpaceStatus status=home.status(); + // m must divide 60-4=56. Dropping 100 exposes a single quotient + // band, which must be reconsidered without an external wakeup. + if (supported) { + if ((status == SS_FAILED) || !m.assigned() || (m.val() != 28)) + return arithmetic_failed + ("sparse modulus must become 28: p, m, y", + IntVarArgs() << p << m << y); + } else if (status != SS_FAILED) { + return arithmetic_failed + ("sparse modulus has no divisor of 56: p, m, y", + IntVarArgs() << p << 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 +2261,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 +2309,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 +2938,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 +2960,55 @@ 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 GcdXYZ("Limits",p,ipls.ipl()); + (void) new GcdXYZ("Limits",p,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); + for (int sw=0; sw<2; sw++) { + (void) new GcdXYX("Nonpositive",Gecode::IntSet(-2,0), + ipls.ipl(),true,sw != 0); + (void) new GcdXYX("Nonpositive",Gecode::IntSet(-2,0), + ipls.ipl(),false,sw != 0); + } + (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 +3158,36 @@ namespace Test { namespace Int { (void) new ArgMinBool(i,1,false); (void) new ArgMinBoolShared(i,false); } + (void) new NumberTheorySparseBounds; + for (int k=0; k<13; k++) + (void) new NumberTheoryAlgebraic(k); + (void) new NumberTheoryLifecycle; + (void) new NumberTheoryLimits; + (void) new NumberTheoryIdentities; + (void) new ProductModInvalidModulus; + (void) new ProductModAlgebraic; + (void) new ProductModNonzero; + (void) new ProductCompactionZero; + for (int r=0; r<2; r++) { + (void) new ProductModFolded(-2,7,r != 0); + (void) new ProductModFolded(6,12,r != 0); + (void) new ProductModFolded(6,15,r != 0); + (void) new ProductModFolded(Gecode::Int::Limits::max,5,r != 0); + } + for (int k=0; k<3; k++) + (void) new ProductInterleaved(k); + (void) new ProductModDispose; + (void) new ProductModCoefficient(false); + (void) new ProductModCoefficient(true); + (void) new ProductModVarBounds; + (void) new ProductModVarSparseModulus(false); + (void) new ProductModVarSparseModulus(true); + (void) new ProductModVarAlgebraic; + (void) new ProductModVarInactive; + (void) new ProductBoundsLarge; + (void) new ProductSimplifySign; + (void) new ProductPowerAlias; + (void) new ProductInverseBounds; } };