#define PI 3.141592653589793238462643
#endif
-#include <iostream.h>
-#include <math.h>
+
+#include <TMath.h>
+#include <Riostream.h>
+//#include <math.h>
#include <stdlib.h>
/////////////////////////// Helper Error Function
-void ComplexError (char *msg)
+void ComplexError (const char *msg)
{
cout << endl << msg << endl;
exit (1);
}
-void ComplexWarning (char *msg)
+void ComplexWarning (const char *msg)
{
cout << endl << msg << endl;
}
////////////////////////////////////////////////
-inline double sqr(double a) { return a*a; }
inline Complex Complex::operator + () const
{
double operator += (double &a, const Complex &b)
{
- if (fabs(b.Im) < small_epsilon) a+=b.Re;
+ if (TMath::Abs(b.Im) < small_epsilon) a+=b.Re;
else
ComplexError ("Error in double+=Complex: Complex is not double number");
return a;
double operator -= (double &a, const Complex &b)
{
- if (fabs(b.Im) < small_epsilon) a-=b.Re;
+ if (TMath::Abs(b.Im) < small_epsilon) a-=b.Re;
else
ComplexError ("Error in double -=Complex: Complex is not double number");
return a;
double operator *= (double &a, const Complex &b)
{
- if (fabs(b.Im) < small_epsilon) a*=b.Re;
+ if (TMath::Abs(b.Im) < small_epsilon) a*=b.Re;
else
ComplexError ("Error in double *=Complex: Complex is not double number");
return a;
double operator /= (double &a, const Complex &b)
{
- if (fabs(b.Im) < small_epsilon) a/=b.Re;
+ if (TMath::Abs(b.Im) < small_epsilon) a/=b.Re;
else
ComplexError ("Error in double /=Complex: Complex is not double number");
return a;
Complex Complex::operator /= (const Complex &a)
{
double t1, t2, temp;
- if (fabs(a.Re) >= fabs(a.Im))
+ if (TMath::Abs(a.Re) >= TMath::Abs(a.Im))
{
t1 = a.Im / a.Re;
t2 = a.Re + a.Im * t1;
inline Complex operator / (const Complex &a, const Complex &b)
{
double t1, t2;
- if (fabs(b.Re) >= fabs(b.Im))
+ if (TMath::Abs(b.Re) >= TMath::Abs(b.Im))
{
t1 = b.Im / b.Re;
t2 = b.Re + b.Im * t1;
inline int operator == (const Complex &a, double b)
{
- return (a.Re==b) && (fabs(a.Im)<small_epsilon);
+ return (a.Re==b) && (TMath::Abs(a.Im)<small_epsilon);
}
inline int operator == (double b, const Complex &a)
{
- return (a.Re==b) && (fabs(a.Im)<small_epsilon);
+ return (a.Re==b) && (TMath::Abs(a.Im)<small_epsilon);
}
inline int operator == (const Complex &a, const Complex &b)
inline int operator != (const Complex &a, double b)
{
- return (a.Re!=b) || (fabs(a.Im)>small_epsilon);
+ return (a.Re!=b) || (TMath::Abs(a.Im)>small_epsilon);
}
inline int operator != (double b, const Complex &a)
{
- return (a.Re!=b) || (fabs(a.Im)>small_epsilon);
+ return (a.Re!=b) || (TMath::Abs(a.Im)>small_epsilon);
}
inline int operator != (const Complex &a, const Complex &b)
}
double abs (const Complex &a)
{
- if (a.Im == 0) return fabs(a.Re);
- if (a.Re == 0) return fabs(a.Im);
- double R = fabs(a.Re), I = fabs(a.Im);
+ if (a.Im == 0) return TMath::Abs(a.Re);
+ if (a.Re == 0) return TMath::Abs(a.Im);
+ double R = TMath::Abs(a.Re), I = TMath::Abs(a.Im);
return (R >= I) ?
- (R * sqrt (1 + sqr(a.Im/a.Re))) :
- (I * sqrt (1 + sqr(a.Re/a.Im))) ;
+ (R * sqrt (1 + a.Im*a.Im/a.Re/a.Re)) :
+ (I * sqrt (1 + a.Re*a.Re/a.Im/a.Im)) ;
}
double Arg (const Complex &a)
Complex sqrt (const Complex &a, int flag)
{
- if ((a.Re>=0) && (fabs(a.Im)<small_epsilon))
+ if ((a.Re>=0) && (TMath::Abs(a.Im)<small_epsilon))
return flag ? -sqrt(a.Re) : sqrt(a.Re);
- double R = fabs(a.Re), I = fabs(a.Im);
+ double R = TMath::Abs(a.Re), I = TMath::Abs(a.Im);
double w = (R >= I) ?
- sqrt (R/2 * ( 1 + sqrt (1 + sqr(a.Im/a.Re)))):
- sqrt (I/2 * (R/I + sqrt (1 + sqr(a.Re/a.Im))));
+ sqrt (R/2 * ( 1 + sqrt (1 + a.Im*a.Im/a.Re/a.Re))):
+ sqrt (I/2 * (R/I + sqrt (1 + a.Re*a.Re/a.Im/a.Im)));
Complex c;
if (a.Re >= 0)
{
t = sqrt((r^2) - (q^3));
a = ((!r * t) >=0.0) ? -((r+t)^(1.0/3)) : -((r-t)^(1.0/3));
b = ((a == zero) ? zero : (q/a));
- z[0] = -(a+b)/2 - a2/3 + ImUnit*sqrt(3)*(a-b)/2;
- z[1] = -(a+b)/2 - a2/3 - ImUnit*sqrt(3)*(a-b)/2;
+ z[0] = -(a+b)/2 - a2/3 + ImUnit*sqrt(3.0)*(a-b)/2;
+ z[1] = -(a+b)/2 - a2/3 - ImUnit*sqrt(3.0)*(a-b)/2;
z[2] = a + b - a2/3;
}