// XGrid3D, Clive Tooth. February 2013

using System.Collections.Generic;
using System.Collections;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.IO;
using System.Windows.Forms;
using System;
using Utils;

namespace MainSpace
{
  public class XGrid3D /////////////////////////////////////////////////////////
  {
    public class Line : IEquatable<Line> ///////////////////////////////////////
    {
      public static readonly Line zero = new Line();

      public Point p;
      public Point q;

      // An object with this class represents the line pq.

      public Line(
        Point p_param,
        Point q_param
      )
      {
        p = p_param;
        q = q_param;
      } // ctor Point Point

      public Line(): this(Point.zero, Point.zero)
      {
        //
      } // ctor

      public static bool operator ==(Line l0, Line l1)
      {
        return l0.Equals(l1);
      } // ==

      public static bool operator !=(Line l0, Line l1)
      {
        return !(l0 == l1);
      } // !=

      public bool Equals(Line l)
      {
 	      return On(l.p) && On(l.q);
      } // Equals(Line)

      public override bool Equals(object obj)
      {
 	      return obj is Line? Equals((Line)obj): false;
      } // Equals(object)

      public override int GetHashCode()
      {
        return p.GetHashCode() ^ q.GetHashCode();
      } // GetHashCode

      public bool On(Point p)
      {
        return XGrid3D.On(this, p);
      } // On

      public override string ToString()
      {
        return "["+p+" "+q+"]";
      } // ToString

    } // Line //////////////////////////////////////////////////////////////////

    public class Point : IEquatable<Point>, IComparable<Point> /////////////////
    {
      public static readonly Point zero = new Point();

      public int x = 0;
      public int y = 0;
      public int z = 0;
      public int t = 0;

      // An object of this class represents the point (x/t, y/t, z/t).
      // t must not be negative.
      // If any of x, y, z, t are non-zero then gcd(x, y, z, t) must be 1.
      // If t = 0 the point is at infinity.

      public Point(
        int x_param,
        int y_param,
        int z_param,
        int t_param
      )
      {
        x = x_param;
        y = y_param;
        z = z_param;
        t = t_param;
        Normalise();
      } // ctor int int int int

      public Point(
        int x_param,
        int y_param,
        int z_param
      )
      : this(x_param, y_param, z_param, 1)
      {
        //
      } // ctor int int int

      public Point(): this(0, 0, 0, 0)
      {
        //
      } // ctor

      public int CompareTo(Point p)
      {
        return
          x*p.t != p.x*t? x*p.t - p.x*t:
          y*p.t != p.y*t? y*p.t - p.y*t:
                          z*p.t - p.z*t;
      } // CompareTo

      public static bool operator <(Point p0, Point p1)
      {
        return p0.CompareTo(p1) < 0;
      } // <

      public static bool operator <=(Point p0, Point p1)
      {
        return !(p0 > p1);
      } // <=

      public static bool operator ==(Point p0, Point p1)
      {
        return p0.Equals(p1);
      } // ==

      public static bool operator !=(Point p0, Point p1)
      {
        return !(p0 == p1);
      } // !=

      public static bool operator >=(Point p0, Point p1)
      {
        return p1 <= p0;
      } // >=

      public static bool operator >(Point p0, Point p1)
      {
        return p1 < p0;
      } // >

      public static Point operator +(Point p0, Point p1)
      {
        return new Point(p0.t*p1.x+p1.t*p0.x, p0.t*p1.y+p1.t*p0.y, p0.t*p1.z+p1.t*p0.z, p0.t*p1.t);
      } // +

      public static Point operator -(Point p)
      {
        return new Point(p.x, p.y, p.z, -p.t);
      } // - Point

      public static Point operator -(Point p0, Point p1)
      {
        return p0 + -p1;
      } // - Point Point

      public static Point operator /(Point p, int m)
      {
        return new Point(p.x, p.y, p.z, m*p.t);
      } // /

      public static Point operator *(Point p, int m)
      {
        return new Point(m*p.x, m*p.y, m*p.z, p.t);
      } // *

      public static Point operator *(int m, Point p)
      {
        return p*m;
      } // *

      public Clive3Dpoint ClivePoint()
      {
        return new Clive3Dpoint(U.D(x, t), U.D(y, t), U.D(z, t));
      } // ClivePoint

      public bool Equals(Point p) 
      {
 	      return x == p.x && y == p.y && z == p.z && t == p.t;
      } // Equals Point

      public override bool Equals(object obj)
      {
        return obj is Point? Equals((Point)obj): false;
      } // Equals object

      public override int GetHashCode()
      {
        return x^y^z^t;
      } // GetHashCode

      public void Normalise()
      {
        if (
          t  < 0 ||
          t == 0 && x  < 0 ||
          t == 0 && x == 0 && y  < 0 ||
          t == 0 && x == 0 && y == 0 && z  < 0 ||
          false
        )
        {
          x = -x;
          y = -y;
          z = -z;
          t = -t;
        }

        int gcd = U.gcd(x, y, z, t);

        if (gcd > 0)
        {
          x /= gcd;
          y /= gcd;
          z /= gcd;
          t /= gcd;
        }
      } // Normalise

      public bool On(Line l)
      {
        return XGrid3D.On(l, this);
      } // On

      public override string ToString()
      {
        return t == 1? "("+x+" "+y+" "+z+")": "("+x+" "+y+" "+z+" "+t+")";
      } // ToString

    } // Point /////////////////////////////////////////////////////////////////    
    
    public List<Line>   allLines = null;
    public List<Point>  allPoints = null;
    public int          n = 0;
    public int          nn = 0;
    public int          nnn = 0;
    public int          numLines = 0;
    public int          numPoints = 0;

    public XGrid3D(
      int n_param
    )
    {
      n = n_param;

      nn = n*n;
      nnn = n*n*n;
      FindLines();
      FindPoints();
    } // ctor

    public static bool Coplanar(Point p0, Point p1, Point p2, Point p3)
    {
      return
        U.Determinant(
          p0.x, p0.y, p0.z, p0.t, 
          p1.x, p1.y, p1.z, p1.t, 
          p2.x, p2.y, p2.z, p2.t, 
          p3.x, p3.y, p3.z, p3.t 
        ) == 0;
    } // Coplanar Point Point Point Point

    public static bool Coplanar(Line l0, Line l1)
    {
      return Coplanar(l0.p, l0.q, l1.p, l1.q);
    } // Coplanar Line Line

    public void FindLines()
    {
      allLines = new List<Line>();
      for (int i = 0; i < nnn; i++)
      {
        Point p0 = PointFromInt(i);

        for (int j = i+1; j < nnn; j++)
        {
          Line l = new Line(p0, PointFromInt(j));

          if (!allLines.Contains(l))
          {
            allLines.Add(l);
          }
        } // for j
      } // for i
      numLines = allLines.Count;
    } // FindLines

    public void FindPoints()
    {
      allPoints = new List<Point>();
      for (int i = 0; i < nnn; i++)
      {
        allPoints.Add(PointFromInt(i));
      }
      for (int i = 0; i < numLines; i++)
      {
        for (int j = i+1; j < numLines; j++)
        {
          Point p = Intersection(allLines[i], allLines[j]);
          if (p.t == 0 || allPoints.Contains(p))
          {
            //
          }
          else
          {
            allPoints.Add(p);
          }
        }
      }
      numPoints = allPoints.Count;
    } // FindPoints

    public static Point Intersection(Line l0, Line l1)
    {
      if (Coplanar(l0, l1) && l0 != l1)
      {
        int a, b;

        a = l1.p.y*l1.q.x - l1.p.x*l1.q.y + l0.q.y*(l1.p.x - l1.q.x) + l0.q.x*(-l1.p.y + l1.q.y);
        b = (l0.p.x - l0.q.x)*(l1.p.y - l1.q.y) - (l0.p.y - l0.q.y)*(l1.p.x - l1.q.x);
        if (a == 0 && b == 0)
        {
          a = l1.p.z*l1.q.y - l1.p.y*l1.q.z + l0.q.z*(l1.p.y - l1.q.y) + l0.q.y*(-l1.p.z + l1.q.z);
          b = (l0.p.y - l0.q.y)*(l1.p.z - l1.q.z) - (l0.p.z - l0.q.z)*(l1.p.y - l1.q.y);
          if (a == 0 && b== 0)
          {
            a = l1.p.x*l1.q.z - l1.p.z*l1.q.x + l0.q.x*(l1.p.z - l1.q.z) + l0.q.z*(-l1.p.x + l1.q.x);
            b = (l0.p.z - l0.q.z)*(l1.p.x - l1.q.x) - (l0.p.x - l0.q.x)*(l1.p.z - l1.q.z);
            if (a == 0 && b== 0)
            {
              U.GiveUp("Intersection problem");
            }
          }
        }
        return (a*l0.p + (b-a)*l0.q)/b;
      }
      return new Point();
    } // Intersection

    public static bool On(Line l, Point r)
    {
      int nx = l.p.t*(l.q.t*r.x - l.q.x*r.t);
      int dx = r.t*(l.p.x*l.q.t - l.p.t*l.q.x);
      int ny = l.p.t*(l.q.t*r.y - l.q.y*r.t);
      int dy = r.t*(l.p.y*l.q.t - l.p.t*l.q.y);
      int nz = l.p.t*(l.q.t*r.z - l.q.z*r.t);
      int dz = r.t*(l.p.z*l.q.t - l.p.t*l.q.z);

      return nx*dy == ny*dx && ny*dz == nz*dy && nz*dx == nx*dz; // All three tests are needed.
    } // On

    private int PointComparer(int x, int y)
    {
      return allPoints[x].CompareTo(allPoints[y]);
    } // PointComparer

    public Point PointFromInt(int a)
    {
      return new Point(a/nn, a/n%n, a%n);
    } // PointFromInt

  } // XGrid3D

} // MainSpace