#include "EXTERN.h"
#include "perl.h"
#include "XSUB.h"

#include <math.h>
#define max(x,y) (((x)>(y)) ? (x) : (y))

MODULE = ST::Xten       PACKAGE = ST::Xten      

# Dieser Code wird zur Initialisierungszeit ausgefuehrt.
# Das ist nuetzlich, um einige Libraries zu initialisieren.
# Wir geben nur eine einfache Meldung aus.
BOOT:
printf("Initializing the ST::Xten Library (faking...)\n");

# t_norm1() liefert die Summennorm eines dreidimensionalen
#     Vektors, der durch drei Koordinaten angegeben wird.
double
t_norm1(x,y,z)
    double  x;
    double  y;
    double  z;
    CODE:
    RETVAL = fabs(x) + fabs(y) + fabs(z);
        OUTPUT:
    RETVAL

# t_norm2() liefert die euklidische Norm eines
#   dreidimensionalen Vektors, der durch
#   drei Koordinaten angegeben wird.
double
t_norm2(x,y,z)
    double  x;
    double  y;
    double  z;
    CODE:
    RETVAL = sqrt(x*x + y*y + z*z);
    OUTPUT:
    RETVAL

# t_norminf() liefert die Maximumnorm eines
#       dreidimensionalen Vektors, der
#       durch drei Koordinaten angegeben wird.
double
t_norminf(x,y,z)
    double  x;
    double  y;
    double  z;
    CODE:
    RETVAL = max(max(fabs(x),fabs(y)),fabs(z));
    OUTPUT:
    RETVAL

# v_norm1() liefert die Summennorm eines n-dimensionalen Vektors,
#     der durch eine Koordinatenliste angegeben wird.
double
v_norm1(x1, ...)
    double      x1;
    PREINIT:
    unsigned int    i;
    double      result;
    CODE:
    result = fabs(x1);
    for (i=1; i<items; i++) {
      result += fabs((double)SvNV(ST(i)));
    }
    RETVAL = result;
    OUTPUT:
    RETVAL

# v_norm2() liefert die euklidische Norm eines n-dimensionalen
#     Vektors, der durch eine Koordinatenliste angegeben wird.
double
v_norm2(x1, ...)
    double      x1;
    PREINIT:
    unsigned int    i;
    double      result;
    CODE:
    result = x1*x1;
    for (i=1; i<items; i++) {
      /* VORSICHT! Das kann sehr gross werden! */
      result += (double)(SvNV(ST(i))) * (double)(SvNV(ST(i)));
    }
    result = sqrt(result);
    RETVAL = result;
    OUTPUT:
    RETVAL

# v_norminf() liefert die Maximumnorm eines n-dimensionalen Vektors,
#       der durch eine Koordinatenliste angegeben wird.
double
v_norminf(x1, ...)
    double      x1;
    PREINIT:
    unsigned int    i;
    double      result;
    CODE:
    result = x1;
    for (i=1; i<items; i++) {
      result = max(result, fabs((double)SvNV(ST(i))));
        }
    RETVAL = result;
    OUTPUT:
    RETVAL

# v_normalize() normalisiert einen n-dimensionalen Vektor
#         bezueglich der euklidischen Norm.
#         Der Vektor ist durch eine Koordinatenliste
#         spezifiziert, und es wird eine Rueckgabeliste
#         zurueckgegeben.
double*
v_normalize(x1, ...)
    double      x1;
    PREINIT:
    unsigned int    i;
    double *    normvect;
    double      norm2;
    PPCODE:

    /* Erst die euklidische Norm ausrechnen.        */
    /* Wir koennen leider nicht v_norm2() aufrufen! */
    norm2 = x1*x1;
    for (i=1; i<items; i++) {
      /* VORSICHT! Das kann sehr gross werden! */
      norm2 += (double)(SvNV(ST(i))) * (double)(SvNV(ST(i)));
    }
    norm2 = sqrt(norm2);
    if (!norm2) croak("can't normalize with 0 norm!\n");

    /* Nun koennen wir normalisieren. */
    /* Hinweis: New() ist das malloc() von Perl */
    New(0, normvect, items, double);
    normvect[0] = x1/norm2;
    for (i=1; i<items; i++) {
      normvect[i] = (double)(SvNV(ST(i)))/norm2;
        }

    /* Wir pushen die Ergebniswerte auf den Stack */
    for (i=0; i<items; i++) {
      XPUSHs(sv_2mortal(newSVnv(normvect[i])));
        }

    /* Und wir raeumen unseren Zwischenvektor wieder auf */
    /* Hinweis: Safefree() ist das free() von Perl */
    Safefree(normvect);
    XSRETURN(items);