Go: package main

import (
        "fmt"
)

func main() {
        var x1, y1, x2, y2, x3, y3, px, py float64
        var oax, oay, obx, oby, opx, opy float64
        var det float64
        var s, t float64

        for {
                _, err := fmt.Scanf("%f %f %f %f %f %f %f %f", &x1, &y1, &x2, &y2, &x3, &y3, &px, &py)
                if err != nil {
                        break
                }
                oax = x2 - x1
                oay = y2 - y1
                obx = x3 - x1
                oby = y3 - y1
                opx = px - x1
                opy = py - y1
                det = oax*oby - obx*oay

                if det == 0 {
                        fmt.Printf("NO\n")
                } else {
                        s = (oby*opx - obx*opy) / det
                        t = (oax*opy - oay*opx) / det

                        if (s+t < 1 && s > 0 && t > 0) {
                                fmt.Printf("YES\n")
                        } else {
                                fmt.Printf("NO\n")
                        }
                }
        }
}

C++: #include <bits/stdc++.h>

using namespace std;

#define rep(i, n) for (int i = 0; i < (int)(n); i++)
#define each(i, c) for (__typeof((c).begin()) i = (c).begin(); i != (c).end(); ++i)
#define chmin(a, b) a = min(a, b)
#define chmax(a, b) a = max(a, b)
#define pb push_back
#define mp make_pair

#define X real()
#define Y imag()

const double EPS = 1e-12;
const double INF = 1e12;
const double PI = acos(-1.0);

double SQ(double a){ return a * a; }
bool EQ(double a, double b){ return abs(a - b) < EPS; }
bool LT(double a, double b){ return a - b < -EPS; }
bool LEQ(double a, double b){ return a - b < EPS; }

double toRad(double t){ return t / 180 * PI; }
double toDeg(double t){ return t * 180 / PI; }

typedef complex<double> Point;
namespace std {
	bool operator < (const Point &a, const Point &b){
		if (a.X != b.X) return a.X < b.X;
		return a.Y < b.Y;
	}
	
	istream& operator >> (istream &is, Point &a){
		double x, y;
		is >> x >> y;
		a = Point(x, y);
		return is;
	}
	
	ostream& operator << (ostream &os, Point &a){
		return os << a.X << " " << a.Y;
	}
}

struct Line {
	Point p1, p2;
};

typedef Line Segment;

struct Circle {
	Point p;
	double r;
};

typedef vector<Point> Polygon;


double norm(Point p){ return SQ(p.X) + SQ(p.Y); }

double dot(Point a, Point b){ return (conj(a) * b).X; }
double cross(Point a, Point b){ return (conj(a) * b).Y; }

enum CCW {
	COUNTER_CLOCKWISE = 1,
	CLOCKWISE = -1,
	ONLINE_BACK = 2,
	ONLINE_FRONT = -2,
	ON_SEGMENT = 0,
};

int ccw(Point p0, Point p1, Point p2){
	Point a = p1 - p0;
	Point b = p2 - p0;
	if (cross(a, b) > EPS) return COUNTER_CLOCKWISE;
	if (cross(a, b) < -EPS) return CLOCKWISE;
	if (dot(a, b) < -EPS) return ONLINE_BACK;
	if (LT(norm(a), norm(b))) return ONLINE_FRONT;
	return ON_SEGMENT;
}

bool orthogonal(Line l1, Line l2){
	return EQ(dot(l1.p2 - l1.p1, l2.p2 - l2.p1), 0.0);
}

bool parallel(Line l1, Line l2){
	return EQ(cross(l1.p2 - l1.p1, l2.p2 - l2.p1), 0.0);
}

Point project(Line l, Point p){
	Point base = l.p2 - l.p1;
	double r = dot(p - l.p1, base) / norm(base);
	return l.p1 + base * r;
}

Point reflect(Line &l, Point &p){
	return p + (project(l, p) - p) * 2.0;
}

bool sameLine(Line l1, Line l2){
	return abs(cross(l1.p2 - l1.p1, l2.p1 - l1.p1)) < EPS;
}

bool intersectLP(Line l, Point p){
	return LEQ(abs(cross(l.p1 - p, l.p2 - p)), 0.0);
}

bool intersectLL(Line l1, Line l2){
	return !parallel(l1, l2) || sameLine(l1, l2);
}

bool intersectLS(Line l, Segment s){
	return cross(l.p2 - l.p1, s.p1 - l.p1) * cross(l.p2 - l.p1, s.p2 - l.p1) < -EPS;
}

bool intersectSP(Segment s, Point p){
	return norm(s.p1 - p) + norm(s.p2 - p) - norm(s.p2 - s.p1) < EPS;
}

bool intersectSS(Segment s1, Segment s2){
	return ccw(s1.p1, s1.p2, s2.p1) * ccw(s1.p1, s1.p2, s2.p2) <= 0
		&& ccw(s2.p1, s2.p2, s1.p1) * ccw(s2.p1, s2.p2, s1.p2) <= 0;
}

double distanceLP(Line l, Point p){
	return abs(cross(l.p2 - l.p1, p - l.p1)) / abs(l.p2 - l.p1);
}

double distanceLL(Line l1, Line l2){
	if (intersectLL(l1, l2)) return 0.0;
	return distanceLP(l1, l2.p1);
}

double distanceLS(Line l, Segment s){
	if (intersectLS(l, s)) return 0.0;
	return min(distanceLP(l, s.p1), distanceLP(l, s.p2));
}

double distanceSP(Segment s, Point p){
	if (dot(s.p2 - s.p1, p - s.p1) < 0.0) return abs(p - s.p1);
	if (dot(s.p1 - s.p2, p - s.p2) < 0.0) return abs(p - s.p2);
	return distanceLP(s, p);
}

double distanceSS(Segment s1, Segment s2){
	if (intersectSS(s1, s2)) return 0.0;
	return min(
		min(distanceSP(s1, s2.p1), distanceSP(s1, s2.p2)),
		min(distanceSP(s2, s1.p1), distanceSP(s2, s1.p2))
	);
}

Point crossPointLL(Line l1, Line l2){
	assert(intersectLL(l1, l2));
	assert(!sameLine(l1, l2));
	Point base = l2.p2 - l2.p1;
	double d1 = abs(cross(base, l1.p1 - l2.p1));
	double d2 = abs(cross(base, l1.p2 - l2.p1));
	double t = d1 / (d1 + d2);
	return l1.p1 + (l1.p2 - l1.p1) * t;
}

Point crossPointLS(Line l, Segment s){
	assert(intersectLS(l, s));
	return crossPointLL(l, s);
}

Point crossPointSS(Segment s1, Segment s2){
	assert(intersectSS(s1, s2));
	return crossPointLL(s1, s2);
}

enum {
	IN = 2,
	ON = 1,
	OUT = 0,
};

int contains(const Polygon &g, Point p){
	int n = g.size();
	bool x = false;
	rep(i, n){
		Point a = g[i] - p, b = g[(i + 1) % n] - p;
		if (abs(cross(a, b)) < EPS && dot(a, b) < EPS) return ON;
		if (a.Y > b.Y) swap(a, b);
		if (a.Y < EPS && EPS < b.Y && cross(a, b) > EPS) x = !x;
	}
	return x ? IN : OUT;
}

Polygon convexHull(Polygon &s){
	Polygon u, l;
	if (s.size() < 3) return s;
	sort(s.begin(), s.end());
	
	u.pb(s[0]); u.pb(s[1]);
	l.pb(s[s.size() - 1]); l.pb(s[s.size() - 2]);
	
	for (int i = 2; i < s.size(); i++){
		int n = u.size();
		while (n >= 2 && ccw(u[n - 2], u[n - 1], s[i]) != CLOCKWISE){
			u.pop_back();
			n--;
		}
		u.pb(s[i]);
	}
	
	for (int i = s.size() - 3; i >= 0; i--){
		int n = l.size();
		while (n >= 2 && ccw(l[n - 2], l[n - 1], s[i]) != CLOCKWISE){
			l.pop_back();
			n--;
		}
		l.pb(s[i]);
	}
	
	reverse(l.begin(), l.end());
	for (int i = u.size() - 2; i >= 1; i--) l.pb(u[i]);
	
	return l;
}

int n;
Polygon g(3);
Point p;

int main()
{
	while (true){
		rep(i, 3) if (!(cin >> g[i])) return 0;
		cin >> p;
		
		puts(contains(g, p) ? "YES" : "NO");
	}
}

C++: #include <iostream>
#include <cmath>
using namespace std;

double cross_prod(double x0, double y0, double x1, double y1) {
	return x0 * y1 - x1 * y0;
}

int main() {
	double x[4], y[4];
	while(cin >> x[0] >> y[0]) {
		for(int i = 1; i < 4; i++)
			cin >> x[i] >> y[i];
		double a, b;
		a = b = 0;
		for(int i = 0; i < 3; i++) {
			double tmp;
			tmp = cross_prod(x[i]-x[3], y[i]-y[3], x[(i+1)%3]-x[3], y[(i+1)%3]-y[3]);
			a += tmp;
			b += abs(tmp);
		}
		a = abs(a);
		if(abs(a-b) < 1e-7)
			cout << "YES" << endl;
		else
			cout << "NO" << endl;
	}
}

C++: #include <iostream>
#include <cmath>
#include <complex>
using namespace std;

typedef double D;
typedef complex<D> P;

D cross(P x, P y){return imag(conj(x)*y);}

bool triContain(P a, P b, P c, P p){
    P ab = b - a;
    P ap = p - a;
    P bc = c - b;
    P bp = p - b;
    P ca = a - c;
    P cp = p - c;
    if(cross(ab, ap) > 0 && cross(bc, bp) > 0 && cross(ca, cp) > 0)return true;
    if(cross(ab, ap) < 0 && cross(bc, bp) < 0 && cross(ca, cp) < 0)return true;
    return false;
}

int main(){

   D ax, ay, bx, by, cx, cy, px, py;
   while(cin>>ax>>ay>>bx>>by>>cx>>cy>>px>>py){
       P a(ax, ay), b(bx, by), c(cx, cy), p(px, py);
       if(triContain(a, b, c, p))cout<<"YES"<<endl;
       else cout<<"NO"<<endl;
   }

    return 0;
}