//---------------------------------------------------------------------------
// created:  20040808  by N.V.Shokhirev
// modified: 200408930  by N.V.Shokhirev
// modified: 20060103 - switched to real
//---------------------------------------------------------------------------
#ifndef MatUtilsH
#define MatUtilsH

#include <math.h>
#define eps0 1.0e-9

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

typedef double real;

// Macheps + 1.0 > 1.0
const real Macheps = 2.22044604925031E-16;  // for double

// SafeMacheps = 1024*Macheps + 1.0 > 1.0
const real SafeMacheps = 2.27373675443231744-13;  // for double

// SqrtMacheps = sqrt(2.22044604925031e-16)
const real SqrtMacheps = 1.4901161193848e-8;

// SafeMinreal > 0.0
const real MinReal = 4.94065645841247E-324;  // for double

const real SqrtMinReal = 2.222758749485E-162;  // for double

// SafeMinreal = 1024*Minreal
const real SafeMinReal = 5.05923221341436928e-321;

// VerySafeMinreal = 1048576*Minreal
const real VerySafeMinReal = 5.18065378653631e-318;

// Pascal square: x**2
real Sqr(real x); // { return (x*x); }

real sqr(real x); //  { return (x*x); }

// FORTRAN abs
real abs(real x); // { return ( (x >= 0.0) ? x : -x); }

int iabs(int x); // { return ( (x >= 0.0) ? x : -x); }

// FORTRAN min/max
real max1(real x1, real x2); // { return ( (x1 > x2) ? x1 : x2); }

real min1(real x1, real x2); // { return ( (x1 < x2) ? x1 : x2); }

int max0(int x1, int x2); // { return ( (x1 > x2) ? x1 : x2); }

int min0(int x1, int x2); // { return ( (x1 < x2) ? x1 : x2); }

real round1(real x); // { return floor(x+0.5);};

int round0(real x); // { return int(round1(x)); }

// Pascal signum
real Sign(real x); // Sign = 1 if x > 0; Sign = 0 if x is 0; Sign = -1 if x < 0;

int iSign (int x);

// FORTRAN signum
real sign(real x1, real x2); // { return ( x2 >= 0.0 ? abs(x1) : -abs(x1) ; }

real sign2(real a, real b);  // sign2 = sign

int isign (int x1, int x2); // { return ( x2 >= 0 ? iabs(x1) : -iabs(x1) ); }

int isign2 (int a, int b);  // isign2 = isign

// other sugnum functions
real sign0(real x); // sign1 = Sign

real sign1(real x); // {  if ( x < 0.0 ) { return -1.0; } else { return 1.0; };}

// safe sqrt: sqroot = sqrt(x + y)
real sqroot1(real x, real y);

// safe sqrt: sqroot2 = sqrt(x*x+y*y); functionally the same as pythag
real sqroot2(real x, real y);

// sqrt(a**2+b**2) without overflow or destructive underflow;
// EISPACK algorithm; functionally the same as sqroot2
real pythag(real a, real b);

// min/max of 3 variables
real max3(real x1, real x2, real x3); // { return max1(max1(x1,x3),x2); }

real min3(real x1, real x2, real x3); // { return min1(min1(x1,x3),x2); }

// if (abs(x2) > SafeMinReal) return x1/x2; else return 0.0;
real SafeDiv(real x1, real x2);

//   if ( abs(x)> eps) return x; else return 0.0;
real Zero(real x, real eps);

// float equal
bool Equal(real x1, real x2, real eps); // { return (fabs(x1-x2) < eps); };

#endif