diff --git a/changelog.in b/changelog.in
index 451776aa77..c20924b2be 100755
--- a/changelog.in
+++ b/changelog.in
@@ -103,6 +103,17 @@ Fix float power and nth-root tests with shared variables by
 initializing their exponent fields and rejecting nth-root test
 assignments whose source interval is not fully non-negative.
 
+[ENTRY]
+Module: int
+What:   performance
+Rank:   minor
+Issue:  202
+Thanks: Fabio Tardivo
+[DESCRIPTION]
+Strengthen bin-packing propagation with lower bounds based on
+dual-feasible functions, improving early detection of infeasible
+packing states.
+
 [ENTRY]
 Module: other
 What:   change
diff --git a/gecode/int/bin-packing.hh b/gecode/int/bin-packing.hh
index 909e00be62..54be725a98 100755
--- a/gecode/int/bin-packing.hh
+++ b/gecode/int/bin-packing.hh
@@ -4,9 +4,11 @@
  *     Christian Schulte <schulte@gecode.dev>
  *
  *  Contributing authors:
+ *     Fabio Tardivo <ftardivo@nmsu.edu>
  *     Stefano Gualandi <stefano.gualandi@gmail.com>
  *
  *  Copyright:
+ *     Fabio Tardivo, 2024
  *     Stefano Gualandi, 2013
  *     Christian Schulte, 2010
  *
@@ -127,12 +129,18 @@ namespace Gecode { namespace Int { namespace BinPacking {
     int operator [](int i) const;
   };
 
+  /// Range of lambda values
+  struct LambdaRange {
+    int min;
+    int max;
+  };
 
   /**
    * \brief Bin-packing propagator
    *
    * The algorithm is taken from:
    *   Paul Shaw. A Constraint for Bin Packing. CP 2004.
+   *   Tardivo et al. CP for Bin Packing with Multi-Core and GPUs. CP 2024.
    *
    * Requires \code #include <gecode/int/bin-packing.hh> \endcode
    *
@@ -175,6 +183,40 @@ namespace Gecode { namespace Int { namespace BinPacking {
     virtual Actor* copy(Space& home);
     /// Destructor
     virtual size_t dispose(Space& home);
+    /// Reductions
+    static int const n_reductions = 3;
+    static void calc_reductions(const ViewArray<Item>& bs,
+                                const ViewArray<OffsetView>& l,
+                                int* weights_base_reduction,
+                                int& capacity_base_reduction,
+                                int* delta_reductions);
+    /// Dual Feasible Functions
+    static int f_ccm1(int w, int l, int c);
+    static int f_mt(int w, int l, int c);
+    static int f_bj1(int w, int l, int c);
+    static int f_vb2_base(int w, int l, int c);
+    static int f_vb2(int w, int l, int c);
+    static int f_fs1(int w, int l, int c);
+    static int f_rad2_base(int w, int l, int c);
+    static int f_rad2(int w, int l, int c);
+    static int const n_lambda_samples = 256;
+    static LambdaRange l_ccm1(int c);
+    static LambdaRange l_mt(int c);
+    static LambdaRange l_bj1(int c);
+    static LambdaRange l_vb2(int c);
+    static LambdaRange l_fs1(int c);
+    static LambdaRange l_rad2(int c);
+    static LambdaRange sanitize_lambda_range(LambdaRange lambda,
+                                             int n_weights, int max_weight);
+    /// Lower bound
+    template<int f(int,int,int)>
+    static int calc_dff_lower_bound_single_lambda(const int* weights,
+                                                  int n_weights,
+                                                  int capacity, int lambda);
+    template<int f(int,int,int), LambdaRange l(int)>
+    static int calc_dff_lower_bound(const int* weights, int n_weights,
+                                    int capacity, int n_not_zero_weights,
+                                    int max_weight, bool sanitize = false);
   };
 
 
diff --git a/gecode/int/bin-packing/propagate.cpp b/gecode/int/bin-packing/propagate.cpp
index 9425fdd460..43c1beb6a2 100755
--- a/gecode/int/bin-packing/propagate.cpp
+++ b/gecode/int/bin-packing/propagate.cpp
@@ -3,7 +3,11 @@
  *  Main authors:
  *     Christian Schulte <schulte@gecode.dev>
  *
+ *  Contributing authors:
+ *     Fabio Tardivo <ftardivo@nmsu.edu>
+ *
  *  Copyright:
+ *     Fabio Tardivo, 2024
  *     Christian Schulte, 2010
  *
  *  This file is part of Gecode, the generic constraint
@@ -278,77 +282,83 @@ namespace Gecode { namespace Int { namespace BinPacking {
     // Perform lower bound checking
     if (n > 0) {
 
-      // Find capacity estimate (we start from bs[0] as it might be
-      // not packable, actually (will be detected later anyway)!
-      int c = bs[0].size();
-      for (int j=0; j<m; j++)
-        c = std::max(c,l[j].max());
-
-      // Count how many items have a certain size (bucket sort)
-      int* n_s = region.alloc<int>(c+1);
-
-      for (int i=0; i<c+1; i++)
-        n_s[i] = 0;
-
-      // Count unpacked items
-      for (int i=0; i<n; i++)
-        n_s[bs[i].size()]++;
-
-      // Number of items and remaining bin load
-      int nm = n;
-
-      // Only count positive remaining bin loads
-      for (int j=0; j<m; j++)
-        if (l[j].max() < 0) {
-          return ES_FAILED;
-        } else if (c > l[j].max()) {
-          n_s[c - l[j].max()]++; nm++;
-        }
-
-      // Sizes of items and remaining bin loads
-      int* s = region.alloc<int>(nm);
-
-      // Setup sorted sizes
       {
-        int k=0;
-        for (int i=c+1; i--; )
-          for (int j=n_s[i]; j--; )
-            s[k++]=i;
-        assert(k == nm);
-      }
-
-      // Items in N1 are from 0 ... n1 - 1
-      int n1 = 0;
-      // Items in N2 are from n1 ... n12 - 1, we count elements in N1 and N2
-      int n12 = 0;
-      // Items in N3 are from n12 ... n3 - 1
-      int n3 = 0;
-      // Free space in N2
-      int f2 = 0;
-      // Total size of items in N3
-      int s3 = 0;
-
-      // Initialize n12 and f2
-      for (; (n12 < nm) && (s[n12] > c/2); n12++)
-        f2 += c - s[n12];
-
-      // Initialize n3 and s3
-      for (n3 = n12; n3 < nm; n3++)
-        s3 += s[n3];
-
-      // Compute lower bounds
-      for (int k=0; k<=c/2; k++) {
-        // Make N1 larger by adding elements and N2 smaller
-        for (; (n1 < nm) && (s[n1] > c-k); n1++)
-          f2 -= c - s[n1];
-        assert(n1 <= n12);
-        // Make N3 smaller by removing elements
-        for (; (s[n3-1] < k) && (n3 > n12); n3--)
-          s3 -= s[n3-1];
-        // Overspill
-        int o = (s3 > f2) ? ((s3 - f2 + c - 1) / c) : 0;
-        if (n12 + o > m)
-          return ES_FAILED;
+        // Allocate auxiliary data
+        int n_bins = l.size();
+        int n_weights_reductions = n_bins + bs.size();
+        int* weights_base_reduction =
+          region.alloc<int>(n_weights_reductions);
+        int* weights_current_reduction =
+          region.alloc<int>(n_weights_reductions);
+        int* delta_reductions = region.alloc<int>(n_reductions);
+
+        // Initialize reductions
+        int capacity_base_reduction = 0;
+        calc_reductions(bs, l, weights_base_reduction,
+                        capacity_base_reduction, delta_reductions);
+
+        for (int r_idx = 0; r_idx < n_reductions; r_idx += 1) {
+          // Calculate reduction parameters
+          int& delta = delta_reductions[r_idx];
+          int capacity = capacity_base_reduction + delta;
+
+          if (capacity > 0) {
+            // Calculate reduction
+            int n_not_zero_weights = 0;
+            int max_weight = 0;
+            int* weights = weights_current_reduction;
+            for (int w_idx = 0; w_idx < n_weights_reductions; w_idx += 1) {
+              int weight =
+                weights_base_reduction[w_idx] + (w_idx < n_bins ? delta : 0);
+              n_not_zero_weights += weight != 0;
+              max_weight = std::max(max_weight, weight);
+              weights[w_idx] = weight;
+            }
+
+            // Check lower bounds
+            int lower_bound =
+              calc_dff_lower_bound<f_ccm1, l_ccm1>
+                (weights, n_weights_reductions, capacity,
+                 n_not_zero_weights, max_weight);
+            if (lower_bound > n_bins)
+              return ES_FAILED;
+
+            lower_bound =
+              calc_dff_lower_bound<f_mt, l_mt>
+                (weights, n_weights_reductions, capacity,
+                 n_not_zero_weights, max_weight);
+            if (lower_bound > n_bins)
+              return ES_FAILED;
+
+            lower_bound =
+              calc_dff_lower_bound<f_bj1, l_bj1>
+                (weights, n_weights_reductions, capacity,
+                 n_not_zero_weights, max_weight);
+            if (lower_bound > n_bins)
+              return ES_FAILED;
+
+            lower_bound =
+              calc_dff_lower_bound<f_vb2, l_vb2>
+                (weights, n_weights_reductions, capacity,
+                 n_not_zero_weights, max_weight, true);
+            if (lower_bound > n_bins)
+              return ES_FAILED;
+
+            lower_bound =
+              calc_dff_lower_bound<f_fs1, l_fs1>
+                (weights, n_weights_reductions, capacity,
+                 n_not_zero_weights, max_weight, true);
+            if (lower_bound > n_bins)
+              return ES_FAILED;
+
+            lower_bound =
+              calc_dff_lower_bound<f_rad2, l_rad2>
+                (weights, n_weights_reductions, capacity,
+                 n_not_zero_weights, max_weight);
+            if (lower_bound > n_bins)
+              return ES_FAILED;
+          }
+        }
       }
       region.free();
     }
@@ -394,7 +404,49 @@ namespace Gecode { namespace Int { namespace BinPacking {
     }
   }
 
+  void
+  Pack::calc_reductions(const ViewArray<Item>& bs,
+                        const ViewArray<OffsetView>& l,
+                        int* weights_base_reduction,
+                        int& capacity_base_reduction,
+                        int* delta_reductions) {
+    // Reset values
+    int n_weights_reductions = l.size() + bs.size();
+    for (int w_idx = 0; w_idx < n_weights_reductions; w_idx += 1)
+      weights_base_reduction[w_idx] = 0;
+
+    // R0
+    int n_bins = l.size();
+    int n_items = bs.size();
+    capacity_base_reduction = 0;
+    for (int b_idx = 0; b_idx < n_bins; b_idx += 1)
+      capacity_base_reduction =
+        std::max(capacity_base_reduction, l[b_idx].max());
+    for (int b_idx = 0; b_idx < n_bins; b_idx += 1)
+      weights_base_reduction[b_idx] =
+        capacity_base_reduction - l[b_idx].max();
+    for (int i_idx = 0; i_idx < n_items; i_idx += 1) {
+      bool i_assigned = bs[i_idx].bin().assigned();
+      int i_weight = bs[i_idx].size();
+      if (i_assigned) {
+        int b_idx = bs[i_idx].bin().val();
+        weights_base_reduction[b_idx] += i_weight;
+      } else {
+        weights_base_reduction[n_bins + i_idx] = i_weight;
+      }
+    }
+
+    // RMin, RMax
+    int smallest_virtual_weight = std::numeric_limits<int>::max();
+    for (int b_idx = 0; b_idx < n_bins; b_idx += 1)
+      smallest_virtual_weight =
+        std::min(smallest_virtual_weight, weights_base_reduction[b_idx]);
+    delta_reductions[0] = -smallest_virtual_weight;
+    delta_reductions[1] = 0;
+    delta_reductions[2] =
+      capacity_base_reduction - 2 * smallest_virtual_weight + 1;
+  }
+
 }}}
 
 // STATISTICS: int-prop
-
diff --git a/gecode/int/bin-packing/propagate.hpp b/gecode/int/bin-packing/propagate.hpp
index 36241a3a0d..1c6e41eb98 100755
--- a/gecode/int/bin-packing/propagate.hpp
+++ b/gecode/int/bin-packing/propagate.hpp
@@ -3,7 +3,11 @@
  *  Main authors:
  *     Christian Schulte <schulte@gecode.dev>
  *
+ *  Contributing authors:
+ *     Fabio Tardivo <ftardivo@nmsu.edu>
+ *
  *  Copyright:
+ *     Fabio Tardivo, 2024
  *     Christian Schulte, 2010
  *
  *  This file is part of Gecode, the generic constraint
@@ -208,7 +212,204 @@ namespace Gecode { namespace Int { namespace BinPacking {
     return nosum(s, a, b, ap, bp);
   }
 
+  forceinline int
+  Pack::f_ccm1(int w, int l, int c) {
+    // Conditions
+    int const c0 = 2 * w > c;  // x > c / 2
+    int const c1 = 2 * w == c; // x == c / 2
+    int const c2 = 2 * w < c;  // x < c / 2
+
+    // Values
+    int const v0 = 2 * ((c / l) - ((c - w) / l));
+    int const v1 = c / l;
+    int const v2 = 2 * (w / l);
+
+    return (c0 * v0) + (c1 * v1) + (c2 * v2);
+  }
+
+  forceinline int
+  Pack::f_mt(int w, int l, int c) {
+    // Conditions
+    int const c0 = w < l;
+    int const c1 = (l <= w) && (w <= c - l);
+    int const c2 = c - l < w;
+
+    // Values
+    int const v0 = 0;
+    int const v1 = w;
+    int const v2 = c;
+
+    return (c0 * v0) + (c1 * v1) + (c2 * v2);
+  }
+
+  forceinline int
+  Pack::f_bj1(int w, int l, int c) {
+    // Auxiliary values
+    int const p = l - (c % l);
+
+    // Conditions
+    int const c0 = w % l <= c % l;
+    int const c1 = w % l > c % l;
+
+    // Values
+    int const v0 = (w / l) * p;
+    int const v1 = (w / l) * p + (w % l) - (c % l);
+
+    return (c0 * v0) + (c1 * v1);
+  }
+
+  forceinline int
+  Pack::f_vb2_base(int w, int l, int c) {
+    int v = ceil_div_pp(l * w, c);
+    return v > 0 ? v - 1 : 0;
+  }
+
+  forceinline int
+  Pack::f_vb2(int w, int l, int c) {
+    int const c0 = 2 * w > c;
+    int const c1 = 2 * w == c;
+    int const c2 = 2 * w < c;
+
+    int const t0 = f_vb2_base(c, l, c);
+    int const t1 = f_vb2_base(w, l, c);
+    int const t2 = f_vb2_base(c - w, l, c);
+
+    int const v0 = 2 * t0 - 2 * t2;
+    int const v1 = t0;
+    int const v2 = 2 * t1;
+
+    return (c0 * v0) + (c1 * v1) + (c2 * v2);
+  }
+
+  forceinline int
+  Pack::f_fs1(int w, int l, int c) {
+    // Conditions
+    int const c0 = w * (l + 1) % c == 0;
+    int const c1 = w * (l + 1) % c != 0;
+
+    // Values
+    int const v0 = w * l;
+    int const v1 = ((w * (l + 1)) / c) * c;
+
+    return (c0 * v0) + (c1 * v1);
+  }
+
+  forceinline int
+  Pack::f_rad2_base(int w, int l, int c) {
+    // Conditions
+    int const c0 = w < l;
+    int const c1 = (l <= w) && (w <= c - 2 * l);
+    int const c2 = (c - 2 * l < w) && (w < 2 * l);
+
+    // Values
+    int const v0 = 0;
+    int const v1 = c / 3;
+    int const v2 = c / 2;
+
+    return (c0 * v0) + (c1 * v1) + (c2 * v2);
+  }
+
+  forceinline int
+  Pack::f_rad2(int w, int l, int c) {
+    // Conditions
+    int const c0 = w < 2 * l;
+    int const c1 = 2 * l <= w;
+
+    // Values
+    int const v0 = f_rad2_base(w, l, c);
+    int const v1 = c - f_rad2_base(c - w, l, c);
+
+    return (c0 * v0) + (c1 * v1);
+  }
+
+  forceinline LambdaRange
+  Pack::l_ccm1(int c) {
+    return {1, c / 2};
+  }
+
+  forceinline LambdaRange
+  Pack::l_mt(int c) {
+    return {0, c / 2};
+  }
+
+  forceinline LambdaRange
+  Pack::l_bj1(int c) {
+    return {1, c};
+  }
+
+  forceinline LambdaRange
+  Pack::l_vb2(int c) {
+    return {2, c};
+  }
+
+  forceinline LambdaRange
+  Pack::l_fs1(int) {
+    return {1, 100};
+  }
+
+  forceinline LambdaRange
+  Pack::l_rad2(int c) {
+    return {c / 4 + 1, c / 3};
+  }
+
+  forceinline LambdaRange
+  Pack::sanitize_lambda_range(LambdaRange lambda_range,
+                              int n_not_zero_weights, int max_weight) {
+    long long int n = static_cast<long long int>(n_not_zero_weights) *
+      static_cast<long long int>(max_weight);
+    if (n != 0) {
+      int l_max =
+        std::min(static_cast<int>(std::numeric_limits<int>::max() / n),
+                 lambda_range.max);
+      return {lambda_range.min, l_max};
+    } else {
+      return {0, -1};
+    }
+  }
+
+  template<int f(int,int,int)>
+  forceinline int
+  Pack::calc_dff_lower_bound_single_lambda(const int* weights, int n_weights,
+                                           int capacity, int lambda) {
+    // Transform and sum the weights
+    int sum_transformed_weights = 0;
+    for (int w_idx = 0; w_idx < n_weights; w_idx += 1)
+      sum_transformed_weights += f(weights[w_idx], lambda, capacity);
+
+    // Transform the capacity
+    int transformed_capacity = f(capacity, lambda, capacity);
+
+    // Lower bound
+    return ceil_div_pp(sum_transformed_weights, transformed_capacity);
+  }
+
+  template<int f(int,int,int), LambdaRange l(int)>
+  forceinline int
+  Pack::calc_dff_lower_bound(const int* weights, int n_weights,
+                             int capacity, int n_not_zero_weights,
+                             int max_weight, bool sanitize) {
+    LambdaRange lambda_range = l(capacity);
+    lambda_range =
+      sanitize ? sanitize_lambda_range(lambda_range, n_not_zero_weights,
+                                       max_weight) :
+                 lambda_range;
+    if (lambda_range.min >= lambda_range.max)
+      return 0;
+
+    int lower_bound = 0;
+    int l_step =
+      ceil_div_pp(lambda_range.max - lambda_range.min + 1,
+                  n_lambda_samples + 1);
+    for (int lambda = lambda_range.min + l_step;
+         lambda < lambda_range.max; lambda += l_step) {
+      int cur_lower_bound =
+        calc_dff_lower_bound_single_lambda<f>(weights, n_weights,
+                                              capacity, lambda);
+      lower_bound = std::max(lower_bound, cur_lower_bound);
+    }
+    return lower_bound;
+  }
+
 }}}
 
 // STATISTICS: int-prop
-
diff --git a/test/int/bin-packing.cpp b/test/int/bin-packing.cpp
index dcf78ada03..f35ca8cd28 100755
--- a/test/int/bin-packing.cpp
+++ b/test/int/bin-packing.cpp
@@ -309,6 +309,37 @@ namespace Test { namespace Int {
       }
     };
 
+    /// Test that DFF lower bounds detect infeasible packing
+    class DFFLowerBound : public Base {
+    protected:
+      /// Simple test space class
+      class TestSpace : public Gecode::Space {
+      public:
+        /// Constructor
+        TestSpace(void) {}
+        /// Copy function
+        virtual Gecode::Space* copy(void) {
+          return nullptr;
+        }
+      };
+    public:
+      /// Constructor
+      DFFLowerBound(void)
+        : Base("Int::BinPacking::DFFLowerBound") {}
+      /// Run the actual test
+      virtual bool run(void) {
+        using namespace Gecode;
+        TestSpace* home = new TestSpace;
+        IntVarArgs l(*home, 3, 0, 5);
+        IntVarArgs b(*home, 6, 0, 2);
+        IntArgs s({4,2,2,2,2,2});
+        binpacking(*home, l, b, s);
+        bool failed = home->status() == SS_FAILED;
+        delete home;
+        return failed;
+      }
+    };
+
     /// Help class to create and register tests
     class Create {
     public:
@@ -392,6 +423,8 @@ namespace Test { namespace Int {
           (void) new CliqueMBPT(c48);
           (void) new CliqueMBPT(c49);
         }
+
+        (void) new DFFLowerBound;
       }
     };
 
@@ -405,4 +438,3 @@ namespace Test { namespace Int {
 
 
 // STATISTICS: test-int
-