//*CMZ :  2.00/07 15/05/98  16.47.13  by  Rene Brun
//*CMZ :  2.00/06 12/05/98  17.48.55  by  Fons Rademakers
//*CMZ :  2.00/01 14/03/98  05.56.51  by  Fons Rademakers
//*CMZ :  1.03/09 16/12/97  17.49.59  by  Fons Rademakers
//*-- Author :    Fons Rademakers   29/07/95

//*KEEP,CopyRight,T=C.
/*************************************************************************
 * Copyright(c) 1995-1998, The ROOT System, All rights reserved.         *
 * Authors: Rene Brun, Nenad Buncic, Valery Fine, Fons Rademakers.       *
 *                                                                       *
 * Permission to use, copy, modify and distribute this software and its  *
 * documentation for non-commercial purposes is hereby granted without   *
 * fee, provided that the above copyright notice appears in all copies   *
 * and that both the copyright notice and this permission notice appear  *
 * in the supporting documentation. The authors make no claims about the *
 * suitability of this software for any purpose.                         *
 * It is provided "as is" without express or implied warranty.           *
 *************************************************************************/
//*KEND.

//////////////////////////////////////////////////////////////////////////
//                                                                      //
// TMath                                                                //
//                                                                      //
// Encapsulate math routines (i.e. provide a kind of namespace).        //
// For the time being avoid templates.                                  //
//                                                                      //
//////////////////////////////////////////////////////////////////////////

//*KEEP,TMath.
#include "TMath.h"
//*KEND.
#include <math.h>

//const Double_t
//   TMath::Pi = 3.14159265358979323846,
//   TMath::E  = 2.7182818284590452354;


ClassImp(TMath)

//______________________________________________________________________________
#if (defined(sun) && !defined(R__I386) && !defined(__SunOS_5_6)) || \
    (defined(__OPTIMIZE__) && \
     (__GNUC__ > 2 || (__GNUC__ == 2 && __GNUC_MINOR__ > 7)))
   extern "C" void sincos(Double_t, Double_t*, Double_t*);
#else
   inline void sincos(Double_t a, Double_t *sinp, Double_t *cosp)
      { *cosp = cos(a); *sinp = sin(a); }
#endif

//______________________________________________________________________________
#if defined(R__MAC) || defined(R__KCC)
Double_t hypot(Double_t x, Double_t y) {
   return sqrt(x*x+y*y);
}
#endif

//______________________________________________________________________________
 Long_t TMath::Sqrt(Long_t x)
{
   return (Long_t) (sqrt((Double_t)x) + 0.5);
}

//______________________________________________________________________________
Long_t TMath::Hypot(Long_t x, Long_t y)     // return sqrt(px*px + py*py);
{
   return (Long_t) (hypot(x, y) + 0.5);
}

//______________________________________________________________________________
 Double_t TMath::Hypot(Double_t x, Double_t y)
{
   return hypot(x, y);
}

//______________________________________________________________________________
 Double_t TMath::ASinH(Double_t x)
{
#if defined(WIN32) || defined(R__KCC)
   return log(x+sqrt(x*x+1));
#else
   return asinh(x);
#endif
}

//______________________________________________________________________________
 Double_t TMath::ACosH(Double_t x)
{
#if defined(WIN32) || defined(R__KCC)
   return log(x+sqrt(x*x-1));
#else
   return acosh(x);
#endif
}

//______________________________________________________________________________
 Double_t TMath::ATanH(Double_t x)
{
#if defined(WIN32) || defined(R__KCC)
   return log((1+x)/(1-x))/2;
#else
   return atanh(x);
#endif
}

//______________________________________________________________________________
 Double_t TMath::Ceil(Double_t x)
{
   return ceil(x);
}

//______________________________________________________________________________
 Double_t TMath::Floor(Double_t x)
{
   return floor(x);
}

//______________________________________________________________________________
 Double_t TMath::Log2(Double_t x)
{
   return log(x)/log(2.0);
}

//______________________________________________________________________________
 Double_t TMath::Atan2(Long_t y, Long_t x)
{
   Double_t dd = (0.75 - (atan2(-(Double_t)y,(Double_t)x) / TMath::Pi() + 1.0)/2.0 ) * 360.0;
   while (dd < 0.0)
      dd += 360.0;
   return dd;
}

//______________________________________________________________________________
 void TMath::Sincos(Double_t a, Double_t *sinp, Double_t *cosp)
{
   sincos(2 * TMath::Pi() * (90.0 - a) / 360.0, sinp, cosp);
   *sinp = -*sinp;
}

//______________________________________________________________________________
 Long_t TMath::NextPrime(Long_t x)
{
   // Return next prime number after x.

   if (x <= 3)
      return 3;

   if (x % 2 == 0)
      x++;

   Long_t sqr = (Long_t) sqrt((Double_t)x) + 1;

   for (;;) {
      Long_t n;
      for (n = 3; (n <= sqr) && ((x % n) != 0); n += 2)
         ;
      if (n > sqr)
         return x;
      x += 2;
   }
}

//______________________________________________________________________________
 Float_t *TMath::Cross(Float_t v1[3],Float_t v2[3],Float_t out[3])
{
   // Calculate the Cross Product of two vectors
   //         out = [v1 x v2]

    out[0] = v1[1] * v2[2] - v1[2] * v2[1];
    out[1] = v1[2] * v2[0] - v1[0] * v2[2];
    out[2] = v1[0] * v2[1] - v1[1] * v2[0];

    return out;
}

//______________________________________________________________________________
 Float_t TMath::Normalize(Float_t v[3])
{
   // Normalize a vector v in place
   // Return:
   //    The norm of the original vector

    Float_t d = Sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
    if (d != 0)
    {
      v[0] /= d;
      v[1] /= d;
      v[2] /= d;
    }
    return d;
}

//______________________________________________________________________________
 Double_t *TMath::Cross(Double_t v1[3],Double_t v2[3],Double_t out[3])
{
   // Calculate the Cross Product of two vectors
   //         out = [v1 x v2]

    out[0] = v1[1] * v2[2] - v1[2] * v2[1];
    out[1] = v1[2] * v2[0] - v1[0] * v2[2];
    out[2] = v1[0] * v2[1] - v1[1] * v2[0];
    return out;
}

//______________________________________________________________________________
 Double_t TMath::Normalize(Double_t v[3])
{
   // Normalize a vector v in place
   //  Return:
   //    The norm of the original vector

    Double_t d = Sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
    if (d != 0)
    {
      v[0] /= d;
      v[1] /= d;
      v[2] /= d;
    }
    return d;
}

//______________________________________________________________________________
 Float_t *TMath::Normal2Plane(Float_t p1[3],Float_t p2[3],Float_t p3[3], Float_t normal[3])
{
   // Calculate a normal vector of a plane
   //
   //  Input:
   //     Float_t *p1,*p2,*p3  -  3 3D points belonged the plane to define it.
   //
   //  Return:
   //     Pointer to 3D normal vector (normalized)
   //

     Float_t v1[3], v2[3];

     v1[0] = p2[0] - p1[0];
     v1[1] = p2[1] - p1[1];
     v1[2] = p2[2] - p1[2];

     v2[0] = p3[0] - p1[0];
     v2[1] = p3[1] - p1[1];
     v2[2] = p3[2] - p1[2];

     NormCross(v1,v2,normal);
     return normal;
}

//______________________________________________________________________________
 Double_t *TMath::Normal2Plane(Double_t p1[3],Double_t p2[3],Double_t p3[3], Double_t normal[3])
{
   // Calculate a normal vector of a plane
   //
   //  Input:
   //     Float_t *p1,*p2,*p3  -  3 3D points belonged the plane to define it.
   //
   //  Return:
   //     Pointer to 3D normal vector (normalized)
   //

     Double_t v1[3], v2[3];

     v1[0] = p2[0] - p1[0];
     v1[1] = p2[1] - p1[1];
     v1[2] = p2[2] - p1[2];

     v2[0] = p3[0] - p1[0];
     v2[1] = p3[1] - p1[1];
     v2[2] = p3[2] - p1[2];

     NormCross(v1,v2,normal);
     return normal;
}