//---------------------------------------------------------------------------
// created:  20040808  by N.V.Shokhirev
// modified: 20050809  by N.V.Shokhirev
// modified: 20060103 - switched to real
//---------------------------------------------------------------------------
#include <vcl.h>
#pragma hdrstop

#include "MatUtils.h"

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

#pragma package(smart_init)

// Pascal square
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)
{
  if ( x > 0.0 ) { return 1.0; }
  else if ( x < 0.0 ) { return -1.0; }
  else { return 0.0; };
}

int iSign (int x)
{
  if (x > 0) { return 1; }
  else if (x == 0)  { return 0; }
  else { return -1; };
}

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

real sign2(real a, real b)   // sign2 = sign
{
  if ( b < 0.0 ) { return -abs(a); }
  else { return abs(a); };
}

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

int isign2 (int a, int b)    // isign2 = isign
{
  if ( b < 0 ) { return -iabs(a); }
  else { return iabs(a); };
}

// other sugnum functions
real sign0(real x)
{
  if ( x > 0.0 ) { return 1.0; }
  else if ( x < 0.0 ) { return -1.0; }
  else { return 0.0; };
}

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

// safe sqrt: sqroot = sqrt(x + y)
real sqroot1(real x, real y)
{
  if ((x + y) <= 0) { return 0.0; }
  else if (x > y) { return sqrt(x)*sqrt(1.0 + y/x); }
  else { return sqrt(y)*sqrt(1.0 + x/y); };
}

// safe sqrt: sqroot2 = sqrt(x*x+y*y); functionally the same as pythag
real sqroot2(real x, real y)
{
  if (abs(x) > abs(y)) { return abs(x)*sqrt(1.0+sqr(y/x)); }
  else if (abs(y) == 0.0) { return 0.0; }
  else { return abs(y)*sqrt(1.0+sqr(x/y)); };
}

// finds dsqrt(a**2+b**2) without overflow or destructive underflow;
// EISPACK algorithm; functionally the same as sqroot2
real pythag(real a, real b)
{
  real p, r;
  p = max1(abs(a),abs(b));
  if (p == 0.0) { return 0.0; };
  r = min1(abs(a),abs(b));
  if (r < SqrtMinReal*p) { return p; }
  else { return p*sqrt(1.0 + sqr(r/p)); };
}

// 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)
{ return ( (abs(x2) > SafeMinReal) ? x1/x2 : 0.0 ); }

//   if ( abs(x)> eps) return x; else return 0.0;
real Zero(real x, real eps)
{ return ( (abs(x)> eps) ? x : 0.0 ); }

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