using System;
using System.IO;
using System.Collections;
using System.Collections.Generic;

namespace MyGeometry
{
	public class SymMatrix3d
	{
		public double a, b, c, d, e, f;

		public SymMatrix3d()
		{
			a = b = c = d = e = f = 0;
		}
		public SymMatrix3d(Vector3d v) : this(v.x, v.y, v.z)
		{
		}
		public SymMatrix3d(double x, double y, double z)
		{
			a = x * x;
			b = x * y;
			c = x * z;
			d = y * y;
			e = y * z;
			f = z * z;
		}
		public SymMatrix3d(double a, double b, double c, double d, double e, double f)
		{
			this.a = a;
			this.b = b;
			this.c = c;
			this.d = d;
			this.e = e;
			this.f = f;
		}

		public void Add(SymMatrix3d other)
		{
			this.a += other.a;
			this.b += other.b;
			this.c += other.c;
			this.d += other.d;
			this.e += other.e;
			this.f += other.f;
		}
		public static SymMatrix3d operator +(SymMatrix3d m, SymMatrix3d other)
		{
			return new SymMatrix3d(
				m.a + other.a,
				m.b + other.b,
				m.c + other.c,
				m.d + other.d,
				m.e + other.e,
				m.f + other.f
				);
		}
		public static SymMatrix3d operator -(SymMatrix3d m, SymMatrix3d other)
		{
			return new SymMatrix3d(
				m.a - other.a,
				m.b - other.b,
				m.c - other.c,
				m.d - other.d,
				m.e - other.e,
				m.f - other.f
				);
		}
		public static SymMatrix3d operator *(SymMatrix3d m, double other)
		{
			return new SymMatrix3d(
				m.a * other,
				m.b * other,
				m.c * other,
				m.d * other,
				m.e * other,
				m.f * other
				);
		}
		public static SymMatrix3d operator /(SymMatrix3d m, double other)
		{
			return new SymMatrix3d(
				m.a / other,
				m.b / other,
				m.c / other,
				m.d / other,
				m.e / other,
				m.f / other
				);
		}
		public static Vector3d operator *(SymMatrix3d m, Vector3d other)
		{
			return new Vector3d(
				m.a * other.x + m.b * other.y + m.c * other.z,
				m.b * other.x + m.d * other.y + m.e * other.z,
				m.c * other.x + m.e * other.y + m.f * other.z);
		}
		public Vector3d GetEigenvalues()
		{
			double a11 = this.a, a12 = this.b, a22 = this.d, a13 = this.c, a23 = this.e, a33 = this.f;
			double c0 = (a11*a22*a33)+(2.0*a12*a13*a23)-(a11*a23*a23)-(a22*a13*a13)-(a33*a12*a12); // Mat. Determinant
			double c1 = a11*a22-a12*a12+a11*a33-a13*a13+a22*a33-a23*a23;
			double c2 = a11+a22+a33;
			double a = (3*c1-c2*c2)/3.0;
			double b = (9*c1*c2-2*c2*c2*c2-27*c0)/27.0;
			double Q = ((b*b)/4.0)+((a*a*a)/27.0);

			if (Q > 1.0e-12) { return new Vector3d(a11, a11, a11); } // Evecs = (1,0,0), (0,1,0), (0,0,1)

			double l1, l2, l3; // Eigenvalues to be computed

			if (Q>=0) {double p=(b>0)?(Math.Pow(b/2, 1.0/3.0)):(0); l1=l2=((c2/3.0)+p); l3=((c2/3.0)-2.0*p);}
			else
			{
			double t = Math.Atan2(Math.Sqrt(-Q), -b/2.0)/3.0, r = Math.Pow(((b*b)/4.0)-Q, 1.0/6.0);
			double cos_t = Math.Cos(t), sin_t = Math.Sin(t);
			double sq3 = Math.Sqrt(3.0);
			l1 = l2 = l3 = (c2/3.0);
			l1 += (2*r*cos_t);
			l2 -= r*(cos_t+sq3*sin_t);
			l3 -= r*(cos_t-sq3*sin_t);
			}

			double L1 = 0, L2 = 0, L3 = 0;
			if (l1<=l2 && l1<=l3) { L1=l1; L2=(l2<l3)?(l2):(l3); L3=(l2<l3)?(l3):(l2); }
			else if (l2<=l1 && l2<=l3) { L1=l2; L2=(l1<l3)?(l1):(l3); L3=(l1<l3)?(l3):(l1); }
			else { L1=l3; L2=(l1<l2)?(l1):(l2); L3=(l1<l2)?(l2):(l1); }
			return new Vector3d(L1, L2, L3);
		}
		public Vector3d GetMaxEigenvalues()
		{
			double a11 = -this.a, a12 = -this.b, a22 = -this.d, a13 = -this.c, a23 = -this.e, a33 = -this.f;
			double c0 = (a11 * a22 * a33) + (2.0 * a12 * a13 * a23) - (a11 * a23 * a23) - (a22 * a13 * a13) - (a33 * a12 * a12); // Mat. Determinant
			double c1 = a11 * a22 - a12 * a12 + a11 * a33 - a13 * a13 + a22 * a33 - a23 * a23;
			double c2 = a11 + a22 + a33;
			double a = (3 * c1 - c2 * c2) / 3.0;
			double b = (9 * c1 * c2 - 2 * c2 * c2 * c2 - 27 * c0) / 27.0;
			double Q = ((b * b) / 4.0) + ((a * a * a) / 27.0);

			if (Q > 1.0e-12) { return new Vector3d(a11, a11, a11); } // Evecs = (1,0,0), (0,1,0), (0,0,1)

			double l1, l2, l3; // Eigenvalues to be computed

			if (Q >= 0) { double p = (b > 0) ? (Math.Pow(b / 2, 1.0 / 3.0)) : (0); l1 = l2 = ((c2 / 3.0) + p); l3 = ((c2 / 3.0) - 2.0 * p); }
			else
			{
				double t = Math.Atan2(Math.Sqrt(-Q), -b / 2.0) / 3.0, r = Math.Pow(((b * b) / 4.0) - Q, 1.0 / 6.0);
				double cos_t = Math.Cos(t), sin_t = Math.Sin(t);
				double sq3 = Math.Sqrt(3.0);
				l1 = l2 = l3 = (c2 / 3.0);
				l1 += (2 * r * cos_t);
				l2 -= r * (cos_t + sq3 * sin_t);
				l3 -= r * (cos_t - sq3 * sin_t);
			}

			double L1 = 0, L2 = 0, L3 = 0;
			if (l1 <= l2 && l1 <= l3) { L1 = l1; L2 = (l2 < l3) ? (l2) : (l3); L3 = (l2 < l3) ? (l3) : (l2); }
			else if (l2 <= l1 && l2 <= l3) { L1 = l2; L2 = (l1 < l3) ? (l1) : (l3); L3 = (l1 < l3) ? (l3) : (l1); }
			else { L1 = l3; L2 = (l1 < l2) ? (l1) : (l2); L3 = (l1 < l2) ? (l2) : (l1); }
			return new Vector3d(L1, L2, L3);
		}
		public Vector3d GetMinEigenvector() 
		{
			double a11 = a, a12 = b, a22 = d, a13 = c, a23 = e, a33 = f;
			double l, l1, l2, l3, c0, c1, c2;
			Vector3d eig = GetEigenvalues();
			l = eig.x;
			l2 = eig.y;
			l3 = eig.z;

			if (l==l3 && l==l2) { return new Vector3d(1,0,0); }

			a11-=l; a22-=l; a33-=l;
			double u11 = a22*a33-a23*a23, u12 = a13*a23-a12*a33, u13 = a12*a23-a13*a22;
			double u22 = a11*a33-a13*a13, u23 = a12*a13-a23*a11, u33 = a11*a22-a12*a12;
			l1 = u11*u11+u12*u12+u13*u13;
			l2 = u12*u12+u22*u22+u23*u23;
			l3 = u13*u13+u23*u23+u33*u33;

			if (l1>=l2 && l1>=l3) {c0=u11; c1=u12; c2=u13; l=l1;}
			else if (l2>=l1 && l2>=l3) {c0=u12; c1=u22; c2=u23; l=l2;}
			else {c0=u13; c1=u23; c2=u33; l=l3;}

			l = Math.Sqrt(l);
			return new Vector3d(c0/l, c1/l, c2/l);
		}
		public Vector3d GetMaxEigenvector()
		{
			double a11 = -a, a12 = -b, a22 = -d, a13 = -c, a23 = -e, a33 = -f;
			double l, l1, l2, l3, c0, c1, c2;
			Vector3d eig = GetMaxEigenvalues();
			l = eig.x;
			l2 = eig.y;
			l3 = eig.z;

			if (l == l3 && l == l2) { return new Vector3d(1, 0, 0); }

			a11 -= l; a22 -= l; a33 -= l;
			double u11 = a22 * a33 - a23 * a23, u12 = a13 * a23 - a12 * a33, u13 = a12 * a23 - a13 * a22;
			double u22 = a11 * a33 - a13 * a13, u23 = a12 * a13 - a23 * a11, u33 = a11 * a22 - a12 * a12;
			l1 = u11 * u11 + u12 * u12 + u13 * u13;
			l2 = u12 * u12 + u22 * u22 + u23 * u23;
			l3 = u13 * u13 + u23 * u23 + u33 * u33;

			if (l1 >= l2 && l1 >= l3) { c0 = u11; c1 = u12; c2 = u13; l = l1; }
			else if (l2 >= l1 && l2 >= l3) { c0 = u12; c1 = u22; c2 = u23; l = l2; }
			else { c0 = u13; c1 = u23; c2 = u33; l = l3; }

			l = Math.Sqrt(l);
			return new Vector3d(c0 / l, c1 / l, c2 / l);
		}
		public bool Invert()
		{
			double det, pos, neg, t;
			double[] ans = new double[6];

			pos = neg = 0.0;
			t = a * d * f; if (t > 0) pos += t; else neg += t;
			t = b * e * c; if (t > 0) pos += t; else neg += t;
			t = c * b * e; if (t > 0) pos += t; else neg += t;
			t = -c * d * c; if (t > 0) pos += t; else neg += t;
			t = -b * b * f; if (t > 0) pos += t; else neg += t;
			t = -a * e * e; if (t > 0) pos += t; else neg += t;
			det = pos+neg;

			t = det/(pos-neg);
			if (Math.Abs(t) >= 1.0e-15)
			{
				ans[0] =  (d * f - e * e) /det;
				ans[1] = -(b * f - e * c) /det;
				ans[2] =  (a * f - c * c) /det;
				ans[3] =  (b * e - d * c) /det;
				ans[4] = -(a * e - b * c) /det;
				ans[5] =  (a * d - b * b) /det;
				a = ans[0]; b = ans[1]; d = ans[2];
				c = ans[3]; e = ans[4]; f = ans[5];
				return true;
			}

			return false;
		}
	}

	public class SymMatrix4d
	{
		public double a, b, c, d, e, f, g, h, i, j;

		public SymMatrix4d()
		{
			a = b = c = d = e = f = g = h = i= j = 0;
		}
		public SymMatrix4d(Vector4d v) : this(v.x, v.y, v.z, v.w)
		{
		}
		public SymMatrix4d(double x, double y, double z, double w)
		{
			a = x * x; b = x * y; c = x * z; d = x * w;
			e = y * y; f = y * z; g = y * w;
			h = z * z; i = z * w;
			j = w * w;
		}
		public SymMatrix4d(
			double a, double b, double c, double d, 
			double e, double f, double g, double h,
			double i, double j)
		{
			this.a = a;
			this.b = b;
			this.c = c;
			this.d = d;
			this.e = e;
			this.f = f;
			this.g = g;
			this.h = h;
			this.i = i;
			this.j = j;
		}

		public void Add(SymMatrix4d other)
		{
			this.a += other.a;
			this.b += other.b;
			this.c += other.c;
			this.d += other.d;
			this.e += other.e;
			this.f += other.f;
			this.g += other.g;
			this.h += other.h;
			this.i += other.i;
			this.j += other.j;
		}

		public static SymMatrix4d operator +(SymMatrix4d m, SymMatrix4d other)
		{
			return new SymMatrix4d(
				m.a + other.a,
				m.b + other.b,
				m.c + other.c,
				m.d + other.d,
				m.e + other.e,
				m.f + other.f,
				m.g + other.g,
				m.h + other.h,
				m.i + other.i,
				m.j + other.j
				);
		}
		public static SymMatrix4d operator *(SymMatrix4d m, double other)
		{
			return new SymMatrix4d(
				m.a * other,
				m.b * other,
				m.c * other,
				m.d * other,
				m.e * other,
				m.f * other,
				m.g * other,
				m.h * other,
				m.i * other,
				m.j * other
				);
		}
		public static Vector4d operator *(SymMatrix4d m, Vector4d other)
		{
			return new Vector4d(
				m.a * other.x + m.b * other.y + m.c * other.z + m.d * other.w,
				m.b * other.x + m.e * other.y + m.f * other.z + m.g * other.w,
				m.c * other.x + m.f * other.y + m.h * other.z + m.i * other.w,
				m.d * other.x + m.g * other.y + m.i * other.z + m.j * other.w);
		}
		public bool Invert()
		{
			 if (this.a <= 0) return false;

			 double d00 = 1.0/this.a; 
			 double L10 = this.b;
			 double l10 = this.b*d00;
			 double L20 = this.c;
			 double l20 = this.c*d00;
			 double L30 = this.d;
			 double l30 = this.d*d00;
			 double d11 = this.e-(L10*l10);

			 if (d11 <= 0) return false; else d11 = 1.0/d11;

			 double L21 = (this.f-(L10*l20));
			 double l21 = L21*d11;
			 double L31 = (this.g-(L10*l30));
			 double l31 = L31*d11;
			 double d22 = this.h-(L20*l20)-(L21*l21);

			 if (d22 <= 0) return false; else d22 = 1.0/d22;

			 double L32 = (this.i-(L20*l30)-(L21*l31));
			 double l32 = L32*d22;
			 double d33 = this.j-(L30*l30)-(L31*l31)-(L32*l32);

			 if (d33 <= 0) return false; else d33 = 1.0/d33;

			 L20 = l10*l21-l20;
			 L31 = l21*l32-l31;
			 L30 = l32*l20-L31*l10-l30;

			 this.a = d00+(-l10)*((-l10)*d11)+L20*(L20*d22)+L30*(L30*d33);
			 this.b = ((-l10)*d11)+L20*((-l21)*d22)+L30*(L31*d33);
			 this.e = d11+(-l21)*((-l21)*d22)+L31*(L31*d33);
			 this.c = (L20*d22)+L30*((-l32)*d33);
			 this.f = ((-l21)*d22)+L31*((-l32)*d33);
			 this.h = d22+(-l32)*((-l32)*d33);
			 this.d = (L30*d33);
			 this.g = (L31*d33);
			 this.i = ((-l32)*d33);
			 this.j = d33;

			 return true;
		}


	}

}