/* $Id: numerics.h,v 1.2 1996/10/04 15:07:13 roitzsch Exp $ */

#ifndef NUMERICS_H
#define NUMERICS_H

#include "general.h"
#include "kvector.h"


const double DoublePrec = 1e-15;
const float  FloatPrec  = 1e-7;

const double MaxDouble  = 1e+300;
const double MinDouble  = 1e-300;
const float  Maxfloat   = 1e+37;
const float  Minfloat   = 1e-37;

inline float  machPrec(float /*dummy*/)  { return FloatPrec;  }
inline double machPrec(double /*dummy*/) { return DoublePrec; }

inline double machMax(double /*dummy*/) { return MaxDouble; }
inline float  machMax(float  /*dummy*/) { return Maxfloat;  }
inline double machMin(double /*dummy*/) { return MinDouble; }
inline float  machMin(float  /*dummy*/) { return Minfloat;  }

//inline int machMax(int dummy)  { return int(pow(2,64)); }
//inline int machMin(int dummy)  { static int i = int(pow(2,64));  return i+1; }

#ifndef MACOS
inline int Round(double x) { return int(x + 0.5); }
#endif

inline float quot(float x, float y)
{ if (!y==0) return x/y; else return (x==0 ? 0 : Maxfloat); }
inline double quot(double x, double y)
{ if (!y==0) return x/y; else return (x==0 ? 0 : MaxDouble); }

//---------------------  Real/Complex arithmetic  -------------------------


inline double Abs(Complex x){ return sqrt(real(x)*real(x) + imag(x)*imag(x)); }

inline double Abs(double x) { return (x < 0.0)?  -x : x; }
inline float  Abs(float  x) { return (x < 0.0)?  -x : x; }
inline short  Abs(short  x) { return (x < 0.0)?  -x : x; }
inline int    Abs(int    x) { return (x < 0.0)?  -x : x; }

inline int sign(double x)   { return (x < 0.0)?  -1 : 1; }
inline int sign(float  x)   { return (x < 0.0)?  -1 : 1; }
inline int sign(int    x)   { return (x < 0.0)?  -1 : 1; }
inline int sign(short  x)   { return (x < 0.0)?  -1 : 1; }

#ifndef __DECCXX
inline Complex sqr(Complex x) { return x*x; }
#endif
inline double  sqr(double x)  { return x*x; }
inline float   sqr(float  x)  { return x*x; }
inline int     sqr(int    x)  { return x*x; }
inline short   sqr(short  x)  { return x*x; }

inline double 	Max(double x, double y) { return (x>y) ? x : y; }
inline float 	Max(float x, float y) 	{ return (x>y) ? x : y; }
inline int 	Max(int x, int y) 	{ return (x>y) ? x : y; }
inline double 	Min(double x, double y) { return (x<y) ? x : y; }
inline float 	Min(float x, float y) 	{ return (x<y) ? x : y; }
inline int 	Min(int x, int y) 	{ return (x<y) ? x : y; }

//-------------------------------------------------------------------------

inline float  conj(float  x) { return x; }	// dummies for overloading
inline double conj(double x) { return x; }

inline double  arg(float  /*x*/) { return 0.0; }	// dummies for overloading
inline double  arg(double /*x*/) { return 0.0; }

inline float  real(float  x) { return x; }	// dummies for overloading
inline double real(double x) { return x; }
inline float  imag(float  /*x*/) { return 0.0; }	// dummies for overloading
inline double imag(double /*x*/) { return 0.0; }

//inline int operator == (Complex c1, Complex c2) { return 0; }	// dummies
inline int operator <  (Complex /*c1*/, Complex /*c2*/) { return 0; }	
inline int operator >  (Complex /*c1*/, Complex /*c2*/) { return 0; }
inline int operator <= (Complex /*c1*/, Complex /*c2*/) { return 0; }
inline int operator >= (Complex /*c1*/, Complex /*c2*/) { return 0; }
inline int equal       (Complex /*c1*/, Complex /*c2*/) { return 0; }

//-------------------------------------------------------------------------
//-------------------------------------------------------------------------

template<class T> inline T dot(const Vector<T>&  x, const Vector<T>&  y)
{
    int i;
    T sum = 0.0;
    for (i=x.low(); i<=x.high(); ++i) sum += x[i] * y[i];
    return sum;
};

template<class T> inline T cdot(const Vector<T>&  x, const Vector<T>&  y)  
								 // conj(x)*y
{
    int i;
    T sum = 0.0;
    for (i=x.low(); i<=x.high(); ++i) sum += conj(x[i]) * y[i];
    return sum;
};
//-------------------------------------------------------------------------

template<class T> inline void CrossProduct(Vector<T>& z, const Vector<T>& x,
					   const Vector<T>& y)
{
    z[1] = x[2]*y[3] - x[3]*y[2];
    z[2] = x[3]*y[1] - x[1]*y[3];
    z[3] = x[1]*y[2] - x[2]*y[1];
};
								    // z = x^y

inline float  Norm(const Vector<float>&   x) { return sqrt(dot(x,x)); }
inline double Norm(const Vector<double>&  x) { return sqrt(dot(x,x)); }
       double Norm(const Vector<Complex>& x);

//-------------------------------------------------------------------------

int equal(double x, double y);
int equal(float  x, float  y);

//   int operator == (double x, double y);	// does not work
//   int operator != (double x, double y);
//   int operator == (float  x, float  y);
//   int operator != (float  x, float  y);

//-------------------------------------------------------------------------

// template<class T>  Norm(T x, T y) { return sqrt(x*x + y*y); }

//-------------------------------------------------------------------------


#endif