diff --git a/gecode/int/bin-packing.hh b/gecode/int/bin-packing.hh
index ef1cc329b9..51a58109c5 100755
--- a/gecode/int/bin-packing.hh
+++ b/gecode/int/bin-packing.hh
@@ -127,12 +127,18 @@ namespace Gecode { namespace Int { namespace BinPacking {
     int operator [](int i) const;
   };
 
+  /// Range of lambda values  
+  struct LambdaRange {int min; int max;};
+
+  /// Integer Dynamic Array
+  using IntDynamicArray = Support::DynamicArray<int,Region>;
 
   /**
    * \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 +181,31 @@ namespace Gecode { namespace Int { namespace BinPacking {
     virtual Actor* copy(Space& home);
     /// Destructor
     virtual size_t dispose(Space& home);
+    /// Reductions
+    static int const nReductions = 3;
+    static void calcReductions(const ViewArray<Item>& bs, const ViewArray<OffsetView>& l, IntDynamicArray & weightsBaseReduction, int & capacityBaseReduction, IntDynamicArray & deltaReductions);
+    /// Dual Feasible Functions
+    static int fCCM1(int w, int l, int c);
+    static int fMT(int w, int l, int c);
+    static int fBJ1(int w, int l, int c);
+    static int fVB2Base(int w, int l, int c);
+    static int fVB2(int w, int l, int c);
+    static int fFS1(int w, int l, int c);
+    static int fRAD2Base(int w, int l, int c);
+    static int fRAD2(int w, int l, int c);
+    static int const nLambdaSamples = 256;
+    static LambdaRange lCCM1(int c);
+    static LambdaRange lMT(int c);
+    static LambdaRange lBJ1(int c);
+    static LambdaRange lVB2(int c);
+    static LambdaRange lFS1(int c);
+    static LambdaRange lRAD2(int c);
+    static LambdaRange sanitizeLambdaRange(LambdaRange lambda, int nWeights, int maxWeight);
+    /// Lower bound
+    template<int f(int,int,int)>
+    static int calcDffLowerboundSingleLambda(const IntDynamicArray & weights, int capacity, int lambda);
+    template<int f(int,int,int), LambdaRange l(int)>
+    static int calcDffLowerbound(const IntDynamicArray & weights, int capacity, int nNotZeroWeights, int maxWeight, bool sanitize = false);
   };
 
 
diff --git a/gecode/int/bin-packing/propagate.cpp b/gecode/int/bin-packing/propagate.cpp
index 080b0cadd2..4b762a1a7f 100755
--- a/gecode/int/bin-packing/propagate.cpp
+++ b/gecode/int/bin-packing/propagate.cpp
@@ -278,77 +278,136 @@ 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());
+      bool useDffLowerbounds = true;
+      if (not useDffLowerbounds)
+      {
+        
+        // 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);
+        // 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;
+        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()]++;
+        // 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;
+        // 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++;
-        }
+        // 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);
+        // 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);
-      }
+        // 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;
+        // 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;
+        }
+      }
+      else
+      {
+        // Allocate auxiliary data
+        int nBins = l.size();
+        int nWeightsReductions = nBins + bs.size();
+        IntDynamicArray weightsBaseReduction(region, nWeightsReductions);
+        IntDynamicArray weightsCurrentReduction(region, nWeightsReductions);
+        IntDynamicArray deltaReductions(region, nReductions);
+
+        // Initialize reductions
+        int capacityBaseReduction = 0;
+        calcReductions(bs, l, weightsBaseReduction, capacityBaseReduction, deltaReductions);
+
+        for (int rIdx = 0; rIdx < nReductions; rIdx += 1)
+        {
+          // Calculate reduction parameters
+          int & delta = deltaReductions[rIdx];
+          int capacity = capacityBaseReduction + delta;
+
+          if (capacity > 0)
+          {
+            // Calculate reduction
+            int nNotZeroWeights = 0;
+            int maxWeight = 0;
+            IntDynamicArray & weights = weightsCurrentReduction;
+            for (int wIdx = 0; wIdx < nWeightsReductions; wIdx += 1)
+            {
+              int weight = weightsBaseReduction[wIdx] + (wIdx < nBins ? delta : 0);
+              nNotZeroWeights += weight != 0;
+              maxWeight = std::max(maxWeight, weight);
+              weights[wIdx] = weight;
+            }
+
+            // Check lowerbounds
+            int lowerbound = calcDffLowerbound<fCCM1, lCCM1>(weights, capacity, nNotZeroWeights, maxWeight);
+            if (lowerbound > nBins) return ES_FAILED;
+
+            lowerbound = calcDffLowerbound<fMT, lMT>(weights, capacity, nNotZeroWeights, maxWeight);
+            if (lowerbound > nBins) return ES_FAILED;
+
+            lowerbound = calcDffLowerbound<fBJ1, lBJ1>(weights, capacity, nNotZeroWeights, maxWeight);
+            if (lowerbound > nBins) return ES_FAILED;
+
+            lowerbound = calcDffLowerbound<fVB2, lVB2>(weights, capacity, nNotZeroWeights, maxWeight, true);
+            if (lowerbound > nBins) return ES_FAILED;
+
+            lowerbound = calcDffLowerbound<fFS1, lFS1>(weights, capacity, nNotZeroWeights, maxWeight, true);
+            if (lowerbound > nBins) return ES_FAILED;
+
+            lowerbound = calcDffLowerbound<fRAD2, lRAD2>(weights, capacity, nNotZeroWeights, maxWeight);
+            if (lowerbound > nBins) return ES_FAILED;
+          }
+        }
       }
       region.free();
     }
@@ -394,6 +453,54 @@ namespace Gecode { namespace Int { namespace BinPacking {
     }
   }
 
+  void 
+  Pack::calcReductions(const ViewArray<Item>& bs, const ViewArray<OffsetView>& l, IntDynamicArray & weightsBaseReduction, int & capacityBaseReduction, IntDynamicArray & deltaReductions)
+  {     
+      // Reset values
+      int nWeightsReductions = weightsBaseReduction.capacity();
+      for (int wIdx = 0; wIdx < nWeightsReductions; wIdx += 1)
+      {
+        weightsBaseReduction[wIdx] = 0;
+      }
+
+      // R0
+      int nBins = l.size();
+      int nItems = bs.size();
+      capacityBaseReduction = 0;
+      for (int bIdx = 0; bIdx < nBins; bIdx += 1)
+      {
+          capacityBaseReduction = std::max(capacityBaseReduction, l[bIdx].max());
+      }
+      for (int bIdx = 0; bIdx < nBins; bIdx += 1)
+      {
+          weightsBaseReduction[bIdx] = capacityBaseReduction - l[bIdx].max();
+      }
+      for (int iIdx = 0; iIdx < nItems; iIdx += 1)
+      {
+          bool iAssigned = bs[iIdx].bin().assigned();
+          int iWeight = bs[iIdx].size();
+          if (iAssigned)
+          {
+            int bIdx = bs[iIdx].bin().val();
+            weightsBaseReduction[bIdx] += iWeight;
+          }
+          else
+          {
+              weightsBaseReduction[nBins + iIdx] = iWeight;
+          }
+      }
+
+      // RMin, RMax
+      int smallestVirtualWeight = std::numeric_limits<int>::max();
+      for (auto bIdx = 0; bIdx < nBins; bIdx += 1)
+      {
+          smallestVirtualWeight = std::min(smallestVirtualWeight, weightsBaseReduction[bIdx]);
+      }
+      deltaReductions[0] = -smallestVirtualWeight;
+      deltaReductions[1] = 0;
+      deltaReductions[2] = capacityBaseReduction - 2 * smallestVirtualWeight + 1;
+  }
+
 }}}
 
 // STATISTICS: int-prop
diff --git a/gecode/int/bin-packing/propagate.hpp b/gecode/int/bin-packing/propagate.hpp
index a8054c1755..eccb7b2bb1 100755
--- a/gecode/int/bin-packing/propagate.hpp
+++ b/gecode/int/bin-packing/propagate.hpp
@@ -208,6 +208,185 @@ namespace Gecode { namespace Int { namespace BinPacking {
     return nosum(s, a, b, ap, bp);
   }
 
+  forceinline int 
+  Pack::fCCM1(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::fMT(int w, int l, int c)
+  {
+    // Conditions
+    int const c0 = w < l;
+    int const c1 = l <= w and 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::fBJ1(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::fVB2Base(int w, int l, int c)
+  {
+      int v = ceil_div_pp(l * w, c);
+      return v > 0 ? v - 1 : 0;
+  }
+
+  forceinline int 
+  Pack::fVB2(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 = fVB2Base(c, l, c);
+    int const t1 = fVB2Base(w, l, c);
+    int const t2 = fVB2Base(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::fFS1(int w, int l, int c)
+  {
+    // Conditions
+    int c0 = w * (l + 1) % c == 0;
+    int c1 = w * (l + 1) % c != 0;
+
+    // Values
+    int v0 = w * l;
+    int v1 = ((w * (l + 1)) / c) * c;
+
+    return (c0 * v0) + (c1 * v1);
+  }
+
+  forceinline int 
+  Pack::fRAD2Base(int w, int l, int c)
+  {
+    // Conditions
+    int const c0 = w < l;
+    int const c1 = l <= w and w <= c - 2 * l;
+    int const c2 = c - 2 * l < w and 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::fRAD2(int w, int l, int c)
+  {
+    // Conditions
+    int const c0 = w < 2 * l;
+    int const c1 = 2 * l <= w;
+
+    // Values
+    int const v0 = fRAD2Base(w, l, c);
+    int const v1 = c - fRAD2Base(c - w, l, c);
+
+    return (c0 * v0) + (c1 * v1);
+  }
+
+  forceinline LambdaRange Pack::lCCM1(int c) {return {1, c / 2};};
+
+  forceinline LambdaRange Pack::lMT(int c)   {return {0, c / 2};}; // The 0 value is included to calculate the L0 bound
+
+  forceinline LambdaRange Pack::lBJ1(int c)  {return {1, c};}
+
+  forceinline LambdaRange Pack::lVB2(int c)  {return {2, c};};
+
+  forceinline LambdaRange Pack::lFS1(int c)  {return {1, 100};};
+
+  forceinline LambdaRange Pack::lRAD2(int c) {return {c / 4 + 1, c / 3};};
+
+  forceinline LambdaRange 
+  Pack::sanitizeLambdaRange(LambdaRange lambdaRange, int nNotZeroWeights, int maxWeight)
+  {
+    if (nNotZeroWeights * maxWeight != 0)
+    {
+        int lMax = std::min(std::numeric_limits<int>::max() / (nNotZeroWeights * maxWeight), lambdaRange.max);
+        return {lambdaRange.min, lMax};
+    }
+    else
+    {
+        return {0,-1};
+    }
+  }
+
+  template<int f(int,int,int)>
+  forceinline int 
+  Pack::calcDffLowerboundSingleLambda(const IntDynamicArray& weights, int capacity, int lambda)
+  {
+      // Transform and sum the weights
+      int sumTransformedWeights = 0;
+      int nWeights = weights.capacity();
+      for (int wIdx = 0; wIdx < nWeights; wIdx += 1)
+      {
+          sumTransformedWeights += f(weights[wIdx], lambda, capacity);
+      }
+
+      // Transform the capacity
+      int transformedCapacity = f(capacity, lambda, capacity);
+
+      // Lowerbound
+      return ceil_div_pp(sumTransformedWeights, transformedCapacity);
+  }
+
+  template<int f(int,int,int), LambdaRange l(int)>
+  forceinline int 
+  Pack::calcDffLowerbound(const IntDynamicArray& weights, int capacity, int nNotZeroWeights, int maxWeight, bool sanitize)
+  {
+      int fLowerbound = 0;
+      LambdaRange lambdaRange = l(capacity);
+      lambdaRange = sanitize ? sanitizeLambdaRange(lambdaRange, nNotZeroWeights, maxWeight) : lambdaRange;
+      int lStep = ceil_div_pp(lambdaRange.max - lambdaRange.min + 1, nLambdaSamples + 1);
+      for (int lambda = lambdaRange.min + lStep; lambda < lambdaRange.max; lambda += lStep)
+      {
+          int lowerbound = calcDffLowerboundSingleLambda<f>(weights, capacity, lambda);
+          fLowerbound = std::max(fLowerbound, lowerbound);
+      }
+      return fLowerbound;
+  }
+
 }}}
 
 // STATISTICS: int-prop
diff --git a/gecode/support/dynamic-array.hpp b/gecode/support/dynamic-array.hpp
index 299fcf03a2..346cadc3c4 100644
--- a/gecode/support/dynamic-array.hpp
+++ b/gecode/support/dynamic-array.hpp
@@ -75,6 +75,9 @@ namespace Gecode { namespace Support {
 
     /// Cast in to pointer of type \a T
     operator T*(void);
+
+    /// Return the capacity of the array
+    int capacity(void) const;
   };
 
 
@@ -154,6 +157,12 @@ namespace Gecode { namespace Support {
     return x;
   }
 
+  template<class T, class A>
+  forceinline int
+  DynamicArray<T,A>::capacity() const {
+    return n;
+  }
+
 }}
 
 // STATISTICS: support-any