Большое переименование и рефакторинг углов / Big renaming and refactoring of angles

BGC/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 BGC
+{
+    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
+    }
+}
+

BGC/BGC.csproj

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

BGC/F32Angle.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 BGC
+{
+    public static class F32Angle
+    {
+        public static float ToRadians(float angle, AngleUnit unit)
+        {
+            if (unit == AngleUnit.DEGREES)
+            {
+                return angle * F32Degrees.RADIANS_IN_DEGREE;
+            }
+
+            if (unit == AngleUnit.TURNS)
+            {
+                return angle * F32Radians.TWO_PI;
+            }
+
+            return angle;
+        }
+
+        public static float ToDegrees(float angle, AngleUnit unit)
+        {
+            if (unit == AngleUnit.RADIANS)
+            {
+                return angle * F32Radians.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 * F32Radians.TURNS_IN_RADIAN;
+            }
+
+            if (unit == AngleUnit.DEGREES)
+            {
+                return angle * F32Degrees.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 F32Radians.TWO_PI;
+        }
+    
+        public static float GetHalfCircle(AngleUnit unit)
+        {
+            if (unit == AngleUnit.DEGREES)
+            {
+                return 180.0f;
+            }
+
+            if (unit == AngleUnit.TURNS)
+            {
+                return 0.5f;
+            }
+
+            return F32Radians.PI;
+        }
+    
+        public static float GetQuarterCircle(AngleUnit unit)
+        {
+            if (unit == AngleUnit.DEGREES)
+            {
+                return 90.0f;
+            }
+
+            if (unit == AngleUnit.TURNS)
+            {
+                return 0.25f;
+            }
+
+            return F32Radians.HALF_OF_PI;
+        }
+
+        public static float Normalize(float angle, AngleUnit unit, AngleRange range)
+        {
+            if (unit == AngleUnit.DEGREES)
+            {
+                return F32Degrees.Normalize(angle, range);
+            }
+
+            if (unit == AngleUnit.TURNS)
+            {
+                return F32Turns.Normalize(angle, range);
+            }
+
+            return F32Radians.Normalize(angle, range);
+        }
+    }
+}

BGC/F32Degrees.cs

@@ -0,0 +1,68 @@
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 18 Nov 2024
+ */
+
+namespace BGC
+{
+    public class F32Degrees
+    {
+        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;
+        }
+    }
+}

BGC/F32Matrix2x2.cs

@@ -0,0 +1,352 @@
+/*
+ * 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 BGC
+{
+    public struct F32Matrix2x2
+    {
+        public float r1c1, r1c2;
+
+        public float r2c1, r2c2;
+
+        public F32Matrix2x2()
+        {
+            this.r1c1 = 0.0f;
+            this.r1c2 = 0.0f;
+
+            this.r2c1 = 0.0f;
+            this.r2c2 = 0.0f;
+        }
+
+        public F32Matrix2x2(float d1, float d2)
+        {
+            this.r1c1 = d1;
+            this.r1c2 = 0.0f;
+
+            this.r2c1 = 0.0f;
+            this.r2c2 = d2;
+        }
+
+        public F32Matrix2x2(in F32Matrix2x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+        }
+
+        public F32Matrix2x2(in F64Matrix2x2 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 -F32Utility.EPSYLON <= determinant && determinant <= F32Utility.EPSYLON;
+        }
+
+        public void Transpose()
+        {
+            (this.r1c2, this.r2c1) = (this.r2c1, this.r1c2);
+        }
+
+        public bool Invert()
+        {
+            float determinant = this.GetDeterminant();
+
+            if (-F32Utility.EPSYLON <= determinant && determinant <= F32Utility.EPSYLON)
+            {
+                return false;
+            }
+
+            float r1c1 =  this.r2c2;
+            float r1c2 = -this.r1c2;
+
+            float r2c1 = -this.r2c1;
+            float r2c2 =  this.r1c1;
+
+            this.r1c1 = r1c1 / determinant;
+            this.r1c2 = r1c2 / determinant;
+
+            this.r2c1 = r2c1 / determinant;
+            this.r2c2 = r2c2 / determinant;
+
+            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 F32Matrix2x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+        }
+
+        public void SetValues(in F64Matrix2x2 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 F32Matrix2x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r2c2 = matrix.r2c2;
+            (this.r1c2, this.r2c1) = (matrix.r2c1, matrix.r1c2);
+        }
+
+        public void SetTransposedOf(in F64Matrix2x2 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 F32Matrix2x2 matrix)
+        {
+            float determinant = matrix.GetDeterminant();
+
+            if (-F32Utility.EPSYLON <= determinant && determinant <= F32Utility.EPSYLON)
+            {
+                return false;
+            }
+
+            float r1c1 =  matrix.r2c2;
+            float r1c2 = -matrix.r1c2;
+
+            float r2c1 = -matrix.r2c1;
+            float r2c2 =  matrix.r1c1;
+
+            this.r1c1 = r1c1 / determinant;
+            this.r1c2 = r1c2 / determinant;
+
+            this.r2c1 = r2c1 / determinant;
+            this.r2c2 = r2c2 / determinant;
+
+            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 static void Add(in F32Matrix2x2 matrix1, in F32Matrix2x2 matrix2, out F32Matrix2x2 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 F32Matrix2x2 minuend, in F32Matrix2x2 subtrahend, out F32Matrix2x2 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 F32Matrix2x2 multiplicand, float multiplier, out F32Matrix2x2 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 F32Matrix2x2 dividend, float divisor, out F32Matrix2x2 quotient)
+        {
+            quotient.r1c1 = dividend.r1c1 / divisor;
+            quotient.r1c2 = dividend.r1c2 / divisor;
+
+            quotient.r2c1 = dividend.r2c1 / divisor;
+            quotient.r2c2 = dividend.r2c2 / divisor;
+        }
+
+        public static void GetWeightedSum2(
+            float weight1, in F32Matrix2x2 matrix1,
+            float weight2, in F32Matrix2x2 matrix2,
+            out F32Matrix2x2 sum)
+        {
+            sum.r1c1 = matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2;
+            sum.r1c2 = matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2;
+
+            sum.r2c1 = matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2;
+            sum.r2c2 = matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2;
+        }
+
+        public static void GetWeightedSum3(
+            float weight1, in F32Matrix2x2 matrix1,
+            float weight2, in F32Matrix2x2 matrix2,
+            float weight3, in F32Matrix2x2 matrix3,
+            out F32Matrix2x2 sum)
+        {
+            sum.r1c1 = matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2 + matrix3.r1c1 * weight3;
+            sum.r1c2 = matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2 + matrix3.r1c2 * weight3;
+
+            sum.r2c1 = matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2 + matrix3.r2c1 * weight3;
+            sum.r2c2 = matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2 + matrix3.r2c2 * weight3;
+        }
+
+        public static void GetWeightedSum4(
+            float weight1, in F32Matrix2x2 matrix1,
+            float weight2, in F32Matrix2x2 matrix2,
+            float weight3, in F32Matrix2x2 matrix3,
+            float weight4, in F32Matrix2x2 matrix4,
+            out F32Matrix2x2 sum)
+        {
+            sum.r1c1 = (matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2) + (matrix3.r1c1 * weight3 + matrix4.r1c1 * weight4);
+            sum.r1c2 = (matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2) + (matrix3.r1c2 * weight3 + matrix4.r1c2 * weight4);
+
+            sum.r2c1 = (matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2) + (matrix3.r2c1 * weight3 + matrix4.r2c1 * weight4);
+            sum.r2c2 = (matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2) + (matrix3.r2c2 * weight3 + matrix4.r2c2 * weight4);
+        }
+
+        public static void GetWeightedSum5(
+            float weight1, in F32Matrix2x2 matrix1,
+            float weight2, in F32Matrix2x2 matrix2,
+            float weight3, in F32Matrix2x2 matrix3,
+            float weight4, in F32Matrix2x2 matrix4,
+            float weight5, in F32Matrix2x2 matrix5,
+            out F32Matrix2x2 sum)
+        {
+            sum.r1c1 = (matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2) + (matrix3.r1c1 * weight3 + matrix4.r1c1 * weight4) + matrix5.r1c1 * weight5;
+            sum.r1c2 = (matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2) + (matrix3.r1c2 * weight3 + matrix4.r1c2 * weight4) + matrix5.r1c2 * weight5;
+
+            sum.r2c1 = (matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2) + (matrix3.r2c1 * weight3 + matrix4.r2c1 * weight4) + matrix5.r2c1 * weight5;
+            sum.r2c2 = (matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2) + (matrix3.r2c2 * weight3 + matrix4.r2c2 * weight4) + matrix5.r2c2 * weight5;
+        }
+
+        public static void GetRightProduct(in F32Matrix2x2 matrix, in F32Vector2 vector, out F32Vector2 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 F32Vector2 vector, in F32Matrix2x2 matrix, out F32Vector2 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 F32Matrix2x2 matrix)
+        {
+            matrix.r1c1 = 0.0f;
+            matrix.r1c2 = 0.0f;
+
+            matrix.r2c1 = 0.0f;
+            matrix.r2c2 = 0.0f;
+        }
+
+        public static void LoadIdentity(out F32Matrix2x2 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 F32Matrix2x2 matrix)
+        {
+            matrix.r1c1 = d1;
+            matrix.r1c2 = 0.0f;
+
+            matrix.r2c1 = 0.0f;
+            matrix.r2c2 = d2;
+        }
+    }
+}

BGC/F32Matrix2x3.cs

@@ -0,0 +1,310 @@
+/*
+ * 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 BGC
+{
+    public struct F32Matrix2x3
+    {
+        public float r1c1, r1c2;
+        public float r2c1, r2c2;
+        public float r3c1, r3c2;
+    
+        public F32Matrix2x3()
+        {
+            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 F32Matrix2x3(in F32Matrix2x3 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 F32Matrix2x3(in F64Matrix2x3 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 F32Matrix2x3(in F32Matrix3x2 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 F32Matrix2x3(in F64Matrix3x2 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 F32Matrix2x3 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 F64Matrix2x3 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 F32Matrix3x2 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 F64Matrix3x2 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 static void Add(in F32Matrix2x3 matrix1, in F32Matrix2x3 matrix2, out F32Matrix2x3 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 F32Matrix2x3 minuend, in F32Matrix2x3 subtrahend, out F32Matrix2x3 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 F32Matrix2x3 multiplicand, float multiplier, out F32Matrix2x3 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 F32Matrix2x3 dividend, float divisor, out F32Matrix2x3 quotient)
+        {
+            quotient.r1c1 = dividend.r1c1 / divisor;
+            quotient.r1c2 = dividend.r1c2 / divisor;
+
+            quotient.r2c1 = dividend.r2c1 / divisor;
+            quotient.r2c2 = dividend.r2c2 / divisor;
+
+            quotient.r3c1 = dividend.r3c1 / divisor;
+            quotient.r3c2 = dividend.r3c2 / divisor;
+        }
+
+        public static void GetWeightedSum2(
+            float weight1, in F32Matrix2x3 matrix1,
+            float weight2, in F32Matrix2x3 matrix2,
+            out F32Matrix2x3 sum)
+        {
+            sum.r1c1 = matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2;
+            sum.r1c2 = matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2;
+
+            sum.r2c1 = matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2;
+            sum.r2c2 = matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2;
+
+            sum.r3c1 = matrix1.r3c1 * weight1 + matrix2.r3c1 * weight2;
+            sum.r3c2 = matrix1.r3c2 * weight1 + matrix2.r3c2 * weight2;
+        }
+
+        public static void GetWeightedSum3(
+            float weight1, in F32Matrix2x3 matrix1,
+            float weight2, in F32Matrix2x3 matrix2,
+            float weight3, in F32Matrix2x3 matrix3,
+            out F32Matrix2x3 sum)
+        {
+            sum.r1c1 = matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2 + matrix3.r1c1 * weight3;
+            sum.r1c2 = matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2 + matrix3.r1c2 * weight3;
+
+            sum.r2c1 = matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2 + matrix3.r2c1 * weight3;
+            sum.r2c2 = matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2 + matrix3.r2c2 * weight3;
+
+            sum.r3c1 = matrix1.r3c1 * weight1 + matrix2.r3c1 * weight2 + matrix3.r3c1 * weight3;
+            sum.r3c2 = matrix1.r3c2 * weight1 + matrix2.r3c2 * weight2 + matrix3.r3c2 * weight3;
+        }
+
+        public static void GetWeightedSum4(
+            float weight1, in F32Matrix2x3 matrix1,
+            float weight2, in F32Matrix2x3 matrix2,
+            float weight3, in F32Matrix2x3 matrix3,
+            float weight4, in F32Matrix2x3 matrix4,
+            out F32Matrix2x3 sum)
+        {
+            sum.r1c1 = (matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2) + (matrix3.r1c1 * weight3 + matrix4.r1c1 * weight4);
+            sum.r1c2 = (matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2) + (matrix3.r1c2 * weight3 + matrix4.r1c2 * weight4);
+
+            sum.r2c1 = (matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2) + (matrix3.r2c1 * weight3 + matrix4.r2c1 * weight4);
+            sum.r2c2 = (matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2) + (matrix3.r2c2 * weight3 + matrix4.r2c2 * weight4);
+
+            sum.r3c1 = (matrix1.r3c1 * weight1 + matrix2.r3c1 * weight2) + (matrix3.r3c1 * weight3 + matrix4.r3c1 * weight4);
+            sum.r3c2 = (matrix1.r3c2 * weight1 + matrix2.r3c2 * weight2) + (matrix3.r3c2 * weight3 + matrix4.r3c2 * weight4);
+        }
+
+        public static void GetWeightedSum5(
+            float weight1, in F32Matrix2x3 matrix1,
+            float weight2, in F32Matrix2x3 matrix2,
+            float weight3, in F32Matrix2x3 matrix3,
+            float weight4, in F32Matrix2x3 matrix4,
+            float weight5, in F32Matrix2x3 matrix5,
+            out F32Matrix2x3 sum)
+        {
+            sum.r1c1 = (matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2) + (matrix3.r1c1 * weight3 + matrix4.r1c1 * weight4) + matrix5.r1c1 * weight5;
+            sum.r1c2 = (matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2) + (matrix3.r1c2 * weight3 + matrix4.r1c2 * weight4) + matrix5.r1c2 * weight5;
+
+            sum.r2c1 = (matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2) + (matrix3.r2c1 * weight3 + matrix4.r2c1 * weight4) + matrix5.r2c1 * weight5;
+            sum.r2c2 = (matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2) + (matrix3.r2c2 * weight3 + matrix4.r2c2 * weight4) + matrix5.r2c2 * weight5;
+
+            sum.r3c1 = (matrix1.r3c1 * weight1 + matrix2.r3c1 * weight2) + (matrix3.r3c1 * weight3 + matrix4.r3c1 * weight4) + matrix5.r3c1 * weight5;
+            sum.r3c2 = (matrix1.r3c2 * weight1 + matrix2.r3c2 * weight2) + (matrix3.r3c2 * weight3 + matrix4.r3c2 * weight4) + matrix5.r3c2 * weight5;
+        }
+
+        public static void GetRightProduct(in F32Matrix2x3 matrix, in F32Vector2 vector, out F32Vector3 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 F32Vector3 vector, in F32Matrix2x3 matrix, out F32Vector2 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;
+        }
+    }
+}

BGC/F32Matrix3x2.cs

@@ -0,0 +1,289 @@
+/*
+ * 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 BGC
+{
+    public struct F32Matrix3x2
+    {
+        public float r1c1, r1c2, r1c3;
+        public float r2c1, r2c2, r2c3;
+
+        public F32Matrix3x2()
+        {
+            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 F32Matrix3x2(in F32Matrix3x2 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 F32Matrix3x2(in F64Matrix3x2 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 F32Matrix3x2(in F32Matrix2x3 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 F32Matrix3x2(in F64Matrix2x3 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 F32Matrix3x2 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 F64Matrix3x2 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 F32Matrix2x3 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 F64Matrix2x3 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 static void Add(in F32Matrix3x2 matrix1, in F32Matrix3x2 matrix2, out F32Matrix3x2 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 F32Matrix3x2 minuend, in F32Matrix3x2 subtrahend, out F32Matrix3x2 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 F32Matrix3x2 multiplicand, float multiplier, out F32Matrix3x2 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 F32Matrix3x2 dividend, float divisor, out F32Matrix3x2 quotient)
+        {
+            quotient.r1c1 = dividend.r1c1 / divisor;
+            quotient.r1c2 = dividend.r1c2 / divisor;
+            quotient.r1c3 = dividend.r1c3 / divisor;
+
+            quotient.r2c1 = dividend.r2c1 / divisor;
+            quotient.r2c2 = dividend.r2c2 / divisor;
+            quotient.r2c3 = dividend.r2c3 / divisor;
+        }
+
+        public static void GetWeightedSum2(
+            float weight1, in F32Matrix3x2 matrix1,
+            float weight2, in F32Matrix3x2 matrix2,
+            out F32Matrix3x2 sum)
+        {
+            sum.r1c1 = matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2;
+            sum.r1c2 = matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2;
+            sum.r1c3 = matrix1.r1c3 * weight1 + matrix2.r1c3 * weight2;
+
+            sum.r2c1 = matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2;
+            sum.r2c2 = matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2;
+            sum.r2c3 = matrix1.r2c3 * weight1 + matrix2.r2c3 * weight2;
+        }
+
+        public static void GetWeightedSum3(
+            float weight1, in F32Matrix3x2 matrix1,
+            float weight2, in F32Matrix3x2 matrix2,
+            float weight3, in F32Matrix3x2 matrix3,
+            out F32Matrix3x2 sum)
+        {
+            sum.r1c1 = matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2 + matrix3.r1c1 * weight3;
+            sum.r1c2 = matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2 + matrix3.r1c2 * weight3;
+            sum.r1c3 = matrix1.r1c3 * weight1 + matrix2.r1c3 * weight2 + matrix3.r1c3 * weight3;
+
+            sum.r2c1 = matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2 + matrix3.r2c1 * weight3;
+            sum.r2c2 = matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2 + matrix3.r2c2 * weight3;
+            sum.r2c3 = matrix1.r2c3 * weight1 + matrix2.r2c3 * weight2 + matrix3.r2c3 * weight3;
+        }
+
+        public static void GetWeightedSum4(
+            float weight1, in F32Matrix3x2 matrix1,
+            float weight2, in F32Matrix3x2 matrix2,
+            float weight3, in F32Matrix3x2 matrix3,
+            float weight4, in F32Matrix3x2 matrix4,
+            out F32Matrix3x2 sum)
+        {
+            sum.r1c1 = (matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2) + (matrix3.r1c1 * weight3 + matrix4.r1c1 * weight4);
+            sum.r1c2 = (matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2) + (matrix3.r1c2 * weight3 + matrix4.r1c2 * weight4);
+            sum.r1c3 = (matrix1.r1c3 * weight1 + matrix2.r1c3 * weight2) + (matrix3.r1c3 * weight3 + matrix4.r1c3 * weight4);
+
+            sum.r2c1 = (matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2) + (matrix3.r2c1 * weight3 + matrix4.r2c1 * weight4);
+            sum.r2c2 = (matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2) + (matrix3.r2c2 * weight3 + matrix4.r2c2 * weight4);
+            sum.r2c3 = (matrix1.r2c3 * weight1 + matrix2.r2c3 * weight2) + (matrix3.r2c3 * weight3 + matrix4.r2c3 * weight4);
+        }
+
+        public static void GetWeightedSum5(
+            float weight1, in F32Matrix3x2 matrix1,
+            float weight2, in F32Matrix3x2 matrix2,
+            float weight3, in F32Matrix3x2 matrix3,
+            float weight4, in F32Matrix3x2 matrix4,
+            float weight5, in F32Matrix3x2 matrix5,
+            out F32Matrix3x2 sum)
+        {
+            sum.r1c1 = (matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2) + (matrix3.r1c1 * weight3 + matrix4.r1c1 * weight4) + matrix5.r1c1 * weight5;
+            sum.r1c2 = (matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2) + (matrix3.r1c2 * weight3 + matrix4.r1c2 * weight4) + matrix5.r1c2 * weight5;
+            sum.r1c3 = (matrix1.r1c3 * weight1 + matrix2.r1c3 * weight2) + (matrix3.r1c3 * weight3 + matrix4.r1c3 * weight4) + matrix5.r1c3 * weight5;
+
+            sum.r2c1 = (matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2) + (matrix3.r2c1 * weight3 + matrix4.r2c1 * weight4) + matrix5.r2c1 * weight5;
+            sum.r2c2 = (matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2) + (matrix3.r2c2 * weight3 + matrix4.r2c2 * weight4) + matrix5.r2c2 * weight5;
+            sum.r2c3 = (matrix1.r2c3 * weight1 + matrix2.r2c3 * weight2) + (matrix3.r2c3 * weight3 + matrix4.r2c3 * weight4) + matrix5.r2c3 * weight5;
+        }
+
+        public static void GetRightProduct(in F32Matrix3x2 matrix, in F32Vector3 vector, out F32Vector2 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 F32Vector2 vector, in F32Matrix3x2 matrix, out F32Vector3 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;
+        }
+    }
+}

BGC/F32Matrix3x3.cs

@@ -0,0 +1,532 @@
+/*
+ * 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 BGC
+{
+    public struct F32Matrix3x3
+    {
+        public float r1c1, r1c2, r1c3;
+
+        public float r2c1, r2c2, r2c3;
+
+        public float r3c1, r3c2, r3c3;
+
+        public F32Matrix3x3()
+        {
+            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 F32Matrix3x3(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.r3c1 = 0.0f;
+            this.r3c2 = 0.0f;
+            this.r3c3 = d3;
+        }
+
+        public F32Matrix3x3(in F32Matrix3x3 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 F32Matrix3x3(in F64Matrix3x3 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 -F32Utility.EPSYLON <= determinant && determinant <= F32Utility.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 (-F32Utility.EPSYLON <= determinant && determinant <= F32Utility.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;
+
+            this.r1c1 = r1c1 / determinant;
+            this.r1c2 = r1c2 / determinant;
+            this.r1c3 = r1c3 / determinant;
+
+            this.r2c1 = r2c1 / determinant;
+            this.r2c2 = r2c2 / determinant;
+            this.r2c3 = r2c3 / determinant;
+
+            this.r3c1 = r3c1 / determinant;
+            this.r3c2 = r3c2 / determinant;
+            this.r3c3 = r3c3 / determinant;
+
+            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 F32Matrix3x3 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 F64Matrix3x3 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 F32Matrix3x3 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 F64Matrix3x3 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 F32Matrix3x3 matrix)
+        {
+            float determinant = matrix.GetDeterminant();
+    
+            if (-F32Utility.EPSYLON <= determinant && determinant <= F32Utility.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;
+
+            this.r1c1 = r1c1 / determinant;
+            this.r1c2 = r1c2 / determinant;
+            this.r1c3 = r1c3 / determinant;
+
+            this.r2c1 = r2c1 / determinant;
+            this.r2c2 = r2c2 / determinant;
+            this.r2c3 = r2c3 / determinant;
+
+            this.r3c1 = r3c1 / determinant;
+            this.r3c2 = r3c2 / determinant;
+            this.r3c3 = r3c3 / determinant;
+
+            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 static void Add(in F32Matrix3x3 matrix1, in F32Matrix3x3 matrix2, out F32Matrix3x3 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 F32Matrix3x3 minuend, in F32Matrix3x3 subtrahend, out F32Matrix3x3 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 F32Matrix3x3 multiplicand, float multiplier, out F32Matrix3x3 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 F32Matrix3x3 dividend, float divisor, out F32Matrix3x3 quotient)
+        {
+            quotient.r1c1 = dividend.r1c1 / divisor;
+            quotient.r1c2 = dividend.r1c2 / divisor;
+            quotient.r1c3 = dividend.r1c3 / divisor;
+
+            quotient.r2c1 = dividend.r2c1 / divisor;
+            quotient.r2c2 = dividend.r2c2 / divisor;
+            quotient.r2c3 = dividend.r2c3 / divisor;
+
+            quotient.r3c1 = dividend.r3c1 / divisor;
+            quotient.r3c2 = dividend.r3c2 / divisor;
+            quotient.r3c3 = dividend.r3c3 / divisor;
+        }
+
+        public static void GetWeightedSum2(
+            float weight1, in F32Matrix3x3 matrix1,
+            float weight2, in F32Matrix3x3 matrix2,
+            out F32Matrix3x3 sum)
+        {
+            sum.r1c1 = matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2;
+            sum.r1c2 = matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2;
+            sum.r1c3 = matrix1.r1c3 * weight1 + matrix2.r1c3 * weight2;
+
+            sum.r2c1 = matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2;
+            sum.r2c2 = matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2;
+            sum.r2c3 = matrix1.r2c3 * weight1 + matrix2.r2c3 * weight2;
+
+            sum.r3c1 = matrix1.r3c1 * weight1 + matrix2.r3c1 * weight2;
+            sum.r3c2 = matrix1.r3c2 * weight1 + matrix2.r3c2 * weight2;
+            sum.r3c3 = matrix1.r3c3 * weight1 + matrix2.r3c3 * weight2;
+        }
+
+        public static void GetWeightedSum3(
+            float weight1, in F32Matrix3x3 matrix1,
+            float weight2, in F32Matrix3x3 matrix2,
+            float weight3, in F32Matrix3x3 matrix3,
+            out F32Matrix3x3 sum)
+        {
+            sum.r1c1 = matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2 + matrix3.r1c1 * weight3;
+            sum.r1c2 = matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2 + matrix3.r1c2 * weight3;
+            sum.r1c3 = matrix1.r1c3 * weight1 + matrix2.r1c3 * weight2 + matrix3.r1c3 * weight3;
+
+            sum.r2c1 = matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2 + matrix3.r2c1 * weight3;
+            sum.r2c2 = matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2 + matrix3.r2c2 * weight3;
+            sum.r2c3 = matrix1.r2c3 * weight1 + matrix2.r2c3 * weight2 + matrix3.r2c3 * weight3;
+
+            sum.r3c1 = matrix1.r3c1 * weight1 + matrix2.r3c1 * weight2 + matrix3.r3c1 * weight3;
+            sum.r3c2 = matrix1.r3c2 * weight1 + matrix2.r3c2 * weight2 + matrix3.r3c2 * weight3;
+            sum.r3c3 = matrix1.r3c3 * weight1 + matrix2.r3c3 * weight2 + matrix3.r3c3 * weight3;
+        }
+
+        public static void GetWeightedSum4(
+            float weight1, in F32Matrix3x3 matrix1,
+            float weight2, in F32Matrix3x3 matrix2,
+            float weight3, in F32Matrix3x3 matrix3,
+            float weight4, in F32Matrix3x3 matrix4,
+            out F32Matrix3x3 sum)
+        {
+            sum.r1c1 = (matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2) + (matrix3.r1c1 * weight3 + matrix4.r1c1 * weight4);
+            sum.r1c2 = (matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2) + (matrix3.r1c2 * weight3 + matrix4.r1c2 * weight4);
+            sum.r1c3 = (matrix1.r1c3 * weight1 + matrix2.r1c3 * weight2) + (matrix3.r1c3 * weight3 + matrix4.r1c3 * weight4);
+
+            sum.r2c1 = (matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2) + (matrix3.r2c1 * weight3 + matrix4.r2c1 * weight4);
+            sum.r2c2 = (matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2) + (matrix3.r2c2 * weight3 + matrix4.r2c2 * weight4);
+            sum.r2c3 = (matrix1.r2c3 * weight1 + matrix2.r2c3 * weight2) + (matrix3.r2c3 * weight3 + matrix4.r2c3 * weight4);
+
+            sum.r3c1 = (matrix1.r3c1 * weight1 + matrix2.r3c1 * weight2) + (matrix3.r3c1 * weight3 + matrix4.r3c1 * weight4);
+            sum.r3c2 = (matrix1.r3c2 * weight1 + matrix2.r3c2 * weight2) + (matrix3.r3c2 * weight3 + matrix4.r3c2 * weight4);
+            sum.r3c3 = (matrix1.r3c3 * weight1 + matrix2.r3c3 * weight2) + (matrix3.r3c3 * weight3 + matrix4.r3c3 * weight4);
+        }
+
+        public static void GetWeightedSum5(
+            float weight1, in F32Matrix3x3 matrix1,
+            float weight2, in F32Matrix3x3 matrix2,
+            float weight3, in F32Matrix3x3 matrix3,
+            float weight4, in F32Matrix3x3 matrix4,
+            float weight5, in F32Matrix3x3 matrix5,
+            out F32Matrix3x3 sum)
+        {
+            sum.r1c1 = (matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2) + (matrix3.r1c1 * weight3 + matrix4.r1c1 * weight4) + matrix5.r1c1 * weight5;
+            sum.r1c2 = (matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2) + (matrix3.r1c2 * weight3 + matrix4.r1c2 * weight4) + matrix5.r1c2 * weight5;
+            sum.r1c3 = (matrix1.r1c3 * weight1 + matrix2.r1c3 * weight2) + (matrix3.r1c3 * weight3 + matrix4.r1c3 * weight4) + matrix5.r1c3 * weight5;
+
+            sum.r2c1 = (matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2) + (matrix3.r2c1 * weight3 + matrix4.r2c1 * weight4) + matrix5.r2c1 * weight5;
+            sum.r2c2 = (matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2) + (matrix3.r2c2 * weight3 + matrix4.r2c2 * weight4) + matrix5.r2c2 * weight5;
+            sum.r2c3 = (matrix1.r2c3 * weight1 + matrix2.r2c3 * weight2) + (matrix3.r2c3 * weight3 + matrix4.r2c3 * weight4) + matrix5.r2c3 * weight5;
+
+            sum.r3c1 = (matrix1.r3c1 * weight1 + matrix2.r3c1 * weight2) + (matrix3.r3c1 * weight3 + matrix4.r3c1 * weight4) + matrix5.r3c1 * weight5;
+            sum.r3c2 = (matrix1.r3c2 * weight1 + matrix2.r3c2 * weight2) + (matrix3.r3c2 * weight3 + matrix4.r3c2 * weight4) + matrix5.r3c2 * weight5;
+            sum.r3c3 = (matrix1.r3c3 * weight1 + matrix2.r3c3 * weight2) + (matrix3.r3c3 * weight3 + matrix4.r3c3 * weight4) + matrix5.r3c3 * weight5;
+        }
+
+        public static void GetRightProduct(in F32Matrix3x3 matrix, in F32Vector3 vector, out F32Vector3 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 F32Vector3 vector, in F32Matrix3x3 matrix, out F32Vector3 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 F32Matrix3x3 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 F32Matrix3x3 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 F32Matrix3x3 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;
+        }
+    }
+}

BGC/F32MatrixProduct.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 BGC
+{
+    public class F32MatrixProduct
+    {
+        public static void Get2x2At2x2(in F32Matrix2x2 left, in F32Matrix2x2 right, out F32Matrix2x2 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 F32Matrix2x2 left, in F32Matrix3x2 right, out F32Matrix3x2 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 F32Matrix2x3 left, in F32Matrix2x2 right, out F32Matrix2x3 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 F32Matrix2x3 left, in F32Matrix3x2 right, out F32Matrix3x3 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 F32Matrix3x2 left, in F32Matrix3x3 right, out F32Matrix3x2 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 F32Matrix3x2 left, in F32Matrix2x3 right, out F32Matrix2x2 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 F32Matrix3x3 left, in F32Matrix2x3 right, out F32Matrix2x3 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 F32Matrix3x3 left, in F32Matrix3x3 right, out F32Matrix3x3 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;
+        }
+    }
+}

BGC/F32Quaternion.cs

@@ -0,0 +1,203 @@
+using System;
+
+namespace BGC
+{
+    public struct F32Quaternion
+    {
+        public float s0, x1, x2, x3;
+
+        public F32Quaternion(float s0, float x1, float x2, float x3)
+        {
+            this.s0 = s0;
+            this.x1 = x1;
+            this.x2 = x2;
+            this.x3 = x3;
+        }
+
+        public F32Quaternion(in F32Quaternion quaternion)
+        {
+            this.s0 = quaternion.s0;
+            this.x1 = quaternion.x1;
+            this.x2 = quaternion.x2;
+            this.x3 = quaternion.x3;
+        }
+
+        public F32Quaternion(in F64Quaternion 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 F32Quaternion quaternion)
+        {
+            this.s0 = quaternion.s0;
+            this.x1 = quaternion.x1;
+            this.x2 = quaternion.x2;
+            this.x3 = quaternion.x3;
+        }
+
+        public void SetValues(in F64Quaternion 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 F32Quaternion quaternion)
+        {
+            this.s0 = quaternion.s0;
+            this.x1 = -quaternion.x1;
+            this.x2 = -quaternion.x2;
+            this.x3 = -quaternion.x3;
+        }
+
+        public readonly void MakeRotationMatrix(out F32Matrix3x3 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 (-F32Utility.EPSYLON <= squareModule && squareModule <= F32Utility.EPSYLON)
+            {
+                F32Matrix3x3.LoadIdentity(out matrix);
+                return;
+            }
+
+            float corrector1;
+            float corrector2;
+
+            if (1.0f - F32Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0f + F32Utility.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 F32Matrix3x3 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 (-F32Utility.EPSYLON <= squareModule && squareModule <= F32Utility.EPSYLON)
+            {
+                F32Matrix3x3.LoadIdentity(out matrix);
+                return;
+            }
+
+            float corrector1;
+            float corrector2;
+
+            if (1.0f - F32Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0f + F32Utility.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 F32Quaternion quaternion1, in F32Quaternion quaternion2, out F32Quaternion 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 F32Quaternion minuend, in F32Quaternion subtrahend, out F32Quaternion 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 F32Quaternion left, in F32Quaternion right, out F32Quaternion 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;
+        }
+    }
+}

BGC/F32Radians.cs

@@ -0,0 +1,76 @@
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 18 Nov 2024
+ */
+
+namespace BGC
+{
+    public class F32Radians
+    {
+        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;
+        }
+    }
+}

BGC/F32Rotation3.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 BGC
+{
+    public struct F32Rotation3
+    {
+        private float angle;
+
+        private F32Vector3 axis;
+
+        public F32Rotation3(F32Rotation3 rotation)
+        {
+            this.angle = rotation.angle;
+            this.axis = rotation.axis;
+        }
+
+        public void Reset()
+        {
+            this.angle = 0.0f;
+            this.axis.Reset();
+        }
+    }
+}

BGC/F32Turns.cs

@@ -0,0 +1,63 @@
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 18 Nov 2024
+ */
+
+namespace BGC
+{
+    public class F32Turns
+    {
+        public static float ToRadians(float turns)
+        {
+            return turns * F32Radians.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 * F32Radians.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;
+        }
+    }
+}

BGC/F32Utility.cs

@@ -0,0 +1,21 @@
+using System;
+
+namespace BGC
+{
+    public class F32Utility
+    {
+        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;
+    }
+}
+

BGC/F32Vector2.cs

@@ -0,0 +1,353 @@
+/*
+ * 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 BGC
+{
+    public struct F32Vector2
+    {
+        public static readonly F32Vector2 ZERO = new F32Vector2(0.0f, 0.0f);
+
+        public float x1;
+        public float x2;
+
+        public F32Vector2(float x1, float x2)
+        {
+            this.x1 = x1;
+            this.x2 = x2;
+        }
+
+        public F32Vector2(in F32Vector2 vector)
+        {
+            this.x1 = vector.x1;
+            this.x2 = vector.x2;
+        }
+
+        public F32Vector2(in F64Vector2 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 - F32Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0f + F32Utility.TWO_EPSYLON)
+            {
+                return 1;
+            }
+
+            if (squareModule <= F32Utility.SQUARE_EPSYLON)
+            {
+                this.Reset();
+                return 0;
+            }
+
+            float module = MathF.Sqrt(squareModule);
+
+            this.x1 /= module;
+            this.x2 /= module;
+
+            return 1;
+        }
+
+        public readonly F32Vector2 GetNormalized()
+        {
+            float squareModule = this.GetSquareModule();
+
+            if (1.0f - F32Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0f + F32Utility.TWO_EPSYLON)
+            {
+                return this;
+            }
+
+            if (squareModule <= F32Utility.SQUARE_EPSYLON)
+            {
+                return ZERO;
+            }
+
+            float module = MathF.Sqrt(squareModule);
+
+            return new F32Vector2(
+                this.x1 / module,
+                this.x2 / module
+            );
+        }
+
+        public void Reverse()
+        {
+            this.x1 = -this.x1;
+            this.x2 = -this.x2;
+        }
+
+        public readonly bool IsZero()
+        {
+            return this.GetSquareModule() <= F32Utility.SQUARE_EPSYLON;
+        }
+
+        public readonly bool IsUnit()
+        {
+            float squareModule = this.GetSquareModule();
+            return 1.0f - F32Utility.TWO_EPSYLON <= squareModule && squareModule <= F32Utility.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 F32Vector2 vector)
+        {
+            this.x1 = vector.x1;
+            this.x2 = vector.x2;
+        }
+
+        public void SetValues(in F64Vector2 vector)
+        {
+            this.x1 = (float)vector.x1;
+            this.x2 = (float)vector.x2;
+        }
+
+        public void SetReverseOf(in F32Vector2 vector)
+        {
+            this.x1 = -vector.x1;
+            this.x2 = -vector.x2;
+        }
+
+        public void SetReverseOf(in F64Vector2 vector)
+        {
+            this.x1 = -(float)vector.x1;
+            this.x2 = -(float)vector.x2;
+        }
+        public readonly override string ToString()
+        {
+            return String.Format("SPVector2({0}, {1})", this.x1, this.x2);
+        }
+
+        public static void Add(in F32Vector2 vector1, in F32Vector2 vector2, out F32Vector2 sum)
+        {
+            sum.x1 = vector1.x1 + vector2.x1;
+            sum.x2 = vector1.x2 + vector2.x2;
+        }
+
+        public static void Subtract(in F32Vector2 minuend, in F32Vector2 subtrahend, out F32Vector2 difference)
+        {
+            difference.x1 = minuend.x1 - subtrahend.x1;
+            difference.x2 = minuend.x2 - subtrahend.x2;
+        }
+
+        public static void Muliply(in F32Vector2 multiplicand, float multiplier, out F32Vector2 product)
+        {
+            product.x1 = multiplicand.x1 * multiplier;
+            product.x2 = multiplicand.x2 * multiplier;
+        }
+
+        public static void Divide(in F32Vector2 dividend, float divisor, out F32Vector2 quotient)
+        {
+            quotient.x1 = dividend.x1 / divisor;
+            quotient.x2 = dividend.x2 / divisor;
+        }
+
+        public static void GetWeightedSum2(
+            float weight1, in F32Vector2 vector1,
+            float weight2, in F32Vector2 vector2,
+            out F32Vector2 sum
+            )
+        {
+            sum.x1 = vector1.x1 * weight1 + vector2.x1 * weight2;
+            sum.x2 = vector1.x2 * weight1 + vector2.x2 * weight2;
+        }
+
+        public static void GetWeightedSum3(
+            float weight1, in F32Vector2 vector1,
+            float weight2, in F32Vector2 vector2,
+            float weight3, in F32Vector2 vector3,
+            out F32Vector2 sum
+            )
+        {
+            sum.x1 = vector1.x1 * weight1 + vector2.x1 * weight2 + vector3.x1 * weight3;
+            sum.x2 = vector1.x2 * weight1 + vector2.x2 * weight2 + vector3.x2 * weight3;
+        }
+
+        public static void GetWeightedSum4(
+            float weight1, in F32Vector2 vector1,
+            float weight2, in F32Vector2 vector2,
+            float weight3, in F32Vector2 vector3,
+            float weight4, in F32Vector2 vector4,
+            out F32Vector2 sum
+            )
+        {
+            sum.x1 = (vector1.x1 * weight1 + vector2.x1 * weight2) + (vector3.x1 * weight3 + vector4.x1 * weight4);
+            sum.x2 = (vector1.x2 * weight1 + vector2.x2 * weight2) + (vector3.x2 * weight3 + vector4.x2 * weight4);
+        }
+
+        public static void GetWeightedSum5(
+            float weight1, in F32Vector2 vector1,
+            float weight2, in F32Vector2 vector2,
+            float weight3, in F32Vector2 vector3,
+            float weight4, in F32Vector2 vector4,
+            float weight5, in F32Vector2 vector5,
+            out F32Vector2 sum
+            )
+        {
+            sum.x1 = (vector1.x1 * weight1 + vector2.x1 * weight2) + (vector3.x1 * weight3 + vector4.x1 * weight4);
+            sum.x2 = (vector1.x2 * weight1 + vector2.x2 * weight2) + (vector3.x2 * weight3 + vector4.x2 * weight4);
+        }
+
+        public static void GetMean2(
+            in F32Vector2 vector1,
+            in F32Vector2 vector2,
+            out F32Vector2 result
+            )
+        {
+            result.x1 = (vector1.x1 + vector2.x1) * 0.5f;
+            result.x2 = (vector1.x2 + vector2.x2) * 0.5f;
+        }
+
+        public static void GetMean3(
+            in F32Vector2 vector1,
+            in F32Vector2 vector2,
+            in F32Vector2 vector3,
+            out F32Vector2 result
+            )
+        {
+            result.x1 = (vector1.x1 + vector2.x1 + vector3.x1) * F32Utility.ONE_THIRD;
+            result.x2 = (vector1.x2 + vector2.x2 + vector3.x2) * F32Utility.ONE_THIRD;
+        }
+
+        public static void GetMean4(
+            in F32Vector2 vector1,
+            in F32Vector2 vector2,
+            in F32Vector2 vector3,
+            in F32Vector2 vector4,
+            out F32Vector2 result
+            )
+        {
+            result.x1 = ((vector1.x1 + vector2.x1) + (vector3.x1 + vector4.x1)) * 0.25f;
+            result.x2 = ((vector1.x2 + vector2.x2) + (vector3.x2 + vector4.x2)) * 0.25f;
+        }
+
+        public static void GetMean5(
+            in F32Vector2 vector1,
+            in F32Vector2 vector2,
+            in F32Vector2 vector3,
+            in F32Vector2 vector4,
+            in F32Vector2 vector5,
+            out F32Vector2 result
+            )
+        {
+            result.x1 = ((vector1.x1 + vector2.x1) + (vector3.x1 + vector4.x1) + vector5.x1) * 0.2f;
+            result.x2 = ((vector1.x2 + vector2.x2) + (vector3.x2 + vector4.x2) + vector5.x2) * 0.2f;
+        }
+
+        public static float GetScalarProduct(in F32Vector2 vector1, in F32Vector2 vector2)
+        {
+            return vector1.x1 * vector2.x1 + vector1.x2 * vector2.x2;
+        }
+ 
+        public static float GetCrossProduct(in F32Vector2 vector1, in F32Vector2 vector2)
+        {
+            return vector1.x1 * vector2.x2 - vector1.x2 * vector2.x1;
+        }
+
+        public static float GetAngle(in F32Vector2 vector1, in F32Vector2 vector2, AngleUnit unit)
+        {
+    		float squareModule1 = vector1.GetSquareModule();
+
+    		if (squareModule1 <= F32Utility.SQUARE_EPSYLON)
+            {
+    			return 0.0f;
+    		}
+
+    		float squareModule2 = vector2.GetSquareModule();
+
+    		if (squareModule2 <= F32Utility.SQUARE_EPSYLON)
+            {
+    			return 0.0f;
+    		}
+
+    		float cosine = F32Vector2.GetScalarProduct(vector1, vector2) / MathF.Sqrt(squareModule1 * squareModule2);
+
+    		if (1.0f - F32Utility.EPSYLON <= cosine)
+            {
+                return 0.0f;
+    		}
+
+    		if (cosine <= -(1.0f - F32Utility.EPSYLON))
+            {
+                return F32Angle.GetHalfCircle(unit);
+		    }
+    
+    		return F32Radians.ToUnits(MathF.Acos(cosine), unit);
+        }
+    
+        public static float GetSquareDistance(in F32Vector2 vector1, in F32Vector2 vector2)
+        {
+            float dx1 = vector1.x1 - vector2.x1;
+            float dx2 = vector1.x2 - vector2.x2;
+        
+            return dx1 * dx1 + dx2 * dx2;
+        }
+    
+        public static float GetDistance(in F32Vector2 vector1, in F32Vector2 vector2)
+        {
+            return MathF.Sqrt(GetSquareDistance(vector1, vector2));
+        }
+    
+        public static bool AreEqual(in F32Vector2 vector1, in F32Vector2 vector2)
+        {
+            float squareModule1 = vector1.GetSquareModule();
+            float squareModule2 = vector2.GetSquareModule();
+            float squareModule3 = GetSquareDistance(vector1, vector2);
+
+            // 2.0f means dimension amount
+            if (squareModule1 < F32Utility.EPSYLON_EFFECTIVENESS_LIMIT || squareModule2 < F32Utility.EPSYLON_EFFECTIVENESS_LIMIT)
+            {
+                return squareModule3 < (2.0f * F32Utility.SQUARE_EPSYLON);
+            }
+
+            if (squareModule1 <= squareModule2)
+            {
+                return squareModule3 <= (2.0f * F32Utility.SQUARE_EPSYLON) * squareModule2;
+            }
+
+            return squareModule3 <= (2.0f * F32Utility.SQUARE_EPSYLON) * squareModule1;
+        }
+    }
+}

BGC/F32Vector3.cs

@@ -0,0 +1,404 @@
+/*
+ * 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 BGC
+{
+    public struct F32Vector3
+    {
+        public static readonly F32Vector3 ZERO = new F32Vector3(0.0f, 0.0f, 0.0f);
+
+        public float x1;
+        public float x2;
+        public float x3;
+
+        public F32Vector3(float x1, float x2, float x3)
+        {
+            this.x1 = x1;
+            this.x2 = x2;
+            this.x3 = x3;
+        }
+
+        public F32Vector3(in F32Vector3 vector)
+        {
+            this.x1 = vector.x1;
+            this.x2 = vector.x2;
+            this.x3 = vector.x3;
+        }
+
+        public F32Vector3(in F64Vector3 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 - F32Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0f + F32Utility.TWO_EPSYLON)
+            {
+                return 1;
+            }
+
+            if (squareModule <= F32Utility.SQUARE_EPSYLON)
+            {
+                this.Reset();
+                return 0;
+            }
+
+            float module = MathF.Sqrt(squareModule);
+
+            this.x1 /= module;
+            this.x2 /= module;
+            this.x3 /= module;
+
+            return 1;
+        }
+
+        public readonly F32Vector3 GetNormalized()
+        {
+            float squareModule = this.GetSquareModule();
+
+            if (1.0f - F32Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0f + F32Utility.TWO_EPSYLON)
+            {
+                return this;
+            }
+
+            if (squareModule <= F32Utility.SQUARE_EPSYLON)
+            {
+                return ZERO;
+            }
+
+            float module = MathF.Sqrt(squareModule);
+
+            return new F32Vector3(
+                this.x1 / module,
+                this.x2 / module,
+                this.x3 / module
+            );
+        }
+
+        public void Reverse()
+        {
+            this.x1 = -this.x1;
+            this.x2 = -this.x2;
+            this.x3 = -this.x3;
+        }
+
+        public readonly bool IsZero()
+        {
+            return this.GetSquareModule() <= F32Utility.SQUARE_EPSYLON;
+        }
+
+        public readonly bool IsUnit()
+        {
+            float squareModule = this.GetSquareModule();
+            return 1.0f - F32Utility.TWO_EPSYLON <= squareModule && squareModule <= F32Utility.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 F64Vector3 vector)
+        {
+            this.x1 = (float)vector.x1;
+            this.x2 = (float)vector.x2;
+            this.x3 = (float)vector.x3;
+        }
+
+        public void SetValues(in F32Vector3 vector)
+        {
+            this.x1 = vector.x1;
+            this.x2 = vector.x2;
+            this.x3 = vector.x3;
+        }
+
+        public void SetReverseOf(in F32Vector3 vector)
+        {
+            this.x1 = -vector.x1;
+            this.x2 = -vector.x2;
+            this.x3 = -vector.x3;
+        }
+
+        public void SetReverseOf(in F64Vector3 vector)
+        {
+            this.x1 = -(float)vector.x1;
+            this.x2 = -(float)vector.x2;
+            this.x3 = -(float)vector.x3;
+        }
+
+        public readonly override string ToString()
+        {
+            return String.Format("SPVector3({0}, {1}, {2})", this.x1, this.x2, this.x3);
+        }
+
+        public static void Add(in F32Vector3 vector1, in F32Vector3 vector2, out F32Vector3 sum)
+        {
+            sum.x1 = vector1.x1 + vector2.x1;
+            sum.x2 = vector1.x2 + vector2.x2;
+            sum.x3 = vector1.x3 + vector2.x3;
+        }
+
+        public static void Subtract(in F32Vector3 minuend, in F32Vector3 subtrahend, out F32Vector3 difference)
+        {
+            difference.x1 = minuend.x1 - subtrahend.x1;
+            difference.x2 = minuend.x2 - subtrahend.x2;
+            difference.x3 = minuend.x3 - subtrahend.x3;
+        }
+
+        public static void Muliply(in F32Vector3 multiplicand, float multiplier, out F32Vector3 product)
+        {
+            product.x1 = multiplicand.x1 * multiplier;
+            product.x2 = multiplicand.x2 * multiplier;
+            product.x3 = multiplicand.x3 * multiplier;
+        }
+
+        public static void Divide(in F32Vector3 dividend, float divisor, out F32Vector3 quotient)
+        {
+            quotient.x1 = dividend.x1 / divisor;
+            quotient.x2 = dividend.x2 / divisor;
+            quotient.x3 = dividend.x3 / divisor;
+        }
+
+        public static void GetWeightedSum2(
+            float weight1, in F32Vector3 vector1,
+            float weight2, in F32Vector3 vector2,
+            out F32Vector3 sum
+            )
+        {
+            sum.x1 = vector1.x1 * weight1 + vector2.x1 * weight2;
+            sum.x2 = vector1.x2 * weight1 + vector2.x2 * weight2;
+            sum.x3 = vector1.x3 * weight1 + vector2.x3 * weight2;
+        }
+
+        public static void GetWeightedSum3(
+            float weight1, in F32Vector3 vector1,
+            float weight2, in F32Vector3 vector2,
+            float weight3, in F32Vector3 vector3,
+            out F32Vector3 sum
+            )
+        {
+            sum.x1 = vector1.x1 * weight1 + vector2.x1 * weight2 + vector3.x1 * weight3;
+            sum.x2 = vector1.x2 * weight1 + vector2.x2 * weight2 + vector3.x2 * weight3;
+            sum.x3 = vector1.x3 * weight1 + vector2.x3 * weight2 + vector3.x3 * weight3;
+        }
+
+        public static void GetWeightedSum4(
+            float weight1, in F32Vector3 vector1,
+            float weight2, in F32Vector3 vector2,
+            float weight3, in F32Vector3 vector3,
+            float weight4, in F32Vector3 vector4,
+            out F32Vector3 sum
+            )
+        {
+            sum.x1 = (vector1.x1 * weight1 + vector2.x1 * weight2) + (vector3.x1 * weight3 + vector4.x1 * weight4);
+            sum.x2 = (vector1.x2 * weight1 + vector2.x2 * weight2) + (vector3.x2 * weight3 + vector4.x2 * weight4);
+            sum.x3 = (vector1.x3 * weight1 + vector2.x3 * weight2) + (vector3.x3 * weight3 + vector4.x3 * weight4);
+        }
+
+        public static void GetWeightedSum5(
+            float weight1, in F32Vector3 vector1,
+            float weight2, in F32Vector3 vector2,
+            float weight3, in F32Vector3 vector3,
+            float weight4, in F32Vector3 vector4,
+            float weight5, in F32Vector3 vector5,
+            out F32Vector3 sum
+            )
+        {
+            sum.x1 = (vector1.x1 * weight1 + vector2.x1 * weight2) + (vector3.x1 * weight3 + vector4.x1 * weight4) + vector5.x1 * weight5;
+            sum.x2 = (vector1.x2 * weight1 + vector2.x2 * weight2) + (vector3.x2 * weight3 + vector4.x2 * weight4) + vector5.x2 * weight5;
+            sum.x3 = (vector1.x3 * weight1 + vector2.x3 * weight2) + (vector3.x3 * weight3 + vector4.x3 * weight4) + vector5.x3 * weight5;
+        }
+
+        public static void GetMean2(
+            in F32Vector3 vector1,
+            in F32Vector3 vector2,
+            out F32Vector3 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 F32Vector3 vector1,
+            in F32Vector3 vector2,
+            in F32Vector3 vector3,
+            out F32Vector3 result
+            )
+        {
+            result.x1 = (vector1.x1 + vector2.x1 + vector3.x1)  * F32Utility.ONE_THIRD;
+            result.x2 = (vector1.x2 + vector2.x2 + vector3.x2)  * F32Utility.ONE_THIRD;
+            result.x3 = (vector1.x3 + vector2.x3 + vector3.x3)  * F32Utility.ONE_THIRD;
+        }
+
+        public static void GetMean4(
+            in F32Vector3 vector1,
+            in F32Vector3 vector2,
+            in F32Vector3 vector3,
+            in F32Vector3 vector4,
+            out F32Vector3 result
+            )
+        {
+            result.x1 = ((vector1.x1 + vector2.x1) + (vector3.x1 + vector4.x1)) * 0.25f;
+            result.x2 = ((vector1.x2 + vector2.x2) + (vector3.x2 + vector4.x2)) * 0.25f;
+            result.x3 = ((vector1.x3 + vector2.x3) + (vector3.x3 + vector4.x3)) * 0.25f;
+        }
+
+        public static void GetMean5(
+            in F32Vector3 vector1,
+            in F32Vector3 vector2,
+            in F32Vector3 vector3,
+            in F32Vector3 vector4,
+            in F32Vector3 vector5,
+            out F32Vector3 result
+            )
+        {
+            result.x1 = ((vector1.x1 + vector2.x1) + (vector3.x1 + vector4.x1) + vector5.x1) * 0.2f;
+            result.x2 = ((vector1.x2 + vector2.x2) + (vector3.x2 + vector4.x2) + vector5.x2) * 0.2f;
+            result.x3 = ((vector1.x3 + vector2.x3) + (vector3.x3 + vector4.x3) + vector5.x3) * 0.2f;
+        }
+
+        public static float GetScalarProduct(in F32Vector3 vector1, in F32Vector3 vector2)
+        {
+            return vector1.x1 * vector2.x1 + vector1.x2 * vector2.x2 + vector1.x3 * vector2.x3;
+        }
+ 
+        public static void GetCrossProduct(in F32Vector3 vector1, in F32Vector3 vector2, out F32Vector3 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 F32Vector3 vector1, in F32Vector3 vector2, in F32Vector3 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 F32Vector3 vector1, in F32Vector3 vector2, in F32Vector3 vector3, out F32Vector3 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 F32Vector3 vector1, in F32Vector3 vector2, AngleUnit unit)
+        {
+    		float squareModule1 = vector1.GetSquareModule();
+
+    		if (squareModule1 <= F32Utility.SQUARE_EPSYLON)
+            {
+    			return 0.0f;
+    		}
+
+    		float squareModule2 = vector2.GetSquareModule();
+
+    		if (squareModule2 <= F32Utility.SQUARE_EPSYLON)
+            {
+    			return 0.0f;
+    		}
+
+    		float cosine = F32Vector3.GetScalarProduct(vector1, vector2) / MathF.Sqrt(squareModule1 * squareModule2);
+
+    		if (1.0f - F32Utility.EPSYLON <= cosine)
+            {
+                return 0.0f;
+    		}
+
+    		if (cosine <= -(1.0f - F32Utility.EPSYLON))
+            {
+                return F32Angle.GetHalfCircle(unit);
+		    }
+    
+    		return F32Radians.ToUnits(MathF.Acos(cosine), unit);
+        }
+    
+        public static float GetSquareDistance(in F32Vector3 vector1, in F32Vector3 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 F32Vector3 vector1, in F32Vector3 vector2)
+        {
+            return MathF.Sqrt(GetSquareDistance(vector1, vector2));
+        }
+    
+        public static bool AreEqual(in F32Vector3 vector1, in F32Vector3 vector2)
+        {
+            float squareModule1 = vector1.GetSquareModule();
+            float squareModule2 = vector2.GetSquareModule();
+            float squareModule3 = GetSquareDistance(vector1, vector2);
+
+            // 3.0f means dimension amount
+            if (squareModule1 < F32Utility.EPSYLON_EFFECTIVENESS_LIMIT || squareModule2 < F32Utility.EPSYLON_EFFECTIVENESS_LIMIT)
+            {
+                return squareModule3 < (3.0f * F32Utility.SQUARE_EPSYLON);
+            }
+
+            if (squareModule1 <= squareModule2)
+            {
+                return squareModule3 <= (3.0f * F32Utility.SQUARE_EPSYLON) * squareModule2;
+            }
+
+            return squareModule3 <= (3.0f * F32Utility.SQUARE_EPSYLON) * squareModule1;
+        }
+    }
+}

BGC/F32Versor.cs

@@ -0,0 +1,366 @@
+/*
+ * 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 BGC
+{
+    public struct F32Versor
+    {
+        public static readonly F32Versor IDLE = new F32Versor(1.0f, 0.0f, 0.0f, 0.0f, 1.0f);
+
+        private float s0;
+        private float x1;
+        private float x2;
+        private float x3;
+        private float corrector;
+
+        private F32Versor(float s0, float x1, float x2, float x3, float corrector)
+        {
+            this.s0 = s0;
+            this.x1 = x1;
+            this.x2 = x2;
+            this.x3 = x3;
+            this.corrector = corrector;
+        }
+
+        public F32Versor()
+        {
+            this.s0 = 1.0f;
+            this.x1 = 0.0f;
+            this.x2 = 0.0f;
+            this.x3 = 0.0f;
+            this.corrector = 1.0f;
+        }
+
+        public F32Versor(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);
+
+            this.corrector = 2.0f - squareModule;
+
+            if (squareModule < 1.0f - F32Utility.TWO_EPSYLON || 1.0f + F32Utility.TWO_EPSYLON < squareModule)
+            {
+                this.Normalize(squareModule);
+            }
+            else if (s0 < -1.0f + F32Utility.EPSYLON || 1.0f - F32Utility.EPSYLON < s0)
+            {
+                this.Reset();
+            }
+        }
+
+        public F32Versor(in F32Versor versor)
+        {
+            this.s0 = versor.s0;
+            this.x1 = versor.x1;
+            this.x2 = versor.x2;
+            this.x3 = versor.x3;
+            this.corrector = versor.corrector;
+        }
+
+        public F32Versor(in F64Versor versor)
+        {
+            this.s0 = (float)versor.GetScalar();
+            this.x1 = (float)versor.GetX1();
+            this.x2 = (float)versor.GetX2();
+            this.x3 = (float)versor.GetX3();
+
+            float squareModule = (this.s0 * this.s0 + this.x1 * this.x1) + (this.x2 * this.x2 + this.x3 * this.x3);
+
+            this.corrector = 2.0f - squareModule;
+
+            if (squareModule < 1.0f - F32Utility.TWO_EPSYLON || 1.0f + F32Utility.TWO_EPSYLON < squareModule)
+            {
+                this.Normalize(squareModule);
+            }
+            else if (s0 < -1.0f + F32Utility.EPSYLON || 1.0f - F32Utility.EPSYLON < s0)
+            {
+                this.Reset();
+            }
+        }
+
+        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 float GetCorrector()
+        {
+            return this.corrector;
+        }
+
+        public readonly bool IsIdle()
+        {
+            return this.s0 <= -(1.0f - F32Utility.EPSYLON) || (1.0f - F32Utility.EPSYLON) <= this.s0;
+        }
+
+        public void Reset()
+        {
+            this.s0 = 1.0f;
+            this.x1 = 0.0f;
+            this.x2 = 0.0f;
+            this.x3 = 0.0f;
+            this.corrector = 1.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 - F32Utility.TWO_EPSYLON) || 1.0f - F32Utility.TWO_EPSYLON <= this.s0) {
+                return 0.0f;
+            }
+
+            if (-F32Utility.EPSYLON <= this.s0 && this.s0 <= F32Utility.EPSYLON)
+            {
+                return F32Angle.GetHalfCircle(unit);
+            }
+
+            return F32Radians.ToUnits(2.0f * MathF.Acos(s0), unit);
+        }
+
+        public readonly void MakeRotationMatrix(out F32Matrix3x3 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 = 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;
+
+            float corrector2 = 2.0f * this.corrector;
+
+            matrix.r1c1 = this.corrector * ((s0s0 + x1x1) - (x2x2 + x3x3));
+            matrix.r2c2 = this.corrector * ((s0s0 + x2x2) - (x1x1 + x3x3));
+            matrix.r3c3 = this.corrector * ((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 F32Matrix3x3 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 = 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;
+
+            float corrector2 = 2.0f * this.corrector;
+
+            matrix.r1c1 = this.corrector * ((s0s0 + x1x1) - (x2x2 + x3x3));
+            matrix.r2c2 = this.corrector * ((s0s0 + x2x2) - (x1x1 + x3x3));
+            matrix.r3c3 = this.corrector * ((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 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 (squareModule < 1.0f - F32Utility.TWO_EPSYLON || 1.0f + F32Utility.TWO_EPSYLON < squareModule)
+            {
+                this.Normalize(squareModule);
+                return;
+            }
+            
+            if (s0 < -1.0f + F32Utility.EPSYLON || 1.0f - F32Utility.EPSYLON < s0) {
+                this.Reset();
+                return;
+            }
+
+            this.corrector = 2.0f - squareModule;
+        }
+
+        public void SetValues(in F32Versor versor)
+        {
+            this.s0 = versor.s0;
+            this.x1 = versor.x1;
+            this.x2 = versor.x2;
+            this.x3 = versor.x3;
+            this.corrector = versor.corrector;
+        }
+
+        public void SetValues(in F64Versor versor)
+        {
+            this.SetValues(
+                (float)versor.GetScalar(),
+                (float)versor.GetX1(),
+                (float)versor.GetX2(),
+                (float)versor.GetX3()
+            );
+        }
+
+        public void SetInverted(in F32Versor versor)
+        {
+            this.s0 = versor.s0;
+            this.x1 = -versor.x1;
+            this.x2 = -versor.x2;
+            this.x3 = -versor.x3;
+            this.corrector = versor.corrector;
+        }
+
+        public void SetInverted(in F64Versor versor)
+        {
+            this.SetValues(
+                (float)versor.GetScalar(),
+                -(float)versor.GetX1(),
+                -(float)versor.GetX2(),
+                -(float)versor.GetX3()
+            );
+        }
+
+        public readonly void Turn(in F32Vector3 vector, out F32Vector3 result)
+        {
+            float multiplier = 2.0f * this.corrector;
+            float tx1 = multiplier * (this.x2 * vector.x3 - this.x3 * vector.x2);
+            float tx2 = multiplier * (this.x3 * vector.x1 - this.x1 * vector.x3);
+            float tx3 = multiplier * (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 F32Vector3 vector, out F32Vector3 result)
+        {
+            float multiplier = 2.0f * this.corrector;
+            float tx1 = multiplier * (this.x2 * vector.x3 - this.x3 * vector.x2);
+            float tx2 = multiplier * (this.x3 * vector.x1 - this.x1 * vector.x3);
+            float tx3 = multiplier * (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 <= F32Utility.SQUARE_EPSYLON || (this.x1 * this.x1 + this.x2 * this.x2 + this.x3 * this.x3) <= F32Utility.SQUARE_EPSYLON * squareModule)
+            {
+                this.Reset();
+                return;
+            }
+
+            float module = MathF.Sqrt(squareModule);
+
+            this.s0 /= module;
+            this.x1 /= module;
+            this.x2 /= module;
+            this.x3 /= module;
+
+            this.corrector = (this.s0 * this.s0 + this.x1 * this.x1) + (this.x2 * this.x2 + this.x3 * this.x3);
+        }
+
+        public static void Combine(in F32Versor second, in F32Versor first, out F32Versor 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.corrector = 2.0f - squareModule;
+            result.s0 = s0;
+            result.x1 = x1;
+            result.x2 = x2;
+            result.x3 = x3;
+
+            if (squareModule < 1.0f - F32Utility.TWO_EPSYLON || 1.0f + F32Utility.TWO_EPSYLON < squareModule)
+            {
+                result.Normalize(squareModule);
+            }
+        }
+        
+        public static void LoadIdle(out F32Versor versor)
+        {
+            versor.corrector = 1.0f;
+            versor.s0 = 1.0f;
+            versor.x1 = 0.0f;
+            versor.x2 = 0.0f;
+            versor.x3 = 0.0f;
+        }
+    }
+}

BGC/F64Angle.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 BGC
+{
+    public static class F64Angle
+    {
+        public static double ToRadians(double angle, AngleUnit unit)
+        {
+            if (unit == AngleUnit.DEGREES)
+            {
+                return angle * F64Degrees.RADIANS_IN_DEGREE;
+            }
+
+            if (unit == AngleUnit.TURNS)
+            {
+                return angle * F64Radians.TWO_PI;
+            }
+
+            return angle;
+        }
+
+        public static double ToDegrees(double angle, AngleUnit unit)
+        {
+            if (unit == AngleUnit.RADIANS)
+            {
+                return angle * F64Radians.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 * F64Radians.TURNS_IN_RADIAN;
+            }
+
+            if (unit == AngleUnit.DEGREES)
+            {
+                return angle * F64Degrees.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 F64Radians.TWO_PI;
+        }
+    
+        public static double GetHalfCircle(AngleUnit unit)
+        {
+            if (unit == AngleUnit.DEGREES)
+            {
+                return 180.0;
+            }
+
+            if (unit == AngleUnit.TURNS)
+            {
+                return 0.5;
+            }
+
+            return F64Radians.PI;
+        }
+    
+        public static double GetQuarterCircle(AngleUnit unit)
+        {
+            if (unit == AngleUnit.DEGREES)
+            {
+                return 90.0;
+            }
+
+            if (unit == AngleUnit.TURNS)
+            {
+                return 0.25;
+            }
+
+            return F64Radians.HALF_OF_PI;
+        }
+
+        public static double Normalize(double angle, AngleUnit unit, AngleRange range)
+        {
+            if (unit == AngleUnit.DEGREES)
+            {
+                return F64Degrees.Normalize(angle, range);
+            }
+
+            if (unit == AngleUnit.TURNS)
+            {
+                return F64Turns.Normalize(angle, range);
+            }
+
+            return F64Radians.Normalize(angle, range);
+        }
+    }
+}

BGC/F64Degrees.cs

@@ -0,0 +1,68 @@
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 18 Nov 2024
+ */
+
+namespace BGC
+{
+    public class F64Degrees
+    {
+        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;
+        }
+    }
+}

BGC/F64Matrix2x2.cs

@@ -0,0 +1,352 @@
+/*
+ * 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 BGC
+{
+    public struct F64Matrix2x2
+    {
+        public double r1c1, r1c2;
+
+        public double r2c1, r2c2;
+
+        public F64Matrix2x2()
+        {
+            this.r1c1 = 0.0;
+            this.r1c2 = 0.0;
+
+            this.r2c1 = 0.0;
+            this.r2c2 = 0.0;
+        }
+
+        public F64Matrix2x2(double d1, double d2)
+        {
+            this.r1c1 = d1;
+            this.r1c2 = 0.0;
+
+            this.r2c1 = 0.0;
+            this.r2c2 = d2;
+        }
+
+        public F64Matrix2x2(in F64Matrix2x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+        }
+
+        public F64Matrix2x2(in F32Matrix2x2 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 -F64Utility.EPSYLON <= determinant && determinant <= F64Utility.EPSYLON;
+        }
+
+        public void Transpose()
+        {
+            (this.r1c2, this.r2c1) = (this.r2c1, this.r1c2);
+        }
+
+        public bool Invert()
+        {
+            double determinant = this.GetDeterminant();
+
+            if (-F64Utility.EPSYLON <= determinant && determinant <= F64Utility.EPSYLON)
+            {
+                return false;
+            }
+
+            double r1c1 =  this.r2c2;
+            double r1c2 = -this.r1c2;
+
+            double r2c1 = -this.r2c1;
+            double r2c2 =  this.r1c1;
+
+            this.r1c1 = r1c1 / determinant;
+            this.r1c2 = r1c2 / determinant;
+
+            this.r2c1 = r2c1 / determinant;
+            this.r2c2 = r2c2 / determinant;
+
+            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 F64Matrix2x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+        }
+
+        public void SetValues(in F32Matrix2x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r1c2;
+
+            this.r2c1 = matrix.r2c1;
+            this.r2c2 = matrix.r2c2;
+        }
+
+        public void SetTransposedOf(in F64Matrix2x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r2c2 = matrix.r2c2;
+            (this.r1c2, this.r2c1) = (matrix.r2c1, matrix.r1c2);
+        }
+
+        public void SetTransposedOf(in F32Matrix2x2 matrix)
+        {
+            this.r1c1 = matrix.r1c1;
+            this.r1c2 = matrix.r2c1;
+            
+            this.r2c1 = matrix.r1c2;
+            this.r2c2 = matrix.r2c2;
+        }
+
+        public bool SetInvertedOf(in F64Matrix2x2 matrix)
+        {
+            double determinant = matrix.GetDeterminant();
+
+            if (-F64Utility.EPSYLON <= determinant && determinant <= F64Utility.EPSYLON)
+            {
+                return false;
+            }
+
+            double r1c1 =  matrix.r2c2;
+            double r1c2 = -matrix.r1c2;
+
+            double r2c1 = -matrix.r2c1;
+            double r2c2 =  matrix.r1c1;
+
+            this.r1c1 = r1c1 / determinant;
+            this.r1c2 = r1c2 / determinant;
+
+            this.r2c1 = r2c1 / determinant;
+            this.r2c2 = r2c2 / determinant;
+
+            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 static void Add(in F64Matrix2x2 matrix1, in F64Matrix2x2 matrix2, out F64Matrix2x2 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 F64Matrix2x2 minuend, in F64Matrix2x2 subtrahend, out F64Matrix2x2 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 F64Matrix2x2 multiplicand, double multiplier, out F64Matrix2x2 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 F64Matrix2x2 dividend, double divisor, out F64Matrix2x2 quotient)
+        {
+            quotient.r1c1 = dividend.r1c1 / divisor;
+            quotient.r1c2 = dividend.r1c2 / divisor;
+
+            quotient.r2c1 = dividend.r2c1 / divisor;
+            quotient.r2c2 = dividend.r2c2 / divisor;
+        }
+
+        public static void GetWeightedSum2(
+            double weight1, in F64Matrix2x2 matrix1,
+            double weight2, in F64Matrix2x2 matrix2,
+            out F64Matrix2x2 sum)
+        {
+            sum.r1c1 = matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2;
+            sum.r1c2 = matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2;
+
+            sum.r2c1 = matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2;
+            sum.r2c2 = matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2;
+        }
+
+        public static void GetWeightedSum3(
+            float weight1, in F32Matrix2x2 matrix1,
+            float weight2, in F32Matrix2x2 matrix2,
+            float weight3, in F32Matrix2x2 matrix3,
+            out F64Matrix2x2 sum)
+        {
+            sum.r1c1 = matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2 + matrix3.r1c1 * weight3;
+            sum.r1c2 = matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2 + matrix3.r1c2 * weight3;
+
+            sum.r2c1 = matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2 + matrix3.r2c1 * weight3;
+            sum.r2c2 = matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2 + matrix3.r2c2 * weight3;
+        }
+
+        public static void GetWeightedSum4(
+            float weight1, in F32Matrix2x2 matrix1,
+            float weight2, in F32Matrix2x2 matrix2,
+            float weight3, in F32Matrix2x2 matrix3,
+            float weight4, in F32Matrix2x2 matrix4,
+            out F64Matrix2x2 sum)
+        {
+            sum.r1c1 = (matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2) + (matrix3.r1c1 * weight3 + matrix4.r1c1 * weight4);
+            sum.r1c2 = (matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2) + (matrix3.r1c2 * weight3 + matrix4.r1c2 * weight4);
+
+            sum.r2c1 = (matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2) + (matrix3.r2c1 * weight3 + matrix4.r2c1 * weight4);
+            sum.r2c2 = (matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2) + (matrix3.r2c2 * weight3 + matrix4.r2c2 * weight4);
+        }
+
+        public static void GetWeightedSum5(
+            double weight1, in F64Matrix2x2 matrix1,
+            double weight2, in F64Matrix2x2 matrix2,
+            double weight3, in F64Matrix2x2 matrix3,
+            double weight4, in F64Matrix2x2 matrix4,
+            double weight5, in F64Matrix2x2 matrix5,
+            out F64Matrix2x2 sum)
+        {
+            sum.r1c1 = (matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2) + (matrix3.r1c1 * weight3 + matrix4.r1c1 * weight4) + matrix5.r1c1 * weight5;
+            sum.r1c2 = (matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2) + (matrix3.r1c2 * weight3 + matrix4.r1c2 * weight4) + matrix5.r1c2 * weight5;
+
+            sum.r2c1 = (matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2) + (matrix3.r2c1 * weight3 + matrix4.r2c1 * weight4) + matrix5.r2c1 * weight5;
+            sum.r2c2 = (matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2) + (matrix3.r2c2 * weight3 + matrix4.r2c2 * weight4) + matrix5.r2c2 * weight5;
+        }
+
+        public static void GetRightProduct(in F64Matrix2x2 matrix, in F64Vector2 vector, out F64Vector2 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 F64Vector2 vector, in F64Matrix2x2 matrix, out F64Vector2 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 F64Matrix2x2 matrix)
+        {
+            matrix.r1c1 = 0.0;
+            matrix.r1c2 = 0.0;
+
+            matrix.r2c1 = 0.0;
+            matrix.r2c2 = 0.0;
+        }
+
+        public static void LoadIdentity(out F64Matrix2x2 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 F64Matrix2x2 matrix)
+        {
+            matrix.r1c1 = d1;
+            matrix.r1c2 = 0.0;
+
+            matrix.r2c1 = 0.0;
+            matrix.r2c2 = d2;
+        }
+    }
+}

BGC/F64Matrix2x3.cs

@@ -0,0 +1,310 @@
+/*
+ * 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 BGC
+{
+    public struct F64Matrix2x3
+    {
+        public double r1c1, r1c2;
+        public double r2c1, r2c2;
+        public double r3c1, r3c2;
+    
+        public F64Matrix2x3()
+        {
+            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 F64Matrix2x3(in F64Matrix2x3 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 F64Matrix2x3(in F32Matrix2x3 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 F64Matrix2x3(in F64Matrix3x2 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 F64Matrix2x3(in F32Matrix3x2 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 F64Matrix2x3 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 F32Matrix2x3 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 F64Matrix3x2 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 F32Matrix3x2 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 static void Add(in F64Matrix2x3 matrix1, in F64Matrix2x3 matrix2, out F64Matrix2x3 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 F64Matrix2x3 minuend, in F64Matrix2x3 subtrahend, out F64Matrix2x3 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 F64Matrix2x3 multiplicand, double multiplier, out F64Matrix2x3 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 F64Matrix2x3 dividend, double divisor, out F64Matrix2x3 quotient)
+        {
+            quotient.r1c1 = dividend.r1c1 / divisor;
+            quotient.r1c2 = dividend.r1c2 / divisor;
+
+            quotient.r2c1 = dividend.r2c1 / divisor;
+            quotient.r2c2 = dividend.r2c2 / divisor;
+
+            quotient.r3c1 = dividend.r3c1 / divisor;
+            quotient.r3c2 = dividend.r3c2 / divisor;
+        }
+
+        public static void GetWeightedSum2(
+            double weight1, in F64Matrix2x3 matrix1,
+            double weight2, in F64Matrix2x3 matrix2,
+            out F64Matrix2x3 sum)
+        {
+            sum.r1c1 = matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2;
+            sum.r1c2 = matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2;
+
+            sum.r2c1 = matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2;
+            sum.r2c2 = matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2;
+
+            sum.r3c1 = matrix1.r3c1 * weight1 + matrix2.r3c1 * weight2;
+            sum.r3c2 = matrix1.r3c2 * weight1 + matrix2.r3c2 * weight2;
+        }
+
+        public static void GetWeightedSum3(
+            double weight1, in F64Matrix2x3 matrix1,
+            double weight2, in F64Matrix2x3 matrix2,
+            double weight3, in F64Matrix2x3 matrix3,
+            out F64Matrix2x3 sum)
+        {
+            sum.r1c1 = matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2 + matrix3.r1c1 * weight3;
+            sum.r1c2 = matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2 + matrix3.r1c2 * weight3;
+
+            sum.r2c1 = matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2 + matrix3.r2c1 * weight3;
+            sum.r2c2 = matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2 + matrix3.r2c2 * weight3;
+
+            sum.r3c1 = matrix1.r3c1 * weight1 + matrix2.r3c1 * weight2 + matrix3.r3c1 * weight3;
+            sum.r3c2 = matrix1.r3c2 * weight1 + matrix2.r3c2 * weight2 + matrix3.r3c2 * weight3;
+        }
+
+        public static void GetWeightedSum4(
+            double weight1, in F64Matrix2x3 matrix1,
+            double weight2, in F64Matrix2x3 matrix2,
+            double weight3, in F64Matrix2x3 matrix3,
+            double weight4, in F64Matrix2x3 matrix4,
+            out F64Matrix2x3 sum)
+        {
+            sum.r1c1 = (matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2) + (matrix3.r1c1 * weight3 + matrix4.r1c1 * weight4);
+            sum.r1c2 = (matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2) + (matrix3.r1c2 * weight3 + matrix4.r1c2 * weight4);
+
+            sum.r2c1 = (matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2) + (matrix3.r2c1 * weight3 + matrix4.r2c1 * weight4);
+            sum.r2c2 = (matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2) + (matrix3.r2c2 * weight3 + matrix4.r2c2 * weight4);
+
+            sum.r3c1 = (matrix1.r3c1 * weight1 + matrix2.r3c1 * weight2) + (matrix3.r3c1 * weight3 + matrix4.r3c1 * weight4);
+            sum.r3c2 = (matrix1.r3c2 * weight1 + matrix2.r3c2 * weight2) + (matrix3.r3c2 * weight3 + matrix4.r3c2 * weight4);
+        }
+
+        public static void GetWeightedSum5(
+            double weight1, in F64Matrix2x3 matrix1,
+            double weight2, in F64Matrix2x3 matrix2,
+            double weight3, in F64Matrix2x3 matrix3,
+            double weight4, in F64Matrix2x3 matrix4,
+            double weight5, in F64Matrix2x3 matrix5,
+            out F64Matrix2x3 sum)
+        {
+            sum.r1c1 = (matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2) + (matrix3.r1c1 * weight3 + matrix4.r1c1 * weight4) + matrix5.r1c1 * weight5;
+            sum.r1c2 = (matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2) + (matrix3.r1c2 * weight3 + matrix4.r1c2 * weight4) + matrix5.r1c2 * weight5;
+
+            sum.r2c1 = (matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2) + (matrix3.r2c1 * weight3 + matrix4.r2c1 * weight4) + matrix5.r2c1 * weight5;
+            sum.r2c2 = (matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2) + (matrix3.r2c2 * weight3 + matrix4.r2c2 * weight4) + matrix5.r2c2 * weight5;
+
+            sum.r3c1 = (matrix1.r3c1 * weight1 + matrix2.r3c1 * weight2) + (matrix3.r3c1 * weight3 + matrix4.r3c1 * weight4) + matrix5.r3c1 * weight5;
+            sum.r3c2 = (matrix1.r3c2 * weight1 + matrix2.r3c2 * weight2) + (matrix3.r3c2 * weight3 + matrix4.r3c2 * weight4) + matrix5.r3c2 * weight5;
+        }
+
+        public static void GetRightProduct(in F64Matrix2x3 matrix, in F64Vector2 vector, out F64Vector3 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 F64Vector3 vector, in F64Matrix2x3 matrix, out F64Vector2 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;
+        }
+    }
+}

BGC/F64Matrix3x2.cs

@@ -0,0 +1,292 @@
+/*
+ * 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 BGC
+{
+    public struct F64Matrix3x2
+    {
+        public double r1c1, r1c2, r1c3;
+        public double r2c1, r2c2, r2c3;
+
+        public F64Matrix3x2()
+        {
+            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 F64Matrix3x2(in F64Matrix3x2 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 F64Matrix3x2(in F32Matrix3x2 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 F64Matrix3x2(in F64Matrix2x3 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 F64Matrix3x2(in F32Matrix2x3 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 F64Matrix3x2 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 F32Matrix3x2 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 F64Matrix2x3 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 F32Matrix2x3 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 static void Add(in F64Matrix3x2 matrix1, in F64Matrix3x2 matrix2, out F64Matrix3x2 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 F64Matrix3x2 minuend, in F64Matrix3x2 subtrahend, out F64Matrix3x2 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 F64Matrix3x2 multiplicand, double multiplier, out F64Matrix3x2 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 F64Matrix3x2 dividend, double divisor, out F64Matrix3x2 quotient)
+        {
+            quotient.r1c1 = dividend.r1c1 / divisor;
+            quotient.r1c2 = dividend.r1c2 / divisor;
+            quotient.r1c3 = dividend.r1c3 / divisor;
+
+            quotient.r2c1 = dividend.r2c1 / divisor;
+            quotient.r2c2 = dividend.r2c2 / divisor;
+            quotient.r2c3 = dividend.r2c3 / divisor;
+        }
+
+        public static void GetWeightedSum2(
+            double weight1, in F64Matrix3x2 matrix1,
+            double weight2, in F64Matrix3x2 matrix2,
+            out F64Matrix3x2 sum)
+        {
+            sum.r1c1 = matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2;
+            sum.r1c2 = matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2;
+            sum.r1c3 = matrix1.r1c3 * weight1 + matrix2.r1c3 * weight2;
+
+            sum.r2c1 = matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2;
+            sum.r2c2 = matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2;
+            sum.r2c3 = matrix1.r2c3 * weight1 + matrix2.r2c3 * weight2;
+        }
+
+        public static void GetWeightedSum3(
+            double weight1, in F64Matrix3x2 matrix1,
+            double weight2, in F64Matrix3x2 matrix2,
+            double weight3, in F64Matrix3x2 matrix3,
+            out F64Matrix3x2 sum)
+        {
+            sum.r1c1 = matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2 + matrix3.r1c1 * weight3;
+            sum.r1c2 = matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2 + matrix3.r1c2 * weight3;
+            sum.r1c3 = matrix1.r1c3 * weight1 + matrix2.r1c3 * weight2 + matrix3.r1c3 * weight3;
+
+            sum.r2c1 = matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2 + matrix3.r2c1 * weight3;
+            sum.r2c2 = matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2 + matrix3.r2c2 * weight3;
+            sum.r2c3 = matrix1.r2c3 * weight1 + matrix2.r2c3 * weight2 + matrix3.r2c3 * weight3;
+        }
+
+        public static void GetWeightedSum4(
+            double weight1, in F64Matrix3x2 matrix1,
+            double weight2, in F64Matrix3x2 matrix2,
+            double weight3, in F64Matrix3x2 matrix3,
+            double weight4, in F64Matrix3x2 matrix4,
+            out F64Matrix3x2 sum)
+        {
+            sum.r1c1 = (matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2) + (matrix3.r1c1 * weight3 + matrix4.r1c1 * weight4);
+            sum.r1c2 = (matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2) + (matrix3.r1c2 * weight3 + matrix4.r1c2 * weight4);
+            sum.r1c3 = (matrix1.r1c3 * weight1 + matrix2.r1c3 * weight2) + (matrix3.r1c3 * weight3 + matrix4.r1c3 * weight4);
+
+            sum.r2c1 = (matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2) + (matrix3.r2c1 * weight3 + matrix4.r2c1 * weight4);
+            sum.r2c2 = (matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2) + (matrix3.r2c2 * weight3 + matrix4.r2c2 * weight4);
+            sum.r2c3 = (matrix1.r2c3 * weight1 + matrix2.r2c3 * weight2) + (matrix3.r2c3 * weight3 + matrix4.r2c3 * weight4);
+        }
+
+        public static void GetWeightedSum5(
+            double weight1, in F64Matrix3x2 matrix1,
+            double weight2, in F64Matrix3x2 matrix2,
+            double weight3, in F64Matrix3x2 matrix3,
+            double weight4, in F64Matrix3x2 matrix4,
+            double weight5, in F64Matrix3x2 matrix5,
+            out F64Matrix3x2 sum)
+        {
+            sum.r1c1 = (matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2) + (matrix3.r1c1 * weight3 + matrix4.r1c1 * weight4) + matrix5.r1c1 * weight5;
+            sum.r1c2 = (matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2) + (matrix3.r1c2 * weight3 + matrix4.r1c2 * weight4) + matrix5.r1c2 * weight5;
+            sum.r1c3 = (matrix1.r1c3 * weight1 + matrix2.r1c3 * weight2) + (matrix3.r1c3 * weight3 + matrix4.r1c3 * weight4) + matrix5.r1c3 * weight5;
+
+            sum.r2c1 = (matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2) + (matrix3.r2c1 * weight3 + matrix4.r2c1 * weight4) + matrix5.r2c1 * weight5;
+            sum.r2c2 = (matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2) + (matrix3.r2c2 * weight3 + matrix4.r2c2 * weight4) + matrix5.r2c2 * weight5;
+            sum.r2c3 = (matrix1.r2c3 * weight1 + matrix2.r2c3 * weight2) + (matrix3.r2c3 * weight3 + matrix4.r2c3 * weight4) + matrix5.r2c3 * weight5;
+        }
+
+        public static void GetRightProduct(in F64Matrix3x2 matrix, in F64Vector3 vector, out F64Vector2 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 F64Vector2 vector, in F64Matrix3x2 matrix, out F64Vector3 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;
+        }
+    }
+}

BGC/F64Matrix3x3.cs

@@ -0,0 +1,528 @@
+/*
+ * 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 BGC
+{
+    public struct F64Matrix3x3
+    {
+        public double r1c1, r1c2, r1c3;
+
+        public double r2c1, r2c2, r2c3;
+
+        public double r3c1, r3c2, r3c3;
+
+        public F64Matrix3x3()
+        {
+            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 F64Matrix3x3(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.r3c1 = 0.0;
+            this.r3c2 = 0.0;
+            this.r3c3 = d3;
+        }
+
+        public F64Matrix3x3(in F64Matrix3x3 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 F64Matrix3x3(in F32Matrix3x3 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 -F64Utility.EPSYLON <= determinant && determinant <= F64Utility.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 (-F64Utility.EPSYLON <= determinant && determinant <= F64Utility.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;
+
+            this.r1c1 = r1c1 / determinant;
+            this.r1c2 = r1c2 / determinant;
+            this.r1c3 = r1c3 / determinant;
+
+            this.r2c1 = r2c1 / determinant;
+            this.r2c2 = r2c2 / determinant;
+            this.r2c3 = r2c3 / determinant;
+
+            this.r3c1 = r3c1 / determinant;
+            this.r3c2 = r3c2 / determinant;
+            this.r3c3 = r3c3 / determinant;
+
+            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(F64Matrix3x3 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(F32Matrix3x3 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 F64Matrix3x3 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 F32Matrix3x3 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 F64Matrix3x3 matrix)
+        {
+            double determinant = matrix.GetDeterminant();
+    
+            if (-F64Utility.EPSYLON <= determinant && determinant <= F64Utility.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;
+
+            this.r1c1 = r1c1 / determinant;
+            this.r1c2 = r1c2 / determinant;
+            this.r1c3 = r1c3 / determinant;
+
+            this.r2c1 = r2c1 / determinant;
+            this.r2c2 = r2c2 / determinant;
+            this.r2c3 = r2c3 / determinant;
+
+            this.r3c1 = r3c1 / determinant;
+            this.r3c2 = r3c2 / determinant;
+            this.r3c3 = r3c3 / determinant;
+
+            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 static void Add(in F64Matrix3x3 matrix1, in F64Matrix3x3 matrix2, out F64Matrix3x3 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 F64Matrix3x3 minuend, in F64Matrix3x3 subtrahend, out F64Matrix3x3 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 F64Matrix3x3 multiplicand, double multiplier, out F64Matrix3x3 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 F64Matrix3x3 dividend, double divisor, out F64Matrix3x3 quotient)
+        {
+            quotient.r1c1 = dividend.r1c1 / divisor;
+            quotient.r1c2 = dividend.r1c2 / divisor;
+            quotient.r1c3 = dividend.r1c3 / divisor;
+
+            quotient.r2c1 = dividend.r2c1 / divisor;
+            quotient.r2c2 = dividend.r2c2 / divisor;
+            quotient.r2c3 = dividend.r2c3 / divisor;
+
+            quotient.r3c1 = dividend.r3c1 / divisor;
+            quotient.r3c2 = dividend.r3c2 / divisor;
+            quotient.r3c3 = dividend.r3c3 / divisor;
+        }
+
+        public static void GetWeightedSum2(
+            double weight1, in F64Matrix3x3 matrix1,
+            double weight2, in F64Matrix3x3 matrix2,
+            out F64Matrix3x3 sum)
+        {
+            sum.r1c1 = matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2;
+            sum.r1c2 = matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2;
+            sum.r1c3 = matrix1.r1c3 * weight1 + matrix2.r1c3 * weight2;
+
+            sum.r2c1 = matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2;
+            sum.r2c2 = matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2;
+            sum.r2c3 = matrix1.r2c3 * weight1 + matrix2.r2c3 * weight2;
+
+            sum.r3c1 = matrix1.r3c1 * weight1 + matrix2.r3c1 * weight2;
+            sum.r3c2 = matrix1.r3c2 * weight1 + matrix2.r3c2 * weight2;
+            sum.r3c3 = matrix1.r3c3 * weight1 + matrix2.r3c3 * weight2;
+        }
+
+        public static void GetWeightedSum3(
+            double weight1, in F64Matrix3x3 matrix1,
+            double weight2, in F64Matrix3x3 matrix2,
+            double weight3, in F64Matrix3x3 matrix3,
+            out F64Matrix3x3 sum)
+        {
+            sum.r1c1 = matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2 + matrix3.r1c1 * weight3;
+            sum.r1c2 = matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2 + matrix3.r1c2 * weight3;
+            sum.r1c3 = matrix1.r1c3 * weight1 + matrix2.r1c3 * weight2 + matrix3.r1c3 * weight3;
+
+            sum.r2c1 = matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2 + matrix3.r2c1 * weight3;
+            sum.r2c2 = matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2 + matrix3.r2c2 * weight3;
+            sum.r2c3 = matrix1.r2c3 * weight1 + matrix2.r2c3 * weight2 + matrix3.r2c3 * weight3;
+
+            sum.r3c1 = matrix1.r3c1 * weight1 + matrix2.r3c1 * weight2 + matrix3.r3c1 * weight3;
+            sum.r3c2 = matrix1.r3c2 * weight1 + matrix2.r3c2 * weight2 + matrix3.r3c2 * weight3;
+            sum.r3c3 = matrix1.r3c3 * weight1 + matrix2.r3c3 * weight2 + matrix3.r3c3 * weight3;
+        }
+
+        public static void GetWeightedSum4(
+            double weight1, in F64Matrix3x3 matrix1,
+            double weight2, in F64Matrix3x3 matrix2,
+            double weight3, in F64Matrix3x3 matrix3,
+            double weight4, in F64Matrix3x3 matrix4,
+            out F64Matrix3x3 sum)
+        {
+            sum.r1c1 = (matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2) + (matrix3.r1c1 * weight3 + matrix4.r1c1 * weight4);
+            sum.r1c2 = (matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2) + (matrix3.r1c2 * weight3 + matrix4.r1c2 * weight4);
+            sum.r1c3 = (matrix1.r1c3 * weight1 + matrix2.r1c3 * weight2) + (matrix3.r1c3 * weight3 + matrix4.r1c3 * weight4);
+
+            sum.r2c1 = (matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2) + (matrix3.r2c1 * weight3 + matrix4.r2c1 * weight4);
+            sum.r2c2 = (matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2) + (matrix3.r2c2 * weight3 + matrix4.r2c2 * weight4);
+            sum.r2c3 = (matrix1.r2c3 * weight1 + matrix2.r2c3 * weight2) + (matrix3.r2c3 * weight3 + matrix4.r2c3 * weight4);
+
+            sum.r3c1 = (matrix1.r3c1 * weight1 + matrix2.r3c1 * weight2) + (matrix3.r3c1 * weight3 + matrix4.r3c1 * weight4);
+            sum.r3c2 = (matrix1.r3c2 * weight1 + matrix2.r3c2 * weight2) + (matrix3.r3c2 * weight3 + matrix4.r3c2 * weight4);
+            sum.r3c3 = (matrix1.r3c3 * weight1 + matrix2.r3c3 * weight2) + (matrix3.r3c3 * weight3 + matrix4.r3c3 * weight4);
+        }
+
+        public static void GetWeightedSum5(
+            double weight1, in F64Matrix3x3 matrix1,
+            double weight2, in F64Matrix3x3 matrix2,
+            double weight3, in F64Matrix3x3 matrix3,
+            double weight4, in F64Matrix3x3 matrix4,
+            double weight5, in F64Matrix3x3 matrix5,
+            out F64Matrix3x3 sum)
+        {
+            sum.r1c1 = (matrix1.r1c1 * weight1 + matrix2.r1c1 * weight2) + (matrix3.r1c1 * weight3 + matrix4.r1c1 * weight4) + matrix5.r1c1 * weight5;
+            sum.r1c2 = (matrix1.r1c2 * weight1 + matrix2.r1c2 * weight2) + (matrix3.r1c2 * weight3 + matrix4.r1c2 * weight4) + matrix5.r1c2 * weight5;
+            sum.r1c3 = (matrix1.r1c3 * weight1 + matrix2.r1c3 * weight2) + (matrix3.r1c3 * weight3 + matrix4.r1c3 * weight4) + matrix5.r1c3 * weight5;
+
+            sum.r2c1 = (matrix1.r2c1 * weight1 + matrix2.r2c1 * weight2) + (matrix3.r2c1 * weight3 + matrix4.r2c1 * weight4) + matrix5.r2c1 * weight5;
+            sum.r2c2 = (matrix1.r2c2 * weight1 + matrix2.r2c2 * weight2) + (matrix3.r2c2 * weight3 + matrix4.r2c2 * weight4) + matrix5.r2c2 * weight5;
+            sum.r2c3 = (matrix1.r2c3 * weight1 + matrix2.r2c3 * weight2) + (matrix3.r2c3 * weight3 + matrix4.r2c3 * weight4) + matrix5.r2c3 * weight5;
+
+            sum.r3c1 = (matrix1.r3c1 * weight1 + matrix2.r3c1 * weight2) + (matrix3.r3c1 * weight3 + matrix4.r3c1 * weight4) + matrix5.r3c1 * weight5;
+            sum.r3c2 = (matrix1.r3c2 * weight1 + matrix2.r3c2 * weight2) + (matrix3.r3c2 * weight3 + matrix4.r3c2 * weight4) + matrix5.r3c2 * weight5;
+            sum.r3c3 = (matrix1.r3c3 * weight1 + matrix2.r3c3 * weight2) + (matrix3.r3c3 * weight3 + matrix4.r3c3 * weight4) + matrix5.r3c3 * weight5;
+        }
+
+        public static void GetRightProduct(in F64Matrix3x3 matrix, in F64Vector3 vector, out F64Vector3 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 F64Vector3 vector, in F64Matrix3x3 matrix, out F64Vector3 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 F64Matrix3x3 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 F64Matrix3x3 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 F64Matrix3x3 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;
+        }
+    }
+}

BGC/F64MatrixProduct.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 BGC
+{
+    public class F64MatrixProduct
+    {
+        public static void Get2x2At2x2(in F64Matrix2x2 left, in F64Matrix2x2 right, out F64Matrix2x2 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 F64Matrix2x2 left, in F64Matrix3x2 right, out F64Matrix3x2 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 F64Matrix2x3 left, in F64Matrix2x2 right, out F64Matrix2x3 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 F64Matrix2x3 left, in F64Matrix3x2 right, out F64Matrix3x3 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 F64Matrix3x2 left, in F64Matrix3x3 right, out F64Matrix3x2 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 F64Matrix3x2 left, in F64Matrix2x3 right, out F64Matrix2x2 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 F64Matrix3x3 left, in F64Matrix2x3 right, out F64Matrix2x3 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 F64Matrix3x3 left, in F64Matrix3x3 right, out F64Matrix3x3 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;
+        }
+    }
+}

BGC/F64Quaternion.cs

@@ -0,0 +1,204 @@
+using System;
+using System.Numerics;
+
+namespace BGC
+{
+    public struct F64Quaternion
+    {
+        public double s0, x1, x2, x3;
+
+        public F64Quaternion(double s0, double x1, double x2, double x3)
+        {
+            this.s0 = s0;
+            this.x1 = x1;
+            this.x2 = x2;
+            this.x3 = x3;
+        }
+
+        public F64Quaternion(in F32Quaternion quaternion)
+        {
+            this.s0 = quaternion.s0;
+            this.x1 = quaternion.x1;
+            this.x2 = quaternion.x2;
+            this.x3 = quaternion.x3;
+        }
+
+        public F64Quaternion(in F64Quaternion 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 F32Quaternion quaternion)
+        {
+            this.s0 = quaternion.s0;
+            this.x1 = quaternion.x1;
+            this.x2 = quaternion.x2;
+            this.x3 = quaternion.x3;
+        }
+
+        public void SetValues(in F64Quaternion quaternion)
+        {
+            this.s0 = quaternion.s0;
+            this.x1 = quaternion.x1;
+            this.x2 = quaternion.x2;
+            this.x3 = quaternion.x3;
+        }
+
+        public void SetConjugateOf(in F64Quaternion quaternion)
+        {
+            this.s0 = quaternion.s0;
+            this.x1 = -quaternion.x1;
+            this.x2 = -quaternion.x2;
+            this.x3 = -quaternion.x3;
+        }
+
+        public readonly void MakeRotationMatrix(out F64Matrix3x3 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 (-F64Utility.EPSYLON <= squareModule && squareModule <= F64Utility.EPSYLON)
+            {
+                F64Matrix3x3.LoadIdentity(out matrix);
+                return;
+            }
+
+            double corrector1;
+            double corrector2;
+
+            if (1.0 - F64Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0 + F64Utility.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 F64Matrix3x3 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 (-F64Utility.EPSYLON <= squareModule && squareModule <= F64Utility.EPSYLON)
+            {
+                F64Matrix3x3.LoadIdentity(out matrix);
+                return;
+            }
+
+            double corrector1;
+            double corrector2;
+
+            if (1.0 - F64Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0 + F64Utility.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 F64Quaternion quaternion1, in F64Quaternion quaternion2, out F64Quaternion 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 F64Quaternion minuend, in F64Quaternion subtrahend, out F64Quaternion 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 F64Quaternion left, in F64Quaternion right, out F64Quaternion 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;
+        }
+    }
+}

BGC/F64Radians.cs

@@ -0,0 +1,75 @@
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 18 Nov 2024
+ */
+
+namespace BGC
+{
+    public class F64Radians
+    {
+        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;
+        }
+    }
+}

BGC/F64Turns.cs

@@ -0,0 +1,63 @@
+
+/*
+ * Author: Andrey Pokidov
+ * Date: 18 Nov 2024
+ */
+
+namespace BGC
+{
+    public class F64Turns
+    {
+        public static double TurnsToRadians(double turns)
+        {
+            return turns * F64Radians.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 * F64Radians.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;
+        }
+    }
+}

BGC/F64Utility.cs

@@ -0,0 +1,20 @@
+using System;
+
+namespace BGC
+{
+    public class F64Utility
+    {
+        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;
+    }
+}

BGC/F64Vector2.cs

@@ -0,0 +1,348 @@
+/*
+ * 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 BGC
+{
+    public struct F64Vector2
+    {
+        public static readonly F64Vector2 ZERO = new F64Vector2(0.0, 0.0);
+
+        public double x1;
+        public double x2;
+
+        public F64Vector2(double x1, double x2)
+        {
+            this.x1 = x1;
+            this.x2 = x2;
+        }
+
+        public F64Vector2(in F64Vector2 vector)
+        {
+            this.x1 = vector.x1;
+            this.x2 = vector.x2;
+        }
+
+        public F64Vector2(in F32Vector2 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 - F64Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0 + F64Utility.TWO_EPSYLON)
+            {
+                return 1;
+            }
+
+            if (squareModule <= F64Utility.SQUARE_EPSYLON)
+            {
+                this.Reset();
+                return 0;
+            }
+
+            double module = Math.Sqrt(squareModule);
+
+            this.x1 /= module;
+            this.x2 /= module;
+
+            return 1;
+        }
+
+        public readonly F64Vector2 GetNormalized()
+        {
+            double squareModule = this.GetSquareModule();
+
+            if (1.0 - F64Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0 + F64Utility.TWO_EPSYLON)
+            {
+                return this;
+            }
+
+            if (squareModule <= F64Utility.SQUARE_EPSYLON)
+            {
+                return ZERO;
+            }
+
+            double module = Math.Sqrt(squareModule);
+
+            return new F64Vector2(
+                this.x1 / module,
+                this.x2 / module
+            );
+        }
+
+        public void Reverse()
+        {
+            this.x1 = -this.x1;
+            this.x2 = -this.x2;
+        }
+
+        public readonly bool IsZero()
+        {
+            return this.GetSquareModule() <= F64Utility.SQUARE_EPSYLON;
+        }
+
+        public readonly bool IsUnit()
+        {
+            double squareModule = this.GetSquareModule();
+            return 1.0 - F64Utility.TWO_EPSYLON <= squareModule && squareModule <= F64Utility.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 F64Vector2 vector)
+        {
+            this.x1 = vector.x1;
+            this.x2 = vector.x2;
+        }
+
+        public void SetValues(in F32Vector2 vector)
+        {
+            this.x1 = vector.x1;
+            this.x2 = vector.x2;
+        }
+
+        public void SetReverseOf(in F64Vector2 vector)
+        {
+            this.x1 = -vector.x1;
+            this.x2 = -vector.x2;
+        }
+
+        public void SetReverseOf(in F32Vector2 vector)
+        {
+            this.x1 = -vector.x1;
+            this.x2 = -vector.x2;
+        }
+
+        public readonly override string ToString()
+        {
+            return String.Format("DPVector2({0}, {1})", this.x1, this.x2);
+        }
+
+        public static void Add(in F64Vector2 vector1, in F64Vector2 vector2, out F64Vector2 sum)
+        {
+            sum.x1 = vector1.x1 + vector2.x1;
+            sum.x2 = vector1.x2 + vector2.x2;
+        }
+
+        public static void Subtract(in F64Vector2 minuend, in F64Vector2 subtrahend, out F64Vector2 difference)
+        {
+            difference.x1 = minuend.x1 - subtrahend.x1;
+            difference.x2 = minuend.x2 - subtrahend.x2;
+        }
+
+        public static void Muliply(in F64Vector2 multiplicand, double multiplier, out F64Vector2 product)
+        {
+            product.x1 = multiplicand.x1 * multiplier;
+            product.x2 = multiplicand.x2 * multiplier;
+        }
+
+        public static void Divide(in F64Vector2 dividend, double divisor, out F64Vector2 quotient)
+        {
+            quotient.x1 = dividend.x1 / divisor;
+            quotient.x2 = dividend.x2 / divisor;
+        }
+
+        public static void GetWeightedSum2(
+            double weight1, in F64Vector2 vector1,
+            double weight2, in F64Vector2 vector2,
+            out F64Vector2 sum
+            )
+        {
+            sum.x1 = vector1.x1 * weight1 + vector2.x1 * weight2;
+            sum.x2 = vector1.x2 * weight1 + vector2.x2 * weight2;
+        }
+
+        public static void GetWeightedSum3(
+            double weight1, in F64Vector2 vector1,
+            double weight2, in F64Vector2 vector2,
+            double weight3, in F64Vector2 vector3,
+            out F64Vector2 sum
+            )
+        {
+            sum.x1 = vector1.x1 * weight1 + vector2.x1 * weight2 + vector3.x1 * weight3;
+            sum.x2 = vector1.x2 * weight1 + vector2.x2 * weight2 + vector3.x2 * weight3;
+        }
+
+        public static void GetWeightedSum4(
+            double weight1, in F64Vector2 vector1,
+            double weight2, in F64Vector2 vector2,
+            double weight3, in F64Vector2 vector3,
+            double weight4, in F64Vector2 vector4,
+            out F64Vector2 sum)
+        {
+            sum.x1 = (vector1.x1 * weight1 + vector2.x1 * weight2) + (vector3.x1 * weight3 + vector4.x1 * weight4);
+            sum.x2 = (vector1.x2 * weight1 + vector2.x2 * weight2) + (vector3.x2 * weight3 + vector4.x2 * weight4);
+        }
+
+        public static void GetWeightedSum5(
+            double weight1, in F64Vector2 vector1,
+            double weight2, in F64Vector2 vector2,
+            double weight3, in F64Vector2 vector3,
+            double weight4, in F64Vector2 vector4,
+            double weight5, in F64Vector2 vector5,
+            out F64Vector2 sum)
+        {
+            sum.x1 = (vector1.x1 * weight1 + vector2.x1 * weight2) + (vector3.x1 * weight3 + vector4.x1 * weight4) + vector5.x1 * weight5;
+            sum.x2 = (vector1.x2 * weight1 + vector2.x2 * weight2) + (vector3.x2 * weight3 + vector4.x2 * weight4) + vector5.x2 * weight5;
+        }
+
+        public static void GetMean2(
+            in F64Vector2 vector1,
+            in F64Vector2 vector2,
+            out F64Vector2 result)
+        {
+            result.x1 = (vector1.x1 + vector2.x1) * 0.5;
+            result.x2 = (vector1.x2 + vector2.x2) * 0.5;
+        }
+
+        public static void GetMean3(
+            in F64Vector2 vector1,
+            in F64Vector2 vector2,
+            in F64Vector2 vector3,
+            out F64Vector2 result)
+        {
+            result.x1 = (vector1.x1 + vector2.x1 + vector3.x1) * F64Utility.ONE_THIRD;
+            result.x2 = (vector1.x2 + vector2.x2 + vector3.x2) * F64Utility.ONE_THIRD;
+        }
+
+        public static void GetMean4(
+            in F64Vector2 vector1,
+            in F64Vector2 vector2,
+            in F64Vector2 vector3,
+            in F64Vector2 vector4,
+            out F64Vector2 result)
+        {
+            result.x1 = ((vector1.x1 + vector2.x1) + (vector3.x1 + vector4.x1)) * 0.25;
+            result.x2 = ((vector1.x2 + vector2.x2) + (vector3.x2 + vector4.x2)) * 0.25;
+        }
+
+        public static void GetMean5(
+            in F64Vector2 vector1,
+            in F64Vector2 vector2,
+            in F64Vector2 vector3,
+            in F64Vector2 vector4,
+            in F64Vector2 vector5,
+            out F64Vector2 result)
+        {
+            result.x1 = ((vector1.x1 + vector2.x1) + (vector3.x1 + vector4.x1) + vector5.x1) * 0.2;
+            result.x2 = ((vector1.x2 + vector2.x2) + (vector3.x2 + vector4.x2) + vector5.x2) * 0.2;
+        }
+
+        public static double GetScalarProduct(in F64Vector2 vector1, in F64Vector2 vector2)
+        {
+            return vector1.x1 * vector2.x1 + vector1.x2 * vector2.x2;
+        }
+ 
+        public static double GetCrossProduct(in F64Vector2 vector1, in F64Vector2 vector2)
+        {
+            return vector1.x1 * vector2.x2 - vector1.x2 * vector2.x1;
+        }
+
+        public static double GetAngle(in F64Vector2 vector1, in F64Vector2 vector2, AngleUnit unit)
+        {
+    		double squareModule1 = vector1.GetSquareModule();
+
+    		if (squareModule1 <= F64Utility.SQUARE_EPSYLON)
+            {
+    			return 0.0;
+    		}
+
+    		double squareModule2 = vector2.GetSquareModule();
+
+    		if (squareModule2 <= F64Utility.SQUARE_EPSYLON)
+            {
+    			return 0.0;
+    		}
+
+    		double cosine = F64Vector2.GetScalarProduct(vector1, vector2) / Math.Sqrt(squareModule1 * squareModule2);
+
+    		if (1.0 - F64Utility.EPSYLON <= cosine)
+            {
+                return 0.0;
+    		}
+
+    		if (cosine <= -(1.0 - F64Utility.EPSYLON))
+            {
+                return F64Angle.GetHalfCircle(unit);
+    		}
+    
+    		return F64Radians.ToUnits(Math.Acos(cosine), unit);
+        }
+
+        public static double GetSquareDistance(in F64Vector2 vector1, in F64Vector2 vector2)
+        {
+            double dx1 = vector1.x1 - vector2.x1;
+            double dx2 = vector1.x2 - vector2.x2;
+        
+            return dx1 * dx1 + dx2 * dx2;
+        }
+    
+        public static double GetDistance(in F64Vector2 vector1, in F64Vector2 vector2)
+        {
+            return Math.Sqrt(GetSquareDistance(vector1, vector2));
+        }
+    
+        public static bool AreEqual(in F64Vector2 vector1, in F64Vector2 vector2)
+        {
+            double squareModule1 = vector1.GetSquareModule();
+            double squareModule2 = vector2.GetSquareModule();
+            double squareModule3 = GetSquareDistance(vector1, vector2);
+
+            // 2.0 means dimension amount
+            if (squareModule1 < F64Utility.EPSYLON_EFFECTIVENESS_LIMIT || squareModule2 < F64Utility.EPSYLON_EFFECTIVENESS_LIMIT)
+            {
+                return squareModule3 < (2.0 * F64Utility.SQUARE_EPSYLON);
+            }
+
+            if (squareModule1 <= squareModule2)
+            {
+                return squareModule3 <= (2.0 * F64Utility.SQUARE_EPSYLON) * squareModule2;
+            }
+
+            return squareModule3 <= (2.0 * F64Utility.SQUARE_EPSYLON) * squareModule1;
+        }
+    }
+}

BGC/F64Vector3.cs

@@ -0,0 +1,397 @@
+/*
+ * 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 BGC
+{
+    public struct F64Vector3
+    {
+        public static readonly F64Vector3 ZERO = new F64Vector3(0.0, 0.0, 0.0);
+
+        public double x1;
+        public double x2;
+        public double x3;
+
+        public F64Vector3(double x1, double x2, double x3)
+        {
+            this.x1 = x1;
+            this.x2 = x2;
+            this.x3 = x3;
+        }
+
+        public F64Vector3(in F64Vector3 vector)
+        {
+            this.x1 = vector.x1;
+            this.x2 = vector.x2;
+            this.x3 = vector.x3;
+        }
+
+        public F64Vector3(in F32Vector3 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 - F64Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0 + F64Utility.TWO_EPSYLON)
+            {
+                return 1;
+            }
+
+            if (squareModule <= F64Utility.SQUARE_EPSYLON)
+            {
+                this.Reset();
+                return 0;
+            }
+
+            double module = Math.Sqrt(squareModule);
+
+            this.x1 /= module;
+            this.x2 /= module;
+            this.x3 /= module;
+
+            return 1;
+        }
+
+        public readonly F64Vector3 GetNormalized()
+        {
+            double squareModule = this.GetSquareModule();
+
+            if (1.0 - F64Utility.TWO_EPSYLON <= squareModule && squareModule <= 1.0 + F64Utility.TWO_EPSYLON)
+            {
+                return this;
+            }
+
+            if (squareModule <= F64Utility.SQUARE_EPSYLON)
+            {
+                return ZERO;
+            }
+
+            double module = Math.Sqrt(squareModule);
+
+            return new F64Vector3(
+                this.x1 / module,
+                this.x2 / module,
+                this.x3 / module
+            );
+        }
+
+        public void Reverse()
+        {
+            this.x1 = -this.x1;
+            this.x2 = -this.x2;
+            this.x3 = -this.x3;
+        }
+
+        public readonly bool IsZero()
+        {
+            return this.GetSquareModule() <= F64Utility.SQUARE_EPSYLON;
+        }
+
+        public readonly bool IsUnit()
+        {
+            double squareModule = this.GetSquareModule();
+            return 1.0 - F64Utility.TWO_EPSYLON <= squareModule && squareModule <= F64Utility.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 F64Vector3 vector)
+        {
+            this.x1 = vector.x1;
+            this.x2 = vector.x2;
+            this.x3 = vector.x3;
+        }
+
+        public void SetValues(in F32Vector3 vector)
+        {
+            this.x1 = vector.x1;
+            this.x2 = vector.x2;
+            this.x3 = vector.x3;
+        }
+
+        public void SetReverseOf(in F64Vector3 vector)
+        {
+            this.x1 = -vector.x1;
+            this.x2 = -vector.x2;
+            this.x3 = -vector.x3;
+        }
+
+        public void SetReverseOf(in F32Vector3 vector)
+        {
+            this.x1 = -vector.x1;
+            this.x2 = -vector.x2;
+            this.x3 = -vector.x3;
+        }
+
+        public readonly override string ToString()
+        {
+            return String.Format("DPVector3({0}, {1}, {2})", this.x1, this.x2, this.x3);
+        }
+
+        public static void Add(in F64Vector3 vector1, in F64Vector3 vector2, out F64Vector3 sum)
+        {
+            sum.x1 = vector1.x1 + vector2.x1;
+            sum.x2 = vector1.x2 + vector2.x2;
+            sum.x3 = vector1.x3 + vector2.x3;
+        }
+
+        public static void Subtract(in F64Vector3 minuend, in F64Vector3 subtrahend, out F64Vector3 difference)
+        {
+            difference.x1 = minuend.x1 - subtrahend.x1;
+            difference.x2 = minuend.x2 - subtrahend.x2;
+            difference.x3 = minuend.x3 - subtrahend.x3;
+        }
+
+        public static void Muliply(in F64Vector3 multiplicand, double multiplier, out F64Vector3 product)
+        {
+            product.x1 = multiplicand.x1 * multiplier;
+            product.x2 = multiplicand.x2 * multiplier;
+            product.x3 = multiplicand.x3 * multiplier;
+        }
+
+        public static void Divide(in F64Vector3 dividend, double divisor, out F64Vector3 quotient)
+        {
+            quotient.x1 = dividend.x1 / divisor;
+            quotient.x2 = dividend.x2 / divisor;
+            quotient.x3 = dividend.x3 / divisor;
+        }
+
+        public static void GetWeightedSum2(
+            double weight1, in F64Vector3 vector1,
+            double weight2, in F64Vector3 vector2,
+            out F64Vector3 sum)
+        {
+            sum.x1 = vector1.x1 * weight1 + vector2.x1 * weight2;
+            sum.x2 = vector1.x2 * weight1 + vector2.x2 * weight2;
+            sum.x3 = vector1.x3 * weight1 + vector2.x3 * weight2;
+        }
+
+        public static void GetWeightedSum3(
+            double weight1, in F64Vector3 vector1,
+            double weight2, in F64Vector3 vector2,
+            double weight3, in F64Vector3 vector3,
+            out F64Vector3 sum)
+        {
+            sum.x1 = vector1.x1 * weight1 + vector2.x1 * weight2 + vector3.x1 * weight3;
+            sum.x2 = vector1.x2 * weight1 + vector2.x2 * weight2 + vector3.x2 * weight3;
+            sum.x3 = vector1.x3 * weight1 + vector2.x3 * weight2 + vector3.x3 * weight3;
+        }
+
+        public static void GetWeightedSum4(
+            double weight1, in F64Vector3 vector1,
+            double weight2, in F64Vector3 vector2,
+            double weight3, in F64Vector3 vector3,
+            double weight4, in F64Vector3 vector4,
+            out F64Vector3 sum)
+        {
+            sum.x1 = (vector1.x1 * weight1 + vector2.x1 * weight2) + (vector3.x1 * weight3 + vector4.x1 * weight4);
+            sum.x2 = (vector1.x2 * weight1 + vector2.x2 * weight2) + (vector3.x2 * weight3 + vector4.x2 * weight4);
+            sum.x3 = (vector1.x3 * weight1 + vector2.x3 * weight2) + (vector3.x3 * weight3 + vector4.x3 * weight4);
+        }
+
+        public static void GetWeightedSum5(
+            double weight1, in F64Vector3 vector1,
+            double weight2, in F64Vector3 vector2,
+            double weight3, in F64Vector3 vector3,
+            double weight4, in F64Vector3 vector4,
+            double weight5, in F64Vector3 vector5,
+            out F64Vector3 sum)
+        {
+            sum.x1 = (vector1.x1 * weight1 + vector2.x1 * weight2) + (vector3.x1 * weight3 + vector4.x1 * weight4) + vector5.x1 * weight5;
+            sum.x2 = (vector1.x2 * weight1 + vector2.x2 * weight2) + (vector3.x2 * weight3 + vector4.x2 * weight4) + vector5.x2 * weight5;
+            sum.x3 = (vector1.x3 * weight1 + vector2.x3 * weight2) + (vector3.x3 * weight3 + vector4.x3 * weight4) + vector5.x3 * weight5;
+        }
+
+        public static void GetMean2(
+            in F64Vector3 vector1,
+            in F64Vector3 vector2,
+            out F64Vector3 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 F64Vector3 vector1,
+            in F64Vector3 vector2,
+            in F64Vector3 vector3,
+            out F64Vector3 result)
+        {
+            result.x1 = (vector1.x1 + vector2.x1 + vector3.x1)  * F64Utility.ONE_THIRD;
+            result.x2 = (vector1.x2 + vector2.x2 + vector3.x2)  * F64Utility.ONE_THIRD;
+            result.x3 = (vector1.x3 + vector2.x3 + vector3.x3)  * F64Utility.ONE_THIRD;
+        }
+
+        public static void GetMean4(
+            in F64Vector3 vector1,
+            in F64Vector3 vector2,
+            in F64Vector3 vector3,
+            in F64Vector3 vector4,
+            out F64Vector3 result)
+        {
+            result.x1 = ((vector1.x1 + vector2.x1) + (vector3.x1 + vector4.x1)) * 0.25;
+            result.x2 = ((vector1.x2 + vector2.x2) + (vector3.x2 + vector4.x2)) * 0.25;
+            result.x3 = ((vector1.x3 + vector2.x3) + (vector3.x3 + vector4.x3)) * 0.25;
+        }
+
+        public static void GetMean5(
+            in F64Vector3 vector1,
+            in F64Vector3 vector2,
+            in F64Vector3 vector3,
+            in F64Vector3 vector4,
+            in F64Vector3 vector5,
+            out F64Vector3 result)
+        {
+            result.x1 = ((vector1.x1 + vector2.x1) + (vector3.x1 + vector4.x1) + vector5.x1) * 0.2;
+            result.x2 = ((vector1.x2 + vector2.x2) + (vector3.x2 + vector4.x2) + vector5.x2) * 0.2;
+            result.x3 = ((vector1.x3 + vector2.x3) + (vector3.x3 + vector4.x3) + vector5.x3) * 0.2;
+        }
+
+        public static double GetScalarProduct(in F64Vector3 vector1, in F64Vector3 vector2)
+        {
+            return vector1.x1 * vector2.x1 + vector1.x2 * vector2.x2 + vector1.x3 * vector2.x3;
+        }
+ 
+        public static void GetCrossProduct(in F64Vector3 vector1, in F64Vector3 vector2, out F64Vector3 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 F64Vector3 vector1, in F64Vector3 vector2, in F64Vector3 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 F64Vector3 vector1, in F64Vector3 vector2, in F64Vector3 vector3, out F64Vector3 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 F64Vector3 vector1, in F64Vector3 vector2, AngleUnit unit)
+        {
+    		double squareModule1 = vector1.GetSquareModule();
+
+    		if (squareModule1 <= F64Utility.SQUARE_EPSYLON)
+            {
+    			return 0.0;
+    		}
+
+    		double squareModule2 = vector2.GetSquareModule();
+
+    		if (squareModule2 <= F64Utility.SQUARE_EPSYLON)
+            {
+    			return 0.0;
+    		}
+
+    		double cosine = F64Vector3.GetScalarProduct(vector1, vector2) / Math.Sqrt(squareModule1 * squareModule2);
+
+    		if (1.0 - F64Utility.EPSYLON <= cosine)
+            {
+                return 0.0;
+    		}
+
+    		if (cosine <= -(1.0 - F64Utility.EPSYLON))
+            {
+                return F64Angle.GetHalfCircle(unit);
+    		}
+    
+    		return F64Radians.ToUnits(Math.Acos(cosine), unit);
+        }
+    
+        public static double GetSquareDistance(in F64Vector3 vector1, in F64Vector3 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 F64Vector3 vector1, in F64Vector3 vector2)
+        {
+            return Math.Sqrt(GetSquareDistance(vector1, vector2));
+        }
+    
+        public static bool AreEqual(in F64Vector3 vector1, in F64Vector3 vector2)
+        {
+            double squareModule1 = vector1.GetSquareModule();
+            double squareModule2 = vector2.GetSquareModule();
+            double squareModule3 = GetSquareDistance(vector1, vector2);
+
+            // 3.0 means dimension amount
+            if (squareModule1 < F64Utility.EPSYLON_EFFECTIVENESS_LIMIT || squareModule2 < F64Utility.EPSYLON_EFFECTIVENESS_LIMIT)
+            {
+                return squareModule3 < (3.0 * F64Utility.SQUARE_EPSYLON);
+            }
+
+            if (squareModule1 <= squareModule2)
+            {
+                return squareModule3 <= (3.0 * F64Utility.SQUARE_EPSYLON) * squareModule2;
+            }
+
+            return squareModule3 <= (3.0 * F64Utility.SQUARE_EPSYLON) * squareModule1;
+        }
+    }
+}
+

BGC/F64Versor.cs

@@ -0,0 +1,359 @@
+/*
+ * 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 BGC
+{
+    public struct F64Versor
+    {
+        public static readonly F64Versor IDLE = new F64Versor(1.0, 0.0, 0.0, 0.0, 1.0);
+
+        private double s0;
+        private double x1;
+        private double x2;
+        private double x3;
+        private double corrector;
+
+        private F64Versor(double s0, double x1, double x2, double x3, double corrector)
+        {
+            this.s0 = s0;
+            this.x1 = x1;
+            this.x2 = x2;
+            this.x3 = x3;
+            this.corrector = corrector;
+        }
+
+        public F64Versor()
+        {
+            this.s0 = 1.0;
+            this.x1 = 0.0;
+            this.x2 = 0.0;
+            this.x3 = 0.0;
+            this.corrector = 1.0;
+        }
+
+        public F64Versor(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);
+
+            this.corrector = 2.0 - squareModule;
+
+            if (squareModule < 1.0 - F64Utility.TWO_EPSYLON || 1.0 + F64Utility.TWO_EPSYLON < squareModule)
+            {
+                this.Normalize(squareModule);
+            }
+        }
+
+        public F64Versor(in F64Versor versor)
+        {
+            this.s0 = versor.s0;
+            this.x1 = versor.x1;
+            this.x2 = versor.x2;
+            this.x3 = versor.x3;
+            this.corrector = versor.corrector;
+        }
+
+        public F64Versor(in F32Versor versor)
+        {
+            this.s0 = versor.GetScalar();
+            this.x1 = versor.GetX1();
+            this.x2 = versor.GetX2();
+            this.x3 = versor.GetX3();
+
+            double squareModule = (this.s0 * this.s0 + this.x1 * this.x1) + (this.x2 * this.x2 + this.x3 * this.x3);
+
+            this.corrector = 2.0 - squareModule;
+
+            if (squareModule < 1.0 - F64Utility.TWO_EPSYLON || 1.0 + F64Utility.TWO_EPSYLON < squareModule)
+            {
+                this.Normalize(squareModule);
+            }
+        }
+
+        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 double GetCorrector()
+        {
+            return this.corrector;
+        }
+
+        public readonly bool IsIdle()
+        {
+            return this.s0 <= -(1.0 - F64Utility.EPSYLON) || (1.0 - F64Utility.EPSYLON) <= this.s0;
+        }
+
+        public void Reset()
+        {
+            this.s0 = 1.0;
+            this.x1 = 0.0;
+            this.x2 = 0.0;
+            this.x3 = 0.0;
+            this.corrector = 1.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 - F64Utility.TWO_EPSYLON) || 1.0 - F64Utility.TWO_EPSYLON <= this.s0) {
+                return 0.0;
+            }
+
+            if (-F64Utility.EPSYLON <= this.s0 && this.s0 <= F64Utility.EPSYLON)
+            {
+                return F64Angle.GetHalfCircle(unit);
+            }
+
+            return F64Radians.ToUnits(2.0 * Math.Acos(s0), unit);
+        }
+
+        public readonly void MakeRotationMatrix(out F64Matrix3x3 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 = 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;
+
+            double corrector2 = 2.0 * this.corrector;
+
+            matrix.r1c1 = this.corrector * ((s0s0 + x1x1) - (x2x2 + x3x3));
+            matrix.r2c2 = this.corrector * ((s0s0 + x2x2) - (x1x1 + x3x3));
+            matrix.r3c3 = this.corrector * ((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 F64Matrix3x3 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 = 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;
+
+            double corrector2 = 2.0 * this.corrector;
+
+            matrix.r1c1 = this.corrector * ((s0s0 + x1x1) - (x2x2 + x3x3));
+            matrix.r2c2 = this.corrector * ((s0s0 + x2x2) - (x1x1 + x3x3));
+            matrix.r3c3 = this.corrector * ((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 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);
+
+            this.corrector = 2.0 - squareModule;
+
+            if (squareModule < 1.0 - F64Utility.TWO_EPSYLON || 1.0 + F64Utility.TWO_EPSYLON < squareModule)
+            {
+                this.Normalize(squareModule);
+            }
+        }
+
+        public void SetValues(in F64Versor versor)
+        {
+            this.s0 = versor.s0;
+            this.x1 = versor.x1;
+            this.x2 = versor.x2;
+            this.x3 = versor.x3;
+            this.corrector = versor.corrector;
+        }
+
+        public void SetValues(in F32Versor versor)
+        {
+            this.s0 = versor.GetScalar();
+            this.x1 = versor.GetX1();
+            this.x2 = versor.GetX2();
+            this.x3 = versor.GetX3();
+
+            double squareModule = (this.s0 * this.s0 + this.x1 * this.x1) + (this.x2 * this.x2 + this.x3 * this.x3);
+
+            this.corrector = 2.0 - squareModule;
+
+            if (squareModule < 1.0 - F64Utility.TWO_EPSYLON || 1.0 + F64Utility.TWO_EPSYLON < squareModule)
+            {
+                this.Normalize(squareModule);
+            }
+        }
+
+        public void SetInverted(in F64Versor versor)
+        {
+            this.s0 = versor.s0;
+            this.x1 = -versor.x1;
+            this.x2 = -versor.x2;
+            this.x3 = -versor.x3;
+            this.corrector = versor.corrector;
+        }
+
+        public void SetInverted(in F32Versor versor)
+        {
+            this.s0 = versor.GetScalar();
+            this.x1 = -versor.GetX1();
+            this.x2 = -versor.GetX2();
+            this.x3 = -versor.GetX3();
+
+            double squareModule = (this.s0 * this.s0 + this.x1 * this.x1) + (this.x2 * this.x2 + this.x3 * this.x3);
+
+            this.corrector = 2.0 - squareModule;
+
+            if (squareModule < 1.0 - F64Utility.TWO_EPSYLON || 1.0 + F64Utility.TWO_EPSYLON < squareModule)
+            {
+                this.Normalize(squareModule);
+            }
+        }
+
+        public readonly void Turn(in F64Vector3 vector, out F64Vector3 result)
+        {
+            double multiplier = 2.0 * this.corrector;
+            double tx1 = multiplier * (this.x2 * vector.x3 - this.x3 * vector.x2);
+            double tx2 = multiplier * (this.x3 * vector.x1 - this.x1 * vector.x3);
+            double tx3 = multiplier * (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 F64Vector3 vector, out F64Vector3 result)
+        {
+            double multiplier = 2.0 * this.corrector;
+            double tx1 = multiplier * (this.x2 * vector.x3 - this.x3 * vector.x2);
+            double tx2 = multiplier * (this.x3 * vector.x1 - this.x1 * vector.x3);
+            double tx3 = multiplier * (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 <= F64Utility.SQUARE_EPSYLON || (this.x1 * this.x1 + this.x2 * this.x2 + this.x3 * this.x3) <= F64Utility.SQUARE_EPSYLON * squareModule)
+            {
+                this.Reset();
+                return;
+            }
+
+            double module = Math.Sqrt(squareModule);
+
+            this.s0 /= module;
+            this.x1 /= module;
+            this.x2 /= module;
+            this.x3 /= module;
+
+            this.corrector = 2.0 - (this.s0 * this.s0 + this.x1 * this.x1) - (this.x2 * this.x2 + this.x3 * this.x3);
+        }
+
+        public static void Combine(in F64Versor second, in F64Versor first, out F64Versor 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);
+
+            result.s0 = s0;
+            result.x1 = x1;
+            result.x2 = x2;
+            result.x3 = x3;
+
+            double squareModule = (s0 * s0 + x1 * x1) + (x2 * x2 + x3 * x3);
+
+            result.corrector = 2.0 - squareModule;
+
+            if (squareModule < 1.0 - F64Utility.TWO_EPSYLON || 1.0 + F64Utility.TWO_EPSYLON < squareModule)
+            {
+                result.Normalize(squareModule);
+            }
+        }
+    }
+}