unit PolyMath;

{$modeswitch autoderef on}

interface

uses
         Math, SysUtils;

const
         Eps = 0.001;

type
         float = single;

         function Interpolate(const a, b: float; const x: float): float;
         function Clamp(const x: float; const a, b: float): float;
         function FloatSign(const x: float): float;

type
         PVec3 = ^TVec3;
         TVec3 = object
         type
           TData = array[0 .. 2] of float;
         var
           data: TData;
           property X: float read data[0] write data[0];
           property Y: float read data[1] write data[1];
           property Z: float read data[2] write data[2];
           function Make(const vx, vy, vz: float): TVec3; static;
           function SqrLength: float;
           function Length: float;
           function Normalized: TVec3;
           function MaybeNormalized: TVec3; // вектор по семантике нормализованный, но мог денормализоваться из-за ошибок округления
           function IsNormalized: boolean;
           function IsZero: boolean;
           function RandomDirection: TVec3; static;
         const
           Zero: TVec3 = (data: (0, 0, 0));
         type
           TArray = array of TVec3;
         var
           function Split(const v: array of float): TArray; static;
         end;
         operator +(const a, b: TVec3): TVec3;
         operator -(const a, b: TVec3): TVec3;
         operator -(const v: TVec3): TVec3;
         operator *(const a, b: TVec3): TVec3;
         operator *(const k: float; const v: TVec3): TVec3;
         operator *(const v: TVec3; const k: float): TVec3;
         operator /(const v: TVec3; const k: float): TVec3;
         operator **(const a, b: TVec3): float; // скалярное произведение
         operator ><(const a, b: TVec3): TVec3; // векторное произведение
         function Interpolate(const a, b: TVec3; const pos: float): TVec3;
         function Distance(const a, b: TVec3): float;
         function Angle(const a, b: TVec3): float;

type
         TMat4 = object
         type
           TData = array[0 .. 3, 0 .. 3] of float;
         var
           data: TData;
           function ExplicitRows(const x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, xa, xb, xc, xd, xe, xf: float): TMat4; static;
           function FrustumProjection(const left, right, bottom, top, zNear, zFar: float): TMat4; static;
           function PerspectiveProjection(const yFov, aspect, zNear, zFar: float): TMat4; static;
           function Translation(const v: TVec3): TMat4; static;
           function Rotation(const angle: float; const axis: TVec3): TMat4; static;
           function Rotation(const a, b: TVec3): TMat4; static;
           function Determinant: float;
           function Inversed: TMat4;

           // Частный случай — матрица поворота по направлениям вправо, вверх и назад.
           function WithBasis(const x, y, z: TVec3): TMat4; static;
         const
           Identity: TMat4 = (data: ((1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 1, 0), (0, 0, 0, 1)));
         end;
         operator *(const mat: TMat4; const vec: TVec3): TVec3;
         operator *(const a, b: TMat4): TMat4;

type
         TPlane = object
         type
           TData = record
           case integer of
             0: (linear: array[0 .. 3] of float);
             1: (v3: TVec3);
           end;
         var
           data: TData;
           property A: float read data.linear[0] write data.linear[0];
           property B: float read data.linear[1] write data.linear[1];
           property C: float read data.linear[2] write data.linear[2];
           property D: float read data.linear[3] write data.linear[3];
           property ABC: TVec3 read data.v3 write data.v3;
           function Make(const pa, pb, pc: TVec3): TPlane; static;
           function Make(const point, normal: TVec3): TPlane; static;
           function AnyPoint: TVec3;
           function IntersectLine(const p, v: TVec3; out int: TVec3): boolean;
           function IntersectSegment(const pa, pb: TVec3; out int: TVec3): boolean;
         end;

         PFace = ^TFace;
         TFace = object
           v: array of integer;
           function Make(const newV: array of integer): TFace; static;
           function Validate(nv: integer): TFace;
           function CyclicV(i: integer): integer;

         type
           TArray = array of TFace;
           function Split(vertsPerFace: integer; const vs: array of integer): TArray; static;
         end;

         PPolyhedron = ^TPolyhedron;
         TPolyhedron = object
           vts: array of TVec3;
           faces: array of TFace;

           // Вершины граней обходятся против часовой.
           function Make(const newVts: array of TVec3; const scale: float; const newFaces: array of TFace): TPolyhedron; static;
           function Make(const newVts: array of TVec3; const newFaces: array of TFace): TPolyhedron; static;
           function Vertex(face: PFace; id: integer): TVec3;
           function Normal(face: PFace): TVec3;
           function AdjacentFace(current: PFace; edge: integer): PFace;
           function EdgeVector(face: PFace; id: integer): TVec3;
           function MaxRadiusVector: float;

         type
           TSingleFaceIntersection = record
             pos: TVec3;
             edge: integer;
           end;
           TFaceIntersection = array of TSingleFaceIntersection;
           function IntersectFace(face: PFace; const plane: TPlane): TFaceIntersection;
           function SameEdge(faceA: PFace; edgeA: integer; faceB: PFace; edgeB: integer): boolean;
           function Icosahedron(const radius: float): TPolyhedron; static;
           function Cube(const w, h, d: float): TPolyhedron; static;
           function TwoSidedHexagon(const radius: float): TPolyhedron; static;
         end;

implementation

         function Interpolate(const a, b: float; const x: float): float;
         begin
           result := a + x * (b - a);
         end;

         function Clamp(const x: float; const a, b: float): float;
         begin
           if x < a then result := a else
             if x > b then result := b else
               result := x;
         end;

         function FloatSign(const x: float): float;
         begin
           if x > 0 then result := 1 else
             if x < 0 then result := -1 else
               result := 0;
         end;

         function TVec3.Make(const vx, vy, vz: float): TVec3;
         begin
           result.x := vx;
           result.y := vy;
           result.z := vz;
         end;

         function TVec3.SqrLength: float;
         begin
           result := sqr(x) + sqr(y) + sqr(z);
         end;

         function TVec3.Length: float;
         begin
           result := sqrt(SqrLength);
         end;

         function TVec3.Normalized: TVec3;
         var
           len: float;
         begin
           len := Length;
           if len = 0.0 then
             result := self
           else
             result := self / len;
           Assert(result.IsNormalized or result.IsZero);
         end;

         function TVec3.MaybeNormalized: TVec3;
         begin
           if IsNormalized then result := self else result := Normalized;
         end;

         function TVec3.IsNormalized: boolean;
         begin
           result := abs(SqrLength - 1) < Eps;
         end;

         function TVec3.IsZero: boolean;
         begin
           result := (abs(x) < Eps) and (abs(y) < Eps) and (abs(z) < Eps);
         end;

         function TVec3.RandomDirection: TVec3;
         var
           dz, len, angle: float;
         begin
           dz := 2 * random - 1;
           len := sqrt(1 - sqr(dz));
           angle := random * (2 * pi);
           result := Make(cos(angle) * len, sin(angle) * len, dz);
         end;

         function TVec3.Split(const v: array of float): TArray;
         var
           i, iv: integer;
         begin
           SetLength(result, System.length(v) div System.length(TData));
           iv := 0;
           for i := 0 to High(result) do
           begin
             result[i] := TVec3.Make(v[iv], v[iv + 1], v[iv + 2]);
             iv += 3;
           end;
         end;

{$define impl :=
         var
           i: integer;
         begin
           for i := 0 to High(result.data) do
             result.data[i] := {$ifdef un} un {$else} a.data[i] op b.data[i] {$endif};
         end; {$undef op} {$undef un}}
         operator +(const a, b: TVec3): TVec3; {$define op := +} impl
         operator -(const a, b: TVec3): TVec3; {$define op := -} impl
         operator -(const v: TVec3): TVec3; {$define un := -v.data[i]} impl
         operator *(const a, b: TVec3): TVec3; {$define op := *} impl
{$undef impl}

         operator *(const k: float; const v: TVec3): TVec3;
         begin
           result := TVec3.Make(k * v.x, k * v.y, k * v.z);
         end;

         operator *(const v: TVec3; const k: float): TVec3;
         begin
           result := k * v;
         end;

         operator /(const v: TVec3; const k: float): TVec3;
         begin
           result := (1.0 / k) * v;
         end;

         operator **(const a, b: TVec3): float;
         begin
           result := a.x * b.x + a.y * b.y + a.z * b.z;
         end;

         operator ><(const a, b: TVec3): TVec3;
         begin
           result.x := a.y * b.z - a.z * b.y;
           result.y := a.z * b.x - a.x * b.z;
           result.z := a.x * b.y - a.y * b.x;
         end;

         function Interpolate(const a, b: TVec3; const pos: float): TVec3;
         begin
           result := a + pos * (b - a);
         end;

         function Distance(const a, b: TVec3): float;
         begin
           result := (b - a).Length;
         end;

         function Angle(const a, b: TVec3): float;
         var
           len: float;
         begin
           len := sqrt(a.SqrLength * b.SqrLength);
           if len < Eps then result := 0 else result := ArcCos((a ** b) / len);
         end;
...