using System;

using System.Collections.Generic;

using System.Linq.Expressions;

namespace EllipseTruncate {

	public static class Geometry {

		public static readonly double precision = 1e-8;

		public static readonly double subPrecision = precision*0.1;

		public static readonly double precisionSqr = precision*precision;

		public static readonly double sqrt2Pi = Math.Sqrt(2.0*Math.PI);

		public static bool isEqual(double a, double b)

			{ return Math.Abs(b - a) <= precision; }

		public static bool isLess(double a, double b)

			{ return a - precision < b; }

		public static bool isGreater(double a, double b)

			{ return b - precision < a; }

		public static bool isLessOrEqual(double a, double b)

			{ return a + precision <= b; }

		public static bool isGreaterOrEqual(double a, double b)

			{ return b + precision <= a; }

		public static double logNormalDistribuitionUnscaled(double x, double x0, double w) {

			return Math.Exp(-0.5*Math.Pow(Math.Log(x/x0)/w, 2.0))/x;

		}

		public static double logNormalDistribuition(double x, double x0, double w) {

			return logNormalDistribuition(x, x0, w)/(w*sqrt2Pi);

		}

		public static void truncateInfiniteLine(Rectangle bounds, ref Point p0, ref Point p1) {

			if (p0.isEqual(p1)) return;

			Point d = p0 - p1;

			if (Math.Abs(d.x)*bounds.height > bounds.width*Math.Abs(d.y)) {

				// horizontal

				double k = d.y/d.x;

				p1 = new Point(bounds.x1, p0.y + k*(bounds.x1 - p0.x));

				p0 = new Point(bounds.x0, p0.y + k*(bounds.x0 - p0.x));

			} else {

				// vertical

				double k = d.x/d.y;

				p1 = new Point(p0.x + k*(bounds.y1 - p0.y), bounds.y1);

				p0 = new Point(p0.x + k*(bounds.y0 - p0.y), bounds.y0);

			}

		}

		public static class Interpolation<t> {</t>

			public delegate T HalfFunc(T a, double b);

			public delegate T FullFunc(T a, T b);

			public static FullFunc add;

			public static FullFunc sub;

			public static HalfFunc mul;

			public static HalfFunc div;

   			static Interpolation() {

   				{ // add

   					ParameterExpression a = Expression.Parameter(typeof(T));

				    ParameterExpression b = Expression.Parameter(typeof(T));

				    add = Expression.Lambda<fullfunc>(Expression.Add(a, b), a, b).Compile();</fullfunc>

				}

   				{ // sub

   					ParameterExpression a = Expression.Parameter(typeof(T));

				    ParameterExpression b = Expression.Parameter(typeof(T));

				    sub = Expression.Lambda<fullfunc>(Expression.Subtract(a, b), a, b).Compile();</fullfunc>

				}

   				{ // mul

   					ParameterExpression a = Expression.Parameter(typeof(T));

				    ParameterExpression b = Expression.Parameter(typeof(double));

				    mul = Expression.Lambda<halffunc>(Expression.Multiply(a, b), a, b).Compile();</halffunc>

				}

   				{ // div

   					ParameterExpression a = Expression.Parameter(typeof(T));

				    ParameterExpression b = Expression.Parameter(typeof(double));

				    div = Expression.Lambda<halffunc>(Expression.Divide(a, b), a, b).Compile();</halffunc>

				}

   			}

			public static T linear(T p0, T p1, double l)

				{ return add(mul(p0, 1.0-l), mul(p1, l)); }

			public static T spline(T p0, T p1, T t0, T t1, double l) {

				double ll = l*l;

				double lll = ll*l;

				return add( add( mul(p0, ( 2.0*lll - 3.0*ll + 1.0)),

				                 mul(p1, (-2.0*lll + 3.0*ll      )) ),

				            add( mul(t0, (     lll - 2.0*ll + l  )),

				                 mul(t1, (     lll - 1.0*ll      )) ));

	        }

			public static T splineTangent(T p0, T p1, T t0, T t1, double l) {

				double ll = l*l;

				return add( mul(sub(p0, p1), 6.0*(ll - l)),

				            add( mul(t0, ( 3.0*ll - 4.0*l + 1.0)),

				                 mul(t1, ( 3.0*ll - 2.0*l      )) ));

	        }

		}

		public static double interpolationLinear(double p0, double p1, double l)

			{ return p0*(1.0 - l) + p1*l; }

		public static double interpolationSpline(double p0, double p1, double t0, double t1, double l) {

			double ll = l*l;

			double lll = ll*l;

			return p0*( 2.0*lll - 3.0*ll + 1.0)

			     + p1*(-2.0*lll + 3.0*ll      )

			     + t0*(     lll - 2.0*ll + l  )

			     + t1*(     lll - 1.0*ll      );

        }

		public static double interpolationSplineTangent(double p0, double p1, double t0, double t1, double l) {

			double ll = l*l;

			return (p0 - p1)*6.0*(ll - l)

			     + t0*( 3.0*ll - 4.0*l + 1.0)

			     + t1*( 3.0*ll - 2.0*l      );

        }

		public static Point interpolationLinear(Point p0, Point p1, double l)

			{ return p0*(1.0 - l) + p1*l; }

		public static Point interpolationSpline(Point p0, Point p1, Point t0, Point t1, double l) {

			double ll = l*l;

			double lll = ll*l;

			return p0*( 2.0*lll - 3.0*ll + 1.0)

			     + p1*(-2.0*lll + 3.0*ll      )

			     + t0*(     lll - 2.0*ll + l  )

			     + t1*(     lll - 1.0*ll      );

        }

		public static Point interpolationSplineTangent(Point p0, Point p1, Point t0, Point t1, double l) {

			double ll = l*l;

			return (p0 - p1)*6.0*(ll - l)

			     + t0*( 3.0*ll - 4.0*l + 1.0)

			     + t1*( 3.0*ll - 2.0*l      );

        }

	}

}