Aprendre a programar 10 - Algorismes fonamentals

Aquesta és la part de PRO1 a la que la gent fa menys cas. De fet, quan ho vaig fer, vaig passar olímpicament d'ella. No obstant, últimament els ha donat per posar alguna pregunta relacionada a l'últim examen, així que potser millor saber-ne una mica. Tampoc és que siguin gaires, i uns quants ja estan explicats a l'article anterior.

Producte ràpid:

Imaginem que tenim un sistema que només pot fer sumes, restes, productes per 2, i divisions per 2. Com podem implementar un algorisme que multipliqui números de manera més o menys ràpida?

Primera idea:

int prod (int a, int b) {
    int r=0;
    for (int i=0;i<a;++i) r+=b;
    return r;
}


Si mirem el seu cost, farà a iteracions. És a dir, temps Θ(a). Lineal respecte a.
Què passaria si la cridem pel cas a=2^31-1, b=1? Doncs que ho podríem haver fet en una iteració, però ho farem en moltes. Podríem pensar de millorar-ho fent



int prod (int a, int b) {
    int r=0;
    if (a<b) for (int i=0;i<a;++i) r+=b;
    else for (int i=0;i<b;++i) r+=a;
    return r;
}


No obstant, això segueix tenint temps lineal. Podem optimitzar-ho?

Versió logarítmica

int prod (int a, int b) {
    if (b==0) return 0;
    if (b%2==0) return prod (a*2,b/2);
    return prod(a,b-1)+a;
}


Anem a mirar això. Si b==0, a*b=0. Si b és parell, a*b=(a*2)*(b/2). I això el que fa és que, sempre que b sigui parell, estarem dividint entre dos.

Anem a analitzar el cost. Quin és el cas millor? El cas millor és quan b sigui 2^n. En aquest cas, necessitarem n+1 crides. És a dir, temps logarítmic. I quin és el cas pitjor? Quan sigui (2^n)-1. Aquí, la meitat de crides faria la resta, i l'altra meitat la divisió. No obstant, això segueix sent temps logarítmic, ja que, malgrat que fem el doble de crides, Θ(2*log b)=Θ(log b). 

Se'n pot fer una versió iterativa, la teniu a l'arxiu d'algorismes fonamentals de PRO1

Exponenciació ràpida:

Aquesta idea la podem extendre a la potència. 

int exp(int a, int n) {

    int r=1;

    for (int i=0;i<n;++i) r*=a;

    return r;

}

Aquest és l'algorisme que faríem intuitivament. No obstant, això té cost Θ(n). En canvi, sabem que a^n=(a^(n/2))^2. Per tant, aprofitant això, podem fer

int exp(int a, int n) {

    if (n==0) return 1;

    int aux=exp(a,n/2);

    if (n%2==0) return aux*aux;

    return aux*aux*a;

}

Això té temps Θ(log n). 

Cerca dicotòmica

Imaginem que ens demanen fer una funció

int cerca(const vector<int>& v, int l, int r, int x);

Que rep un vector ordenat creixentment, i mira si entre l i r hi ha l'element x. Si no hi és, retorna -1. Com ho podem fer?

int cerca(const vector<int>& v, int l, int r, int x) {

    for (int i=l;i<=r;++i) {

        if (v[i]==x) return i;

    }

    return -1;

}

Ara, això ens obliga a recórrer tot el subvector. Però sabem que el vector està ordenat, per tant arribarà un moment on les iteracions seran inútils, ja que l'element no pot estar allà. Podem fer

int cerca(const vector<int>& v, int l, int r, int x) {

    for (int i=l;i<=r;++i) {

        if (v[i]==x) return i;

        if (v[i]>x) return -1;

    }

    return -1;

}

Això el que fa és que, a la que ens hem passat la posició per on hauria d'estar, retornem -1, perquè no hi pot estar. No obstant, això segueix sent lineal, no? No podem fer per dividir els casos d'alguna manera?

Podem aprofitar-nos d'una cosa. Si volem mirar entre l i r, podem calcular el punt mig ((l+r)/2), i mirar si l'element està allà. Si hi és, ja hem acabat. Si no, mirem si és menor o major. Si l'element que hi ha a aquesta posició és menor que el que busquem, vol dir que està a la dreta. Si no, a l'esquerra. Una obvietat, no? Però, si ens hi fixem, hem dividit per dos la mida del subvector on volem buscar.

És a dir, per trobar x en un vector de mida n, tindrem:

x està en vector de mida n?
x està en vector de mida n/2?
x està en vector de mida n/4?
x està en vector de mida n/8?
...
x està en vector de mida 1?

És a dir, si no hi és, farem log2(n)+2 accessos al vector. En canvi, amb una cerca lineal, faríem n accessos. Per un vector de mida 100.000, amb la cerca dicotòmica necessitarem, com a molt, 19 accessos amb la dicotòmica. Si ho fem amb la lineal, podem veure quants en farem. No cal parlar ja de si hi ha 1.000.000, o 1.000.000.000 d'elements (42 accessos amb la dicotòmica). La mida màxima d'un enter és 2.147.483.647, amb 33 accessos podem comprovar si un element hi és

Podem fer el codi. Queda una cosa similar

int cerca(const vector<int>& v, int l, int r, int x) {

    if (l>r) return -1;

    int mig=(l+r)/2;

    if (v[mig]==x) return mig;

    if (v[mig]<x) return cerca(v,mig+1,r,x);

    return cerca(v,l,mig-1,x);

}

De fet, podem analitzar la complexitat d'això. El cas pitjor serà quan no hi sigui, que el que farà és anar dividint, fins que arribi un punt on l>r, i retorni -1. És a dir, la distància entre l i r s'anirà dividint entre 2, fins que sigui 0, una crida més, i ja estarà. Això és temps Θ(log n)

Suma de vectors dispersos

Quan tenim un vector on molts dels seus elements són 0, utilitzar un vector vol dir malgastar molt espai. Una cosa que podem fer és crear un vector dispers. En aquest vector, tots els elements són una parella de dos elements. El primer membre és la posició, i el segon el valor. Així, el primer vector es convertiria en el segon

0 0 0 0 0 30 0 0 0 0

{5,30}

Com implementem un algorisme que sumi dos d'aquests vectors?

vector<pair<int,int> > suma_dispersa (const vector<pair<int,int> >& v1, const vector<pair<int,int> >&v2) {

    int n, m;

    n=v1.size();

    m=v2.size();

    int i=0,j=0;

    vector<pair<int,int> > r;

    while (i<n and j<m) {

        if (v1[i].first<v2[j].second) {

            r.push_back(v1[i++]);

        }

        else if (v1[i].first>v2[j].second) {

            r.push_back(v2[j++]);

        }

        else {

            r.push_back(pair<int,int>(v1[i].first,v1[i].second+v2[i].second));

            ++i;

            ++j;

        }

    }

    return r;

}

Suposo que s'entén el que faig, no?

Producte de matrius

typedef vector<int> Fila;

typedef vector<Fila> Matriu;

void producte(const Matriu& A, const Matriu& B, Matriu& C) {

    C=Matriu(A.size(),Fila(B.size(),0));

    for (int i=0;i<int(A.size());++i) {

        for (int j=0;j<int(B[0].size());++j) {

            for (int k=0;k<int(B.size());++k) {

                C[i][j]+=A[i][k]*B[k][j];

            }

        }

    }

}

Aquest algorisme seria el que se'ns vindria al cap més intuitivament. Bàsicament és programar el mateix que fem per multiplicar dues matrius a mà. No obstant, no és la millor versió per multiplicar matrius. Si volem millorar-ho una mica sense gaire esforç, podem fer-ho així:

void producte(const Matriu& A, const Matriu& B, Matriu& C) {

    C=Matriu(A.size(),Fila(B.size(),0));

    for (int k=0;k<int(B.size());++k) {

        for (int i=0;i<int(A.size());++i) {

            for (int j=0;j<int(B[0].size());++j) {

                C[i][j]+=A[i][k]*B[k][j];

            }

        }

    }

}

Bàsicament perquè, per temes de memòria cache, és molt més eficient recórrer una matriu per files que per columnes. I com podem veure, a la primera versió, per cada execució del bucle intern, estem recorrent A per files, B per columnes, i C no canvia. En canvi, amb la segona versió, recorem C i B per files, i A no canvia. 

No obstant, les dues versions triguen temps cúbic. Si cal una multiplicació més ràpida, caldria utilitzar l'algorisme de Strassen

Trobar l'arrel d'una funció continua 

Suposo que tothom coneix el teorema de Bolzano. És aquell que diu

Sigui f:[a,b] -> Reals tal que f(a)*f(b)<0, existeix c a l'interval (a,b) tal que f(c)=0

Cosa que, expressada en un llenguatge més col·loquial, ve a ser que si una funció continua en un interval és positiva en un punt de l'interval, i negativa a un altre, en algun moment val 0. Una obvietat, potser, però una obvietat important. Perquè permet trobar una solució d'una funció amb molta precisió. 

double arrel(double a, double b, double epsilon) {

    if (b-a<=epsilon) return a;

    double c=(a+b)/2;

    if (f(a)*f(c)<=0) return arrel(a,c,epsilon);

    return arrel(c,b,epsilon);

}

És a dir, si la diferència entre els dos és menor que epsilon, ja hem trobat una solució. Si no, mirem què passa al mig, i en funció d'això, anem fent.