Laborator 4: Programare Dinamică (continuare)

$n = 3$

Răspuns: 64 (înmulțiri scalare)

Explicație: Avem 3 matrice:

• A de dimensiuni (2, 3)
• B de (3, 4)
• C de (4, 5)

În funcție de ordinea efectuării înmulțirilor matriceale, numărul total de înmulțiri scalare poate să fie foarte diferit:

• $(AB)C$ ⇒ $24 + 40 = 64$ de înmulțiri
  • explicație: $X = (AB)$ generează $2 * 3 * 4 = 24$ înmulțiri, $(XC)$ generează $2 * 4 * 5 = 40$ de înmulțiri
• $A(BC)$ ⇒ $60 + 30 = 90$ de înmulțiri
  • explicație: $X =(BC)$ generează $3 * 4 * 5 = 60$ înmulțiri, $(AX)$ generează $2 * 3 * 5 = 30$ de înmulțiri

Rezultatul optim se obține pentru cea de a treia parantezare: $(AB)C$.

$n = 4$

Răspuns: 48 (înmulțiri scalare)

Explicație: Avem 4 matrice:

• A de dimensiuni (2, 3)
• B de (3, 4)
• C de (4, 2)
• D de (2, 3)

În funcție de ordinea efectuării înmulțirilor matriceale, numărul total de înmulțiri scalare poate să fie foarte diferit:

• $(AB)C)D$ ⇒ $24 + 16 + 12 = 52$ înmulțiri
  • explicație: $X = (AB)$ generează $2 * 3 *4 = 24$ înmulțiri scalare, $Y = (XC)$ generează $2 * 4 * 2 = 16$ înmulțiri scalare, $Z = YD$ generează $2 * 2 *3 = 12$ înmulțiri scalare
• $(A(BC))D$ ⇒ $24 + 12 + 12 = 48$ înmulțiri
  • explicație: $X = (BC)$ generează $3 * 4 * 2 = 24$ înmulțiri scalare, $Y = (AX)$ generează $ 2 * 3 * 2= 12$ înmulțiri scalare, $Z = YD$ generează 2 * 2 * 3$ = 12$ înmulțiri scalare
• $(AB)(CD)$ ⇒ $ = $ inmulțiri
  • explicație: $X = (AB)$ generează $2 * 3 * 4 = 24$ înmulțiri scalare, $Y = (CD)$ generează $4 * 2 * 3 = 24$ înmulțiri scalare, $Z = XY$ generează $2 * 4 * 3 = 24$ înmulțiri scalare
• $A((BC)D)$ ⇒ $24 + 18 + 27 = 69$ înmulțiri
  • explicație: $X = (BC)$ generează $3 * 4 * 2 = 24$ înmulțiri scalare, $Y = (XD)$ generează $3 * 2 * 3 = 18$ înmulțiri scalare, $Z = AY$ generează $3 * 3 * 3 = 27$ înmulțiri scalare
• $A(B(CD))$ ⇒ $24 + 36 + 18 = 78$ înmulțiri
  • explicație: $X = (CD)$ generează $4 * 2 * 3 = 24$ înmulțiri scalare, $Y = (BX)$ generează $3 * 4 * 3 = 36$ înmulțiri scalare, $Z = AY$ generează $2 * 3 * 3 = 18$ înmulțiri scalare

Rezultatul optim se obține pentru cea de a treia parantezare: $((A(BC))D)$.

$n = 4$

Răspuns: 2856 (înmulțiri scalare)

Explicație: Avem 4 matrice:

• A de dimensiuni (13, 5)
• B de (5, 89)
• C de (89, 3)
• D de (3, 34)

În funcție de ordinea efectuării înmulțirilor matriciale, numărul total de înmulțiri scalare poate să fie foarte diferit:

• $((AB)C)D$ ⇒ 10582 înmulțiri
• $(AB)(CD)$ ⇒ 54201 înmulțiri
• $(A(BC))D$ ⇒ 2856 înmulțiri
• $A((BC)D)$ ⇒ 4055 înmulțiri
• …

Rezultatul optim se obține pentru cea de a treia parantezare: $(A(BC))D$.

// kInf este valoarea maximă - "infinitul" nostru
const unsigned long long kInf = std::numeric_limits<unsigned long long>::max();
 
// T = O(n ^ 3) 
// S = O(n ^ 2) - stocăm n x n întregi în tabloul dp
 unsigned long long solve_podm(int n, const vector<int> &d) {
    // dp[i][j] = numărul MINIM înmulțiri scalare cu codare, poate fi calculat produsul
    //            matriceal M_i * M_i+1 * ... * M_j
    vector<vector<unsigned long long>>  dp(n + 1, vector<unsigned long long> (n + 1, kInf));
 
    // Cazul de bază 1: nu am ce înmulți 
    for (int i = 1; i <= n; ++i) {
        dp[i][i] = 0ULL;  // 0 pe unsigned long long (voi folosi mai încolo și 1ULL)
    }
 
    // Cazul de bază 2: matrice d[i - 1] x d[i] înmulțită cu matrice d[i] x d[i + 1] 
    // (matrice pe poziții consecutive)
    for (int i = 1; i < n; ++i) {
        dp[i][i + 1] =  1ULL * d[i - 1] * d[i] * d[i + 1];  
    }
 
    // Cazul general:
    // dp[i][j] = min(dp[i][k] + dp[k + 1][j] + d[i - 1] * d[k] * d[j]), k = i : j - 1
    for (int len = 2; len <= n; ++len) {            // fixăm lungimea intervalului (2, 3, 4, ...)
        for (int i = 1; i + len - 1 <= n; ++i) {    // fixăm capătul din stânga: i
            int j = i + len - 1;                    // capătul din dreapta se deduce: j
 
            // Iterăm prin indicii dintre capete, spărgând șirul de înmulțiri in două (paranteze).
            for (int k = i; k < j; ++k) {
                // M_i * ... M_j  = (M_i * .. * M_k) * (M_k+1 *... * M_j)
                unsigned long long new_sol = dp[i][k] + dp[k + 1][j] + 1ULL * d[i - 1] * d[k] * d[j];
 
                // actualizăm soluția dacă este mai bună
                dp[i][j] = min(dp[i][j], new_sol); 
            }
        }
    }
 
    // Rezultatul se află în dp[1][n]: Numărul MINIM de inmultiri scalare
    // pe care trebuie să le facem pentru a obține produsul M_1 * ... * M_n
    return dp[1][n];
 
}

• definiție : b este inversul modular al lui a în raport cu MOD dacă $ a * b = 1 (modulo \ MOD)$
• utilizare : $ \frac{a}{b} \ \% \ MOD = ((a \ \% \ MOD) * (invers(b) \ \% \ MOD)) \ \% \ MOD $
• calculare : deoarece la PA această discuție are sens doar în contextul posibilității implementării unei recurențe DP în care folosim resturile doar pentru a evita overflow/imposibilitatea de a reține rezultatul pe tipurile standard de tip int (adică nu ne interesează să dăm o metoda generală pentru invers modular), vom simplifica problema - MOD este prim!!!
• Mica teoremă a lui Fermat: Dacă p este un număr prim și a este un număr întreg care nu este multiplu al lui p, atunci $a^{p-1} = 1 (modulo \ p)$.
  • din definiția inversului modular, reiese că a și b nu sunt multipli ai lui MOD
  • introducând notațiile noastre în teoremă și prelucrând obținem
    • $a ^ {MOD - 1} = 1 (modulo \ MOD) <=> a * (a ^{MOD-2}) = 1 (modulo \ MOD)$
    • deci, inversul modular al lui a (în aceste condiții specifice) este $b = a ^ {MOD -2 }$

$n = 1$ sau $n = 2$ sau $n = 3$

Răspuns: 1 (un singur gard)

Explicație: Se poate forma un singur gard în fiecare caz, după cum este ilustrat și în figura Garduri_123.

$n = 4$

Răspuns: 2

Explicație: Se pot forma 2 garduri, în funcție de cum așezăm piesele, după cum este ilustrat și în figura Garduri_4. Observăm că de fiecare dată când punem o piesă în poziție orizontală, de fapt suntem obligați să punem 4 piese, una peste alta!

$n = 5$

Răspuns: 3

Explicație: Se pot forma 3 garduri, în funcție de cum așezăm piesele, după cum este ilustrat și în figura Garduri_5.

• dacă dorim ca acest gard să se termine cu 4 piese în poziție orizontală (una peste alta - marcat cu roșu), atunci la stânga mai ramane de completat un subgard de lungime 1, în toate modurile posibile
• dacă dorim ca acest gard să se termine cu o piesă în poziție verticală (marcat cu roșu), atunci la stânga mai rămâne de completat un subgard de lungime 4, în toate modurile posibile

$n = 6$

Răspuns: 4

Explicație: Se pot forma 4 garduri, în funcție de cum așezăm piesele, după cum este ilustrat și în figura Garduri_6.

• dacă dorim ca acest gard să se termine cu o piesă în poziție verticală (marcat cu roșu), atunci la stânga mai rămâne de completat un subgard de lungime 5, în toate modurile posibile
• dacă dorim ca acest gard să se termine cu 4 piese în poziție orizontală (una peste alta - marcat cu roșu), atunci la stânga mai ramane de completat un subgard de lungime 2, în toate modurile posibile

#define MOD 1009
int gardurile_lui_Gigel(int n) {
    // cazurile de bază
    if (n <= 3) return 1;
    if (n == 4) return 2;
 
    vector<int> dp(n + 1); // păstrez indexarea de la 1 ca în explicații
 
    // cazurile de bază
    dp[1] = dp[2] = dp[3] = 1;
    dp[4] = 2;
 
    // cazul general
    for (int i = 5; i <= n; ++i) {
        dp[i] = (dp[i - 1] + dp[i - 4]) % MOD;
    }
 
    return dp[n];
}

Menționez că am folosit expresia $dp[i] = (dp[i - 1] + dp[i - 4]) \ \% \ MOD$ în loc de $dp[i] = ((dp[i - 1] \ \% \ MOD) + (dp[i - 4] \ \% \ MOD)) \ \% \ MOD$, deoarece, pe valorile anterior calculate în dp, a fost deja aplicată operația $%$.

Am plecat cu numerele $1, 1, 1, 2$ și, la fiecare pas, rezultatul stocat este $\ \% \ MOD$, deci, tot ce este stocat deja în dp este un rest în raport cu MOD. NU mai era nevoie, deci, să aplicăm % și pe termenii din paranteză.

Presupunand ca nu mai exista alte specificatii ale problemei si ca avand cele KMAX cazuri de baza (primele KMAX valori ar trebui stiute/deduse prin alte reguli), atunci un algoritm poate implementa recurenta de mai sus folosind 2 cicluri de tip for (for i = 1 : n, for k = 1 : KMAX …).

• complexitatea temporala : $ T =O(n * KMAX) = O(n) $
  • reamintim ca acea valoarea KMAX este o constanta foarte mica in compartie cu n (ex. KMAX < 100)
• complexitatea spatiala : $ S = O(n) $

• am presupus ca avem nevoie sa retinem doar tabloul dp

• ultima coloana contine toti coeficientii $c_k$ intrucat $dp[i] = \sum_{k=1}^{KMAX} c_k * dp[i - k]$
• celelalte coloane contin doar cate o valoare nenula
  • pe coloana j vom avea valoarea 1 pe linia $j+1$ ($j = 1 : KMAX - 1$)
    • cum obtinem, de exemplu, $dp[i - 1]$?
    • pai avem $dp[i-1]$ chiar si in starea $S_{i-1}$, deci trebuie sa il copiam in starea $S_i$
      • copierea se realizeaza prin inmultirea cu 1
      • daca $dp[i-1]$ era pe ultima pozitiei (pozitia k) in starea $S_{i-1}$, in noua stare $S_i$ este pe penultima pozitie (pozitia $k-1$)
        • deci s-a deplasat la stanga cu o pozitie!
  • in noua stare, noua pozitie este deplasata cu o unitate la stanga fata de starea precedenta
    • de aceea pe coloana $j$, vrem sa avem elementul 1 pe linia $j + 1$ ($j = 1 : KMAX - 1$)
    • cand inmultim $S_{i-1}$ cu coloana $C_j$ dorim sa
      • ce copiam?
        • valoarea $dp[i - KMAX + j]$ din $S_{i-1}$ in $S_{i}$
        • adica sa copiam a j-a valoare de pe linie
      • unde copiam?
        • de pe pozitia $j + 1$ pe pozitia $j$

Mai jos se afla o implementare simplista in C++ care cuprinde toate etapele pe care trebuie sa le realizati in cod, dupa ce stiti cum arata recurenta sub forma matriceala.

#define MOD  1009
#define KMAX 4
 
// C = A * B
void multiply_matrix(int A[KMAX][KMAX], int B[KMAX][KMAX], int C[KMAX][KMAX]) {
    int tmp[KMAX][KMAX];
 
    // tmp = A * B
    for (int i = 0; i < KMAX; ++i) {
        for (int j = 0; j < KMAX; ++j) {
            unsigned long long sum = 0; // presupun ca suma  incape pe 64 de biti
 
            for (int k = 0; k < KMAX; ++k) {
                sum += 1LL * A[i][k] * B[k][j];
            }
 
            tmp[i][j] = sum % MOD;
        }
    }
 
    //  C = tmp
    memcpy(C, tmp, sizeof(tmp));
}
 
// R = C^p
void power_matrix(int C[KMAX][KMAX], int p, int R[KMAX][KMAX]) {
    // tmp = I (matricea identitate)
    int tmp[KMAX][KMAX];
    for (int i = 0; i < KMAX; ++i) {
        for (int j = 0; j < KMAX; ++j) {
            tmp[i][j] = (i == j) ? 1 : 0;
        }
    }
 
    while (p != 1) {
        if  (p % 2 == 0) {
            multiply_matrix(C, C, C);     // C = C*C
            p /= 2;                       // ramane de calculat C^(p/2)
        } else {
            // reduc la cazul anterior:
            multiply_matrix(tmp, C, tmp); // tmp = tmp*C
            --p;                          // ramane de calculat C^(p-1)
        }
    }
 
    // avem o parte din rezultat in C si o parte in tmp
    multiply_matrix(C, tmp, R);           // rezultat = tmp * C
}
 
int garduri_rapide(int n) {
    // cazurile de baza
    if (n <= 3) return 1;
    if (n == 4) return 2;
 
    // construiesc matricea C
    int C[KMAX][KMAX] = { {0, 0, 0, 1},
                          {1, 0, 0, 0},
                          {0, 1, 0, 0},
                          {0, 0, 1, 1}};
   // vreau sa aplic formula S_n = S_4 * C^(n-4)
 
   // C = C^(n-4)
   power_matrix(C, n - 4, C);
 
   // sol = S_4 * C = dp[n] (se afla pe ultima pozitie din S_n,
   // deci voi folosi ultima coloana din C)
   int sol = 1 * C[0][3] + 1 * C[1][3] + 1 * C[2][3] + 2 * C[3][3];
   return sol % MOD;
}

In arhiva demo-lab04.zip gasiti o sursa completa in care se realizeaza:

• o verificare a faptului ca cele 2 implementari ( gardurile_lui_Gigel si garduri_rapide) produc aceleasi rezultate
• un benchmark in care cele 2 implementari sunt comparate
  • pe sistem uzual (laptop) s-au obtinut urmatoarele rezulate:

test case: varianta simpla
n = 100000000 sol = 119; time = 0.984545 s
test case: varianta rapida
n = 100000000 sol = 119; time = 0.000021 s
 
 
test case: varianta simpla
n = 1000000000 sol = 812; time = 9.662377 s
test case: varianta rapida
n = 1000000000 sol = 812; time = 0.000022 s

• se observa clar diferenta intre cele 2 solutii (am confirmat ceea ce spunea si teoria: $O(n) $ vs $O(log(n))$); aceasta tehnica imbunatateste drastic o solutie gasita relativ usor.

$n = 3$

Raspuns: $7$

Explicatie: Toate subsirurile posibile sunt

• $[2]$
• $[2, 6]$
• $[2, 6, 4]$
• $[2, 4]$
• $[6]$
• $[6, 4]$
• $[4]$

Toate subsirurile de mai sus au suma para.

$n = 3$

Raspuns: $3$

Explicatie: Toate subsirurile posibile sunt

• $[2]$
• $[2, 1]$
• $[2, 1, 3]$
• $[2, 3]$
• $[1]$
• $[1, 3]$
• $[3]$

Subsirurile cu suma para sunt: $[2]$, $[2, 1, 3]$, $[1, 3]$.

$n = 3$

Raspuns: $3$

Explicatie: Toate subsirurile posibile sunt

• $[3]$
• $[3, 2]$
• $[3, 2, 1]$
• $[3, 1]$
• $[2]$
• $[2, 1]$
• $[1]$

Subsirurile cu suma para sunt: $[3, 2, 1]$, $[3, 1]$, $[2]$.

Cate subsiruri au suma impara?

$n = 5$ si $expr = ['T', '&', 'F', '^', 'T']$ (expr = [ true and false xor true])

Raspuns: $2$

Explicatie: Exista 2 moduri corecte de a paranteza expresia astfel incat sa obtinem rezultatul true (1).

• $ T&(F^T) $
• $ (T&F)^T $

  Complexitate temporală dorita este $O(n ^ 3)$.
  
  Optional, se pot defini functii ajutatoare precum **is_operand**, **is_operator**, **evaluate**.

Recomandam sa NU fie una din cele 3 probleme de la Test PA 2017. Recomandam sa le incercati dupa ce recapitulati acasa DP1 si DP2, pentru a verifica daca cunostintele acumulate sunt la nivelul asteptat.

Rezolvati problema extraterestrii de la Test PA 2017.

Rezolvati problema Secvente de la Test PA 2017.

Rezolvati problema PA Country de la Test PA 2017.

Rezolvati pe infoarena problema iepuri.

Hint: Exponentiere logaritmica pe matrice

Solutie:

• $dp[0] = X; dp[1] = Y; dp[0] = Z; $
• $dp[i] = (A * dp[i-1] + B * dp[i-2] + C * dp[i-3]) \ \% \ 666013$

Pentru punctaj maxim, pentru fiecare test se foloseste ecuatia matriceala atasata. Complexitate: $O(T * log(n))$.

Rezolvati pe leetcode problema Minimum Path Sum.

Rezolvati pe infoarena problema Lacusta.

Rezolvati pe infoarena problema Suma4.

Rezolvati pe infoarena problema subsir.

Rezolvati pe infoarena problema 2sah.

Hint: Exponentiere logaritmica pe matrice

O descrie detaliata se afla in arhiva OJI 2015.

Articolul de pe leetcode conține o listă cu diverse tipuri de probleme de programare dinamică, din toate categoriile discutate la PA.