Laborator 4: Programare Dinamică (continuare)

$n = 3$

Raspuns: 64 (inmultiri scalare)

Explicatie: Avem 3 matrici:

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

In functie de ordinea efectuarii inmultirilor matriciale , numarul total de inmultiri scalare poate sa fie foarte diferit:

• $(AB)C$ ⇒ $24 + 40 = 64$ de inmultiri
  • explicatie: $X = (AB)$ genereaza $2 * 3 * 4 = 24$ inmultiri, $(XC)$ genereaza $2 * 4 * 5 = 40$ de inmultiri
• $A(BC)$ ⇒ $60 + 30 = 90$ de inmultiri
  • explicatie: $X =(BC)$ genereaza $3 * 4 * 5 = 60$ inmultiri, $(AX)$ genereaza $2 * 3 * 5 = 30$ de inmultiri

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

$n = 4$

Raspuns: 48 (inmultiri scalare)

Explicatie: Avem 4 matrici:

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

In functie de ordinea efectuarii inmultirilor matriciale, numarul total de inmultiri scalare poate sa fie foarte diferit:

• $(AB)C)D$ ⇒ $24 + 16 + 12 = 52$ inmultiri
  • explicatie: $X = (AB)$ genereaza $2 * 3 *4 = 24$ inmultiri scalare, $Y = (XC)$ genereaza $2 * 4 * 2 = 16$ inmultiri scalare, $Z = YD$ genereaza $2 * 2 *3 = 12$ inmultiri scalare
• $(A(BC))D$ ⇒ $24 + 12 + 12 = 48$ inmultiri
  • explicatie: $X = (BC)$ genereaza $3 * 4 * 2 = 24$ inmultiri scalare, $Y = (AX)$ genereaza $ 2 * 3 * 2= 12$ inmultiri scalare, $Z = YD$ genereaza 2 * 2 * 3$ = 12$ inmultiri scalare
• $(AB)(CD)$ ⇒ $ = $ inmultiri
  • explicatie: $X = (AB)$ genereaza $2 * 3 * 4 = 24$ inmultiri scalare, $Y = (CD)$ genereaza $4 * 2 * 3 = $24 inmultiri scalare, $Z = XY$ genereaza $2 * 4 * 3 = 24$ inmultiri scalare
• $A((BC)D)$ ⇒ $24 + 18 + 27 = 69$ inmultiri
  • explicatie: $X = (BC)$ genereaza $3 * 4 * 2 = 24$ inmultiri scalare, $Y = (XD)$ genereaza $3 * 2 * 3 = 18$ inmultiri scalare, $Z = AY$ genereaza $3 * 3 * 3 = 27$ inmultiri scalare
• $A(B(CD))$ ⇒ $24 + 36 + 18 = 78$ inmultiri
  • explicatie: $X = (CD)$ genereaza $4 * 2 * 3 = 24$ inmultiri scalare, $Y = (BX)$ genereaza $3 * 4 * 3 =36 $ inmultiri scalare, $Z = AY$ genereaza $2 * 3 * 3 = 18$ inmultiri scalare

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

$n = 4$

Raspuns: 2856 (inmultiri scalare)

Explicatie: Avem 4 matrici:

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

In functie de ordinea efectuarii inmultirilor matriciale , numarul total de inmultiri scalare poate sa fie foarte diferit:

• $((AB)C)D$ ⇒ 10582 inmultiri
• $(AB)(CD)$ ⇒ 54201 inmultiri
• $(A(BC))D$ ⇒ 2856 inmultiri
• $A((BC)D)$ ⇒ 4055 inmultiri
• …

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

// kInf este valoarea maxima - "infinitul" nostru
const unsigned long long kInf = std::numeric_limits<unsigned long long>::max();
 
// T = O(n ^ 3) 
// S = O(n ^ 2) - stocam n x n intregi in tabloul dp
 unsigned long long solve_podm(int n, const vector<int> &d) {
    // dp[i][j] = numarul MINIM inmultiri scalare cu codare poate fi calculat produsul
    //            matricial M_i * M_i+1 * ... * M_j
    vector<vector<unsigned long long>>  dp(n + 1, vector<unsigned long long> (n + 1, kInf));
 
    // Cazul de baza 1: nu am ce inmulti 
    for (int i = 1; i <= n; ++i) {
        dp[i][i] = 0ULL;  // 0 pe unsigned long long (voi folosi mai incolo si 1ULL)
    }
 
    // Cazul de baza 2: matrice d[i - 1] x d[i] inmultia cu matrice d[i] x d[i + 1] 
    // (matrice pe pozitii 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) {            // fixam lungimea intervalului (2, 3, 4, ...)
        for (int i = 1; i + len - 1 <= n; ++i) {    // fixam capatul din stanga: i
            int j = i + len - 1;                    // capatul din dreapta se deduce: j
 
            // Iteram prin indicii dintre capete, spargand sirul de inmultiri in doua (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];
 
                // actualizam solutia daca este mai buna
                dp[i][j] = min(dp[i][j], new_sol); 
            }
        }
    }
 
    // Rezultatul se afla in dp[1][n]: Numarul MINIM de inmultiri scalare
    // pe care trebuie sa le facem pentru a obtine produsul M_1 * ... * M_n
    return dp[1][n];
 
}

• definitie : b este inversul modular a lui a in raport cu MOD daca $ a * b = 1 (modulo \ MOD)$
• utilizare : $ \frac{a}{b} \ \% \ MOD = ((a \ \% \ MOD) * (invers(b) \ \% \ MOD)) \ \% \ MOD $
• calculare : deoarece la PA aceasta discutie are sens doar in contextul posibilitatii implementarii unei recurente DP in care folosim resturile doar pentru a evita overflow/imposibilitatea de a retine rezultatul pe tipurile standard de tip int (adica nu ne intereseaza sa dam o metoda generala pentru invers modular), vom simplifica problema - MOD este prim!!!
• Mica teorema a lui Fermat: Daca p este un numar prim si a este un număr intreg care nu este multiplu al lui p, atunci $a^{p-1} = 1 (modulo \ p)$.
  • din definitia inversului modular, reiese ca a si b nu sunt multipli ai lui MOD
  • introducand notatiile noastre in teorema si prelucrand obtinem
    • $a ^ {MOD - 1} = 1 (modulo \ MOD) <=> a * (a ^{MOD-2}) = 1 (modulo \ MOD)$
    • deci inversul modular al lui a (in aceste conditii specifice) este $b = a ^ {MOD -2 }$

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

Raspuns: 1 (un singur gard)

Explicatie: Se poate forma un singur gard in fiecare caz, dupa cum este ilustrat si in figura Garduri_123.

$n = 4$

Raspuns: 2

Explicatie: Se pot forma 2 garduri, in functie de cum asezam piesele, dupa cum este ilustrat si in figura Garduri_4. Observam ca de fiecare daca cand punem o piesa in pozitie orizontala, de fapt suntem obligati sa punem 4 piese, una peste alta!

$n = 5$

Raspuns: 3

Explicatie: Se pot forma 3 garduri, in functie de cum asezam piesele, dupa cum este ilustrat si in figura Garduri_5.

• daca dorim ca acest gard sa se termine cu 4 piese in pozitie orizontala (una peste alta - marcat cu rosu), atunci la stanga mai ramane de completat un subgard de lungime 1, in toate modurile posibile
• daca dorim ca acest gard sa se termine cu o piesa in pozitie verticala (marcat cu rosu), atunci la stanga mai ramane de completat un subgard de lungime 4, in toate modurile posibile

$n = 6$

Raspuns: 4

Explicatie: Se pot forma 4 garduri, in functie de cum asezam piesele, dupa cum este ilustrat si in figura Garduri_6.

• daca dorim ca acest gard sa se termine cu o piesa in pozitie verticala (marcat cu rosu), atunci la stanga mai ramane de completat un subgard de lungime 5, in toate modurile posibile
• daca dorim ca acest gard sa se termine cu 4 piese in pozitie orizontala (una peste alta - marcat cu rosu), atunci la stanga mai ramane de completat un subgard de lungime 2, in toate modurile posibile

#define MOD 1009
int gardurile_lui_Gigel(int n) {
    // cazurile de baza
    if (n <= 3) return 1;
    if (n == 4) return 2;
 
    vector<int> dp(n + 1); // pastrez indexarea de la 1 ca in explicatii
 
    // cazurile de baza
    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];
}

Mentionez ca am folosit expresia $dp[i] = (dp[i - 1] + dp[i - 4]) \ \% \ MOD$ in loc de $dp[i] = ((dp[i - 1] \ \% \ MOD) + (dp[i - 4] \ \% \ MOD)) \ \% \ MOD$, deoarece [e valorile anterior calcule in dp a fost deja aplicata operatia $%$.

Am plecat cu numerele $1, 1, 1, 2$ si la fiecare pas rezultatul stocat este $\ \% \ MOD$, deci tot ce este stocat deja in dp este un rest in raport cu MOD. NU mai era nevoie deci sa aplica % si pe termenii din paranteza.

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 subsirule 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 subsirule 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 subsirule 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]$.

Dar cate subsiruri au suma impara?

Problema este preluata de aici. Solutia se gaseste aici.

$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.