Переименование проектов и оптимизация вычислений / Renaming of projects and optimization of computations

BasicGeometry/Angle.cs

@@ -0,0 +1,48 @@
+/*
+ * Copyright 2019-2025 Andrey Pokidov <andrey.pokidov@gmail.com>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 1 Feb 2019
+ */
+
+namespace BasicGeometry
+{
+    public enum AngleUnit
+    {
+        RADIANS = 1,
+        DEGREES = 2,
+        TURNS = 3,
+
+    }
+
+    public enum AngleRange
+    {
+        /// <summary>
+        /// The measure of an angle with a range of:
+        /// [0, 360) degrees, [0, 2xPI) radians, [0, 1) turns, [0, 400) gradians
+        /// </summary>
+        UNSIGNED_RANGE = 1,
+
+        /// <summary>
+        /// The measure of an angle with a range of:
+        /// (-180, 180] degrees, (-PI, PI] radians, (-0.5, 0.5] turns, (-200, 200] gradians
+        /// </summary>
+        SIGNED_RANGE = 2
+    }
+}
+

BasicGeometry/BasicGeometry.csproj

@@ -0,0 +1,10 @@
+<Project Sdk="Microsoft.NET.Sdk">
+
+  <PropertyGroup>
+    <TargetFramework>net6.0</TargetFramework>
+    <ImplicitUsings>enable</ImplicitUsings>
+    <Nullable>enable</Nullable>
+  </PropertyGroup>
+
+
+</Project>

BasicGeometry/FP32Angle.cs

@@ -0,0 +1,133 @@
+/*
+ * Copyright 2019-2025 Andrey Pokidov <andrey.pokidov@gmail.com>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+using System;
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 1 Feb 2019
+ */
+
+namespace BasicGeometry
+{
+    public static class FP32Angle
+    {
+        public static float ToRadians(float angle, AngleUnit unit)
+        {
+            if (unit == AngleUnit.DEGREES)
+            {
+                return angle * FP32Degrees.RADIANS_IN_DEGREE;
+            }
+
+            if (unit == AngleUnit.TURNS)
+            {
+                return angle * FP32Radians.TWO_PI;
+            }
+
+            return angle;
+        }
+
+        public static float ToDegrees(float angle, AngleUnit unit)
+        {
+            if (unit == AngleUnit.RADIANS)
+            {
+                return angle * FP32Radians.DEGREES_IN_RADIAN;
+            }
+
+            if (unit == AngleUnit.TURNS)
+            {
+                return angle * 360.0f;
+            }
+
+            return angle;
+        }
+
+        public static float ToTurns(float angle, AngleUnit unit)
+        {
+            if (unit == AngleUnit.RADIANS)
+            {
+                return angle * FP32Radians.TURNS_IN_RADIAN;
+            }
+
+            if (unit == AngleUnit.DEGREES)
+            {
+                return angle * FP32Degrees.TURNS_IN_DEGREE;
+            }
+
+            return angle;
+        }
+
+        public static float GetFullCircle(AngleUnit unit)
+        {
+            if (unit == AngleUnit.DEGREES)
+            {
+                return 360.0f;
+            }
+
+            if (unit == AngleUnit.TURNS)
+            {
+                return 1.0f;
+            }
+
+            return FP32Radians.TWO_PI;
+        }
+    
+        public static float GetHalfCircle(AngleUnit unit)
+        {
+            if (unit == AngleUnit.DEGREES)
+            {
+                return 180.0f;
+            }
+
+            if (unit == AngleUnit.TURNS)
+            {
+                return 0.5f;
+            }
+
+            return FP32Radians.PI;
+        }
+    
+        public static float GetQuarterCircle(AngleUnit unit)
+        {
+            if (unit == AngleUnit.DEGREES)
+            {
+                return 90.0f;
+            }
+
+            if (unit == AngleUnit.TURNS)
+            {
+                return 0.25f;
+            }
+
+            return FP32Radians.HALF_OF_PI;
+        }
+
+        public static float Normalize(float angle, AngleUnit unit, AngleRange range)
+        {
+            if (unit == AngleUnit.DEGREES)
+            {
+                return FP32Degrees.Normalize(angle, range);
+            }
+
+            if (unit == AngleUnit.TURNS)
+            {
+                return FP32Turns.Normalize(angle, range);
+            }
+
+            return FP32Radians.Normalize(angle, range);
+        }
+    }
+}

BasicGeometry/FP32Degrees.cs

@@ -0,0 +1,68 @@
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 18 Nov 2024
+ */
+
+namespace BasicGeometry
+{
+    public class FP32Degrees
+    {
+        public const float RADIANS_IN_DEGREE = 1.745329252E-2f;
+        public const float TURNS_IN_DEGREE = 2.7777777778E-3f;
+
+        public static float ToRadians(float degrees)
+        {
+            return degrees * RADIANS_IN_DEGREE;
+        }
+
+        public static float ToTurns(float degrees)
+        {
+            return degrees * TURNS_IN_DEGREE;
+        }
+
+        public static float ToUnits(float degrees, AngleUnit toUnit)
+        {
+            if (toUnit == AngleUnit.RADIANS)
+            {
+                return degrees * RADIANS_IN_DEGREE;
+            }
+
+            if (toUnit == AngleUnit.TURNS)
+            {
+                return degrees * TURNS_IN_DEGREE;
+            }
+
+            return degrees;
+        }
+
+        public static float Normalize(float radians, AngleRange range)
+        {
+            if (range == AngleRange.UNSIGNED_RANGE)
+            {
+                if (0.0f <= radians && radians < 360.0f)
+                {
+                    return radians;
+                }
+            }
+            else
+            {
+                if (-180.0f < radians && radians <= 180.0f)
+                {
+                    return radians;
+                }
+            }
+
+            float turns = radians * TURNS_IN_DEGREE;
+
+            turns -= MathF.Floor(turns);
+
+            if (range == AngleRange.SIGNED_RANGE && turns > 0.5f)
+            {
+                turns -= 1.0f;
+            }
+
+            return turns * 360.0f;
+        }
+    }
+}

BasicGeometry/FP32Matrix2x2.cs

@@ -0,0 +1,294 @@
+/*
+ * Copyright 2019-2025 Andrey Pokidov <andrey.pokidov@gmail.com>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+using System;
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 10 Feb 2019
+ */
+
+namespace BasicGeometry
+{
+    public struct FP32Matrix2x2
+    {
+        public float r1c1 = 0.0f, r1c2 = 0.0f;
+        public float r2c1 = 0.0f, r2c2 = 0.0f;
+
+        public FP32Matrix2x2(float d1, float d2)
+        {
+            this.r1c1 = d1;
+            this.r2c2 = d2;
+        }
+
+        public FP32Matrix2x2(in FP32Matrix2x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+        }
+
+        public FP32Matrix2x2(in FP64Matrix2x2 matrix)
+        {
+            this.r1c1 = (float)matrix.r1c1;
+            this.r1c2 = (float)matrix.r1c2;
+
+            this.r2c1 = (float)matrix.r2c1;
+            this.r2c2 = (float)matrix.r2c2;
+        }
+
+        public readonly float GetDeterminant()
+        {
+            return this.r1c1 * this.r2c2 - this.r1c2 * this.r2c1;
+        }
+
+        public readonly bool IsSingular()
+        {
+            float determinant = this.GetDeterminant();
+            return -FP32Utility.EPSYLON <= determinant && determinant <= FP32Utility.EPSYLON;
+        }
+
+        public void Transpose()
+        {
+            (this.r1c2, this.r2c1) = (this.r2c1, this.r1c2);
+        }
+
+        public bool Invert()
+        {
+            float determinant = this.GetDeterminant();
+
+            if (-FP32Utility.EPSYLON <= determinant && determinant <= FP32Utility.EPSYLON)
+            {
+                return false;
+            }
+
+            float r1c1 =  this.r2c2;
+            float r1c2 = -this.r1c2;
+
+            float r2c1 = -this.r2c1;
+            float r2c2 =  this.r1c1;
+
+            float multiplier = 1.0f / determinant;
+
+            this.r1c1 = r1c1 * multiplier;
+            this.r1c2 = r1c2 * multiplier;
+
+            this.r2c1 = r2c1 * multiplier;
+            this.r2c2 = r2c2 * multiplier;
+
+            return true;
+        }
+
+        public void Reset()
+        {
+            this.r1c1 = 0.0f;
+            this.r1c2 = 0.0f;
+
+            this.r2c1 = 0.0f;
+            this.r2c2 = 0.0f;
+        }
+
+        public void SetToIdentity()
+        {
+            this.r1c1 = 1.0f;
+            this.r1c2 = 0.0f;
+
+            this.r2c1 = 0.0f;
+            this.r2c2 = 1.0f;
+        }
+
+        public void SetToDiagonal(float d1, float d2)
+        {
+            this.r1c1 = d1;
+            this.r1c2 = 0.0f;
+
+            this.r2c1 = 0.0f;
+            this.r2c2 = d2;
+        }
+
+        public void SetValues(in FP32Matrix2x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+        }
+
+        public void SetValues(in FP64Matrix2x2 matrix)
+        {
+            this.r1c1 = (float)matrix.r1c1;
+            this.r1c2 = (float)matrix.r1c2;
+
+            this.r2c1 = (float)matrix.r2c1;
+            this.r2c2 = (float)matrix.r2c2;
+        }
+
+        public void SetTransposedOf(in FP32Matrix2x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r2c2 = matrix.r2c2;
+            (this.r1c2, this.r2c1) = (matrix.r2c1, matrix.r1c2);
+        }
+
+        public void SetTransposedOf(in FP64Matrix2x2 matrix)
+        {
+            this.r1c1 = (float)matrix.r1c1;
+            this.r1c2 = (float)matrix.r2c1;
+
+            this.r2c1 = (float)matrix.r1c2;
+            this.r2c2 = (float)matrix.r2c2;
+        }
+
+        public bool SetInvertedOf(in FP32Matrix2x2 matrix)
+        {
+            float determinant = matrix.GetDeterminant();
+
+            if (-FP32Utility.EPSYLON <= determinant && determinant <= FP32Utility.EPSYLON)
+            {
+                return false;
+            }
+
+            float r1c1 =  matrix.r2c2;
+            float r1c2 = -matrix.r1c2;
+
+            float r2c1 = -matrix.r2c1;
+            float r2c2 =  matrix.r1c1;
+
+            float multiplier = 1.0f / determinant;
+
+            this.r1c1 = r1c1 * multiplier;
+            this.r1c2 = r1c2 * multiplier;
+
+            this.r2c1 = r2c1 * multiplier;
+            this.r2c2 = r2c2 * multiplier;
+
+            return true;
+        }
+
+        public void SetRow1(float c1, float c2)
+        {
+            this.r1c1 = c1;
+            this.r1c2 = c2;
+        }
+
+        public void SetRow2(float c1, float c2)
+        {
+            this.r2c1 = c1;
+            this.r2c2 = c2;
+        }
+
+        public void SetColumn1(float r1, float r2)
+        {
+            this.r1c1 = r1;
+            this.r2c1 = r2;
+        }
+
+        public void SetColumn2(float r1, float r2)
+        {
+            this.r1c2 = r1;
+            this.r2c2 = r2;
+        }
+
+        public void AppendScaled(in FP32Matrix2x2 matrix, float scale)
+        {
+            this.r1c1 += matrix.r1c1 * scale;
+            this.r1c2 += matrix.r1c2 * scale;
+
+            this.r2c1 += matrix.r2c1 * scale;
+            this.r2c2 += matrix.r2c2 * scale;
+        }
+
+        public static void Add(in FP32Matrix2x2 matrix1, in FP32Matrix2x2 matrix2, out FP32Matrix2x2 result)
+        {
+            result.r1c1 = matrix1.r1c1 + matrix2.r1c1;
+            result.r1c2 = matrix1.r1c2 + matrix2.r1c2;
+
+            result.r2c1 = matrix1.r2c1 + matrix2.r2c1;
+            result.r2c2 = matrix1.r2c2 + matrix2.r2c2;
+        }
+
+        public static void Subtract(in FP32Matrix2x2 minuend, in FP32Matrix2x2 subtrahend, out FP32Matrix2x2 difference)
+        {
+            difference.r1c1 = minuend.r1c1 - subtrahend.r1c1;
+            difference.r1c2 = minuend.r1c2 - subtrahend.r1c2;
+
+            difference.r2c1 = minuend.r2c1 - subtrahend.r2c1;
+            difference.r2c2 = minuend.r2c2 - subtrahend.r2c2;
+        }
+
+        public static void Multiply(in FP32Matrix2x2 multiplicand, float multiplier, out FP32Matrix2x2 product)
+        {
+            product.r1c1 = multiplicand.r1c1 * multiplier;
+            product.r1c2 = multiplicand.r1c2 * multiplier;
+
+            product.r2c1 = multiplicand.r2c1 * multiplier;
+            product.r2c2 = multiplicand.r2c2 * multiplier;
+        }
+
+        public static void Divide(in FP32Matrix2x2 dividend, float divisor, out FP32Matrix2x2 quotient)
+        {
+            Multiply(dividend, 1.0f / divisor, out quotient);
+        }
+
+        public static void GetRightProduct(in FP32Matrix2x2 matrix, in FP32Vector2 vector, out FP32Vector2 result)
+        {
+            float x1 = matrix.r1c1 * vector.x1 + matrix.r1c2 * vector.x2;
+            float x2 = matrix.r2c1 * vector.x1 + matrix.r2c2 * vector.x2;
+
+            result.x1 = x1;
+            result.x2 = x2;
+        }
+
+        public static void GetLeftProduct(in FP32Vector2 vector, in FP32Matrix2x2 matrix, out FP32Vector2 result)
+        {
+            float x1 = vector.x1 * matrix.r1c1 + vector.x2 * matrix.r2c1;
+            float x2 = vector.x1 * matrix.r1c2 + vector.x2 * matrix.r2c2;
+
+            result.x1 = x1;
+            result.x2 = x2;
+        }
+
+        public static void LoadZero(out FP32Matrix2x2 matrix)
+        {
+            matrix.r1c1 = 0.0f;
+            matrix.r1c2 = 0.0f;
+
+            matrix.r2c1 = 0.0f;
+            matrix.r2c2 = 0.0f;
+        }
+
+        public static void LoadIdentity(out FP32Matrix2x2 matrix)
+        {
+            matrix.r1c1 = 1.0f;
+            matrix.r1c2 = 0.0f;
+
+            matrix.r2c1 = 0.0f;
+            matrix.r2c2 = 1.0f;
+        }
+
+        public static void LoadDiagonal(float d1, float d2, out FP32Matrix2x2 matrix)
+        {
+            matrix.r1c1 = d1;
+            matrix.r1c2 = 0.0f;
+
+            matrix.r2c1 = 0.0f;
+            matrix.r2c2 = d2;
+        }
+    }
+}

BasicGeometry/FP32Matrix2x3.cs

@@ -0,0 +1,237 @@
+/*
+ * Copyright 2019-2025 Andrey Pokidov <andrey.pokidov@gmail.com>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+using System;
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 11 Nov 2024
+ */
+namespace BasicGeometry
+{
+    public struct FP32Matrix2x3
+    {
+        public float r1c1 = 0.0f, r1c2 = 0.0f;
+        public float r2c1 = 0.0f, r2c2 = 0.0f;
+        public float r3c1 = 0.0f, r3c2 = 0.0f;
+    
+        public FP32Matrix2x3(in FP32Matrix2x3 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+        
+            this.r3c1 = matrix.r3c1;
+            this.r3c2 = matrix.r3c2;
+        }
+    
+        public FP32Matrix2x3(in FP64Matrix2x3 matrix)
+        {
+            this.r1c1 = (float) matrix.r1c1;
+            this.r1c2 = (float) matrix.r1c2;
+
+            this.r2c1 = (float) matrix.r2c1;
+            this.r2c2 = (float) matrix.r2c2;
+        
+            this.r3c1 = (float) matrix.r3c1;
+            this.r3c2 = (float) matrix.r3c2;
+        }
+    
+        public FP32Matrix2x3(in FP32Matrix3x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r2c1;
+
+            this.r2c1 = matrix.r1c2;
+            this.r2c2 = matrix.r2c2;
+        
+            this.r3c1 = matrix.r1c3;
+            this.r3c2 = matrix.r2c3;
+        }
+    
+        public FP32Matrix2x3(in FP64Matrix3x2 matrix)
+        {
+            this.r1c1 = (float) matrix.r1c1;
+            this.r1c2 = (float) matrix.r2c1;
+
+            this.r2c1 = (float) matrix.r1c2;
+            this.r2c2 = (float) matrix.r2c2;
+        
+            this.r3c1 = (float) matrix.r1c3;
+            this.r3c2 = (float) matrix.r2c3;
+        }
+    
+        public void Reset()
+        {
+            this.r1c1 = 0.0f;
+            this.r1c2 = 0.0f;
+        
+            this.r2c1 = 0.0f;
+            this.r2c2 = 0.0f;
+        
+            this.r3c1 = 0.0f;
+            this.r3c2 = 0.0f;
+        }
+    
+        public void SetValues(in FP32Matrix2x3 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+        
+            this.r3c1 = matrix.r3c1;
+            this.r3c2 = matrix.r3c2;
+        }
+    
+        public void SetValues(in FP64Matrix2x3 matrix)
+        {
+            this.r1c1 = (float) matrix.r1c1;
+            this.r1c2 = (float) matrix.r1c2;
+
+            this.r2c1 = (float) matrix.r2c1;
+            this.r2c2 = (float) matrix.r2c2;
+        
+            this.r3c1 = (float) matrix.r3c1;
+            this.r3c2 = (float) matrix.r3c2;
+        }
+    
+        public void SetTransposed(in FP32Matrix3x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r2c1;
+
+            this.r2c1 = matrix.r1c2;
+            this.r2c2 = matrix.r2c2;
+        
+            this.r3c1 = matrix.r1c3;
+            this.r3c2 = matrix.r2c3;
+        }
+    
+        public void SetTransposed(in FP64Matrix3x2 matrix)
+        {
+            this.r1c1 = (float) matrix.r1c1;
+            this.r1c2 = (float) matrix.r2c1;
+
+            this.r2c1 = (float) matrix.r1c2;
+            this.r2c2 = (float) matrix.r2c2;
+        
+            this.r3c1 = (float) matrix.r1c3;
+            this.r3c2 = (float) matrix.r2c3;
+        }
+
+        public void SetRow1(float c1, float c2)
+        {
+            this.r1c1 = c1;
+            this.r1c2 = c2;
+        }
+
+        public void SetRow2(float c1, float c2)
+        {
+            this.r2c1 = c1;
+            this.r2c2 = c2;
+        }
+
+        public void SetRow3(float c1, float c2)
+        {
+            this.r3c1 = c1;
+            this.r3c2 = c2;
+        }
+
+        public void SetColumn1(float r1, float r2, float r3)
+        {
+            this.r1c1 = r1;
+            this.r2c1 = r2;
+            this.r3c1 = r3;
+        }
+
+        public void SetColumn2(float r1, float r2, float r3)
+        {
+            this.r1c2 = r1;
+            this.r2c2 = r2;
+            this.r3c2 = r3;
+        }
+
+        public void AppendScaled(in FP32Matrix2x3 matrix, float scale)
+        {
+            this.r1c1 += matrix.r1c1 * scale;
+            this.r1c2 += matrix.r1c2 * scale;
+
+            this.r2c1 += matrix.r2c1 * scale;
+            this.r2c2 += matrix.r2c2 * scale;
+
+            this.r3c1 += matrix.r3c1 * scale;
+            this.r3c2 += matrix.r3c2 * scale;
+        }
+
+        public static void Add(in FP32Matrix2x3 matrix1, in FP32Matrix2x3 matrix2, out FP32Matrix2x3 sum)
+        {
+            sum.r1c1 = matrix1.r1c1 + matrix2.r1c1;
+            sum.r1c2 = matrix1.r1c2 + matrix2.r1c2;
+
+            sum.r2c1 = matrix1.r2c1 + matrix2.r2c1;
+            sum.r2c2 = matrix1.r2c2 + matrix2.r2c2;
+
+            sum.r3c1 = matrix1.r3c1 + matrix2.r3c1;
+            sum.r3c2 = matrix1.r3c2 + matrix2.r3c2;
+        }
+
+        public static void Subtract(in FP32Matrix2x3 minuend, in FP32Matrix2x3 subtrahend, out FP32Matrix2x3 difference)
+        {
+            difference.r1c1 = minuend.r1c1 - subtrahend.r1c1;
+            difference.r1c2 = minuend.r1c2 - subtrahend.r1c2;
+
+            difference.r2c1 = minuend.r2c1 - subtrahend.r2c1;
+            difference.r2c2 = minuend.r2c2 - subtrahend.r2c2;
+
+            difference.r3c1 = minuend.r3c1 - subtrahend.r3c1;
+            difference.r3c2 = minuend.r3c2 - subtrahend.r3c2;
+        }
+
+        public static void Multiply(in FP32Matrix2x3 multiplicand, float multiplier, out FP32Matrix2x3 product)
+        {
+            product.r1c1 = multiplicand.r1c1 * multiplier;
+            product.r1c2 = multiplicand.r1c2 * multiplier;
+
+            product.r2c1 = multiplicand.r2c1 * multiplier;
+            product.r2c2 = multiplicand.r2c2 * multiplier;
+
+            product.r3c1 = multiplicand.r3c1 * multiplier;
+            product.r3c2 = multiplicand.r3c2 * multiplier;
+        }
+
+        public static void Divide(in FP32Matrix2x3 dividend, float divisor, out FP32Matrix2x3 quotient)
+        {
+            Multiply(dividend, 1.0f / divisor, out quotient);
+        }
+
+        public static void GetRightProduct(in FP32Matrix2x3 matrix, in FP32Vector2 vector, out FP32Vector3 result)
+        {
+            result.x1 = matrix.r1c1 * vector.x1 + matrix.r1c2 * vector.x2;
+            result.x2 = matrix.r2c1 * vector.x1 + matrix.r2c2 * vector.x2;
+            result.x3 = matrix.r3c1 * vector.x1 + matrix.r3c2 * vector.x2;
+        }
+
+        public static void GetLeftProduct(in FP32Vector3 vector, in FP32Matrix2x3 matrix, out FP32Vector2 result)
+        {
+            result.x1 = vector.x1 * matrix.r1c1 + vector.x2 * matrix.r2c1 + vector.x3 * matrix.r3c1;
+            result.x2 = vector.x1 * matrix.r1c2 + vector.x2 * matrix.r2c2 + vector.x3 * matrix.r3c2;
+        }
+    }
+}

BasicGeometry/FP32Matrix3x2.cs

@@ -0,0 +1,221 @@
+/*
+ * Copyright 2019-2025 Andrey Pokidov <andrey.pokidov@gmail.com>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 11 Nov 2024
+ */
+namespace BasicGeometry
+{
+    public struct FP32Matrix3x2
+    {
+        public float r1c1 = 0.0f, r1c2 = 0.0f, r1c3 = 0.0f;
+        public float r2c1 = 0.0f, r2c2 = 0.0f, r2c3 = 0.0f;
+
+        public FP32Matrix3x2(in FP32Matrix3x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+            this.r1c3 = matrix.r1c3;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+            this.r2c3 = matrix.r2c3;
+        }
+
+        public FP32Matrix3x2(in FP64Matrix3x2 matrix)
+        {
+            this.r1c1 = (float) matrix.r1c1;
+            this.r1c2 = (float) matrix.r1c2;
+            this.r1c3 = (float) matrix.r1c3;
+
+            this.r2c1 = (float) matrix.r2c1;
+            this.r2c2 = (float) matrix.r2c2;
+            this.r2c3 = (float) matrix.r2c3;
+        }
+
+        public FP32Matrix3x2(in FP32Matrix2x3 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r2c1;
+            this.r1c3 = matrix.r3c1;
+
+            this.r2c1 = matrix.r1c2;
+            this.r2c2 = matrix.r2c2;
+            this.r2c3 = matrix.r3c2;
+        }
+
+        public FP32Matrix3x2(in FP64Matrix2x3 matrix)
+        {
+            this.r1c1 = (float) matrix.r1c1;
+            this.r1c2 = (float) matrix.r2c1;
+            this.r1c3 = (float) matrix.r3c1;
+
+            this.r2c1 = (float) matrix.r1c2;
+            this.r2c2 = (float) matrix.r2c2;
+            this.r2c3 = (float) matrix.r3c2;
+        }
+
+        public void Reset()
+        {
+            this.r1c1 = 0.0f;
+            this.r1c2 = 0.0f;
+            this.r1c3 = 0.0f;
+
+            this.r2c1 = 0.0f;
+            this.r2c2 = 0.0f;
+            this.r2c3 = 0.0f;
+        }
+    
+        public void SetValues(in FP32Matrix3x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+            this.r1c3 = matrix.r1c3;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+            this.r2c3 = matrix.r2c3;
+        }
+
+        public void SetValues(in FP64Matrix3x2 matrix)
+        {
+            this.r1c1 = (float) matrix.r1c1;
+            this.r1c2 = (float) matrix.r1c2;
+            this.r1c3 = (float) matrix.r1c3;
+
+            this.r2c1 = (float) matrix.r2c1;
+            this.r2c2 = (float) matrix.r2c2;
+            this.r2c3 = (float) matrix.r2c3;
+        }
+    
+        public void SetTransposed(in FP32Matrix2x3 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r2c1;
+            this.r1c3 = matrix.r3c1;
+
+            this.r2c1 = matrix.r1c2;
+            this.r2c2 = matrix.r2c2;
+            this.r2c3 = matrix.r3c2;
+        }
+
+        public void SetTransposed(in FP64Matrix2x3 matrix)
+        {
+            this.r1c1 = (float) matrix.r1c1;
+            this.r1c2 = (float) matrix.r2c1;
+            this.r1c3 = (float) matrix.r3c1;
+
+            this.r2c1 = (float) matrix.r1c2;
+            this.r2c2 = (float) matrix.r2c2;
+            this.r2c3 = (float) matrix.r3c2;
+        }
+
+        public void SetRow1(float c1, float c2, float c3)
+        {
+            this.r1c1 = c1;
+            this.r1c2 = c2;
+            this.r1c3 = c3;
+        }
+
+        public void SetRow2(float c1, float c2, float c3)
+        {
+            this.r2c1 = c1;
+            this.r2c2 = c2;
+            this.r2c3 = c3;
+        }
+
+        public void SetColumn1(float r1, float r2)
+        {
+            this.r1c1 = r1;
+            this.r2c1 = r2;
+        }
+
+        public void SetColumn2(float r1, float r2)
+        {
+            this.r1c2 = r1;
+            this.r2c2 = r2;
+        }
+
+        public void SetColumn3(float r1, float r2)
+        {
+            this.r1c3 = r1;
+            this.r2c3 = r2;
+        }
+
+        public void AppendScaled(in FP32Matrix3x2 matrix, float scale)
+        {
+            this.r1c1 += matrix.r1c1 * scale;
+            this.r1c2 += matrix.r1c2 * scale;
+            this.r1c3 += matrix.r1c3 * scale;
+
+            this.r2c1 += matrix.r2c1 * scale;
+            this.r2c2 += matrix.r2c2 * scale;
+            this.r2c3 += matrix.r2c3 * scale;
+        }
+
+        public static void Add(in FP32Matrix3x2 matrix1, in FP32Matrix3x2 matrix2, out FP32Matrix3x2 sum)
+        {
+            sum.r1c1 = matrix1.r1c1 + matrix2.r1c1;
+            sum.r1c2 = matrix1.r1c2 + matrix2.r1c2;
+            sum.r1c3 = matrix1.r1c3 + matrix2.r1c3;
+
+            sum.r2c1 = matrix1.r2c1 + matrix2.r2c1;
+            sum.r2c2 = matrix1.r2c2 + matrix2.r2c2;
+            sum.r2c3 = matrix1.r2c3 + matrix2.r2c3;
+        }
+
+        public static void Subtract(in FP32Matrix3x2 minuend, in FP32Matrix3x2 subtrahend, out FP32Matrix3x2 difference)
+        {
+            difference.r1c1 = minuend.r1c1 - subtrahend.r1c1;
+            difference.r1c2 = minuend.r1c2 - subtrahend.r1c2;
+            difference.r1c3 = minuend.r1c3 - subtrahend.r1c3;
+
+            difference.r2c1 = minuend.r2c1 - subtrahend.r2c1;
+            difference.r2c2 = minuend.r2c2 - subtrahend.r2c2;
+            difference.r2c3 = minuend.r2c3 - subtrahend.r2c3;
+        }
+
+        public static void Multiply(in FP32Matrix3x2 multiplicand, float multiplier, out FP32Matrix3x2 product)
+        {
+            product.r1c1 = multiplicand.r1c1 * multiplier;
+            product.r1c2 = multiplicand.r1c2 * multiplier;
+            product.r1c3 = multiplicand.r1c3 * multiplier;
+
+            product.r2c1 = multiplicand.r2c1 * multiplier;
+            product.r2c2 = multiplicand.r2c2 * multiplier;
+            product.r2c3 = multiplicand.r2c3 * multiplier;
+        }
+
+        public static void Divide(in FP32Matrix3x2 dividend, float divisor, out FP32Matrix3x2 quotient)
+        {
+            Multiply(dividend, 1.0f / divisor, out quotient);
+        }
+
+        public static void GetRightProduct(in FP32Matrix3x2 matrix, in FP32Vector3 vector, out FP32Vector2 result)
+        {
+            result.x1 = matrix.r1c1 * vector.x1 + matrix.r1c2 * vector.x2 + matrix.r1c3 * vector.x3;
+            result.x2 = matrix.r2c1 * vector.x1 + matrix.r2c2 * vector.x2 + matrix.r2c3 * vector.x3;
+        }
+
+        public static void GetLeftProduct(in FP32Vector2 vector, in FP32Matrix3x2 matrix, out FP32Vector3 result)
+        {
+            result.x1 = vector.x1 * matrix.r1c1 + vector.x2 * matrix.r2c1;
+            result.x2 = vector.x1 * matrix.r1c2 + vector.x2 * matrix.r2c2;
+            result.x3 = vector.x1 * matrix.r1c3 + vector.x2 * matrix.r2c3;
+        }
+    }
+}

BasicGeometry/FP32Matrix3x3.cs

@@ -0,0 +1,438 @@
+/*
+ * Copyright 2019-2025 Andrey Pokidov <andrey.pokidov@gmail.com>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+using System;
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 10 Feb 2019
+ */
+
+namespace BasicGeometry
+{
+    public struct FP32Matrix3x3
+    {
+        public float r1c1 = 0.0f, r1c2 = 0.0f, r1c3 = 0.0f;
+        public float r2c1 = 0.0f, r2c2 = 0.0f, r2c3 = 0.0f;
+        public float r3c1 = 0.0f, r3c2 = 0.0f, r3c3 = 0.0f;
+
+        public FP32Matrix3x3(float d1, float d2, float d3)
+        {
+            this.r1c1 = d1;
+            this.r2c2 = d2;
+            this.r3c3 = d3;
+        }
+
+        public FP32Matrix3x3(in FP32Matrix3x3 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+            this.r1c3 = matrix.r1c3;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+            this.r2c3 = matrix.r2c3;
+
+            this.r3c1 = matrix.r3c1;
+            this.r3c2 = matrix.r3c2;
+            this.r3c3 = matrix.r3c3;
+        }
+
+        public FP32Matrix3x3(in FP64Matrix3x3 matrix)
+        {
+            this.r1c1 = (float)matrix.r1c1;
+            this.r1c2 = (float)matrix.r1c2;
+            this.r1c3 = (float)matrix.r1c3;
+
+            this.r2c1 = (float)matrix.r2c1;
+            this.r2c2 = (float)matrix.r2c2;
+            this.r2c3 = (float)matrix.r2c3;
+
+            this.r3c1 = (float)matrix.r3c1;
+            this.r3c2 = (float)matrix.r3c2;
+            this.r3c3 = (float)matrix.r3c3;
+        }
+
+        public readonly float GetDeterminant()
+        {
+            return this.r1c1 * (this.r2c2 * this.r3c3 - this.r2c3 * this.r3c2)
+                 + this.r1c2 * (this.r2c3 * this.r3c1 - this.r2c1 * this.r3c3)
+                 + this.r1c3 * (this.r2c1 * this.r3c2 - this.r2c2 * this.r3c1);
+        }
+
+        public readonly bool IsSingular()
+        {
+            float determinant = this.GetDeterminant();
+            return -FP32Utility.EPSYLON <= determinant && determinant <= FP32Utility.EPSYLON;
+        }
+
+        public void Transpose()
+        {
+            (this.r1c2, this.r2c1) = (this.r2c1, this.r1c2);
+            (this.r1c3, this.r3c1) = (this.r3c1, this.r1c3);
+            (this.r2c3, this.r3c2) = (this.r3c2, this.r2c3);
+        }
+
+        public bool Invert()
+        {
+            float determinant = this.GetDeterminant();
+    
+            if (-FP32Utility.EPSYLON <= determinant && determinant <= FP32Utility.EPSYLON) {
+                return false;
+            }
+
+            float r1c1 = this.r2c2 * this.r3c3 - this.r2c3 * this.r3c2;
+            float r1c2 = this.r1c3 * this.r3c2 - this.r1c2 * this.r3c3;
+            float r1c3 = this.r1c2 * this.r2c3 - this.r1c3 * this.r2c2;
+
+            float r2c1 = this.r2c3 * this.r3c1 - this.r2c1 * this.r3c3;
+            float r2c2 = this.r1c1 * this.r3c3 - this.r1c3 * this.r3c1;
+            float r2c3 = this.r1c3 * this.r2c1 - this.r1c1 * this.r2c3;
+
+            float r3c1 = this.r2c1 * this.r3c2 - this.r2c2 * this.r3c1;
+            float r3c2 = this.r1c2 * this.r3c1 - this.r1c1 * this.r3c2;
+            float r3c3 = this.r1c1 * this.r2c2 - this.r1c2 * this.r2c1;
+
+            float mutiplier = 1.0f / determinant;
+
+            this.r1c1 = r1c1 * mutiplier;
+            this.r1c2 = r1c2 * mutiplier;
+            this.r1c3 = r1c3 * mutiplier;
+
+            this.r2c1 = r2c1 * mutiplier;
+            this.r2c2 = r2c2 * mutiplier;
+            this.r2c3 = r2c3 * mutiplier;
+
+            this.r3c1 = r3c1 * mutiplier;
+            this.r3c2 = r3c2 * mutiplier;
+            this.r3c3 = r3c3 * mutiplier;
+
+            return true;
+        }
+
+        public void Reset()
+        {
+            this.r1c1 = 0.0f;
+            this.r1c2 = 0.0f;
+            this.r1c3 = 0.0f;
+
+            this.r2c1 = 0.0f;
+            this.r2c2 = 0.0f;
+            this.r2c3 = 0.0f;
+
+            this.r3c1 = 0.0f;
+            this.r3c2 = 0.0f;
+            this.r3c3 = 0.0f;
+        }
+
+        public void SetToIdentity()
+        {
+            this.r1c1 = 1.0f;
+            this.r1c2 = 0.0f;
+            this.r1c3 = 0.0f;
+
+            this.r2c1 = 0.0f;
+            this.r2c2 = 1.0f;
+            this.r2c3 = 0.0f;
+
+            this.r3c1 = 0.0f;
+            this.r3c2 = 0.0f;
+            this.r3c3 = 1.0f;
+        }
+
+        public void SetToDiagonal(float d1, float d2, float d3)
+        {
+            this.r1c1 = d1;
+            this.r1c2 = 0.0f;
+            this.r1c3 = 0.0f;
+
+            this.r2c1 = 0.0f;
+            this.r2c2 = d2;
+            this.r2c3 = 0.0f;
+
+            this.r2c1 = 0.0f;
+            this.r2c2 = 0.0f;
+            this.r2c3 = d3;
+        }
+
+        public void SetValues(in FP32Matrix3x3 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+            this.r1c3 = matrix.r1c3;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+            this.r2c3 = matrix.r2c3;
+
+            this.r3c1 = matrix.r3c1;
+            this.r3c2 = matrix.r3c2;
+            this.r3c3 = matrix.r3c3;
+        }
+
+        public void SetValues(in FP64Matrix3x3 matrix)
+        {
+            this.r1c1 = (float)matrix.r1c1;
+            this.r1c2 = (float)matrix.r1c2;
+            this.r1c3 = (float)matrix.r1c3;
+
+            this.r2c1 = (float)matrix.r2c1;
+            this.r2c2 = (float)matrix.r2c2;
+            this.r2c3 = (float)matrix.r2c3;
+
+            this.r3c1 = (float)matrix.r3c1;
+            this.r3c2 = (float)matrix.r3c2;
+            this.r3c3 = (float)matrix.r3c3;
+        }
+
+        public void SetTransposedOf(in FP32Matrix3x3 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r2c2 = matrix.r2c2;
+            this.r3c3 = matrix.r3c3;
+
+            (this.r1c2, this.r2c1) = (matrix.r2c1, matrix.r1c2);
+            (this.r1c3, this.r3c1) = (matrix.r3c1, matrix.r1c3);
+            (this.r2c3, this.r3c2) = (matrix.r3c2, matrix.r2c3);
+        }
+
+        public void SetTransposedOf(in FP64Matrix3x3 matrix)
+        {
+            this.r1c1 = (float)matrix.r1c1;
+            this.r1c2 = (float)matrix.r2c1;
+            this.r1c3 = (float)matrix.r3c1;
+
+            this.r2c1 = (float)matrix.r1c2;
+            this.r2c2 = (float)matrix.r2c2;
+            this.r2c3 = (float)matrix.r3c2;
+
+            this.r3c1 = (float)matrix.r1c3;
+            this.r3c2 = (float)matrix.r2c3;
+            this.r3c3 = (float)matrix.r3c3;
+        }
+
+        public bool SetInvertedOf(in FP32Matrix3x3 matrix)
+        {
+            float determinant = matrix.GetDeterminant();
+    
+            if (-FP32Utility.EPSYLON <= determinant && determinant <= FP32Utility.EPSYLON) {
+                return false;
+            }
+
+            float r1c1 = matrix.r2c2 * matrix.r3c3 - matrix.r2c3 * matrix.r3c2;
+            float r1c2 = matrix.r1c3 * matrix.r3c2 - matrix.r1c2 * matrix.r3c3;
+            float r1c3 = matrix.r1c2 * matrix.r2c3 - matrix.r1c3 * matrix.r2c2;
+
+            float r2c1 = matrix.r2c3 * matrix.r3c1 - matrix.r2c1 * matrix.r3c3;
+            float r2c2 = matrix.r1c1 * matrix.r3c3 - matrix.r1c3 * matrix.r3c1;
+            float r2c3 = matrix.r1c3 * matrix.r2c1 - matrix.r1c1 * matrix.r2c3;
+
+            float r3c1 = matrix.r2c1 * matrix.r3c2 - matrix.r2c2 * matrix.r3c1;
+            float r3c2 = matrix.r1c2 * matrix.r3c1 - matrix.r1c1 * matrix.r3c2;
+            float r3c3 = matrix.r1c1 * matrix.r2c2 - matrix.r1c2 * matrix.r2c1;
+
+            float mutiplier = 1.0f / determinant;
+
+            this.r1c1 = r1c1 * mutiplier;
+            this.r1c2 = r1c2 * mutiplier;
+            this.r1c3 = r1c3 * mutiplier;
+
+            this.r2c1 = r2c1 * mutiplier;
+            this.r2c2 = r2c2 * mutiplier;
+            this.r2c3 = r2c3 * mutiplier;
+
+            this.r3c1 = r3c1 * mutiplier;
+            this.r3c2 = r3c2 * mutiplier;
+            this.r3c3 = r3c3 * mutiplier;
+
+            return true;
+        }
+
+        public void SetRow1(float c1, float c2, float c3)
+        {
+            this.r1c1 = c1;
+            this.r1c2 = c2;
+            this.r1c3 = c3;
+        }
+
+        public void SetRow2(float c1, float c2, float c3)
+        {
+            this.r2c1 = c1;
+            this.r2c2 = c2;
+            this.r2c3 = c3;
+        }
+
+        public void SetRow3(float c1, float c2, float c3)
+        {
+            this.r3c1 = c1;
+            this.r3c2 = c2;
+            this.r3c3 = c3;
+        }
+
+        public void SetColumn1(float r1, float r2, float r3)
+        {
+            this.r1c1 = r1;
+            this.r2c1 = r2;
+            this.r3c1 = r3;
+        }
+
+        public void SetColumn2(float r1, float r2, float r3)
+        {
+            this.r1c2 = r1;
+            this.r2c2 = r2;
+            this.r3c2 = r3;
+        }
+
+        public void SetColumn3(float r1, float r2, float r3)
+        {
+            this.r1c3 = r1;
+            this.r2c3 = r2;
+            this.r3c3 = r3;
+        }
+
+        public void AppendScaled(in FP32Matrix3x3 matrix, float scale)
+        {
+            this.r1c1 += matrix.r1c1* scale;
+            this.r1c2 += matrix.r1c2* scale;
+            this.r1c3 += matrix.r1c3* scale;
+
+            this.r2c1 += matrix.r2c1* scale;
+            this.r2c2 += matrix.r2c2* scale;
+            this.r2c3 += matrix.r2c3* scale;
+
+            this.r3c1 += matrix.r3c1* scale;
+            this.r3c2 += matrix.r3c2* scale;
+            this.r3c3 += matrix.r3c3* scale;
+        }
+
+        public static void Add(in FP32Matrix3x3 matrix1, in FP32Matrix3x3 matrix2, out FP32Matrix3x3 sum)
+        {
+            sum.r1c1 = matrix1.r1c1 + matrix2.r1c1;
+            sum.r1c2 = matrix1.r1c2 + matrix2.r1c2;
+            sum.r1c3 = matrix1.r1c3 + matrix2.r1c3;
+
+            sum.r2c1 = matrix1.r2c1 + matrix2.r2c1;
+            sum.r2c2 = matrix1.r2c2 + matrix2.r2c2;
+            sum.r2c3 = matrix1.r2c3 + matrix2.r2c3;
+
+            sum.r3c1 = matrix1.r3c1 + matrix2.r3c1;
+            sum.r3c2 = matrix1.r3c2 + matrix2.r3c2;
+            sum.r3c3 = matrix1.r3c3 + matrix2.r3c3;
+        }
+
+        public static void Subtract(in FP32Matrix3x3 minuend, in FP32Matrix3x3 subtrahend, out FP32Matrix3x3 difference)
+        {
+            difference.r1c1 = minuend.r1c1 - subtrahend.r1c1;
+            difference.r1c2 = minuend.r1c2 - subtrahend.r1c2;
+            difference.r1c3 = minuend.r1c3 - subtrahend.r1c3;
+
+            difference.r2c1 = minuend.r2c1 - subtrahend.r2c1;
+            difference.r2c2 = minuend.r2c2 - subtrahend.r2c2;
+            difference.r2c3 = minuend.r2c3 - subtrahend.r2c3;
+
+            difference.r3c1 = minuend.r3c1 - subtrahend.r3c1;
+            difference.r3c2 = minuend.r3c2 - subtrahend.r3c2;
+            difference.r3c3 = minuend.r3c3 - subtrahend.r3c3;
+        }
+
+        public static void Multiply(in FP32Matrix3x3 multiplicand, float multiplier, out FP32Matrix3x3 product)
+        {
+            product.r1c1 = multiplicand.r1c1 * multiplier;
+            product.r1c2 = multiplicand.r1c2 * multiplier;
+            product.r1c3 = multiplicand.r1c3 * multiplier;
+
+            product.r2c1 = multiplicand.r2c1 * multiplier;
+            product.r2c2 = multiplicand.r2c2 * multiplier;
+            product.r2c3 = multiplicand.r2c3 * multiplier;
+
+            product.r3c1 = multiplicand.r3c1 * multiplier;
+            product.r3c2 = multiplicand.r3c2 * multiplier;
+            product.r3c3 = multiplicand.r3c3 * multiplier;
+        }
+
+        public static void Divide(in FP32Matrix3x3 dividend, float divisor, out FP32Matrix3x3 quotient)
+        {
+            Multiply(dividend, 1.0f / divisor, out quotient);
+        }
+
+        public static void GetRightProduct(in FP32Matrix3x3 matrix, in FP32Vector3 vector, out FP32Vector3 result)
+        {
+            float x1 = matrix.r1c1 * vector.x1 + matrix.r1c2 * vector.x2 + matrix.r1c3 * vector.x3;
+            float x2 = matrix.r2c1 * vector.x1 + matrix.r2c2 * vector.x2 + matrix.r2c3 * vector.x3;
+            float x3 = matrix.r3c1 * vector.x1 + matrix.r3c2 * vector.x2 + matrix.r3c3 * vector.x3;
+
+            result.x1 = x1;
+            result.x2 = x2;
+            result.x3 = x3;
+        }
+
+        public static void GetLeftProduct(in FP32Vector3 vector, in FP32Matrix3x3 matrix, out FP32Vector3 result)
+        {
+            float x1 = vector.x1 * matrix.r1c1 + vector.x2 * matrix.r2c1 + vector.x3 * matrix.r3c1;
+            float x2 = vector.x1 * matrix.r1c2 + vector.x2 * matrix.r2c2 + vector.x3 * matrix.r3c2;
+            float x3 = vector.x1 * matrix.r1c3 + vector.x2 * matrix.r2c3 + vector.x3 * matrix.r3c3;
+
+            result.x1 = x1;
+            result.x2 = x2;
+            result.x3 = x3;
+        }
+
+        public static void LoadZero(out FP32Matrix3x3 matrix)
+        {
+            matrix.r1c1 = 0.0f;
+            matrix.r1c2 = 0.0f;
+            matrix.r1c3 = 0.0f;
+
+            matrix.r2c1 = 0.0f;
+            matrix.r2c2 = 0.0f;
+            matrix.r2c3 = 0.0f;
+
+            matrix.r3c1 = 0.0f;
+            matrix.r3c2 = 0.0f;
+            matrix.r3c3 = 0.0f;
+        }
+
+        public static void LoadIdentity(out FP32Matrix3x3 matrix)
+        {
+            matrix.r1c1 = 1.0f;
+            matrix.r1c2 = 0.0f;
+            matrix.r1c3 = 0.0f;
+
+            matrix.r2c1 = 0.0f;
+            matrix.r2c2 = 1.0f;
+            matrix.r2c3 = 0.0f;
+
+            matrix.r3c1 = 0.0f;
+            matrix.r3c2 = 0.0f;
+            matrix.r3c3 = 1.0f;
+        }
+
+        public static void LoadDiagonal(float d1, float d2, float d3, out FP32Matrix3x3 matrix)
+        {
+            matrix.r1c1 = d1;
+            matrix.r1c2 = 0.0f;
+            matrix.r1c3 = 0.0f;
+
+            matrix.r2c1 = 0.0f;
+            matrix.r2c2 = d2;
+            matrix.r2c3 = 0.0f;
+
+            matrix.r3c1 = 0.0f;
+            matrix.r3c2 = 0.0f;
+            matrix.r3c3 = d3;
+        }
+    }
+}

BasicGeometry/FP32MatrixProduct.cs

@@ -0,0 +1,173 @@
+/*
+ * Copyright 2019-2025 Andrey Pokidov <andrey.pokidov@gmail.com>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+using System;
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 11 Nov 2024
+ */
+namespace BasicGeometry
+{
+    public class FP32MatrixProduct
+    {
+        public static void Get2x2At2x2(in FP32Matrix2x2 left, in FP32Matrix2x2 right, out FP32Matrix2x2 product)
+        {
+            float r1c1 = left.r1c1 * right.r1c1 + left.r1c2 * right.r2c1;
+            float r1c2 = left.r1c1 * right.r1c2 + left.r1c2 * right.r2c2;
+
+            float r2c1 = left.r2c1 * right.r1c1 + left.r2c2 * right.r2c1;
+            float r2c2 = left.r2c1 * right.r1c2 + left.r2c2 * right.r2c2;
+
+            product.r1c1 = r1c1;
+            product.r1c2 = r1c2;
+
+            product.r2c1 = r2c1;
+            product.r2c2 = r2c2;
+        }
+
+        public static void Get2x2At3x2(in FP32Matrix2x2 left, in FP32Matrix3x2 right, out FP32Matrix3x2 product)
+        {
+            float r1c1 = left.r1c1 * right.r1c1 + left.r1c2 * right.r2c1;
+            float r1c2 = left.r1c1 * right.r1c2 + left.r1c2 * right.r2c2;
+            float r1c3 = left.r1c1 * right.r1c3 + left.r1c2 * right.r2c3;
+
+            float r2c1 = left.r2c1 * right.r1c1 + left.r2c2 * right.r2c1;
+            float r2c2 = left.r2c1 * right.r1c2 + left.r2c2 * right.r2c2;
+            float r2c3 = left.r2c1 * right.r1c3 + left.r2c2 * right.r2c3;
+
+            product.r1c1 = r1c1;
+            product.r1c2 = r1c2;
+            product.r1c3 = r1c3;
+
+            product.r2c1 = r2c1;
+            product.r2c2 = r2c2;
+            product.r2c3 = r2c3;
+        }
+    
+        public static void Get2x3At2x2(in FP32Matrix2x3 left, in FP32Matrix2x2 right, out FP32Matrix2x3 product)
+        {
+            float r1c1 = left.r1c1 * right.r1c1 + left.r1c2 * right.r2c1;
+            float r1c2 = left.r1c1 * right.r1c2 + left.r1c2 * right.r2c2;
+        
+            float r2c1 = left.r2c1 * right.r1c1 + left.r2c2 * right.r2c1;
+            float r2c2 = left.r2c1 * right.r1c2 + left.r2c2 * right.r2c2;
+        
+            float r3c1 = left.r3c1 * right.r1c1 + left.r3c2 * right.r2c1;
+            float r3c2 = left.r3c1 * right.r1c2 + left.r3c2 * right.r2c2;
+        
+            product.r1c1 = r1c1;
+            product.r1c2 = r1c2;
+        
+            product.r2c1 = r2c1;
+            product.r2c2 = r2c2;
+        
+            product.r3c1 = r3c1;
+            product.r3c2 = r3c2;
+        }
+
+        public static void Get2x3At3x2(in FP32Matrix2x3 left, in FP32Matrix3x2 right, out FP32Matrix3x3 product)
+        {
+            product.r1c1 = left.r1c1 * right.r1c1 + left.r1c2 * right.r2c1;
+            product.r1c2 = left.r1c1 * right.r1c2 + left.r1c2 * right.r2c2;
+            product.r1c3 = left.r1c1 * right.r1c3 + left.r1c2 * right.r2c3;
+
+            product.r2c1 = left.r2c1 * right.r1c1 + left.r2c2 * right.r2c1;
+            product.r2c2 = left.r2c1 * right.r1c2 + left.r2c2 * right.r2c2;
+            product.r2c3 = left.r2c1 * right.r1c3 + left.r2c2 * right.r2c3;
+
+            product.r3c1 = left.r3c1 * right.r1c1 + left.r3c2 * right.r2c1;
+            product.r3c2 = left.r3c1 * right.r1c2 + left.r3c2 * right.r2c2;
+            product.r3c3 = left.r3c1 * right.r1c2 + left.r3c2 * right.r2c3;
+        }
+
+        public static void Get3x2At3x3(in FP32Matrix3x2 left, in FP32Matrix3x3 right, out FP32Matrix3x2 product)
+        {
+            float r1c1 = left.r1c1 * right.r1c1 + left.r1c2 * right.r2c1 + left.r1c3 * right.r3c1;
+            float r1c2 = left.r1c1 * right.r1c2 + left.r1c2 * right.r2c2 + left.r1c3 * right.r3c2;
+            float r1c3 = left.r1c1 * right.r1c3 + left.r1c2 * right.r2c3 + left.r1c3 * right.r3c3;
+
+            float r2c1 = left.r2c1 * right.r1c1 + left.r2c2 * right.r2c1 + left.r2c3 * right.r3c1;
+            float r2c2 = left.r2c1 * right.r1c2 + left.r2c2 * right.r2c2 + left.r2c3 * right.r3c2;
+            float r2c3 = left.r2c1 * right.r1c3 + left.r2c2 * right.r2c3 + left.r2c3 * right.r3c3;
+
+            product.r1c1 = r1c1;
+            product.r1c2 = r1c2;
+            product.r1c3 = r1c3;
+
+            product.r2c1 = r2c1;
+            product.r2c2 = r2c2;
+            product.r2c3 = r2c3;
+        }
+
+        public static void Get3x2At2x3(in FP32Matrix3x2 left, in FP32Matrix2x3 right, out FP32Matrix2x2 product)
+        {
+            product.r1c1 = left.r1c1 * right.r1c1 + left.r1c2 * right.r2c1 + left.r1c3 * right.r3c1;
+            product.r1c2 = left.r1c1 * right.r1c2 + left.r1c2 * right.r2c2 + left.r1c3 * right.r3c2;
+
+            product.r2c1 = left.r2c1 * right.r1c1 + left.r2c2 * right.r2c1 + left.r2c3 * right.r3c1;
+            product.r2c2 = left.r2c1 * right.r1c2 + left.r2c2 * right.r2c2 + left.r2c3 * right.r3c2;
+        }
+
+        public static void Get3x3At2x3(in FP32Matrix3x3 left, in FP32Matrix2x3 right, out FP32Matrix2x3 product)
+        {
+            float r1c1 = left.r1c1 * right.r1c1 + left.r1c2 * right.r2c1 + left.r1c3 * right.r3c1;
+            float r1c2 = left.r1c1 * right.r1c2 + left.r1c2 * right.r2c2 + left.r1c3 * right.r3c2;
+
+            float r2c1 = left.r2c1 * right.r1c1 + left.r2c2 * right.r2c1 + left.r2c3 * right.r3c1;
+            float r2c2 = left.r2c1 * right.r1c2 + left.r2c2 * right.r2c2 + left.r2c3 * right.r3c2;
+
+            float r3c1 = left.r3c1 * right.r1c1 + left.r3c2 * right.r2c1 + left.r3c3 * right.r3c1;
+            float r3c2 = left.r3c1 * right.r1c2 + left.r3c2 * right.r2c2 + left.r3c3 * right.r3c2;
+
+            product.r1c1 = r1c1;
+            product.r1c2 = r1c2;
+
+            product.r2c1 = r2c1;
+            product.r2c2 = r2c2;
+
+            product.r3c1 = r3c1;
+            product.r3c2 = r3c2;
+        }
+
+        public static void Get3x3At3x3(in FP32Matrix3x3 left, in FP32Matrix3x3 right, out FP32Matrix3x3 product)
+        {
+            float r1c1 = left.r1c1 * right.r1c1 + left.r1c2 * right.r2c1 + left.r1c3 * right.r3c1;
+            float r1c2 = left.r1c1 * right.r1c2 + left.r1c2 * right.r2c2 + left.r1c3 * right.r3c2;
+            float r1c3 = left.r1c1 * right.r1c3 + left.r1c2 * right.r2c3 + left.r1c3 * right.r3c3;
+
+            float r2c1 = left.r2c1 * right.r1c1 + left.r2c2 * right.r2c1 + left.r2c3 * right.r3c1;
+            float r2c2 = left.r2c1 * right.r1c2 + left.r2c2 * right.r2c2 + left.r2c3 * right.r3c2;
+            float r2c3 = left.r2c1 * right.r1c3 + left.r2c2 * right.r2c3 + left.r2c3 * right.r3c3;
+
+            float r3c1 = left.r3c1 * right.r1c1 + left.r3c2 * right.r2c1 + left.r3c3 * right.r3c1;
+            float r3c2 = left.r3c1 * right.r1c2 + left.r3c2 * right.r2c2 + left.r3c3 * right.r3c2;
+            float r3c3 = left.r3c1 * right.r1c3 + left.r3c2 * right.r2c3 + left.r3c3 * right.r3c3;
+
+            product.r1c1 = r1c1;
+            product.r1c2 = r1c2;
+            product.r1c3 = r1c3;
+
+            product.r2c1 = r2c1;
+            product.r2c2 = r2c2;
+            product.r2c3 = r2c3;
+
+            product.r3c1 = r3c1;
+            product.r3c2 = r3c2;
+            product.r3c3 = r3c3;
+        }
+    }
+}

BasicGeometry/FP32Quaternion.cs

@@ -0,0 +1,203 @@
+using System;
+
+namespace BasicGeometry
+{
+    public struct FP32Quaternion
+    {
+        public float s0 = 0.0f, x1 = 0.0f, x2 = 0.0f, x3 = 0.0f;
+
+        public FP32Quaternion(float s0, float x1, float x2, float x3)
+        {
+            this.s0 = s0;
+            this.x1 = x1;
+            this.x2 = x2;
+            this.x3 = x3;
+        }
+
+        public FP32Quaternion(in FP32Quaternion quaternion)
+        {
+            this.s0 = quaternion.s0;
+            this.x1 = quaternion.x1;
+            this.x2 = quaternion.x2;
+            this.x3 = quaternion.x3;
+        }
+
+        public FP32Quaternion(in FP64Quaternion quaternion)
+        {
+            this.s0 = (float)quaternion.s0;
+            this.x1 = (float)quaternion.x1;
+            this.x2 = (float)quaternion.x2;
+            this.x3 = (float)quaternion.x3;
+        }
+
+        public void Reset()
+        {
+            this.s0 = 0.0f;
+            this.x1 = 0.0f;
+            this.x2 = 0.0f;
+            this.x3 = 0.0f;
+        }
+
+        public void Conjugate()
+        {
+            this.x1 = -this.x1;
+            this.x2 = -this.x2;
+            this.x3 = -this.x3;
+        }
+
+        public void SetValues(float s0, float x1, float x2, float x3)
+        {
+            this.s0 = s0;
+            this.x1 = x1;
+            this.x2 = x2;
+            this.x3 = x3;
+        }
+
+        public void SetValues(in FP32Quaternion quaternion)
+        {
+            this.s0 = quaternion.s0;
+            this.x1 = quaternion.x1;
+            this.x2 = quaternion.x2;
+            this.x3 = quaternion.x3;
+        }
+
+        public void SetValues(in FP64Quaternion quaternion)
+        {
+            this.s0 = (float)quaternion.s0;
+            this.x1 = (float)quaternion.x1;
+            this.x2 = (float)quaternion.x2;
+            this.x3 = (float)quaternion.x3;
+        }
+
+        public void SetConjugateOf(in FP32Quaternion quaternion)
+        {
+            this.s0 = quaternion.s0;
+            this.x1 = -quaternion.x1;
+            this.x2 = -quaternion.x2;
+            this.x3 = -quaternion.x3;
+        }
+
+        public readonly void MakeRotationMatrix(out FP32Matrix3x3 matrix)
+        {
+            float s0s0 = this.s0 * this.s0;
+            float x1x1 = this.x1 * this.x1;
+            float x2x2 = this.x2 * this.x2;
+            float x3x3 = this.x3 * this.x3;
+
+            float squareModule = (s0s0 + x1x1) + (x2x2 + x3x3);
+
+            if (-FP32Utility.EPSYLON <= squareModule && squareModule <= FP32Utility.EPSYLON)
+            {
+                FP32Matrix3x3.LoadIdentity(out matrix);
+                return;
+            }
+
+            float corrector1;
+            float corrector2;
+
+            if (1.0f - FP32Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0f + FP32Utility.TWO_EPSYLON) {
+                corrector1 = 2.0f - squareModule;
+                corrector2 = 2.0f * corrector1;
+            }
+            else {
+                corrector1 = 1.0f / squareModule;
+                corrector2 = 2.0f / squareModule;
+            }
+
+            float s0x1 = this.s0 * this.x1;
+            float s0x2 = this.s0 * this.x2;
+            float s0x3 = this.s0 * this.x3;
+            float x1x2 = this.x1 * this.x2;
+            float x1x3 = this.x1 * this.x3;
+            float x2x3 = this.x2 * this.x3;
+
+            matrix.r1c1 = corrector1 * ((s0s0 + x1x1) - (x2x2 + x3x3));
+            matrix.r2c2 = corrector1 * ((s0s0 + x2x2) - (x1x1 + x3x3));
+            matrix.r3c3 = corrector1 * ((s0s0 + x3x3) - (x1x1 + x2x2));
+
+            matrix.r1c2 = corrector2 * (x1x2 - s0x3);
+            matrix.r2c3 = corrector2 * (x2x3 - s0x1);
+            matrix.r3c1 = corrector2 * (x1x3 - s0x2);
+
+            matrix.r2c1 = corrector2 * (x1x2 + s0x3);
+            matrix.r3c2 = corrector2 * (x2x3 + s0x1);
+            matrix.r1c3 = corrector2 * (x1x3 + s0x2);
+        }
+
+        public readonly void MakeReverseMatrix(out FP32Matrix3x3 matrix)
+        {
+            float s0s0 = this.s0 * this.s0;
+            float x1x1 = this.x1 * this.x1;
+            float x2x2 = this.x2 * this.x2;
+            float x3x3 = this.x3 * this.x3;
+
+            float squareModule = (s0s0 + x1x1) + (x2x2 + x3x3);
+
+            if (-FP32Utility.EPSYLON <= squareModule && squareModule <= FP32Utility.EPSYLON)
+            {
+                FP32Matrix3x3.LoadIdentity(out matrix);
+                return;
+            }
+
+            float corrector1;
+            float corrector2;
+
+            if (1.0f - FP32Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0f + FP32Utility.TWO_EPSYLON) {
+                corrector1 = 2.0f - squareModule;
+                corrector2 = 2.0f * corrector1;
+            }
+            else {
+                corrector1 = 1.0f / squareModule;
+                corrector2 = 2.0f / squareModule;
+            }
+
+            float s0x1 = this.s0 * this.x1;
+            float s0x2 = this.s0 * this.x2;
+            float s0x3 = this.s0 * this.x3;
+            float x1x2 = this.x1 * this.x2;
+            float x1x3 = this.x1 * this.x3;
+            float x2x3 = this.x2 * this.x3;
+
+            matrix.r1c1 = corrector1 * ((s0s0 + x1x1) - (x2x2 + x3x3));
+            matrix.r2c2 = corrector1 * ((s0s0 + x2x2) - (x1x1 + x3x3));
+            matrix.r3c3 = corrector1 * ((s0s0 + x3x3) - (x1x1 + x2x2));
+
+            matrix.r1c2 = corrector2 * (x1x2 + s0x3);
+            matrix.r2c3 = corrector2 * (x2x3 + s0x1);
+            matrix.r3c1 = corrector2 * (x1x3 + s0x2);
+
+            matrix.r2c1 = corrector2 * (x1x2 - s0x3);
+            matrix.r3c2 = corrector2 * (x2x3 - s0x1);
+            matrix.r1c3 = corrector2 * (x1x3 - s0x2);
+        }
+
+        public static void Add(in FP32Quaternion quaternion1, in FP32Quaternion quaternion2, out FP32Quaternion sum)
+        {
+            sum.s0 = quaternion1.s0 + quaternion2.s0;
+            sum.x1 = quaternion1.x1 + quaternion2.x1;
+            sum.x2 = quaternion1.x2 + quaternion2.x2;
+            sum.x3 = quaternion1.x3 + quaternion2.x3;
+        }
+
+        public static void Subtract(in FP32Quaternion minuend, in FP32Quaternion subtrahend, out FP32Quaternion difference)
+        {
+            difference.s0 = minuend.s0 - subtrahend.s0;
+            difference.x1 = minuend.x1 - subtrahend.x1;
+            difference.x2 = minuend.x2 - subtrahend.x2;
+            difference.x3 = minuend.x3 - subtrahend.x3;
+        }
+
+        public static void Multiply(in FP32Quaternion left, in FP32Quaternion right, out FP32Quaternion product)
+        {
+            float s0 = (left.s0 * right.s0 - left.x1 * right.x1) - (left.x2 * right.x2 + left.x3 * right.x3);
+            float x1 = (left.x1 * right.s0 + left.s0 * right.x1) - (left.x3 * right.x2 - left.x2 * right.x3);
+            float x2 = (left.x2 * right.s0 + left.s0 * right.x2) - (left.x1 * right.x3 - left.x3 * right.x1);
+            float x3 = (left.x3 * right.s0 + left.s0 * right.x3) - (left.x2 * right.x1 - left.x1 * right.x2);
+
+            product.s0 = s0;
+            product.x1 = x1;
+            product.x2 = x2;
+            product.x3 = x3;
+        }
+    }
+}

BasicGeometry/FP32Radians.cs

@@ -0,0 +1,76 @@
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 18 Nov 2024
+ */
+
+namespace BasicGeometry
+{
+    public class FP32Radians
+    {
+        public const float PI = 3.1415926536f;
+        public const float TWO_PI = 6.2831853072f;
+        public const float HALF_OF_PI = 1.5707963268f;
+        public const float THIRD_OF_PI = 1.0471975512f;
+        public const float FOURTH_OF_PI = 0.7853981634f;
+        public const float SIXTH_OF_PI = 0.5235987756f;
+
+        public const float DEGREES_IN_RADIAN = 57.295779513f;
+        public const float TURNS_IN_RADIAN = 0.1591549431f;
+
+        public static float ToDegrees(float radians)
+        {
+            return radians * DEGREES_IN_RADIAN;
+        }
+
+        public static float ToTurns(float radians)
+        {
+            return radians * TURNS_IN_RADIAN;
+        }
+
+        public static float ToUnits(float radians, AngleUnit toUnit)
+        {
+            if (toUnit == AngleUnit.DEGREES)
+            {
+                return radians * DEGREES_IN_RADIAN;
+            }
+
+            if (toUnit == AngleUnit.TURNS)
+            {
+                return radians * TURNS_IN_RADIAN;
+            }
+
+            return radians;
+        }
+
+        public static float Normalize(float radians, AngleRange range)
+        {
+            if (range == AngleRange.UNSIGNED_RANGE)
+            {
+                if (0.0 <= radians && radians < TWO_PI)
+                {
+                    return radians;
+                }
+
+            }
+            else
+            {
+                if (-PI < radians && radians <= PI)
+                {
+                    return radians;
+                }
+            }
+
+            float turns = radians * TURNS_IN_RADIAN;
+
+            turns -= MathF.Floor(turns);
+
+            if (range == AngleRange.SIGNED_RANGE && turns > 0.5f)
+            {
+                turns -= 1.0f;
+            }
+
+            return turns * TWO_PI;
+        }
+    }
+}

BasicGeometry/FP32Rotation3.cs

@@ -0,0 +1,44 @@
+/*
+ * Copyright 2019-2025 Andrey Pokidov <andrey.pokidov@gmail.com>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+using System;
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 2 Feb 2019
+ */
+
+namespace BasicGeometry
+{
+    public struct FP32Rotation3
+    {
+        private float angle = 0.0f;
+
+        private FP32Vector3 axis;
+
+        public FP32Rotation3(FP32Rotation3 rotation)
+        {
+            this.angle = rotation.angle;
+            this.axis = rotation.axis;
+        }
+
+        public void Reset()
+        {
+            this.angle = 0.0f;
+            this.axis.Reset();
+        }
+    }
+}

BasicGeometry/FP32Turns.cs

@@ -0,0 +1,63 @@
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 18 Nov 2024
+ */
+
+namespace BasicGeometry
+{
+    public class FP32Turns
+    {
+        public static float ToRadians(float turns)
+        {
+            return turns * FP32Radians.TWO_PI;
+        }
+
+        public static float ToDegrees(float turns)
+        {
+            return turns * 360.0f;
+        }
+
+        public static float ToUnits(float turns, AngleUnit toUnit)
+        {
+            if (toUnit == AngleUnit.RADIANS)
+            {
+                return turns * FP32Radians.TWO_PI;
+            }
+
+            if (toUnit == AngleUnit.DEGREES)
+            {
+                return turns * 360.0f;
+            }
+
+            return turns;
+        }
+
+        public static float Normalize(float turns, AngleRange range)
+        {
+            if (range == AngleRange.UNSIGNED_RANGE)
+            {
+                if (0.0f <= turns && turns < 1.0f)
+                {
+                    return turns;
+                }
+            }
+            else
+            {
+                if (-0.5f < turns && turns <= 0.5f)
+                {
+                    return turns;
+                }
+            }
+
+            float rest = turns - MathF.Floor(turns);
+
+            if (range == AngleRange.SIGNED_RANGE && rest > 0.5f)
+            {
+                rest -= 1.0f;
+            }
+
+            return rest;
+        }
+    }
+}

BasicGeometry/FP32Utility.cs

@@ -0,0 +1,21 @@
+using System;
+
+namespace BasicGeometry
+{
+    public class FP32Utility
+    {
+        public const float EPSYLON = 5E-7f;
+        public const float TWO_EPSYLON = 1E-6f;
+        public const float SQUARE_EPSYLON = 2.5E-13f;
+        
+        public const float EPSYLON_EFFECTIVENESS_LIMIT = 1.0f;
+
+        public const float ONE_THIRD = 0.333333333f;
+        public const float ONE_SIXTH = 0.166666667f;
+        public const float ONE_NINETH = 0.111111111f;
+
+        public const float GOLDEN_RATIO_HIGH = 1.618034f;
+        public const float GOLDEN_RATIO_LOW = 0.618034f;
+    }
+}
+

BasicGeometry/FP32Vector2.cs

@@ -0,0 +1,257 @@
+/*
+ * Copyright 2019-2025 Andrey Pokidov <andrey.pokidov@gmail.com>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+using System;
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 1 Feb 2019
+ */
+
+namespace BasicGeometry
+{
+    public struct FP32Vector2
+    {
+        public static readonly FP32Vector2 ZERO = new FP32Vector2(0.0f, 0.0f);
+
+        public float x1 = 0.0f;
+        public float x2 = 0.0f;
+
+        public FP32Vector2(float x1, float x2)
+        {
+            this.x1 = x1;
+            this.x2 = x2;
+        }
+
+        public FP32Vector2(in FP32Vector2 vector)
+        {
+            this.x1 = vector.x1;
+            this.x2 = vector.x2;
+        }
+
+        public FP32Vector2(in FP64Vector2 vector)
+        {
+            this.x1 = (float)vector.x1;
+            this.x2 = (float)vector.x2;
+        }
+
+        public readonly float GetSquareModule()
+        {
+            return this.x1 * this.x1 + this.x2 * this.x2;
+        }
+
+        public readonly float GetModule()
+        {
+            return MathF.Sqrt(this.GetSquareModule());
+        }
+
+        public int Normalize()
+        {
+            float squareModule = this.GetSquareModule();
+
+            if (1.0f - FP32Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0f + FP32Utility.TWO_EPSYLON)
+            {
+                return 1;
+            }
+
+            if (squareModule <= FP32Utility.SQUARE_EPSYLON)
+            {
+                this.Reset();
+                return 0;
+            }
+
+            float module = MathF.Sqrt(squareModule);
+
+            this.x1 /= module;
+            this.x2 /= module;
+
+            return 1;
+        }
+
+        public void Reverse()
+        {
+            this.x1 = -this.x1;
+            this.x2 = -this.x2;
+        }
+
+        public readonly bool IsZero()
+        {
+            return this.GetSquareModule() <= FP32Utility.SQUARE_EPSYLON;
+        }
+
+        public readonly bool IsUnit()
+        {
+            float squareModule = this.GetSquareModule();
+            return 1.0f - FP32Utility.TWO_EPSYLON <= squareModule && squareModule <= FP32Utility.EPSYLON;
+        }
+
+        public void Reset()
+        {
+            this.x1 = 0.0f;
+            this.x2 = 0.0f;
+        }
+
+        public void SetValues(float x1, float x2)
+        {
+            this.x1 = x1;
+            this.x2 = x2;
+        }
+
+        public void SetValues(in FP32Vector2 vector)
+        {
+            this.x1 = vector.x1;
+            this.x2 = vector.x2;
+        }
+
+        public void SetValues(in FP64Vector2 vector)
+        {
+            this.x1 = (float)vector.x1;
+            this.x2 = (float)vector.x2;
+        }
+
+        public void SetReverseOf(in FP32Vector2 vector)
+        {
+            this.x1 = -vector.x1;
+            this.x2 = -vector.x2;
+        }
+
+        public void SetReverseOf(in FP64Vector2 vector)
+        {
+            this.x1 = -(float)vector.x1;
+            this.x2 = -(float)vector.x2;
+        }
+
+        public void AppendScaled(FP32Vector2 summand, float scale)
+        {
+            this.x1 += summand.x1 * scale;
+            this.x2 += summand.x2 * scale;
+        }
+        
+        public readonly override string ToString()
+        {
+            return String.Format("SPVector2({0}, {1})", this.x1, this.x2);
+        }
+
+        public static void Add(in FP32Vector2 vector1, in FP32Vector2 vector2, out FP32Vector2 sum)
+        {
+            sum.x1 = vector1.x1 + vector2.x1;
+            sum.x2 = vector1.x2 + vector2.x2;
+        }
+
+        public static void Subtract(in FP32Vector2 minuend, in FP32Vector2 subtrahend, out FP32Vector2 difference)
+        {
+            difference.x1 = minuend.x1 - subtrahend.x1;
+            difference.x2 = minuend.x2 - subtrahend.x2;
+        }
+
+        public static void Multiply(in FP32Vector2 multiplicand, float multiplier, out FP32Vector2 product)
+        {
+            product.x1 = multiplicand.x1 * multiplier;
+            product.x2 = multiplicand.x2 * multiplier;
+        }
+
+        public static void Divide(in FP32Vector2 dividend, float divisor, out FP32Vector2 quotient)
+        {
+            Multiply(dividend, 1.0f / divisor, out quotient);
+        }
+
+        public static void GetMean2(in FP32Vector2 vector1, in FP32Vector2 vector2, out FP32Vector2 result)
+        {
+            result.x1 = (vector1.x1 + vector2.x1) * 0.5f;
+            result.x2 = (vector1.x2 + vector2.x2) * 0.5f;
+        }
+
+        public static void GetMean3(in FP32Vector2 vector1, in FP32Vector2 vector2, in FP32Vector2 vector3, out FP32Vector2 result)
+        {
+            result.x1 = (vector1.x1 + vector2.x1 + vector3.x1) * FP32Utility.ONE_THIRD;
+            result.x2 = (vector1.x2 + vector2.x2 + vector3.x2) * FP32Utility.ONE_THIRD;
+        }
+
+        public static float GetScalarProduct(in FP32Vector2 vector1, in FP32Vector2 vector2)
+        {
+            return vector1.x1 * vector2.x1 + vector1.x2 * vector2.x2;
+        }
+ 
+        public static float GetCrossProduct(in FP32Vector2 vector1, in FP32Vector2 vector2)
+        {
+            return vector1.x1 * vector2.x2 - vector1.x2 * vector2.x1;
+        }
+
+        public static float GetAngle(in FP32Vector2 vector1, in FP32Vector2 vector2, AngleUnit unit)
+        {
+            float squareModule1 = vector1.GetSquareModule();
+
+            if (squareModule1 <= FP32Utility.SQUARE_EPSYLON)
+            {
+                return 0.0f;
+            }
+
+            float squareModule2 = vector2.GetSquareModule();
+
+            if (squareModule2 <= FP32Utility.SQUARE_EPSYLON)
+            {
+                return 0.0f;
+            }
+
+            float cosine = FP32Vector2.GetScalarProduct(vector1, vector2) / MathF.Sqrt(squareModule1 * squareModule2);
+
+            if (1.0f - FP32Utility.EPSYLON <= cosine)
+            {
+                return 0.0f;
+            }
+
+            if (cosine <= -(1.0f - FP32Utility.EPSYLON))
+            {
+                return FP32Angle.GetHalfCircle(unit);
+            }
+    
+            return FP32Radians.ToUnits(MathF.Acos(cosine), unit);
+        }
+    
+        public static float GetSquareDistance(in FP32Vector2 vector1, in FP32Vector2 vector2)
+        {
+            float dx1 = vector1.x1 - vector2.x1;
+            float dx2 = vector1.x2 - vector2.x2;
+        
+            return dx1 * dx1 + dx2 * dx2;
+        }
+    
+        public static float GetDistance(in FP32Vector2 vector1, in FP32Vector2 vector2)
+        {
+            return MathF.Sqrt(GetSquareDistance(vector1, vector2));
+        }
+    
+        public static bool AreEqual(in FP32Vector2 vector1, in FP32Vector2 vector2)
+        {
+            float squareModule1 = vector1.GetSquareModule();
+            float squareModule2 = vector2.GetSquareModule();
+            float squareModule3 = GetSquareDistance(vector1, vector2);
+
+            // 2.0f means dimension amount
+            if (squareModule1 < FP32Utility.EPSYLON_EFFECTIVENESS_LIMIT || squareModule2 < FP32Utility.EPSYLON_EFFECTIVENESS_LIMIT)
+            {
+                return squareModule3 < (2.0f * FP32Utility.SQUARE_EPSYLON);
+            }
+
+            if (squareModule1 <= squareModule2)
+            {
+                return squareModule3 <= (2.0f * FP32Utility.SQUARE_EPSYLON) * squareModule2;
+            }
+
+            return squareModule3 <= (2.0f * FP32Utility.SQUARE_EPSYLON) * squareModule1;
+        }
+    }
+}

BasicGeometry/FP32Vector3.cs

@@ -0,0 +1,300 @@
+/*
+ * Copyright 2019-2025 Andrey Pokidov <andrey.pokidov@gmail.com>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+using System;
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 1 Feb 2019
+ */
+
+namespace BasicGeometry
+{
+    public struct FP32Vector3
+    {
+        public static readonly FP32Vector3 ZERO = new FP32Vector3(0.0f, 0.0f, 0.0f);
+
+        public float x1 = 0.0f;
+        public float x2 = 0.0f;
+        public float x3 = 0.0f;
+
+        public FP32Vector3(float x1, float x2, float x3)
+        {
+            this.x1 = x1;
+            this.x2 = x2;
+            this.x3 = x3;
+        }
+
+        public FP32Vector3(in FP32Vector3 vector)
+        {
+            this.x1 = vector.x1;
+            this.x2 = vector.x2;
+            this.x3 = vector.x3;
+        }
+
+        public FP32Vector3(in FP64Vector3 vector)
+        {
+            this.x1 = (float)vector.x1;
+            this.x2 = (float)vector.x2;
+            this.x3 = (float)vector.x3;
+        }
+
+        public readonly float GetSquareModule()
+        {
+            return this.x1 * this.x1 + this.x2 * this.x2 + this.x3 * this.x3;
+        }
+
+        public readonly float GetModule()
+        {
+            return MathF.Sqrt(this.GetSquareModule());
+        }
+
+        public int Normalize()
+        {
+            float squareModule = this.GetSquareModule();
+
+            if (1.0f - FP32Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0f + FP32Utility.TWO_EPSYLON)
+            {
+                return 1;
+            }
+
+            if (squareModule <= FP32Utility.SQUARE_EPSYLON)
+            {
+                this.Reset();
+                return 0;
+            }
+
+            float module = MathF.Sqrt(squareModule);
+
+            this.x1 /= module;
+            this.x2 /= module;
+            this.x3 /= module;
+
+            return 1;
+        }
+
+        public void Reverse()
+        {
+            this.x1 = -this.x1;
+            this.x2 = -this.x2;
+            this.x3 = -this.x3;
+        }
+
+        public readonly bool IsZero()
+        {
+            return this.GetSquareModule() <= FP32Utility.SQUARE_EPSYLON;
+        }
+
+        public readonly bool IsUnit()
+        {
+            float squareModule = this.GetSquareModule();
+            return 1.0f - FP32Utility.TWO_EPSYLON <= squareModule && squareModule <= FP32Utility.EPSYLON;
+        }
+
+        public void Reset()
+        {
+            this.x1 = 0.0f;
+            this.x2 = 0.0f;
+            this.x3 = 0.0f;
+        }
+
+        public void SetValues(float x1, float x2, float x3)
+        {
+            this.x1 = x1;
+            this.x2 = x2;
+            this.x3 = x3;
+        }
+
+        public void SetValues(in FP64Vector3 vector)
+        {
+            this.x1 = (float)vector.x1;
+            this.x2 = (float)vector.x2;
+            this.x3 = (float)vector.x3;
+        }
+
+        public void SetValues(in FP32Vector3 vector)
+        {
+            this.x1 = vector.x1;
+            this.x2 = vector.x2;
+            this.x3 = vector.x3;
+        }
+
+        public void SetReverseOf(in FP32Vector3 vector)
+        {
+            this.x1 = -vector.x1;
+            this.x2 = -vector.x2;
+            this.x3 = -vector.x3;
+        }
+
+        public void SetReverseOf(in FP64Vector3 vector)
+        {
+            this.x1 = -(float)vector.x1;
+            this.x2 = -(float)vector.x2;
+            this.x3 = -(float)vector.x3;
+        }
+
+        public void AppendScaled(FP32Vector3 summand, float scale)
+        {
+            this.x1 += summand.x1 * scale;
+            this.x2 += summand.x2 * scale;
+            this.x3 += summand.x3 * scale;
+        }
+
+        public readonly override string ToString()
+        {
+            return String.Format("SPVector3({0}, {1}, {2})", this.x1, this.x2, this.x3);
+        }
+
+        public static void Add(in FP32Vector3 vector1, in FP32Vector3 vector2, out FP32Vector3 sum)
+        {
+            sum.x1 = vector1.x1 + vector2.x1;
+            sum.x2 = vector1.x2 + vector2.x2;
+            sum.x3 = vector1.x3 + vector2.x3;
+        }
+
+        public static void Subtract(in FP32Vector3 minuend, in FP32Vector3 subtrahend, out FP32Vector3 difference)
+        {
+            difference.x1 = minuend.x1 - subtrahend.x1;
+            difference.x2 = minuend.x2 - subtrahend.x2;
+            difference.x3 = minuend.x3 - subtrahend.x3;
+        }
+
+        public static void Multiply(in FP32Vector3 multiplicand, float multiplier, out FP32Vector3 product)
+        {
+            product.x1 = multiplicand.x1 * multiplier;
+            product.x2 = multiplicand.x2 * multiplier;
+            product.x3 = multiplicand.x3 * multiplier;
+        }
+
+        public static void Divide(in FP32Vector3 dividend, float divisor, out FP32Vector3 quotient)
+        {
+            Multiply(dividend, 1.0f / divisor, out quotient);
+        }
+
+        public static void GetMean2(in FP32Vector3 vector1, in FP32Vector3 vector2, out FP32Vector3 result)
+        {
+            result.x1 = (vector1.x1 + vector2.x1) * 0.5f;
+            result.x2 = (vector1.x2 + vector2.x2) * 0.5f;
+            result.x3 = (vector1.x3 + vector2.x3) * 0.5f;
+        }
+
+        public static void GetMean3(in FP32Vector3 vector1, in FP32Vector3 vector2, in FP32Vector3 vector3, out FP32Vector3 result)
+        {
+            result.x1 = (vector1.x1 + vector2.x1 + vector3.x1)  * FP32Utility.ONE_THIRD;
+            result.x2 = (vector1.x2 + vector2.x2 + vector3.x2)  * FP32Utility.ONE_THIRD;
+            result.x3 = (vector1.x3 + vector2.x3 + vector3.x3)  * FP32Utility.ONE_THIRD;
+        }
+
+        public static float GetScalarProduct(in FP32Vector3 vector1, in FP32Vector3 vector2)
+        {
+            return vector1.x1 * vector2.x1 + vector1.x2 * vector2.x2 + vector1.x3 * vector2.x3;
+        }
+ 
+        public static void GetCrossProduct(in FP32Vector3 vector1, in FP32Vector3 vector2, out FP32Vector3 result)
+        {
+            float x1 = vector1.x2 * vector2.x3 - vector1.x3 * vector2.x2;
+            float x2 = vector1.x3 * vector2.x1 - vector1.x1 * vector2.x3;
+            float x3 = vector1.x1 * vector2.x2 - vector1.x2 * vector2.x1;
+
+            result.x1 = x1;
+            result.x2 = x2;
+            result.x3 = x3;
+        }
+
+        public static float GetTripleProduct(in FP32Vector3 vector1, in FP32Vector3 vector2, in FP32Vector3 vector3)
+        {
+            return vector1.x1 * (vector2.x2 * vector3.x3 - vector2.x3 * vector3.x2)
+                 + vector1.x2 * (vector2.x3 * vector3.x1 - vector2.x1 * vector3.x3)
+                 + vector1.x3 * (vector2.x1 * vector3.x2 - vector2.x2 * vector3.x1);
+        }
+ 
+        public static void GetDoubleCrossProduct(in FP32Vector3 vector1, in FP32Vector3 vector2, in FP32Vector3 vector3, out FP32Vector3 result)
+        {
+            // [a x [b x c]] = b * (a, c) - c * (a, b)
+            float ac = GetScalarProduct(vector1, vector3);
+            float ab = GetScalarProduct(vector1, vector2);
+
+            result.x1 = ac * vector2.x1 - ab * vector3.x1;
+            result.x2 = ac * vector2.x2 - ab * vector3.x2;
+            result.x3 = ac * vector2.x3 - ab * vector3.x3;
+        }
+
+        public static float GetAngle(in FP32Vector3 vector1, in FP32Vector3 vector2, AngleUnit unit)
+        {
+            float squareModule1 = vector1.GetSquareModule();
+
+            if (squareModule1 <= FP32Utility.SQUARE_EPSYLON)
+            {
+                return 0.0f;
+            }
+
+            float squareModule2 = vector2.GetSquareModule();
+
+            if (squareModule2 <= FP32Utility.SQUARE_EPSYLON)
+            {
+                return 0.0f;
+            }
+
+            float cosine = FP32Vector3.GetScalarProduct(vector1, vector2) / MathF.Sqrt(squareModule1 * squareModule2);
+
+            if (1.0f - FP32Utility.EPSYLON <= cosine)
+            {
+                return 0.0f;
+            }
+
+            if (cosine <= -(1.0f - FP32Utility.EPSYLON))
+            {
+                return FP32Angle.GetHalfCircle(unit);
+            }
+    
+            return FP32Radians.ToUnits(MathF.Acos(cosine), unit);
+        }
+    
+        public static float GetSquareDistance(in FP32Vector3 vector1, in FP32Vector3 vector2)
+        {
+            float dx1 = vector1.x1 - vector2.x1;
+            float dx2 = vector1.x2 - vector2.x2;
+            float dx3 = vector1.x3 - vector2.x3;
+        
+            return dx1 * dx1 + dx2 * dx2 + dx3 * dx3;
+        }
+    
+        public static float GetDistance(in FP32Vector3 vector1, in FP32Vector3 vector2)
+        {
+            return MathF.Sqrt(GetSquareDistance(vector1, vector2));
+        }
+    
+        public static bool AreEqual(in FP32Vector3 vector1, in FP32Vector3 vector2)
+        {
+            float squareModule1 = vector1.GetSquareModule();
+            float squareModule2 = vector2.GetSquareModule();
+            float squareModule3 = GetSquareDistance(vector1, vector2);
+
+            // 3.0f means dimension amount
+            if (squareModule1 < FP32Utility.EPSYLON_EFFECTIVENESS_LIMIT || squareModule2 < FP32Utility.EPSYLON_EFFECTIVENESS_LIMIT)
+            {
+                return squareModule3 < (3.0f * FP32Utility.SQUARE_EPSYLON);
+            }
+
+            if (squareModule1 <= squareModule2)
+            {
+                return squareModule3 <= (3.0f * FP32Utility.SQUARE_EPSYLON) * squareModule2;
+            }
+
+            return squareModule3 <= (3.0f * FP32Utility.SQUARE_EPSYLON) * squareModule1;
+        }
+    }
+}

BasicGeometry/FP32Versor.cs

@@ -0,0 +1,304 @@
+/*
+ * Copyright 2019-2025 Andrey Pokidov <andrey.pokidov@gmail.com>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/*
+ * Author: Andrey Pokidov
+* Date: 20 Oct 2024
+ */
+
+namespace BasicGeometry
+{
+    public struct FP32Versor
+    {
+        private float s0 = 1.0f;
+        private float x1 = 0.0f;
+        private float x2 = 0.0f;
+        private float x3 = 0.0f;
+
+        public FP32Versor(float s0, float x1, float x2, float x3)
+        {
+            LoadValues(s0, x1, x2, x3, out this);
+        }
+
+        public FP32Versor(in FP32Versor versor)
+        {
+            this.s0 = versor.s0;
+            this.x1 = versor.x1;
+            this.x2 = versor.x2;
+            this.x3 = versor.x3;
+        }
+
+        public FP32Versor(in FP64Versor versor)
+        {
+            this.s0 = (float)versor.GetScalar();
+            this.x1 = (float)versor.GetX1();
+            this.x2 = (float)versor.GetX2();
+            this.x3 = (float)versor.GetX3();
+        }
+
+        public readonly float GetScalar()
+        {
+            return this.s0;
+        }
+
+        public readonly float GetX1()
+        {
+            return this.x1;
+        }
+
+        public readonly float GetX2()
+        {
+            return this.x2;
+        }
+
+        public readonly float GetX3()
+        {
+            return this.x3;
+        }
+
+        public readonly bool IsIdle()
+        {
+            return this.s0 <= -(1.0f - FP32Utility.EPSYLON) || (1.0f - FP32Utility.EPSYLON) <= this.s0;
+        }
+
+        public void Reset()
+        {
+            this.s0 = 1.0f;
+            this.x1 = 0.0f;
+            this.x2 = 0.0f;
+            this.x3 = 0.0f;
+        }
+
+        public void Invert()
+        {
+            this.x1 = -this.x1;
+            this.x2 = -this.x2;
+            this.x3 = -this.x3;
+        }
+
+        public readonly float GetAngle(AngleUnit unit)
+        {
+            if (this.s0 <= -(1.0f - FP32Utility.TWO_EPSYLON) || 1.0f - FP32Utility.TWO_EPSYLON <= this.s0) {
+                return 0.0f;
+            }
+
+            if (-FP32Utility.EPSYLON <= this.s0 && this.s0 <= FP32Utility.EPSYLON)
+            {
+                return FP32Angle.GetHalfCircle(unit);
+            }
+
+            return FP32Radians.ToUnits(2.0f * MathF.Acos(s0), unit);
+        }
+
+        public readonly void MakeRotationMatrix(out FP32Matrix3x3 matrix)
+        {
+            float s0s0 = this.s0 * this.s0;
+            float x1x1 = this.x1 * this.x1;
+            float x2x2 = this.x1 * this.x2;
+            float x3x3 = this.x1 * this.x3;
+
+            float s0x1 = 2.0f * this.s0 * this.x1;
+            float s0x2 = 2.0f * this.s0 * this.x2;
+            float s0x3 = 2.0f * this.s0 * this.x3;
+
+            float x1x2 = 2.0f * this.x1 * this.x2;
+            float x1x3 = 2.0f * this.x1 * this.x3;
+            float x2x3 = 2.0f * this.x2 * this.x3;
+
+            matrix.r1c1 = (s0s0 + x1x1) - (x2x2 + x3x3);
+            matrix.r2c2 = (s0s0 + x2x2) - (x1x1 + x3x3);
+            matrix.r3c3 = (s0s0 + x3x3) - (x1x1 + x2x2);
+
+            matrix.r1c2 = x1x2 - s0x3;
+            matrix.r2c3 = x2x3 - s0x1;
+            matrix.r3c1 = x1x3 - s0x2;
+
+            matrix.r2c1 = x1x2 + s0x3;
+            matrix.r3c2 = x2x3 + s0x1;
+            matrix.r1c3 = x1x3 + s0x2;
+        }
+
+        public readonly void MakeReverseMatrix(out FP32Matrix3x3 matrix)
+        {
+            float s0s0 = this.s0 * this.s0;
+            float x1x1 = this.x1 * this.x1;
+            float x2x2 = this.x1 * this.x2;
+            float x3x3 = this.x1 * this.x3;
+
+            float s0x1 = 2.0f * this.s0 * this.x1;
+            float s0x2 = 2.0f * this.s0 * this.x2;
+            float s0x3 = 2.0f * this.s0 * this.x3;
+
+            float x1x2 = 2.0f * this.x1 * this.x2;
+            float x1x3 = 2.0f * this.x1 * this.x3;
+            float x2x3 = 2.0f * this.x2 * this.x3;
+
+            matrix.r1c1 = (s0s0 + x1x1) - (x2x2 + x3x3);
+            matrix.r2c2 = (s0s0 + x2x2) - (x1x1 + x3x3);
+            matrix.r3c3 = (s0s0 + x3x3) - (x1x1 + x2x2);
+
+            matrix.r1c2 = x1x2 + s0x3;
+            matrix.r2c3 = x2x3 + s0x1;
+            matrix.r3c1 = x1x3 + s0x2;
+
+            matrix.r2c1 = x1x2 - s0x3;
+            matrix.r3c2 = x2x3 - s0x1;
+            matrix.r1c3 = x1x3 - s0x2;
+        }
+
+        public void SetValues(float s0, float x1, float x2, float x3)
+        {
+            this.s0 = s0;
+            this.x1 = x1;
+            this.x2 = x2;
+            this.x3 = x3;
+
+            float squareModule = (s0 * s0 + x1 * x1) + (x2 * x2 + x3 * x3);
+            
+            if (1.0f - FP32Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0f + FP32Utility.TWO_EPSYLON)
+            {
+                return;
+            }
+
+            this.Normalize(squareModule);
+        }
+
+        public void SetValues(in FP32Versor versor)
+        {
+            this.s0 = versor.s0;
+            this.x1 = versor.x1;
+            this.x2 = versor.x2;
+            this.x3 = versor.x3;
+        }
+
+        public void SetValues(in FP64Versor versor)
+        {
+            this.s0 = (float) versor.GetScalar();
+            this.x1 = (float) versor.GetX1();
+            this.x2 = (float) versor.GetX2();
+            this.x3 = (float) versor.GetX3();
+        }
+
+        public void SetInverted(in FP32Versor versor)
+        {
+            this.s0 = versor.s0;
+            this.x1 = -versor.x1;
+            this.x2 = -versor.x2;
+            this.x3 = -versor.x3;
+        }
+
+        public void SetInverted(in FP64Versor versor)
+        {
+            this.s0 = (float) versor.GetScalar();
+            this.x1 = (float) -versor.GetX1();
+            this.x2 = (float) -versor.GetX2();
+            this.x3 = (float) -versor.GetX3();
+        }
+
+        public readonly void Turn(in FP32Vector3 vector, out FP32Vector3 result)
+        {
+            float tx1 = 2.0f * (this.x2 * vector.x3 - this.x3 * vector.x2);
+            float tx2 = 2.0f * (this.x3 * vector.x1 - this.x1 * vector.x3);
+            float tx3 = 2.0f * (this.x1 * vector.x2 - this.x2 * vector.x1);
+
+            float x1 = (vector.x1 + tx1 * this.s0) + (this.x2 * tx3 - this.x3 * tx2);
+            float x2 = (vector.x2 + tx2 * this.s0) + (this.x3 * tx1 - this.x1 * tx3);
+            float x3 = (vector.x3 + tx3 * this.s0) + (this.x1 * tx2 - this.x2 * tx1);
+
+            result.x1 = x1;
+            result.x2 = x2;
+            result.x3 = x3;
+        }
+
+        public readonly void TurnBack(in FP32Vector3 vector, out FP32Vector3 result)
+        {
+            float tx1 = 2.0f * (this.x2 * vector.x3 - this.x3 * vector.x2);
+            float tx2 = 2.0f * (this.x3 * vector.x1 - this.x1 * vector.x3);
+            float tx3 = 2.0f * (this.x1 * vector.x2 - this.x2 * vector.x1);
+
+            float x1 = (vector.x1 - tx1 * this.s0) + (this.x2 * tx3 - this.x3 * tx2);
+            float x2 = (vector.x2 - tx2 * this.s0) + (this.x3 * tx1 - this.x1 * tx3);
+            float x3 = (vector.x3 - tx3 * this.s0) + (this.x1 * tx2 - this.x2 * tx1);
+
+            result.x1 = x1;
+            result.x2 = x2;
+            result.x3 = x3;
+        }
+
+        private void Normalize(float squareModule)
+        {
+            if (squareModule <= FP32Utility.SQUARE_EPSYLON)
+            {
+                this.Reset();
+                return;
+            }
+
+            float multiplier = MathF.Sqrt(1.0f / squareModule);
+
+            this.s0 *= multiplier;
+            this.x1 *= multiplier;
+            this.x2 *= multiplier;
+            this.x3 *= multiplier;
+        }
+
+        public static void Combine(in FP32Versor second, in FP32Versor first, out FP32Versor result)
+        {
+            float s0 = (second.s0 * first.s0 - second.x1 * first.x1) - (second.x2 * first.x2 + second.x3 * first.x3);
+            float x1 = (second.x1 * first.s0 + second.s0 * first.x1) - (second.x3 * first.x2 - second.x2 * first.x3);
+            float x2 = (second.x2 * first.s0 + second.s0 * first.x2) - (second.x1 * first.x3 - second.x3 * first.x1);
+            float x3 = (second.x3 * first.s0 + second.s0 * first.x3) - (second.x2 * first.x1 - second.x1 * first.x2);
+
+            float squareModule = (s0 * s0 + x1 * x1) + (x2 * x2 + x3 * x3);
+
+            result.s0 = s0;
+            result.x1 = x1;
+            result.x2 = x2;
+            result.x3 = x3;
+
+            if (1.0f - FP32Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0f + FP32Utility.TWO_EPSYLON)
+            {
+                return;
+            }
+
+            result.Normalize(squareModule);
+        }
+
+        public static void LoadIdle(out FP32Versor versor)
+        {
+            versor.s0 = 1.0f;
+            versor.x1 = 0.0f;
+            versor.x2 = 0.0f;
+            versor.x3 = 0.0f;
+        }
+        
+        public static void LoadValues(float s0, float x1, float x2, float x3, out FP32Versor versor)
+        {
+            versor.s0 = s0;
+            versor.x1 = x1;
+            versor.x2 = x2;
+            versor.x3 = x3;
+
+            float squareModule = (s0 * s0 + x1 * x1) + (x2 * x2 + x3 * x3);
+            
+            if (1.0f - FP32Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0f + FP32Utility.TWO_EPSYLON)
+            {
+                return;
+            }
+
+            versor.Normalize(squareModule);
+        }
+    }
+}

BasicGeometry/FP64Angle.cs

@@ -0,0 +1,133 @@
+/*
+ * Copyright 2019-2025 Andrey Pokidov <andrey.pokidov@gmail.com>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+using System;
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 1 Feb 2019
+ */
+
+namespace BasicGeometry
+{
+    public static class FP64Angle
+    {
+        public static double ToRadians(double angle, AngleUnit unit)
+        {
+            if (unit == AngleUnit.DEGREES)
+            {
+                return angle * FP64Degrees.RADIANS_IN_DEGREE;
+            }
+
+            if (unit == AngleUnit.TURNS)
+            {
+                return angle * FP64Radians.TWO_PI;
+            }
+
+            return angle;
+        }
+
+        public static double ToDegrees(double angle, AngleUnit unit)
+        {
+            if (unit == AngleUnit.RADIANS)
+            {
+                return angle * FP64Radians.DEGREES_IN_RADIAN;
+            }
+
+            if (unit == AngleUnit.TURNS)
+            {
+                return angle * 360.0;
+            }
+
+            return angle;
+        }
+
+        public static double ToTurns(double angle, AngleUnit unit)
+        {
+            if (unit == AngleUnit.RADIANS)
+            {
+                return angle * FP64Radians.TURNS_IN_RADIAN;
+            }
+
+            if (unit == AngleUnit.DEGREES)
+            {
+                return angle * FP64Degrees.TURNS_IN_DEGREE;
+            }
+
+            return angle;
+        }
+
+        public static double GetFullCircle(AngleUnit unit)
+        {
+            if (unit == AngleUnit.DEGREES)
+            {
+                return 360.0;
+            }
+
+            if (unit == AngleUnit.TURNS)
+            {
+                return 1.0;
+            }
+
+            return FP64Radians.TWO_PI;
+        }
+    
+        public static double GetHalfCircle(AngleUnit unit)
+        {
+            if (unit == AngleUnit.DEGREES)
+            {
+                return 180.0;
+            }
+
+            if (unit == AngleUnit.TURNS)
+            {
+                return 0.5;
+            }
+
+            return FP64Radians.PI;
+        }
+    
+        public static double GetQuarterCircle(AngleUnit unit)
+        {
+            if (unit == AngleUnit.DEGREES)
+            {
+                return 90.0;
+            }
+
+            if (unit == AngleUnit.TURNS)
+            {
+                return 0.25;
+            }
+
+            return FP64Radians.HALF_OF_PI;
+        }
+
+        public static double Normalize(double angle, AngleUnit unit, AngleRange range)
+        {
+            if (unit == AngleUnit.DEGREES)
+            {
+                return FP64Degrees.Normalize(angle, range);
+            }
+
+            if (unit == AngleUnit.TURNS)
+            {
+                return FP64Turns.Normalize(angle, range);
+            }
+
+            return FP64Radians.Normalize(angle, range);
+        }
+    }
+}

BasicGeometry/FP64Degrees.cs

@@ -0,0 +1,68 @@
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 18 Nov 2024
+ */
+
+namespace BasicGeometry
+{
+    public class FP64Degrees
+    {
+        public const double RADIANS_IN_DEGREE = 1.74532925199432958E-2;
+        public const double TURNS_IN_DEGREE = 2.77777777777777778E-3;
+
+        public static double ToRadians(double degrees)
+        {
+            return degrees * RADIANS_IN_DEGREE;
+        }
+
+        public static double ToTurns(double degrees)
+        {
+            return degrees * TURNS_IN_DEGREE;
+        }
+
+        public static double ToUnits(double degrees, AngleUnit toUnit)
+        {
+            if (toUnit == AngleUnit.RADIANS)
+            {
+                return degrees * RADIANS_IN_DEGREE;
+            }
+
+            if (toUnit == AngleUnit.TURNS)
+            {
+                return degrees * TURNS_IN_DEGREE;
+            }
+
+            return degrees;
+        }
+
+        public static double Normalize(double degrees, AngleRange range)
+        {
+            if (range == AngleRange.UNSIGNED_RANGE)
+            {
+                if (0.0 <= degrees && degrees < 360.0)
+                {
+                    return degrees;
+                }
+            }
+            else
+            {
+                if (-180.0 < degrees && degrees <= 180.0)
+                {
+                    return degrees;
+                }
+            }
+
+            double turns = degrees * TURNS_IN_DEGREE;
+
+            turns -= Math.Floor(turns);
+
+            if (range == AngleRange.SIGNED_RANGE && turns > 0.5)
+            {
+                turns -= 1.0;
+            }
+
+            return turns * 360.0;
+        }
+    }
+}

BasicGeometry/FP64Matrix2x2.cs

@@ -0,0 +1,294 @@
+/*
+ * Copyright 2019-2025 Andrey Pokidov <andrey.pokidov@gmail.com>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+using System;
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 10 Feb 2019
+ */
+
+namespace BasicGeometry
+{
+    public struct FP64Matrix2x2
+    {
+        public double r1c1 = 0.0, r1c2 = 0.0;
+        public double r2c1 = 0.0, r2c2 = 0.0;
+
+        public FP64Matrix2x2(double d1, double d2)
+        {
+            this.r1c1 = d1;
+            this.r2c2 = d2;
+        }
+
+        public FP64Matrix2x2(in FP64Matrix2x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+        }
+
+        public FP64Matrix2x2(in FP32Matrix2x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+        }
+
+        public readonly double GetDeterminant()
+        {
+            return this.r1c1 * this.r2c2 - this.r1c2 * this.r2c1;
+        }
+
+        public readonly bool IsSingular()
+        {
+            double determinant = this.GetDeterminant();
+            return -FP64Utility.EPSYLON <= determinant && determinant <= FP64Utility.EPSYLON;
+        }
+
+        public void Transpose()
+        {
+            (this.r1c2, this.r2c1) = (this.r2c1, this.r1c2);
+        }
+
+        public bool Invert()
+        {
+            double determinant = this.GetDeterminant();
+
+            if (-FP64Utility.EPSYLON <= determinant && determinant <= FP64Utility.EPSYLON)
+            {
+                return false;
+            }
+
+            double r1c1 =  this.r2c2;
+            double r1c2 = -this.r1c2;
+
+            double r2c1 = -this.r2c1;
+            double r2c2 =  this.r1c1;
+
+            double multiplier = 1.0 / determinant;
+
+            this.r1c1 = r1c1 * multiplier;
+            this.r1c2 = r1c2 * multiplier;
+
+            this.r2c1 = r2c1 * multiplier;
+            this.r2c2 = r2c2 * multiplier;
+
+            return true;
+        }
+
+        public void Reset()
+        {
+            this.r1c1 = 0.0;
+            this.r1c2 = 0.0;
+
+            this.r2c1 = 0.0;
+            this.r2c2 = 0.0;
+        }
+
+        public void SetToIdentity()
+        {
+            this.r1c1 = 1.0;
+            this.r1c2 = 0.0;
+
+            this.r2c1 = 0.0;
+            this.r2c2 = 1.0;
+        }
+
+        public void SetToDiagonal(double d1, double d2)
+        {
+            this.r1c1 = d1;
+            this.r1c2 = 0.0;
+
+            this.r2c1 = 0.0;
+            this.r2c2 = d2;
+        }
+
+        public void SetValues(in FP64Matrix2x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+        }
+
+        public void SetValues(in FP32Matrix2x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+        }
+
+        public void SetTransposedOf(in FP64Matrix2x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r2c2 = matrix.r2c2;
+            (this.r1c2, this.r2c1) = (matrix.r2c1, matrix.r1c2);
+        }
+
+        public void SetTransposedOf(in FP32Matrix2x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r2c1;
+            
+            this.r2c1 = matrix.r1c2;
+            this.r2c2 = matrix.r2c2;
+        }
+
+        public bool SetInvertedOf(in FP64Matrix2x2 matrix)
+        {
+            double determinant = matrix.GetDeterminant();
+
+            if (-FP64Utility.EPSYLON <= determinant && determinant <= FP64Utility.EPSYLON)
+            {
+                return false;
+            }
+
+            double r1c1 =  matrix.r2c2;
+            double r1c2 = -matrix.r1c2;
+
+            double r2c1 = -matrix.r2c1;
+            double r2c2 =  matrix.r1c1;
+
+            double multiplier = 1.0 / determinant;
+
+            this.r1c1 = r1c1 * multiplier;
+            this.r1c2 = r1c2 * multiplier;
+
+            this.r2c1 = r2c1 * multiplier;
+            this.r2c2 = r2c2 * multiplier;
+
+            return true;
+        }
+
+        public void SetRow1(double c1, double c2)
+        {
+            this.r1c1 = c1;
+            this.r1c2 = c2;
+        }
+
+        public void SetRow2(double c1, double c2)
+        {
+            this.r2c1 = c1;
+            this.r2c2 = c2;
+        }
+
+        public void SetColumn1(double r1, double r2)
+        {
+            this.r1c1 = r1;
+            this.r2c1 = r2;
+        }
+
+        public void SetColumn2(double r1, double r2)
+        {
+            this.r1c2 = r1;
+            this.r2c2 = r2;
+        }
+
+        public void AppendScaled(in FP64Matrix2x2 matrix, double scale)
+        {
+            this.r1c1 += matrix.r1c1 * scale;
+            this.r1c2 += matrix.r1c2 * scale;
+
+            this.r2c1 += matrix.r2c1 * scale;
+            this.r2c2 += matrix.r2c2 * scale;
+        }
+
+        public static void Add(in FP64Matrix2x2 matrix1, in FP64Matrix2x2 matrix2, out FP64Matrix2x2 sum)
+        {
+            sum.r1c1 = matrix1.r1c1 + matrix2.r1c1;
+            sum.r1c2 = matrix1.r1c2 + matrix2.r1c2;
+
+            sum.r2c1 = matrix1.r2c1 + matrix2.r2c1;
+            sum.r2c2 = matrix1.r2c2 + matrix2.r2c2;
+        }
+
+        public static void Subtract(in FP64Matrix2x2 minuend, in FP64Matrix2x2 subtrahend, out FP64Matrix2x2 difference)
+        {
+            difference.r1c1 = minuend.r1c1 - subtrahend.r1c1;
+            difference.r1c2 = minuend.r1c2 - subtrahend.r1c2;
+
+            difference.r2c1 = minuend.r2c1 - subtrahend.r2c1;
+            difference.r2c2 = minuend.r2c2 - subtrahend.r2c2;
+        }
+
+        public static void Multiply(in FP64Matrix2x2 multiplicand, double multiplier, out FP64Matrix2x2 product)
+        {
+            product.r1c1 = multiplicand.r1c1 * multiplier;
+            product.r1c2 = multiplicand.r1c2 * multiplier;
+
+            product.r2c1 = multiplicand.r2c1 * multiplier;
+            product.r2c2 = multiplicand.r2c2 * multiplier;
+        }
+
+        public static void Divide(in FP64Matrix2x2 dividend, double divisor, out FP64Matrix2x2 quotient)
+        {
+            Multiply(dividend, 1.0 / divisor, out quotient);
+        }
+
+        public static void GetRightProduct(in FP64Matrix2x2 matrix, in FP64Vector2 vector, out FP64Vector2 result)
+        {
+            double x1 = matrix.r1c1 * vector.x1 + matrix.r1c2 * vector.x2;
+            double x2 = matrix.r2c1 * vector.x1 + matrix.r2c2 * vector.x2;
+
+            result.x1 = x1;
+            result.x2 = x2;
+        }
+
+        public static void GetLeftProduct(in FP64Vector2 vector, in FP64Matrix2x2 matrix, out FP64Vector2 result)
+        {
+            double x1 = vector.x1 * matrix.r1c1 + vector.x2 * matrix.r2c1;
+            double x2 = vector.x1 * matrix.r1c2 + vector.x2 * matrix.r2c2;
+
+            result.x1 = x1;
+            result.x2 = x2;
+        }
+
+        public static void LoadZero(out FP64Matrix2x2 matrix)
+        {
+            matrix.r1c1 = 0.0;
+            matrix.r1c2 = 0.0;
+
+            matrix.r2c1 = 0.0;
+            matrix.r2c2 = 0.0;
+        }
+
+        public static void LoadIdentity(out FP64Matrix2x2 matrix)
+        {
+            matrix.r1c1 = 1.0;
+            matrix.r1c2 = 0.0;
+
+            matrix.r2c1 = 0.0;
+            matrix.r2c2 = 1.0;
+        }
+
+        public static void LoadDiagonal(double d1, double d2, out FP64Matrix2x2 matrix)
+        {
+            matrix.r1c1 = d1;
+            matrix.r1c2 = 0.0;
+
+            matrix.r2c1 = 0.0;
+            matrix.r2c2 = d2;
+        }
+    }
+}

BasicGeometry/FP64Matrix2x3.cs

@@ -0,0 +1,237 @@
+/*
+ * Copyright 2019-2025 Andrey Pokidov <andrey.pokidov@gmail.com>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+using System;
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 11 Nov 2024
+ */
+namespace BasicGeometry
+{
+    public struct FP64Matrix2x3
+    {
+        public double r1c1 = 0.0, r1c2 = 0.0;
+        public double r2c1 = 0.0, r2c2 = 0.0;
+        public double r3c1 = 0.0, r3c2 = 0.0;
+    
+        public FP64Matrix2x3(in FP64Matrix2x3 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+        
+            this.r3c1 = matrix.r3c1;
+            this.r3c2 = matrix.r3c2;
+        }
+    
+        public FP64Matrix2x3(in FP32Matrix2x3 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+        
+            this.r3c1 = matrix.r3c1;
+            this.r3c2 = matrix.r3c2;
+        }
+    
+        public FP64Matrix2x3(in FP64Matrix3x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r2c1;
+
+            this.r2c1 = matrix.r1c2;
+            this.r2c2 = matrix.r2c2;
+        
+            this.r3c1 = matrix.r1c3;
+            this.r3c2 = matrix.r2c3;
+        }
+    
+        public FP64Matrix2x3(in FP32Matrix3x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r2c1;
+
+            this.r2c1 = matrix.r1c2;
+            this.r2c2 = matrix.r2c2;
+        
+            this.r3c1 = matrix.r1c3;
+            this.r3c2 = matrix.r2c3;
+        }
+    
+        public void Reset()
+        {
+            this.r1c1 = 0.0;
+            this.r1c2 = 0.0;
+        
+            this.r2c1 = 0.0;
+            this.r2c2 = 0.0;
+        
+            this.r3c1 = 0.0;
+            this.r3c2 = 0.0;
+        }
+    
+        public void SetValues(in FP64Matrix2x3 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+        
+            this.r3c1 = matrix.r3c1;
+            this.r3c2 = matrix.r3c2;
+        }
+    
+        public void SetValues(in FP32Matrix2x3 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+        
+            this.r3c1 = matrix.r3c1;
+            this.r3c2 = matrix.r3c2;
+        }
+    
+        public void SetTransposed(in FP64Matrix3x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r2c1;
+
+            this.r2c1 = matrix.r1c2;
+            this.r2c2 = matrix.r2c2;
+        
+            this.r3c1 = matrix.r1c3;
+            this.r3c2 = matrix.r2c3;
+        }
+    
+        public void SetTransposed(in FP32Matrix3x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r2c1;
+
+            this.r2c1 = matrix.r1c2;
+            this.r2c2 = matrix.r2c2;
+        
+            this.r3c1 = matrix.r1c3;
+            this.r3c2 = matrix.r2c3;
+        }
+
+        public void SetRow1(double c1, double c2)
+        {
+            this.r1c1 = c1;
+            this.r1c2 = c2;
+        }
+
+        public void SetRow2(double c1, double c2)
+        {
+            this.r2c1 = c1;
+            this.r2c2 = c2;
+        }
+
+        public void SetRow3(double c1, double c2)
+        {
+            this.r3c1 = c1;
+            this.r3c2 = c2;
+        }
+
+        public void SetColumn1(double r1, double r2, double r3)
+        {
+            this.r1c1 = r1;
+            this.r2c1 = r2;
+            this.r3c1 = r3;
+        }
+
+        public void SetColumn2(double r1, double r2, double r3)
+        {
+            this.r1c2 = r1;
+            this.r2c2 = r2;
+            this.r3c2 = r3;
+        }
+
+        public void AppendScaled(in FP64Matrix2x3 matrix, double scale)
+        {
+            this.r1c1 += matrix.r1c1 * scale;
+            this.r1c2 += matrix.r1c2 * scale;
+
+            this.r2c1 += matrix.r2c1 * scale;
+            this.r2c2 += matrix.r2c2 * scale;
+
+            this.r3c1 += matrix.r3c1 * scale;
+            this.r3c2 += matrix.r3c2 * scale;
+        }
+
+        public static void Add(in FP64Matrix2x3 matrix1, in FP64Matrix2x3 matrix2, out FP64Matrix2x3 sum)
+        {
+            sum.r1c1 = matrix1.r1c1 + matrix2.r1c1;
+            sum.r1c2 = matrix1.r1c2 + matrix2.r1c2;
+
+            sum.r2c1 = matrix1.r2c1 + matrix2.r2c1;
+            sum.r2c2 = matrix1.r2c2 + matrix2.r2c2;
+
+            sum.r3c1 = matrix1.r3c1 + matrix2.r3c1;
+            sum.r3c2 = matrix1.r3c2 + matrix2.r3c2;
+        }
+
+        public static void Subtract(in FP64Matrix2x3 minuend, in FP64Matrix2x3 subtrahend, out FP64Matrix2x3 difference)
+        {
+            difference.r1c1 = minuend.r1c1 - subtrahend.r1c1;
+            difference.r1c2 = minuend.r1c2 - subtrahend.r1c2;
+
+            difference.r2c1 = minuend.r2c1 - subtrahend.r2c1;
+            difference.r2c2 = minuend.r2c2 - subtrahend.r2c2;
+
+            difference.r3c1 = minuend.r3c1 - subtrahend.r3c1;
+            difference.r3c2 = minuend.r3c2 - subtrahend.r3c2;
+        }
+
+        public static void Multiply(in FP64Matrix2x3 multiplicand, double multiplier, out FP64Matrix2x3 product)
+        {
+            product.r1c1 = multiplicand.r1c1 * multiplier;
+            product.r1c2 = multiplicand.r1c2 * multiplier;
+
+            product.r2c1 = multiplicand.r2c1 * multiplier;
+            product.r2c2 = multiplicand.r2c2 * multiplier;
+
+            product.r3c1 = multiplicand.r3c1 * multiplier;
+            product.r3c2 = multiplicand.r3c2 * multiplier;
+        }
+
+        public static void Divide(in FP64Matrix2x3 dividend, double divisor, out FP64Matrix2x3 quotient)
+        {
+            Multiply(dividend, 1.0 / divisor, out quotient);
+        }
+
+        public static void GetRightProduct(in FP64Matrix2x3 matrix, in FP64Vector2 vector, out FP64Vector3 result)
+        {
+            result.x1 = matrix.r1c1 * vector.x1 + matrix.r1c2 * vector.x2;
+            result.x2 = matrix.r2c1 * vector.x1 + matrix.r2c2 * vector.x2;
+            result.x3 = matrix.r3c1 * vector.x1 + matrix.r3c2 * vector.x2;
+        }
+
+        public static void GetLeftProduct(in FP64Vector3 vector, in FP64Matrix2x3 matrix, out FP64Vector2 result)
+        {
+            result.x1 = vector.x1 * matrix.r1c1 + vector.x2 * matrix.r2c1 + vector.x3 * matrix.r3c1;
+            result.x2 = vector.x1 * matrix.r1c2 + vector.x2 * matrix.r2c2 + vector.x3 * matrix.r3c2;
+        }
+    }
+}

BasicGeometry/FP64Matrix3x2.cs

@@ -0,0 +1,224 @@
+/*
+ * Copyright 2019-2025 Andrey Pokidov <andrey.pokidov@gmail.com>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+using System;
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 11 Nov 2024
+ */
+namespace BasicGeometry
+{
+    public struct FP64Matrix3x2
+    {
+        public double r1c1 = 0.0, r1c2 = 0.0, r1c3 = 0.0;
+        public double r2c1 = 0.0, r2c2 = 0.0, r2c3 = 0.0;
+
+        public FP64Matrix3x2(in FP64Matrix3x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+            this.r1c3 = matrix.r1c3;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+            this.r2c3 = matrix.r2c3;
+        }
+
+        public FP64Matrix3x2(in FP32Matrix3x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+            this.r1c3 = matrix.r1c3;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+            this.r2c3 = matrix.r2c3;
+        }
+
+        public FP64Matrix3x2(in FP64Matrix2x3 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r2c1;
+            this.r1c3 = matrix.r3c1;
+
+            this.r2c1 = matrix.r1c2;
+            this.r2c2 = matrix.r2c2;
+            this.r2c3 = matrix.r2c2;
+        }
+
+        public FP64Matrix3x2(in FP32Matrix2x3 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r2c1;
+            this.r1c3 = matrix.r3c1;
+
+            this.r2c1 = matrix.r1c2;
+            this.r2c2 = matrix.r2c2;
+            this.r2c3 = matrix.r2c2;
+        }
+
+        public void Reset()
+        {
+            this.r1c1 = 0.0;
+            this.r1c2 = 0.0;
+            this.r1c3 = 0.0;
+
+            this.r2c1 = 0.0;
+            this.r2c2 = 0.0;
+            this.r2c3 = 0.0;
+        }
+    
+
+        public void SetValues(in FP64Matrix3x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+            this.r1c3 = matrix.r1c3;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+            this.r2c3 = matrix.r2c3;
+        }
+
+        public void SetValues(in FP32Matrix3x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+            this.r1c3 = matrix.r1c3;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+            this.r2c3 = matrix.r2c3;
+        }
+
+        public void SetTransposed(in FP64Matrix2x3 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r2c1;
+            this.r1c3 = matrix.r3c1;
+
+            this.r2c1 = matrix.r1c2;
+            this.r2c2 = matrix.r2c2;
+            this.r2c3 = matrix.r2c2;
+        }
+
+        public void SetTransposed(in FP32Matrix2x3 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r2c1;
+            this.r1c3 = matrix.r3c1;
+
+            this.r2c1 = matrix.r1c2;
+            this.r2c2 = matrix.r2c2;
+            this.r2c3 = matrix.r2c2;
+        }
+    
+        public void SetRow1(double c1, double c2, double c3)
+        {
+            this.r1c1 = c1;
+            this.r1c2 = c2;
+            this.r1c3 = c3;
+        }
+
+        public void SetRow2(double c1, double c2, double c3)
+        {
+            this.r2c1 = c1;
+            this.r2c2 = c2;
+            this.r2c3 = c3;
+        }
+
+        public void SetColumn1(double r1, double r2)
+        {
+            this.r1c1 = r1;
+            this.r2c1 = r2;
+        }
+
+        public void SetColumn2(double r1, double r2)
+        {
+            this.r1c2 = r1;
+            this.r2c2 = r2;
+        }
+
+        public void SetColumn3(double r1, double r2)
+        {
+            this.r1c3 = r1;
+            this.r2c3 = r2;
+        }
+
+        public void AppendScaled(in FP64Matrix3x2 matrix, double scale)
+        {
+            this.r1c1 += matrix.r1c1 * scale;
+            this.r1c2 += matrix.r1c2 * scale;
+            this.r1c3 += matrix.r1c3 * scale;
+
+            this.r2c1 += matrix.r2c1 * scale;
+            this.r2c2 += matrix.r2c2 * scale;
+            this.r2c3 += matrix.r2c3 * scale;
+        }
+
+        public static void Add(in FP64Matrix3x2 matrix1, in FP64Matrix3x2 matrix2, out FP64Matrix3x2 sum)
+        {
+            sum.r1c1 = matrix1.r1c1 + matrix2.r1c1;
+            sum.r1c2 = matrix1.r1c2 + matrix2.r1c2;
+            sum.r1c3 = matrix1.r1c3 + matrix2.r1c3;
+
+            sum.r2c1 = matrix1.r2c1 + matrix2.r2c1;
+            sum.r2c2 = matrix1.r2c2 + matrix2.r2c2;
+            sum.r2c3 = matrix1.r2c3 + matrix2.r2c3;
+        }
+
+        public static void Subtract(in FP64Matrix3x2 minuend, in FP64Matrix3x2 subtrahend, out FP64Matrix3x2 difference)
+        {
+            difference.r1c1 = minuend.r1c1 - subtrahend.r1c1;
+            difference.r1c2 = minuend.r1c2 - subtrahend.r1c2;
+            difference.r1c3 = minuend.r1c3 - subtrahend.r1c3;
+
+            difference.r2c1 = minuend.r2c1 - subtrahend.r2c1;
+            difference.r2c2 = minuend.r2c2 - subtrahend.r2c2;
+            difference.r2c3 = minuend.r2c3 - subtrahend.r2c3;
+        }
+
+        public static void Multiply(in FP64Matrix3x2 multiplicand, double multiplier, out FP64Matrix3x2 product)
+        {
+            product.r1c1 = multiplicand.r1c1 * multiplier;
+            product.r1c2 = multiplicand.r1c2 * multiplier;
+            product.r1c3 = multiplicand.r1c3 * multiplier;
+
+            product.r2c1 = multiplicand.r2c1 * multiplier;
+            product.r2c2 = multiplicand.r2c2 * multiplier;
+            product.r2c3 = multiplicand.r2c3 * multiplier;
+        }
+
+        public static void Divide(in FP64Matrix3x2 dividend, double divisor, out FP64Matrix3x2 quotient)
+        {
+            Multiply(dividend, 1.0 / divisor, out quotient);
+        }
+
+        public static void GetRightProduct(in FP64Matrix3x2 matrix, in FP64Vector3 vector, out FP64Vector2 result)
+        {
+            result.x1 = matrix.r1c1 * vector.x1 + matrix.r1c2 * vector.x2 + matrix.r1c3 * vector.x3;
+            result.x2 = matrix.r2c1 * vector.x1 + matrix.r2c2 * vector.x2 + matrix.r2c3 * vector.x3;
+        }
+
+        public static void GetLeftProduct(in FP64Vector2 vector, in FP64Matrix3x2 matrix, out FP64Vector3 result)
+        {
+            result.x1 = vector.x1 * matrix.r1c1 + vector.x2 * matrix.r2c1;
+            result.x2 = vector.x1 * matrix.r1c2 + vector.x2 * matrix.r2c2;
+            result.x3 = vector.x1 * matrix.r1c3 + vector.x2 * matrix.r2c3;
+        }
+    }
+}

BasicGeometry/FP64Matrix3x3.cs

@@ -0,0 +1,434 @@
+/*
+ * Copyright 2019-2025 Andrey Pokidov <andrey.pokidov@gmail.com>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+using System;
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 10 Feb 2019
+ */
+
+namespace BasicGeometry
+{
+    public struct FP64Matrix3x3
+    {
+        public double r1c1 = 0.0, r1c2 = 0.0, r1c3 = 0.0;
+        public double r2c1 = 0.0, r2c2 = 0.0, r2c3 = 0.0;
+        public double r3c1 = 0.0, r3c2 = 0.0, r3c3 = 0.0;
+
+        public FP64Matrix3x3(double d1, double d2, double d3)
+        {
+            this.r1c1 = d1;
+            this.r2c2 = d2;
+            this.r3c3 = d3;
+        }
+
+        public FP64Matrix3x3(in FP64Matrix3x3 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+            this.r1c3 = matrix.r1c3;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+            this.r2c3 = matrix.r2c3;
+
+            this.r3c1 = matrix.r3c1;
+            this.r3c2 = matrix.r3c2;
+            this.r3c3 = matrix.r3c3;
+        }
+
+        public FP64Matrix3x3(in FP32Matrix3x3 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+            this.r1c3 = matrix.r1c3;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+            this.r2c3 = matrix.r2c3;
+
+            this.r3c1 = matrix.r3c1;
+            this.r3c2 = matrix.r3c2;
+            this.r3c3 = matrix.r3c3;
+        }
+
+        public readonly double GetDeterminant()
+        {
+            return this.r1c1 * (this.r2c2 * this.r3c3 - this.r2c3 * this.r3c2)
+                 + this.r1c2 * (this.r2c3 * this.r3c1 - this.r2c1 * this.r3c3)
+                 + this.r1c3 * (this.r2c1 * this.r3c2 - this.r2c2 * this.r3c1);
+        }
+
+        public readonly bool IsSingular()
+        {
+            double determinant = this.GetDeterminant();
+            return -FP64Utility.EPSYLON <= determinant && determinant <= FP64Utility.EPSYLON;
+        }
+
+        public void Transpose()
+        {
+            (this.r1c2, this.r2c1) = (this.r2c1, this.r1c2);
+            (this.r1c3, this.r3c1) = (this.r3c1, this.r1c3);
+            (this.r2c3, this.r3c2) = (this.r3c2, this.r2c3);
+        }
+
+        public bool Invert()
+        {
+            double determinant = this.GetDeterminant();
+    
+            if (-FP64Utility.EPSYLON <= determinant && determinant <= FP64Utility.EPSYLON) {
+                return false;
+            }
+
+            double r1c1 = this.r2c2 * this.r3c3 - this.r2c3 * this.r3c2;
+            double r1c2 = this.r1c3 * this.r3c2 - this.r1c2 * this.r3c3;
+            double r1c3 = this.r1c2 * this.r2c3 - this.r1c3 * this.r2c2;
+
+            double r2c1 = this.r2c3 * this.r3c1 - this.r2c1 * this.r3c3;
+            double r2c2 = this.r1c1 * this.r3c3 - this.r1c3 * this.r3c1;
+            double r2c3 = this.r1c3 * this.r2c1 - this.r1c1 * this.r2c3;
+
+            double r3c1 = this.r2c1 * this.r3c2 - this.r2c2 * this.r3c1;
+            double r3c2 = this.r1c2 * this.r3c1 - this.r1c1 * this.r3c2;
+            double r3c3 = this.r1c1 * this.r2c2 - this.r1c2 * this.r2c1;
+
+            double mutiplier = 1.0 / determinant;
+
+            this.r1c1 = r1c1 * mutiplier;
+            this.r1c2 = r1c2 * mutiplier;
+            this.r1c3 = r1c3 * mutiplier;
+
+            this.r2c1 = r2c1 * mutiplier;
+            this.r2c2 = r2c2 * mutiplier;
+            this.r2c3 = r2c3 * mutiplier;
+
+            this.r3c1 = r3c1 * mutiplier;
+            this.r3c2 = r3c2 * mutiplier;
+            this.r3c3 = r3c3 * mutiplier;
+
+            return true;
+        }
+
+        public void Reset()
+        {
+            this.r1c1 = 0.0;
+            this.r1c2 = 0.0;
+            this.r1c3 = 0.0;
+
+            this.r2c1 = 0.0;
+            this.r2c2 = 0.0;
+            this.r2c3 = 0.0;
+
+            this.r3c1 = 0.0;
+            this.r3c2 = 0.0;
+            this.r3c3 = 0.0;
+        }
+
+        public void SetToIdentity()
+        {
+            this.r1c1 = 1.0;
+            this.r1c2 = 0.0;
+            this.r1c3 = 0.0;
+
+            this.r2c1 = 0.0;
+            this.r2c2 = 1.0;
+            this.r2c3 = 0.0;
+
+            this.r3c1 = 0.0;
+            this.r3c2 = 0.0;
+            this.r3c3 = 1.0;
+        }
+
+        public void SetToDiagonal(double d1, double d2, double d3)
+        {
+            this.r1c1 = d1;
+            this.r1c2 = 0.0;
+            this.r1c3 = 0.0;
+
+            this.r2c1 = 0.0;
+            this.r2c2 = d2;
+            this.r2c3 = 0.0;
+
+            this.r2c1 = 0.0;
+            this.r2c2 = 0.0;
+            this.r2c3 = d3;
+        }
+
+        public void SetValues(FP64Matrix3x3 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+            this.r1c3 = matrix.r1c3;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+            this.r2c3 = matrix.r2c3;
+
+            this.r3c1 = matrix.r3c1;
+            this.r3c2 = matrix.r3c2;
+            this.r3c3 = matrix.r3c3;
+        }
+
+        public void SetValues(FP32Matrix3x3 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+            this.r1c3 = matrix.r1c3;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+            this.r2c3 = matrix.r2c3;
+
+            this.r3c1 = matrix.r3c1;
+            this.r3c2 = matrix.r3c2;
+            this.r3c3 = matrix.r3c3;
+        }
+
+        public void SetTransposedOf(in FP64Matrix3x3 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r2c2 = matrix.r2c2;
+            this.r3c3 = matrix.r3c3;
+
+            (this.r1c2, this.r2c1) = (matrix.r2c1, matrix.r1c2);
+            (this.r1c3, this.r3c1) = (matrix.r3c1, matrix.r1c3);
+            (this.r2c3, this.r3c2) = (matrix.r3c2, matrix.r2c3);
+        }
+
+        public void SetTransposedOf(in FP32Matrix3x3 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r2c2 = matrix.r2c2;
+            this.r3c3 = matrix.r3c3;
+
+            (this.r1c2, this.r2c1) = (matrix.r2c1, matrix.r1c2);
+            (this.r1c3, this.r3c1) = (matrix.r3c1, matrix.r1c3);
+            (this.r2c3, this.r3c2) = (matrix.r3c2, matrix.r2c3);
+        }
+
+        public bool SetInvertedOf(in FP64Matrix3x3 matrix)
+        {
+            double determinant = matrix.GetDeterminant();
+    
+            if (-FP64Utility.EPSYLON <= determinant && determinant <= FP64Utility.EPSYLON) {
+                return false;
+            }
+
+            double r1c1 = matrix.r2c2 * matrix.r3c3 - matrix.r2c3 * matrix.r3c2;
+            double r1c2 = matrix.r1c3 * matrix.r3c2 - matrix.r1c2 * matrix.r3c3;
+            double r1c3 = matrix.r1c2 * matrix.r2c3 - matrix.r1c3 * matrix.r2c2;
+
+            double r2c1 = matrix.r2c3 * matrix.r3c1 - matrix.r2c1 * matrix.r3c3;
+            double r2c2 = matrix.r1c1 * matrix.r3c3 - matrix.r1c3 * matrix.r3c1;
+            double r2c3 = matrix.r1c3 * matrix.r2c1 - matrix.r1c1 * matrix.r2c3;
+
+            double r3c1 = matrix.r2c1 * matrix.r3c2 - matrix.r2c2 * matrix.r3c1;
+            double r3c2 = matrix.r1c2 * matrix.r3c1 - matrix.r1c1 * matrix.r3c2;
+            double r3c3 = matrix.r1c1 * matrix.r2c2 - matrix.r1c2 * matrix.r2c1;
+
+            double mutiplier = 1.0 / determinant;
+
+            this.r1c1 = r1c1 * mutiplier;
+            this.r1c2 = r1c2 * mutiplier;
+            this.r1c3 = r1c3 * mutiplier;
+
+            this.r2c1 = r2c1 * mutiplier;
+            this.r2c2 = r2c2 * mutiplier;
+            this.r2c3 = r2c3 * mutiplier;
+
+            this.r3c1 = r3c1 * mutiplier;
+            this.r3c2 = r3c2 * mutiplier;
+            this.r3c3 = r3c3 * mutiplier;
+
+            return true;
+        }
+
+        public void SetRow1(double c1, double c2, double c3)
+        {
+            this.r1c1 = c1;
+            this.r1c2 = c2;
+            this.r1c3 = c3;
+        }
+
+        public void SetRow2(double c1, double c2, double c3)
+        {
+            this.r2c1 = c1;
+            this.r2c2 = c2;
+            this.r2c3 = c3;
+        }
+
+        public void SetRow3(double c1, double c2, double c3)
+        {
+            this.r3c1 = c1;
+            this.r3c2 = c2;
+            this.r3c3 = c3;
+        }
+
+        public void SetColumn1(double r1, double r2, double r3)
+        {
+            this.r1c1 = r1;
+            this.r2c1 = r2;
+            this.r3c1 = r3;
+        }
+
+        public void SetColumn2(double r1, double r2, double r3)
+        {
+            this.r1c2 = r1;
+            this.r2c2 = r2;
+            this.r3c2 = r3;
+        }
+
+        public void SetColumn3(double r1, double r2, double r3)
+        {
+            this.r1c3 = r1;
+            this.r2c3 = r2;
+            this.r3c3 = r3;
+        }
+
+        public void AppendScaled(in FP64Matrix3x3 matrix, double scale)
+        {
+            this.r1c1 += matrix.r1c1 * scale;
+            this.r1c2 += matrix.r1c2 * scale;
+            this.r1c3 += matrix.r1c3 * scale;
+
+            this.r2c1 += matrix.r2c1 * scale;
+            this.r2c2 += matrix.r2c2 * scale;
+            this.r2c3 += matrix.r2c3 * scale;
+
+            this.r3c1 += matrix.r3c1 * scale;
+            this.r3c2 += matrix.r3c2 * scale;
+            this.r3c3 += matrix.r3c3 * scale;
+        }
+
+        public static void Add(in FP64Matrix3x3 matrix1, in FP64Matrix3x3 matrix2, out FP64Matrix3x3 sum)
+        {
+            sum.r1c1 = matrix1.r1c1 + matrix2.r1c1;
+            sum.r1c2 = matrix1.r1c2 + matrix2.r1c2;
+            sum.r1c3 = matrix1.r1c3 + matrix2.r1c3;
+
+            sum.r2c1 = matrix1.r2c1 + matrix2.r2c1;
+            sum.r2c2 = matrix1.r2c2 + matrix2.r2c2;
+            sum.r2c3 = matrix1.r2c3 + matrix2.r2c3;
+
+            sum.r3c1 = matrix1.r3c1 + matrix2.r3c1;
+            sum.r3c2 = matrix1.r3c2 + matrix2.r3c2;
+            sum.r3c3 = matrix1.r3c3 + matrix2.r3c3;
+        }
+
+        public static void Subtract(in FP64Matrix3x3 minuend, in FP64Matrix3x3 subtrahend, out FP64Matrix3x3 difference)
+        {
+            difference.r1c1 = minuend.r1c1 - subtrahend.r1c1;
+            difference.r1c2 = minuend.r1c2 - subtrahend.r1c2;
+            difference.r1c3 = minuend.r1c3 - subtrahend.r1c3;
+
+            difference.r2c1 = minuend.r2c1 - subtrahend.r2c1;
+            difference.r2c2 = minuend.r2c2 - subtrahend.r2c2;
+            difference.r2c3 = minuend.r2c3 - subtrahend.r2c3;
+
+            difference.r3c1 = minuend.r3c1 - subtrahend.r3c1;
+            difference.r3c2 = minuend.r3c2 - subtrahend.r3c2;
+            difference.r3c3 = minuend.r3c3 - subtrahend.r3c3;
+        }
+
+        public static void Multiply(in FP64Matrix3x3 multiplicand, double multiplier, out FP64Matrix3x3 product)
+        {
+            product.r1c1 = multiplicand.r1c1 * multiplier;
+            product.r1c2 = multiplicand.r1c2 * multiplier;
+            product.r1c3 = multiplicand.r1c3 * multiplier;
+
+            product.r2c1 = multiplicand.r2c1 * multiplier;
+            product.r2c2 = multiplicand.r2c2 * multiplier;
+            product.r2c3 = multiplicand.r2c3 * multiplier;
+
+            product.r3c1 = multiplicand.r3c1 * multiplier;
+            product.r3c2 = multiplicand.r3c2 * multiplier;
+            product.r3c3 = multiplicand.r3c3 * multiplier;
+        }
+
+        public static void Divide(in FP64Matrix3x3 dividend, double divisor, out FP64Matrix3x3 quotient)
+        {
+            Multiply(dividend, 1.0 / divisor, out quotient);
+        }
+
+        public static void GetRightProduct(in FP64Matrix3x3 matrix, in FP64Vector3 vector, out FP64Vector3 result)
+        {
+            double x1 = matrix.r1c1 * vector.x1 + matrix.r1c2 * vector.x2 + matrix.r1c3 * vector.x3;
+            double x2 = matrix.r2c1 * vector.x1 + matrix.r2c2 * vector.x2 + matrix.r2c3 * vector.x3;
+            double x3 = matrix.r3c1 * vector.x1 + matrix.r3c2 * vector.x2 + matrix.r3c3 * vector.x3;
+
+            result.x1 = x1;
+            result.x2 = x2;
+            result.x3 = x3;
+        }
+
+        public static void GetLeftProduct(in FP64Vector3 vector, in FP64Matrix3x3 matrix, out FP64Vector3 result)
+        {
+            double x1 = vector.x1 * matrix.r1c1 + vector.x2 * matrix.r2c1 + vector.x3 * matrix.r3c1;
+            double x2 = vector.x1 * matrix.r1c2 + vector.x2 * matrix.r2c2 + vector.x3 * matrix.r3c2;
+            double x3 = vector.x1 * matrix.r1c3 + vector.x2 * matrix.r2c3 + vector.x3 * matrix.r3c3;
+
+            result.x1 = x1;
+            result.x2 = x2;
+            result.x3 = x3;
+        }
+
+        public static void LoadZero(out FP64Matrix3x3 matrix)
+        {
+            matrix.r1c1 = 0.0;
+            matrix.r1c2 = 0.0;
+            matrix.r1c3 = 0.0;
+
+            matrix.r2c1 = 0.0;
+            matrix.r2c2 = 0.0;
+            matrix.r2c3 = 0.0;
+
+            matrix.r3c1 = 0.0;
+            matrix.r3c2 = 0.0;
+            matrix.r3c3 = 0.0;
+        }
+
+        public static void LoadIdentity(out FP64Matrix3x3 matrix)
+        {
+            matrix.r1c1 = 1.0;
+            matrix.r1c2 = 0.0;
+            matrix.r1c3 = 0.0;
+
+            matrix.r2c1 = 0.0;
+            matrix.r2c2 = 1.0;
+            matrix.r2c3 = 0.0;
+
+            matrix.r3c1 = 0.0;
+            matrix.r3c2 = 0.0;
+            matrix.r3c3 = 1.0;
+        }
+
+        public static void LoadDiagonal(double d1, double d2, double d3, out FP64Matrix3x3 matrix)
+        {
+            matrix.r1c1 = d1;
+            matrix.r1c2 = 0.0;
+            matrix.r1c3 = 0.0;
+
+            matrix.r2c1 = 0.0;
+            matrix.r2c2 = d2;
+            matrix.r2c3 = 0.0;
+
+            matrix.r3c1 = 0.0;
+            matrix.r3c2 = 0.0;
+            matrix.r3c3 = d3;
+        }
+    }
+}

BasicGeometry/FP64MatrixProduct.cs

@@ -0,0 +1,173 @@
+/*
+ * Copyright 2019-2025 Andrey Pokidov <andrey.pokidov@gmail.com>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+using System;
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 11 Nov 2024
+ */
+namespace BasicGeometry
+{
+    public class FP64MatrixProduct
+    {
+        public static void Get2x2At2x2(in FP64Matrix2x2 left, in FP64Matrix2x2 right, out FP64Matrix2x2 product)
+        {
+            double r1c1 = left.r1c1 * right.r1c1 + left.r1c2 * right.r2c1;
+            double r1c2 = left.r1c1 * right.r1c2 + left.r1c2 * right.r2c2;
+
+            double r2c1 = left.r2c1 * right.r1c1 + left.r2c2 * right.r2c1;
+            double r2c2 = left.r2c1 * right.r1c2 + left.r2c2 * right.r2c2;
+
+            product.r1c1 = r1c1;
+            product.r1c2 = r1c2;
+
+            product.r2c1 = r2c1;
+            product.r2c2 = r2c2;
+        }
+
+        public static void Get2x2At3x2(in FP64Matrix2x2 left, in FP64Matrix3x2 right, out FP64Matrix3x2 product)
+        {
+            double r1c1 = left.r1c1 * right.r1c1 + left.r1c2 * right.r2c1;
+            double r1c2 = left.r1c1 * right.r1c2 + left.r1c2 * right.r2c2;
+            double r1c3 = left.r1c1 * right.r1c3 + left.r1c2 * right.r2c3;
+
+            double r2c1 = left.r2c1 * right.r1c1 + left.r2c2 * right.r2c1;
+            double r2c2 = left.r2c1 * right.r1c2 + left.r2c2 * right.r2c2;
+            double r2c3 = left.r2c1 * right.r1c3 + left.r2c2 * right.r2c3;
+
+            product.r1c1 = r1c1;
+            product.r1c2 = r1c2;
+            product.r1c3 = r1c3;
+
+            product.r2c1 = r2c1;
+            product.r2c2 = r2c2;
+            product.r2c3 = r2c3;
+        }
+    
+        public static void Get2x3At2x2(in FP64Matrix2x3 left, in FP64Matrix2x2 right, out FP64Matrix2x3 product)
+        {
+            double r1c1 = left.r1c1 * right.r1c1 + left.r1c2 * right.r2c1;
+            double r1c2 = left.r1c1 * right.r1c2 + left.r1c2 * right.r2c2;
+        
+            double r2c1 = left.r2c1 * right.r1c1 + left.r2c2 * right.r2c1;
+            double r2c2 = left.r2c1 * right.r1c2 + left.r2c2 * right.r2c2;
+        
+            double r3c1 = left.r3c1 * right.r1c1 + left.r3c2 * right.r2c1;
+            double r3c2 = left.r3c1 * right.r1c2 + left.r3c2 * right.r2c2;
+        
+            product.r1c1 = r1c1;
+            product.r1c2 = r1c2;
+        
+            product.r2c1 = r2c1;
+            product.r2c2 = r2c2;
+        
+            product.r3c1 = r3c1;
+            product.r3c2 = r3c2;
+        }
+
+        public static void Get2x3At3x2(in FP64Matrix2x3 left, in FP64Matrix3x2 right, out FP64Matrix3x3 product)
+        {
+            product.r1c1 = left.r1c1 * right.r1c1 + left.r1c2 * right.r2c1;
+            product.r1c2 = left.r1c1 * right.r1c2 + left.r1c2 * right.r2c2;
+            product.r1c3 = left.r1c1 * right.r1c3 + left.r1c2 * right.r2c3;
+        
+            product.r2c1 = left.r2c1 * right.r1c1 + left.r2c2 * right.r2c1;
+            product.r2c2 = left.r2c1 * right.r1c2 + left.r2c2 * right.r2c2;
+            product.r2c3 = left.r2c1 * right.r1c3 + left.r2c2 * right.r2c3;
+        
+            product.r3c1 = left.r3c1 * right.r1c1 + left.r3c2 * right.r2c1;
+            product.r3c2 = left.r3c1 * right.r1c2 + left.r3c2 * right.r2c2;
+            product.r3c3 = left.r3c1 * right.r1c2 + left.r3c2 * right.r2c3;
+        }
+
+        public static void Get3x2At3x3(in FP64Matrix3x2 left, in FP64Matrix3x3 right, out FP64Matrix3x2 product)
+        {
+            double r1c1 = left.r1c1 * right.r1c1 + left.r1c2 * right.r2c1 + left.r1c3 * right.r3c1;
+            double r1c2 = left.r1c1 * right.r1c2 + left.r1c2 * right.r2c2 + left.r1c3 * right.r3c2;
+            double r1c3 = left.r1c1 * right.r1c3 + left.r1c2 * right.r2c3 + left.r1c3 * right.r3c3;
+
+            double r2c1 = left.r2c1 * right.r1c1 + left.r2c2 * right.r2c1 + left.r2c3 * right.r3c1;
+            double r2c2 = left.r2c1 * right.r1c2 + left.r2c2 * right.r2c2 + left.r2c3 * right.r3c2;
+            double r2c3 = left.r2c1 * right.r1c3 + left.r2c2 * right.r2c3 + left.r2c3 * right.r3c3;
+
+            product.r1c1 = r1c1;
+            product.r1c2 = r1c2;
+            product.r1c3 = r1c3;
+
+            product.r2c1 = r2c1;
+            product.r2c2 = r2c2;
+            product.r2c3 = r2c3;
+        }
+
+        public static void Get3x2At2x3(in FP64Matrix3x2 left, in FP64Matrix2x3 right, out FP64Matrix2x2 product)
+        {
+            product.r1c1 = left.r1c1 * right.r1c1 + left.r1c2 * right.r2c1 + left.r1c3 * right.r3c1;
+            product.r1c2 = left.r1c1 * right.r1c2 + left.r1c2 * right.r2c2 + left.r1c3 * right.r3c2;
+
+            product.r2c1 = left.r2c1 * right.r1c1 + left.r2c2 * right.r2c1 + left.r2c3 * right.r3c1;
+            product.r2c2 = left.r2c1 * right.r1c2 + left.r2c2 * right.r2c2 + left.r2c3 * right.r3c2;
+        }
+
+        public static void Get3x3At2x3(in FP64Matrix3x3 left, in FP64Matrix2x3 right, out FP64Matrix2x3 product)
+        {
+            double r1c1 = left.r1c1 * right.r1c1 + left.r1c2 * right.r2c1 + left.r1c3 * right.r3c1;
+            double r1c2 = left.r1c1 * right.r1c2 + left.r1c2 * right.r2c2 + left.r1c3 * right.r3c2;
+
+            double r2c1 = left.r2c1 * right.r1c1 + left.r2c2 * right.r2c1 + left.r2c3 * right.r3c1;
+            double r2c2 = left.r2c1 * right.r1c2 + left.r2c2 * right.r2c2 + left.r2c3 * right.r3c2;
+
+            double r3c1 = left.r3c1 * right.r1c1 + left.r3c2 * right.r2c1 + left.r3c3 * right.r3c1;
+            double r3c2 = left.r3c1 * right.r1c2 + left.r3c2 * right.r2c2 + left.r3c3 * right.r3c2;
+
+            product.r1c1 = r1c1;
+            product.r1c2 = r1c2;
+
+            product.r2c1 = r2c1;
+            product.r2c2 = r2c2;
+
+            product.r3c1 = r3c1;
+            product.r3c2 = r3c2;
+        }
+
+        public static void Get3x3At3x3(in FP64Matrix3x3 left, in FP64Matrix3x3 right, out FP64Matrix3x3 product)
+        {
+            double r1c1 = left.r1c1 * right.r1c1 + left.r1c2 * right.r2c1 + left.r1c3 * right.r3c1;
+            double r1c2 = left.r1c1 * right.r1c2 + left.r1c2 * right.r2c2 + left.r1c3 * right.r3c2;
+            double r1c3 = left.r1c1 * right.r1c3 + left.r1c2 * right.r2c3 + left.r1c3 * right.r3c3;
+
+            double r2c1 = left.r2c1 * right.r1c1 + left.r2c2 * right.r2c1 + left.r2c3 * right.r3c1;
+            double r2c2 = left.r2c1 * right.r1c2 + left.r2c2 * right.r2c2 + left.r2c3 * right.r3c2;
+            double r2c3 = left.r2c1 * right.r1c3 + left.r2c2 * right.r2c3 + left.r2c3 * right.r3c3;
+
+            double r3c1 = left.r3c1 * right.r1c1 + left.r3c2 * right.r2c1 + left.r3c3 * right.r3c1;
+            double r3c2 = left.r3c1 * right.r1c2 + left.r3c2 * right.r2c2 + left.r3c3 * right.r3c2;
+            double r3c3 = left.r3c1 * right.r1c3 + left.r3c2 * right.r2c3 + left.r3c3 * right.r3c3;
+
+            product.r1c1 = r1c1;
+            product.r1c2 = r1c2;
+            product.r1c3 = r1c3;
+
+            product.r2c1 = r2c1;
+            product.r2c2 = r2c2;
+            product.r2c3 = r2c3;
+
+            product.r3c1 = r3c1;
+            product.r3c2 = r3c2;
+            product.r3c3 = r3c3;
+        }
+    }
+}

BasicGeometry/FP64Quaternion.cs

@@ -0,0 +1,204 @@
+using System;
+using System.Numerics;
+
+namespace BasicGeometry
+{
+    public struct FP64Quaternion
+    {
+        public double s0 = 0.0, x1 = 0.0, x2 = 0.0, x3 = 0.0;
+
+        public FP64Quaternion(double s0, double x1, double x2, double x3)
+        {
+            this.s0 = s0;
+            this.x1 = x1;
+            this.x2 = x2;
+            this.x3 = x3;
+        }
+
+        public FP64Quaternion(in FP32Quaternion quaternion)
+        {
+            this.s0 = quaternion.s0;
+            this.x1 = quaternion.x1;
+            this.x2 = quaternion.x2;
+            this.x3 = quaternion.x3;
+        }
+
+        public FP64Quaternion(in FP64Quaternion quaternion)
+        {
+            this.s0 = quaternion.s0;
+            this.x1 = quaternion.x1;
+            this.x2 = quaternion.x2;
+            this.x3 = quaternion.x3;
+        }
+
+        public void Reset()
+        {
+            this.s0 = 0.0;
+            this.x1 = 0.0;
+            this.x2 = 0.0;
+            this.x3 = 0.0;
+        }
+
+        public void Conjugate()
+        {
+            this.x1 = -this.x1;
+            this.x2 = -this.x2;
+            this.x3 = -this.x3;
+        }
+
+        public void SetValues(double s0, double x1, double x2, double x3)
+        {
+            this.s0 = s0;
+            this.x1 = x1;
+            this.x2 = x2;
+            this.x3 = x3;
+        }
+
+        public void SetValues(in FP32Quaternion quaternion)
+        {
+            this.s0 = quaternion.s0;
+            this.x1 = quaternion.x1;
+            this.x2 = quaternion.x2;
+            this.x3 = quaternion.x3;
+        }
+
+        public void SetValues(in FP64Quaternion quaternion)
+        {
+            this.s0 = quaternion.s0;
+            this.x1 = quaternion.x1;
+            this.x2 = quaternion.x2;
+            this.x3 = quaternion.x3;
+        }
+
+        public void SetConjugateOf(in FP64Quaternion quaternion)
+        {
+            this.s0 = quaternion.s0;
+            this.x1 = -quaternion.x1;
+            this.x2 = -quaternion.x2;
+            this.x3 = -quaternion.x3;
+        }
+
+        public readonly void MakeRotationMatrix(out FP64Matrix3x3 matrix)
+        {
+            double s0s0 = this.s0 * this.s0;
+            double x1x1 = this.x1 * this.x1;
+            double x2x2 = this.x2 * this.x2;
+            double x3x3 = this.x3 * this.x3;
+
+            double squareModule = (s0s0 + x1x1) + (x2x2 + x3x3);
+
+            if (-FP64Utility.EPSYLON <= squareModule && squareModule <= FP64Utility.EPSYLON)
+            {
+                FP64Matrix3x3.LoadIdentity(out matrix);
+                return;
+            }
+
+            double corrector1;
+            double corrector2;
+
+            if (1.0 - FP64Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0 + FP64Utility.TWO_EPSYLON) {
+                corrector1 = 2.0 - squareModule;
+                corrector2 = 2.0 * corrector1;
+            }
+            else {
+                corrector1 = 1.0 / squareModule;
+                corrector2 = 2.0 / squareModule;
+            }
+
+            double s0x1 = this.s0 * this.x1;
+            double s0x2 = this.s0 * this.x2;
+            double s0x3 = this.s0 * this.x3;
+            double x1x2 = this.x1 * this.x2;
+            double x1x3 = this.x1 * this.x3;
+            double x2x3 = this.x2 * this.x3;
+
+            matrix.r1c1 = corrector1 * ((s0s0 + x1x1) - (x2x2 + x3x3));
+            matrix.r2c2 = corrector1 * ((s0s0 + x2x2) - (x1x1 + x3x3));
+            matrix.r3c3 = corrector1 * ((s0s0 + x3x3) - (x1x1 + x2x2));
+
+            matrix.r1c2 = corrector2 * (x1x2 - s0x3);
+            matrix.r2c3 = corrector2 * (x2x3 - s0x1);
+            matrix.r3c1 = corrector2 * (x1x3 - s0x2);
+
+            matrix.r2c1 = corrector2 * (x1x2 + s0x3);
+            matrix.r3c2 = corrector2 * (x2x3 + s0x1);
+            matrix.r1c3 = corrector2 * (x1x3 + s0x2);
+        }
+
+        public readonly void MakeReverseMatrix(out FP64Matrix3x3 matrix)
+        {
+            double s0s0 = this.s0 * this.s0;
+            double x1x1 = this.x1 * this.x1;
+            double x2x2 = this.x2 * this.x2;
+            double x3x3 = this.x3 * this.x3;
+
+            double squareModule = (s0s0 + x1x1) + (x2x2 + x3x3);
+
+            if (-FP64Utility.EPSYLON <= squareModule && squareModule <= FP64Utility.EPSYLON)
+            {
+                FP64Matrix3x3.LoadIdentity(out matrix);
+                return;
+            }
+
+            double corrector1;
+            double corrector2;
+
+            if (1.0 - FP64Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0 + FP64Utility.TWO_EPSYLON) {
+                corrector1 = 2.0 - squareModule;
+                corrector2 = 2.0 * corrector1;
+            }
+            else {
+                corrector1 = 1.0 / squareModule;
+                corrector2 = 2.0 / squareModule;
+            }
+
+            double s0x1 = this.s0 * this.x1;
+            double s0x2 = this.s0 * this.x2;
+            double s0x3 = this.s0 * this.x3;
+            double x1x2 = this.x1 * this.x2;
+            double x1x3 = this.x1 * this.x3;
+            double x2x3 = this.x2 * this.x3;
+
+            matrix.r1c1 = corrector1 * ((s0s0 + x1x1) - (x2x2 + x3x3));
+            matrix.r2c2 = corrector1 * ((s0s0 + x2x2) - (x1x1 + x3x3));
+            matrix.r3c3 = corrector1 * ((s0s0 + x3x3) - (x1x1 + x2x2));
+
+            matrix.r1c2 = corrector2 * (x1x2 + s0x3);
+            matrix.r2c3 = corrector2 * (x2x3 + s0x1);
+            matrix.r3c1 = corrector2 * (x1x3 + s0x2);
+
+            matrix.r2c1 = corrector2 * (x1x2 - s0x3);
+            matrix.r3c2 = corrector2 * (x2x3 - s0x1);
+            matrix.r1c3 = corrector2 * (x1x3 - s0x2);
+        }
+
+        public static void Add(in FP64Quaternion quaternion1, in FP64Quaternion quaternion2, out FP64Quaternion sum)
+        {
+            sum.s0 = quaternion1.s0 + quaternion2.s0;
+            sum.x1 = quaternion1.x1 + quaternion2.x1;
+            sum.x2 = quaternion1.x2 + quaternion2.x2;
+            sum.x3 = quaternion1.x3 + quaternion2.x3;
+        }
+
+        public static void Subtract(in FP64Quaternion minuend, in FP64Quaternion subtrahend, out FP64Quaternion difference)
+        {
+            difference.s0 = minuend.s0 - subtrahend.s0;
+            difference.x1 = minuend.x1 - subtrahend.x1;
+            difference.x2 = minuend.x2 - subtrahend.x2;
+            difference.x3 = minuend.x3 - subtrahend.x3;
+        }
+
+        public static void Multiply(in FP64Quaternion left, in FP64Quaternion right, out FP64Quaternion product)
+        {
+            double s0 = (left.s0 * right.s0 - left.x1 * right.x1) - (left.x2 * right.x2 + left.x3 * right.x3);
+            double x1 = (left.x1 * right.s0 + left.s0 * right.x1) - (left.x3 * right.x2 - left.x2 * right.x3);
+            double x2 = (left.x2 * right.s0 + left.s0 * right.x2) - (left.x1 * right.x3 - left.x3 * right.x1);
+            double x3 = (left.x3 * right.s0 + left.s0 * right.x3) - (left.x2 * right.x1 - left.x1 * right.x2);
+
+            product.s0 = s0;
+            product.x1 = x1;
+            product.x2 = x2;
+            product.x3 = x3;
+        }
+    }
+}

BasicGeometry/FP64Radians.cs

@@ -0,0 +1,75 @@
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 18 Nov 2024
+ */
+
+namespace BasicGeometry
+{
+    public class FP64Radians
+    {
+        public const double PI = 3.14159265358979324;
+        public const double TWO_PI = 6.28318530717958648;
+        public const double HALF_OF_PI = 1.57079632679489662;
+        public const double THIRD_OF_PI = 1.04719755119659775;
+        public const double FOURTH_OF_PI = 0.78539816339744831;
+        public const double SIXTH_OF_PI = 0.523598775598298873;
+
+        public const double DEGREES_IN_RADIAN = 57.2957795130823209;
+        public const double TURNS_IN_RADIAN = 0.159154943091895336;
+
+        public static double ToDegrees(double radians)
+        {
+            return radians * DEGREES_IN_RADIAN;
+        }
+
+        public static double ToTurns(double radians)
+        {
+            return radians * TURNS_IN_RADIAN;
+        }
+
+        public static double ToUnits(double radians, AngleUnit toUnit)
+        {
+            if (toUnit == AngleUnit.DEGREES)
+            {
+                return radians * DEGREES_IN_RADIAN;
+            }
+
+            if (toUnit == AngleUnit.TURNS)
+            {
+                return radians * TURNS_IN_RADIAN;
+            }
+
+            return radians;
+        }
+
+        public static double Normalize(double radians, AngleRange range)
+        {
+            if (range == AngleRange.UNSIGNED_RANGE)
+            {
+                if (0.0 <= radians && radians < TWO_PI)
+                {
+                    return radians;
+                }
+            }
+            else
+            {
+                if (-PI < radians && radians <= PI)
+                {
+                    return radians;
+                }
+            }
+
+            double turns = radians * TURNS_IN_RADIAN;
+
+            turns -= Math.Floor(turns);
+
+            if (range == AngleRange.SIGNED_RANGE && turns > 0.5)
+            {
+                turns -= 1.0;
+            }
+
+            return turns * TWO_PI;
+        }
+    }
+}

BasicGeometry/FP64Turns.cs

@@ -0,0 +1,63 @@
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 18 Nov 2024
+ */
+
+namespace BasicGeometry
+{
+    public class FP64Turns
+    {
+        public static double TurnsToRadians(double turns)
+        {
+            return turns * FP64Radians.TWO_PI;
+        }
+
+        public static double TurnsToDegrees(double turns)
+        {
+            return turns * 360.0;
+        }
+
+        public static double ToUnits(double turns, AngleUnit toUnit)
+        {
+            if (toUnit == AngleUnit.RADIANS)
+            {
+                return turns * FP64Radians.TWO_PI;
+            }
+
+            if (toUnit == AngleUnit.DEGREES)
+            {
+                return turns * 360.0;
+            }
+
+            return turns;
+        }
+
+        public static double Normalize(double turns, AngleRange range)
+        {
+            if (range == AngleRange.UNSIGNED_RANGE)
+            {
+                if (0.0 <= turns && turns < 1.0)
+                {
+                    return turns;
+                }
+            }
+            else
+            {
+                if (-0.5 < turns && turns <= 0.5)
+                {
+                    return turns;
+                }
+            }
+
+            double rest = turns - Math.Floor(turns);
+
+            if (range == AngleRange.SIGNED_RANGE && rest > 0.5)
+            {
+                rest -= 1.0;
+            }
+
+            return rest;
+        }
+    }
+}

BasicGeometry/FP64Utility.cs

@@ -0,0 +1,20 @@
+using System;
+
+namespace BasicGeometry
+{
+    public class FP64Utility
+    {
+        public const double EPSYLON = 5E-14;
+        public const double TWO_EPSYLON = 1E-13;
+        public const double SQUARE_EPSYLON = 2.5E-27;
+        
+        public const double EPSYLON_EFFECTIVENESS_LIMIT = 1.0;
+
+        public const double ONE_THIRD = 0.333333333333333333;
+        public const double ONE_SIXTH = 0.166666666666666667;
+        public const double ONE_NINETH = 0.111111111111111111;
+
+        public const double GOLDEN_RATIO_HIGH = 1.61803398874989485;
+        public const double GOLDEN_RATIO_LOW = 0.61803398874989485;
+    }
+}

BasicGeometry/FP64Vector2.cs

@@ -0,0 +1,257 @@
+/*
+ * Copyright 2019-2025 Andrey Pokidov <andrey.pokidov@gmail.com>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+using System;
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 1 Feb 2019
+ */
+
+namespace BasicGeometry
+{
+    public struct FP64Vector2
+    {
+        public static readonly FP64Vector2 ZERO = new FP64Vector2(0.0, 0.0);
+
+        public double x1 = 0.0;
+        public double x2 = 0.0;
+
+        public FP64Vector2(double x1, double x2)
+        {
+            this.x1 = x1;
+            this.x2 = x2;
+        }
+
+        public FP64Vector2(in FP64Vector2 vector)
+        {
+            this.x1 = vector.x1;
+            this.x2 = vector.x2;
+        }
+
+        public FP64Vector2(in FP32Vector2 vector)
+        {
+            this.x1 = vector.x1;
+            this.x2 = vector.x2;
+        }
+
+        public readonly double GetSquareModule()
+        {
+            return this.x1 * this.x1 + this.x2 * this.x2;
+        }
+
+        public readonly double GetModule()
+        {
+            return Math.Sqrt(this.GetSquareModule());
+        }
+
+        public int Normalize()
+        {
+            double squareModule = this.GetSquareModule();
+
+            if (1.0 - FP64Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0 + FP64Utility.TWO_EPSYLON)
+            {
+                return 1;
+            }
+
+            if (squareModule <= FP64Utility.SQUARE_EPSYLON)
+            {
+                this.Reset();
+                return 0;
+            }
+
+            double module = Math.Sqrt(squareModule);
+
+            this.x1 /= module;
+            this.x2 /= module;
+
+            return 1;
+        }
+
+        public void Reverse()
+        {
+            this.x1 = -this.x1;
+            this.x2 = -this.x2;
+        }
+
+        public readonly bool IsZero()
+        {
+            return this.GetSquareModule() <= FP64Utility.SQUARE_EPSYLON;
+        }
+
+        public readonly bool IsUnit()
+        {
+            double squareModule = this.GetSquareModule();
+            return 1.0 - FP64Utility.TWO_EPSYLON <= squareModule && squareModule <= FP64Utility.EPSYLON;
+        }
+
+        public void Reset()
+        {
+            this.x1 = 0.0;
+            this.x2 = 0.0;
+        }
+
+        public void SetValues(double x1, double x2)
+        {
+            this.x1 = x1;
+            this.x2 = x2;
+        }
+
+        public void SetValues(in FP64Vector2 vector)
+        {
+            this.x1 = vector.x1;
+            this.x2 = vector.x2;
+        }
+
+        public void SetValues(in FP32Vector2 vector)
+        {
+            this.x1 = vector.x1;
+            this.x2 = vector.x2;
+        }
+
+        public void SetReverseOf(in FP64Vector2 vector)
+        {
+            this.x1 = -vector.x1;
+            this.x2 = -vector.x2;
+        }
+
+        public void SetReverseOf(in FP32Vector2 vector)
+        {
+            this.x1 = -vector.x1;
+            this.x2 = -vector.x2;
+        }
+
+        public void AppendScaled(FP64Vector2 summand, double scale)
+        {
+            this.x1 += summand.x1 * scale;
+            this.x2 += summand.x2 * scale;
+        }
+
+        public readonly override string ToString()
+        {
+            return String.Format("DPVector2({0}, {1})", this.x1, this.x2);
+        }
+
+        public static void Add(in FP64Vector2 vector1, in FP64Vector2 vector2, out FP64Vector2 sum)
+        {
+            sum.x1 = vector1.x1 + vector2.x1;
+            sum.x2 = vector1.x2 + vector2.x2;
+        }
+
+        public static void Subtract(in FP64Vector2 minuend, in FP64Vector2 subtrahend, out FP64Vector2 difference)
+        {
+            difference.x1 = minuend.x1 - subtrahend.x1;
+            difference.x2 = minuend.x2 - subtrahend.x2;
+        }
+
+        public static void Multiply(in FP64Vector2 multiplicand, double multiplier, out FP64Vector2 product)
+        {
+            product.x1 = multiplicand.x1 * multiplier;
+            product.x2 = multiplicand.x2 * multiplier;
+        }
+
+        public static void Divide(in FP64Vector2 dividend, double divisor, out FP64Vector2 quotient)
+        {
+            Multiply(dividend, 1.0 / divisor, out quotient);
+        }
+
+        public static void GetMean2(in FP64Vector2 vector1, in FP64Vector2 vector2, out FP64Vector2 result)
+        {
+            result.x1 = (vector1.x1 + vector2.x1) * 0.5;
+            result.x2 = (vector1.x2 + vector2.x2) * 0.5;
+        }
+
+        public static void GetMean3(in FP64Vector2 vector1, in FP64Vector2 vector2, in FP64Vector2 vector3, out FP64Vector2 result)
+        {
+            result.x1 = (vector1.x1 + vector2.x1 + vector3.x1) * FP64Utility.ONE_THIRD;
+            result.x2 = (vector1.x2 + vector2.x2 + vector3.x2) * FP64Utility.ONE_THIRD;
+        }
+
+        public static double GetScalarProduct(in FP64Vector2 vector1, in FP64Vector2 vector2)
+        {
+            return vector1.x1 * vector2.x1 + vector1.x2 * vector2.x2;
+        }
+ 
+        public static double GetCrossProduct(in FP64Vector2 vector1, in FP64Vector2 vector2)
+        {
+            return vector1.x1 * vector2.x2 - vector1.x2 * vector2.x1;
+        }
+
+        public static double GetAngle(in FP64Vector2 vector1, in FP64Vector2 vector2, AngleUnit unit)
+        {
+            double squareModule1 = vector1.GetSquareModule();
+
+            if (squareModule1 <= FP64Utility.SQUARE_EPSYLON)
+            {
+                return 0.0;
+            }
+
+            double squareModule2 = vector2.GetSquareModule();
+
+            if (squareModule2 <= FP64Utility.SQUARE_EPSYLON)
+            {
+                return 0.0;
+            }
+
+            double cosine = FP64Vector2.GetScalarProduct(vector1, vector2) / Math.Sqrt(squareModule1 * squareModule2);
+
+            if (1.0 - FP64Utility.EPSYLON <= cosine)
+            {
+                return 0.0;
+            }
+
+            if (cosine <= -(1.0 - FP64Utility.EPSYLON))
+            {
+                return FP64Angle.GetHalfCircle(unit);
+            }
+    
+            return FP64Radians.ToUnits(Math.Acos(cosine), unit);
+        }
+
+        public static double GetSquareDistance(in FP64Vector2 vector1, in FP64Vector2 vector2)
+        {
+            double dx1 = vector1.x1 - vector2.x1;
+            double dx2 = vector1.x2 - vector2.x2;
+        
+            return dx1 * dx1 + dx2 * dx2;
+        }
+    
+        public static double GetDistance(in FP64Vector2 vector1, in FP64Vector2 vector2)
+        {
+            return Math.Sqrt(GetSquareDistance(vector1, vector2));
+        }
+    
+        public static bool AreEqual(in FP64Vector2 vector1, in FP64Vector2 vector2)
+        {
+            double squareModule1 = vector1.GetSquareModule();
+            double squareModule2 = vector2.GetSquareModule();
+            double squareModule3 = GetSquareDistance(vector1, vector2);
+
+            // 2.0 means dimension amount
+            if (squareModule1 < FP64Utility.EPSYLON_EFFECTIVENESS_LIMIT || squareModule2 < FP64Utility.EPSYLON_EFFECTIVENESS_LIMIT)
+            {
+                return squareModule3 < (2.0 * FP64Utility.SQUARE_EPSYLON);
+            }
+
+            if (squareModule1 <= squareModule2)
+            {
+                return squareModule3 <= (2.0 * FP64Utility.SQUARE_EPSYLON) * squareModule2;
+            }
+
+            return squareModule3 <= (2.0 * FP64Utility.SQUARE_EPSYLON) * squareModule1;
+        }
+    }
+}

BasicGeometry/FP64Vector3.cs

@@ -0,0 +1,300 @@
+/*
+ * Copyright 2019-2025 Andrey Pokidov <andrey.pokidov@gmail.com>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+using System;
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 1 Feb 2019
+ */
+
+namespace BasicGeometry
+{
+    public struct FP64Vector3
+    {
+        public static readonly FP64Vector3 ZERO = new FP64Vector3(0.0, 0.0, 0.0);
+
+        public double x1 = 0.0;
+        public double x2 = 0.0;
+        public double x3 = 0.0;
+
+        public FP64Vector3(double x1, double x2, double x3)
+        {
+            this.x1 = x1;
+            this.x2 = x2;
+            this.x3 = x3;
+        }
+
+        public FP64Vector3(in FP64Vector3 vector)
+        {
+            this.x1 = vector.x1;
+            this.x2 = vector.x2;
+            this.x3 = vector.x3;
+        }
+
+        public FP64Vector3(in FP32Vector3 vector)
+        {
+            this.x1 = vector.x1;
+            this.x2 = vector.x2;
+            this.x3 = vector.x3;
+        }
+
+        public readonly double GetSquareModule()
+        {
+            return this.x1 * this.x1 + this.x2 * this.x2 + this.x3 * this.x3;
+        }
+
+        public readonly double GetModule()
+        {
+            return Math.Sqrt(this.GetSquareModule());
+        }
+
+        public int Normalize()
+        {
+            double squareModule = this.GetSquareModule();
+
+            if (1.0 - FP64Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0 + FP64Utility.TWO_EPSYLON)
+            {
+                return 1;
+            }
+
+            if (squareModule <= FP64Utility.SQUARE_EPSYLON)
+            {
+                this.Reset();
+                return 0;
+            }
+
+            double module = Math.Sqrt(squareModule);
+
+            this.x1 /= module;
+            this.x2 /= module;
+            this.x3 /= module;
+
+            return 1;
+        }
+
+        public void Reverse()
+        {
+            this.x1 = -this.x1;
+            this.x2 = -this.x2;
+            this.x3 = -this.x3;
+        }
+
+        public readonly bool IsZero()
+        {
+            return this.GetSquareModule() <= FP64Utility.SQUARE_EPSYLON;
+        }
+
+        public readonly bool IsUnit()
+        {
+            double squareModule = this.GetSquareModule();
+            return 1.0 - FP64Utility.TWO_EPSYLON <= squareModule && squareModule <= FP64Utility.EPSYLON;
+        }
+
+        public void Reset()
+        {
+            this.x1 = 0.0;
+            this.x2 = 0.0;
+            this.x3 = 0.0;
+        }
+
+        public void SetValues(double x1, double x2, double x3)
+        {
+            this.x1 = x1;
+            this.x2 = x2;
+            this.x3 = x3;
+        }
+
+        public void SetValues(in FP64Vector3 vector)
+        {
+            this.x1 = vector.x1;
+            this.x2 = vector.x2;
+            this.x3 = vector.x3;
+        }
+
+        public void SetValues(in FP32Vector3 vector)
+        {
+            this.x1 = vector.x1;
+            this.x2 = vector.x2;
+            this.x3 = vector.x3;
+        }
+
+        public void SetReverseOf(in FP64Vector3 vector)
+        {
+            this.x1 = -vector.x1;
+            this.x2 = -vector.x2;
+            this.x3 = -vector.x3;
+        }
+
+        public void SetReverseOf(in FP32Vector3 vector)
+        {
+            this.x1 = -vector.x1;
+            this.x2 = -vector.x2;
+            this.x3 = -vector.x3;
+        }
+
+        public void AppendScaled(FP64Vector3 summand, double scale)
+        {
+            this.x1 += summand.x1 * scale;
+            this.x2 += summand.x2 * scale;
+            this.x3 += summand.x3 * scale;
+        }
+
+        public readonly override string ToString()
+        {
+            return String.Format("DPVector3({0}, {1}, {2})", this.x1, this.x2, this.x3);
+        }
+
+        public static void Add(in FP64Vector3 vector1, in FP64Vector3 vector2, out FP64Vector3 sum)
+        {
+            sum.x1 = vector1.x1 + vector2.x1;
+            sum.x2 = vector1.x2 + vector2.x2;
+            sum.x3 = vector1.x3 + vector2.x3;
+        }
+
+        public static void Subtract(in FP64Vector3 minuend, in FP64Vector3 subtrahend, out FP64Vector3 difference)
+        {
+            difference.x1 = minuend.x1 - subtrahend.x1;
+            difference.x2 = minuend.x2 - subtrahend.x2;
+            difference.x3 = minuend.x3 - subtrahend.x3;
+        }
+
+        public static void Multiply(in FP64Vector3 multiplicand, double multiplier, out FP64Vector3 product)
+        {
+            product.x1 = multiplicand.x1 * multiplier;
+            product.x2 = multiplicand.x2 * multiplier;
+            product.x3 = multiplicand.x3 * multiplier;
+        }
+
+        public static void Divide(in FP64Vector3 dividend, double divisor, out FP64Vector3 quotient)
+        {
+            Multiply(dividend, 1.0 / divisor, out quotient);
+        }
+
+        public static void GetMean2(in FP64Vector3 vector1, in FP64Vector3 vector2, out FP64Vector3 result)
+        {
+            result.x1 = (vector1.x1 + vector2.x1) * 0.5;
+            result.x2 = (vector1.x2 + vector2.x2) * 0.5;
+            result.x3 = (vector1.x3 + vector2.x3) * 0.5;
+        }
+
+        public static void GetMean3(in FP64Vector3 vector1, in FP64Vector3 vector2, in FP64Vector3 vector3, out FP64Vector3 result)
+        {
+            result.x1 = (vector1.x1 + vector2.x1 + vector3.x1)  * FP64Utility.ONE_THIRD;
+            result.x2 = (vector1.x2 + vector2.x2 + vector3.x2)  * FP64Utility.ONE_THIRD;
+            result.x3 = (vector1.x3 + vector2.x3 + vector3.x3)  * FP64Utility.ONE_THIRD;
+        }
+
+        public static double GetScalarProduct(in FP64Vector3 vector1, in FP64Vector3 vector2)
+        {
+            return vector1.x1 * vector2.x1 + vector1.x2 * vector2.x2 + vector1.x3 * vector2.x3;
+        }
+ 
+        public static void GetCrossProduct(in FP64Vector3 vector1, in FP64Vector3 vector2, out FP64Vector3 result)
+        {
+            double x1 = vector1.x2 * vector2.x3 - vector1.x3 * vector2.x2;
+            double x2 = vector1.x3 * vector2.x1 - vector1.x1 * vector2.x3;
+            double x3 = vector1.x1 * vector2.x2 - vector1.x2 * vector2.x1;
+
+            result.x1 = x1;
+            result.x2 = x2;
+            result.x3 = x3;
+        }
+
+        public static double GetTripleProduct(in FP64Vector3 vector1, in FP64Vector3 vector2, in FP64Vector3 vector3)
+        {
+            return vector1.x1 * (vector2.x2 * vector3.x3 - vector2.x3 * vector3.x2)
+                 + vector1.x2 * (vector2.x3 * vector3.x1 - vector2.x1 * vector3.x3)
+                 + vector1.x3 * (vector2.x1 * vector3.x2 - vector2.x2 * vector3.x1);
+        }
+ 
+        public static void GetDoubleCrossProduct(in FP64Vector3 vector1, in FP64Vector3 vector2, in FP64Vector3 vector3, out FP64Vector3 result)
+        {
+            // [a x [b x c]] = b * (a, c) - c * (a, b)
+            double ac = GetScalarProduct(vector1, vector3);
+            double ab = GetScalarProduct(vector1, vector2);
+
+            result.x1 = ac * vector2.x1 - ab * vector3.x1;
+            result.x2 = ac * vector2.x2 - ab * vector3.x2;
+            result.x3 = ac * vector2.x3 - ab * vector3.x3;
+        }
+
+        public static double GetAngle(in FP64Vector3 vector1, in FP64Vector3 vector2, AngleUnit unit)
+        {
+            double squareModule1 = vector1.GetSquareModule();
+
+            if (squareModule1 <= FP64Utility.SQUARE_EPSYLON)
+            {
+                return 0.0;
+            }
+
+            double squareModule2 = vector2.GetSquareModule();
+
+            if (squareModule2 <= FP64Utility.SQUARE_EPSYLON)
+            {
+                return 0.0;
+            }
+
+            double cosine = FP64Vector3.GetScalarProduct(vector1, vector2) / Math.Sqrt(squareModule1 * squareModule2);
+
+            if (1.0 - FP64Utility.EPSYLON <= cosine)
+            {
+                return 0.0;
+            }
+
+            if (cosine <= -(1.0 - FP64Utility.EPSYLON))
+            {
+                return FP64Angle.GetHalfCircle(unit);
+            }
+    
+            return FP64Radians.ToUnits(Math.Acos(cosine), unit);
+        }
+    
+        public static double GetSquareDistance(in FP64Vector3 vector1, in FP64Vector3 vector2)
+        {
+            double dx1 = vector1.x1 - vector2.x1;
+            double dx2 = vector1.x2 - vector2.x2;
+            double dx3 = vector1.x3 - vector2.x3;
+        
+            return dx1 * dx1 + dx2 * dx2 + dx3 * dx3;
+        }
+    
+        public static double GetDistance(in FP64Vector3 vector1, in FP64Vector3 vector2)
+        {
+            return Math.Sqrt(GetSquareDistance(vector1, vector2));
+        }
+    
+        public static bool AreEqual(in FP64Vector3 vector1, in FP64Vector3 vector2)
+        {
+            double squareModule1 = vector1.GetSquareModule();
+            double squareModule2 = vector2.GetSquareModule();
+            double squareModule3 = GetSquareDistance(vector1, vector2);
+
+            // 3.0 means dimension amount
+            if (squareModule1 < FP64Utility.EPSYLON_EFFECTIVENESS_LIMIT || squareModule2 < FP64Utility.EPSYLON_EFFECTIVENESS_LIMIT)
+            {
+                return squareModule3 < (3.0 * FP64Utility.SQUARE_EPSYLON);
+            }
+
+            if (squareModule1 <= squareModule2)
+            {
+                return squareModule3 <= (3.0 * FP64Utility.SQUARE_EPSYLON) * squareModule2;
+            }
+
+            return squareModule3 <= (3.0 * FP64Utility.SQUARE_EPSYLON) * squareModule1;
+        }
+    }
+}

BasicGeometry/FP64Versor.cs

@@ -0,0 +1,312 @@
+/*
+ * Copyright 2019-2025 Andrey Pokidov <andrey.pokidov@gmail.com>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 20 Oct 2024
+ */
+
+namespace BasicGeometry
+{
+    public struct FP64Versor
+    {
+        private double s0 = 1.0;
+        private double x1 = 0.0;
+        private double x2 = 0.0;
+        private double x3 = 0.0;
+
+        public FP64Versor(double s0, double x1, double x2, double x3)
+        {
+            LoadValues(s0, x1, x2, x3, out this);
+        }
+
+        public FP64Versor(in FP64Versor versor)
+        {
+            this.s0 = versor.s0;
+            this.x1 = versor.x1;
+            this.x2 = versor.x2;
+            this.x3 = versor.x3;
+        }
+
+        public FP64Versor(in FP32Versor versor)
+        {
+            LoadValues(
+                versor.GetScalar(),
+                versor.GetX1(),
+                versor.GetX2(),
+                versor.GetX3(),
+                out this);
+        }
+
+        public readonly double GetScalar()
+        {
+            return this.s0;
+        }
+
+        public readonly double GetX1()
+        {
+            return this.x1;
+        }
+
+        public readonly double GetX2()
+        {
+            return this.x2;
+        }
+
+        public readonly double GetX3()
+        {
+            return this.x3;
+        }
+
+        public readonly bool IsIdle()
+        {
+            return this.s0 <= -(1.0 - FP64Utility.EPSYLON) || (1.0 - FP64Utility.EPSYLON) <= this.s0;
+        }
+
+        public void Reset()
+        {
+            this.s0 = 1.0;
+            this.x1 = 0.0;
+            this.x2 = 0.0;
+            this.x3 = 0.0;
+        }
+
+        public void Invert()
+        {
+            this.x1 = -this.x1;
+            this.x2 = -this.x2;
+            this.x3 = -this.x3;
+        }
+
+        public readonly double GetAngle(AngleUnit unit)
+        {
+            if (this.s0 <= -(1.0 - FP64Utility.TWO_EPSYLON) || 1.0 - FP64Utility.TWO_EPSYLON <= this.s0) {
+                return 0.0;
+            }
+
+            if (-FP64Utility.EPSYLON <= this.s0 && this.s0 <= FP64Utility.EPSYLON)
+            {
+                return FP64Angle.GetHalfCircle(unit);
+            }
+
+            return FP64Radians.ToUnits(2.0 * Math.Acos(s0), unit);
+        }
+
+        public readonly void MakeRotationMatrix(out FP64Matrix3x3 matrix)
+        {
+            double s0s0 = this.s0 * this.s0;
+            double x1x1 = this.x1 * this.x1;
+            double x2x2 = this.x1 * this.x2;
+            double x3x3 = this.x1 * this.x3;
+
+
+            double s0x1 = 2.0 * this.s0 * this.x1;
+            double s0x2 = 2.0 * this.s0 * this.x2;
+            double s0x3 = 2.0 * this.s0 * this.x3;
+
+            double x1x2 = 2.0 * this.x1 * this.x2;
+            double x1x3 = 2.0 * this.x1 * this.x3;
+            double x2x3 = 2.0 * this.x2 * this.x3;
+
+            matrix.r1c1 = (s0s0 + x1x1) - (x2x2 + x3x3);
+            matrix.r2c2 = (s0s0 + x2x2) - (x1x1 + x3x3);
+            matrix.r3c3 = (s0s0 + x3x3) - (x1x1 + x2x2);
+
+            matrix.r1c2 = x1x2 - s0x3;
+            matrix.r2c3 = x2x3 - s0x1;
+            matrix.r3c1 = x1x3 - s0x2;
+
+            matrix.r2c1 = x1x2 + s0x3;
+            matrix.r3c2 = x2x3 + s0x1;
+            matrix.r1c3 = x1x3 + s0x2;
+        }
+
+        public readonly void MakeReverseMatrix(out FP64Matrix3x3 matrix)
+        {
+            double s0s0 = this.s0 * this.s0;
+            double x1x1 = this.x1 * this.x1;
+            double x2x2 = this.x1 * this.x2;
+            double x3x3 = this.x1 * this.x3;
+
+            double s0x1 = 2.0 * this.s0 * this.x1;
+            double s0x2 = 2.0 * this.s0 * this.x2;
+            double s0x3 = 2.0 * this.s0 * this.x3;
+
+            double x1x2 = 2.0 * this.x1 * this.x2;
+            double x1x3 = 2.0 * this.x1 * this.x3;
+            double x2x3 = 2.0 * this.x2 * this.x3;
+
+            matrix.r1c1 = (s0s0 + x1x1) - (x2x2 + x3x3);
+            matrix.r2c2 = (s0s0 + x2x2) - (x1x1 + x3x3);
+            matrix.r3c3 = (s0s0 + x3x3) - (x1x1 + x2x2);
+
+            matrix.r1c2 = x1x2 + s0x3;
+            matrix.r2c3 = x2x3 + s0x1;
+            matrix.r3c1 = x1x3 + s0x2;
+
+            matrix.r2c1 = x1x2 - s0x3;
+            matrix.r3c2 = x2x3 - s0x1;
+            matrix.r1c3 = x1x3 - s0x2;
+        }
+
+        public void SetValues(double s0, double x1, double x2, double x3)
+        {
+            this.s0 = s0;
+            this.x1 = x1;
+            this.x2 = x2;
+            this.x3 = x3;
+            
+            double squareModule = (s0 * s0 + x1 * x1) + (x2 * x2 + x3 * x3);
+            
+            if (1.0 - FP64Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0 + FP64Utility.TWO_EPSYLON)
+            {
+                return;
+            }
+
+            this.Normalize(squareModule);
+        }
+
+        public void SetValues(in FP64Versor versor)
+        {
+            this.s0 = versor.s0;
+            this.x1 = versor.x1;
+            this.x2 = versor.x2;
+            this.x3 = versor.x3;
+        }
+
+        public void SetValues(in FP32Versor versor)
+        {
+            LoadValues(
+                versor.GetScalar(),
+                versor.GetX1(),
+                versor.GetX2(),
+                versor.GetX3(),
+                out this);
+        }
+
+        public void SetInverted(in FP64Versor versor)
+        {
+            this.s0 = versor.s0;
+            this.x1 = -versor.x1;
+            this.x2 = -versor.x2;
+            this.x3 = -versor.x3;
+        }
+
+        public void SetInverted(in FP32Versor versor)
+        {
+            LoadValues(
+                versor.GetScalar(),
+                versor.GetX1(),
+                versor.GetX2(),
+                versor.GetX3(),
+                out this);
+        }
+
+        public readonly void Turn(in FP64Vector3 vector, out FP64Vector3 result)
+        {
+            double tx1 = 2.0 * (this.x2 * vector.x3 - this.x3 * vector.x2);
+            double tx2 = 2.0 * (this.x3 * vector.x1 - this.x1 * vector.x3);
+            double tx3 = 2.0 * (this.x1 * vector.x2 - this.x2 * vector.x1);
+
+            double x1 = (vector.x1 + tx1 * this.s0) + (this.x2 * tx3 - this.x3 * tx2);
+            double x2 = (vector.x2 + tx2 * this.s0) + (this.x3 * tx1 - this.x1 * tx3);
+            double x3 = (vector.x3 + tx3 * this.s0) + (this.x1 * tx2 - this.x2 * tx1);
+
+            result.x1 = x1;
+            result.x2 = x2;
+            result.x3 = x3;
+        }
+
+        public readonly void TurnBack(in FP64Vector3 vector, out FP64Vector3 result)
+        {
+            double tx1 = 2.0 * (this.x2 * vector.x3 - this.x3 * vector.x2);
+            double tx2 = 2.0 * (this.x3 * vector.x1 - this.x1 * vector.x3);
+            double tx3 = 2.0 * (this.x1 * vector.x2 - this.x2 * vector.x1);
+
+            double x1 = (vector.x1 - tx1 * this.s0) + (this.x2 * tx3 - this.x3 * tx2);
+            double x2 = (vector.x2 - tx2 * this.s0) + (this.x3 * tx1 - this.x1 * tx3);
+            double x3 = (vector.x3 - tx3 * this.s0) + (this.x1 * tx2 - this.x2 * tx1);
+
+            result.x1 = x1;
+            result.x2 = x2;
+            result.x3 = x3;
+        }
+
+        private void Normalize(double squareModule)
+        {
+            if (squareModule <= FP64Utility.SQUARE_EPSYLON || (this.x1 * this.x1 + this.x2 * this.x2 + this.x3 * this.x3) <= FP64Utility.SQUARE_EPSYLON * squareModule)
+            {
+                this.Reset();
+                return;
+            }
+
+            double module = Math.Sqrt(squareModule);
+
+            this.s0 /= module;
+            this.x1 /= module;
+            this.x2 /= module;
+            this.x3 /= module;
+        }
+
+        public static void Combine(in FP64Versor second, in FP64Versor first, out FP64Versor result)
+        {
+            double s0 = (second.s0 * first.s0 - second.x1 * first.x1) - (second.x2 * first.x2 + second.x3 * first.x3);
+            double x1 = (second.x1 * first.s0 + second.s0 * first.x1) - (second.x3 * first.x2 - second.x2 * first.x3);
+            double x2 = (second.x2 * first.s0 + second.s0 * first.x2) - (second.x1 * first.x3 - second.x3 * first.x1);
+            double x3 = (second.x3 * first.s0 + second.s0 * first.x3) - (second.x2 * first.x1 - second.x1 * first.x2);
+
+            double squareModule = (s0 * s0 + x1 * x1) + (x2 * x2 + x3 * x3);
+
+            result.s0 = s0;
+            result.x1 = x1;
+            result.x2 = x2;
+            result.x3 = x3;
+
+            if (1.0 - FP64Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0 + FP64Utility.TWO_EPSYLON)
+            {
+                return;
+            }
+
+            result.Normalize(squareModule);
+        }
+
+        public static void LoadIdle(out FP64Versor versor)
+        {
+            versor.s0 = 1.0;
+            versor.x1 = 0.0;
+            versor.x2 = 0.0;
+            versor.x3 = 0.0;
+        }
+        
+        public static void LoadValues(double s0, double x1, double x2, double x3, out FP64Versor versor)
+        {
+            versor.s0 = s0;
+            versor.x1 = x1;
+            versor.x2 = x2;
+            versor.x3 = x3;
+
+            double squareModule = (s0 * s0 + x1 * x1) + (x2 * x2 + x3 * x3);
+            
+            if (1.0 - FP64Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0 + FP64Utility.TWO_EPSYLON)
+            {
+                return;
+            }
+
+            versor.Normalize(squareModule);
+        }
+    }
+}