Source: classlib/include/nmath.h

#ifndef NMATH__
#define NMATH__

#include 
#include 
#include 
#include 

#include 

#define EPS 0.001
#define DPI ( M_PI * 2.0 )
#define EquEps( a, b, eps )	(fabs((a)-(b)) < (eps))
#define ZeroEps( a, eps )	(fabs(a) < (eps))
#define EquEPS( a, b )	(fabs((a)-(b)) < EPS)
#define ZeroEPS( a )	(fabs(a) < EPS)
#define DimMat (sizeof( double ) * 16)
#define Freccia 1.0
#define MAXNPTI 1000
#define BIG_FLOAT	10e7

// conversione da gradi a radianti
#define G2R(x)	(( ((double)(x)) * M_PI) / 180.0)
// conversione radianti in gradi
#define R2G(x)	(( ((double)(x)) * 180.0) / M_PI)

// informazioni per funzioni di elaborazione archi
struct ArcInfo
{
  float cx, cy, raggio;
  double angStart, angStop;
};

#ifdef __cplusplus

#ifndef SWAPFUN__
#define SWAPFUN__

// funzione template di swapping
template  inline void Swap( T &t1, T &t2 )
{
	T tmp;
	tmp = t1;
	t1  = t2;
	t2  = tmp;
}

// in t1 sempre il maggiore dei due
template  inline void maxSwap( T &t1, T &t2 )
{
	if( t2 > t1 ) {
		T tmp;
		tmp = t1;
		t1  = t2;
		t2  = tmp;
	}
}

// in t1 sempre il minore dei due
template  inline void minSwap( T &t1, T &t2 )
{
	if( t2 < t1 ) {
		T tmp;
		tmp = t1;
		t1  = t2;
		t2  = tmp;
	}
}

#endif // SWAPFUN__

class Linea2D;
class Linea3D;

class Matrix3
{
public:
    double m[4][4];

    Matrix3();
    Matrix3( const double ma[4][4] ) { copymatrix( ma ); }
    Matrix3( const float ma[4][4] ) { copymatrix( m, ma ); }
    Matrix3( const Matrix3 &ma ) { copymatrix( ma.m ); }
    virtual ~Matrix3() {}
    void MultMatrix( Matrix3& );
    void MultMatrix( double ma[4][4] );
    void ReadMatrix( void *m1 ) { memcpy( m1, m, DimMat ); }
    void WriteMatrix( void *m1 ) { memcpy( m, m1, DimMat ); }
    void Clear0() { memset( m, 0, DimMat ); }
    void ClearI();
    void ScalaMatrix( double x, double y, double z );
    void TraslaMatrix( double x, double y, double z );
    void RotateXMatrix( double angolo );
    void RotateYMatrix( double angolo );
    void RotateYMatrixCart( double angolo );
    void RotateZMatrix( double angolo );
    void RotateXMatrix( double s, double c );
    void RotateYMatrix( double s, double c );
    void RotateYMatrixCart( double s, double c );
    void RotateZMatrix( double s, double c );
    void MultPoint( double &x, double &y, double &z ) const;
    void MultPoint( float &x, float &y, float &z ) const;
    void MultPoint( int &x, int &y, int &z ) const;
    void MultPoint( double p[4] ) const;
    void GetIdentita( void *mi );

    void operator = ( const Matrix3& m2 ) { copymatrix( m2.m ); }
    void operator = ( const double ma[4][4] ) { copymatrix( ma ); }
    double& operator() ( int a, int b ) { return m[a][b]; }
   	int operator == ( const Matrix3& m ) const {return &m == this;}

protected:
    void copymatrix( double m1[4][4], const double m2[4][4] ) { memcpy( m1, m2, DimMat ); }
    void copymatrix( double m1[4][4], const float m2[4][4] );
    void copymatrix( const double m1[4][4] ) { memcpy( m, m1, DimMat ); }

private:
    static const double Identita[4][4];
};

class Punto2D
{
friend class Linea2D;
public:
    float x, y;

    Punto2D() { x=y=0.0f; }
    Punto2D( float xx, float yy ) { x = xx; y = yy; }
    Punto2D( const Punto2D &p ) { x = p.x; y = p.y; }
    virtual ~Punto2D() {}
    double GetDistance( const Punto2D& p ) const;
	Punto2D& AccumulaMin( const Punto2D& p );
	Punto2D& AccumulaMax( const Punto2D& p );
	Punto2D& AccumulaMin( float x, float y );
	Punto2D& AccumulaMax( float x, float y );
    void MultMatrixXY( Matrix3& m );
    void MultMatrixXZ( Matrix3& m );
    Punto2D& operator = ( const Punto2D &p ) { x=p.x; y=p.y; return *this; }
    Punto2D  operator - ( const Punto2D& p2 ) const;
    Punto2D  operator + ( const Punto2D& p2 ) const;
    Punto2D  operator * ( const double scalare ) const;
    Punto2D& operator -= ( const Punto2D& p2 );
    Punto2D& operator += ( const Punto2D& p2 );
    int operator == ( const Punto2D& p2 ) const;
    int operator != ( const Punto2D& p2 ) const;

    static void SetEpsilon( float eps );
    static void GetEpsilon( float&, float& );

private:
    static float epsilonX;
	static float epsilonY;
};


class Punto3D
{
friend class Linea3D;
public:
    float x, y, z;

    Punto3D() { x=y=z=0.0f; }
    Punto3D( float xx, float yy, float zz ) { x = xx; y = yy; z = zz; }
    Punto3D( const Punto3D &p ) { x = p.x; y = p.y; z = p.z; }
    virtual ~Punto3D() {}
    double GetDistance( const Punto3D& p ) const;
    double GetMagnitude( const Punto3D& p ) const;
	Punto3D& AccumulaMin( const Punto3D& p );
	Punto3D& AccumulaMax( const Punto3D& p );
	Punto3D& AccumulaMin( float x, float y, float z );
	Punto3D& AccumulaMax( float x, float y, float z );
    void MultMatrix( Matrix3& m );
    Punto3D& operator = ( const Punto3D &p ) { x=p.x; y=p.y; z=p.z; return *this; }
    Punto3D  operator - ( const Punto3D& p2 ) const;
    Punto3D  operator + ( const Punto3D& p2 ) const;
    Punto3D  operator * ( const double scalare ) const;
    Punto3D  operator * ( const Punto3D &p ) const;
    Punto3D& operator -= ( const Punto3D& p2 );
    Punto3D& operator += ( const Punto3D& p2 );
    int operator == ( const Punto3D& p2 ) const;
    int operator != ( const Punto3D& p2 ) const;

    static void SetEpsilon( float eps );
    static void GetEpsilon( float&, float&, float& );

private:
    static float epsilonX;
	static float epsilonY;
	static float epsilonZ;
};

class Linea2D
{
private:
    static float epsilonAng;
protected:
    Punto2D p1, p2, b1, b2;
    double coeff;
    char Tipo;

public:
    Linea2D();
    Linea2D( const Punto2D& pa, const Punto2D& pb );
    Linea2D( const Punto3D& pa, const Punto3D& pb );
    Linea2D( const Punto3D& pa, const Punto3D& pb, int );
    Linea2D( const Linea2D &l );
    virtual ~Linea2D() {}

	// informazioni sulla linea
    int InfoLinea();
    int InSide( const Punto2D& p ) const;
    int InLinea( const Linea2D &l ) const;
    int InLineaFine( Linea2D &l ) const;
    int IsPuntoRetta( const Punto2D& p ) const;
    int IsPuntoLinea( const Punto2D& p ) const;
    int GetTipo() const { return (int)Tipo; }
    double GetCoeff() const { return coeff; }
    int HalfPlane( const Punto2D &p ) const;
    int IsNull() const { return Tipo == 3; }
    double GetAngolo() const;
	// distanza
    double GetDistance( const Punto2D& p ) const;
    double GetDistance( const Linea2D& l ) const;
	// intersezione
    int IntersecaRette( const Linea2D &l, Punto2D &intr ) const;
    int IntersecaLinee( const Linea2D &l, Punto2D &intr ) const;
    int IntersecaPerpendicolare( Punto2D &p, Punto2D &intr ) const;
	// perpendicolari/parallele
    Linea2D CostruisciPerpendicolare( int ptflag ) const;
	Linea2D ParallelaPerPunto(Punto2D& pto) const;
	Linea2D PerpendicolarePerPunto(Punto2D& pto) const;
	// operatori
    Linea2D& operator = ( const Linea2D& l );
    Punto2D operator[] ( int idx ) const;
    int  operator == ( const Linea2D& l2 ) const;
    int  operator != ( const Linea2D& l2 ) const;
    Linea2D operator & ( const Linea2D& l2 ) const;
	int operator | ( const Linea2D& l2 ) const;
	int operator || ( const Linea2D& l2 ) const;
	// epsilon
    static void SetEpsilonAng( float eps );
    static void GetEpsilonAng( float& );

protected:
    void InitLinea();
};

class Linea3D
{
protected:
    Punto3D p1, p2;

public:
    Linea3D() {}
    Linea3D( Punto3D& pa, Punto3D& pb ) { p1=pa; p2=pb; }
    Linea3D( Linea3D &l ) { p1 = l.p1; p2 = l.p2; }
    virtual ~Linea3D() {}
    double GetDistance( const Punto3D& p ) const;
    Linea3D& operator = ( const Linea3D& l ) { p1 = l.p1; p2 = l.p2; return *this; }
    int  operator == ( const Linea3D& l2 ) const;
    int  operator != ( const Linea3D& l2 ) const;
    int  operator & ( const Linea3D& l2 ) const;

protected:
    void InitLinea();
};

typedef TArray ArrayMatrix3;
typedef TArray ArrayPunto2D;
typedef TArray ArrayPunto3D;
typedef TArray ArrayLinea2D;
typedef TArray ArrayLinea3D;

extern "C" double CoeffAng( double xa, double ya, double xb, double yb );

int arcpointlist( Punto2D centro,
	double start, double stop, float raggio, ArrayPunto2D &punti );
int arcpointlist( Punto2D centro,
	double start, double stop, float raggio, ArrayPunto3D &punti );
int ellipseinlist( Punto2D centro, double a, double b, ArrayPunto2D &punti );
int ellipseinlist( Punto2D centro, double a, double b, ArrayPunto3D &punti );

// per programmi che non sono il C++
extern "C"
{
  int arcpointlist( float x, float y,
			double start, double stop, float raggio,
			float *xlist, float *ylist, int maxpunti );
  int ellipseinlist( float x, float y,
			double a, double b,
			float *xlist, float *ylist, int maxpunti );
}

int ArcoPerTrePunti( Punto2D estremo1, Punto2D estremo2, Punto2D p,	ArcInfo &ai );

int Cerchio3Punti( Punto2D estremo1, Punto2D estremo2,
			Punto2D PtMedio, Punto2D& Centro, double& raggio );
int Arco3Punti( Punto2D estremo1, Punto2D estremo2, Punto2D PtMedio,
			Punto2D& Centro, double& raggio, double& AngStart, double& AngStop );

// per programmi che non sono il C++
extern "C"
{
  int ArcoPerTrePunti(
		float p1x, float p1y,
		float p2x, float p2y,
		float p3x, float p3y,
		ArcInfo *ai );
}

#else  // __cplusplus

// eventuali programmi solo C non C++ includono quest sezione

double CoeffAng( double xa, double ya, double xb, double yb );

// funzioni di serializzazione arco in punti
int arcpointlist( float x, float y,
	double start, double stop, float raggio,
	float *xlist, float *ylist, int maxpunti );
int ellipseinlist( float x, float y,
	double a, double b,
	float *xlist, float *ylist, int maxpunti );

// calcolo di un arco che passi per tre punti dati
int ArcoPerTrePunti(
		float p1x, float p1y,
		float p2x, float p2y,
		float p3x, float p3y,
		float* pcx, float* pcy,
		float* raggio
		);

#endif // __cplusplus

#endif //  NMATH__