#include <math.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <stdarg.h>
#include <limits.h>
#include <assert.h>
#include "iqm.h"

#define ASSERT(c) if(c) {}

#ifdef NULL
#undef NULL
#endif
#define NULL 0

#ifdef _WIN32
#ifndef M_PI
#define M_PI 3.1415926535897932384626433832795
#endif
#define strcasecmp _stricmp
#define strncasecmp _strnicmp
#endif

typedef unsigned char uchar;
typedef unsigned short ushort;
typedef unsigned int uint;
typedef signed long long int llong;
typedef unsigned long long int ullong;

inline void *operator new(size_t size)
{
	void *p = malloc(size);
	if(!p) abort();
	return p;
}
inline void *operator new[](size_t size)
{
	void *p = malloc(size);
	if(!p) abort();
	return p;
}
inline void operator delete(void *p) { if(p) free(p); }
inline void operator delete[](void *p) { if(p) free(p); }

inline void *operator new(size_t, void *p) { return p; }
inline void *operator new[](size_t, void *p) { return p; }
inline void operator delete(void *, void *) {}
inline void operator delete[](void *, void *) {}

#ifdef swap
#undef swap
#endif
template<class T>
static inline void swap(T &a, T &b)
{
	T t = a;
	a = b;
	b = t;
}
#ifdef max
#undef max
#endif
#ifdef min
#undef min
#endif
template<class T>
static inline T max(T a, T b)
{
	return a > b ? a : b;
}
template<class T>
static inline T min(T a, T b)
{
	return a < b ? a : b;
}

#define clamp(a,b,c) (max(b, min(a, c)))

#define loop(v,m) for(int v = 0; v<int(m); v++)
#define loopi(m) loop(i,m)
#define loopj(m) loop(j,m)
#define loopk(m) loop(k,m)
#define loopl(m) loop(l,m)

#ifdef WIN32
#ifdef M_PI
#undef M_PI
#endif
#define M_PI 3.14159265

#ifndef __GNUC__
#pragma warning (3: 4189)       // local variable is initialized but not referenced
#pragma warning (disable: 4244) // conversion from 'int' to 'float', possible loss of data
#pragma warning (disable: 4267) // conversion from 'size_t' to 'int', possible loss of data
#pragma warning (disable: 4355) // 'this' : used in base member initializer list
#pragma warning (disable: 4996) // 'strncpy' was declared deprecated
#endif

#define strcasecmp _stricmp
#define PATHDIV '\\'
#else
#define __cdecl
#define _vsnprintf vsnprintf
#define PATHDIV '/'
#endif

// easy safe strings

#define MAXSTRLEN 260
typedef char string[MAXSTRLEN];

inline void vformatstring(char *d, const char *fmt, va_list v, int len = MAXSTRLEN) { _vsnprintf(d, len, fmt, v); d[len-1] = 0; }
inline char *copystring(char *d, const char *s, size_t len = MAXSTRLEN)
{ 
	size_t slen = min(strlen(s)+1, len);
	memcpy(d, s, slen);
	d[slen-1] = 0;
	return d;
}
inline char *concatstring(char *d, const char *s) { size_t len = strlen(d); return copystring(d+len, s, MAXSTRLEN-len); }

template<size_t N> inline void formatstring(char (&d)[N], const char *fmt, ...)
{
	va_list v;
	va_start(v, fmt);
	vformatstring(d, fmt, v, int(N));
	va_end(v);
}

#define defformatstring(d,...) string d; formatstring(d, __VA_ARGS__)
#define defvformatstring(d,last,fmt) string d; { va_list ap; va_start(ap, last); vformatstring(d, fmt, ap); va_end(ap); }

inline char *newstring(size_t l)                { return new char[l+1]; }
inline char *newstring(const char *s, size_t l) { return copystring(newstring(l), s, l+1); }
inline char *newstring(const char *s)           { size_t l = strlen(s); char *d = newstring(l); memcpy(d, s, l+1); return d; }

#define loopv(v)    for(int i = 0; i<(v).length(); i++)
#define loopvj(v)   for(int j = 0; j<(v).length(); j++)
#define loopvk(v)   for(int k = 0; k<(v).length(); k++)
#define loopvrev(v) for(int i = (v).length()-1; i>=0; i--)

template <class T> struct vector
{
	static const int MINSIZE = 8;

	T *buf;
	int alen, ulen;

	vector() : buf(NULL), alen(0), ulen(0)
	{
	}

	vector(const vector &v) : buf(NULL), alen(0), ulen(0)
	{
		*this = v;
	}

	~vector() { setsize(0); if(buf) delete[] (uchar *)buf; }

	vector<T> &operator=(const vector<T> &v)
	{
		setsize(0);
		if(v.length() > alen) growbuf(v.length());
		loopv(v) add(v[i]);
		return *this;
	}

	T &add(const T &x)
	{
		if(ulen==alen) growbuf(ulen+1);
		new (&buf[ulen]) T(x);
		return buf[ulen++];
	}

	T &add()
	{
		if(ulen==alen) growbuf(ulen+1);
		new (&buf[ulen]) T;
		return buf[ulen++];
	}

	T &dup()
	{
		if(ulen==alen) growbuf(ulen+1);
		new (&buf[ulen]) T(buf[ulen-1]);
		return buf[ulen++];
	}

	bool inrange(uint i) const { return i<uint(ulen); }
	bool inrange(int i) const { return i>=0 && i<ulen; }

	T &pop() { return buf[--ulen]; }
	T &last() { return buf[ulen-1]; }
	void drop() { ulen--; buf[ulen].~T(); }

	bool empty() const { return ulen==0; }
	int capacity() const { return alen; }
	int length() const { return ulen; }
	T &operator[](int i) { ASSERT(i>=0 && i<ulen); return buf[i]; }
	const T &operator[](int i) const { ASSERT(i >= 0 && i<ulen); return buf[i]; }

	void setsize(int i) { ASSERT(i <= ulen); ulen = i; }

	void swap(vector<T> &v)
	{
		::swap(buf, v.buf);
		::swap(ulen, v.ulen);
		::swap(alen, v.alen);
	}

	T *getbuf() { return buf; }
	const T *getbuf() const { return buf; }
	bool inbuf(const T *e) const { return e >= buf && e < &buf[ulen]; }

	void growbuf(int sz)
	{
		int olen = alen;
		if(!alen) alen = max(MINSIZE, sz);
		else while(alen < sz) alen *= 2;
		if(alen <= olen) return;
		uchar *newbuf = new uchar[alen*sizeof(T)];
		if(olen > 0)
		{
			memcpy(newbuf, buf, olen*sizeof(T));
			delete[] (uchar *)buf;
		}
		buf = (T *)newbuf;
	}

	T *reserve(int sz)
	{
		if(ulen+sz > alen) growbuf(ulen+sz);
		return &buf[ulen];
	}

	void advance(int sz)
	{
		ulen += sz;
	}

	void put(const T *v, int n)
	{
		memcpy(reserve(n), v, n*sizeof(T));
		advance(n);
	}
};

static inline uint hthash(const char *key)
{
	uint h = 5381;
	for(int i = 0, k; (k = key[i]); i++) h = ((h<<5)+h)^k;    // bernstein k=33 xor
	return h;
}

static inline bool htcmp(const char *x, const char *y)
{
	return !strcmp(x, y);
}

static inline uint hthash(int key)
{
	return key;
}

static inline bool htcmp(int x, int y)
{
	return x==y;
}

static inline bool htcmp(double x, double y)
{
	return x == y;
}

static inline uint hthash(double k)
{
	union { double f; uint h[sizeof(double)/sizeof(uint)]; } conv;
	conv.f = k;
	uint hash = conv.h[0];
	for(size_t i = 1; i < sizeof(conv.h)/sizeof(uint); i++) hash ^= conv.h[i];
	return hash;
}

template <class K, class T> struct hashtable
{
	typedef K key;
	typedef const K const_key;
	typedef T value;
	typedef const T const_value;

	enum { CHUNKSIZE = 64 };

	struct chain      { T data; K key; chain *next; };
	struct chainchunk { chain chains[CHUNKSIZE]; chainchunk *next; };

	int size;
	int numelems;
	chain **table;

	chainchunk *chunks;
	chain *unused;

	hashtable(int size = 1<<10)
	  : size(size)
	{
		numelems = 0;
		chunks = NULL;
		unused = NULL;
		table = new chain *[size];
		loopi(size) table[i] = NULL;
	}

	~hashtable()
	{
		if(table) delete[] table;
		deletechunks();
	}

	chain *insert(const K &key, uint h)
	{
		if(!unused)
		{
			chainchunk *chunk = new chainchunk;
			chunk->next = chunks;
			chunks = chunk;
			loopi(CHUNKSIZE-1) chunk->chains[i].next = &chunk->chains[i+1];
			chunk->chains[CHUNKSIZE-1].next = unused;
			unused = chunk->chains;
		}
		chain *c = unused;
		unused = unused->next;
		c->key = key;
		c->next = table[h]; 
		table[h] = c;
		numelems++;
		return c;
	}

	#define HTFIND(success, fail) \
		uint h = hthash(key)&(size-1); \
		for(chain *c = table[h]; c; c = c->next) \
		{ \
			if(htcmp(key, c->key)) return (success); \
		} \
		return (fail);

	template<class L>
	T *access(const L &key)
	{
		HTFIND(&c->data, NULL);
	}

	template<class L>
	T &access(const L &key, const T &data)
	{
		HTFIND(c->data, insert(key, h)->data = data);
	}

	template<class L>
	const T &find(const L &key, const T &notfound)
	{
		HTFIND(c->data, notfound);
	}

	template<class L>
	T &operator[](const L &key)
	{
		HTFIND(c->data, insert(key, h)->data);
	}

	#undef HTFIND
   
	template<class L>
	bool remove(const L &key)
	{
		uint h = hthash(key)&(size-1); 
		for(chain **p = &table[h], *c = table[h]; c; p = &c->next, c = c->next)
		{
			if(htcmp(key, c->key))
			{
				*p = c->next;
				c->data.~T();
				c->key.~K();
				new (&c->data) T;
				new (&c->key) K;
				c->next = unused;
				unused = c;
				numelems--;
				return true;
			}
		}
		return false;
	}

	void deletechunks()
	{
		for(chainchunk *nextchunk; chunks; chunks = nextchunk)
		{
			nextchunk = chunks->next;
			delete chunks;
		}
	}

	void clear()
	{
		if(!numelems) return;
		loopi(size) table[i] = NULL;
		numelems = 0;
		unused = NULL;
		deletechunks();
	}
};

#define enumerate(ht,k,e,t,f,b) loopi((ht).size)  for(hashtable<k,t>::chain *enumc = (ht).table[i]; enumc;) { hashtable<k,t>::const_key &e = enumc->key; t &f = enumc->data; enumc = enumc->next; b; }

template<class T>
struct unionfind
{
	struct ufval
	{
		int rank, next;
		T val;

		ufval(const T &val) : rank(0), next(-1), val(val) {}
	};

	vector<ufval> ufvals;

	void clear()
	{
		ufvals.setsize(0);
	}

	const T &find(int k, const T &noval, const T &initval)
	{
		if(k>=ufvals.length()) return initval;
		while(ufvals[k].next>=0) k = ufvals[k].next;
		if(ufvals[k].val == noval) ufvals[k].val = initval;
		return ufvals[k].val;
	}

	int compressfind(int k)
	{
		if(ufvals[k].next<0) return k;
		return ufvals[k].next = compressfind(ufvals[k].next);
	}

	void unite (int x, int y, const T &noval)
	{
		while(ufvals.length() <= max(x, y)) ufvals.add(ufval(noval));
		x = compressfind(x);
		y = compressfind(y);
		if(x==y) return;
		ufval &xval = ufvals[x], &yval = ufvals[y];
		if(xval.rank < yval.rank) xval.next = y;
		else
		{
			yval.next = x;
			if(xval.rank==yval.rank) yval.rank++;
		}
	}
};

template<class T>
struct listnode
{
	T *prev, *next;
};

template<class T>
struct list
{
	typedef listnode<T> node;

	int size;
	listnode<T> nodes;

	list() { clear(); }

	bool empty() const { return nodes.prev == nodes.next; }
	bool notempty() const { return nodes.prev != nodes.next; }

	T *first() const { return nodes.next; }
	T *last() const { return nodes.prev; }
	T *end() const { return (T *)&nodes; }

	void clear()
	{
		size = 0;
		nodes.prev = nodes.next = (T *)&nodes;
	}

	T *remove(T *node)
	{
		size--;
		node->prev->next = node->next;
		node->next->prev = node->prev;
		return node;
	}

	T *insertafter(T *node, T *pos)
	{
		size++;
		node->next = pos->next;
		node->next->prev = node;
		node->prev = pos;
		pos->next = node;
		return node;
	}

	T *insertbefore(T *node, T *pos)
	{
		size++;
		node->prev = pos->prev;
		node->prev->next = node;
		node->next = pos;
		pos->prev = node;
		return node;
	} 

	T *insertfirst(T *node) { return insertafter(node, end()); }
	T *insertlast(T *node) { return insertbefore(node, end()); }

	T *removefirst() { return remove(first()); }
	T *removelast() { return remove(last()); }
};

static inline bool islittleendian() { union { int i; uchar b[sizeof(int)]; } conv; conv.i = 1; return conv.b[0] != 0; }
inline ushort endianswap16(ushort n) { return (n<<8) | (n>>8); }
inline uint endianswap32(uint n) { return (n<<24) | (n>>24) | ((n>>8)&0xFF00) | ((n<<8)&0xFF0000); }
inline ullong endianswap64(ullong n) { return endianswap32(uint(n >> 32)) | ((ullong)endianswap32(uint(n)) << 32); }
template<class T> inline T endianswap(T n) { union { T t; uint i; } conv; conv.t = n; conv.i = endianswap32(conv.i); return conv.t; }
template<> inline uchar endianswap<uchar>(uchar n) { return n; }
template<> inline char endianswap<char>(char n) { return n; }
template<> inline ushort endianswap<ushort>(ushort n) { return endianswap16(n); }
template<> inline short endianswap<short>(short n) { return endianswap16(n); }
template<> inline uint endianswap<uint>(uint n) { return endianswap32(n); }
template<> inline int endianswap<int>(int n) { return endianswap32(n); }
template<> inline ullong endianswap<ullong>(ullong n) { return endianswap64(n); }
template<> inline llong endianswap<llong>(llong n) { return endianswap64(n); }
template<> inline double endianswap<double>(double n) { union { double t; uint i; } conv; conv.t = n; conv.i = endianswap64(conv.i); return conv.t; }
template<class T> inline void endianswap(T *buf, int len) { for(T *end = &buf[len]; buf < end; buf++) *buf = endianswap(*buf); }
template<class T> inline T endiansame(T n) { return n; }
template<class T> inline void endiansame(T *buf, int len) {}
template<class T> inline T lilswap(T n) { return islittleendian() ? n : endianswap(n); }
template<class T> inline void lilswap(T *buf, int len) { if(!islittleendian()) endianswap(buf, len); }
template<class T> inline T bigswap(T n) { return islittleendian() ? endianswap(n) : n; }
template<class T> inline void bigswap(T *buf, int len) { if(islittleendian()) endianswap(buf, len); }

/* workaround for some C platforms that have these two functions as macros - not used anywhere */
#ifdef getchar
#undef getchar
#endif
#ifdef putchar
#undef putchar
#endif

struct stream
{
	virtual ~stream() {}
	virtual void close() = 0;
	virtual bool end() = 0;
	virtual long tell() { return -1; }
	virtual bool seek(long offset, int whence = SEEK_SET) { return false; }
	virtual long size();
	virtual size_t read(void *buf, size_t len) { return 0; }
	virtual size_t write(const void *buf, size_t len) { return 0; }
	virtual int getchar() { uchar c; return read(&c, 1) == 1 ? c : -1; }
	virtual bool putchar(int n) { uchar c = n; return write(&c, 1) == 1; }
	virtual bool getline(char *str, size_t len);
	virtual bool putstring(const char *str) { size_t len = strlen(str); return write(str, len) == len; }
	virtual bool putline(const char *str) { return putstring(str) && putchar('\n'); }
	virtual int printf(const char *fmt, ...) { return -1; }

	template<class T> bool put(T n) { return write(&n, sizeof(n)) == sizeof(n); }
	template<class T> bool putlil(T n) { return put<T>(lilswap(n)); }
	template<class T> bool putbig(T n) { return put<T>(bigswap(n)); }

	template<class T> T get() { T n; return read(&n, sizeof(n)) == sizeof(n) ? n : 0; }
	template<class T> T getlil() { return lilswap(get<T>()); }
	template<class T> T getbig() { return bigswap(get<T>()); }
};

long stream::size()
{
	long pos = tell(), endpos;
	if(pos < 0 || !seek(0, SEEK_END)) return -1;
	endpos = tell();
	return pos == endpos || seek(pos, SEEK_SET) ? endpos : -1;
}

bool stream::getline(char *str, size_t len)
{
	loopi(len-1)
	{
		if(read(&str[i], 1) != 1) { str[i] = '\0'; return i > 0; }
		else if(str[i] == '\n') { str[i+1] = '\0'; return true; }
	}
	if(len > 0) str[len-1] = '\0';
	return true;
}

struct filestream : stream
{
	FILE *file;

	filestream() : file(NULL) {}
	~filestream() { close(); }

	bool open(const char *name, const char *mode)
	{
		if(file) return false;
		file = fopen(name, mode);
		return file!=NULL;
	}

	void close()
	{
		if(file) { fclose(file); file = NULL; }
	}

	bool end() { return feof(file)!=0; }
	long tell() { return ftell(file); }
	bool seek(long offset, int whence) { return fseek(file, offset, whence) >= 0; }
	size_t read(void *buf, size_t len) { return fread(buf, 1, len, file); }
	size_t write(const void *buf, size_t len) { return fwrite(buf, 1, len, file); }
	int getchar() { return fgetc(file); }
	bool putchar(int c) { return fputc(c, file)!=EOF; }
	bool getline(char *str, int len) { return fgets(str, len, file)!=NULL; }
	bool putstring(const char *str) { return fputs(str, file)!=EOF; }

	int printf(const char *fmt, ...)
	{
		va_list v;
		va_start(v, fmt);
		int result = vfprintf(file, fmt, v);
		va_end(v);
		return result;
	}
};

char *path(char *s)
{
	for(char *curpart = s;;)
	{
		char *endpart = strchr(curpart, '&');
		if(endpart) *endpart = '\0';
		if(curpart[0]=='<')
		{
			char *file = strrchr(curpart, '>');
			if(!file) return s;
			curpart = file+1;
		}
		for(char *t = curpart; (t = strpbrk(t, "/\\")); *t++ = PATHDIV);
		for(char *prevdir = NULL, *curdir = s;;)
		{
			prevdir = curdir[0]==PATHDIV ? curdir+1 : curdir;
			curdir = strchr(prevdir, PATHDIV);
			if(!curdir) break;
			if(prevdir+1==curdir && prevdir[0]=='.')
			{
				memmove(prevdir, curdir+1, strlen(curdir+1)+1);
				curdir = prevdir;
			}
			else if(curdir[1]=='.' && curdir[2]=='.' && curdir[3]==PATHDIV)
			{
				if(prevdir+2==curdir && prevdir[0]=='.' && prevdir[1]=='.') continue;
				memmove(prevdir, curdir+4, strlen(curdir+4)+1);
				curdir = prevdir;
			}
		}
		if(endpart)
		{
			*endpart = '&';
			curpart = endpart+1;
		}
		else break;
	}
	return s;
}

char *path(const char *s, bool copy)
{
	static string tmp;
	copystring(tmp, s);
	path(tmp);
	return tmp;
}

const char *parentdir(const char *directory)
{
	const char *p = directory + strlen(directory);
	while(p > directory && *p != '/' && *p != '\\') p--;
	static string parent;
	size_t len = p-directory+1;
	copystring(parent, directory, len);
	return parent;
}

stream *openfile(const char *filename, const char *mode)
{
	filestream *file = new filestream;
	if(!file->open(path(filename, true), mode)) { delete file; return NULL; }
	return file;
}

struct Vec4;

struct Vec3
{
	union
	{
		struct { double x, y, z; };
		double v[3];
		uint h[3*sizeof(double)/sizeof(uint)];
	};

	Vec3() {}
	Vec3(double x, double y, double z) : x(x), y(y), z(z) {}
	explicit Vec3(const double *v) : x(v[0]), y(v[1]), z(v[2]) {}
	explicit Vec3(const Vec4 &v);

	double &operator[](int i) { return v[i]; }
	double operator[](int i) const { return v[i]; }

	bool operator==(const Vec3 &o) const { return x == o.x && y == o.y && z == o.z; }
	bool operator!=(const Vec3 &o) const { return x != o.x || y != o.y || z != o.z; }
	bool operator<(const Vec3 &o) const { return x < o.x || y < o.y || z < o.z; }
	bool operator>(const Vec3 &o) const { return x > o.x || y > o.y || z > o.z; }

	Vec3 operator+(const Vec3 &o) const { return Vec3(x+o.x, y+o.y, z+o.z); }
	Vec3 operator-(const Vec3 &o) const { return Vec3(x-o.x, y-o.y, z-o.z); }
	Vec3 operator+(double k) const { return Vec3(x+k, y+k, z+k); }
	Vec3 operator-(double k) const { return Vec3(x-k, y-k, z-k); }
	Vec3 operator-() const { return Vec3(-x, -y, -z); }
	Vec3 operator*(const Vec3 &o) const { return Vec3(x*o.x, y*o.y, z*o.z); }
	Vec3 operator/(const Vec3 &o) const { return Vec3(x/o.x, y/o.y, z/o.z); }
	Vec3 operator*(double k) const { return Vec3(x*k, y*k, z*k); }
	Vec3 operator/(double k) const { return Vec3(x/k, y/k, z/k); }

	Vec3 &operator+=(const Vec3 &o) { x += o.x; y += o.y; z += o.z; return *this; }
	Vec3 &operator-=(const Vec3 &o) { x -= o.x; y -= o.y; z -= o.z; return *this; }
	Vec3 &operator+=(double k) { x += k; y += k; z += k; return *this; }
	Vec3 &operator-=(double k) { x -= k; y -= k; z -= k; return *this; }
	Vec3 &operator*=(const Vec3 &o) { x *= o.x; y *= o.y; z *= o.z; return *this; }
	Vec3 &operator/=(const Vec3 &o) { x /= o.x; y /= o.y; z /= o.z; return *this; }
	Vec3 &operator*=(double k) { x *= k; y *= k; z *= k; return *this; }
	Vec3 &operator/=(double k) { x /= k; y /= k; z /= k; return *this; }

	double dot(const Vec3 &o) const { return x*o.x + y*o.y + z*o.z; }
	double magnitude() const { return sqrt(dot(*this)); }
	double squaredlen() const { return dot(*this); }
	double dist(const Vec3 &o) const { return (*this - o).magnitude(); }
	Vec3 normalize() const { return *this * (1.0 / magnitude()); }
	Vec3 cross(const Vec3 &o) const { return Vec3(y*o.z-z*o.y, z*o.x-x*o.z, x*o.y-y*o.x); }
	Vec3 reflect(const Vec3 &n) const { return *this - n*2.0*dot(n); }
	Vec3 project(const Vec3 &n) const { return *this - n*dot(n); }

	Vec3 zxy() const { return Vec3(z, x, y); }
	Vec3 zyx() const { return Vec3(z, y, x); }
	Vec3 yxz() const { return Vec3(y, x, z); }
	Vec3 yzx() const { return Vec3(y, z, x); }
	Vec3 xzy() const { return Vec3(x, z, y); }
};

static inline bool htcmp(const Vec3 &x, const Vec3 &y)
{
	return x == y;
}

static inline uint hthash(const Vec3 &k)
{
	uint hash = k.h[0];
	for(size_t i = 1; i < sizeof(k.h)/sizeof(uint); i++) hash ^= k.h[i];
	return hash;
}

struct Vec4
{
	union
	{
		struct { double x, y, z, w; };
		double v[4];
		uint h[4*sizeof(double)/sizeof(uint)];
	};

	Vec4() {}
	Vec4(double x, double y, double z, double w) : x(x), y(y), z(z), w(w) {}
	explicit Vec4(const Vec3 &p, double w = 0) : x(p.x), y(p.y), z(p.z), w(w) {}
	explicit Vec4(const double *v) : x(v[0]), y(v[1]), z(v[2]), w(v[3]) {}

	double &operator[](int i)       { return v[i]; }
	double  operator[](int i) const { return v[i]; }

	bool operator==(const Vec4 &o) const { return x == o.x && y == o.y && z == o.z && w == o.w; }
	bool operator!=(const Vec4 &o) const { return x != o.x || y != o.y || z != o.z || w != o.w; }
	bool operator<(const Vec4 &o) const { return x < o.x || y < o.y || z < o.z || w < o.w; }
	bool operator>(const Vec4 &o) const { return x > o.x || y > o.y || z > o.z || w > o.w; }

	Vec4 operator+(const Vec4 &o) const { return Vec4(x+o.x, y+o.y, z+o.z, w+o.w); }
	Vec4 operator-(const Vec4 &o) const { return Vec4(x-o.x, y-o.y, z-o.z, w-o.w); }
	Vec4 operator+(double k) const { return Vec4(x+k, y+k, z+k, w+k); }
	Vec4 operator-(double k) const { return Vec4(x-k, y-k, z-k, w-k); }
	Vec4 operator-() const { return Vec4(-x, -y, -z, -w); }
	Vec4 operator*(double k) const { return Vec4(x*k, y*k, z*k, w*k); }
	Vec4 operator/(double k) const { return Vec4(x/k, y/k, z/k, w/k); }
	Vec4 addw(double f) const { return Vec4(x, y, z, w + f); }

	Vec4 &operator+=(const Vec4 &o) { x += o.x; y += o.y; z += o.z; w += o.w; return *this; }
	Vec4 &operator+=(const Vec3 &o) { x += o.x; y += o.y; z += o.z; return * this; }
	Vec4 &operator-=(const Vec4 &o) { x -= o.x; y -= o.y; z -= o.z; w -= o.w; return *this; }
	Vec4 &operator-=(const Vec3 &o) { x -= o.x; y -= o.y; z -= o.z; return * this; }
	Vec4 &operator+=(double k) { x += k; y += k; z += k; w += k; return *this; }
	Vec4 &operator-=(double k) { x -= k; y -= k; z -= k; w -= k; return *this; }
	Vec4 &operator*=(double k) { x *= k; y *= k; z *= k; w *= k; return *this; }
	Vec4 &operator/=(double k) { x /= k; y /= k; z /= k; w /= k; return *this; }

	double dot3(const Vec4 &o) const { return x*o.x + y*o.y + z*o.z; }
	double dot3(const Vec3 &o) const { return x*o.x + y*o.y + z*o.z; }
	double dot(const Vec4 &o) const { return dot3(o) + w*o.w; }
	double dot(const Vec3 &o) const  { return x*o.x + y*o.y + z*o.z + w; }
	double magnitude() const  { return sqrt(dot(*this)); }
	double magnitude3() const { return sqrt(dot3(*this)); }
	Vec4 normalize() const { return *this * (1.0 / magnitude()); }
	Vec3 cross3(const Vec4 &o) const { return Vec3(y*o.z-z*o.y, z*o.x-x*o.z, x*o.y-y*o.x); }
	Vec3 cross3(const Vec3 &o) const { return Vec3(y*o.z-z*o.y, z*o.x-x*o.z, x*o.y-y*o.x); }

	void setxyz(const Vec3 &o) { x = o.x; y = o.y; z = o.z; }
};

inline Vec3::Vec3(const Vec4 &v) : x(v.x), y(v.y), z(v.z) {}

static inline bool htcmp(const Vec4 &x, const Vec4 &y)
{
	return x == y;
}

static inline uint hthash(const Vec4 &k)
{
	uint hash = k.h[0];
	for(size_t i = 1; i < sizeof(k.h)/sizeof(uint); i++) hash ^= k.h[i];
	return hash;
}

struct Matrix3x3;
struct Matrix3x4;

struct Quat : Vec4
{
	Quat() {}
	Quat(double x, double y, double z, double w) : Vec4(x, y, z, w) {}
	Quat(double angle, const Vec3 &axis)
	{
		double s = sin(0.5*angle);
		x = s*axis.x;
		y = s*axis.y;
		z = s*axis.z;
		w = cos(0.5*angle);
	}
	explicit Quat(const Vec3 &v) : Vec4(v.x, v.y, v.z, -sqrt(max(1.0 - v.squaredlen(), 0.0))) {}
	explicit Quat(const Matrix3x3 &m) { convertmatrix(m); }
	explicit Quat(const Matrix3x4 &m) { convertmatrix(m); }

	void restorew()
	{
		w = -sqrt(max(1.0 - dot3(*this), 0.0));
	}

	Quat operator*(const Quat &o) const
	{
		return Quat(w*o.x + x*o.w + y*o.z - z*o.y,
					w*o.y - x*o.z + y*o.w + z*o.x,
					w*o.z + x*o.y - y*o.x + z*o.w,
					w*o.w - x*o.x - y*o.y - z*o.z);
	}
	Quat &operator*=(const Quat &o) { return (*this = *this * o); }

	Quat operator+(const Vec4 &o) const { return Quat(x+o.x, y+o.y, z+o.z, w+o.w); }
	Quat &operator+=(const Vec4 &o) { return (*this = *this + o); }
	Quat operator-(const Vec4 &o) const { return Quat(x-o.x, y-o.y, z-o.z, w-o.w); }
	Quat &operator-=(const Vec4 &o) { return (*this = *this - o); }

	Quat operator-() const { return Quat(-x, -y, -z, w); }

	void flip() { x = -x; y = -y; z = -z; w = -w; }

	Vec3 transform(const Vec3 &p) const
	{
		return p + cross3(cross3(p) + p*w)*2.0;
	}

	template<class M>
	void convertmatrix(const M &m)
	{
		double trace = m.a.x + m.b.y + m.c.z;
		if(trace>0)
		{
			double r = sqrt(1 + trace), inv = 0.5/r;
			w = 0.5*r;
			x = (m.c.y - m.b.z)*inv;
			y = (m.a.z - m.c.x)*inv;
			z = (m.b.x - m.a.y)*inv;
		}
		else if(m.a.x > m.b.y && m.a.x > m.c.z)
		{
			double r = sqrt(1 + m.a.x - m.b.y - m.c.z), inv = 0.5/r;
			x = 0.5*r;
			y = (m.b.x + m.a.y)*inv;
			z = (m.a.z + m.c.x)*inv;
			w = (m.c.y - m.b.z)*inv;
		}
		else if(m.b.y > m.c.z)
		{
			double r = sqrt(1 + m.b.y - m.a.x - m.c.z), inv = 0.5/r;
			x = (m.b.x + m.a.y)*inv;
			y = 0.5*r;
			z = (m.c.y + m.b.z)*inv;
			w = (m.a.z - m.c.x)*inv;
		}
		else
		{
			double r = sqrt(1 + m.c.z - m.a.x - m.b.y), inv = 0.5/r;
			x = (m.a.z + m.c.x)*inv;
			y = (m.c.y + m.b.z)*inv;
			z = 0.5*r;
			w = (m.b.x - m.a.y)*inv;
		}
	}

	static Quat fromangles(const Vec3 &rot)
	{
		double cx = cos(rot.x/2), sx = sin(rot.x/2),
			   cy = cos(rot.y/2), sy = sin(rot.y/2),
			   cz = cos(rot.z/2), sz = sin(rot.z/2);
		Quat q(sx*cy*cz - cx*sy*sz,
			   cx*sy*cz + sx*cy*sz,
			   cx*cy*sz - sx*sy*cz,
			   cx*cy*cz + sx*sy*sz);
		if(q.w > 0) q.flip();
		return q;
	}

	static Quat fromdegrees(const Vec3 &rot) { return fromangles(rot * (M_PI / 180)); }
};

struct Matrix3x3
{
	Vec3 a, b, c;

	Matrix3x3() {}
	Matrix3x3(const Vec3 &a, const Vec3 &b, const Vec3 &c) : a(a), b(b), c(c) {}
	explicit Matrix3x3(const Quat &q) { convertquat(q); }
	explicit Matrix3x3(const Quat &q, const Vec3 &scale)
	{
		convertquat(q);
		a *= scale;
		b *= scale;
		c *= scale;
	}

	void convertquat(const Quat &q)
	{
		double x = q.x, y = q.y, z = q.z, w = q.w,
			   tx = 2*x, ty = 2*y, tz = 2*z,
			   txx = tx*x, tyy = ty*y, tzz = tz*z,
			   txy = tx*y, txz = tx*z, tyz = ty*z,
			   twx = w*tx, twy = w*ty, twz = w*tz;
		a = Vec3(1 - (tyy + tzz), txy - twz, txz + twy);
		b = Vec3(txy + twz, 1 - (txx + tzz), tyz - twx);
		c = Vec3(txz - twy, tyz + twx, 1 - (txx + tyy));
	}

	Matrix3x3 operator*(const Matrix3x3 &o) const
	{
		return Matrix3x3(
			o.a*a.x + o.b*a.y + o.c*a.z,
			o.a*b.x + o.b*b.y + o.c*b.z,
			o.a*c.x + o.b*c.y + o.c*c.z);
	}
	Matrix3x3 &operator*=(const Matrix3x3 &o) { return (*this = *this * o); }

	void transpose(const Matrix3x3 &o)
	{
		a = Vec3(o.a.x, o.b.x, o.c.x);
		b = Vec3(o.a.y, o.b.y, o.c.y);
		c = Vec3(o.a.z, o.b.z, o.c.z);
	}
	void transpose() { transpose(Matrix3x3(*this)); }

	Vec3 transform(const Vec3 &o) const { return Vec3(a.dot(o), b.dot(o), c.dot(o)); }

	float determinant()
	{
		return
			a.x * b.y * c.z +
			a.y * b.z * c.x +
			a.z * b.x * c.y -
			a.z * b.y * c.x -
			a.y * b.x * c.z -
			a.x * b.z * c.y;
	}
};

struct Matrix3x4
{
	Vec4 a, b, c;

	Matrix3x4() {}
	Matrix3x4(const Vec4 &a, const Vec4 &b, const Vec4 &c) : a(a), b(b), c(c) {}
	explicit Matrix3x4(const Matrix3x3 &rot, const Vec3 &trans)
		: a(Vec4(rot.a, trans.x)), b(Vec4(rot.b, trans.y)), c(Vec4(rot.c, trans.z))
	{
	}
	explicit Matrix3x4(const Quat &rot, const Vec3 &trans)
	{
		*this = Matrix3x4(Matrix3x3(rot), trans);
	}
	explicit Matrix3x4(const Quat &rot, const Vec3 &trans, const Vec3 &scale)
	{
		*this = Matrix3x4(Matrix3x3(rot, scale), trans);
	}

	Matrix3x4 operator*(float k) const { return Matrix3x4(*this) *= k; }
	Matrix3x4 &operator*=(float k)
	{
		a *= k;
		b *= k;
		c *= k;
		return *this;
	}

	Matrix3x4 operator+(const Matrix3x4 &o) const { return Matrix3x4(*this) += o; }
	Matrix3x4 &operator+=(const Matrix3x4 &o)
	{
		a += o.a;
		b += o.b;
		c += o.c;
		return *this;
	}
	Matrix3x4 operator+(const Vec3 &o) const { return Matrix3x4(*this) += o; }
	Matrix3x4 &operator+=(const Vec3 &o)
	{
		a[3] += o[0];
		b[3] += o[1];
		c[3] += o[2];
		return *this;
	}

	void invert(const Matrix3x4 &o)
	{
		Matrix3x3 invrot(Vec3(o.a.x, o.b.x, o.c.x), Vec3(o.a.y, o.b.y, o.c.y), Vec3(o.a.z, o.b.z, o.c.z));
		invrot.a /= invrot.a.squaredlen();
		invrot.b /= invrot.b.squaredlen();
		invrot.c /= invrot.c.squaredlen();
		Vec3 trans(o.a.w, o.b.w, o.c.w);
		a = Vec4(invrot.a, -invrot.a.dot(trans));
		b = Vec4(invrot.b, -invrot.b.dot(trans));
		c = Vec4(invrot.c, -invrot.c.dot(trans));
	}
	void invert() { invert(Matrix3x4(*this)); }

	Matrix3x4 operator*(const Matrix3x4 &o) const
	{
		return Matrix3x4(
			(o.a*a.x + o.b*a.y + o.c*a.z).addw(a.w),
			(o.a*b.x + o.b*b.y + o.c*b.z).addw(b.w),
			(o.a*c.x + o.b*c.y + o.c*c.z).addw(c.w));
	}
	Matrix3x4 &operator*=(const Matrix3x4 &o) { return (*this = *this * o); }

	Vec3 transform(const Vec3 &o) const { return Vec3(a.dot(o), b.dot(o), c.dot(o)); }
};

void conoutf(const char *s, ...)
{
	defvformatstring(msg,s,s);
	printf("%s\n", msg);
}

void fatal(const char *s, ...)    // failure exit
{
	defvformatstring(msg,s,s);
	fprintf(stderr, "%s\n", msg);

	exit(EXIT_FAILURE);
}