typedef struct tagPlane
{
   float a, b, c, d;
} Plane_t;

typedef float vec_t;
typedef vec_t vec3_t[ 3 ];
typedef vec_t mat4x4_t[ 16 ];

enum
{
   plane_left = 0,
   plane_right,
   plane_top,
   plane_bottom,
   plane_near,
   plane_far,
};

Plane_t m_pPlanes[ 6 ];

void NormalizePlane( Plane_t &p )
{
   float _nMag = ( float )sqrt( p.a * p.a + p.b * p.b + p.c * p.c );

   p.a /= _nMag;
   p.b /= _nMag;
   p.c /= _nMag;   
   p.d /= _nMag;
}

float PointPlaneDistance( const Plane_t &p, const vec3_t &v )
{
   return ( p.a * v[ 0 ] + p.b * v[ 1 ] + p.c * v[ 2 ] + p.d );
}

bool PointInFrustum( const vec3_t &v )
{
   float d = 0;
   
   for ( int i = 0; i < 6; i++ )
   {
      d = PointPlaneDistance( m_pPlanes[ i ], v );

      if ( d <= 0 )
         return false;
   }

   return true;
}

void CalculateFrustum()
{
   mat4x4_t m1, m2, m3;

   glGetFloatv( GL_MODELVIEW_MATRIX, m1 );
   glGetFloatv( GL_PROJECTION_MATRIX, m2 );

   m3[ 0 ]  = m1[ 0 ]  * m2[ 0 ] + m1[ 1 ]  * m2[ 4 ] + m1[ 2 ]  * m2[ 8 ]  + m1[ 3 ]  * m2[ 12 ];
   m3[ 1 ]  = m1[ 0 ]  * m2[ 1 ] + m1[ 1 ]  * m2[ 5 ] + m1[ 2 ]  * m2[ 9 ]  + m1[ 3 ]  * m2[ 13 ];
   m3[ 2 ]  = m1[ 0 ]  * m2[ 2 ] + m1[ 1 ]  * m2[ 6 ] + m1[ 2 ]  * m2[ 10 ] + m1[ 3 ]  * m2[ 14 ];
   m3[ 3 ]  = m1[ 0 ]  * m2[ 3 ] + m1[ 1 ]  * m2[ 7 ] + m1[ 2 ]  * m2[ 11 ] + m1[ 3 ]  * m2[ 15 ];

   m3[ 4 ]  = m1[ 4 ]  * m2[ 0 ] + m1[ 5 ]  * m2[ 4 ] + m1[ 6 ]  * m2[ 8 ]  + m1[ 7 ]  * m2[ 12 ];
   m3[ 5 ]  = m1[ 4 ]  * m2[ 1 ] + m1[ 5 ]  * m2[ 5 ] + m1[ 6 ]  * m2[ 9 ]  + m1[ 7 ]  * m2[ 13 ];
   m3[ 6 ]  = m1[ 4 ]  * m2[ 2 ] + m1[ 5 ]  * m2[ 6 ] + m1[ 6 ]  * m2[ 10 ] + m1[ 7 ]  * m2[ 14 ];
   m3[ 7 ]  = m1[ 4 ]  * m2[ 3 ] + m1[ 5 ]  * m2[ 7 ] + m1[ 6 ]  * m2[ 11 ] + m1[ 7 ]  * m2[ 15 ];

   m3[ 8 ]  = m1[ 8 ]  * m2[ 0 ] + m1[ 9 ]  * m2[ 4 ] + m1[ 10 ] * m2[ 8 ]  + m1[ 11 ] * m2[ 12 ];
   m3[ 9 ]  = m1[ 8 ]  * m2[ 1 ] + m1[ 9 ]  * m2[ 5 ] + m1[ 10 ] * m2[ 9 ]  + m1[ 11 ] * m2[ 13 ];
   m3[ 10 ] = m1[ 8 ]  * m2[ 2 ] + m1[ 9 ]  * m2[ 6 ] + m1[ 10 ] * m2[ 10 ] + m1[ 11 ] * m2[ 14 ];
   m3[ 11 ] = m1[ 8 ]  * m2[ 3 ] + m1[ 9 ]  * m2[ 7 ] + m1[ 10 ] * m2[ 11 ] + m1[ 11 ] * m2[ 15 ];

   m3[ 12 ] = m1[ 12 ] * m2[ 0 ] + m1[ 13 ] * m2[ 4 ] + m1[ 14 ] * m2[ 8 ]  + m1[ 15 ] * m2[ 12 ];
   m3[ 13 ] = m1[ 12 ] * m2[ 1 ] + m1[ 13 ] * m2[ 5 ] + m1[ 14 ] * m2[ 9 ]  + m1[ 15 ] * m2[ 13 ];
   m3[ 14 ] = m1[ 12 ] * m2[ 2 ] + m1[ 13 ] * m2[ 6 ] + m1[ 14 ] * m2[ 10 ] + m1[ 15 ] * m2[ 14 ];
   m3[ 15 ] = m1[ 12 ] * m2[ 3 ] + m1[ 13 ] * m2[ 7 ] + m1[ 14 ] * m2[ 11 ] + m1[ 15 ] * m2[ 15 ];

   m_pPlanes[ plane_left ].a   = m3[ 3 ]  + m3[ 0 ];
   m_pPlanes[ plane_left ].b   = m3[ 7 ]  + m3[ 4 ];
   m_pPlanes[ plane_left ].c   = m3[ 11 ] + m3[ 8 ];
   m_pPlanes[ plane_left ].d   = m3[ 15 ] + m3[ 12 ];

   m_pPlanes[ plane_right ].a  = m3[ 3 ]  - m3[ 0 ];
   m_pPlanes[ plane_right ].b  = m3[ 7 ]  - m3[ 4 ];
   m_pPlanes[ plane_right ].c  = m3[ 11 ] - m3[ 8 ];
   m_pPlanes[ plane_right ].d  = m3[ 15 ] - m3[ 12 ];   

   m_pPlanes[ plane_bottom ].a = m3[ 3 ]  + m3[ 1 ];
   m_pPlanes[ plane_bottom ].b = m3[ 7 ]  + m3[ 5 ];
   m_pPlanes[ plane_bottom ].c = m3[ 11 ] + m3[ 9 ];
   m_pPlanes[ plane_bottom ].d = m3[ 15 ] + m3[ 13 ];

   m_pPlanes[ plane_top ].a   = m3[ 3 ]  - m3[ 1 ];
   m_pPlanes[ plane_top ].b   = m3[ 7 ]  - m3[ 5 ];
   m_pPlanes[ plane_top ].c   = m3[ 11 ] - m3[ 9 ];
   m_pPlanes[ plane_top ].d   = m3[ 15 ] - m3[ 13 ];

   m_pPlanes[ plane_near ].a   = m3[ 3 ]  + m3[ 2 ];
   m_pPlanes[ plane_near ].b   = m3[ 7 ]  + m3[ 6 ];
   m_pPlanes[ plane_near ].c   = m3[ 11 ] + m3[ 10 ];
   m_pPlanes[ plane_near ].d   = m3[ 15 ] + m3[ 14 ];

   m_pPlanes[ plane_far ].a   = m3[ 3 ]  - m3[ 2 ];
   m_pPlanes[ plane_far ].b   = m3[ 7 ]  - m3[ 6 ];
   m_pPlanes[ plane_far ].c   = m3[ 11 ] - m3[ 10 ];
   m_pPlanes[ plane_far ].d   = m3[ 15 ] - m3[ 14 ];

   for ( int i = 0; i < 6; i++ )
      NormalizePlane( m_pPlanes[ i ] );

   if ( PointInFrustum( m_pPoint ) )
      printf( "point in frustum" );
   else
      printf( "poit not in frustum" );
}