Startseite C / C++ Borland C++ Builder Archive nach Kategorien "Mathematik"

Kategorie: Mathematik

[BC++] Dritte Wurzel aus einer negativen Zahl (pow domain error) ziehen

Sicherlich ist es schon vielen so ergangen: Excel rechnet es, der Taschenrechner auch, bloß im C++-Code ist der Wurm drin. Beim Ziehen der dritten Wurzel aus einer negativen Zahl mittels pow-Funktion (x^1/3) poppt ein pow domain error auf:

#include <math.h>

double x = -27.0;
double y = pow(x, (1.0 / 3.0)); // <-- pow domain error

Die Lösung ist rel. einfach. Für ungerade Wurzel-Exponenten kann man die Wurzel aus dem Absolutwert der Basis ziehen und das Ergebnis dann mit dem Vorzeichenwert nochmal multiplizieren:

#include <math.h>

double x = -27.0;
double ax = fabs(x);
double y = pow(ax, (1.0 / 3.0)) * x / ax;

Das Problem besteht in der Umsetzung der pow-Funktion über Logarithmen:

pow(x,y) = (e^ln(x))^y = e^(y*ln(x))

d.h.

pow(-1,y) = e^(y*ln(-1)) --> der nat. Logarithmus akzeptiert keine neg. Argumente oder 0 --> pow domain error

Potenzieren 

#include <math.h>

...

void __fastcall TMain::Button1Click(TObject *Sender)
{
     double dbl;
     dbl = pow (2,4); // 2^4 ausrechnen
     Application->MessageBox(AnsiString(dbl).c_str(),"",MB_OK);
}

Matrizeninversion mit dem Gauss-Jordan-Algorithmus

  • Aufrufbeispiel:
#include "matinvert.h"
.
.
.
void __fastcall TForm1::Button1Click(TObject *Sender)
{
    int iSize = 3;
    double Determinante;

    double Matrix[iSize][iSize]={{5,5,4},
                                 {9,3,2},
                                 {3,8,7}};

    // invertierte Matrix steht nach der Inversion wieder in 'Matrix'
    // zusätzlich wird die Determinante der Matrix zurückgegeben

    Determinante = invert_matrix(*Matrix, iSize);
}
  • matinvert.h
#ifndef MATINVERT_H
#define MATINVERT_H

#define MAT_AT(A, COLW, J, I) (A[(J)*(COLW)+(I)])
double invert_matrix(double *A, int dim);

#endif
  • matinvert.cpp
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "matinvert.h"

double invert_matrix(double *A, int dim)
{
    int i, j, k;
    int *pivot_i, *pivot_j; /* Ort des Pivot */
    double pivot; /* das Pivot selbst */
    double tmp;
    double det = 1.0;

    pivot_i = (int*)malloc(dim*sizeof(int));
    pivot_j = (int*)malloc(dim*sizeof(int));

    for (k=0; k<dim; k++)
    {
        /* k-tes Pivot finden */
        pivot = MAT_AT(A, dim, k, k); /* suche init. */
        pivot_i[k] = k;
        pivot_j[k] = k;
        for (i=k; i<dim; i++)
            for (j=k; j<dim; j++)
                if (fabs(MAT_AT(A, dim, i, j)) -> fabs(pivot))
                {

                   pivot_i[k] = i;
                   pivot_j[k] = j;
                   pivot = MAT_AT(A, dim, i, j);

                }

        /* aufmultiplizierte Pivots ergeben spaeter die Determinante */
        det *= pivot;

        if (det == 0.0)
        {
            /* singulaere Matrix */
            free(pivot_i);
            free(pivot_j);
            return 0.0;
        }

        /* Zeilen vertauschen, Vorzeichen beachten */
        i = pivot_i[k];
        if (i!=k) /* wenn sich die Zeilen unterscheiden ... */
           for (j=0; j<dim; j++)
           {
               tmp = -MAT_AT(A, dim, k, j);
               MAT_AT(A, dim, k, j) = MAT_AT(A, dim, i, j);
               MAT_AT(A, dim, i, j) = tmp;
           }

        /* Spalten vertauschen */
        j = pivot_j[k];
        if (j!=k) /* wenn sich die Spalten unterscheiden ... */
           for (i=0; i<dim; i++)
           {
               tmp = -MAT_AT(A, dim, i, k);
               MAT_AT(A, dim, i, k) = MAT_AT(A, dim, i, j);
               MAT_AT(A, dim, i, j) = tmp;
           }

        /* Spalte durch -pivot teilen */
        for (i=0; i<dim; i++)
            if (i!=k) /* Pivot nicht einbeziehen */
               MAT_AT(A, dim, i, k) /= -pivot;

        /* Matrix reduzieren ... */
        for (i=0; i<dim; i++)
        {
            tmp = MAT_AT(A, dim, i, k);
            for (j=0; j<dim; j++)
                if (i!=k && j!=k) /* Pivot nicht einbeziehen */
                   MAT_AT(A, dim, i, j) += tmp * MAT_AT(A, dim, k, j);
        }

        /* Zeile durch Pivot teilen */
        for (j=0; j<dim; j++)
            if (j!=k) /* Pivot nicht einbeziehen */
               MAT_AT(A, dim, k, j) /= pivot;

        /* Pivot durch seinen Kehrwert ersetzen */
        MAT_AT(A, dim, k, k) = 1./pivot;

    }

    /* In einem weiteren Durchlauf Zeilen/Spalten vertauschen */
    for (k = dim-2; k ->= 0; k--)
    {
        i = pivot_j[k]; /* mit Pivot-Spalte korrespondierende Zeilen */
        if (i!=k) /* wenn sich die Zeilen unterscheiden ... */
           for (j = 0; j < dim; j++)
           {
               tmp = MAT_AT(A, dim, k, j);
               MAT_AT(A, dim, k, j) = -MAT_AT(A, dim, i, j);
               MAT_AT(A, dim, i, j) = tmp;
           }

        j = pivot_i[k]; /* mit Pivot-Zeile korrespondierende Spalten */
        if (j!=k) /* wenn sich die Zeilen unterscheiden ... */
           for (i=0; i<dim; i++)
           {
               tmp = MAT_AT(A, dim, i, k);
               MAT_AT(A, dim, i, k) = -MAT_AT(A, dim, i, j);
               MAT_AT(A, dim, i, j) = tmp;
           }
    }

    free(pivot_i);
    free(pivot_j);
    return det; /* Determinante zurückgeben */
}

Schneller rechnen

  • durch Verwendung der in fastmath.h definierten Rountinen
  • der Geschwindigkeitsvorteil ergibt sich aus den fehlenden Fehlerabfragen und Fehlerbehandlungsroutinen
  • einbinden mit #include “fastmath.h”
  • Aufruf der Funktionen mit der vorangestellten Zeichenfolge _fm_
  • Beispiel: int z = _fm_tan(Wert);