momente şi schiţe de informatică şi matematică
To attain knowledge, write. To attain wisdom, rewrite.

În jurul unei probleme de echilibru (I)

C++ | STL | combinatorică | partiţii | vector
2011 may

Vom nota prin [n] mulţimea {1, 2, ..., n} (mai zicem "segmentul" natural, de la 1 la n). Pentru mulţimi finite de numere naturale X si Y, notăm prin s(X) suma numerelor din mulţimea X şi dacă s(X) ≥ s(Y), notăm d(X, Y) = s(X) - s(Y).

Vom viza descompunerile mulţimii [n] în două submulţimi A şi B = [n] - A, pentru care valoarea d(A, B) este cea mai mică posibilă. Pentru o referire comodă, să numim s-partiţie (citim "espartiţie") o astfel de descompunere.

Putem contura ideile de lucru punând întâi o întrebare mai generală: există o partiţie [n] = (A, B) încât d(A, B) să aibă o valoare fixată în prealabil? Exemplificarea următoare sugerează în bună măsură dezvoltările care ar fi de făcut:

{1, 2, 3, 4, 5} nu admite nicio descompunere (A, B) încât d(A, B) să fie 2 - fiindcă avem s([5]) = 15 = s(A) + s(B) şi atunci s(A) - s(B) = 2 ar duce la 2*s(A) = 17 (absurd); în schimb, cu d(A, B) = 3 găsim 2*s(A) = 18 - iar aceasta pare suficient pentru a şi evidenţia o descompunere: scădem cel mai mare element, s(A) - 5 = 9 - 5 = 4 şi deja găsim A = {5, 4} şi B = [5] - A = {1, 2, 3}.

Exerciţii
    1. Motivaţi: descompunerea indicată mai sus (pentru [5] şi d(A, B) = 3) este unică.
    2. [40] admite o partiţie (A, B) încât s(A) = s(B) (adică având d(A, B) = 0)?
    3. Găsiţi o partiţie (A, B) pentru [10], cu d(A, B) = 5. Câte astfel de partiţii există?

Arătăm întâi că minimul d(A, B) este 0 sau 1 - în funcţie de valoarea n % 4 - şi deducem (inclusiv pe seama unui mic algoritm de construcţie)există s-partiţii, pentru orice n.

Existenţa şi construcţia unei s-partiţii

Avem s([n]) = 1 + 2 + ... + n = n*(n+1)/2 = T(n) - numărul triunghiular de rangul n (alăturat sugerăm geometria implicată: cele două triunghiuri de câte 1+2+...+n steluţe constituie împreună un dreptunghi de n*(n+1) elemente).

Avem s(A) + s(B) = s([n]) = T(n) şi pe de altă parte, s-a cerut ca s(A) - s(B) = d(A, B).
Însumând, rezultă că 2*s(A) = T(n) + d(A, B) (deci T(n) şi d(A, B) trebuie să aibă aceeaşi paritate).

Pentru ca d(A, B) să fie minim, deducem că ar trebui să luăm d(A, B) = 0 dacă T(n) este par şi respectiv, d(A, B) = 1 dacă T(n) este impar; dar - de aceea am şi folosit "ar trebui" - va fi totuşi de dovedit că există A (altfel ar fi de analizat prima valoare d(A, B) mai mare ca 1 - deci 2, 4, etc. dacă T(n) este par şi respectiv 3, 5, etc. dacă T(n) este impar - pentru care A poate fi construită).

T(n) = n*(n+1)/2 este par numai dacă n sau (n+1) este multiplu de 4. Deci trebuie să luăm d(A, B) = 0 pentru n % 4 = 0 sau n % 4 = 3 şi respectiv d(A, B) = 1 pentru n % 4 = 1 sau n % 4 = 2.

În oricare caz, următorul algoritm construieşte mulţimea căutată A: scădem din s([n]) şi apoi din suma rămasă, numerele n, n-1, n-2, etc. până când rezultă o valoare mai mică decât scăzătorul curent; mulţimea A este formată din scăzătorii respectivi, împreună cu valoarea ultimei scăderi.

Exemplificare, pentru n = 14:

T(14) = 105, 14 % 4 = 2 ==> s(14) = 53.
53 - 14 = 39 > 14; A = {14}
          39 - 13 = 26 > 13; A = {14, 13}
                    26 - 12 = 14 > 12; A = {14, 13, 12} 
                              14 - 11 = 3 < 11; A = {14, 13, 12, 11},
A = {14, 13, 12, 11, 3}.

Exerciţii
    4. Derivaţi din mulţimea A rezultată mai sus, alte s-partiţii ale segmentului natural [14].
    5. Mulţimea A calculată mai sus are 5 elemente; există s-partiţii pentru n = 14, în care A să aibă mai puţin de 5 elemente?
    6. Implementaţi (C++, javaScript, sau alt limbaj) algoritmul descris mai sus.

Observaţii

1) Mulţimea A (şi deci o s-partiţie pentru [n]) este unic determinată prin algoritmul de mai sus şi intuim că ea poate servi ca "bază" pentru obţinerea tuturor celorlalte s-partiţii (de exemplu, înlocuind mai sus 14 = 10 + 4 obţinem o altă s-partiţie pentru [14]).
Fiindcă această "bază" s-a obţinut colectând succesiv cel mai mare dintre numerele posibile în condiţiile specificate de algoritm, rezultă pe de o parte că "baza" are cardinalul minim posibil şi pe de altă parte, că "baza" este cea mai mare într-o anumită ordine lexicografică (văzând mulţimile A din s-partiţii ca mulţimi ordonate descrescător şi ajustând la comparaţii numerice "ordinea lexicografică" obişnuită).
Dar nu prea ajunge să intuim… şi încă nu putem specula coerent aceste observaţii.

2) Punctul de plecare a acestui articol despre s-partiţii este o problemă pe care am "tratat-o" anterior (doar că în manieră infernal-didactică) în Exerciţii de programare în C, C++: adună sau scade între ele numerele 1, 2, ..., n încât să rezulte 1. Din cele de mai sus, rezultă că dacă n % 4 este 1 sau 2, iar (A, B) este o s-partiţie a mulţimii [n] atunci avem într-adevăr s(A) - s(B) = 1, cum cerea problema menţionată (în secvenţa 1, 2, ..., n punem plus numerelor conţinute de A şi minus celor din B).

Este adevărat că puteam argumenta simplu (fără nici un calcul, sau construcţie)există s-partiţii: şirul tuturor diferenţelor s(X) - s(Y), când X parcurge mulţimea părţilor lui [n], iar Y este complementara lui X este finit - prin urmare are un cel mai mic termen d(A, B). Dar (spre deosebire de calculul şi algoritmul descris mai sus) acest raţionament nu poate evidenţia şi care este acest d(A, B).

Numărul de s-partiţii

Mai sus am găsit 53 = 14 + 13 + 12 + 11 +3 ceea ce este o partiţie întreagă a numărului 53; la fel, oricare mulţime A dintr-o s-partiţie a mulţimii [14], reprezintă o "partiţie întreagă" în părţi distincte de valoare maximă 14, a lui 53 (dar şi 53 = 25 + 25 + 1 + 1 + 1 este o "partiţie întreagă" a lui 53).

Din cele stabilite mai sus rezultă că s-partiţiile segmentului [n] reprezintă exact partiţiile întregi cu părţi distincte de valoare cel mult n ale numărului T(n)/2 (dacă n%4 este 0 sau 3), sau ale numărului (T(n) + 1)/2 (dacă n%4 este 1 sau 2). Aceasta ar permite calcularea numărului de s-partiţii exploatând formule sau tabele pentru partiţiile întregi (a vedea A000041).

Dar aici vom proceda direct, ilustrând o metodă specifică analizei combinatorice: pentru a număra anumite configuraţii se caută o funcţie generatoare pentru numerele respective, adică o functie care conduce la o serie de puteri ∑ak xk ai cărei coeficienţi ak sunt numerele căutate (de exemplu, (1+x)n este o funcţie generatoare pentru coeficienţii binomiali: în dezvoltarea corespunzătoare, ak este numărul combinărilor de n elemente luate câte k).

Produsul (1 + x)(1 + x2)(1 + x3)...(1 + xn) este un polinom de gradul 1 + 2 + 3 + ... + n = T(n), fiindcă monomul de grad maxim al dezvoltării rezultă din x1x2x3...xn.

Monomul care conţine x6 va rezulta însumând monoamele asemenea x6 (obţinut înmulţind x6 din a 6-a paranteză cu termenii 1 din celelalte), x5x1 (înmulţind x5 din a 5-a paranteză, x1 din prima şi termenii 1 din celelalte), x4x2 şi în sfârşit, x3x2x1 (desigur, am presupus că n este suficient de mare; de exemplu dacă am înmulţi numai primele 4 paranteze - adică am lua n = 4 - atunci coeficientul lui x6 ar fi doar 2).

Deci în final, x6 va avea coeficientul 4; dar să vedem şi faptul că monoamele asemenea intermediare "reprezintă" exact partiţiile întregi stricte (cu termeni distincţi) ale lui 6 (anume 6, 5 + 1, 4 + 2 şi 3 + 2 + 1), în număr de 4.

În general, coeficientul lui xk din produsul considerat mai sus ne va da numărul de partiţii întregi stricte ale numărului k. Prin urmare - având în vedere legătura evidenţiată mai sus între numărul de s-partiţii şi numărul de partiţii întregi stricte - obţinem:

Teoremă. Numărul de s-partiţii ale segmentului natural [n] (cu n ≥ 3) este dat de coeficientul lui

xT(n)/2 dacă n % 4 ∈ {0, 3}, sau

x(T(n) + 1)/2 dacă n % 4 ∈ {1, 2}

din dezvoltarea polinomială a produsului (1 + x)(1 + x2)(1 + x3)...(1 + xn).

Să notăm însă că produsul indicat nu este "pe deplin" o funcţie generatoare pentru valorile dorite aici, fiindcă nu toţi coeficienţii dezvoltării sunt relevanţi - ci numai cel de rang T(n)/2 sau (T(n)+1)/2.

Tabelul cardinalelor s-partiţiilor (C++)

Cu alte cuvinte, dacă am vrea să constituim un tabel al numărului de s-partiţii pentru n = 3..100 va trebui să dezvoltăm P1 = (1+x)(1+x2)(1+x3) şi să reţinem în tabel coeficientul lui x3 (fiindcă T(3)/2 = 3), apoi să calculăm P2 = P1*(1+x4) şi să înscriem în tabel coeficientul lui x5 (T(4)/2 = 5), apoi P3 = P2*(1+x5) înscriind coeficientul lui x8 (fiindcă (T(5)+1)/2 = 8), etc.

Angajând C++, alegem să modelăm polinoamele ca obiecte de tip vector (#include <vector>); cu "vectori" se lucrează la fel ca şi cu tablouri C obişnuite, însă cu avantajul că nu mai trebuie să ne preocupăm de alocarea memoriei şi nici să menţinem separat "dimensiunea curentă". De exemplu, vector<int> V(5, -1); alocă un tablou pentru 5 întregi cu valoarea iniţială -1; apoi, V[1]=5 şi V.push_back(22) transformă tabloul în (-1,5,-1,-1,-1,22); dimensiunea curentă este dată de V.size().

Polinoamele pot să fie implicate în fel de fel de probleme, încât este firesc să constituim un modul separat în care să grupăm operaţiile cu polinoame. În cazul nostru avem nevoie numai de operaţia de înmulţire şi putem considera următorul fişier "header":

/* polinom.h  introduce Polinom şi operaţii cu polinoame */
#include <vector>
typedef unsigned long long coeficient; // 0..2^64 - 1 
typedef std::vector<coeficient> Polinom;
Polinom produs(Polinom&, Polinom&);

Coeficienţii care ar fi de obţinut aici sunt mari; GCC 4.5 (pe Linux) prevede şi tipul long long, extinzând reprezentarea întregilor la 8 octeţi (64 de biţi) - astfel că (folosind modificatorul unsigned) vom putea opera cu valori întregi de până la 264 - 1 = 18,446,744,073,709,551,615.

/* polinom.cpp  Presupune că polinoamele sunt astfel încât:
rangul în vector = gradul monomului asociat */
#include "polinom.h"

Polinom produs(Polinom& a, Polinom& b) {
    int a_sz = a.size(), // = grad + 1
        b_sz = b.size();
    int c_sz = a_sz + b_sz - 1;
    Polinom c(c_sz, 0); // produsul iniţial (0,0,...,0)
    for(int i = 0; i < a_sz; i++)
        for(int j = 0; j < b_sz; j++)
            c[i + j] += a[i] * b[j];
    return c;
}

În programul "principal" trebuie să includem "headerul" polinom.h (şi nu "polinom.cpp"):

// count_s-partitions.cpp
/* Numărul de s-partiţii ale mulţimii [1..n] este egal cu
coeficientul de rang T(n)/2 sau cel de rang (T(n) + 1)/2 
din dezvoltarea (1+x)(1+x^2)(1+x^3)...(1+x^n) */
#include "polinom.h"
#include <iostream>
using namespace std;

int main() {
    Polinom p(2, 1); // polinomul iniţial 1 + x
    for(int n = 2; n <= 100; n++) { // coeficienţi > 2^64 dacă n > 72
        Polinom next(n + 1, 0);
        next[0] = next[n] = 1; // (1,0,0,...,0,0,1)
        p = produs(p, next);
        int Tn = n*(n+1) / 2;  
        switch(n % 4) {
            case 1: case 2: Tn += 1; // fără 'break'!
            case 0: case 3: cout << n << ": " << p[Tn / 2] << '\n'; 
        }
    }
}

Pentru a obţine programul executabil, invocăm g++ de pe linia de comandă:

vb@vb:~/docere/doc/Prolog$ g++ -o cnt-parts  polinom.cpp  count_s-partitions.cpp

şi apoi, pentru a obţine rezultatele într-un fişier (şi nu pe ecran, cum s-a prevăzut iniţial):

vb@vb:~/docere/doc/Prolog$ ./cnt-parts  >  coef-part.txt

Consultând "coef-part.txt" constatăm privitor la ordinul de mărime al coeficienţilor calculaţi:

...
68: 714733339229024336    coloana ultimei cifre: 23
69: 1398852939566205174
70: 2738645100853765060
71: 5363288299585278800
72: 10506331021814142340  coloana ultimei cifre: 25
73: 2140221165180875491   coloana ultimei cifre: 24

Evident, coeficienţii vizaţi cresc odată cu n; deci pe fiecare rând, coloana ultimei cifre ar trebui să fie cel puţin egală cu coloana ultimei cifre de pe rândul precedent. Constatăm că aceasta este valabil până la linia n = 72 (valoarea înscrisă la n = 73 este mai mică decât cea pentru n = 72, ceea ce este absurd).

Deducem că valoarea corectă pentru n = 73 deja depăşeşte domeniul valorilor "long long"; deci valoarea corectă pentru n = 73 ar trebui să fie cu 264 mai mare decât cea reţinută (fiindcă în cazul depăşirii unui domeniu standard de valori întregi se "reţine" restul faţă de limita maximă):

73:      2,140,221,165,180,875,491  rezultatul înscris (trunchiat) de program
2^64:   18,446,744,073,709,551,616
suma:   20,586,965,238,890,427,107  valoarea corectă pentru n = 73

Deci programul de mai sus poate furniza rezultatele dorite doar până la n = 72 (valorile afişate începând cu n = 73 sunt corecte doar modulo 264). Însă putem modifica programul redat încât să obţinem un tabel mai amplu, de exemplu până la n = 500 - anume astfel: includem gmp.h (implicând biblioteca GMP) şi pentru simplitate, renunţăm să folosim clasa "vector" (revenind la tablourile C clasice - acum tablouri de valori de tip mpz_t, introdus de GMP pentru reprezentarea întregilor "cu oricâte cifre"; am dat un exemplu de folosire a bibliotecii GMP în "Hello World!").

Verificări, reluări şi simplificări

Este firesc în general, obiceiul de a testa sau a verifica rezultatele furnizate de un program sau altul. Dar dacă programul respectiv exprimă direct şi fără "şiretlicuri" de programare un rezultat teoretic deja demonstrat, atunci chiar nu există motive să ne îndoim de valabilitatea rezultatelor obţinute!

Programul de mai sus exprimă aproape liniar teorema pe care am demonstrat-o anterior, aşa că… rezultatele furnizate (până la n = 72) sunt corecte!

Simpla verificare directă a unui rezultat este mai degrabă neinteresantă (cu excepţii de genul cazului evidenţiat mai sus pentru n = 73) şi de obicei încercăm să corelăm verificarea manuală cu aspecte de context - în cazul de faţă, cu un mecanism de generare a configuraţiilor numărate de program.

Să probăm manual că rezultatul 14 înscris de program pentru n = 8 este într-adevăr corect. Să încercăm să enumerăm "ordonat" s-partiţiile mulţimii [8], reluând algoritmul de construcţie pe care l-am evidenţiat la început.

Avem T(8)/2 = 18; selectăm cel mai mare element 8 şi repetăm pentru 18 - 8 = 10 (având 10 ≥ 7): selectăm 7 şi rămâne 3 - deci "cea mai mare" într-o anumită ordine lexicografică a s-partiţiilor este (8, 7, 3). Pentru a obţine pe următoarea în această ordine lexicografică, ar trebui să modificăm precedenta s-partiţie începând de la sfârşitul ei (menţinând pe cât se poate, valorile înscrise la începutul ei) - deci să-l descompunem pe 3, găsind astfel (8, 7, 2, 1). Mai departe, 1 şi 2 nu se pot descompune, deci revenim la valoarea 7 şi la "baza" 10 din care scăzusem iniţial pe 7; selectăm 6 în loc de 7 şi 10 - 6 = 4, obţinând "următoarea" s-partiţie (8, 6, 4). Apoi descompunem 4 găsind (8, 6, 3, 1), revenim la 6 şi înlocuim cu 5 obţinând (8, 5, 4, 1), etc.:

[8, 7, 3]             18-8=10 ≥ 7, 10-7 = 3
[8, 7, 2, 1]          descompune pe 3
[8, 6, 4]             înlocuieşte 7, de la "baza" 10: 10-6 = 4
[8, 6, 3, 1]          descompune pe 4
[8, 5, 4, 1]          înlocuieşte 6, de la baza 10: 10-5=5 ≥ 5, 5-4 = 1
[8, 5, 3, 2]          4 nu se poate descompune! Revenim la 5: 5-3 = 2
[8, 4, 3, 2, 1]       înlocuieşte 5, de la baza 10: 10-4=6 ≥ 3, 6-3 = 3

Fiecare listă este aici ordonată descrescător şi împreună, listele redate mai sus reprezintă toate părţile "A" din s-partiţiile mulţimii [8] care conţin ("încep") cu 8 (în număr de 7 s-partiţii). Întrucât pentru n = 8 (şi în general, pentru n%4 = 0 sau 3) avem s(A) = s(B), rezultă că mai departe vom obţine de fapt complementarele mulţimilor redate mai sus (altfel spus, obţinem a doua oară fiecare s-partiţie dintre cele deja reprezentate mai sus):

[7, 6, 5]             Înlocuieşte 8, de la baza 18: 18-7=11 ≥ 6, 11-6 = 5
[7, 6, 4, 1]          Complementară cu [8, 5, 3, 2] din lista de mai sus…
[7, 6, 3, 2]
[7, 5, 4, 2]
[7, 5, 3, 2, 1]
[6, 5, 4, 3]
[6, 5, 4, 2, 1]        Complementară cu [8, 7, 3]

Se validează astfel că sunt exact 14 s-partiţii pentru [8], dar am reuşit să relevăm mai mult: algoritmul de construcţie a unei s-partiţii redat la început poate fi modificat (aparent fără prea mare efort) încât să furnizeze în serie, toate s-partiţiile lui [n], în ordine lexicografică descrescătoare.

Pe de altă parte, tocmai verificarea manuală întreprinsă mai sus ne-a dezvăluit un "defect" pe care nu l-am observat pe parcursul "teoriei" elaborate: în cazul când n%4 este 0 sau 3 (şi numai în acest caz) fiecare s-partiţie (ca mulţime constituită din două părţi disjuncte) apare de două ori; în acest caz, numărul care s-a obţinut prin Teorema demonstrată mai sus este de fapt nu numărul s-partiţiilor (care acum este doar jumătatea numărului obţinut), ci "numărul de moduri" în care putem exprima numărul T(n)/2 ca sumă de termeni distincţi din [n].

Dacă vrem, putem "îndrepta" lucrurile impunând explicit din start că (A, B) ≠ (B, A), pentru oricare s-partiţie (A, B). Dar aceasta (deşi are sens, ba chiar este asumată - dar în mod tacit - în cele de mai sus) este mai degrabă o reparaţie conjuncturală: în realitate nu există 14 s-partiţii distincte pentru [8], ci numai 7; poate fi util să le enumerăm pe toate cele 14, dar de numărat sunt doar 7.

Faţă de aceste constatări, "reparaţia" justă constă probabil în rescrierea lucrurilor. Vom rescrie în întregime, păstrând însă această primă versiune aşa cum este: aici am sugerat şi am explorat treptat căile şi demersurile necesare în jurul temei propuse.

vezi Cărţile mele (de programare)

docerpro | Prev | Next