Tangens hiperboliczny

    -- Sebastian Pawlak, 2002.


Projekt dotyczy zaprogramowania współbieżnego numerycznego algorytmu wyznaczania wartości tangensa hiperbolicznego tanh(x) dla zadanej wartości argumentu x.



Aparat matematyczny

Funkcję tanh(x) da się wyrazić wzorem postaci

wzor 1

Z kolei wartości funkcji sinh oraz cosh w danym punkcie można wyznaczyć korzystając z szeregu Taylora.

Jeżeli funkcja f(x) ma n-tą pochodną f^n(x) w pewnym przedziale domkniętym zawierającym punkt a, wówczas dla każdego x z tego przedziału ma miejsce następujący wzór Taylora:

wzor 2

gdzie a < c_n < x przy x > a i x < c_n < a przy x < a zwany szeregiem Taylora.

Rozwinięcie funkcji sinh w szereg Taylora jest postaci

wzor 3

Rozwinięcie funkcji cosh w szereg Taylora jest postaci

wzor 4

Dla obu funkcji za zmienną a podstawiamy zero, gdyż punkt a=0 należy do dziedziny tych funkcji.



Metoda sekwencyjna

Algorytm sekwencyjny polega na sumowaniu n początkowych wyrazów szeregu Taylora dla funkcji sinusa hiperbolicznego oraz cosinusa hiperbolicznego, a następnie podzieleniu sumy n wyrazów sinh przez sumę n wyrazów cosh. Od liczby zsumowanych wyrazów n zależy dokładność otrzymanego wyniku.



Implementacja algorytmu sekwencyjnego

Fakt, iż kolejne wyrazy szeregów są podobne do swoich poprzedników umożliwia zredukowanie liczby operacji arytmetycznych.


Kod źródłowy pliku "sekwencyjny.c":

/* Wyznaczanie wartosci tangensa hiperbolicznego.
 * Wersja sekwencyjna.
 *
 * Sebastian Pawlak, s1222.
 */
 
#define DOKLADNOSC 100

int main(void)
{
    long double x = 1.37; /* `x', dla ktorego liczymy */
    int n;
    long double sinus = 0, cosinus = 0;
    long double silnia = 1;
    long double xx = x;

    /* Wyznaczanie wartosci sinusa hiperbolicznego */
    for(n = 0; n < DOKLADNOSC; n++) {
        sinus += xx / silnia;
        if(n < DOKLADNOSC - 1) {
            silnia *= (2 * n + 2) * (2 * n + 3);
            xx *= x * x;
        }
    }

    /* Wyznaczanie wartosci cosinusa hiperbolicznego */
    xx = 1;
    silnia = 1;
    for(n = 0; n < DOKLADNOSC; n++) {
        cosinus += xx / silnia;
        if(n < DOKLADNOSC - 1) {
            silnia *= (2 * n + 1) * (2 * n + 2);
            xx *= x * x;
        }
    }

    printf("sinus(%.3Lg) = %Lg\n", x, sinus);
    printf("cosinus(%.3Lg) = %Lg\n", x, cosinus);
    printf("tangens(%.3Lg) = %Lg\n", x, sinus / cosinus);
       
    return 0;
}



Optymalizacja algorytmu sekwencyjnego

Można zauważyć, że odpowiadające sobie wyrazy szeregu sinh oraz cosh są do siebie podobne. Wyraz szeregu sinh może być wyznaczony na podstawie odpowiadającego wyrazu szeregu cosh poprzez wymnożenie go przez x oraz podzielenie przez (2n+1). W implementacji rozwiązanie to jest nieco zmodyfikowane, choć idea pozostaje ta sama.


Kod źródłowy pliku "sekwencyjny2.c":

/* Wyznaczanie wartosci tangensa hiperbolicznego.
 * Wersja sekwencyjna.
 * Wersja zoptymalizowana.
 *
 * Sebastian Pawlak, s1222.
 */
 
#define DOKLADNOSC 100

int main(void)
{
    long double x = 1.37; /* `x', dla ktorego liczymy */
    int n;
    long double sinus = 0, cosinus = 0;
    long double silnia = 1;
    long double xx = 1;

    /* Wyznaczanie wartosci cosinusa i sinusa
     * hiperbolicznego.
     * Fakt, ze odpowiadajace sobie wyrazy obu szeregow
     * sa podobne mozna wykorzystac do zmniejszenia
     * liczby operacji arytmetycznych.
     */
    for(n = 0; n < DOKLADNOSC; n++) {
        cosinus += xx / silnia;
        sinus += xx * x / (silnia * (2 * n + 1));
        if(n < DOKLADNOSC - 1) {
            silnia *= (2 * n + 1) * (2 * n + 2);
            xx *= x * x;
        }
    }

    printf("sinus(%.3Lg) = %Lg\n", x, sinus);
    printf("cosinus(%.3Lg) = %Lg\n", x, cosinus);
    printf("tangens(%.3Lg) = %Lg\n", x, sinus / cosinus);
       
    return 0;
}



Metoda równoległa oparta na "zoptymalizowanym algorytmie sekwencyjnym"

Przedstawiony tu algorytm wyznaczania wartości tanh na n procesach polega na podzieleniu szeregów sinh i cosh na, w miarę możliwości, n równych przedziałów. Każdy z procesów wyznacza sumę wyrazów szeregu sinh oraz cosh w zadanym przedziale. Te sumy częściowe są następnie przekazywane (wraz z sumowaniem funkcją MPI_Reduce do procesu master). Komunikację między procesorami została zrealizowana za pomocą funkcji dostępnych w bibliotece MPI-CH.


Kod źródłowy pliku "rownolegly.c":

/* Wyznaczanie wartosci tangensa hiperbolicznego.
 * Wersja zrownoleglona.
 *
 * Sebastian Pawlak, s1222.
 */

#include "mpi.h"

#define DOKLADNOSC 100
#define MAX_LICZBA_PROCESOW 16

int dzielenie(int, int);
long double obliczSilnie(int);
long double obliczPotege(long double, int);

int main(int argc, char **argv)
{
    long double x = 1.37; /* `x', dla ktorego liczymy */
    int n;
    long double silnia;
    long double xx;

    int rank, size; /* numer procesu, liczba procesow */
    double czasStart;
    /* przedzialy, w ktorych beda liczyc procesy */
    struct przedzial {
        int nMin, nMax;
    } przedzialyTab[MAX_LICZBA_PROCESOW] = { 0 };
    /* wyniki czesciowe, suma wynikow czesciowych */
    struct {
        long double sinus, cosinus;
    } wynikCzesciowy = { 0 }, sumaWynikow = { 0 };

    MPI_Init(&argc, &argv);
    MPI_Comm_size(MPI_COMM_WORLD, &size);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);

    czasStart = MPI_Wtime(); /* pomiar czasu wykonywania */

    if(size > DOKLADNOSC) {
        printf("Przydzielono zbyt wiele procesow !\n");
        MPI_Abort(MPI_COMM_WORLD, -1);
    }

    /* MASTER ustala przedzialy dla wszystkich procesow */
    if(rank == 0) {
        przedzialyTab[0].nMin = 0;
        przedzialyTab[0].nMax = dzielenie(DOKLADNOSC, size) - 1;
        for(n = 1; n < size; n++) {
            przedzialyTab[n].nMin = przedzialyTab[n - 1].nMax + 1;
            przedzialyTab[n].nMax = 
              dzielenie(DOKLADNOSC - przedzialyTab[n - 1].nMax - 1,
                        size - n) + przedzialyTab[n - 1]. nMax;
        }
    }

    /* rozsylanie informacji o przedzialach do poszczegolnych
     * procesow 
     */
    MPI_Scatter(przedzialyTab, sizeof(struct przedzial), MPI_BYTE,
                przedzialyTab, sizeof(struct przedzial), MPI_BYTE,
                0, MPI_COMM_WORLD);

    /* Wyznaczanie wartosci cosinusa i sinusa hiperbolicznego.
     * Fakt, ze odpowiadajace sobie wyrazy obu szeregow sa 
     * podobne mozna wykorzystac do zmniejszenia liczby 
     * operacji arytmetycznych.
     */
    xx = obliczPotege(x, 2 * przedzialyTab[0].nMin);
    silnia = obliczSilnie(2 * przedzialyTab[0].nMin);

    for(n = przedzialyTab[0].nMin; n <= przedzialyTab[0].nMax; n++) {
        wynikCzesciowy.cosinus += xx / silnia;
        wynikCzesciowy.sinus += xx * x / (silnia * (2 * n + 1));
        if(n < przedzialyTab[0].nMax) {
            silnia *= (2 * n + 1) * (2 * n + 2);
            xx *= x * x;
        }
    }

    /* kazdy proces wypisuje wyniki czesciowe */
    printf("%d : <%d,%d> - sinh = %Lg, cosh = %Lg\n", rank,
           przedzialyTab[0].nMin, przedzialyTab[0].nMax,
           wynikCzesciowy.sinus, wynikCzesciowy.cosinus);

    /* wysylanie wynikow czesciowych do MASTERa,
     * ktore sa sumowane w sumaWynikow 
     */
    MPI_Reduce(&wynikCzesciowy, &sumaWynikow,
               2, MPI_LONG_DOUBLE, MPI_SUM,
               0, MPI_COMM_WORLD);
    
    /* MASTER wypisuje wyniki */
    if(rank == 0) {
        printf("sinus(%.3Lg) = %Lg\n", x, sumaWynikow.sinus);
        printf("cosinus(%.3Lg) = %Lg\n",
               x, sumaWynikow.cosinus);
        printf("tangens(%.3Lg) = %Lg\n",
               x, sumaWynikow.sinus / sumaWynikow.cosinus);
        printf("czas dzialania programu: %g sekund\n",
               MPI_Wtime() - czasStart);
    }

    MPI_Finalize();
    return 0;
}

/* Funkcja zwraca, zaokraglony do liczby calkowitej, wynik
 * dzielenia dwoch liczb.
 * Gdy liczba po przecinku jest >= 0.5 to wynik zaokraglany
 * jest w gore.
 */
int dzielenie(int dzielna, int dzielnik)
{
    return (double)dzielna / dzielnik - 
           (double)(dzielna / dzielnik) >= 0.5
	   ? dzielna / dzielnik + 1 : dzielna / dzielnik;
}

/* Funkcja zwraca silnie z liczby n.
 */
long double obliczSilnie(int n)
{
    long double wynik = 1.0;

    if(n < 0) {
        printf("Blad w wysolaniu funkcji obliczSilnie !\n");
        MPI_Abort(MPI_COMM_WORLD, -1);
    }
    
    if(n == 0)
        return wynik;

    for( ; n > 0; n--)
        wynik *= n;

    return wynik;
}

/* Funkcja zwraca x podniesiony do potegi n.
 */
long double obliczPotege(long double x, int n)
{
    long double wynik = 1;

    if(n < 0) {
        printf("Blad w wywolaniu funkcji obliczPotege !\n");
        MPI_Abort(MPI_COMM_WORLD, -1);
    }

    for( ; n > 0; n--)
        wynik *= x;

    return wynik;
}



Problemy występujące w implementacji

Opisany wcześniej algorytm współbieżny polega na dekompozycji n początkowych wyrazów szeregów na, w miarę możliwości, równe części. Każda taka część jest przetwarzana przez osobny proces. Podczas obliczeń, wartości kolejnych wyrazów szeregów wyznaczane są na podstawie poprzedniego wyrazu. W momencie, gdy szereg dzielony jest na części, powstaje potrzeba obliczenia pierwszego elementu każdego nowego szeregu. Jest to zasadniczy problem uniemożliwiający uzyskanie wprost proporcjonalnego przyspieszenia, wraz ze zwiększeniem liczby procesorów.

W programie istnieje potrzeba wyznaczania wartości silni. Ponieważ funkcję tę cechuje uzyskiwanie bardzo duży wartości dla nawet niewielkich argumentów, to w obawie o niedopuszczenie do przekroczenia zakresu wartości zmiennych należy ogranicznyć liczbę sumowanych wyrazów ciągu, co implikuje zmniejszenie dokładność (choć w tym przypadku jest ona wystarczająca).

w3cw3c
automatyka przemysłowa