Uitwerking van opgave 8 van paragraaf 6.18 uit het boek "De programmeertaal C", 4de editie van Kelley en Pohl.

Deze pagina is bestemd voor studenten van de Haagse Hogeschool - Academie voor Technology, Innovation & Society Delft.

Gegeven is een polynoom van graad n: f(x) = a₀ + a₁ x + a₂ x² + a₃ x³ + ... + a_n xⁿ waarbij de coëfficiënten a₀, a₁, ... , a_n reële getallen voorstellen en a_n ongelijk aan 0 is. In een C programma kunnen we de coëfficiënten opslaan in een array met (bijvoorbeeld) de naam p. Element p[0] bevat dan de coëfficiënt a₀ en element p[1] bevat coëfficiënt a₁, enz.

Gevraagd wordt om een functie genaamd eval te schrijven die de waarde van een polynoom (met graad n, waarvan de coëfficiënten opgeslagen zijn in een array p) in het punt x berekent. Het prototype van de functie (de functiedeclaratie) is gegeven:

double eval(double p[], double x, int n);

We beginnen met het schrijven van een testprogramma:

#include <stdio.h>

/* © 2012 Harry Broeders 
   Uitwerking van opgave 8 van paragraaf 6.18 uit het boek "De programmeertaal C", 4de editie van Kelley en Pohl */

double eval(double p[], double x, int n) {
    /* dummy */
    return 0;
}

void printErrorAlsOngelijk(double x, double y) {
    /* if (x != y) werkt niet vanwege afrondingsfouten ! */
    double verschil = x - y;
    double absoluutVerschil;
    if (verschil > 0) {
        absoluutVerschil = verschil;
    }
    else {
        absoluutVerschil = -verschil;
    }
    if (absoluutVerschil >= 1E-10) {
        printf("Error: verwachte waarde = %.10f maar werkelijke waarde is %.10f\n", y, x);
    }
}

int main(void) {
    double p0[] = {-3};                     /* f(x) = -3                                         */
    double p1[] = {4, -4};                  /* f(x) =  4 +   -4*x                                */
    double p2[] = {-2, -1, 1};              /* f(x) = -2 +   -1*x +    1*x^2                     */
    double p3[] = {-2, -1.5, 0.75, 0.25};   /* f(x) = -2 + -1.5*x + 0.75*x^2 +  0.25*x^3         */
    double p4[] = {3, 0.3, -0.5, 6.125, 2}; /* f(x) =  3 +  0.3*x + -0.5*x^2 + 6.125*x^3 + 2*x^4 */
    double y;

    y = eval(p0, 3.2, 0);
    printf("Als f(x) = -3 dan geldt f(3.2) = %f\n", y);
    printErrorAlsOngelijk(y, -3);
    y = eval(p1, 3.2, 1);
    printf("Als f(x) = 4 - 4*x dan geldt f(3.2) = %f\n", y);
    printErrorAlsOngelijk(y, -8.8);
    y = eval(p2, 3.2, 2);
    printf("Als f(x) = -2 - 1*x + x^2 dan geldt f(3.2) = %f\n", y);
    printErrorAlsOngelijk(y, 5.04);
    y = eval(p3, 3.2, 3);
    printf("Als f(x) = -2 - 1.5*x + 0.75*x^2 + 0.25*x^3 dan geldt f(3.2) = %f\n", y);
    printErrorAlsOngelijk(y, 9.072);
    y = eval(p4, 3.2, 4);
    printf("Als f(x) = 3 + 0.3*x - 0.5*x^2 + 6.125*x^3 + 2*x^4 dan geldt f(3.2) = %f\n", y);
    printErrorAlsOngelijk(y, 409.2592);

    fflush(stdin);
    getchar();
    return 0;
}

Om te kijken of het antwoord correct is moeten we de berekende waarde vergelijken met de verwachte waarde (die we op een rekenmachine hebben berekend). Omdat beide getallen van het type double zijn kunnen we beide getallen niet met behulp van de operator != vergelijken. Als er namelijk in de berekening een heel klein afrondfoutje is opgetreden (en dat zal zeker gebeuren) dan zal de operator != aangeven dat het antwoord niet correct is terwijl het in feite wel correct is (op een klein, niet te vermijden, afrondfoutje na). Dit probleem wordt opgelost door de berekende waarde goed te rekenen als:
| berekende waarde - verwachte waarde | < 10^-10.

Het bovenstaande programma geeft de volgende uitvoer:

Als f(x) = -3 dan geldt f(3.2) = 0.000000 Error: verwachte waarde = -3.0000000000 maar werkelijke waarde is 0.0000000000 Als f(x) = 4 - 4*x dan geldt f(3.2) = 0.000000 Error: verwachte waarde = -8.8000000000 maar werkelijke waarde is 0.0000000000 Als f(x) = -2 - 1*x + x^2 dan geldt f(3.2) = 0.000000 Error: verwachte waarde = 5.0400000000 maar werkelijke waarde is 0.0000000000 Als f(x) = -2 - 1.5*x + 0.75*x^2 + 0.25*x^3 dan geldt f(3.2) = 0.000000 Error: verwachte waarde = 9.0720000000 maar werkelijke waarde is 0.0000000000 Als f(x) = 3 + 0.3*x - 0.5*x^2 + 6.125*x^3 + 2*x^4 dan geldt f(3.2) = 0.000000 Error: verwachte waarde = 409.2592000000 maar werkelijke waarde is 0.0000000000

De functie eval kan als volgt geïmplementeerd worden:

double eval(double p[], double x, int n) {
    int i, j;
    double res = p[0];
    for (i = 1 ; i <= n; i = i + 1) {
        /* bereken x^i */
        double xTotDeMachti = x;
        for (j = 1; j < i; j = j + 1) {
            xTotDeMachti = xTotDeMachti * x;
        }
        res = res + p[i] * xTotDeMachti;
    }
    return res;
}

In het boek wordt dit de "naïeve" aanpak genoemd. De functie eval kan echter efficiënter geïmplementeerd worden door gebruik te maken van de regel van Horner. Volgens deze regel kunnen we f(x) = a₀ + a₁ x + a₂ x² + a₃ x³ + ... + a_n xⁿ ook schrijven als f(x) = a₀ + x ( a₁ + x (a₂ + x ( a₃ ... + x (a_n) ... ))). Merk op dat we in dit geval de coëfficiënten van hoog naar laag gebruiken, omdat we de haakjes van binnen naar buiten moeten uitwerken. De functie evalHorner kan als volgt geïmplementeerd worden:

double evalHorner(double p[], double x, int n) {
    int i;
    double res = p[n];
    for (i = n - 1 ; i >= 0; i = i - 1) {
        res = p[i] + x * res;
    }
    return res;
}

Het hele programma kun je hier vinden: 6.18 opgave 8.c.

Tot slot wordt in het boek de vraag gesteld: Hoeveel optellingen en vermenigvuldigingen komen er voor in elk van de twee versies van de functie (eval en evalHorner).

In de functie eval wordt f(x) berekend met de formule: f(x) = a₀ + a₁ x + a₂ x² + a₃ x³ + ... + a_n xⁿ. Dit kunnen we ook schrijven (door de machten uit te schrijven) als f(x) = a₀ + a₁ x + a₂ x x + a₃ x x x + ... + a_n x x x ... x. We zien dat het aantal optellingen bij de naieve methode gelijk is aan de graad van de polynoom (n). Het aatal vermenigvuldigingen bij de naieve methode is gelijk aan 0 + 1 + 2 + 3 + ... + n (waarbij n de graad van de polynoom is). Deze reeks is gelijk aan 0.5 * n * (n + 1). Voor een 4^de graads polynoom is het aantal vermenigvuldigingen dus 0.5 * 4 * 5 = 10.

In de functie evalHorner wordt f(x) berekend met de formule: f(x) = a₀ + x ( a₁ + x (a₂ + x ( a₃ ... + x (a_n) ... ))). We zien dat het aantal optellingen bij de methode van Horner gelijk is aan de graad van de polynoom (n). Het aatal vermenigvuldigingen bij de methode van Horner is ook gelijk aan de graad van de polynoom (n). Voor een 4^de graads polynoom is het aantal vermenigvuldigingen dus 4.

In het programma 6.18 opgave 8b.czijn instructies toegevoegd om het aantal optellingen en vermenigvuldigingen te tellen. De uitvoer van dit programma is zoals verwacht:

Naief graad = 0 : aantalOptellingen = 0, aantalVermenigvuldigingen = 0 graad = 1 : aantalOptellingen = 1, aantalVermenigvuldigingen = 1 graad = 2 : aantalOptellingen = 2, aantalVermenigvuldigingen = 3 graad = 3 : aantalOptellingen = 3, aantalVermenigvuldigingen = 6 graad = 4 : aantalOptellingen = 4, aantalVermenigvuldigingen = 10 Horner graad = 0 : aantalOptellingen = 0, aantalVermenigvuldigingen = 0 graad = 1 : aantalOptellingen = 1, aantalVermenigvuldigingen = 1 graad = 2 : aantalOptellingen = 2, aantalVermenigvuldigingen = 2 graad = 3 : aantalOptellingen = 3, aantalVermenigvuldigingen = 3 graad = 4 : aantalOptellingen = 4, aantalVermenigvuldigingen = 4

De implementatie met behulp van de regel van Horner is dus sneller omdat er (veel) minder vermenigvuldigingen nodig zijn. Dit zal vooral bij grote waarden van n opvallen. Voor n = 1000000 geld dat bij de naïeve implementatie 500000500000 vermenigvuldigingen nodig zijn terwijl bij de implementatie die gebruik maakt van de regel van Horner slechts 1000000 vermenigvuldigingen nodig zijn. Een PC met een rekenkracht van 20 GFlops doet ongeveer 25 seconden over de vermenigvuldigingen in de naïeve implementatie en 0,00005 seconden over de vermenigvuldigingen in de Horner implementatie.