/*
 $Id: numerics.cc,v 1.1.1.1 1996/10/02 10:35:43 roitzsch Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "numerics.h"


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


double Norm(const Vector<Complex>& x)
{
    int i;
    double sum = 0.0;
    for (i=x.low(); i<=x.high(); ++i) 
      		 sum += (real(x[i])*real(x[i]) + imag(x[i])*imag(x[i]));
    return sqrt(sum);
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


//-------------------------------------------------------------------------
/*
   float machPrec(float dummy)
   {
   static int calls = 0;
   static float eps = 1.0;
   
   if (calls==0)
   {
   do { eps *= 0.5; }  while (float(1.0+eps) != float(1.0));
   eps *= 2.0;
   calls = 1;
   }
   // cout << "\n 1+eps      = "  << form("%1.9f",float(1.0 + eps)) << "\n";
   // cout << "\n 1+0.5*eps  = "  << form("%1.9f",float(1.0 + 0.5*eps)) << "\n";
   // cout << "\n eps = "  << eps  << "\n";
   return eps;
   }
   
   double machPrec(double dummy)
   {
   static int  calls = 0;
   static double eps = 1.0;
   
   if (calls==0)
   {
   do { eps *= 0.5; }  while (double(1.0+eps) != double(1.0));
   eps *= 2.0;
   calls = 1;
   }
   return eps;
   }
   //-------------------------------------------------------------------------
   
   int machMax(int dummy)  { return int(pow(2,64)); }
   int machMin(int dummy)  { static int i = int(pow(2,64));  return i+1; }
   
   double machMax(double dummy)
   {
   static int calls = 0;
   static double max = 1e20;
   double maxp1 = 10.*max, maxp2 = 100.*max;
   
   if (calls == 0)
   {
   do { maxp2 *= 10; maxp1 *= 10; max *= 10; }  while (maxp1 < maxp2); 
   calls = 1;
   }
   return max;
   }
   float machMax(float dummy)
   {
   static int calls = 0;
   static float max = 1e20;
   float maxp1 = 10.*max, maxp2 = 100.*max;
   
   if (calls == 0)
   {
   do { maxp2 *= 10; maxp1 *= 10; max *= 10; }  while (maxp1 < maxp2); 
   calls = 1;
   }
   return max;
   }
   
   double machMin(double dummy)
   {
   static int calls = 0;
   static double min = 1e-20;
   double minp1 = 0.1*min, minp2 = 0.01*min;
   
   if (calls == 0)
   {
   do { minp2 *= 0.1; minp1 *= 0.1; min *= 0.1; }  while (minp1 > minp2); 
   calls = 1;
   }
   return min;
   }
   float machMin(float dummy)
   {
   static int calls = 0;
   static float min = 1e-20;
   float minp1 = 0.1*min, minp2 = 0.01*min;
   
   if (calls == 0)
   {
   do { minp2 *= 0.1; minp1 *= 0.1; min *= 0.1; }  while (minp1 > minp2); 
   calls = 1;
   }
   return min;
   }
*/
//-------------------------------------------------------------------------

int equal(double x, double y) 
{
    static double tiny = machMin (double(0.0));
    static double prec = machPrec(double(0.0));
    double xAbs = Abs(x);

    if (xAbs <= tiny) return (Abs(y)   <= tiny);  	// x approx.== 0 
    else  	      return (Abs(x-y) <= prec*xAbs);  
}

int equal(float x, float y) 
{
    static float tiny = machMin (float(0.0));
    static float prec = machPrec(float(0.0));
    float xAbs = Abs(x);

    if (xAbs <= tiny) return (Abs(y)   <= tiny);  	// x approx.== 0 
    else  	      return (Abs(x-y) <= prec*xAbs);  
}