Algorithmic Problem Solving: Dynamic Programming and Optimization Techniques

Circular Dynamic Programming

Circular Transportation Problem

In circular DP problems, we can transform a circular road into a linear one by breaking it at the connection between points N and 1, then duplicating the segment to create a path from 1 to 2×n. For each original position i, we now have two positions: i and n+i. This transformation ensures that all point pairs are calculated correctly, and the distance between points i and j (where j ∈ [i-n/2, i] and i ∈ [1, 2×n]) equals either i-j or n+i-j, eliminating the need for case-by-case analysis.

The cost between positions i and j is given by Ai + Aj + i - j. Since Ai + i remains constant, we need to maximize Aj - j. This can be efficiently solved using a monotonic queue optimization: first remove invalid indices from the queue, then use the front element to update the answer, add the current Ai - i to the queue, and maintain the monotonic property.

for(int i = 1; i <= n; i++)  
    values[i] = values[n+i] = readValue();
int doubledSize = n << 1;
queueFront = queueBack = 1;
queue[queueBack] = 0;
for(int i = 1; i <= doubledSize; i++){
    while(queueFront <= queueBack && queue[queueFront] < i - n/2)  
        queueFront++;
    maxResult = max(maxResult, values[i] + values[queue[queueFront]] + i - queue[queueFront]);
    while(values[i] - i >= values[queue[queueBack]] - queue[queueBack] && queueFront <= queueBack) 
        queueBack--;
    queue[++queueBack] = i; 
}
output(maxResult);

Gaussian Elimination

Back-Substitution Gaussian Elimination

This method transforms the matrix into upper triangular form and then performs back-substitution. We must determine whether the system has no solution or infinite solutions by counting the number of pivot elements (r). If any remaining free variables have a non-zero constant term, the system has no solution. Otherwise, it has infinite solutions. Note that "no solution" takes precedence over "infinite solutions" in our analysis.

#define EPSILON 1e-12
int performGaussianElimination(){
    int pivotCount = 1;
    for(int i = 1; i <= n; i++){
        int maxRow = pivotCount;
        for(int j = pivotCount + 1; j <= n; j++)
            if(fabs(matrix[maxRow][i]) < fabs(matrix[j][i]))  
                maxRow = j;
        if(pivotCount != maxRow)  
            swap(matrix[pivotCount], matrix[maxRow]);
        double divisor = matrix[pivotCount][i];
        if(fabs(divisor) < EPSILON)  
            continue;
        for(int j = i; j <= n + 1; j++)  
            matrix[pivotCount][j] /= divisor;
        for(int j = pivotCount + 1; j <= n; j++){
            divisor = matrix[j][i];
            for(int k = i; k <= n + 1; k++)
                matrix[j][k] -= matrix[pivotCount][k] * divisor;
        }
        ++pivotCount;
    }
    for(int i = pivotCount; i <= n; i++){
        if(fabs(matrix[i][n + 1]) > EPSILON)  
            return -1;
    }
    if(pivotCount <= n)  
        return 0;
    solution[n] = matrix[n][n + 1];
    for(int i = n - 1; i >= 1; --i){
        solution[i] = matrix[i][n + 1];
        for(int j = i + 1; j <= n; j++){
            solution[i] -= matrix[i][j] * solution[j];
        }
    }
    return 1;
}

Jordan Elimination

The basic approach is similar to back-substitution, but instead of back-substituting, we transform the matrix into diagonal form, where the diagonal elements directly give the solution. This is achieved by eliminating entire columns during the elimination process.

#define EPSILON 1e-12
int performJordanElimination(){
    int pivotCount = 1;
    for(int i = 1; i <= n; i++){
        int maxRow = pivotCount;
        for(int j = pivotCount + 1; j <= n; j++)
            if(fabs(matrix[maxRow][i]) < fabs(matrix[j][i]))  
                maxRow = j;
        if(maxRow != pivotCount)  
            swap(matrix[maxRow], matrix[pivotCount]);
        if(fabs(matrix[pivotCount][i]) < EPSILON)  
            continue;
        double divisor = matrix[pivotCount][i];
        for(int j = i; j <= n + 1; j++)  
            matrix[pivotCount][j] /= divisor;
        for(int j = 1; j <= n; j++){
            if(j == pivotCount)  
                continue;
            divisor = matrix[j][i];
            for(int k = i; k <= n + 1; k++){
                matrix[j][k] -= matrix[pivotCount][k] * divisor;
            }
        }
        pivotCount++;
    }
    for(int i = pivotCount; i <= n; i++)
        if(fabs(matrix[i][n + 1]) > EPSILON)
            return -1;
    if(pivotCount <= n)  
        return 0;
    for(int i = 1; i <= n; i++)  
        solution[i] = matrix[i][n + 1] / matrix[i][i];
    return 1;
}

Example: Broken Robot

This is a DP problem with after-effects that requires Gaussian elimination for solution. Since the robot can move left or right, the states in each row influence each other, making recursive solution impossible. We use Gaussian elimination to solve the system of equations.

Let dp[i][j] represent the expected number of steps to reach the last row from position (i, j). The base case is dp[n][j] = 0. The DP transitions are:

For j = 1:
dp[i][1] = 1 + (1/3)dp[i][1] + (1/3)dp[i][2] + (1/3)dp[i+1][1]

For 1 < j < m:
dp[i][j] = 1 + (1/4)dp[i][j] + (1/4)dp[i][j+1] + (1/4)dp[i][j-1] + (1/4)dp[i+1][j]

For j = m:
dp[i][m] = 1 + (1/3)dp[i][m] + (1/3)dp[i][m-1] + (1/3)dp[i+1][m]

Rearranging these equations, we get a system of m linear equations. Since row i+1 has already been computed, we treat those terms as constants. This results in a banded matrix, which allows for O(m×d²) complexity where d is the bandwidth (treated as a constant), giving an overall complexiyt of O(m×n).

#include<iostream>
#include<cmath>
#include<algorithm>
using namespace std;
#define EPSILON 1e-7
const int BANDWIDTH = 2;
int rows, cols, startX, startY;
double dp[1010][1010];
double equation[1010][1010];

void solveLinearSystem(int currentRow){
    for(int i = 1; i <= cols; i++){
        if(fabs(equation[i][i]) < EPSILON){
            for(int j = i + 1; j <= min(cols, i + BANDWIDTH); j++)
                if(fabs(equation[j][i]) > fabs(equation[i][i])){
                    swap(equation[j], equation[i]);
                    break;
                }
        }
        if(fabs(equation[i][i]) < EPSILON)  
            continue;
        double divisor = equation[i][i];
        for(int j = i; j <= min(i + 2 * BANDWIDTH, cols); j++)   
            equation[i][j] /= divisor;
        equation[i][cols + 1] /= divisor;
        for(int j = i + 1; j <= min(i + BANDWIDTH, cols); j++){
            divisor = equation[j][i];
            for(int k = i; k <= min(i + 2 * BANDWIDTH, cols); k++)
                equation[j][k] -= equation[i][k] * divisor;
            equation[j][cols + 1] -= equation[i][cols + 1] * divisor;
        }
    }
    dp[currentRow][cols] = equation[cols][cols + 1];
    for(int i = cols - 1; i >= 1; --i){
        dp[currentRow][i] = equation[i][cols + 1];
        for(int j = i + 1; j <= min(i + 2 * BANDWIDTH, cols); j++)
            dp[currentRow][i] -= dp[currentRow][j] * equation[i][j];
    }
}

int main(){
    cin >> rows >> cols >> startX >> startY;
    if(cols == 1){
        printf("%.2lf", (rows - startX) * 2.0);
        return 0;
    }
    for(int i = rows - 1; i >= 1; i--){
        equation[1][1] = equation[cols][cols] = 2.0 / 3;
        equation[1][2] = equation[cols][cols - 1] = -1.0 / 3;
        equation[1][cols + 1] = 1 + dp[i + 1][1] / 3;
        equation[cols][cols + 1] = 1 + dp[i + 1][cols] / 3;
        for(int j = 2; j <= cols - 1; j++){
            equation[j][j - 1] = equation[j][j + 1] = -1.0 / 4;
            equation[j][j] = 3.0 / 4;
            equation[j][cols + 1] = 1 + dp[i + 1][j] / 4;
        }
        solveLinearSystem(i);
    }
    printf("%.10lf", dp[startX][startY]);
    return 0;
}

Monotonic Queue Optimization for Dynamic Programming

This technique significantly reduces complexity by handling decision varible bounds with monotonic queues.

Stock Trading Problem

Let dp[i][j] represent the maximum profit achievable by the end of day i with j stocks in hand. For each day i, we have several options:

  1. Initial purchase: dp[i][j] = -price[i] × j
  2. No action: dp[i][j] = max(dp[i][j], dp[i-1][j])
  3. Buy additional stocks (if i > w): dp[i][j] = max(dp[i][j], dp[i-w-1][k] - (j-k) × price[i]), where k < j and j-k ≤ limit[i]
  4. Sell stocks: Similar to buying but with j enumerated in reverse order
#include<iostream>
#include<algorithm>
using namespace std;

const int MAX_STOCKS = 2010;
int days, maxStocks, cooldown;
int buyPrice[MAX_STOCKS], sellPrice[MAX_STOCKS], buyLimit[MAX_STOCKS], sellLimit[MAX_STOCKS];
int dp[MAX_STOCKS][MAX_STOCKS];
int queue[MAX_STOCKS << 1], queueFront, queueBack;
int maxProfit = -1e9;

int main(){
    cin >> days >> maxStocks >> cooldown;
    for(int i = 1; i <= days; i++)
        cin >> buyPrice[i] >> sellPrice[i] >> buyLimit[i] >> sellLimit[i];
    
    memset(dp, 128, sizeof(dp));
    for(int i = 1; i <= days; i++)
        for(int j = 0; j <= buyLimit[i]; j++)
            dp[i][j] = -buyPrice[i] * j;
            
    for(int i = 1; i <= days; i++){
        for(int j = 0; j <= maxStocks; j++)
            dp[i][j] = max(dp[i][j], dp[i-1][j]);
            
        if(i < cooldown + 1)  
            continue;
            
        queueFront = 1, queueBack = 0;
        for(int j = 0; j <= maxStocks; j++){
            while(queueFront <= queueBack && queue[queueFront] < j - buyLimit[i])  
                queueFront++;
            if(queueFront <= queueBack)  
                dp[i][j] = max(dp[i][j], dp[i-cooldown-1][queue[queueFront]] + queue[queueFront] * buyPrice[i] - j * buyPrice[i]);
            while(queueFront <= queueBack && dp[i-cooldown-1][j] + j * buyPrice[i] >= dp[i-cooldown-1][queue[queueBack]] + queue[queueBack] * buyPrice[i])  
                queueBack--;
            queue[++queueBack] = j;
        }
        
        queueFront = 1, queueBack = 0;
        for(int j = maxStocks; j >= 0; j--){
            while(queueFront <= queueBack && queue[queueFront] > j + sellLimit[i])  
                queueFront++;
            if(queueFront <= queueBack)  
                dp[i][j] = max(dp[i][j], dp[i-cooldown-1][queue[queueFront]] + queue[queueFront] * sellPrice[i] - j * sellPrice[i]);
            while(queueFront <= queueBack && dp[i-cooldown-1][j] + j * sellPrice[i] >= dp[i-cooldown-1][queue[queueBack]] + queue[queueBack] * sellPrice[i])  
                queueBack--;
            queue[++queueBack] = j;
        }
    }
    
    for(int j = 0; j <= maxStocks; j++)  
        maxProfit = max(maxProfit, dp[days][j]);
        
    cout << maxProfit << endl;
    return 0;
}

Buying Feed Problem

This problem demonstrates monotonic queue optimization for the knapsack problem. Let dp[i][j] represent the minimum cost to reach the i-th store with j tons of feed.

  • If no purchase at store i: dp[i][j] = dp[i-1][j] + (x[i] - x[i-1]) × j²
  • If purchasing at store i: dp[i][j] = min(dp[i-1][k] + (x[i] - x[i-1]) × k² + (j-k) × cost[i]), where k < j and j-k ≤ inventory[i]
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;

const int MAX_STORES = 510;
ll targetAmount, homePosition, numStores;
struct Store {
    ll position, inventory, cost;
} stores[MAX_STORES];
ll dp[MAX_STORES][10010];
ll queue[MAX_STORES << 1], queueFront, queueBack;

ll calculateValue(int storeIdx, int prevAmount) {
    return dp[storeIdx-1][prevAmount] + 
           (stores[storeIdx].position - stores[storeIdx-1].position) * prevAmount * prevAmount - 
           prevAmount * stores[storeIdx].cost;
}

int main(){
    cin >> targetAmount >> homePosition >> numStores;
    for(int i = 1; i <= numStores; i++)  
        cin >> stores[i].position >> stores[i].inventory >> stores[i].cost;
    
    sort(stores + 1, stores + numStores + 1, [](Store a, Store b) {
        return a.position < b.position;
    });
    
    memset(dp, 0x3f, sizeof(dp));
    for(int j = 0; j <= stores[1].inventory; j++)
        dp[1][j] = stores[1].cost * j;
        
    for(int i = 2; i <= numStores; i++){
        dp[i][0] = 0;
        queueFront = 1, queueBack = 0;
        queue[++queueBack] = 0;
        for(int j = 1; j <= targetAmount; j++){
            dp[i][j] = min(dp[i][j], dp[i-1][j] + 
                          (stores[i].position - stores[i-1].position) * j * j);
            while(queueFront <= queueBack && queue[queueFront] < j - stores[i].inventory)  
                queueFront++;
            if(queueFront <= queueBack)  
                dp[i][j] = min(dp[i][j], calculateValue(i, queue[queueFront]) + j * stores[i].cost);
            while(queueFront <= queueBack && calculateValue(i, queue[queueBack]) >= calculateValue(i, j))  
                queueBack--;
            queue[++queueBack] = j;
        }
    }
    
    cout << dp[numStores][targetAmount] + (homePosition - stores[numStores].position) * targetAmount * targetAmount << endl;
    return 0;
}

Beautiful Waltz Problem

This interesting problem can be solved using monotonic queue optimization. The O(N×M×T) approach considers using or not using magic for each time period, with separate handling for each of the 4 directions. However, we can optimize this to O(N×M×K) by recognizing that using magic throughout an entire time period is equivalent to using it for individual times within that period.

Let dp[i][j][t] represent the maximum sliding distance when the piano slides to position (i, j) at the end of time period t. For upward sliding (d[t] = 1), we have:

dp[i][j][t] = max(dp[i][j][t-1], dp[k][j][t-1] + k - j), where k > i and k - i ≤ s[t]

Similar formulas apply for other directions. We can use monotonic queues to optimize these transitions.

#include<iostream>
#include<algorithm>
using namespace std;

const int MAX_SIZE = 210;
int rows, cols, startX, startY, numPeriods;
bool grid[MAX_SIZE][MAX_SIZE];
int direction[MAX_SIZE], duration[MAX_SIZE];
int dp[MAX_SIZE][MAX_SIZE][MAX_SIZE];
int queue[MAX_SIZE << 1], queueFront, queueBack;
int maxDistance;

int main(){
    cin >> rows >> cols >> startX >> startY >> numPeriods;
    char c;
    for(int i = 1; i <= rows; i++){
        for(int j = 1; j <= cols; j++){
            cin >> c;
            grid[i][j] = (c == '.' ? 0 : 1);
        }
        cin >> c;
    }
    for(int start, end, i = 1; i <= numPeriods; i++){
        cin >> start >> end >> direction[i];
        duration[i] = end - start + 1;
    }
    
    memset(dp, 128, sizeof(dp));
    for(int t = 0; t <= numPeriods; t++)  
        dp[startX][startY][t] = 0;
        
    for(int t = 1; t <= numPeriods; t++){
        if(direction[t] == 1){ // Upward
            for(int j = 1; j <= cols; j++){
                queueFront = 1, queueBack = 0;
                for(int i = rows; i >= 1; --i){
                    if(grid[i][j]){ queueFront = 1, queueBack = 0; continue; }
                    while(queueFront <= queueBack && queue[queueFront] > i + duration[t])  
                        queueFront++;
                    dp[i][j][t] = max(dp[i][j][t], dp[i][j][t-1]);
                    if(queueFront <= queueBack)  
                        dp[i][j][t] = max(dp[i][j][t], dp[queue[queueFront]][j][t-1] + queue[queueFront] - i);
                    while(queueFront <= queueBack && dp[i][j][t-1] + i >= dp[queue[queueBack]][j][t-1] + queue[queueBack]) 
                        queueBack--;
                    queue[++queueBack] = i;
                    maxDistance = max(maxDistance, dp[i][j][t]);
                }
            }
        }
        if(direction[t] == 2){ // Downward
            for(int j = 1; j <= cols; j++){
                queueFront = 1, queueBack = 0;
                for(int i = 1; i <= rows; i++){
                    if(grid[i][j]){ queueFront = 1, queueBack = 0; continue; }
                    while(queueFront <= queueBack && queue[queueFront] < i - duration[t])  
                        queueFront++;
                    dp[i][j][t] = max(dp[i][j][t], dp[i][j][t-1]);
                    if(queueFront <= queueBack)  
                        dp[i][j][t] = max(dp[i][j][t], dp[queue[queueFront]][j][t-1] + i - queue[queueFront]);
                    while(queueFront <= queueBack && dp[i][j][t-1] - i >= dp[queue[queueBack]][j][t-1] - queue[queueBack])  
                        queueBack--;
                    queue[++queueBack] = i;
                    maxDistance = max(maxDistance, dp[i][j][t]);
                }
            }
        }
        if(direction[t] == 3){ // Left
            for(int i = 1; i <= rows; i++){
                queueFront = 1, queueBack = 0;
                for(int j = cols; j >= 1; --j){
                    if(grid[i][j]){ queueFront = 1, queueBack = 0; continue; }
                    while(queueFront <= queueBack && queue[queueFront] > j + duration[t])  
                        queueFront++;
                    dp[i][j][t] = max(dp[i][j][t], dp[i][j][t-1]);
                    if(queueFront <= queueBack)  
                        dp[i][j][t] = max(dp[i][j][t], dp[i][queue[queueFront]][t-1] + queue[queueFront] - j);
                    while(queueFront <= queueBack && dp[i][j][t-1] + j >= dp[i][queue[queueBack]][t-1] + queue[queueBack])  
                        queueBack--;
                    queue[++queueBack] = j;
                    maxDistance = max(maxDistance, dp[i][j][t]);
                }
            }
        }
        if(direction[t] == 4){ // Right
            for(int i = 1; i <= rows; i++){
                queueFront = 1, queueBack = 0;
                for(int j = 1; j <= cols; j++){
                    if(grid[i][j]){ queueFront = 1, queueBack = 0; continue; }
                    while(queueFront <= queueBack && queue[queueFront] < j - duration[t])  
                        queueFront++;
                    dp[i][j][t] = max(dp[i][j][t], dp[i][j][t-1]);
                    if(queueFront <= queueBack)  
                        dp[i][j][t] = max(dp[i][j][t], dp[i][queue[queueFront]][t-1] + j - queue[queueFront]);
                    while(queueFront <= queueBack && dp[i][j][t-1] - j >= dp[i][queue[queueBack]][t-1] - queue[queueBack])  
                        queueBack--;
                    queue[++queueBack] = j;
                    maxDistance = max(maxDistance, dp[i][j][t]);
                }
            }
        }
    }
    
    cout << maxDistance << endl;
    return 0;
}

Data Structure Optimizaton for Dynamic Programming

This category of DP problems is relatively straightforward: after deriving the brute-force recurrence, we select an appropriate data structure for optimization.

The Battle of Chibi Problem

This problem requires counting the number of increasing subsequences of length m. Let's first consider the naive DP approach: let dp[i][j] represent the number of increasing subsequences ending at position i with length j. The recurrence is:

dp[i][j] = Σ[k=1 to i-1] [a[k] < a[i]] × dp[k][j-1]

This is a 2-dimensional partial order problem that can be efficiently solved using a Fenwick tree (Binary Indexed Tree) for optimization.

#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;

#define lowbit(x) ((x) & -(x))
const int MAX_SIZE = 1010;
const int MOD = 1e9 + 7;
ll dp[MAX_SIZE][MAX_SIZE];
int sequenceLength, subseqLength;
int values[MAX_SIZE];
int sortedValues[MAX_SIZE], uniqueCount;
ll fenwickTree[MAX_SIZE];

void updateFenwick(int index, ll value) {
    for(int i = index; i <= uniqueCount; i += lowbit(i))
        fenwickTree[i] = (fenwickTree[i] + value) % MOD;
}

ll queryFenwick(int index) {
    ll result = 0;
    for(int i = index; i > 0; i -= lowbit(i))
        result = (result + fenwickTree[i]) % MOD;
    return result;
}

int main(){
    int testCases;
    cin >> testCases;
    for(int testCase = 1; testCase <= testCases; testCase++){
        ll totalSubsequences = 0;
        memset(dp, 0, sizeof(dp));
        cin >> sequenceLength >> subseqLength;
        for(int i = 1; i <= sequenceLength; i++)  
            sortedValues[i] = values[i] = [&](){
                int val; cin >> val; return val;
            }();
            
        sort(sortedValues + 1, sortedValues + sequenceLength + 1);
        uniqueCount = unique(sortedValues + 1, sortedValues + sequenceLength + 1) - sortedValues - 1;
        
        for(int i = 1; i <= sequenceLength; i++)  
            values[i] = lower_bound(sortedValues + 1, sortedValues + uniqueCount + 1, values[i]) - sortedValues;
            
        for(int i = 1; i <= sequenceLength; i++)
            dp[i][1] = 1;
            
        for(int length = 2; length <= subseqLength; length++){
            memset(fenwickTree, 0, sizeof(fenwickTree));
            for(int i = length - 1; i <= sequenceLength; i++){
                dp[i][length] = queryFenwick(values[i] - 1) % MOD;
                updateFenwick(values[i], dp[i][length - 1]);
            }
        }
        
        for(int i = subseqLength; i <= sequenceLength; i++)  
            totalSubsequences = (totalSubsequences + dp[i][subseqLength]) % MOD;
            
        printf("Case #%d: %lld\n", testCase, totalSubsequences);
    }
    return 0;
}

Doubling Optimization for Dynamic Programming

This technique is typically used for problems that can be arbitrarily partitioned, employing binary decomposition to fill values.

Driving Travel Problem

For this problem, the number of days can be arbitrarily partitioned. We first preprocess the next destinations for both drivers A and B at each city using a set data structure. Then we define:

  • dp[i][j][k]: city reached after driving 2^i days from city j, with driver k starting
  • da[i][j][k]: distance driven by A after 2^i days from city j, with driver k starting
  • db[i][j][k]: distance driven by B after 2^i days from city j, with driver k starting

For the first question, we enumerate each city and use a binary approach similar to finding LCA (Lowest Common Ancestor) to determine the optimal starting city. The second question follows a similar approach.

#include<iostream>
#include<set>
#include<algorithm>
using namespace std;

typedef pair<int, int> pii;
#define make_pair mp
const int MAX_CITIES = 1e5 + 10;
int numCities, numQueries;
long long elevation[MAX_CITIES];
int nextA[MAX_CITIES], nextB[MAX_CITIES];

void preprocessNextCities() {
    set<pii> citySet;
    elevation[0] = -2e9 - 1e8;
    elevation[numCities + 1] = 2e9 + 1e8;
    set<pii>::iterator it1, it2, it3;
    citySet.insert(mp(elevation[numCities + 1], numCities + 1));
    citySet.insert(mp(elevation[numCities], numCities));
    citySet.insert(mp(elevation[numCities - 1], numCities - 1));
    nextB[numCities - 1] = numCities;
    
    for(int i = numCities - 2; i >= 1; --i) {
        it1 = it2 = citySet.upper_bound(mp(elevation[i], i));
        it1--;
        int prev1 = (*it1).second;
        int next1 = (*it2).second;
        it1--; it2++;
        int prev2 = (*it1).second;
        int next2 = (*it2).second;
        
        if(elevation[i] - elevation[prev1] <= elevation[next1] - elevation[i] && prev1) {
            nextB[i] = prev1; 
            nextA[i] = (elevation[i] - elevation[prev2] <= elevation[next1] - elevation[i]) ? prev2 : next1;
        }
        else {
            nextB[i] = next1;
            nextA[i] = (elevation[i] - elevation[prev1] <= elevation[next2] - elevation[i]) ? prev1 : next2;
        }
        citySet.insert(mp(elevation[i], i));
    }
}

int dp[25][MAX_CITIES][2];
int distA[25][MAX_CITIES][2], distB[25][MAX_CITIES][2];
int bestStartCity;
double minRatio = 1e9;

int main(){
    cin >> numCities;
    for(int i = 1; i <= numCities; i++)  
        cin >> elevation[i];
    preprocessNextCities();
    
    for(int j = 1; j <= numCities; j++) {
        dp[0][j][0] = nextA[j];
        dp[0][j][1] = nextB[j];
    }
    
    for(int j = 1; j <= numCities; j++)
        for(int k = 0; k <= 1; k++)
            dp[1][j][k] = dp[0][dp[0][j][k]][k ^ 1];
            
    for(int i = 2; i <= 20; i++) {
        for(int j = 1; j <= numCities; j++) {
            for(int k = 0; k <= 1; k++)
                dp[i][j][k] = dp[i - 1][dp[i - 1][j][k]][k];
        }
    }
    
    for(int j = 1; j <= numCities; j++) {
        if(nextA[j])
            distA[0][j][0] = abs(elevation[j] - elevation[nextA[j]]);
        if(nextB[j])
            distB[0][j][1] = abs(elevation[j] - elevation[nextB[j]]);
    }
    
    for(int j = 1; j <= numCities; j++) {
        for(int k = 0; k <= 1; k++) { 
            distA[1][j][k] = distA[0][j][k] + distA[0][dp[0][j][k]][k ^ 1];
            distB[1][j][k] = distB[0][j][k] + distB[0][dp[0][j][k]][k ^ 1];
        }
    }
    
    for(int i = 2; i <= 20; i++) {
        for(int j = 1; j <= numCities; j++) {
            for(int k = 0; k <= 1; k++) {
                distA[i][j][k] = distA[i - 1][j][k] + distA[i - 1][dp[i - 1][j][k]][k];
                distB[i][j][k] = distB[i - 1][j][k] + distB[i - 1][dp[i - 1][j][k]][k];
            }
        }
    }
    
    int maxDays;
    cin >> maxDays;
    for(int startCity = 1; startCity <= numCities; startCity++) {
        int currentCity = startCity;
        int remainingDays = maxDays;
        int totalDistA = 0, totalDistB = 0;
        
        for(int k = 18; k >= 0; --k) {
            if(!dp[k][currentCity][0])  
                continue;
            if(maxDays >= totalDistA + totalDistB + distA[k][currentCity][0] + distB[k][currentCity][0]) {
                totalDistA += distA[k][currentCity][0];
                totalDistB += distB[k][currentCity][0];
                currentCity = dp[k][currentCity][0];
            }
        }
        
        if(totalDistB != 0) {
            double currentRatio = (double)totalDistA / totalDistB;
            if(currentRatio == minRatio)  
                bestStartCity = (elevation[startCity] > elevation[bestStartCity]) ? startCity : bestStartCity;
            if(currentRatio < minRatio)  
                bestStartCity = startCity, minRatio = currentRatio;
        }
    }
    
    cout << bestStartCity << endl;
    
    for(numQueries; numQueries; --numQueries) {
        int startCity, maxDays;
        cin >> startCity >> maxDays;
        int totalDistA = 0, totalDistB = 0;
        
        for(int k = 18; k >= 0; --k) {
            if(!dp[k][startCity][0])  
                continue;
            if(maxDays >= totalDistA + totalDistB + distA[k][startCity][0] + distB[k][startCity][0]) {
                totalDistA += distA[k][startCity][0];
                totalDistB += distB[k][startCity][0];
                startCity = dp[k][startCity][0];
            }
        }
        
        cout << totalDistA << " " << totalDistB << endl;
    }
    return 0;
}

Slope Optimization for Dynamic Programming

Task Scheduling Problem

The first naive O(n²) approach is straightforward. For the O(n) solution, we consider a different DP formulation: let dp[i] represent the minimum cost for the first i tasks. The recurrence becomes:

dp[i] = min(dp[j] + T[i] × (C[i] - C[j]) + s × (C[n] - C[j]))

where T[i] and C[i] are prefix sums. This formulation accounts for the startup cost s for all subsequent groups. Rearranging the equation:

dp[i] = min(dp[j] - (T[i] + s) × C[j]) + T[i] × C[i] + s × C[n]

Assuming we find an optimal decision j, we have:

dp[j] = (T[i] + s) × C[j] + dp[i] - T[i] × C[i] - s × C[n]

This forms a linear equation y = kx + b, where y = dp[j], x = C[j], and k = T[i] + s. We can maintain a convex hull using a monotonic queue, and since the slope T[i] + s is monotonically increasing, we only need to maintain the lower-left endpoint of the convex hull.

#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;

const int MAX_TASKS = 3e5 + 10;
int numTasks, startupCost;
int taskTime[MAX_TASKS], timePrefix[MAX_TASKS];
int taskCost[MAX_TASKS], costPrefix[MAX_TASKS];
ll dp[MAX_TASKS];
int queue[MAX_TASKS << 1], queueFront, queueBack;

int main(){
    cin >> numTasks >> startupCost;
    for(int i = 1; i <= numTasks; i++) {
        cin >> taskTime[i] >> taskCost[i];
        timePrefix[i] = timePrefix[i - 1] + taskTime[i];
        costPrefix[i] = costPrefix[i - 1] + taskCost[i];
    }
    
    memset(dp, 0x3f, sizeof(dp));
    queueFront = 1, queueBack = 0;
    dp[0] = 0;
    queue[++queueBack] = 0;
    
    for(int i = 1; i <= numTasks; i++) {
        while(queueFront < queueBack && 
              (double)(dp[queue[queueFront]] - dp[queue[queueFront + 1]]) / 
              (costPrefix[queue[queueFront]] - costPrefix[queue[queueFront + 1]]) < timePrefix[i] + startupCost)  
            queueFront++;
            
        dp[i] = dp[queue[queueFront]] - (timePrefix[i] + startupCost) * costPrefix[queue[queueFront]] + 
                timePrefix[i] * costPrefix[i] + startupCost * costPrefix[numTasks];
                
        while(queueFront < queueBack && 
              (double)(dp[queue[queueBack]] - dp[queue[queueBack - 1]]) / 
              (costPrefix[queue[queueBack]] - costPrefix[queue[queueBack - 1]]) >= 
              (double)(dp[i] - dp[queue[queueBack]]) / (costPrefix[i] - costPrefix[queue[queueBack]]))  
            queueBack--;
            
        queue[++queueBack] = i;
    }
    
    cout << dp[numTasks] << endl;
    return 0;
}

Cats Transport Problem

For each cat i, we define A[i] as the earliest time a person must depart to pick up cat i. This can be calculated as A[i] = T[i] - Σ[j=1 to h[i]] D[j]. After sorting A to ensure monotonicity and computing prefix sums S, we define dp[i][j] as the minimum waiting time for the first i people to pick up j cats. The recurrence is:

dp[i][j] = min(dp[i-1][k] + Σ[p=k+1 to j] A[j] - A[p])

Simplifying the summation and rearranging, we get:

dp[i-1][k] + S[k] = A[j] × k + dp[i][j] - A[j] × j + S[j]

Since we've sorted the values, the slope A[j] is monotonically increasing, allowing us to use the same monotonic queue approach as in the previous problem.

#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;

const int MAX_CATS = 2e5 + 10;
int numCats, numPeople, numGroups;
int distances[MAX_CATS];
int hillPositions[MAX_CATS], arrivalTimes[MAX_CATS];
ll earliestDeparture[MAX_CATS], prefixSum[MAX_CATS];
ll dp[110][MAX_CATS];
int queue[MAX_CATS << 1], queueFront, queueBack;

double calculateSlope(int groupIndex, int x, int y) {
    double numerator = ((dp[groupIndex-1][x] + prefixSum[x]) - (dp[groupIndex-1][y] + prefixSum[y]));
    double denominator = x - y;
    if(denominator == 0) return numerator > 0 ? 1e18 : -1e18;
    return numerator / denominator;
}

int main(){
    cin >> numCats >> numPeople >> numGroups;
    for(int i = 2; i <= numCats; i++) {
        int d;
        cin >> d;
        distances[i] = distances[i - 1] + d;
    }
    for(int i = 1; i <= numPeople; i++) {
        cin >> hillPositions[i] >> arrivalTimes[i];
        earliestDeparture[i] = arrivalTimes[i] - distances[hillPositions[i]];
    }
    
    sort(earliestDeparture + 1, earliestDeparture + numPeople + 1);
    for(int i = 1; i <= numPeople; i++)
        prefixSum[i] = prefixSum[i - 1] + earliestDeparture[i];
        
    memset(dp, 0x3f, sizeof(dp));
    for(int j = 1; j <= numPeople; j++) {
        dp[1][j] = earliestDeparture[j] * j - prefixSum[j];
    }
    
    for(int i = 2; i <= numGroups; i++) {
        dp[i][0] = 0;
        queueFront = 1, queueBack = 0;
        queue[++queueBack] = 0;
        for(int j = 1; j <= numPeople; j++) {
            while(queueFront < queueBack && 
                  calculateSlope(i, queue[queueFront + 1], queue[queueFront]) <= earliestDeparture[j])  
                queueFront++;
                
            int k = queue[queueFront];
            dp[i][j] = dp[i-1][k] + prefixSum[k] - earliestDeparture[j] * k + 
                      earliestDeparture[j] * j - prefixSum[j];
                      
            while(queueFront < queueBack && 
                  calculateSlope(i, queue[queueBack], queue[queueBack - 1]) >= 
                  calculateSlope(i, j, queue[queueBack]))  
                queueBack--;
                
            queue[++queueBack] = j;
        }
    }
    
    cout << dp[numGroups][numPeople] << endl;
    return 0;
}

Posted on Sat, 15 Aug 2026 16:08:49 +0000 by Shaun