From 42ba33d229a164f4e34b636044c415eb60288ada Mon Sep 17 00:00:00 2001 From: hzhang Date: Fri, 21 Jun 2024 21:48:07 +0800 Subject: [PATCH] M2 --- .gitignore | 5 + Skeleton.csproj | 9 + Skeleton.sln | 16 + global.json | 7 + .../Extensions/ComplexMatrixExtension.cs | 192 +++ .../ComplexNormalMatrixExtension.cs | 42 + .../Extensions/RealDiagonalMatrixExtension.cs | 20 + .../UnitaryInvariantMatrixExtension.cs | 23 + .../IUnitaryInvariantMatrix.cs | 11 + src/Algebra/AffineSpaces/AffineSpace.cs | 163 ++ src/Algebra/AffineSpaces/FAffineSpace.cs | 19 + src/Algebra/CategoryOf.cs | 29 + src/Algebra/DimensionProviders/IDim2.cs | 6 + src/Algebra/DimensionProviders/IDim3.cs | 6 + src/Algebra/DimensionProviders/IDim4.cs | 6 + src/Algebra/DimensionProviders/OfDimension.cs | 18 + src/Algebra/Extensions/GeneralExt.cs | 256 ++++ src/Algebra/FixableObjects/FixableTensor.cs | 50 + .../PostConstructionFixExtension.cs | 11 + src/Algebra/Groups/FUnitaryGroup.cs | 21 + src/Algebra/Groups/Group.cs | 46 + src/Algebra/Matrices/FMatrix.cs | 206 +++ src/Algebra/Matrices/Matrix.cs | 1312 +++++++++++++++++ .../Matrices/MatrixProperty/MatrixProperty.cs | 102 ++ .../NormalMatrices/FDiagonalMatrix.cs | 79 + .../NormalMatrices/FHermitianMatrix.cs | 99 ++ .../NormalMatrices/FLieUnitaryMatrix.cs | 158 ++ .../Matrices/NormalMatrices/FNormalMatrix.cs | 106 ++ .../FSpecialLieUnitaryMatrix.cs | 111 ++ .../NormalMatrices/FSpecialUnitaryMatrix.cs | 94 ++ .../Matrices/NormalMatrices/FUnitaryMatrix.cs | 119 ++ src/Algebra/Morphisms/Homomorphism.cs | 13 + src/Algebra/Morphisms/Isomorphism.cs | 11 + src/Algebra/Morphisms/Morphism.cs | 11 + .../ComplexFieldStructure.cs | 121 ++ .../ScalarFieldStructure/FieldStructure.cs | 325 ++++ .../RealFieldStructure.cs | 86 ++ .../ScalarFieldStructureProvider.cs | 19 + .../ScalarFieldStructure/StaticAccess.cs | 23 + src/Algebra/VectorSpaces/FVectorSpace.cs | 25 + src/Algebra/VectorSpaces/VectorSpace.cs | 319 ++++ src/Algebra/Vectors/FVector.cs | 125 ++ src/Algebra/Vectors/Vector.cs | 499 +++++++ .../Arithmetics/AnalyticAddition.cs | 39 + .../Arithmetics/AnalyticApply.cs | 39 + .../Arithmetics/AnalyticConstant.cs | 46 + .../Arithmetics/AnalyticIdentity.cs | 25 + .../Arithmetics/AnalyticNegation.cs | 36 + .../Arithmetics/AnalyticProduction.cs | 39 + src/Analysis/AnalyticFunctions/IAnalytic.cs | 19 + .../AnalyticFunctions/Implements/Analytic.cs | 63 + .../Polynomials/IPolynomial.cs | 29 + .../Polynomials/Implements/C0Polynomial.cs | 61 + .../Polynomials/Implements/C1Polynomial.cs | 54 + .../Polynomials/Implements/C2Polynomial.cs | 74 + .../Polynomials/Implements/C3Polynomial.cs | 105 ++ .../Polynomials/Implements/C4Polynomial.cs | 140 ++ .../Implements/GeneralPolynomial.cs | 29 + .../Polynomials/Implements/Polynomial.cs | 262 ++++ .../Polynomials/Implements/Xd.cs | 13 + .../BivariantContinuousAddition.cs | 30 + .../BivariantContinuousConstant.cs | 36 + .../BivariantContinuousIdentityX.cs | 17 + .../BivariantContinuousIdentityY.cs | 17 + .../BivariantContinuousNegation.cs | 22 + .../BivariantContinuousProduction.cs | 30 + .../Real/Arithmetics/BivariantContinuousX2.cs | 15 + .../Real/Arithmetics/BivariantContinuousX3.cs | 20 + .../Real/Arithmetics/BivariantContinuousX4.cs | 20 + .../Real/Arithmetics/BivariantContinuousY2.cs | 20 + .../Real/Arithmetics/BivariantContinuousY3.cs | 20 + .../Real/Arithmetics/BivariantContinuousY4.cs | 20 + .../Real/IBivariantContinuous.cs | 23 + .../Real/Implements/BivariantContinuous.cs | 156 ++ .../Real/PolynomialHelpers/PolyCHelper.cs | 435 ++++++ .../Arithmetics/ContinuousAddition.cs | 27 + .../Arithmetics/ContinuousApply.cs | 27 + .../Arithmetics/ContinuousConstant.cs | 33 + .../Arithmetics/ContinuousNegation.cs | 25 + .../Arithmetics/ContinuousProduction.cs | 27 + .../ContinuousFunctions/IContinuous.cs | 18 + .../Implements/Continuous.cs | 58 + src/Constants/AlgebraConstant.cs | 100 ++ src/Constants/CC.cs | 76 + src/DataStructure/CacheItem.cs | 82 ++ src/DataStructure/IsomorphicMap.cs | 60 + src/DataStructure/Link/Link.cs | 236 +++ src/DataStructure/Link/LinkHead.cs | 38 + src/DataStructure/Link/LinkNode.cs | 164 +++ src/DataStructure/Link/LinkTail.cs | 38 + src/DataStructure/Packs/EigenPair.cs | 26 + .../Packs/LinearSystems/LinearSystem.cs | 37 + .../MatrixDecompositions/GeneralUHPair.cs | 22 + .../Packs/MatrixDecompositions/JDPair.cs | 21 + .../Packs/MatrixDecompositions/PLUPack.cs | 37 + .../Packs/MatrixDecompositions/QRPair.cs | 28 + src/DataStructure/Packs/Pack.cs | 124 ++ src/DataStructure/Packs/PlainPack.cs | 78 + .../Packs/PolynomialDivisionResult.cs | 19 + src/DataStructure/Packs/QuadraticRoots.cs | 21 + .../Packs/ReducedRowEchelonForm.cs | 27 + src/Samples/C22Samples.cs | 39 + src/Samples/C33Samples.cs | 51 + src/Samples/C44Samples.cs | 42 + src/Samples/ComplexSamples.cs | 48 + src/Samples/DiagC2Samples.cs | 60 + src/Samples/DiagC3Samples.cs | 63 + src/Samples/DiagC4Samples.cs | 63 + src/Samples/LieSU2Samples.cs | 51 + src/Samples/LieSU3Samples.cs | 51 + src/Samples/LieSU4Samples.cs | 58 + src/Samples/R2Samples.cs | 24 + src/Samples/R3Samples.cs | 24 + src/Samples/R4Samples.cs | 24 + src/Samples/RealSamples.cs | 133 ++ src/Samples/SU2Samples.cs | 82 ++ src/Samples/SU3Samples.cs | 82 ++ src/Samples/SU4Samples.cs | 83 ++ src/Samples/U2Samples.cs | 39 + src/Samples/U3Samples.cs | 39 + src/Samples/U4Samples.cs | 39 + src/Utils/Approximations/ErfApproximation.cs | 105 ++ src/Utils/Helpers/AlgebraHelper.cs | 152 ++ src/Utils/Helpers/Decorators.cs | 54 + src/Utils/Helpers/LinqHelper.cs | 564 +++++++ src/Utils/InverseSampling/InverseSampling.cs | 40 + src/Utils/RandomEngines/Normal.cs | 35 + src/Utils/RandomEngines/RandomSource.cs | 18 + src/Utils/Utils.cs | 326 ++++ 129 files changed, 10711 insertions(+) create mode 100644 .gitignore create mode 100644 Skeleton.csproj create mode 100644 Skeleton.sln create mode 100644 global.json create mode 100644 src/Algebra/AdditionalProperties/Extensions/ComplexMatrixExtension.cs create mode 100644 src/Algebra/AdditionalProperties/Extensions/ComplexNormalMatrixExtension.cs create mode 100644 src/Algebra/AdditionalProperties/Extensions/RealDiagonalMatrixExtension.cs create mode 100644 src/Algebra/AdditionalProperties/Extensions/UnitaryInvariantMatrixExtension.cs create mode 100644 src/Algebra/AdditionalProperties/IUnitaryInvariantMatrix.cs create mode 100644 src/Algebra/AffineSpaces/AffineSpace.cs create mode 100644 src/Algebra/AffineSpaces/FAffineSpace.cs create mode 100644 src/Algebra/CategoryOf.cs create mode 100644 src/Algebra/DimensionProviders/IDim2.cs create mode 100644 src/Algebra/DimensionProviders/IDim3.cs create mode 100644 src/Algebra/DimensionProviders/IDim4.cs create mode 100644 src/Algebra/DimensionProviders/OfDimension.cs create mode 100644 src/Algebra/Extensions/GeneralExt.cs create mode 100644 src/Algebra/FixableObjects/FixableTensor.cs create mode 100644 src/Algebra/FixableObjects/PostConstructionFixExtension.cs create mode 100644 src/Algebra/Groups/FUnitaryGroup.cs create mode 100644 src/Algebra/Groups/Group.cs create mode 100644 src/Algebra/Matrices/FMatrix.cs create mode 100644 src/Algebra/Matrices/Matrix.cs create mode 100644 src/Algebra/Matrices/MatrixProperty/MatrixProperty.cs create mode 100644 src/Algebra/Matrices/NormalMatrices/FDiagonalMatrix.cs create mode 100644 src/Algebra/Matrices/NormalMatrices/FHermitianMatrix.cs create mode 100644 src/Algebra/Matrices/NormalMatrices/FLieUnitaryMatrix.cs create mode 100644 src/Algebra/Matrices/NormalMatrices/FNormalMatrix.cs create mode 100644 src/Algebra/Matrices/NormalMatrices/FSpecialLieUnitaryMatrix.cs create mode 100644 src/Algebra/Matrices/NormalMatrices/FSpecialUnitaryMatrix.cs create mode 100644 src/Algebra/Matrices/NormalMatrices/FUnitaryMatrix.cs create mode 100644 src/Algebra/Morphisms/Homomorphism.cs create mode 100644 src/Algebra/Morphisms/Isomorphism.cs create mode 100644 src/Algebra/Morphisms/Morphism.cs create mode 100644 src/Algebra/ScalarFieldStructure/ComplexFieldStructure.cs create mode 100644 src/Algebra/ScalarFieldStructure/FieldStructure.cs create mode 100644 src/Algebra/ScalarFieldStructure/RealFieldStructure.cs create mode 100644 src/Algebra/ScalarFieldStructure/ScalarFieldStructureProvider.cs create mode 100644 src/Algebra/ScalarFieldStructure/StaticAccess.cs create mode 100644 src/Algebra/VectorSpaces/FVectorSpace.cs create mode 100644 src/Algebra/VectorSpaces/VectorSpace.cs create mode 100644 src/Algebra/Vectors/FVector.cs create mode 100644 src/Algebra/Vectors/Vector.cs create mode 100644 src/Analysis/AnalyticFunctions/Arithmetics/AnalyticAddition.cs create mode 100644 src/Analysis/AnalyticFunctions/Arithmetics/AnalyticApply.cs create mode 100644 src/Analysis/AnalyticFunctions/Arithmetics/AnalyticConstant.cs create mode 100644 src/Analysis/AnalyticFunctions/Arithmetics/AnalyticIdentity.cs create mode 100644 src/Analysis/AnalyticFunctions/Arithmetics/AnalyticNegation.cs create mode 100644 src/Analysis/AnalyticFunctions/Arithmetics/AnalyticProduction.cs create mode 100644 src/Analysis/AnalyticFunctions/IAnalytic.cs create mode 100644 src/Analysis/AnalyticFunctions/Implements/Analytic.cs create mode 100644 src/Analysis/AnalyticFunctions/Polynomials/IPolynomial.cs create mode 100644 src/Analysis/AnalyticFunctions/Polynomials/Implements/C0Polynomial.cs create mode 100644 src/Analysis/AnalyticFunctions/Polynomials/Implements/C1Polynomial.cs create mode 100644 src/Analysis/AnalyticFunctions/Polynomials/Implements/C2Polynomial.cs create mode 100644 src/Analysis/AnalyticFunctions/Polynomials/Implements/C3Polynomial.cs create mode 100644 src/Analysis/AnalyticFunctions/Polynomials/Implements/C4Polynomial.cs create mode 100644 src/Analysis/AnalyticFunctions/Polynomials/Implements/GeneralPolynomial.cs create mode 100644 src/Analysis/AnalyticFunctions/Polynomials/Implements/Polynomial.cs create mode 100644 src/Analysis/AnalyticFunctions/Polynomials/Implements/Xd.cs create mode 100644 src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousAddition.cs create mode 100644 src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousConstant.cs create mode 100644 src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousIdentityX.cs create mode 100644 src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousIdentityY.cs create mode 100644 src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousNegation.cs create mode 100644 src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousProduction.cs create mode 100644 src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousX2.cs create mode 100644 src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousX3.cs create mode 100644 src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousX4.cs create mode 100644 src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousY2.cs create mode 100644 src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousY3.cs create mode 100644 src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousY4.cs create mode 100644 src/Analysis/BivariantFunctions/Real/IBivariantContinuous.cs create mode 100644 src/Analysis/BivariantFunctions/Real/Implements/BivariantContinuous.cs create mode 100644 src/Analysis/BivariantFunctions/Real/PolynomialHelpers/PolyCHelper.cs create mode 100644 src/Analysis/ContinuousFunctions/Arithmetics/ContinuousAddition.cs create mode 100644 src/Analysis/ContinuousFunctions/Arithmetics/ContinuousApply.cs create mode 100644 src/Analysis/ContinuousFunctions/Arithmetics/ContinuousConstant.cs create mode 100644 src/Analysis/ContinuousFunctions/Arithmetics/ContinuousNegation.cs create mode 100644 src/Analysis/ContinuousFunctions/Arithmetics/ContinuousProduction.cs create mode 100644 src/Analysis/ContinuousFunctions/IContinuous.cs create mode 100644 src/Analysis/ContinuousFunctions/Implements/Continuous.cs create mode 100644 src/Constants/AlgebraConstant.cs create mode 100644 src/Constants/CC.cs create mode 100644 src/DataStructure/CacheItem.cs create mode 100644 src/DataStructure/IsomorphicMap.cs create mode 100644 src/DataStructure/Link/Link.cs create mode 100644 src/DataStructure/Link/LinkHead.cs create mode 100644 src/DataStructure/Link/LinkNode.cs create mode 100644 src/DataStructure/Link/LinkTail.cs create mode 100644 src/DataStructure/Packs/EigenPair.cs create mode 100644 src/DataStructure/Packs/LinearSystems/LinearSystem.cs create mode 100644 src/DataStructure/Packs/MatrixDecompositions/GeneralUHPair.cs create mode 100644 src/DataStructure/Packs/MatrixDecompositions/JDPair.cs create mode 100644 src/DataStructure/Packs/MatrixDecompositions/PLUPack.cs create mode 100644 src/DataStructure/Packs/MatrixDecompositions/QRPair.cs create mode 100644 src/DataStructure/Packs/Pack.cs create mode 100644 src/DataStructure/Packs/PlainPack.cs create mode 100644 src/DataStructure/Packs/PolynomialDivisionResult.cs create mode 100644 src/DataStructure/Packs/QuadraticRoots.cs create mode 100644 src/DataStructure/Packs/ReducedRowEchelonForm.cs create mode 100644 src/Samples/C22Samples.cs create mode 100644 src/Samples/C33Samples.cs create mode 100644 src/Samples/C44Samples.cs create mode 100644 src/Samples/ComplexSamples.cs create mode 100644 src/Samples/DiagC2Samples.cs create mode 100644 src/Samples/DiagC3Samples.cs create mode 100644 src/Samples/DiagC4Samples.cs create mode 100644 src/Samples/LieSU2Samples.cs create mode 100644 src/Samples/LieSU3Samples.cs create mode 100644 src/Samples/LieSU4Samples.cs create mode 100644 src/Samples/R2Samples.cs create mode 100644 src/Samples/R3Samples.cs create mode 100644 src/Samples/R4Samples.cs create mode 100644 src/Samples/RealSamples.cs create mode 100644 src/Samples/SU2Samples.cs create mode 100644 src/Samples/SU3Samples.cs create mode 100644 src/Samples/SU4Samples.cs create mode 100644 src/Samples/U2Samples.cs create mode 100644 src/Samples/U3Samples.cs create mode 100644 src/Samples/U4Samples.cs create mode 100644 src/Utils/Approximations/ErfApproximation.cs create mode 100644 src/Utils/Helpers/AlgebraHelper.cs create mode 100644 src/Utils/Helpers/Decorators.cs create mode 100644 src/Utils/Helpers/LinqHelper.cs create mode 100644 src/Utils/InverseSampling/InverseSampling.cs create mode 100644 src/Utils/RandomEngines/Normal.cs create mode 100644 src/Utils/RandomEngines/RandomSource.cs create mode 100644 src/Utils/Utils.cs diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..add57be --- /dev/null +++ b/.gitignore @@ -0,0 +1,5 @@ +bin/ +obj/ +/packages/ +riderModule.iml +/_ReSharper.Caches/ \ No newline at end of file diff --git a/Skeleton.csproj b/Skeleton.csproj new file mode 100644 index 0000000..eb2460e --- /dev/null +++ b/Skeleton.csproj @@ -0,0 +1,9 @@ + + + + net6.0 + enable + enable + + + diff --git a/Skeleton.sln b/Skeleton.sln new file mode 100644 index 0000000..9ad9771 --- /dev/null +++ b/Skeleton.sln @@ -0,0 +1,16 @@ + +Microsoft Visual Studio Solution File, Format Version 12.00 +Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "Skeleton", "Skeleton.csproj", "{C073F473-21EF-4B03-8CED-253124F292A1}" +EndProject +Global + GlobalSection(SolutionConfigurationPlatforms) = preSolution + Debug|Any CPU = Debug|Any CPU + Release|Any CPU = Release|Any CPU + EndGlobalSection + GlobalSection(ProjectConfigurationPlatforms) = postSolution + {C073F473-21EF-4B03-8CED-253124F292A1}.Debug|Any CPU.ActiveCfg = Debug|Any CPU + {C073F473-21EF-4B03-8CED-253124F292A1}.Debug|Any CPU.Build.0 = Debug|Any CPU + {C073F473-21EF-4B03-8CED-253124F292A1}.Release|Any CPU.ActiveCfg = Release|Any CPU + {C073F473-21EF-4B03-8CED-253124F292A1}.Release|Any CPU.Build.0 = Release|Any CPU + EndGlobalSection +EndGlobal diff --git a/global.json b/global.json new file mode 100644 index 0000000..1bcf6c0 --- /dev/null +++ b/global.json @@ -0,0 +1,7 @@ +{ + "sdk": { + "version": "6.0.0", + "rollForward": "latestMinor", + "allowPrerelease": false + } +} \ No newline at end of file diff --git a/src/Algebra/AdditionalProperties/Extensions/ComplexMatrixExtension.cs b/src/Algebra/AdditionalProperties/Extensions/ComplexMatrixExtension.cs new file mode 100644 index 0000000..f44a753 --- /dev/null +++ b/src/Algebra/AdditionalProperties/Extensions/ComplexMatrixExtension.cs @@ -0,0 +1,192 @@ +using System.Numerics; +using Skeleton.Algebra.DimensionProviders; + +namespace Skeleton.Algebra.AdditionalProperties.Extensions; +using C22 = CategoryOf.OnField.FMatrix; +using C33 = CategoryOf.OnField.FMatrix; + +/// +/// extensions for complex matrix +/// +public static class ComplexMatrixExtension +{ + + private struct C22B + { + public C22B(Complex a, Complex b, Complex c, Complex d) + { + X11 = a; + X12 = b; + X21 = c; + X22 = d; + } + public Complex X11; + public Complex X12; + public Complex X21; + public Complex X22; + + public static C22B operator +(C22B a, C22B b) => + new(a.X11 + b.X11, a.X12 + b.X12, a.X21 + b.X21, a.X22 + b.X22); + public static C22B operator -(C22B a, C22B b) => + new(a.X11 - b.X11, a.X12 - b.X12, a.X21 - b.X21, a.X22 - b.X22); + + public string PythonString => new C22(X11, X12, X21, X22).PythonRepresentation; + } + + internal static C22 FastMul(C22 self, C22 other) + { + C22 res = new(); + Complex a = self[1, 1]; + Complex b = self[1, 2]; + Complex c = self[2, 1]; + Complex d = self[2, 2]; + Complex A = other[1, 1]; + Complex B = other[2, 1]; + Complex C = other[1, 2]; + Complex D = other[2, 2]; + + Complex CsD = C - D; + Complex cSa = c - a; + Complex t = a * A; + Complex u = cSa * CsD; + Complex v = (c + d) * (C - A); + Complex w = t + (cSa + d) * (A - CsD); + Complex wPv = w + v; + res[1, 1] = t + b* B; + res[1, 2] = wPv + (b - d - cSa) * D; + Complex wPu = w + u; + Complex r21f = B - A + CsD; + res[2, 1] = wPu + d * r21f; + res[2, 2] = wPu + v; + return res; + } + + private static C22B FastFullMul(C22B a, C22B b) + { + Complex CsD = b.X12 - b.X22; + Complex cSa = a.X21 - a.X11; + Complex t = a.X11 * b.X11; + Complex u = cSa * CsD; + Complex v = (a.X21 + a.X22) * (b.X12 - b.X11); + Complex w = t + (cSa + a.X22) * (b.X11 - CsD); + Complex wPv = w + v; + Complex wPu = w + u; + Complex r21f = b.X21 - b.X11 + CsD; + return new(t + a.X12 * b.X21, wPv + (a.X12 - a.X22 - cSa) * b.X22, wPu + a.X22 * r21f, wPu + v); + } + + private static C22B FastMulP1(C22B a, C22B b) => + new(a.X11 * b.X11 + a.X12 * b.X21, 0, a.X21 * b.X11 + a.X22 * b.X21, 0); + + private static C22B FastMulP2(C22B a, C22B b) => + new(a.X11 * b.X11 + a.X12 * b.X21, a.X11 * b.X12 + a.X12 * b.X22, 0, 0); + + internal static C33 FastMul(C33 self, C33 other) + { + C22B a = new(self[1, 1], self[1, 2], self[2, 1], self[2, 2]); + C22B A = new(other[1, 1], other[1, 2], other[2, 1], other[2, 2]); + C22B CsD = new(other[1, 3] - other[3, 3], 0, other[2, 3], 0); + C22B cSa = new(self[3, 1] - self[1, 1], self[3, 2] - self[1, 2], -self[2, 1], -self[2, 2]); + C22B CsA = new(other[1, 3] - other[1, 1], -other[1, 2], other[2, 3] - other[2, 1], -other[2, 2]); + C22B cPd = new(self[3, 1] + self[3, 3], self[3, 2], 0, 0); + C22B Ds_CsA = new(other[3, 3] - other[1, 3] + other[1, 1], other[1, 2], other[2, 1] - other[2, 3], other[2, 2]); + C22B t = FastFullMul(a, A); + C22B u = FastMulP1(cSa, CsD); + C22B v = FastMulP2(cPd, CsA); + C22B w = t + FastFullMul(cPd - a, Ds_CsA); + C22B wPv = w + v; + + C33 res = new(); + res[1, 1] = t.X11 + self[1, 3] * other[3, 1]; + res[1, 2] = t.X12 + self[1, 3] * other[3, 2]; + res[2, 1] = t.X21 + self[2, 3] * other[3, 1]; + res[2, 2] = t.X22 + self[2, 3] * other[3, 2]; + + res[1, 3] = wPv.X11 + (self[1, 3] - self[3, 3] - cSa.X11) * other[3, 3]; + res[2, 3] = wPv.X21 + (self[2, 3] - cSa.X21) * other[3, 3]; + C22B wPu = w + u; + res[3, 1] = wPu.X11 + self[3, 3] * (other[3, 1] - other[1, 1] + CsD.X11); + res[3, 2] = wPu.X12 + self[3, 3] * (other[3, 2] - other[1, 2] + CsD.X12); + res[3, 3] = wPu.X11 + v.X11; + return res; + } + + internal static void FastSelfMul(C33 self, C33 other) + { + C22B a = new(self[1, 1], self[1, 2], self[2, 1], self[2, 2]); + C22B A = new(other[1, 1], other[1, 2], other[2, 1], other[2, 2]); + C22B CsD = new(other[1, 3] - other[3, 3], 0, other[2, 3], 0); + C22B cSa = new(self[3, 1] - self[1, 1], self[3, 2] - self[1, 2], -self[2, 1], -self[2, 2]); + C22B CsA = new(other[1, 3] - other[1, 1], -other[1, 2], other[2, 3] - other[2, 1], -other[2, 2]); + C22B cPd = new(self[3, 1] + self[3, 3], self[3, 2], 0, 0); + C22B Ds_CsA = new(other[3, 3] - other[1, 3] + other[1, 1], other[1, 2], other[2, 1] - other[2, 3], other[2, 2]); + C22B t = FastFullMul(a, A); + C22B u = FastMulP1(cSa, CsD); + C22B v = FastMulP2(cPd, CsA); + C22B w = t + FastFullMul(cPd - a, Ds_CsA); + C22B wPv = w + v; + self[1, 1] = t.X11 + self[1, 3] * other[3, 1]; + self[1, 2] = t.X12 + self[1, 3] * other[3, 2]; + self[2, 1] = t.X21 + self[2, 3] * other[3, 1]; + self[2, 2] = t.X22 + self[2, 3] * other[3, 2]; + + self[1, 3] = wPv.X11 + (self[1, 3] - self[3, 3] - cSa.X11) * other[3, 3]; + self[2, 3] = wPv.X21 + (self[2, 3] - cSa.X21) * other[3, 3]; + C22B wPu = w + u; + self[3, 1] = wPu.X11 + self[3, 3] * (other[3, 1] - other[1, 1] + CsD.X11); + self[3, 2] = wPu.X12 + self[3, 3] * (other[3, 2] - other[1, 2] + CsD.X12); + self[3, 3] = wPu.X11 + v.X11; + } + + internal static void FastSelfLMul(C33 self, C33 other) + { + C22B a = new(self[1, 1], self[1, 2], self[2, 1], self[2, 2]); + C22B A = new(other[1, 1], other[1, 2], other[2, 1], other[2, 2]); + C22B CsD = new(other[1, 3] - other[3, 3], 0, other[2, 3], 0); + C22B cSa = new(self[3, 1] - self[1, 1], self[3, 2] - self[1, 2], -self[2, 1], -self[2, 2]); + C22B CsA = new(other[1, 3] - other[1, 1], -other[1, 2], other[2, 3] - other[2, 1], -other[2, 2]); + C22B cPd = new(self[3, 1] + self[3, 3], self[3, 2], 0, 0); + C22B Ds_CsA = new(other[3, 3] - other[1, 3] + other[1, 1], other[1, 2], other[2, 1] - other[2, 3], other[2, 2]); + C22B t = FastFullMul(a, A); + C22B u = FastMulP1(cSa, CsD); + C22B v = FastMulP2(cPd, CsA); + C22B w = t + FastFullMul(cPd - a, Ds_CsA); + C22B wPv = w + v; + + other[1, 1] = t.X11 + self[1, 3] * other[3, 1]; + other[1, 2] = t.X12 + self[1, 3] * other[3, 2]; + other[2, 1] = t.X21 + self[2, 3] * other[3, 1]; + other[2, 2] = t.X22 + self[2, 3] * other[3, 2]; + + other[1, 3] = wPv.X11 + (self[1, 3] - self[3, 3] - cSa.X11) * other[3, 3]; + other[2, 3] = wPv.X21 + (self[2, 3] - cSa.X21) * other[3, 3]; + C22B wPu = w + u; + other[3, 1] = wPu.X11 + self[3, 3] * (other[3, 1] - A.X11 + CsD.X11); + other[3, 2] = wPu.X12 + self[3, 3] * (other[3, 2] - A.X12 + CsD.X12); + other[3, 3] = wPu.X11 + v.X11; + } + + internal static void UnitaryFix(this CategoryOf.OnField.FMatrix a) + { + if (a.Property.IsDiagonal) + { + for (int i = 1; i <= a.Dim; i++) + a[i, i] /= (a[i, i] * Complex.Conjugate(a[i, i])); + + } + else + { + a[1] = a[1].Normalize; + for (int i = 2; i <= a.Dim; i++) + { + for (int j = 1; j < i; j++) + a[i] -= a[i].ProjectionOn(a[j]); + a[i] = a[i].Normalize; + } + } + } + + + +} + diff --git a/src/Algebra/AdditionalProperties/Extensions/ComplexNormalMatrixExtension.cs b/src/Algebra/AdditionalProperties/Extensions/ComplexNormalMatrixExtension.cs new file mode 100644 index 0000000..4520572 --- /dev/null +++ b/src/Algebra/AdditionalProperties/Extensions/ComplexNormalMatrixExtension.cs @@ -0,0 +1,42 @@ +using System.Numerics; +using Skeleton.DataStructure.Packs.MatrixDecompositions; + +namespace Skeleton.Algebra.AdditionalProperties.Extensions; + +/// +/// additional methods for complex normal matrix +/// +public static class ComplexNormalMatrixExtension +{ + /// + /// unitary decomposition of complex normal matrix + /// + /// + /// + /// + public static JDPair.OnField.FMatrix> UnitaryDecomposition + (this CategoryOf.OnField.FNormalMatrix a, bool order = false) + + { + JDPair.OnField.FMatrix> res = a.InternalUnitaryDecomposition(); + if (!order) + return res; + bool sorted = false; + while (!sorted) + { + sorted = true; + for (int i = 1; i < a.Dim; i++) + { + if (res.D[i, i].Imaginary > res.D[i + 1, i + 1].Imaginary) + { + sorted = false; + res.D.SwapRow(i, i + 1); + res.D.SwapColumn(i, i + 1); + res.J.SwapColumn(i, i + 1); + } + } + } + + return res; + } +} \ No newline at end of file diff --git a/src/Algebra/AdditionalProperties/Extensions/RealDiagonalMatrixExtension.cs b/src/Algebra/AdditionalProperties/Extensions/RealDiagonalMatrixExtension.cs new file mode 100644 index 0000000..5f431c3 --- /dev/null +++ b/src/Algebra/AdditionalProperties/Extensions/RealDiagonalMatrixExtension.cs @@ -0,0 +1,20 @@ +using System.Numerics; + +namespace Skeleton.Algebra.AdditionalProperties.Extensions; + +/// +/// additional methods for real diagonal matrices +/// +public static class RealDiagonalMatrixExtension +{ + /// + /// rayleigh quotient of real diagonal matrix + /// + /// + public static double Rayleigh + (this CategoryOf.OnField.FDiagonalMatrix a, CategoryOf.OnField.FVector b) => + a + .ComplexCast.OnField.FVector, CategoryOf.OnField.FMatrix>() + .MatrixCast.FHermitianMatrix>() + .Rayleigh(b); +} \ No newline at end of file diff --git a/src/Algebra/AdditionalProperties/Extensions/UnitaryInvariantMatrixExtension.cs b/src/Algebra/AdditionalProperties/Extensions/UnitaryInvariantMatrixExtension.cs new file mode 100644 index 0000000..ac04d14 --- /dev/null +++ b/src/Algebra/AdditionalProperties/Extensions/UnitaryInvariantMatrixExtension.cs @@ -0,0 +1,23 @@ +using System.Numerics; + +namespace Skeleton.Algebra.AdditionalProperties.Extensions; + +/// +/// additional methods for unitary invariant matrices +/// +public static class UnitaryInvariantMatrixExtension +{ + /// + /// Conjugated by + /// + /// + /// + /// + /// + /// + public static TInv + ConjugatedBy(this IUnitaryInvariantMatrix a, CategoryOf.FUnitaryMatrix b) + where TInv : CategoryOf.OnField.FMatrix, IUnitaryInvariantMatrix, new() + => b.ConjugateOn((a as TInv)!); + +} \ No newline at end of file diff --git a/src/Algebra/AdditionalProperties/IUnitaryInvariantMatrix.cs b/src/Algebra/AdditionalProperties/IUnitaryInvariantMatrix.cs new file mode 100644 index 0000000..1604713 --- /dev/null +++ b/src/Algebra/AdditionalProperties/IUnitaryInvariantMatrix.cs @@ -0,0 +1,11 @@ +namespace Skeleton.Algebra.AdditionalProperties; + +/// +/// matrix property invariant under unitary similarity +/// +public interface IUnitaryInvariantMatrix { } + +/// +public interface IUnitaryInvariantMatrix : + IUnitaryInvariantMatrix + where TUnitaryInvariant : IUnitaryInvariantMatrix,new() { } diff --git a/src/Algebra/AffineSpaces/AffineSpace.cs b/src/Algebra/AffineSpaces/AffineSpace.cs new file mode 100644 index 0000000..e923f11 --- /dev/null +++ b/src/Algebra/AffineSpaces/AffineSpace.cs @@ -0,0 +1,163 @@ +using Skeleton.Algebra.FixableObjects; +using Skeleton.Algebra.Matrices; +using Skeleton.Algebra.Vectors; +using Skeleton.Algebra.VectorSpaces; +using Skeleton.DataStructure.Packs; +using Skeleton.DataStructure.Packs.LinearSystems; + +namespace Skeleton.Algebra.AffineSpaces; + +/// +/// general affine space +/// +public abstract class AffineSpace : FixableTensor +{ +} + +/// +public abstract class AffineSpace : AffineSpace + where TAffine : AffineSpace, new() +{ + /// + /// intersection of two affine space + /// + /// + /// + public abstract TAffine Intersection(TAffine otherSpace); + +} + +/// +public abstract class AffineSpace : AffineSpace + where TSpace : VectorSpace, new() + where TAffine : AffineSpace, new() +{ + + + /// + /// underlying vector space + /// + public TSpace BaseSpace { get; set; } = new(); +} + +/// +public abstract class AffineSpace : AffineSpace + where TMatrix : Matrix, new() + where TSpace : VectorSpace, new() + where TAffine : AffineSpace, new() +{ + +} + +/// +public abstract class AffineSpace : AffineSpace + where TVector : Vector, new() + where TMatrix : Matrix, new() + where TSpace : VectorSpace, new() + where TAffine : AffineSpace, new() +{ + + /// + /// bias vector + /// + public TVector Bias { get; set; } = new(); + /// + /// if a point is container in the space + /// + /// + /// + public bool ContainsPoint(TVector p) => BaseSpace.ContainsVector(p - Bias); + + /// + /// construct + /// + /// + /// + public void Fill(TSpace s, TVector b) + { + BaseSpace = s; + Bias = b; + } + + /// + /// generate affine space from a space and a bias vector + /// + /// + /// + /// + public static TAffine Dispatch(TSpace s, TVector b) + { + TAffine res = new (); + res.Fill(s, b); + return res; + } + + /// + /// generate affine space from a set of basis and a bias vector + /// + /// + /// + /// + public static TAffine Dispatch(IEnumerable basis, TVector b) + { + TSpace s = VectorSpace.Dispatch(basis); + return Dispatch(s, b); + } + +} + +/// +public abstract class AffineSpace : AffineSpace + where TVector : Vector, new() + where TMatrix : Matrix, new() + where TSpace : VectorSpace, new() + where TAffine : AffineSpace, new() + where TScalar : notnull +{ + + /// + /// + public LinearSystem InverseSystem { + get + { + TSpace otc = BaseSpace.OrthogonalComplementSpace; + List bs = otc.ToList(); + while(bs.Count < Dim) + bs.Add(new TVector()); + TMatrix mA = Matrix.Dispatch(bs); + TVector rs = mA * Bias; + return new LinearSystem(mA, rs); + } + } + + /// + public override TAffine Intersection(TAffine otherSpace) + { + LinearSystem sa = InverseSystem; + LinearSystem sb = otherSpace.InverseSystem; + + ReducedRowEchelonForm brref = sb.A.RowReduce(); + TMatrix sam = sa.A; + TVector say = sa.Y; + TVector sby = brref.RowOperation * sb.Y; + TMatrix sbm = brref.ReducedForm; + for (int i = 1; i <= sbm.Dim; i++) + if (sbm[i].IsEqualApprox(sbm[i].VarZero)) + { + sam[sam.Dim] = sbm[i]; + say[sam.Dim] = sby[i]; + } + else + { + TVector x = sam.Inv() * say; + if ((sbm * x).IsEqualApprox(sby)) + break; + return new TAffine(); + } + + ReducedRowEchelonForm rref = sam.RowReduce(); + say = rref.RowOperation * say; + sam = rref.ReducedForm; + return sam.SolutionSpace(say); + } +} \ No newline at end of file diff --git a/src/Algebra/AffineSpaces/FAffineSpace.cs b/src/Algebra/AffineSpaces/FAffineSpace.cs new file mode 100644 index 0000000..cc7bc66 --- /dev/null +++ b/src/Algebra/AffineSpaces/FAffineSpace.cs @@ -0,0 +1,19 @@ +using Skeleton.Algebra.AffineSpaces; + +namespace Skeleton.Algebra; + +public partial class CategoryOf +{ + public static partial class OnField + { + /// + public class FAffineSpace : AffineSpace + { + /// + public FAffineSpace() { } + + /// + public override int Dim => FDim; + } + } +} \ No newline at end of file diff --git a/src/Algebra/CategoryOf.cs b/src/Algebra/CategoryOf.cs new file mode 100644 index 0000000..d23932c --- /dev/null +++ b/src/Algebra/CategoryOf.cs @@ -0,0 +1,29 @@ +using System.Reflection; +using Skeleton.Algebra.DimensionProviders; +using Skeleton.Algebra.ScalarFieldStructure; + +namespace Skeleton.Algebra; + +/// +/// +/// +/// +public static partial class CategoryOf +{ + private static readonly int FDim = typeof(TDimension).GetCustomAttribute()!.Dim; + + /// + /// + /// + /// + public static partial class OnField where TScalar : notnull + { + /// + /// Field Structure + /// + public static readonly FieldStructure FieldStructure = FieldStructure.Dispatch(); + + } + +} + diff --git a/src/Algebra/DimensionProviders/IDim2.cs b/src/Algebra/DimensionProviders/IDim2.cs new file mode 100644 index 0000000..098252b --- /dev/null +++ b/src/Algebra/DimensionProviders/IDim2.cs @@ -0,0 +1,6 @@ +namespace Skeleton.Algebra.DimensionProviders; +[OfDimension(2)] +public interface IDim2 +{ + +} \ No newline at end of file diff --git a/src/Algebra/DimensionProviders/IDim3.cs b/src/Algebra/DimensionProviders/IDim3.cs new file mode 100644 index 0000000..493c7f3 --- /dev/null +++ b/src/Algebra/DimensionProviders/IDim3.cs @@ -0,0 +1,6 @@ +namespace Skeleton.Algebra.DimensionProviders; +[OfDimension(3)] +public interface IDim3 +{ + +} \ No newline at end of file diff --git a/src/Algebra/DimensionProviders/IDim4.cs b/src/Algebra/DimensionProviders/IDim4.cs new file mode 100644 index 0000000..6d5fd06 --- /dev/null +++ b/src/Algebra/DimensionProviders/IDim4.cs @@ -0,0 +1,6 @@ +namespace Skeleton.Algebra.DimensionProviders; +[OfDimension(4)] +public interface IDim4 +{ + +} \ No newline at end of file diff --git a/src/Algebra/DimensionProviders/OfDimension.cs b/src/Algebra/DimensionProviders/OfDimension.cs new file mode 100644 index 0000000..3a1ad16 --- /dev/null +++ b/src/Algebra/DimensionProviders/OfDimension.cs @@ -0,0 +1,18 @@ +namespace Skeleton.Algebra.DimensionProviders; + +/// +/// provide dimension information for objects +/// +[AttributeUsage(AttributeTargets.Class | AttributeTargets.Interface)] +public class OfDimension : Attribute +{ + /// + public OfDimension(int dim) + { + Dim = dim; + } + /// + /// dimension of this class + /// + public int Dim { get; } +} diff --git a/src/Algebra/Extensions/GeneralExt.cs b/src/Algebra/Extensions/GeneralExt.cs new file mode 100644 index 0000000..f60c4d9 --- /dev/null +++ b/src/Algebra/Extensions/GeneralExt.cs @@ -0,0 +1,256 @@ +using System.Numerics; +using Skeleton.Algebra.ScalarFieldStructure; +using Skeleton.Utils.Helpers; + +namespace Skeleton.Algebra.Extensions; + +/// +/// +public static class GeneralExt +{ + + + /// + /// + private static readonly ThreadLocal> TolCache = new(() => new Stack()); + + [ThreadStatic] private static double TolValue; + + private static double Tol => TolCache.Value!.TryPeek(out TolValue) ? TolValue : 1E-8D; + + /// + /// + /// + public static void SetTol(double t) + { + TolCache.Value!.Push(t); + } + + /// + /// + public static void ResetTol() + { + TolCache.Value!.Pop(); + } + + private static bool CZero(this float x) + { + return Math.Abs(x) < Tol; + } + + private static bool CZero(this double x) + { + return Math.Abs(x) < Tol; + } + + /// + /// + /// + /// + public static bool CZero(this Complex x) + { + return x.Real.CZero() && x.Imaginary.CZero(); + } + + /// + /// + /// + /// + /// + public static bool IsEqualApprox(this float l, float r) + { + return (l - r).CZero(); + } + + /// + /// + /// + /// + /// + /// + public static bool IsSignificantlyGreaterThen(this double l, double r) + { + if (l.IsEqualApprox(r)) + return false; + return l > r; + } + + /// + /// + /// + /// + /// + /// + public static bool IsSignificantlyLessThen(this double l, double r) + { + if (l.IsEqualApprox(r)) + return false; + return l < r; + } + + /// + /// + /// + /// + /// + public static bool IsEqualApprox(this double l, double r) + { + return (l - r).CZero(); + } + + /// + /// + /// + /// + /// + public static bool IsEqualApprox(this Complex l, Complex r) + => l.Real.IsEqualApprox(r.Real) && l.Imaginary.IsEqualApprox(r.Imaginary); + + + /// + /// + /// + /// + public static Complex LSum(this IEnumerable l) + { + Complex res = Complex.Zero; + foreach (Complex x in l) + res += x; + return res; + } + + /// + /// + /// + /// + public static float LSum(this IEnumerable l) + { + float res = 0; + foreach (float x in l) + res += x; + return res; + } + + /// + /// + /// + /// + public static double LSum(this IEnumerable l) + { + double res = 0; + foreach (double x in l) + res += x; + return res; + } + + /// + /// + /// + /// + public static Complex LProd(this IEnumerable l) + { + Complex res = Complex.One; + foreach (Complex x in l) + res *= x; + return res; + } + + /// + /// + /// + /// + /// + public static bool ApproxContains(this IEnumerable l, Complex k) + { + foreach (Complex x in l) + if (x.IsEqualApprox(k)) + return true; + return false; + } + + /// + /// + /// + /// + /// + public static bool ApproxContains(this IEnumerable l, float k) + { + foreach (float x in l) + if (x.IsEqualApprox(k)) + return true; + return false; + } + + /// + /// + /// + /// + /// + public static bool ApproxContains(this IEnumerable l, double k) + { + foreach (double x in l) + if (x.IsEqualApprox(k)) + return true; + return false; + } + + /// + /// + /// + /// + public static Dictionary RoughGroup(this IEnumerable r) + { + List> clusters = new List>(); + foreach (Complex c in r) + { + bool added = false; + foreach (List cluster in clusters) + { + double distance = cluster.Min(x => (c - x).ManhattanError()); + if (distance < 1E-6d) + { + cluster.Add(c); + added = true; + } + } + + if (!added) + clusters.Add(new List { c }); + } + + Dictionary res = new Dictionary(); + foreach (List cluster in clusters) + { + Complex center = cluster.Aggregate((a, b) => a + b) / cluster.Count; + res[center] = cluster.Count; + } + + return res; + } + + /// + /// group complex numbers by count + /// + /// + /// + public static Dictionary Group(this IEnumerable l) + where TScalar : notnull + { + Dictionary res = new (); + FieldStructure structure = FieldStructure.Dispatch(); + foreach (TScalar val in l) + { + bool exists = false; + foreach (TScalar u in res.Keys) + if (structure.IsEqualApprox(val, u)) + { + exists = true; + res[u] += 1; + break; + } + if (!exists) + res[val] = 1; + } + return res; + } +} diff --git a/src/Algebra/FixableObjects/FixableTensor.cs b/src/Algebra/FixableObjects/FixableTensor.cs new file mode 100644 index 0000000..67e8bd6 --- /dev/null +++ b/src/Algebra/FixableObjects/FixableTensor.cs @@ -0,0 +1,50 @@ +namespace Skeleton.Algebra.FixableObjects; + +/// +/// error in tensor might increase as the increased number of calculations +/// provide methods to fix this +/// +public abstract class FixableTensor +{ + /// + /// Dimension of the tensor + /// + public abstract int Dim { get; } + + /// + /// + protected FixableTensor() { } + + /// + /// if the tensor need to be fixed + /// + internal virtual bool NeedFix => false; + + /// + /// fix the tensor + /// + internal virtual void Fix() + { + + } + + /// + /// Call abstract/virtual methods in constructor to kill warning + /// most derived class's implementation will be called + /// + /// + /// + /// + protected static T CancelAbstract(T a) => a; + + + + internal TFixable FixByPass() + where TFixable : FixableTensor + { + Fix(); + return (this as TFixable)!; + } + +} + diff --git a/src/Algebra/FixableObjects/PostConstructionFixExtension.cs b/src/Algebra/FixableObjects/PostConstructionFixExtension.cs new file mode 100644 index 0000000..0d42478 --- /dev/null +++ b/src/Algebra/FixableObjects/PostConstructionFixExtension.cs @@ -0,0 +1,11 @@ +namespace Skeleton.Algebra.FixableObjects; + +internal static class PostConstructionFixExtension +{ + internal static void PostConstructionFix(this TFixable obj) + where TFixable : FixableTensor + { + if (obj.GetType() == typeof(TFixable) && obj.NeedFix) + obj.Fix(); + } +} \ No newline at end of file diff --git a/src/Algebra/Groups/FUnitaryGroup.cs b/src/Algebra/Groups/FUnitaryGroup.cs new file mode 100644 index 0000000..6bd6298 --- /dev/null +++ b/src/Algebra/Groups/FUnitaryGroup.cs @@ -0,0 +1,21 @@ +using Skeleton.Algebra.Groups; + +namespace Skeleton.Algebra; + +public static partial class CategoryOf +{ + /// + /// unitary group + /// + public class FUnitaryGroup : Group + { + /// + public override FUnitaryMatrix GroupUnit => FUnitaryMatrix.One; + + /// + public override FUnitaryMatrix GroupInv(FUnitaryMatrix element) => element.Inv(); + + /// + public override FUnitaryMatrix GroupOperation(FUnitaryMatrix a, FUnitaryMatrix b) => a * b; + } +} \ No newline at end of file diff --git a/src/Algebra/Groups/Group.cs b/src/Algebra/Groups/Group.cs new file mode 100644 index 0000000..22aa9bf --- /dev/null +++ b/src/Algebra/Groups/Group.cs @@ -0,0 +1,46 @@ +namespace Skeleton.Algebra.Groups; + +/// +/// a set of elements closed under group operation +/// +/// +public abstract class Group +{ + /// + /// unit of the group
+ /// unit * x = x * unit = x for all x in group + /// + public abstract TElement GroupUnit { get; } + /// + /// inverse of the element under group operation + /// + /// + /// + public abstract TElement GroupInv(TElement element); + /// + /// + /// + /// + /// + /// + public abstract TElement GroupOperation(TElement a, TElement b); + + /// + /// Commutator of group + /// + /// a + /// b + /// + public TElement GroupCommutator(TElement a, TElement b) => + GroupOperation(GroupOperation(a, b), GroupOperation(GroupInv(a), GroupInv(b))); + /// + /// a conjugate of b => aba^{-1} + /// + /// a + /// b + /// a b a^{-1} + public TElement Conjugation(TElement a, TElement b) => + GroupOperation(a, GroupOperation(b, GroupInv(a))); + + +} diff --git a/src/Algebra/Matrices/FMatrix.cs b/src/Algebra/Matrices/FMatrix.cs new file mode 100644 index 0000000..8a4a885 --- /dev/null +++ b/src/Algebra/Matrices/FMatrix.cs @@ -0,0 +1,206 @@ +using System.Numerics; +using Skeleton.Algebra.AdditionalProperties.Extensions; +using Skeleton.Algebra.DimensionProviders; +using Skeleton.Algebra.FixableObjects; +using Skeleton.Algebra.Matrices; +using Skeleton.DataStructure.Packs.MatrixDecompositions; + +namespace Skeleton.Algebra; + +public partial class CategoryOf +{ + public static partial class OnField + { + /// + public class FMatrix : Matrix + { + /// + public FMatrix() { } + + /// + public FMatrix(params TScalar[] es) : base(es) => this.PostConstructionFix(); + + /// + public FMatrix(IEnumerable es, bool row = true) : base(es, row) => + this.PostConstructionFix(); + + /// + public override int Dim => FDim; + + /// + public override int Rank() => Im().Rank; + + /// + public override JDPair EigenDecomposition() => + InternalEigenDecomposition(); + + /// + /// real part of the matrix
+ /// homomorphism from field to complex has to be defined + /// + public OnField.FMatrix Real => + new(this.Select(row => row.Real)); + + /// + /// imaginary part of the matrix
+ /// homomorphism from field to complex has to be defined + /// + public OnField.FMatrix Imaginary => + new(this.Select(row => row.Imaginary)); + + /// + /// map the matrix to a complex matrix
+ /// homomorphism from field to complex has to be defined + /// + /// + public OnField.FMatrix AsComplex() => + new(this.Select(x => x.AsComplex())); + + /// + /// cast from matrix over K to a matrix over R
+ /// morphism K -> R has to be defined + /// + /// + /// + public static explicit operator OnField.FMatrix(FMatrix a) => + new(a.Select(x => (OnField.FVector)x)); + + /// + /// cast from matrix over K to a matrix over C
+ /// morphism K -> C has to be defined + /// + /// + /// + public static explicit operator OnField.FMatrix(FMatrix a) => + new(a.Select(x => (OnField.FVector)x)); + + /// + /// cast from matrix over R to a matrix over K
+ /// morphism R -> K has to be defined + /// + /// + /// + public static explicit operator FMatrix(OnField.FMatrix a) => + new(a.Select(x => (FVector)x)); + + /// + /// cast from matrix over C to a matrix over K
+ /// morphism C -> K has to be defined + /// + /// + /// + public static explicit operator FMatrix(OnField.FMatrix a) => + new(a.Select(x => (FVector)x)); + + /// + /// kernel space + /// + public FVectorSpace Ker => Ker(); + + /// + /// average of column vectors + /// + public FVector ColumnAverage + { + get + { + FVector res = new(); + for (int i = 1; i <= FDim; i++) + res += this[i, false]; + return res / FDim; + } + } + + /// + /// + /// + public FVector ColumnAverageBivariant + { + get + { + double t = 0; + FVector res = new(); + for (int i = 1; i <= FDim; i++) + { + double p = Math.Abs(Dim * 0.5 - i); + res += this[i, false] * p; + t += p; + } + return res / t; + } + } + + internal sealed override void ResetOne() + { + for (int i = 1; i <= FDim; i++) + { + this[i, i] = FieldStructure.MultiplicationUnit; + for (int j = 1; j < i; j++) + { + this[i, j] = FieldStructure.AdditionUnit; + this[j, i] = FieldStructure.AdditionUnit; + } + } + } + + internal sealed override void ResetZero() + { + for (int i = 1; i <= FDim; i++) + { + this[i, i] = FieldStructure.AdditionUnit; + for (int j = 1; j < i; j++) + { + this[i, j] = FieldStructure.AdditionUnit; + this[j, i] = FieldStructure.AdditionUnit; + } + } + } + + internal override void SelfMul(FMatrix other) + { + if (this is CategoryOf.OnField.FMatrix c33) + { + ComplexMatrixExtension.FastSelfMul(c33, (other as CategoryOf.OnField.FMatrix)!); + return; + } + + TScalar[] r1 = new TScalar[FDim]; + TScalar[] r2 = new TScalar[FDim]; + for (int r = 1; r <= FDim + 1; r++) + { + for (int c = 1; c <= FDim; c++) + { + if (r <= FDim) + r1[c-1] = this[r] * other[c, false]; + if (r > 1) + this[r - 1, c] = r2[c-1]; + } + (r1, r2) = (r2, r1); + } + } + + internal override void SelfLMul(FMatrix other) + { + if (this is CategoryOf.OnField.FMatrix c33) + { + ComplexMatrixExtension.FastSelfLMul((other as CategoryOf.OnField.FMatrix)!, c33); + return; + } + TScalar[] c1 = new TScalar[FDim]; + TScalar[] c2 = new TScalar[FDim]; + for (int c = 1; c <= FDim + 1; c++) + { + for (int r = 1; r <= FDim; r++) + { + if (c <= FDim) + c1[r - 1] = this[c, false] * other[r]; + if (c > 1) + this[r, c - 1] = c2[r - 1]; + } + (c1, c2) = (c2, c1); + } + } + } + + } +} \ No newline at end of file diff --git a/src/Algebra/Matrices/Matrix.cs b/src/Algebra/Matrices/Matrix.cs new file mode 100644 index 0000000..0bebf8d --- /dev/null +++ b/src/Algebra/Matrices/Matrix.cs @@ -0,0 +1,1312 @@ +using System.Collections; +using System.Numerics; +using Skeleton.Algebra.AdditionalProperties.Extensions; +using Skeleton.Algebra.AffineSpaces; +using Skeleton.Algebra.DimensionProviders; +using Skeleton.Algebra.Extensions; +using Skeleton.Algebra.FixableObjects; +using Skeleton.Algebra.Matrices.MatrixProperty; +using Skeleton.Algebra.ScalarFieldStructure; +using Skeleton.Algebra.Vectors; +using Skeleton.Algebra.VectorSpaces; +using Skeleton.DataStructure.Packs; +using Skeleton.DataStructure.Packs.MatrixDecompositions; + +namespace Skeleton.Algebra.Matrices; + +/// +/// matrices +/// +public abstract class Matrix : FixableTensor +{ + /// + /// string representation of this matrix + /// + public abstract string Representation { get; } + + /// + /// + /// + public abstract string PythonRepresentation { get; } + + /// + /// + /// + public abstract string CSharpRepresentation { get; } + + /// + /// + /// + public abstract string LatexRepresentation { get; } + + /// + /// dimension of null space + /// + /// + public int Nullity() => Dim - Rank(); + + /// + /// dimension of image space + /// + /// + public abstract int Rank(); +} + +/// +/// meta matrices +/// +/// +public abstract class Matrix< TMatrix> : Matrix + where TMatrix : Matrix< TMatrix>, new() +{ + /// + /// zero matrix as variable + /// + public TMatrix VarZero => new(); + + /// + /// id matrix as variable + /// + internal abstract TMatrix VarOne { get; } + + /// + /// process row reduction to the matrix + /// + /// + public abstract ReducedRowEchelonForm RowReduce(); + + /// + /// inverse of the matrix + /// + /// + public TMatrix Inv() => RowReduce().RowOperation; + + internal abstract TMatrix Addition(TMatrix other); + internal abstract TMatrix Subtraction(TMatrix other); + internal abstract TMatrix Negation(); + internal abstract TMatrix Mul(TMatrix other); + internal abstract void SelfMul(TMatrix other); + internal abstract void SelfLMul(TMatrix other); + /// + /// Add two matrix of the same dim over same field + /// + /// + /// + /// + public static TMatrix operator +(Matrix self, TMatrix other) => self.Addition(other); + + /// + /// negation of the matrix + /// + /// + /// + public static TMatrix operator -(Matrix< TMatrix> self) => self.Negation(); + + /// + /// subtract on matrix from another + /// + /// left operand + /// right operand + /// + public static TMatrix operator -(Matrix< TMatrix> self, TMatrix other) => self.Subtraction(other); + /// + /// multiply two matrix of the same dimension and field + /// + /// left operand + /// right operand + /// + public static TMatrix operator *(Matrix< TMatrix> self, TMatrix other) => self.Mul(other); + + /// + /// divide one matrix by another + /// + /// left operand + /// right operand + /// + public static TMatrix operator /(Matrix< TMatrix> self, TMatrix other) => self * other.Inv(); + + /// + /// copy of matrix + /// + /// + public abstract TMatrix Copy(); + internal abstract void ResetZero(); + internal abstract void ResetOne(); + + /// + /// commutator of matrices [a, b] = ab - ba + /// + /// right operand + /// left operand + /// + public static TMatrix Commutator(TMatrix a, TMatrix b) => a * b - b * a; + + /// + /// anti commutator of matrices {a, b} = ab + ba + /// + /// + /// + /// + public static TMatrix AntiCommutator(TMatrix a, TMatrix b) => a * b + b * a; + + /// + /// Transpose of the matrix. Mij = Nji + /// + /// + public abstract TMatrix Transpose(); + + /// + /// transpose current matrix, modify on current matrix without return a new one + /// + public abstract void SelfTranspose(); + + /// + /// Eigen decomposition of the matrix, M = JD inv(J)
+ /// will throw exception if the field is not complete + /// + /// + public abstract JDPair EigenDecomposition(); +} + +/// +public abstract class Matrix< TVector, TMatrix> : Matrix< TMatrix>, IEnumerable + where TVector : Vectors.Vector< TVector>, new() + where TMatrix : Matrix< TVector, TMatrix>, new() +{ + /// + protected Matrix() => Elements = GetNewArray(); + + internal TVector[] GetNewArray() + { + TVector[] res = new TVector[Dim]; + for (int i = 0; i < Dim; i++) + res[i] = new TVector(); + return res; + } + + /// + /// Construct by vectors + /// + /// collection of vectors + /// true if given vectors are row + protected Matrix(IEnumerable rows, bool byRow = true) + { + Elements = rows.ToArray(); + if (!byRow) + CancelAbstract(SelfTranspose)(); + } + internal TVector[] Elements { get; set; } + + /// + /// get/set row/column vector + /// + /// row/column number + /// if true, set, get by row. otherwise by column + public abstract TVector this[int i, bool row = true] { get; set; } + + /// + /// matrix as collection of rows + /// + /// + public IEnumerator GetEnumerator() + { + for (int i = 0; i < Dim; i++) + yield return Elements[i]; + } + IEnumerator IEnumerable.GetEnumerator() => GetEnumerator(); + internal override void ResetZero() + { + Elements = GetNewArray(); + } + + internal override void Fix() + { + base.Fix(); + for (int i = 0; i < Dim; i++) + Elements[i].Fix(); + } + + /// + /// determine if the matrix is close to another matrix + /// + /// + /// + public bool IsEqualApprox(TMatrix other) + { + for(int i = 1; i<=Dim;i++) + if(!this[i].IsEqualApprox(other[i])) + return false; + return true; + } + /// + /// Left Multiplication + /// + /// + /// + internal abstract TVector LMul(TVector other); + + /// + /// left multiplication + /// + /// matrix, left operand + /// vector, right operand + /// + public static TVector operator *(Matrix< TVector, TMatrix> self, TVector other) + => self.LMul(other); + + /// + /// right multiplication + /// + /// + /// + /// + public static TVector operator *(TVector other, Matrix< TVector, TMatrix> self) => self.Transpose() * other; + + /// + /// kernel basis of the image space + /// + /// + public abstract IEnumerable KerBasis(); + + /// + /// Conj of kernel basis + /// + /// + public IEnumerable ConjKerBasis() => KerBasis().Select(x => x.Conj()); + /// + /// + /// + /// + /// + public TRes MatrixCast() + where TRes : TMatrix, new() + { + TRes res = new TRes(); + res.Elements = Elements; + res.Fix(); + return res; + } + internal static TMatrix Dispatch(IEnumerable es, bool row = true) + { + TMatrix res = new () + { + Elements = es.ToArray() + }; + return res; + } + + /// + /// forget additional structures + /// + public TMatrix AsMatrix => (this as TMatrix)!; +} + +/// +/// General idea of square matrix(2x2 to 4x4) +/// +/// Type, double or complex +/// compatible vector +/// the matrix +public abstract class Matrix : Matrix + where TVector : Vector, new() + where TMatrix : Matrix, new() + where TScalar : notnull +{ + /// + /// id matrix + /// + public static readonly TMatrix One = new TMatrix().VarOne; + + /// + /// zero matrix + /// + public static readonly TMatrix Zero = new(); + + /// + protected Matrix() => + Property = new MatrixProperty((this as TMatrix)!); + + /// + protected Matrix(params TScalar[] es) + { + Property = new MatrixProperty((this as TMatrix)!); + ResetZero(); + for (int i = 1; i <= Dim; i++) + for (int j = 1; j <= Dim; j++) + this[j, i] = es[i - 1 + (j - 1) * Dim]; + } + + /// + /// properties of matrix + /// + public MatrixProperty Property { get; } + + /// + protected Matrix(IEnumerable rows, bool byRow = true) : base(rows, byRow) => + Property = new MatrixProperty((this as TMatrix)!); + + /// + protected Matrix(IEnumerable elements) + { + Elements = new TVector[Dim]; + TScalar[] ac = elements.ToArray(); + for (int i = 1; i <= Dim; i++) + { + TVector rowI = new (); + for (int j = 1; j <= Dim; j++) + rowI[j] = ac[(i - 1) * Dim + j - 1]; + Elements[i - 1] = rowI; + } + Property = new MatrixProperty((this as TMatrix)!); + } + + /// + /// + public double MaxError => this.Select(v => v.MaxError).Max(); + + /// + /// get/set element by row and column + /// + /// row number + /// column number + public virtual TScalar this[int r, int c] + { + get => Elements[r - 1][c]; + set => Elements[r - 1][c] = value; + } + + /// + /// get/set row or column vector + /// + /// row/column number + /// if true, set, get by row. otherwise by column + public override TVector this[int i, bool row = true] + { + get + { + if (row) + return Elements[i - 1]; + TVector res = new TVector(); + for (int j = 1; j <= Dim; j++) + res[j] = this[j, i]; + return res; + } + set + { + if (row) + for (int j = 1; j <= Dim; j++) + this[i, j] = value[j]; + else + for (int j = 1; j <= Dim; j++) + this[j, i] = value[j]; + } + } + + + + /// + /// elements on the diagonal + /// + public IEnumerable Diagonals + { + get + { + for (int i = 1; i <= Dim; i++) + yield return this[i, i]; + } + } + + /// + /// + public IEnumerable OffDiagonals + { + get + { + for (int i = 2; i <= Dim; i++) + for (int j = 1; j < i; j++) + { + yield return this[i, j]; + yield return this[j, i]; + } + } + } + + internal override TMatrix VarOne + { + get + { + TMatrix res = VarZero; + res.ResetOne(); + return res; + } + } + + /// + public override void SelfTranspose() + { + for (int i = 1; i <= Dim; i++) + for (int j = 1; j < i; j++) + (this[i, j], this[j, i]) = (this[j, i], this[i, j]); + } + + /// + public override TMatrix Transpose() + { + TMatrix res = VarZero; + for (int i = 1; i <= Dim; i++) + { + res[i, i] = this[i, i]; + for (int j = 1; j < i; j++) + (res[i, j], res[j, i]) = (this[j, i], this[i, j]); + } + + return res; + } + + /// + /// Det of the matrix + /// + /// + public virtual TScalar Det() => PLUDecomposition().Det; + + /// + /// trace of the matrix(sum of diagonal) + /// + /// + public virtual TScalar Trace() + { + TScalar res = StaticAccess.AdditionUnit; + for (int i = 1; i <= Dim; i++) + res = Structure.Addition(res, this[i, i]); + return res; + } + + internal override TVector LMul(TVector other) + { + TVector v = new TVector(); + for (int i = 1; i <= Dim; i++) + v[i] = this[i] * other; + return v; + } + + /// + public override IEnumerable KerBasis() + { + ReducedRowEchelonForm rref = RowReduce(); + TMatrix e = rref.ReducedForm; + TMatrix rbs = e.VarZero; + int row = 1; + for (int col = 1; col <= Dim; col++) + { + if (rref.Pivots.Contains(col)) + { + for (int c = 1; c <= Dim; c++) + if (c != col) + rbs[col, c] = Structure.AdditionInverse(e[row, c]); + row += 1; + } + + if (!rref.Pivots.Contains(col)) + rbs[col, col] = Structure.MultiplicationUnit; + } + + return rbs.Transpose().Where(c => !c.IsEqualApprox(c.VarZero)).ToArray(); + } + + + /// + /// Kernel space of the matrix + /// v | M v = 0 + /// + /// + /// + /// + public TVectorSpace Ker() + where TVectorSpace : VectorSpace, new() + { + TVectorSpace res = new(); + res.Fill(KerBasis()); + return res; + } + + /// + /// + /// + /// + /// + public TVectorSpace ConjKer() + where TVectorSpace : VectorSpace, new() + { + TVectorSpace res = new TVectorSpace(); + res.Fill(ConjKerBasis()); + return res; + } + + + /// + /// Image space of the matrix + /// vector space span by col vectors + /// Im(M) = otc of Ker(M*) + /// + /// + /// + /// + public TVectorSpace Im() + where TVectorSpace : VectorSpace, new() + => VectorSpace.Dispatch(Transpose()); + + /// + /// solution space to the system Mx = y + /// + /// + /// + /// + /// + /// + public TAffineSpace SolutionSpace(TVector y) + where TVectorSpace : VectorSpace, new() + where TAffineSpace : AffineSpace, new() + { + ReducedRowEchelonForm rref = RowReduce(); + TAffineSpace res = new TAffineSpace(); + TVector ty = rref.RowOperation * y; + TVector ry = y.VarZero; + TVectorSpace s = Ker(); + if (ty.EndingZeros < s.Rank) + throw new Exception("No Solution"); // TODO + foreach (int pivot in rref.Pivots) + for (int r = 1; r <= Dim; r++) + if (Structure.IsEqualApprox(rref.ReducedForm[r, pivot], Structure.MultiplicationUnit)) + { + ry[pivot] = ty[r]; + break; + } + res.Fill(s, ry); + return res; + } + + /// + /// Scalar Multiplication + /// + /// + /// + internal virtual TMatrix ScalarMul(TScalar other) + { + TMatrix res = new (); + for (int i = 1; i <= Dim; i++) + for (int j = 1; j <= Dim; j++) + res[i, j] = Structure.Multiplication(this[i, j], other); + return res; + } + + internal virtual TMatrix ScalarDiv(TScalar other) => + ScalarMul(Structure.MultiplicationInverse(other)); + + internal override TMatrix Addition(TMatrix other) + { + TMatrix res = new TMatrix(); + for (int i = 1; i <= Dim; i++) + for (int j = 1; j <= Dim; j++) + res[i, j] = Structure.Addition(this[i, j], other[i, j]); + return res; + } + + internal override TMatrix Negation() + { + TMatrix res = new(); + for (int i = 1; i <= Dim; i++) + for (int j = 1; j <= Dim; j++) + res[i, j] = Structure.AdditionInverse(this[i, j]); + return res; + } + + internal override TMatrix Subtraction(TMatrix other) + { + TMatrix res = new TMatrix(); + for (int i = 1; i <= Dim; i++) + for (int j = 1; j <= Dim; j++) + res[i, j] = Structure.Subtraction(this[i, j], other[i, j]); + return res; + } + + + /// + /// scalar multiplication + /// + /// + /// + /// + public static TMatrix operator *(Matrix self, TScalar other) + => self.ScalarMul(other); + + + /// + /// scalar multiplication + /// + /// + /// + /// + public static TMatrix operator *(TScalar other, Matrix self) + => self.ScalarMul(other); + + /// + /// + /// + /// + /// + /// + public static TMatrix operator *(Matrix self, double other) => + self.ScalarMul(Structure.FromReal(other)); + + /// + /// + /// + /// + /// + /// + public static TMatrix operator *(double other, Matrix self) => + self * other; + /// + /// scalar division + /// + /// + /// + /// + public static TMatrix operator /(Matrix self, TScalar other) + => self.ScalarDiv(other); + + internal TComplexMatrix ComplexCast() + where TComplexVector : Vector, new() + where TComplexMatrix : Matrix, new() + => Matrix + .Dispatch(Elements.Select(x => x.ComplexCast())); + + + internal IEnumerable> EigenInfo() + where TVectorSpace : VectorSpace, new() + { + Dictionary eigenValues = EigenValues().Group(); + foreach (TScalar eigenValue in eigenValues.Keys) + { + TMatrix h = this - eigenValue * VarOne; + TMatrix hN = h.Copy(); + int algebraicMultiplicity = eigenValues[eigenValue]; + List rs = new(); + int rx = -1; + List oS = new(); + while (rx < algebraicMultiplicity) + { + oS.Add(hN); + rx = hN.Nullity(); + hN = hN * h; + rs.Add(rx); + } + + List ss = new List { rs[0] }; + for (int idx = 1; idx < rs.Count; idx++) + ss.Add(rs[idx] - rs[idx - 1]); + List ms = new(); + for (int idx = 1; idx < ss.Count; idx++) + ms.Add(ss[idx - 1] - ss[idx]); + ms.Add(ss.Last()); + List geometricMultiplicities = new(); + for (int ic = ms.Count - 1; ic >= 0; ic--) + for (int ix = 0; ix < ms[ic]; ix++) + geometricMultiplicities.Add(ic + 1); + TVectorSpace w = new(); + foreach (int geometricMultiplicity in geometricMultiplicities) + { + List generalizedEigenVectors = new(); + foreach (TVector eigenVector in oS[geometricMultiplicity - 1].Ker()) + { + if (eigenVector.IsEqualApprox(eigenVector.VarZero)) + break; + if ( + !w.ContainsVector(eigenVector) + && + ( + geometricMultiplicity <= 1 + || + !oS[geometricMultiplicity - 2].Ker().ContainsVector(eigenVector) + ) + ) + { + generalizedEigenVectors.Add(eigenVector); + w.AddBasis(eigenVector); + for (int i = geometricMultiplicity - 2; i < oS.Count - 1 && i >= 0; i++) + { + TVector vh = oS[i] * eigenVector; + if (vh.IsEqualApprox(vh.VarZero)) + break; + generalizedEigenVectors.Add(vh); + w.AddBasis(vh); + } + + break; + } + } + + yield return new EigenPair(eigenValue, generalizedEigenVectors.ToArray()); + } + } + } + + private EigenPair[] Rebalance(EigenPair[] a) + { + List> res = a.ToList(); + double valueTolerance = 1E-5d; + double vectorTolerance = 1E-3d; + while (res.Sum(p => p.EigenVectors.Length) > Dim) + { + int min1 = res.Count + 1; + int min2 = res.Count + 1; + double minValueErr = valueTolerance; + double minVectorErr = vectorTolerance; + for (int i = 0; i < res.Count; i++) + { + for (int j = 0; j < i; j++) + { + if (res[i].EigenVectors.Length != res[j].EigenVectors.Length) + continue; + double valueErr = Structure.MaxError(Structure.Subtraction(res[i].EigenValue, res[j].EigenValue)); + double vectorError = res[i].EigenVectors + .Zip(res[j].EigenVectors) + .Select(pair => (pair.First - pair.Second).MaxError) + .Max(); + if ( + (vectorError < minVectorErr && valueErr < minValueErr + valueTolerance) || + (vectorError < minValueErr + vectorTolerance && valueErr < minValueErr) + ) + { + min1 = i; + min2 = j; + minValueErr = valueErr; + minVectorErr = vectorError; + } + } + } + + TScalar val2 = Structure.Addition(Structure.MultiplicationUnit, Structure.MultiplicationUnit); + res[min1].EigenValue = + Structure.Division(Structure.Addition(res[min1].EigenValue, res[min2].EigenValue), val2); + for (int i = 0; i < res[min1].EigenVectors.Length; i++) + res[min1].EigenVectors[i] = (res[min1].EigenVectors[i] + res[min2].EigenVectors[i]) / 2d; + res.RemoveAt(min2); + } + + return res.ToArray(); + } + + internal JDPair InternalEigenDecomposition() + where TVectorSpace : VectorSpace, new() + { + if (Property.IsDiagonal) + return new JDPair( + VarOne, + this.Copy() + ); + TMatrix p = VarZero; + TMatrix jordan = VarZero; + int h = 1; + EigenPair[] bls = EigenInfo().ToArray(); + if (bls.Length > Dim) + bls = Rebalance(bls); + foreach (EigenPair block in bls) + { + int geoMultiplicity = block.EigenVectors.Length; + for (int idx = geoMultiplicity - 1; idx >= 0; idx--) + { + jordan[h, h] = block.EigenValue; + if (idx != geoMultiplicity - 1) + jordan[h - 1, h] = Structure.MultiplicationUnit; + p[h, false] = block.EigenVectors[idx].Normalize; + h += 1; + } + } + + return new JDPair(p, jordan); + } + + /// + /// + /// + /// + /// + /// + /// + /// + public TComplexMatrix Pow(Complex o) + { + throw new Exception("Should be overloaded"); + } + + /// + /// swap two rows + /// + /// + /// + public void SwapRow(int r1, int r2) + { + for (int i = 1; i <= Dim; i++) + (this[r1, i], this[r2, i]) = (this[r2, i], this[r1, i]); + } + + /// + /// swap two columns + /// + /// + /// + public void SwapColumn(int c1, int c2) + { + for (int i = 1; i <= Dim; i++) + (this[i, c1], this[i, c2]) = (this[i, c2], this[i, c1]); + } + + internal void Fill(TScalar[] es) + { + for (int i = 1; i <= Dim; i++) + for (int j = 1; j <= Dim; j++) + this[i, j] = es[j - 1 + (i - 1) * Dim]; + + } + + internal override TMatrix Mul(TMatrix other) + { + if (this is CategoryOf.OnField.FMatrix self2) + return (ComplexMatrixExtension.FastMul( + self2, + (other as CategoryOf.OnField.FMatrix)! + ) as TMatrix)!; + if (this is CategoryOf.OnField.FMatrix self3) + + return (ComplexMatrixExtension.FastMul( + self3, + (other as CategoryOf.OnField.FMatrix)! + ) as TMatrix)!; + TMatrix res = new(); + for (int r = 1; r <= Dim; r++) + for (int c = 1; c <= Dim; c++) + res[r, c] = this[r] * other[c, false]; + + return res; + } + + /// + public sealed override ReducedRowEchelonForm RowReduce() + { + TMatrix res = Copy(); + TMatrix conj = VarOne; + List pivotCol = new List(); + bool ordered = false; + while (!ordered) + { + ordered = true; + for (int r = 1; r < Dim; r++) + if (res[r].LeadingZeros > res[r + 1].LeadingZeros) + { + ordered = false; + res.SwapRow(r, r + 1); + conj.SwapRow(r, r + 1); + } + } + + for (int r = 1; r <= Dim; r++) + { + int p = 1; + while (p < Dim && Structure.IsEqualApprox(res[r, p], StaticAccess.AdditionUnit)) + p += 1; + if (!Structure.IsEqualApprox(res[r, p], StaticAccess.AdditionUnit)) + { + pivotCol.Add(p); + conj[r] = conj[r] / res[r, p]; + res[r] = res[r] / res[r, p]; + for (int pr = 1; pr <= Dim; pr++) + { + if (r == pr) + continue; + conj[pr] = conj[pr] - conj[r] * res[pr, p]; + res[pr] = res[pr] - res[r] * res[pr, p]; + } + } + } + + ordered = false; + while (!ordered) + { + ordered = true; + for (int r = 1; r < Dim; r++) + if (res[r].LeadingZeros > res[r + 1].LeadingZeros) + { + ordered = false; + res.SwapRow(r, r + 1); + conj.SwapRow(r, r + 1); + } + } + + return new ReducedRowEchelonForm(res, conj, pivotCol.ToArray()); + } + + /// + /// Conjugate of the matrix + /// + /// + public TMatrix Conj() + { + TMatrix res = VarZero; + for (int i = 1; i <= Dim; i++) + for (int j = 1; j <= Dim; j++) + res[i, j] = Structure.Conj(this[i, j]); + return res; + } + + /// + /// the hermitian(transpose + conjugate) of a matrix + /// + /// the hermitian of the matrix + public TMatrix Hermitian() + { + TMatrix res = VarZero; + for (int i = 1; i <= Dim; i++) + for (int j = 1; j <= Dim; j++) + res[i, j] = Structure.Conj(this[j, i]); + return res; + } + + /// + public sealed override TMatrix Copy() + { + TMatrix res = new(); + res.Elements = Elements.Select(x => x.Copy()).ToArray(); + return res; + } + + /// + /// + /// + public IEnumerable FlatElements + { + get + { + for (int i = 1; i <= Dim; i++) + for (int j = 1; j <= Dim; j++) + yield return this[i, j]; + } + } + + /// + public override string Representation => string.Join("\n", Elements.Select(x => x.Representation)); + + /// + public sealed override string PythonRepresentation => "np.matrix([" + string.Join + ( + ",", this.Select(v => + "[" + string.Join(',', v.Select(x => Structure.PythonRepresentation(x))) + "]" + ) + ) + "])"; + + /// + public override string CSharpRepresentation + => $"new {GetType().Name} (" + + string.Join(",", + FlatElements.Select(x => Structure.CSharpRepresentation(x)) + ) + ")"; + + /// + public override string LatexRepresentation + => "\\begin{pmatrix}" + + string.Join("\\\\", + this.Select(v => + string.Join('&', v.Select(x => Structure.Representation(x))))) + + "\\end{pmatrix}"; + + /// + /// save matrix to string + /// + /// + public string DumpString + => "" + + string.Join("", Elements.Select(x => x.DumpString)) + + ""; + + /// + /// load matrix from string + /// + /// + /// + public static TMatrix Resolve(string dumpString) + { + TMatrix res = new(); + res.Elements = dumpString + .Replace("", "") + .Replace("", "") + .Split("") + .Select(Vector.Resolve) + .ToArray(); + return res; + } + + /// + /// M = PLU
+ /// P permutation matrix
+ /// L lower triangle matrix
+ /// U upper triangle matrix
+ /// + /// + public PLUPack PLUDecomposition() + { + TMatrix p = VarOne; + TMatrix u = Copy(); + TMatrix l = VarOne; + TScalar det = Structure.MultiplicationUnit; + + void ESwapRow(int i, int j) + { + p.SwapColumn(i, j); + l.SwapRow(i, j); + l.SwapColumn(i, j); + u.SwapRow(i, j); + det = Structure.AdditionInverse(det); + } + + void ScaleToUnit(int i, int j) + { + TScalar scale = u[i, j]; + if (Structure.IsEqualApprox(scale, StaticAccess.AdditionUnit)) + return; + u[i] /= scale; + l[i, false] *= scale; + det = Structure.Multiplication(det, scale); + } + + void Elimination(int i) + { + bool pivotNotFound = Structure.IsEqualApprox( + u[i, i], + StaticAccess.AdditionUnit + ); + for (int j = i; j <= Dim; j++) + { + ScaleToUnit(j, i); + if ( + pivotNotFound && + !Structure.IsEqualApprox( + u[j, i], + StaticAccess.AdditionUnit + ) + ) + { + ESwapRow(i, j); + pivotNotFound = false; + } + } + + if (pivotNotFound) + det = StaticAccess.AdditionUnit; + for (int j = i + 1; j <= Dim; j++) + { + if (Structure.IsEqualApprox(u[j, i], StaticAccess.AdditionUnit)) + continue; + u[j] -= u[i]; + l[i, false] += l[j, false]; + } + } + + for (int i = 1; i <= Dim; i++) + Elimination(i); + return new PLUPack(p, l, u, det); + } + + /// + /// Decomposition the matrix into product of unitary and upper triangle matrices + /// + /// + public QRPair QRDecomposition() + { + TMatrix q = VarOne; + TMatrix qia = Copy(); + TMatrix id = VarOne; + for (int i = 1; i <= Dim; i++) + { + TVector ei = id[i]; + TVector x = qia[i, false].UpTruncateBy(i - 1); + TScalar alpha = + Structure.IsEqualApprox(Structure.Abs(x[i]), StaticAccess.AdditionUnit) + ? Structure.MultiplicationUnit + : Structure.Division(x[i], Structure.Abs(x[i])); + alpha = Structure.RealMul(alpha, x.Magnitude); + TVector ui = x + alpha * ei; + TVector vi = ui.Normalize; + TScalar dm = vi.Conj() * x; + TScalar omg = Structure.Division(x.Conj() * vi, dm); + if (Structure.IsEqualApprox(dm, StaticAccess.AdditionUnit)) + omg = Structure.MultiplicationUnit; + TMatrix qi = + id - Structure.Addition( + Structure.MultiplicationUnit, + omg + ) * OuterProduct(vi, vi); + for (int j = 1; j < i; j++) + qi[j, j] = Structure.MultiplicationUnit; + qia.SelfLMul(qi); + q.SelfMul(qi.Hermitian()); + + } + + return new QRPair(q, qia); + } + /// + /// structure of scalar field + /// + protected static readonly FieldStructure Structure = FieldStructure.Dispatch(); + /// + /// M = v1 v2* + /// + /// + /// + /// + public static TMatrix OuterProduct(TVector v1, TVector v2) + { + TMatrix res = new(); + TVector h2 = v2.Conj(); + int dim = v1.Dim; + for (int i = 1; i <= dim; i++) + for (int j = 1; j <= dim; j++) + res[i, j] = Structure.Multiplication(v1[i], h2[j]); + return res; + } + + + + /// + /// make the matrix a (dim - d)x(dim - d) matrix + /// M[i,j] = 0 for all i<=d or j<d + /// + /// + /// + public TMatrix UpTruncateBy(int d) + { + TMatrix res = Copy(); + for (int i = 1; i <= d; i++) + for (int j = 1; j <= d; j++) + res[i, j] = StaticAccess.AdditionUnit; + return res; + } + + /// + /// make the matrix a (dim - d) x (dim - d) matrix + /// M[i,j] = 0 for all i>=d and j>=d + /// + /// + /// + public TMatrix DownTruncateBy(int d) + { + TMatrix res = Copy(); + int dim = Dim; + for (int i = d; i <= dim; i++) + for (int j = d; j <= dim; j++) + res[i, j] = StaticAccess.AdditionUnit; + return res; + } + + private JDPair InternalUpperTriangleUnitaryDecomposition() + { + TMatrix u = new(); + TMatrix d = new(); + for (int i = 1; i <= Dim; i++) + { + TVector col = new(); + d[i, i] = this[i, i]; + for (int r = i; r >= 1; r--) + { + if (Structure.IsEqualApprox(this[i, i], this[r, r])) + col[r] = StaticAccess.MultiplicationUnit; + else + col[r] = Structure.Division( + Structure.AdditionInverse(col * this[r]), + Structure.Subtraction(this[r, r], this[i, i]) + ); + } + + u[i, false] = col.Normalize; + } + + return new JDPair(u, d); + } + + internal JDPair InternalUnitaryDecomposition(bool triangle = false) + { + int maxIters = 5000; + double firstTol = 1E-6; + double secondTol = 1E-10; + + TMatrix shift = VarOne * this[Dim, Dim]; + QRPair qr = (this - shift).QRDecomposition(); + TMatrix a = qr.R * qr.Q + shift; + TMatrix u = qr.Q.Hermitian(); + int wIndex = Dim; + double diff = Structure.MaxError(this[wIndex, wIndex - 1]); + double errorDecrease = Structure.MaxError(a[wIndex, wIndex - 1]) - diff; + for (int i = 0; i < maxIters; i++) + { + if (i == maxIters - 10) + { + Console.WriteLine("HERE"); + } + + diff = Structure.MaxError(a[wIndex, wIndex - 1]); + TScalar s = diff < firstTol + ? a[wIndex, wIndex] + : Structure.WilkinsonShift( + a[wIndex - 1, wIndex - 1], + a[wIndex - 1, wIndex], + a[wIndex, wIndex - 1], + a[wIndex, wIndex] + ); + for (int j = 1; j <= Dim; j++) + shift[j, j] = s; + qr = (a - shift).QRDecomposition(); + a = qr.R * qr.Q + shift; + u = qr.Q.Hermitian() * u; + errorDecrease = Structure.MaxError(a[wIndex, wIndex - 1]) - diff; + diff = Structure.MaxError(a[wIndex, wIndex - 1]); + if (a.Property.IsUpperTriangle) + { + if (triangle || a.Property.IsDiagonal) + return new JDPair(u.Hermitian(), a); + JDPair tud = a.InternalUpperTriangleUnitaryDecomposition(); + return new JDPair(u.Hermitian() * tud.J, tud.D); + + } + + if (diff < secondTol || errorDecrease < firstTol) + { + int cIndex = wIndex; + double cMaxErr = diff; + for (int j = 2; j <= Dim; j++) + { + if (Structure.MaxError(a[j, j - 1]) > cMaxErr) + { + cIndex = j; + cMaxErr = Structure.MaxError(a[j, j - 1]); + } + } + + wIndex = cIndex; + } + } + + Console.WriteLine(PythonRepresentation); + Console.WriteLine($"--{wIndex} | {diff} | {errorDecrease}"); + Console.WriteLine(a.PythonRepresentation); + throw new Exception("CONVERGE MIGHT FAILED"); + } + + /// + /// eigen values of the matrix + /// + public virtual IEnumerable EigenValues() => InternalUnitaryDecomposition(true).D.Diagonals; + + +} + + diff --git a/src/Algebra/Matrices/MatrixProperty/MatrixProperty.cs b/src/Algebra/Matrices/MatrixProperty/MatrixProperty.cs new file mode 100644 index 0000000..f74ee09 --- /dev/null +++ b/src/Algebra/Matrices/MatrixProperty/MatrixProperty.cs @@ -0,0 +1,102 @@ +using Skeleton.Algebra.ScalarFieldStructure; +using Skeleton.Algebra.Vectors; + +namespace Skeleton.Algebra.Matrices.MatrixProperty; + +/// +/// test properties of matrix +/// +/// +/// +/// +public class MatrixProperty + where TVector : Vector, new() + where TMatrix : Matrix, new() + where TScalar : notnull +{ + private TMatrix Mat { get; } + private static readonly FieldStructure Structure = FieldStructure.Dispatch(); + internal MatrixProperty(TMatrix m) => Mat = m; + + /// + /// M[i, j] = 0 if i!= j + /// + public bool IsDiagonal => Mat.OffDiagonals + .All(x => Structure.IsEqualApprox(x, StaticAccess.AdditionUnit)); + + /// + /// exist of inv(M) + /// + /// + public bool IsInvertible => !Structure.IsEqualApprox(Mat.Det(), StaticAccess.AdditionUnit); + + /// + /// M M* = M* M + /// + public bool IsNormal => (Mat * Mat.Hermitian()).IsEqualApprox(Mat.Hermitian() * Mat); + + /// + /// M = M* + /// + public bool IsHermitian => Mat.IsEqualApprox(Mat.Hermitian()); + + /// + /// M[i, j] = 0 for all j < i + /// + public bool IsUpperTriangle + { + get + { + for (int i = 2; i <= Mat.Dim; i++) + for (int j = 1; j < i; j++) + if (!Structure.IsEqualApprox(Mat[i, j], StaticAccess.AdditionUnit)) + return false; + return true; + } + } + + /// + /// M[i, j] = 0 for all j > i + /// + public bool IsLowerTriangle + { + get + { + for (int i = 1; i <= Mat.Dim; i++) + for (int j = i + 1; j <= Mat.Dim; j++) + if (!Structure.IsEqualApprox(Mat[i, j], StaticAccess.AdditionUnit)) + return false; + return true; + } + } + + /// + /// U*U = UU* = I + /// + public bool IsUnitary => (Mat * Mat.Hermitian()).IsEqualApprox(Mat.VarOne); + + /// + /// U + U* = 0 + /// + public bool IsSkewHermitian => (Mat + Mat.Hermitian()).IsEqualApprox(Mat.VarZero); + + /// + /// tr(M) = 0 + /// + public bool IsTraceless => Structure.IsEqualApprox(Mat.Trace(), StaticAccess.AdditionUnit); + + /// + /// M M^T = I + /// + public bool IsOrthogonal => (Mat * Mat.Transpose()).IsEqualApprox(Mat.VarOne); + + /// + /// M = M^T + /// + public bool IsSymmetric => Mat.IsEqualApprox(Mat.Transpose()); + + /// + /// M = -M^T + /// + public bool IsSkewSymmetric => Mat.IsEqualApprox(-Mat.Transpose()); +} diff --git a/src/Algebra/Matrices/NormalMatrices/FDiagonalMatrix.cs b/src/Algebra/Matrices/NormalMatrices/FDiagonalMatrix.cs new file mode 100644 index 0000000..b7a6bb4 --- /dev/null +++ b/src/Algebra/Matrices/NormalMatrices/FDiagonalMatrix.cs @@ -0,0 +1,79 @@ +using Skeleton.Algebra.FixableObjects; +using Skeleton.Algebra.ScalarFieldStructure; + +namespace Skeleton.Algebra; + +public partial class CategoryOf +{ + public static partial class OnField + { + /// + public class FDiagonalMatrix : FNormalMatrix + { + /// + public FDiagonalMatrix() + { + U = One; + D = AsMatrix; + } + + /// + public FDiagonalMatrix(IEnumerable es, bool row = true) : base(es, row) => + this.PostConstructionFix(); + + /// + public FDiagonalMatrix(params TScalar[] es) + { + U = One; + if (es.Length == FDim) + for (int i = 1; i <= FDim; i++) + this[i, i] = es[i - 1]; + else + for (int i = 1; i <= FDim; i++) + for (int j = 1; j <= FDim; j++) + this[i, j] = es[i - 1 + (j - 1) * FDim]; + this.PostConstructionFix(); + } + internal override bool NeedFix => base.NeedFix || !Property.IsDiagonal; + internal override void Fix() + { + for (int i = 1; i <= Dim; i++) + { + for (int j = i + 1; j <= Dim; j++) + { + this[i, j] = StaticAccess.AdditionUnit; + this[j, i] = StaticAccess.AdditionUnit; + } + } + base.Fix(); + } + + /// + /// log of the normal matrix + /// + /// + public FDiagonalMatrix Log() + { + FDiagonalMatrix res = new(); + for (int i = 1; i <= Dim; i++) + res[i, i] = FieldStructure.Log(this[i, i]); + return res; + } + + /// + /// exp of the normal matrix + /// + /// + public FDiagonalMatrix Exp() + { + FDiagonalMatrix res = new(); + for (int i = 1; i <= Dim; i++) + res[i, i] = FieldStructure.Exp(this[i, i]); + return res; + } + + internal override void UnitaryFix(FMatrix a) { } + + } + } +} diff --git a/src/Algebra/Matrices/NormalMatrices/FHermitianMatrix.cs b/src/Algebra/Matrices/NormalMatrices/FHermitianMatrix.cs new file mode 100644 index 0000000..0ade2c0 --- /dev/null +++ b/src/Algebra/Matrices/NormalMatrices/FHermitianMatrix.cs @@ -0,0 +1,99 @@ +using System.Numerics; +using Skeleton.Algebra.AdditionalProperties; +using Skeleton.Algebra.AdditionalProperties.Extensions; +using Skeleton.Algebra.FixableObjects; + +namespace Skeleton.Algebra; + +public static partial class CategoryOf +{ + /// + public class FHermitianMatrix : + OnField.FNormalMatrix, + IUnitaryInvariantMatrix + { + /// + /// id hermitian matrix + /// + public new static readonly FHermitianMatrix One = new(OnField.FMatrix.One); + /// + /// zero hermitian matrix + /// + public new static readonly FHermitianMatrix Zero = new(); + /// + public FHermitianMatrix() { } + + /// + public FHermitianMatrix(IEnumerable.FVector> es, bool row = true) : base(es, row) => + this.PostConstructionFix(); + + /// + /// addition + /// + /// + /// + /// + public static FHermitianMatrix operator +(FHermitianMatrix a, FHermitianMatrix b) => new(a.AsMatrix + b); + + /// + /// scalar multiplication + /// + /// hermitian matrix + /// scalar + /// + public static FHermitianMatrix operator *(FHermitianMatrix a, double b) => new(a.AsMatrix * b); + /// + /// scalar multiplication + /// + /// hermitian matrix + /// scalar + /// + public static FHermitianMatrix operator *(double b, FHermitianMatrix a) => a * b; + + /// + /// rayleigh quotient + /// + /// + /// + /// + public static Complex operator /(FHermitianMatrix a, OnField.FVector b) => a.Rayleigh(b); + + /// + /// Rayleigh quotient + /// + /// + /// + public double Rayleigh(OnField.FVector v) + { + Complex t = 0; + Complex w = 0; + for (int i = 1; i <= Dim; i++) + { + w += v[i] * Complex.Conjugate(v[i]); + for (int j = 1; j <= Dim; j++) + t += Complex.Conjugate(v[i]) * v[j] * this[j, i]; + } + + return (t / w).Real; + } + internal override void Fix() + { + for (int i = 1; i <= Dim; i++) + { + for (int j = i+1; j <= Dim; j++) + { + Complex w = (this[i, j] + Complex.Conjugate(this[j, i])) / 2d; + this[i, j] = w; + this[j, i] = Complex.Conjugate(w); + } + } + } + internal override bool NeedFix => base.NeedFix || !Hermitian().IsEqualApprox(this); + + /*/// + public override JDPair.FMatrix> EigenDecomposition() + => this.UnitaryDecomposition();*/ + + internal override void UnitaryFix(OnField.FMatrix a) => a.UnitaryFix(); + } +} diff --git a/src/Algebra/Matrices/NormalMatrices/FLieUnitaryMatrix.cs b/src/Algebra/Matrices/NormalMatrices/FLieUnitaryMatrix.cs new file mode 100644 index 0000000..09412a7 --- /dev/null +++ b/src/Algebra/Matrices/NormalMatrices/FLieUnitaryMatrix.cs @@ -0,0 +1,158 @@ +using System.Numerics; +using Skeleton.Algebra.AdditionalProperties; +using Skeleton.Algebra.AdditionalProperties.Extensions; +using Skeleton.Algebra.Extensions; +using Skeleton.Algebra.FixableObjects; + +namespace Skeleton.Algebra; + +public static partial class CategoryOf +{ + /// + /// equivalence class of Skew hermitian matrices
+ /// A~B if eigenvalues of A = eigenvalues of B + 2 k pi i + /// + public class FLieUnitaryMatrix : + OnField.FNormalMatrix, + IUnitaryInvariantMatrix + { + /// + public FLieUnitaryMatrix() { } + /// + public FLieUnitaryMatrix(IEnumerable.FVector> es, bool row = true) : base(es, row) => + this.PostConstructionFix(); + + /// + /// Construct by unitary matrix and diagonal matrix + /// + /// + /// + public FLieUnitaryMatrix(OnField.FMatrix u, OnField.FMatrix d) : base(u, d) + => this.PostConstructionFix(); + /// + public FLieUnitaryMatrix(params Complex[] es) : base(es) => this.PostConstructionFix(); + + /// + public FLieUnitaryMatrix(OnField.FNormalMatrix n) : base(n) => this.PostConstructionFix(); + + /// + /// zero matrix as skew hermitian matrix + /// + public new static readonly FLieUnitaryMatrix Zero = new(); + /// + /// exp of the matrix + /// + /// + public FUnitaryMatrix Exp() + { + //JDPair.FMatrix> jd = EigenDecomposition(); + //FUnitaryMatrix u = jd.J.MatrixCast(); + FUnitaryMatrix expd = new(); + for (int i = 1; i <= Dim; i++) + expd[i, i] = Complex.Exp(D[i, i]); + return new(U, expd); + } + private bool LieNeedFix + { + get + { + if (U == Zero) + return true; + return + !(this + Hermitian()).IsEqualApprox(Zero) || + !Trace().Real.IsEqualApprox(0) || + D.Diagonals.Any(x => x.Imaginary < -Math.PI || x.Imaginary > Math.PI); + } + } + internal override bool NeedFix => base.NeedFix || LieNeedFix; + internal static void EigenValueFix(Complex[] eigenValues) + { + for (int i = 0; i < eigenValues.Length; i++) + { + double posDiff = eigenValues[i].Imaginary - Math.PI; + double negDiff = -eigenValues[i].Imaginary - Math.PI; + if (posDiff > 0) + eigenValues[i] = + (eigenValues[i].Imaginary - + (1 + Math.Floor(posDiff / (2 * Math.PI))) * 2 * Math.PI) * Complex.ImaginaryOne; + else if (negDiff > 0) + eigenValues[i] = + (eigenValues[i].Imaginary + + (1 + Math.Floor(negDiff / (2 * Math.PI))) * 2 * Math.PI) * Complex.ImaginaryOne; + } + } + internal override void UnitaryFix(OnField.FMatrix a) => a.UnitaryFix(); + internal override void InternalFix() + { + Complex[] eigenValues = D.Diagonals.ToArray(); + EigenValueFix(eigenValues); + for (int i = 0; i < eigenValues.Length; i++) + D[i + 1, i + 1] = eigenValues[i].Imaginary * Complex.ImaginaryOne; + } + /// + /// scalar multiplication + /// + /// matrix + /// scalar + /// + public static FLieUnitaryMatrix operator *(FLieUnitaryMatrix a, double b) => + new(a.AsMatrix * b); + /// + /// scalar multiplication + /// + /// scalar + /// matrix + /// + public static FLieUnitaryMatrix operator *(double b, FLieUnitaryMatrix a) => + new(a.AsMatrix * b); + /// + /// Cayley transformation (I + A)/(I - A) + /// + /// + public FUnitaryMatrix CayleyPositive() => new((One + this) / (One - this)); + /// + /// Cayley transformation (I - A)/(I + A) + /// + /// + public FUnitaryMatrix CayleyNegative() => new((One - this) / (One + this)); + /// + /// addition + /// + /// + /// + /// + public static FLieUnitaryMatrix operator +(FLieUnitaryMatrix a, FLieUnitaryMatrix b) => + new(a.AsMatrix + b); + /// + /// subtraction + /// + /// + /// + /// + public static FLieUnitaryMatrix operator -(FLieUnitaryMatrix a, FLieUnitaryMatrix b) => + new(a.AsMatrix - b); + /// + /// commutator + /// + /// + /// + /// + public static FLieUnitaryMatrix operator ^(FLieUnitaryMatrix a, FLieUnitaryMatrix b) => + new(Commutator(a, b)); + /// + /// anti commutator + /// + /// + /// + /// + public static FHermitianMatrix operator |(FLieUnitaryMatrix a, FLieUnitaryMatrix b) => + new(AntiCommutator(a, b)); + /* + /// + public override JDPair.FMatrix> EigenDecomposition() => + this.UnitaryDecomposition(true); + */ + + } +} + diff --git a/src/Algebra/Matrices/NormalMatrices/FNormalMatrix.cs b/src/Algebra/Matrices/NormalMatrices/FNormalMatrix.cs new file mode 100644 index 0000000..7f6f996 --- /dev/null +++ b/src/Algebra/Matrices/NormalMatrices/FNormalMatrix.cs @@ -0,0 +1,106 @@ +using Skeleton.Algebra.FixableObjects; +using Skeleton.Algebra.ScalarFieldStructure; +using Skeleton.DataStructure.Packs.MatrixDecompositions; + +namespace Skeleton.Algebra; + +public partial class CategoryOf +{ + public static partial class OnField where TScalar : notnull + { + + /// + /// matrices that commutative with its hermitian + /// NN* = N*N + /// + public class FNormalMatrix : FMatrix + { + /// + /// Unitary matrix + /// + public FMatrix U { get; set; } = Zero; + /// + /// Diagonal matrix + /// + public FMatrix D { get; set; } = Zero; + + /// + public FNormalMatrix() { } + + /// + public FNormalMatrix(IEnumerable es, bool row = true) : base(es, row) => + this.PostConstructionFix(); + + /// + public FNormalMatrix(params TScalar[] es) : base(es) => + this.PostConstructionFix(); + + /// + /// construct by unitary matrix and diagonal matrix + /// + /// unitary matrix + /// diagonal matrix + public FNormalMatrix(FMatrix u, FMatrix d) + { + U = u; + D = d; + ConstructFromUnitaryAndDiagonal(); + } + + /// + public FNormalMatrix(FNormalMatrix n) + { + U = n.U; + D = n.D; + ConstructFromUnitaryAndDiagonal(); + } + + private void ConstructFromUnitaryAndDiagonal() + { + ResetOne(); + SelfMul(U); + SelfMul(D); + SelfMul(U.Hermitian()); + } + internal override bool NeedFix => U == Zero; + internal virtual void UnitaryFix(FMatrix a) { } + + internal override void Fix() + { + if (U == Zero) + { + JDPair jd = InternalUnitaryDecomposition(); + U = jd.J; + D = jd.D; + } + + if (!(U * U.Hermitian()).IsEqualApprox(One)) + UnitaryFix(U); + for (int i = 1; i <= FDim; i++) + { + for (int j = 1 + 1; j <= FDim; j++) + { + this[i, j] = StaticAccess.AdditionUnit; + this[j, i] = StaticAccess.AdditionUnit; + } + } + InternalFix(); + ConstructFromUnitaryAndDiagonal(); + } + + /// + /// + /// + /// + public new FNormalMatrix Hermitian() => new(U, D.Hermitian()); + internal virtual void InternalFix() { } + + /// + public override JDPair EigenDecomposition() => new (U, D); + + /// + public override IEnumerable EigenValues() => D.Diagonals; + } + } +} + diff --git a/src/Algebra/Matrices/NormalMatrices/FSpecialLieUnitaryMatrix.cs b/src/Algebra/Matrices/NormalMatrices/FSpecialLieUnitaryMatrix.cs new file mode 100644 index 0000000..0c222c8 --- /dev/null +++ b/src/Algebra/Matrices/NormalMatrices/FSpecialLieUnitaryMatrix.cs @@ -0,0 +1,111 @@ +using System.Numerics; +using Skeleton.Algebra.AdditionalProperties; +using Skeleton.Algebra.Extensions; +using Skeleton.Algebra.FixableObjects; + +namespace Skeleton.Algebra; + +public static partial class CategoryOf +{ + /// + /// equivalence class of skew hermitian matrices with trace = 2 k pi i + /// by A~B if A eigenvalues of A = eigenvalues of B + 2 k pi i + /// + public class FSpecialLieUnitaryMatrix : + FLieUnitaryMatrix, + IUnitaryInvariantMatrix + { + /// + /// zero matrix as skew hermitian matrix with trace = 2 k pi i + /// + public new static readonly FSpecialLieUnitaryMatrix Zero = new(); + + /// + public FSpecialLieUnitaryMatrix() { } + + /// + public FSpecialLieUnitaryMatrix(IEnumerable.FVector> es, bool row = true) : + base(es, row) => this.PostConstructionFix(); + + /// + public FSpecialLieUnitaryMatrix(params Complex[] es) : base(es) => this.PostConstructionFix(); + + /// + public FSpecialLieUnitaryMatrix(OnField.FMatrix u, OnField.FMatrix d) : base(u, d) => + this.PostConstructionFix(); + + /// + public FSpecialLieUnitaryMatrix(OnField.FNormalMatrix n) : base(n) => this.PostConstructionFix(); + + internal override bool NeedFix => base.NeedFix || !Complex.Exp(Trace()).IsEqualApprox(1); + + /// + /// addition of matrices + /// + /// left operand + /// right operand + /// + public static FSpecialLieUnitaryMatrix operator +(FSpecialLieUnitaryMatrix a, FSpecialLieUnitaryMatrix b) => + new(a.AsMatrix + b); + + /// + /// scalar multiplication
+ /// LieSU x R -> LieSU + /// + /// + /// + /// + public static FSpecialLieUnitaryMatrix operator *(FSpecialLieUnitaryMatrix a, double b) => + new(a.U, a.D * b); + + /// + /// scalar multiplication
+ /// LieSU x R -> LieSU + /// + /// + /// + /// + public static FSpecialLieUnitaryMatrix operator *(double b, FSpecialLieUnitaryMatrix a) => a * b; + + /// + /// exp of the matrix + /// + /// + public new FSpecialUnitaryMatrix Exp() + { + OnField.FMatrix expd = new(); + for (int i = 1; i <= FDim; i++) + expd[i, i] = Complex.Exp(D[i, i]); + return new(U, expd); + } + + internal override void InternalFix() + { + base.InternalFix(); + Complex e = D.Trace(); + double eMin = e.Imaginary; + for (int i = 1; i <= Dim; i++) + { + double errP = e.Imaginary - i * Math.PI; + double errN = e.Imaginary + i * Math.PI; + if (Math.Abs(errP) < Math.Abs(eMin)) + eMin = errP; + if (Math.Abs(errN) < Math.Abs(eMin)) + eMin = errN; + } + for (int i = 1; i <= Dim; i++) + D[i, i] = (D[i, i].Imaginary - eMin / Dim) * Complex.ImaginaryOne; + } + + /// + /// commutator + /// + /// + /// + /// + public static FSpecialLieUnitaryMatrix operator ^(FSpecialLieUnitaryMatrix a, FSpecialLieUnitaryMatrix b) => + new(Commutator(a, b)); + + } +} + diff --git a/src/Algebra/Matrices/NormalMatrices/FSpecialUnitaryMatrix.cs b/src/Algebra/Matrices/NormalMatrices/FSpecialUnitaryMatrix.cs new file mode 100644 index 0000000..72d6a7e --- /dev/null +++ b/src/Algebra/Matrices/NormalMatrices/FSpecialUnitaryMatrix.cs @@ -0,0 +1,94 @@ +using System.Numerics; +using Skeleton.Algebra.AdditionalProperties; +using Skeleton.Algebra.FixableObjects; + +namespace Skeleton.Algebra; + +public static partial class CategoryOf +{ + + /// + /// Unitary matrices with det = 1 + /// + public class FSpecialUnitaryMatrix : + FUnitaryMatrix, + IUnitaryInvariantMatrix + { + /// + public FSpecialUnitaryMatrix() => ResetOne(); + + /// + public FSpecialUnitaryMatrix(IEnumerable.FVector> es, bool row = true) : base(es, row) => + this.PostConstructionFix(); + + /// + public FSpecialUnitaryMatrix(params Complex[] es) : base(es) => this.PostConstructionFix(); + + /// + public FSpecialUnitaryMatrix(OnField.FNormalMatrix n) : base(n) => this.PostConstructionFix(); + + /// + public FSpecialUnitaryMatrix(OnField.FMatrix u, OnField.FMatrix d) : base(u, d) => + this.PostConstructionFix(); + + /// + /// id matrix as special unitary matrix + /// + public new static readonly FSpecialUnitaryMatrix One = new(); + + /// + /// + /// + /// + public new FSpecialLieUnitaryMatrix Log() + { + OnField.FMatrix logd = new(); + for (int i = 1; i <= FDim; i++) + logd[i, i] = Complex.Log(D[i, i]); + return new(U, logd); + } + + /// + /// + /// + /// + public new FSpecialUnitaryMatrix Inv() => Hermitian(); + + internal override void InternalFix() + { + base.InternalFix(); + Complex det = 1; + for (int i = 1; i <= FDim; i++) + det *= D[i, i]; + Complex d = Complex.Pow(1d / det, 1d / Dim); + for (int i = 1; i <= FDim; i++) + D[i, i] *= d; + } + + /// + /// transpose + conjugate + /// + /// + public new FSpecialUnitaryMatrix Hermitian() => new(U, D.Hermitian()); + + /// + /// multiplication + /// + /// + /// + /// + public static FSpecialUnitaryMatrix operator *(FSpecialUnitaryMatrix a, FSpecialUnitaryMatrix b) => + new(a.AsMatrix * b); + + /// + /// a/b => a b* + /// + /// + /// + /// + public static FSpecialUnitaryMatrix operator /(FSpecialUnitaryMatrix a, FSpecialUnitaryMatrix b) => + a * b.Hermitian(); + + } + +} diff --git a/src/Algebra/Matrices/NormalMatrices/FUnitaryMatrix.cs b/src/Algebra/Matrices/NormalMatrices/FUnitaryMatrix.cs new file mode 100644 index 0000000..85a6f79 --- /dev/null +++ b/src/Algebra/Matrices/NormalMatrices/FUnitaryMatrix.cs @@ -0,0 +1,119 @@ +using System.Numerics; +using Skeleton.Algebra.AdditionalProperties; +using Skeleton.Algebra.AdditionalProperties.Extensions; +using Skeleton.Algebra.Extensions; +using Skeleton.Algebra.FixableObjects; + +namespace Skeleton.Algebra; + +public static partial class CategoryOf +{ + /// + /// Matrices whose inverse is its hermitian + /// AA* = I = A*A + /// + public class FUnitaryMatrix : + OnField.FNormalMatrix, + IUnitaryInvariantMatrix + { + /// + public FUnitaryMatrix() => ResetOne(); + + /// + public FUnitaryMatrix(IEnumerable.FVector> es, bool row = true) : base(es, row) => + this.PostConstructionFix(); + + /// + public FUnitaryMatrix(params Complex[] es) : base(es) => this.PostConstructionFix(); + + /// + public FUnitaryMatrix(OnField.FMatrix u, OnField.FMatrix d) : base(u, d) => + this.PostConstructionFix(); + + /// + public FUnitaryMatrix(OnField.FNormalMatrix n) : base(n) => + this.PostConstructionFix(); + /// + /// id matrix as unitary matrix + /// + public new static readonly FUnitaryMatrix One = new(); + + /// + /// inverse of the matrix + /// + /// + public new FUnitaryMatrix Inv() => Hermitian(); + + /// + /// multiplication between unitary matrices + /// + /// + /// + /// + public static FUnitaryMatrix operator *(FUnitaryMatrix self, FUnitaryMatrix other) => + new(self.AsMatrix * other); + + /// + /// transpose + conjugate + /// + /// + public new FUnitaryMatrix Hermitian() => new(U, D.Hermitian()); + + internal override void InternalFix() + { + for (int i = 1; i <= FDim; i++) + D[i, i] /= D[i, i].Magnitude; + } + + private bool UnitaryNeedFix => + !D.Diagonals.Aggregate((a, b) => a * b).Magnitude.IsEqualApprox(1); + internal override bool NeedFix => base.NeedFix || UnitaryNeedFix; + + /// + /// a conjugate on b => a b inv(a) + /// + /// + /// + /// + public TUnitaryInvariant ConjugateOn(TUnitaryInvariant a) + where TUnitaryInvariant : OnField.FMatrix, IUnitaryInvariantMatrix, new() + { + TUnitaryInvariant res = new(); + res.Elements = (this * a * Hermitian()).ToArray(); + res.Fix(); + return res; + } + + /// + /// conjugate of real is itself
+ /// so if homomorphism h: R -> K kept the conj
+ /// h(x) = h(conj(x)) = conj(h(x))
+ /// h will map real diag matrix to hermitian matrix in K + /// + /// + /// + public FHermitianMatrix ConjugateOn(OnField.FDiagonalMatrix a) + { + FHermitianMatrix g = a.AsComplex().MatrixCast(); + return ConjugateOn(g); + } + + /// + /// log of unitary matrix + /// + /// + public FLieUnitaryMatrix Log() + { + OnField.FMatrix logd = new(); + for (int i = 1; i <= FDim; i++) + logd[i, i] = Complex.Log(D[i, i]); + return new(U, logd); + } + + /*/// + public override JDPair.FMatrix> EigenDecomposition() => this.UnitaryDecomposition();*/ + + internal override void UnitaryFix(OnField.FMatrix a) => a.UnitaryFix(); + } +} + diff --git a/src/Algebra/Morphisms/Homomorphism.cs b/src/Algebra/Morphisms/Homomorphism.cs new file mode 100644 index 0000000..f1a3eef --- /dev/null +++ b/src/Algebra/Morphisms/Homomorphism.cs @@ -0,0 +1,13 @@ +namespace Skeleton.Algebra.Morphisms; + +/// +/// a morphism K -> D that preserve some structures of K +/// +/// +/// +/// property ketp by morphism +public class Homomorphism : Morphism +where TProperty : Morphism +{ + +} \ No newline at end of file diff --git a/src/Algebra/Morphisms/Isomorphism.cs b/src/Algebra/Morphisms/Isomorphism.cs new file mode 100644 index 0000000..dcf1f6a --- /dev/null +++ b/src/Algebra/Morphisms/Isomorphism.cs @@ -0,0 +1,11 @@ +namespace Skeleton.Algebra.Morphisms; + +/// +/// morphism K -> D that preserve all structures of K +/// +/// +/// +public class Isomorphism : Morphism +{ + +} \ No newline at end of file diff --git a/src/Algebra/Morphisms/Morphism.cs b/src/Algebra/Morphisms/Morphism.cs new file mode 100644 index 0000000..694de68 --- /dev/null +++ b/src/Algebra/Morphisms/Morphism.cs @@ -0,0 +1,11 @@ +namespace Skeleton.Algebra.Morphisms; + +/// +/// a map from K to D +/// +/// +/// +public class Morphism +{ + +} \ No newline at end of file diff --git a/src/Algebra/ScalarFieldStructure/ComplexFieldStructure.cs b/src/Algebra/ScalarFieldStructure/ComplexFieldStructure.cs new file mode 100644 index 0000000..d7749d4 --- /dev/null +++ b/src/Algebra/ScalarFieldStructure/ComplexFieldStructure.cs @@ -0,0 +1,121 @@ +using System.Numerics; + +namespace Skeleton.Algebra.ScalarFieldStructure; + +/// +public class ComplexFieldStructure : FieldStructure +{ + /// + /// + /// + public static readonly ComplexFieldStructure Structure = new(); + + /// + public override Complex Addition(Complex self, Complex other) => self + other; + + /// + public override Complex AdditionUnit => 0d; + + /// + public override Complex Multiplication(Complex self, Complex other) => self * other; + + /// + public override Complex MultiplicationUnit => 1d; + + /// + public override Complex AdditionInverse(Complex self) => -self; + + /// + public override Complex MultiplicationInverse(Complex self) => 1d / self; + + /// + public override Complex Subtraction(Complex self, Complex other) => self - other; + + /// + public override Complex Division(Complex self, Complex other) => self / other; + + /// + public override bool IsEqualApprox(Complex self, Complex other, double? absTol = null, double? relTol = null) + + => RealFieldStructure.Structure.IsEqualApprox(self.Real, other.Real, absTol, relTol) && + RealFieldStructure.Structure.IsEqualApprox(self.Imaginary, other.Imaginary, absTol, relTol); + + /// + public override Complex Fix(Complex x) + => base.Fix(new Complex( + RealFieldStructure.Structure.Fix(x.Real), + RealFieldStructure.Structure.Fix(x.Imaginary) + )); + + /// + public override Complex Conj(Complex self) => Complex.Conjugate(self); + + /// + public override double AsReal(Complex x) => x.Real; + + /// + public override Complex FromReal(double x) => x; + + private const string PlusSign = "+"; + private const string EmptyString = ""; + + /// + public override string RawRepresentation(Complex x) + => RealFieldStructure.Structure.IsEqualApprox(x.Imaginary, 0d) + ? RealFieldStructure.Structure.Representation(x.Real) + : $"{x.Real}{(x.Imaginary > 0 ? PlusSign : EmptyString)}{x.Imaginary}i"; + + /// + public override string RawPythonRepresentation(Complex x) + => RealFieldStructure.Structure.IsEqualApprox(x.Imaginary, 0d) + ? RealFieldStructure.Structure.Representation(x.Real) + : $"{x.Real}{(x.Imaginary > 0 ? PlusSign : EmptyString)}{x.Imaginary}j"; + + /// + public override string RawCSharpRepresentation(Complex x) + => RealFieldStructure.Structure.IsEqualApprox(x.Imaginary, 0d) + ? RealFieldStructure.Structure.Representation(x.Real) + : $"{x.Real}{(x.Imaginary >= 0 ? PlusSign : EmptyString)}{x.Imaginary}*AlgebraConstant.I"; + + /// + public override string DumpString(Complex x) + => "" + + RealFieldStructure.Structure.DumpString(x.Real) + + "" + + RealFieldStructure.Structure.DumpString(x.Imaginary) + + ""; + + /// + public override Complex Resolve(string dumpString) + { + double[] restored = dumpString + .Replace("", "") + .Replace("", "") + .Split("") + .Select(x => RealFieldStructure.Structure.Resolve(x)) + .ToArray(); + return new Complex(restored[0], restored[1]); + } + + /// + public override double MaxError(Complex a) + => Math.Max(Math.Abs(a.Real), Math.Abs(a.Imaginary)); + + /// + public override Complex ComplexCast(Complex a) => a; + + /// + public override Complex RealMul(Complex a, double b) => a * b; + + /// + public override Complex SquareRoot(Complex a) => Complex.Sqrt(a); + + /// + public override Complex FromComplex(Complex a) => a; + + /// + public override Complex Log(Complex x) => Complex.Log(x); + + /// + public override Complex Exp(Complex x) => Complex.Exp(x); +} diff --git a/src/Algebra/ScalarFieldStructure/FieldStructure.cs b/src/Algebra/ScalarFieldStructure/FieldStructure.cs new file mode 100644 index 0000000..c26783f --- /dev/null +++ b/src/Algebra/ScalarFieldStructure/FieldStructure.cs @@ -0,0 +1,325 @@ +using System.Numerics; +using Skeleton.Algebra.Extensions; +using Skeleton.DataStructure.Packs; + +namespace Skeleton.Algebra.ScalarFieldStructure; + +/// +/// operations / structure of the field +/// +public abstract class FieldStructure +{ + +} + +/// +public abstract class FieldStructure : FieldStructure +{ + + /// + /// return the singleton base on type + /// + /// + public static FieldStructure Dispatch() + { + if (typeof(TScalar) == typeof(double)) + return (RealFieldStructure.Structure as FieldStructure)!; + if (typeof(TScalar) == typeof(Complex)) + return (ComplexFieldStructure.Structure as FieldStructure)!; + return ScalarFieldStructureProvider.GetOperation(); + } + + /// + public FieldStructure() + { + StaticAccess.AdditionUnit = AdditionUnit; + StaticAccess.MultiplicationUnit = MultiplicationUnit; + StaticAccess.Structure = this; + } + + /// + /// addition in field + /// + /// + /// + /// + public abstract TScalar Addition(TScalar self, TScalar other); + /// + /// addition unit of the field
+ /// a + unit = a for any a in field + /// + public abstract TScalar AdditionUnit { get; } + /// + /// multiplication in field + /// + /// + /// + /// + public abstract TScalar Multiplication(TScalar self, TScalar other); + /// + /// multiplication unit of the field
+ /// unit * a = a for any a in the field except the addition unit + /// + public abstract TScalar MultiplicationUnit { get; } + /// + /// inverse of element under addition
+ /// a + AddInv(a) = unit(addition) for any a in field + /// + /// + /// + public abstract TScalar AdditionInverse(TScalar self); + /// + /// inverse of element under multiplication
+ /// a * MulInv(a) = unit(multiplication) for any a in field except the addition unit + /// + /// + /// + public abstract TScalar MultiplicationInverse(TScalar self); + /// + /// a + AddInv(b) + /// + /// + /// + /// + public virtual TScalar Subtraction(TScalar self, TScalar other) + => Addition(self, AdditionInverse(other)); + + /// + /// a * MulInv(b) + /// + /// + /// + /// + public virtual TScalar Division(TScalar self, TScalar other) + => Multiplication(self, MultiplicationInverse(other)); + /// + /// conjugate of the scalar + /// + /// + /// + public abstract TScalar Conj(TScalar self); + /// + /// if two elements are close to each other in the field + /// + /// + /// + /// + /// + /// + public abstract bool IsEqualApprox(TScalar self, TScalar other, double? absTol=null, double? relTol = null); + + /// + /// tolerance to check if two scalars are close + /// + public double AbsoluteTolerance { get; set; } = 1E-8d; + + /// + /// + /// + public double RelativeTolerance { get; set; } = 1E-5d; + + /// + /// if calculation in field won't introduce error, override this and leave it empty + /// otherwise, provide method to reduce the error here + /// + /// + /// + public virtual TScalar Fix(TScalar x) + { + if (IsEqualApprox(x, AdditionUnit)) + return AdditionUnit; + if (IsEqualApprox(x, MultiplicationUnit)) + return MultiplicationUnit; + return x; + } + + /// + /// homomorphism from scalar field to real
+ /// ignore if not applicable + /// + /// + /// + public virtual double AsReal(TScalar x) => IsEqualApprox(x, AdditionUnit) ? 0 : 1; + /// + /// conjugated square xx* when field is real/complex + /// + /// + /// + public double ConjX2(TScalar x) => AsReal(Multiplication(x, Conj(x))); + + /// + /// homomorphism from real to scalar
+ /// ignore if not applicable + /// + /// + /// + public virtual TScalar FromReal(double x) => x.IsEqualApprox(0) ? AdditionUnit : MultiplicationUnit; + /// + /// readable string representation of data + /// + /// + /// + public string Representation(TScalar x) => RawRepresentation(Fix(x)); + + /// + /// representation of the original data, without fix + /// + /// + /// + public virtual string RawRepresentation(TScalar x) => $"{x}"; + + /// + /// code to construct object in python + /// + /// + /// + public string PythonRepresentation(TScalar x) => RawPythonRepresentation(Fix(x)); + + /// + /// code to construct scalar object in python
+ /// ignore if not applicable + /// + /// + /// + public virtual string RawPythonRepresentation(TScalar x) => ""; + /// + /// code to construct the object in c# + /// + /// + /// + public string CSharpRepresentation(TScalar x) => RawCSharpRepresentation(Fix(x)); + /// + /// C# string for the original data + /// + /// + /// + public abstract string RawCSharpRepresentation(TScalar x); + + /// + /// dump the object as string which can be restored exactly
+ /// ignore if not applicable
+ /// if override can not contain following string as substring
+ /// <D> </D> <C> </C> <CSPLIT/>
+ /// <VECTOR> <VSPLIT/> </VECTOR>
+ /// <MATRIX> <MSPLIT/> </MSPLIT> + /// + /// + /// + public virtual string DumpString(TScalar x) => ""; + + /// + /// method to restore a value from dump string
+ /// ignore if not applicable + /// + /// + /// + public virtual TScalar Resolve(string dumpString) => AdditionUnit; + /// + /// max diff from zero + /// + /// + /// + public abstract double MaxError(TScalar a); + + /// + /// homomorphism from field to complex
+ /// ignore if not applicable + /// + /// + /// + public virtual Complex ComplexCast(TScalar a) => IsEqualApprox(a, AdditionUnit) ? 0 : 1; + + + /// + /// K x R -> K
+ /// ignore if not applicable + /// + /// + /// + /// + public virtual TScalar RealMul(TScalar a, double b) => a; + + /// + /// for y in field, find x such that x*x = y + /// + /// + /// + public abstract TScalar SquareRoot(TScalar a); + + /// + /// Absolute value + /// + /// + /// + public TScalar Abs(TScalar a) => SquareRoot(Multiplication(a, Conj(a))); + /// + /// a x^2 + b x + c = unit(add) + /// + /// + /// + /// + /// + public QuadraticRoots QuadraticRoot(TScalar a, TScalar b, TScalar c) + { + TScalar invB = AdditionInverse(b); + TScalar b2 = Multiplication(b, b); + TScalar ax2 = Addition(a, a); + TScalar ac = Multiplication(a, c); + TScalar acx2 = Addition(ac, ac); + TScalar acx4 = Addition(acx2, acx2); + TScalar x1 = Division(Addition(invB, SquareRoot(Subtraction(b2, acx4))), ax2); + TScalar x2 = Division(Subtraction(invB, SquareRoot(Subtraction(b2, acx4))), ax2); + return new QuadraticRoots(x1, x2); + } + + internal TScalar WilkinsonShift(TScalar a1, TScalar b1, TScalar b2, TScalar a2) + { + QuadraticRoots es = QuadraticRoot( + AdditionUnit, + AdditionInverse(Addition(a1, a2)), + Subtraction(Multiplication(a1, a2), Multiplication(b1, b2)) + ); + TScalar k = a2; + TScalar s = Addition(Addition(Abs(a1), Abs(b1)), Addition(Abs(b2), Abs(a2))); + if (IsEqualApprox(s, AdditionUnit)) + return k; + TScalar q = Multiplication(Division(b1, s), Division(b2, s)); + + if (IsEqualApprox(q, AdditionUnit)) + return k; + TScalar p = Division( + Subtraction(Division(a1, s), Division(a2, s)), + Addition(MultiplicationUnit, MultiplicationUnit) + ); + TScalar r = SquareRoot(Addition(Multiplication(p, p), q)); + Complex cp = ComplexCast(p); + Complex cr = ComplexCast(r); + if (cp.Real * cr.Real + cp.Imaginary * cr.Imaginary < 0) + r = AdditionInverse(r); + k = Subtraction(k, Multiplication(s, Division(q, Addition(p, r)))); + return k; + } + + /// + /// homomorphism from complex to this field + /// + /// + /// + public virtual TScalar FromComplex(Complex a) => AdditionUnit; + + + /// + /// log of the element in the field + /// + /// + /// + public virtual TScalar Log(TScalar x) => AdditionUnit; + + /// + /// exp of the element in the field + /// + /// + /// + public virtual TScalar Exp(TScalar x) => MultiplicationUnit; +} + diff --git a/src/Algebra/ScalarFieldStructure/RealFieldStructure.cs b/src/Algebra/ScalarFieldStructure/RealFieldStructure.cs new file mode 100644 index 0000000..a8de10e --- /dev/null +++ b/src/Algebra/ScalarFieldStructure/RealFieldStructure.cs @@ -0,0 +1,86 @@ +using System.Numerics; + +namespace Skeleton.Algebra.ScalarFieldStructure; + +/// +public class RealFieldStructure : FieldStructure +{ + /// + /// + /// + public static readonly RealFieldStructure Structure = new(); + + /// + public override double Addition(double self, double other) => self + other; + + /// + public override double AdditionUnit => 0d; + + /// + public override double Multiplication(double self, double other) => self * other; + + /// + public override double MultiplicationUnit => 1d; + + /// + public override double AdditionInverse(double self) => -self; + + /// + public override double MultiplicationInverse(double self) => 1d / self; + + /// + public override double Division(double self, double other) => self / other; + + /// + public override double Subtraction(double self, double other) => self - other; + + /// + public override bool IsEqualApprox(double self, double other, double? absTol = null, double? relTol = null) + { + double absTolerance = absTol ?? AbsoluteTolerance; + double relTolerance = relTol ?? RelativeTolerance; + return Math.Abs(self - other) < absTolerance + relTolerance * Math.Max(Math.Abs(self), Math.Abs(other)); + } + /// + public override double Conj(double self) => self; + /// + public override double AsReal(double x) => x; + /// + public override double FromReal(double x) => x; + /// + public override string RawRepresentation(double x) => $"{x}"; + /// + public override string RawPythonRepresentation(double x) => RawRepresentation(x); + /// + public override string RawCSharpRepresentation(double x) => RawRepresentation(x); + + /// + public override string DumpString(double x) + => "" + string.Join(",", BitConverter.GetBytes(x)) + ""; + + /// + public override double Resolve(string dumpString) + => BitConverter.ToDouble( + dumpString + .Replace("", "") + .Replace("", "") + .Split(",") + .Select(x => Convert.ToByte(x)) + .ToArray() + ); + + /// + public override double MaxError(double a) => Math.Abs(a); + + /// + public override Complex ComplexCast(double a) => a; + + /// + public override double RealMul(double a, double b) => a * b; + + /// + public override double SquareRoot(double a) => Math.Sqrt(a); + + /// + public override double FromComplex(Complex a) => a.Real; +} diff --git a/src/Algebra/ScalarFieldStructure/ScalarFieldStructureProvider.cs b/src/Algebra/ScalarFieldStructure/ScalarFieldStructureProvider.cs new file mode 100644 index 0000000..4b7ba0c --- /dev/null +++ b/src/Algebra/ScalarFieldStructure/ScalarFieldStructureProvider.cs @@ -0,0 +1,19 @@ +namespace Skeleton.Algebra.ScalarFieldStructure; + +/// +/// +/// +public static class ScalarFieldStructureProvider +{ + /// + /// Register customer field structure + /// + /// + /// + public static void Register(FieldStructure structure) + => Cache[typeof(TScalar)] = structure; + + internal static FieldStructure GetOperation() + => (Cache[typeof(TScalar)] as FieldStructure)!; + private static readonly Dictionary Cache = new (); +} \ No newline at end of file diff --git a/src/Algebra/ScalarFieldStructure/StaticAccess.cs b/src/Algebra/ScalarFieldStructure/StaticAccess.cs new file mode 100644 index 0000000..cfce9d0 --- /dev/null +++ b/src/Algebra/ScalarFieldStructure/StaticAccess.cs @@ -0,0 +1,23 @@ +namespace Skeleton.Algebra.ScalarFieldStructure; + +/// +/// +/// +/// +public static class StaticAccess +{ + /// + /// + /// + public static T AdditionUnit { get; set; } + /// + /// + /// + public static T MultiplicationUnit { get; set; } + + /// + /// + /// + public static FieldStructure Structure { get; set; } + +} \ No newline at end of file diff --git a/src/Algebra/VectorSpaces/FVectorSpace.cs b/src/Algebra/VectorSpaces/FVectorSpace.cs new file mode 100644 index 0000000..e13b1ff --- /dev/null +++ b/src/Algebra/VectorSpaces/FVectorSpace.cs @@ -0,0 +1,25 @@ +using Skeleton.Algebra.ScalarFieldStructure; +using Skeleton.Algebra.VectorSpaces; + +namespace Skeleton.Algebra; + +public partial class CategoryOf +{ + public static partial class OnField where TScalar : notnull + { + + /// + public class FVectorSpace : VectorSpace + { + /// + public override int Dim => FDim; + internal override FieldStructure Structure => FieldStructure; + + /// + public FVectorSpace() { } + + /// + public FVectorSpace(IEnumerable es) : base(es) { } + } + } +} \ No newline at end of file diff --git a/src/Algebra/VectorSpaces/VectorSpace.cs b/src/Algebra/VectorSpaces/VectorSpace.cs new file mode 100644 index 0000000..f5abe1c --- /dev/null +++ b/src/Algebra/VectorSpaces/VectorSpace.cs @@ -0,0 +1,319 @@ +using System.Collections; +using Skeleton.Algebra.FixableObjects; +using Skeleton.Algebra.Matrices; +using Skeleton.Algebra.ScalarFieldStructure; +using Skeleton.Algebra.Vectors; +using Skeleton.DataStructure.Packs; + +namespace Skeleton.Algebra.VectorSpaces; + +/// +/// meta vector space +/// +public abstract class VectorSpace : FixableTensor + +{ + /// + /// number of independent basis + /// + public abstract int Rank { get; } + + /// + /// dimension of complementary space + /// + public int Nullity => Dim - Rank; +} +/// +public abstract class VectorSpace< TVectorSpace> : VectorSpace + + where TVectorSpace : VectorSpace< TVectorSpace>, new() +{ + /// + /// OTC of the space + /// + public abstract TVectorSpace OrthogonalComplementSpace { get; } + /// + /// the union vector space + /// + /// + /// + public abstract TVectorSpace UnionWith(TVectorSpace oth); +} +/// +public abstract class VectorSpace< TMatrix, TVectorSpace> : VectorSpace< TVectorSpace> + + where TMatrix : Matrix< TMatrix>, new() + where TVectorSpace : VectorSpace< TMatrix, TVectorSpace>, new() +{ + /// + /// Matrix that transform inner basis to basis + /// + protected TMatrix I2OTransformation { get; set; } = new(); + /// + /// Matrix that transform basis to inner basis + /// + protected TMatrix O2ITransformation { get; set; } = new(); + /// + /// I forgot what is this matrix + /// + protected TMatrix LTransformation { get; set; } = new(); + /// + /// union to vector spaces + /// + /// + /// + /// + public static TVectorSpace operator +(VectorSpace self, TVectorSpace other) + => self.UnionWith(other); + /// + /// intersection + /// + /// + /// + /// + public static TVectorSpace operator *(VectorSpace< TMatrix, TVectorSpace> self, TVectorSpace other) + => (self.OrthogonalComplementSpace + other.OrthogonalComplementSpace).OrthogonalComplementSpace; +} + +/// +public abstract class VectorSpace : VectorSpace, IEnumerable + + where TVector : Vector, new() + where TMatrix : Matrix, new() + where TVectorSpace : VectorSpace< TVector, TMatrix, TVectorSpace>, new() +{ + /// + public VectorSpace() { } + + /// + public VectorSpace(IEnumerable es) => Fill(es); + + /// + /// + /// + /// + /// + public abstract bool ContainsVector(TVector v); + /// + /// as iterator + /// + /// + public TVector this[int idx] => Basis[idx]; + /// + /// basis + /// + protected List Basis { get; set; } = new(); + /// + /// Standard basis + /// + protected List InnerBasis { get; set; } = new(); + /// + public sealed override int Rank => Basis.Count; + /// + /// init from numerable of vectors + /// + /// + public void Fill(IEnumerable bs) + { + InnerBasis = new List(); + Basis = new List(); + foreach (TVector b in bs) + AddBasis(b); + } + + /// + /// add a basis to the space + /// + /// + public abstract void AddBasis(TVector basis); + + + /// + /// union with another vector space + /// + /// + /// + public override TVectorSpace UnionWith(TVectorSpace oth) + { + TVectorSpace res = new TVectorSpace(); + foreach (TVector b in Basis) + res.AddBasis(b); + foreach (TVector b in oth.Basis) + res.AddBasis(b); + return res; + } + /// + public IEnumerator GetEnumerator() + { + foreach (TVector b in Basis) + yield return b; + } + IEnumerator IEnumerable.GetEnumerator() => GetEnumerator(); + + /// + /// if vector belongs to this space + /// + /// + /// + /// + public static bool operator >(VectorSpace< TVector, TMatrix, TVectorSpace> self, TVector oth) + => self.ContainsVector(oth); + /// + /// + /// + /// + /// + /// + public static bool operator <(VectorSpace< TVector, TMatrix, TVectorSpace> self, TVector oth) + => throw new Exception("TODO"); + /// + /// + /// + /// + /// + /// + public static bool operator >(TVector oth, VectorSpace< TVector, TMatrix, TVectorSpace> self) + => throw new Exception("TODO"); + /// + /// if this space contains the vector + /// + /// + /// + /// + public static bool operator <(TVector oth, VectorSpace< TVector, TMatrix, TVectorSpace> self) + => self.ContainsVector(oth); + /// + /// add a basis to the space + /// + /// + /// + /// + public static TVectorSpace operator +(VectorSpace< TVector, TMatrix, TVectorSpace> self, TVector other) + { + self.AddBasis(other); + return (self as TVectorSpace)!; + } + + /// + /// generate vector space from basis + /// + /// + /// + public static TVectorSpace Dispatch(IEnumerable basis) + { + TVectorSpace res = new(); + res.Fill(basis); + return res; + } +} + +/// +/// meta vector space +/// +/// +/// +/// +/// +public abstract class VectorSpace< TScalar, TVector, TMatrix, TVectorSpace> : + VectorSpace< TVector, TMatrix, TVectorSpace> + where TScalar : notnull + where TVector : Vector, new() + where TMatrix : Matrix, new() + where TVectorSpace : VectorSpace, new() + +{ + + /// + /// constructor + /// + protected VectorSpace() + { + InnerBasis = new List(); + I2OTransformation = new TMatrix(); + O2ITransformation = new TMatrix(); + LTransformation = new TMatrix(); + Basis = new List(); + } + + /// + protected VectorSpace(IEnumerable es) : base(es) { } + + /// + /// check if vector is in the vector space + /// + /// + /// + public override bool ContainsVector(TVector v) + { + if (Basis.Count == Dim) + return true; + if (v.IsEqualApprox(v.VarZero)) + return true; + TVector iv = v.Copy(); + foreach (TVector b in InnerBasis) + { + int pivot = b.LeadingZeros + 1; + iv -= b * iv[pivot] / b[pivot]; + iv[pivot] = default!; + } + + return iv.IsEqualApprox(iv.VarZero); + } + + /// + public sealed override void AddBasis(TVector basis) + { + if (Basis.Count == Dim) + return; + if (ContainsVector(basis)) + return; + InnerBasis = new List(); + Basis.Add(basis); + List basisPaddedWithZero = Basis.Select(x => x.Copy()).ToList(); + while (basisPaddedWithZero.Count < Dim) + basisPaddedWithZero.Add(basis.VarZero); + TMatrix m = Matrix< TVector, TMatrix>.Dispatch(basisPaddedWithZero); + ReducedRowEchelonForm rowReduceInfo = m.RowReduce(); + List fullSpaceBasis = Basis.Select(x => x).ToList(); + foreach (TVector kernelBase in m.KerBasis()) + fullSpaceBasis.Add(kernelBase.Conj()); + if (fullSpaceBasis.Count != Dim) + throw new Exception("Unknown Err"); + TMatrix changeBasis = Matrix< TVector, TMatrix>.Dispatch(fullSpaceBasis); + O2ITransformation = changeBasis.Inv() * rowReduceInfo.ReducedForm; + + LTransformation = rowReduceInfo.RowOperation; + int[] np = Utils.Utils + .XRange(1, Dim) + .Where(x => !rowReduceInfo.Pivots.Contains(x)) + .ToArray(); + int kernelDim = 0; + for (int r = 1; r <= Dim; r++) + if (rowReduceInfo.ReducedForm[r].IsEqualApprox(basis.VarZero)) + { + rowReduceInfo.ReducedForm[r, np[kernelDim]] = Structure.MultiplicationUnit; + kernelDim += 1; + } + else + { + InnerBasis.Add(rowReduceInfo.ReducedForm[r]); + } + + I2OTransformation = rowReduceInfo.ReducedForm.Inv() * m; + } + + internal abstract FieldStructure Structure { get; } + + /// + public sealed override TVectorSpace OrthogonalComplementSpace + { + get + { + List fullSpaceBasis = Basis.Select(x => x).ToList(); + while (fullSpaceBasis.Count < Dim) + fullSpaceBasis.Add(new TVector()); + return Matrix< TVector, TMatrix>.Dispatch(fullSpaceBasis) + .ConjKer(); + } + } +} diff --git a/src/Algebra/Vectors/FVector.cs b/src/Algebra/Vectors/FVector.cs new file mode 100644 index 0000000..e663ea9 --- /dev/null +++ b/src/Algebra/Vectors/FVector.cs @@ -0,0 +1,125 @@ +using System.Numerics; +using Skeleton.Algebra.FixableObjects; +using Skeleton.Algebra.ScalarFieldStructure; +using Skeleton.Algebra.Vectors; + +namespace Skeleton.Algebra; + +public partial class CategoryOf +{ + public static partial class OnField + where TScalar : notnull + { + /// + public class FVector : Vector + { + /// + /// zero vector + /// + public static readonly FVector Zero = new(); + /// + public FVector() { } + + /// + public FVector(params TScalar[] es) : base(es) => this.PostConstructionFix(); + + /// + public FVector(IEnumerable es) : base(es) => this.PostConstructionFix(); + + /// + public override int Dim => FDim; + + /// + /// real part of the vector
+ /// homomorphism from field to complex has to be defined + /// + public OnField.FVector Real => + new(this.Select(x => FieldStructure.ComplexCast(x).Real)); + + /// + /// imaginary part of the vector
+ /// homomorphism from field to complex has to be defined + /// + public OnField.FVector Imaginary => + new(this.Select(x => Structure.ComplexCast(x).Imaginary)); + + /// + /// map the vector to a complex vector
+ /// homomorphism from field to complex has to be defined + /// + /// + public OnField.FVector AsComplex() => + new(this.Select(x => Structure.ComplexCast(x))); + + /// + /// outer product of two vector
+ /// M = v v* + /// + /// + /// + public FMatrix OuterProduct(FVector other) + { + FMatrix res = new(); + for (int i = 1; i <= FDim; i++) + for (int j = 1; j <= FDim; j++) + res[i, j] = Structure.Multiplication(this[i], other[j]); + res.Fix(); + return res; + } + /// + /// cast from vector over K to a vector over R
+ /// morphism K -> R has to be defined + /// + /// + /// + public static explicit operator OnField.FVector(FVector a) => a.Real; + /// + /// cast from vector over K to a vector over C
+ /// morphism K -> C has to be defined + /// + /// + /// + public static explicit operator OnField.FVector(FVector a) => + new(a.Select(x => FieldStructure.ComplexCast(x))); + + /// + /// cast from vector over R to a vector over K
+ /// morphism R -> K has to be defined + /// + /// + /// + public static explicit operator FVector(OnField.FVector a) => + new(a.Select(x => FieldStructure.FromReal(x))); + + /// + /// cast from vector over C to a vector over C
+ /// morphism C -> K has to be defined + /// + /// + /// + public static explicit operator FVector(OnField.FVector a) => + new(a.Select(x => FieldStructure.FromComplex(x))); + + /// + protected override FieldStructure Structure => + FieldStructure; + + /// + public override FVector Conj() + { + TScalar[] ele = new TScalar[FDim]; + for (int i = 1; i <= FDim; i++) + ele[i-1] = FieldStructure.Conj(this[i]); + return new(ele); + } + + protected override TScalar Dot(FVector other) + { + TScalar res = StaticAccess.AdditionUnit; + for (int i = 1; i <= FDim; i++) + res = FieldStructure.Addition(res, FieldStructure.Multiplication(this[i], other[i])); + return res; + } + } + } +} diff --git a/src/Algebra/Vectors/Vector.cs b/src/Algebra/Vectors/Vector.cs new file mode 100644 index 0000000..c7c5ff7 --- /dev/null +++ b/src/Algebra/Vectors/Vector.cs @@ -0,0 +1,499 @@ +using System.Collections; +using System.Numerics; +using Skeleton.Algebra.Extensions; +using Skeleton.Algebra.FixableObjects; +using Skeleton.Algebra.ScalarFieldStructure; + +namespace Skeleton.Algebra.Vectors; + +/// +/// meta for 2-4 dim vectors +/// +public abstract class Vector : FixableTensor + +{ + /// + protected Vector() { } + /// + /// length of the vector + /// + public abstract double Magnitude { get; } + + /// + /// count leading zeros in the vector + /// + public abstract int LeadingZeros { get; } + + /// + /// count ending zeros in the vector + /// + public abstract int EndingZeros { get; } + + /// + /// reset the vector to zero vector + /// + protected abstract void ResetZero(); + + /// + /// fix vector with additional action + /// + /// + protected void WithFix(Action f) + { + f(); + Fix(); + } + + /// + /// + /// + public abstract string Representation { get; } + + /// + /// + /// + public abstract string PythonRepresentation { get; } + + /// + /// + /// + public abstract string CSharpRepresentation { get; } + + /// + /// string representation of the matrix to create a new instance with the same data + /// debug use only + /// + public abstract string LatexRepresentation { get; } +} + +/// +/// meta for 2-4 dim vectors +/// +/// +public abstract class Vector : Vector + where TVector : Vector +{ + /// + /// normalize the vector so that length is 1 + /// + public abstract TVector Normalize { get; } + + /// + /// get a variable zero vector + /// + public abstract TVector VarZero { get; } + + /// + /// get a copy of this vector + /// + /// + public abstract TVector Copy(); + + /// + /// + /// + /// + protected abstract TVector Addition(TVector other); + + /// + /// + /// + protected abstract TVector Negation(); + + /// + /// + /// + /// + /// + public static TVector operator +(Vector< TVector> self, TVector other) => self.Addition(other); + + /// + /// + /// + /// + public static TVector operator -(Vector< TVector> self) => self.Negation(); + + /// + /// OPT Required + /// + public static TVector operator -(Vector< TVector> self, TVector other) + { + return self + -other; + } + + /// + /// the projection of this vector on another vector + /// + /// + /// + public abstract TVector ProjectionOn(TVector b); + + /// + /// the rejection of this vector on another vector + /// + /// + /// + public TVector RejectionOn(TVector b) => this - ProjectionOn(b); + + /// + /// if two vectors are close to each other numerically + /// + /// + /// + public abstract bool IsEqualApprox(TVector other); + + /// + /// Conjugate of this vector + /// + /// + public abstract TVector Conj(); + +} + +/// +/// meta for 2-4 dim vector +/// +/// +/// +public abstract class Vector : Vector, IEnumerable + where TVector : Vector, new() +{ + + /// + protected Vector(IEnumerable elements) + { + Elements = elements.ToArray(); + } + + /// + protected Vector() => Elements = new TScalar[Dim]; + + /// + /// + protected TScalar[] Elements { get; set; } + /// + /// + /// + public TScalar this[int index] + { + get => Elements[index - 1]; + set => Elements[index - 1] = value; + } + + /// + /// + /// + public IEnumerator GetEnumerator() + { + for (int i = 0; i < Dim; i++) + yield return Elements[i]; + } + + IEnumerator IEnumerable.GetEnumerator() => GetEnumerator(); + + /// + protected override void ResetZero() + { + Elements = new TScalar[Dim]; + for (int i = 1; i <= Dim; i++) + Elements[i] = StaticAccess.AdditionUnit; + } + + + /// + /// + /// + /// + protected abstract TScalar Dot(TVector other); + + + /// + /// + /// + /// + /// + public static TVector operator *(Vector< TScalar, TVector> self, TScalar other) + { + TVector res = new (); + for (int i = 1; i <= self.Dim; i++) + res[i] = res.Structure.Multiplication(self[i], other); + return res; + } + + /// + /// + /// + /// + /// + public static TVector operator *(TScalar other, Vector< TScalar, TVector> self) => self * other!; + + /// + /// + /// + /// + /// + public static TVector operator /(Vector< TScalar, TVector> self, TScalar other) + { + TVector res = new TVector(); + for (int i = 1; i <= self.Dim; i++) + res[i] = res.Structure.Division(self[i], other); + return res; + } + + /// + /// + /// + /// + /// + /// + public static TVector operator /(Vector< TScalar, TVector> self, double other) + => self / self.Structure.FromReal(other)!; + /// + /// + /// + /// + /// + /// + public static TVector operator *(Vector< TScalar, TVector> self, double other) + => self * self.Structure.FromReal(other)!; + /// + /// + /// + /// + /// + public static TScalar operator *(Vector< TScalar, TVector> self, TVector other) => self.Dot(other); + + /// + /// if this vector is linear dependent with another vector v + /// + /// + /// + public bool IsLinearDependantWith(TVector v) + { + bool started = false; + TScalar comDiv = StaticAccess.AdditionUnit; + for (int i = 1; i <= v.Dim; i++) + { + if ( + Structure.IsEqualApprox(v[i], StaticAccess.AdditionUnit) && + Structure.IsEqualApprox(this[i], StaticAccess.AdditionUnit) + ) + continue; + if (Structure.IsEqualApprox(v[i], StaticAccess.AdditionUnit)) + return false; + TScalar ratio = Structure.Division(this[i], v[i]); + if (!started) + comDiv = ratio; + else if (!Structure.IsEqualApprox(comDiv, ratio)) + return false; + } + + return true; + } + + /// + /// field structure + /// + protected abstract FieldStructure Structure { get; } + + /// + public sealed override bool IsEqualApprox(TVector other) + { + for (int i = 1; i <= Dim; i++) + if (!Structure.IsEqualApprox(this[i], other[i])) + return false; + return true; + } + /// + protected sealed override TVector Negation() + { + TVector res = new TVector(); + for (int i = 1; i <= Dim; i++) + res[i] = Structure.AdditionInverse(this[i]); + return res; + } + + /// + protected sealed override TVector Addition(TVector other) + { + TVector res = new(); + for (int i = 1; i <= Dim; i++) + res[i] = Structure.Addition(this[i], other[i]); + return res; + } + + /// + public sealed override TVector Copy() + { + TVector res = new(); + for (int i = 1; i <= Dim; i++) + res[i] = this[i]; + return res; + } + + internal override void Fix() + { + for (int i = 1; i <= Dim; i++) + this[i] = Structure.Fix(this[i]); + base.Fix(); + } + + /// + public sealed override int EndingZeros + { + get + { + for (int i = 0; i <= Dim; i++) + if (!Structure.IsEqualApprox(Elements[^(i + 1)], StaticAccess.AdditionUnit)) + return i; + return Dim; + } + } + + /// + public sealed override int LeadingZeros { + get + { + for(int i = 0;i.AdditionUnit)) + return i; + return Dim; + } + } + + /*/// + public override TVector Conj() + { + TVector res = new(); + for (int i = 1; i <= Dim; i++) + res[i] = Structure.Conj(this[i]); + return res; + }*/ + + /// + public sealed override TVector ProjectionOn(TVector b) + { + TVector bconj = b.Conj(); + return Structure.Division(this * bconj, b * bconj)! * b; + } + /// + public sealed override double Magnitude + => Math.Sqrt(Elements.Sum(x => Structure.ConjX2(x))); + + /// + public sealed override TVector Normalize => ( + Magnitude.IsEqualApprox(0) ? this as TVector : this / Magnitude + )!; + + /// + public override string Representation => + string.Join(",\t", Elements.Select(x => Structure.Representation(x))); + + /// + public override string PythonRepresentation + => "np.matrix([" + string.Join(",", Elements.Select(x => Structure.PythonRepresentation(x))) + "])"; + + /// + public override string CSharpRepresentation + => $"new {GetType().Name} (" + + string.Join(",", Elements.Select(x => Structure.CSharpRepresentation(x))) + + ");"; + + /// + public override string LatexRepresentation + => "\\begin{pmatrix}" + + string.Join('&', Elements.Select(x => Structure.Representation(x))) + + "\\end{pmatrix}"; + + /// + /// generate vector from elements + /// + /// + /// + public static TVector Dispatch(IEnumerable elements) + { + TVector res = new () + { + Elements = elements.ToArray() + }; + res.PostConstructionFix(); + return res; + } + + /// + /// save object to string + /// + public string DumpString + => "" + + string.Join("", Elements.Select(x => Structure.DumpString(x))) + + ""; + + /// + /// + /// + /// + /// + public static TVector Resolve(string dumpString) + { + TVector res = new (); + res.Elements = dumpString + .Replace("", "") + .Replace("", "") + .Split("") + .Select(x => FieldStructure.Dispatch().Resolve(x)) + .ToArray(); + return res; + } + + + /// + /// + /// + public double MaxError => this.Select(x => Structure.MaxError(x)).Max(); + + internal TComplexVector ComplexCast() + where TComplexVector : Vector, new() + => Vector< Complex, TComplexVector>.Dispatch(Elements.Select(Structure.ComplexCast)); + + /// + public sealed override TVector VarZero + { + get + { + TVector res = new(); + res.Elements = Enumerable.Range(0, Dim).Select(_ => StaticAccess.AdditionUnit).ToArray(); + return res; + } + } + + + /// + /// make the vector a (dim - d) vector + /// V[i] = 0 for all i<=d + /// + /// + /// + public TVector UpTruncateBy(int d) + { + TVector res = Copy(); + for (int i = 1; i <= d; i++) + res[i] = StaticAccess.AdditionUnit; + return res; + } + + /// + /// make the vector a d-1 dimensional vector + /// V[i] = 0 for all i>=d + /// + /// + /// + public TVector DownTruncateBy(int d) + { + TVector res = Copy(); + for (int i = d; i < Dim; i++) + res[i] = StaticAccess.AdditionUnit; + return res; + } + +} + diff --git a/src/Analysis/AnalyticFunctions/Arithmetics/AnalyticAddition.cs b/src/Analysis/AnalyticFunctions/Arithmetics/AnalyticAddition.cs new file mode 100644 index 0000000..c674844 --- /dev/null +++ b/src/Analysis/AnalyticFunctions/Arithmetics/AnalyticAddition.cs @@ -0,0 +1,39 @@ +using System.Numerics; +using Skeleton.Analysis.AnalyticFunctions.Implements; + +namespace Skeleton.Analysis.AnalyticFunctions.Arithmetics; + +/// +/// addition +/// +public class AnalyticAddition : Analytic +{ + /// + /// constructor + /// + /// + /// + public AnalyticAddition(Analytic left, Analytic right) + { + Left = left; + Right = right; + } + + private Analytic Left { get; } + private Analytic Right { get; } + + /// + /// d (f+g)/dx + /// + public override Analytic Derivative => Left.Derivative + Right.Derivative; + + /// + /// [f+g](x) + /// + /// + /// + public override Complex Evaluate(Complex x) + { + return Left.Evaluate(x) + Right.Evaluate(x); + } +} \ No newline at end of file diff --git a/src/Analysis/AnalyticFunctions/Arithmetics/AnalyticApply.cs b/src/Analysis/AnalyticFunctions/Arithmetics/AnalyticApply.cs new file mode 100644 index 0000000..e61a9ec --- /dev/null +++ b/src/Analysis/AnalyticFunctions/Arithmetics/AnalyticApply.cs @@ -0,0 +1,39 @@ +using System.Numerics; +using Skeleton.Analysis.AnalyticFunctions.Implements; + +namespace Skeleton.Analysis.AnalyticFunctions.Arithmetics; + +/// +/// f(x),g(x) -> f(g(x)) +/// +public class AnalyticApply : Analytic +{ + /// + /// Constructor + /// + /// + /// + public AnalyticApply(Analytic outer, Analytic inner) + { + Outer = outer; + Inner = inner; + } + + private Analytic Outer { get; } + private Analytic Inner { get; } + + /// + /// d(f(g(x)))/dx + /// + public override Analytic Derivative => Inner.Derivative * Outer.Derivative[Inner]; + + /// + /// f(g(x)) + /// + /// + /// + public override Complex Evaluate(Complex x) + { + return Outer.Evaluate(Inner.Evaluate(x)); + } +} \ No newline at end of file diff --git a/src/Analysis/AnalyticFunctions/Arithmetics/AnalyticConstant.cs b/src/Analysis/AnalyticFunctions/Arithmetics/AnalyticConstant.cs new file mode 100644 index 0000000..34c92f8 --- /dev/null +++ b/src/Analysis/AnalyticFunctions/Arithmetics/AnalyticConstant.cs @@ -0,0 +1,46 @@ +using System.Numerics; +using Skeleton.Analysis.AnalyticFunctions.Implements; + +namespace Skeleton.Analysis.AnalyticFunctions.Arithmetics; + +/// +/// C +/// +public class AnalyticConstant : Analytic +{ + /// + /// zero + /// + public static readonly AnalyticConstant Zero = new(0); + + /// + /// one + /// + public static readonly AnalyticConstant One = new(1); + + /// + /// constructor + /// + /// + public AnalyticConstant(Complex c) + { + Constant = c; + } + + private Complex Constant { get; } + + /// + /// dc/dx = 0 + /// + public override Analytic Derivative => Zero; + + /// + /// c(x) = c + /// + /// + /// + public override Complex Evaluate(Complex x) + { + return Constant; + } +} \ No newline at end of file diff --git a/src/Analysis/AnalyticFunctions/Arithmetics/AnalyticIdentity.cs b/src/Analysis/AnalyticFunctions/Arithmetics/AnalyticIdentity.cs new file mode 100644 index 0000000..eaf75c2 --- /dev/null +++ b/src/Analysis/AnalyticFunctions/Arithmetics/AnalyticIdentity.cs @@ -0,0 +1,25 @@ +using System.Numerics; +using Skeleton.Analysis.AnalyticFunctions.Implements; + +namespace Skeleton.Analysis.AnalyticFunctions.Arithmetics; + +/// +/// f(x) = x +/// +public class AnalyticIdentity : Analytic +{ + /// + /// dx/dx = 1 + /// + public override Analytic Derivative => AnalyticConstant.One; + + /// + /// x = x + /// + /// + /// + public override Complex Evaluate(Complex x) + { + return x; + } +} \ No newline at end of file diff --git a/src/Analysis/AnalyticFunctions/Arithmetics/AnalyticNegation.cs b/src/Analysis/AnalyticFunctions/Arithmetics/AnalyticNegation.cs new file mode 100644 index 0000000..465cd4e --- /dev/null +++ b/src/Analysis/AnalyticFunctions/Arithmetics/AnalyticNegation.cs @@ -0,0 +1,36 @@ +using System.Numerics; +using Skeleton.Analysis.AnalyticFunctions.Implements; + +namespace Skeleton.Analysis.AnalyticFunctions.Arithmetics; + +/// +/// -f +/// +public class AnalyticNegation : Analytic +{ + /// + /// constructor + /// + /// + public AnalyticNegation(Analytic element) + { + Element = element; + } + + private Analytic Element { get; } + + /// + /// d[-f]/dx = -df/dx + /// + public override Analytic Derivative => -Element.Derivative; + + /// + /// [-f](x) = -x + /// + /// + /// + public override Complex Evaluate(Complex x) + { + return -Element.Evaluate(x); + } +} \ No newline at end of file diff --git a/src/Analysis/AnalyticFunctions/Arithmetics/AnalyticProduction.cs b/src/Analysis/AnalyticFunctions/Arithmetics/AnalyticProduction.cs new file mode 100644 index 0000000..4afceb4 --- /dev/null +++ b/src/Analysis/AnalyticFunctions/Arithmetics/AnalyticProduction.cs @@ -0,0 +1,39 @@ +using System.Numerics; +using Skeleton.Analysis.AnalyticFunctions.Implements; + +namespace Skeleton.Analysis.AnalyticFunctions.Arithmetics; + +/// +/// f(x), g(x) -> f(x)g(x) +/// +public class AnalyticProduction : Analytic +{ + /// + /// Constructor + /// + /// + /// + public AnalyticProduction(Analytic left, Analytic right) + { + Left = left; + Right = right; + } + + private Analytic Left { get; } + private Analytic Right { get; } + + /// + /// dfg/dx = fdg/dx + gdf/dx + /// + public override Analytic Derivative => Left.Derivative * Right + Left * Right.Derivative; + + /// + /// f(x)*g(x) + /// + /// + /// + public override Complex Evaluate(Complex x) + { + return Left.Evaluate(x) * Right.Evaluate(x); + } +} \ No newline at end of file diff --git a/src/Analysis/AnalyticFunctions/IAnalytic.cs b/src/Analysis/AnalyticFunctions/IAnalytic.cs new file mode 100644 index 0000000..bcca257 --- /dev/null +++ b/src/Analysis/AnalyticFunctions/IAnalytic.cs @@ -0,0 +1,19 @@ +using System.Numerics; +using Skeleton.Analysis.AnalyticFunctions.Implements; + +namespace Skeleton.Analysis.AnalyticFunctions; + +/// +/// +public interface IAnalytic +{ + /// + /// + Analytic Derivative { get; } + + /// + /// + /// + /// + Complex Evaluate(Complex x); +} \ No newline at end of file diff --git a/src/Analysis/AnalyticFunctions/Implements/Analytic.cs b/src/Analysis/AnalyticFunctions/Implements/Analytic.cs new file mode 100644 index 0000000..6ce7bd5 --- /dev/null +++ b/src/Analysis/AnalyticFunctions/Implements/Analytic.cs @@ -0,0 +1,63 @@ +using System.Numerics; +using Skeleton.Analysis.AnalyticFunctions.Arithmetics; + +namespace Skeleton.Analysis.AnalyticFunctions.Implements; + +/// +/// analytic functions +/// +public abstract class Analytic : IAnalytic +{ + /// + /// f(g)(x) -> f(g(x)) + /// + /// + public Analytic this[Analytic inner] => new AnalyticApply(this, inner); + + /// + /// df/dx + /// + public abstract Analytic Derivative { get; } + + /// + /// f(x) + /// + /// + /// + public abstract Complex Evaluate(Complex x); + + /// + /// [f+g](x) = f(x) + g(x) + /// + /// + /// + /// + public static Analytic operator +(Analytic a, Analytic b) + => new AnalyticAddition(a, b); + + /// + /// [-f](x) = -f(x) + /// + /// + /// + public static Analytic operator -(Analytic a) + => new AnalyticNegation(a); + + /// + /// [f-g](x) = f(x) - g(x) + /// + /// + /// + /// + public static Analytic operator -(Analytic a, Analytic b) + => -b + a; + + /// + /// [fg](x) = f(x)g(x) + /// + /// + /// + /// + public static Analytic operator *(Analytic a, Analytic b) + => new AnalyticProduction(a, b); +} \ No newline at end of file diff --git a/src/Analysis/AnalyticFunctions/Polynomials/IPolynomial.cs b/src/Analysis/AnalyticFunctions/Polynomials/IPolynomial.cs new file mode 100644 index 0000000..0474da9 --- /dev/null +++ b/src/Analysis/AnalyticFunctions/Polynomials/IPolynomial.cs @@ -0,0 +1,29 @@ +using System.Numerics; +using Skeleton.Analysis.BivariantFunctions.Real.Implements; + +namespace Skeleton.Analysis.AnalyticFunctions.Polynomials; + +/// +/// +public interface IPolynomial : IAnalytic +{ + /// + /// + int Dim { get; } + + /// + /// + int Degree { get; } + + /// + /// + Complex[] Coefficients { get; set; } + + /// + /// + BivariantContinuous Re { get; } + + /// + /// + BivariantContinuous Im { get; } +} \ No newline at end of file diff --git a/src/Analysis/AnalyticFunctions/Polynomials/Implements/C0Polynomial.cs b/src/Analysis/AnalyticFunctions/Polynomials/Implements/C0Polynomial.cs new file mode 100644 index 0000000..92824bd --- /dev/null +++ b/src/Analysis/AnalyticFunctions/Polynomials/Implements/C0Polynomial.cs @@ -0,0 +1,61 @@ +using System.Numerics; +using Skeleton.Analysis.BivariantFunctions.Real.Arithmetics; +using Skeleton.Analysis.BivariantFunctions.Real.Implements; + +namespace Skeleton.Analysis.AnalyticFunctions.Polynomials.Implements; + +/// +/// p = a0 +/// +public class C0Polynomial : Polynomial +{ + /// + /// constructor + /// + /// + public C0Polynomial(Complex a0) + { + Coefficients = new[] { a0 }; + Fix(); + } + + /// + /// a0 + /// + public Complex A0 => Coefficients[0]; + + /// + /// 0 + /// + public override Polynomial Derivative => Zero; + + /// + /// real part + /// + public override BivariantContinuous Re => new BivariantContinuousConstant(A0.Real); + + /// + /// imaginary part + /// + public override BivariantContinuous Im => new BivariantContinuousConstant(A0.Imaginary); + + /// + /// string for debug + /// + public override string Representation => $"{A0.Real} + {A0.Imaginary}i"; + + /// + public override int Dim => 0; + + /// + public override IEnumerable Roots => Array.Empty(); + + /// + /// + /// + /// + public override Complex Evaluate(Complex x) + { + return A0; + } +} \ No newline at end of file diff --git a/src/Analysis/AnalyticFunctions/Polynomials/Implements/C1Polynomial.cs b/src/Analysis/AnalyticFunctions/Polynomials/Implements/C1Polynomial.cs new file mode 100644 index 0000000..127e8e7 --- /dev/null +++ b/src/Analysis/AnalyticFunctions/Polynomials/Implements/C1Polynomial.cs @@ -0,0 +1,54 @@ +using System.Numerics; +using Skeleton.Analysis.BivariantFunctions.Real.Implements; + +namespace Skeleton.Analysis.AnalyticFunctions.Polynomials.Implements; + +/// +/// p(x) = a0 + a1 x +/// +public class C1Polynomial : Polynomial +{ + private static readonly BivariantContinuous X = BivariantContinuous.X; + private static readonly BivariantContinuous Y = BivariantContinuous.Y; + + /// + public C1Polynomial() + { + Coefficients = new Complex[2]; + Fix(); + } + + /// + public C1Polynomial(Complex a0, Complex a1) + { + Coefficients = new[] { a0, a1 }; + Fix(); + } + + private Complex A0 => Coefficients[0]; + private Complex A1 => Coefficients[1]; + + /// + public override string Representation => $"{A0.Real} + {A0.Imaginary}i + ({A1.Real} + {A1.Imaginary}i)x"; + + /// + public override int Dim => 1; + + /// + public override Polynomial Derivative => new C0Polynomial(A1); + + /// + public override IEnumerable Roots => new[] { -A0 / A1 }; + + /// + public override BivariantContinuous Re => A0.Real + A1.Real * X - A1.Imaginary * Y; + + /// + public override BivariantContinuous Im => A0.Imaginary + A1.Real * Y + A1.Imaginary * X; + + /// + public override Complex Evaluate(Complex x) + { + return A0 + x * A1; + } +} \ No newline at end of file diff --git a/src/Analysis/AnalyticFunctions/Polynomials/Implements/C2Polynomial.cs b/src/Analysis/AnalyticFunctions/Polynomials/Implements/C2Polynomial.cs new file mode 100644 index 0000000..6b5633e --- /dev/null +++ b/src/Analysis/AnalyticFunctions/Polynomials/Implements/C2Polynomial.cs @@ -0,0 +1,74 @@ +using System.Numerics; +using Skeleton.Algebra.Extensions; +using Skeleton.Analysis.BivariantFunctions.Real.Implements; +using Skeleton.Analysis.BivariantFunctions.Real.PolynomialHelpers; +using Skeleton.Utils.Helpers; + +namespace Skeleton.Analysis.AnalyticFunctions.Polynomials.Implements; + +/// +/// +public class C2Polynomial : Polynomial +{ + /// + public C2Polynomial() + { + Coefficients = new Complex[3]; + Fix(); + } + + /// + public C2Polynomial(Complex a0, Complex a1, Complex a2) + { + Coefficients = new[] { a0, a1, a2 }; + Fix(); + } + + /// + /// + public Complex A0 => Coefficients[0]; + + /// + /// + public Complex A1 => Coefficients[1]; + + /// + /// + public Complex A2 => Coefficients[2]; + + /// + public override string Representation => + $"{A0.Real} + {A0.Imaginary}i + ({A1.Real} + {A1.Imaginary}i)x + ({A2.Real} + {A2.Imaginary}i)x2"; + + /// + public override int Dim => 2; + + /// + public override Polynomial Derivative => new C1Polynomial(A1, A2 * 2); + + /// + public override IEnumerable Roots => Decorators.WithTol(1E-20, + () => + { + if (A2.IsEqualApprox(0)) + return new C1Polynomial(A0, A1).Roots; + if (A0.IsEqualApprox(0)) + return new C1Polynomial(A1, A2).Roots.Concat(new[] { Complex.Zero }); + Complex x1 = (-A1 + Complex.Sqrt(A1 * A1 - 4 * A2 * A0)) / (2 * A2); + Complex x2 = (-A1 - Complex.Sqrt(A1 * A1 - 4 * A2 * A0)) / (2 * A2); + return new[] { x1, x2 }; + } + ); + + /// + public override BivariantContinuous Re => new PolyCHelper.PolyC2Re(A0, A1, A2); + + /// + public override BivariantContinuous Im => new PolyCHelper.PolyC2Im(A0, A1, A2); + + /// + public override Complex Evaluate(Complex x) + { + return A0 + x * (A1 + x * A2); + } +} \ No newline at end of file diff --git a/src/Analysis/AnalyticFunctions/Polynomials/Implements/C3Polynomial.cs b/src/Analysis/AnalyticFunctions/Polynomials/Implements/C3Polynomial.cs new file mode 100644 index 0000000..132c36c --- /dev/null +++ b/src/Analysis/AnalyticFunctions/Polynomials/Implements/C3Polynomial.cs @@ -0,0 +1,105 @@ +using System.Numerics; +using Skeleton.Algebra.Extensions; +using Skeleton.Analysis.BivariantFunctions.Real.Implements; +using Skeleton.Analysis.BivariantFunctions.Real.PolynomialHelpers; +using Skeleton.Utils.Helpers; + +namespace Skeleton.Analysis.AnalyticFunctions.Polynomials.Implements; + +/// +/// +public class C3Polynomial : Polynomial +{ + /// + public C3Polynomial() + { + Coefficients = new Complex[4]; + Fix(); + } + + /// + public C3Polynomial(Complex a0, Complex a1, Complex a2, Complex a3) + { + Coefficients = new[] { a0, a1, a2, a3 }; + Fix(); + } + + /// + /// + public Complex A0 => Coefficients[0]; + + /// + /// + public Complex A1 => Coefficients[1]; + + /// + /// + public Complex A2 => Coefficients[2]; + + /// + /// + public Complex A3 => Coefficients[3]; + + /// + public override int Dim => 3; + + /// + public override Polynomial Derivative => new C2Polynomial(A1, A2 * 2, A3 * 3); + + /// + public override IEnumerable Roots => + Decorators.WithTol(1E-20, () => + { + if (A3.IsEqualApprox(0)) + return new C2Polynomial(A0, A1, A2).Roots; + if (A0.IsEqualApprox(0)) + return new C2Polynomial(A1, A2, A3).Roots.Concat(new[] { Complex.Zero }); + Polynomial gcd = GreatestCommonDivisor(Derivative); + if (!gcd.IsConstant) + { + Polynomial div = LongDivision(gcd).Quotient; + return div.Roots.Concat(gcd.Roots); + } + + Complex a = A3; + Complex b = A2; + Complex c = A1; + Complex d = A0; + Complex D0 = b * b - 3d * a * c; + Complex D1 = 2d * b * b * b - 9d * a * b * c + 27d * a * a * d; + if (D0.IsEqualApprox(0) && D1.IsEqualApprox(0)) + { + Complex root = -b / (a * 3d); + return new[] { root, root, root }; + } + + Complex C = Complex.Pow((D1 + Complex.Sqrt(D1 * D1 - 4d * D0 * D0 * D0)) / 2d, 1d / 3d); + if (C.IsEqualApprox(0)) + C = Complex.Pow((D1 - Complex.Sqrt(D1 * D1 - 4d * D0 * D0 * D0)) / 2d, 1d / 3d); + if (C.IsEqualApprox(0)) + throw new Exception($"Unknown situation with {a}, {b}, {c}, {d}"); //TODO + Complex ee = (-1d + Complex.Sqrt(-3d)) / 2d; + return new[] + { + -1d / (3d * a) * (b + C + D0 / C), + -1d / (3d * a) * (b + ee * C + D0 / (ee * C)), + -1d / (3d * a) * (b + ee * ee * C + D0 / (ee * ee * C)) + }; + }); + + /// + public override BivariantContinuous Re => new PolyCHelper.PolyC3Re(A0, A1, A2, A3); + + /// + public override BivariantContinuous Im => new PolyCHelper.PolyC3Im(A0, A1, A2, A3); + + /// + public override string Representation => + $"{A0.Real} + {A0.Imaginary}i + ({A1.Real} + {A1.Imaginary}i)x + ({A2.Real} + {A2.Imaginary}i)x2 + ({A3.Real} + {A3.Imaginary}i)X3"; + + /// + public override Complex Evaluate(Complex x) + { + return A0 + x * (A1 + x * (A2 + x * A3)); + } +} \ No newline at end of file diff --git a/src/Analysis/AnalyticFunctions/Polynomials/Implements/C4Polynomial.cs b/src/Analysis/AnalyticFunctions/Polynomials/Implements/C4Polynomial.cs new file mode 100644 index 0000000..3a427ef --- /dev/null +++ b/src/Analysis/AnalyticFunctions/Polynomials/Implements/C4Polynomial.cs @@ -0,0 +1,140 @@ +using System.Numerics; +using Skeleton.Algebra.Extensions; +using Skeleton.Analysis.BivariantFunctions.Real.Implements; +using Skeleton.Analysis.BivariantFunctions.Real.PolynomialHelpers; +using Skeleton.Utils.Helpers; + +namespace Skeleton.Analysis.AnalyticFunctions.Polynomials.Implements; + +/// +/// +public class C4Polynomial : Polynomial +{ + /// + public C4Polynomial() + { + Coefficients = new Complex[5]; + Fix(); + } + + /// + public C4Polynomial(Complex a0, Complex a1, Complex a2, Complex a3, Complex a4) + { + Coefficients = new[] { a0, a1, a2, a3, a4 }; + Fix(); + } + + /// + /// + public Complex A0 => Coefficients[0]; + + /// + /// + public Complex A1 => Coefficients[1]; + + /// + /// + public Complex A2 => Coefficients[2]; + + /// + /// + public Complex A3 => Coefficients[3]; + + /// + /// + public Complex A4 => Coefficients[4]; + + /// + public override int Dim => 4; + + /// + public override Polynomial Derivative => new C3Polynomial(A1, A2 * 2, A3 * 3, A4 * 4); + + /// + public override IEnumerable Roots => + Decorators.WithTol(1e-20, () => + { + if (A4.IsEqualApprox(Complex.Zero)) + return new C3Polynomial(A0, A1, A2, A3).Roots; + if (A0.IsEqualApprox(Complex.Zero)) + return new C3Polynomial(A1, A2, A3, A4).Roots.Concat(new[] { Complex.Zero }); + if (A3.IsEqualApprox(Complex.Zero) && A1.IsEqualApprox(Complex.Zero)) + return new C2Polynomial(A0, A2, A4).Roots + .SelectMany(x => new[] { Complex.Sqrt(x), -Complex.Sqrt(x) }); + + + Polynomial gcd = GreatestCommonDivisor(Derivative); + if (!gcd.IsConstant) + { + Polynomial div = LongDivision(gcd).Quotient; + return div.Roots.Concat(gcd.Roots); + } + + Complex a = A4; + Complex b = A3; + Complex c = A2; + Complex d = A1; + Complex e = A0; + Complex P = (c * c + 12d * a * e - 3d * b * d) / 9d; + Complex Q = (27d * a * d * d + 2d * c * c * c + 27d * b * b * e - 72d * a * c * e - 9d * b * c * d) / 54d; + Complex D = Complex.Sqrt(Q * Q - P * P * P); + Complex ua = Complex.Pow(Q + D, 1d / 3d); + Complex ub = Complex.Pow(Q - D, 1d / 3d); + Complex u = Complex.Abs(ua) > Complex.Abs(ub) ? ua : ub; + Complex v = u.IsEqualApprox(0d) ? 0d : P / u; + Complex omega = -0.5d + Complex.Sqrt(3d) * Complex.ImaginaryOne / 2d; + + Complex m(int k) + { + return Complex.Sqrt(b * b - 8d / 3d * a * c + + 4d * a * (Utils.Utils.IterPow(omega, k - 1) * u + + Utils.Utils.IterPow(omega, 4 - k) * v)); + } + + Complex S(int k) + { + return 2d * b * b - 16d / 3d * a * c - + 4d * a * (Utils.Utils.IterPow(omega, k - 1) * u + Utils.Utils.IterPow(omega, 4 - k) * v); + } + + int i = new[] { 1, 2, 3 }.MaxBy(x => m(x).Magnitude); + Complex mm = m(i); + Complex SS = S(i); + Complex T = 0; + if (mm.Magnitude.IsEqualApprox(0)) + { + mm = 0; + SS = b * b - 8d / 3d * a * c; + T = 0; + } + else + { + T = (8d * a * b * c - 16d * a * a * d - 2d * b * b * b) / mm; + } + + Complex x(int i) + { + return (-b + Utils.Utils.SPow(Math.Ceiling(i / 2f)) * mm + Utils.Utils.SPow(i + 1) * + Complex.Sqrt(SS + Utils.Utils.SPow(Math.Ceiling(i / 2f)) * T)) / + (4d * a); + } + + return new[] { x(1), x(2), x(3), x(4) }; + }); + + + /// + public override BivariantContinuous Re => new PolyCHelper.PolyC4Re(A0, A1, A2, A3, A4); + + /// + public override BivariantContinuous Im => new PolyCHelper.PolyC4Im(A0, A1, A2, A3, A4); + + /// + public override string Representation => + $"{A0.Real} + {A0.Imaginary}i + ({A1.Real} + {A1.Imaginary}i)x + ({A2.Real} + {A2.Imaginary}i)x2 + ({A3.Real} + {A3.Imaginary}i)X3 + ({A4.Real} + {A4.Imaginary}i)x4"; + + /// + public override Complex Evaluate(Complex x) + => A0 + x * (A1 + x * (A2 + x * (A3 + x * A4))); + +} \ No newline at end of file diff --git a/src/Analysis/AnalyticFunctions/Polynomials/Implements/GeneralPolynomial.cs b/src/Analysis/AnalyticFunctions/Polynomials/Implements/GeneralPolynomial.cs new file mode 100644 index 0000000..b5ae451 --- /dev/null +++ b/src/Analysis/AnalyticFunctions/Polynomials/Implements/GeneralPolynomial.cs @@ -0,0 +1,29 @@ +using System.Numerics; +using Skeleton.Analysis.BivariantFunctions.Real.Implements; + +namespace Skeleton.Analysis.AnalyticFunctions.Polynomials.Implements; + +/// +public class GeneralPolynomial : Polynomial +{ + /// + public GeneralPolynomial(int dim, Complex[] coefficients) + { + PDim = dim; + Coefficients = new Complex[dim + 1]; + for (int i = 0; i <= Dim; i++) + Coefficients[i] = coefficients[i]; + } + + private int PDim { get; } + public override IEnumerable Roots { get; } + + /// + public override int Dim => PDim; + + /// + public override BivariantContinuous Re { get; } + + /// + public override BivariantContinuous Im { get; } +} \ No newline at end of file diff --git a/src/Analysis/AnalyticFunctions/Polynomials/Implements/Polynomial.cs b/src/Analysis/AnalyticFunctions/Polynomials/Implements/Polynomial.cs new file mode 100644 index 0000000..9dd3538 --- /dev/null +++ b/src/Analysis/AnalyticFunctions/Polynomials/Implements/Polynomial.cs @@ -0,0 +1,262 @@ +using System.Numerics; +using Skeleton.Algebra.Extensions; +using Skeleton.Algebra.ScalarFieldStructure; +using Skeleton.Analysis.AnalyticFunctions.Implements; +using Skeleton.Analysis.BivariantFunctions.Real.Implements; +using Skeleton.DataStructure.Packs; +using Skeleton.Utils.Helpers; + +namespace Skeleton.Analysis.AnalyticFunctions.Polynomials.Implements; + +/// +/// +public abstract class Polynomial : Analytic, IPolynomial +{ + /// + /// + protected static readonly Polynomial Zero = new C0Polynomial(0); + + private static readonly Polynomial One = new C0Polynomial(1); + private static readonly Polynomial X = new C1Polynomial(0, 1); + private static readonly Polynomial X2 = new C2Polynomial(0, 0, 1); + private static readonly Polynomial X3 = new C3Polynomial(0, 0, 0, 1); + private static readonly Polynomial X4 = new C4Polynomial(0, 0, 0, 0, 1); + private bool IsZero => Decorators.WithTol(1e-30, () => Coefficients.All(x => x.IsEqualApprox(Complex.Zero))); + + /// + /// + public bool IsConstant => + Decorators.WithTol(1E-30, () => Coefficients.Skip(1).All(x => x.IsEqualApprox(Complex.Zero))); + + /// + public override Polynomial Derivative + { + get + { + Complex[] coeff = new Complex[Dim]; + for (int i = 1; i <= Dim; i++) + coeff[i - 1] = Coefficients[i] * i; + return Dispatch(coeff); + } + } + + /// + /// + public virtual string Representation => + string.Join( + " + ", + Coefficients + .Zip(Enumerable.Range(0, Degree + 1)) + .Select( + x + => "(" + + ComplexFieldStructure.Structure.Representation(x.First) + + $")x^{x.Second}" + ) + ); + + /// + /// + public abstract IEnumerable Roots { get; } + + /// + public abstract int Dim { get; } + + /// + public int Degree + { + get + { + for (int i = 0; i <= Dim; i++) + if (Coefficients[^(i + 1)] != Complex.Zero) + return Dim - i; + return 0; + } + } + + /// + public Complex[] Coefficients { get; set; } = Array.Empty(); + + /// + public abstract BivariantContinuous Re { get; } + + /// + public abstract BivariantContinuous Im { get; } + + /// + /// + /// + /// + /// + public static Polynomial Dispatch(Complex[] coefficients) + { + return coefficients.Length switch + { + 1 => new C0Polynomial(coefficients[0]), + 2 => new C1Polynomial(coefficients[0], coefficients[1]), + 3 => new C2Polynomial(coefficients[0], coefficients[1], coefficients[2]), + 4 => new C3Polynomial(coefficients[0], coefficients[1], coefficients[2], coefficients[3]), + 5 => new C4Polynomial(coefficients[0], coefficients[1], coefficients[2], coefficients[3], coefficients[4]), + _ => new GeneralPolynomial(coefficients.Length - 1, coefficients) + }; + } + + /// + /// + /// + /// + /// + public static Polynomial operator +(Polynomial a, Complex b) => a + new C0Polynomial(b); + + /// + /// + /// + /// + /// + public static Polynomial operator +(Complex b, Polynomial a) => a + b; + + /// + /// + /// + /// + /// + public static Polynomial operator -(Polynomial a, Complex b) => a + new C0Polynomial(-b); + + /// + /// + /// + /// + /// + public static Polynomial operator +(Polynomial a, Polynomial b) + { + int d = Math.Max(a.Degree, b.Degree); + int md = Math.Min(a.Degree, b.Degree); + Polynomial lp = new[] { a, b }.MaxBy(p => p.Degree)!; + Complex[] coefficients = new Complex[d + 1]; + for (int i = 0; i <= d; i++) + if (i <= md) + coefficients[i] = a.Coefficients[i] + b.Coefficients[i]; + else + coefficients[i] = lp.Coefficients[i]; + return Dispatch(coefficients); + } + + /// + /// + /// + /// + public static Polynomial operator -(Polynomial a) + => Dispatch(a.Coefficients.Select(x => -x).ToArray()); + + /// + /// + /// + /// + /// + public static Polynomial operator -(Polynomial a, Polynomial b) => -b + a; + + /// + /// + /// + /// + /// + public static Polynomial operator *(Polynomial a, Complex b) + => Dispatch(a.Coefficients.Select(x => x * b).ToArray()); + + /// + /// + /// + /// + /// + public static Polynomial operator *(Complex b, Polynomial a) => a * b; + + /// + /// + /// + /// + /// + public static Polynomial operator /(Polynomial a, Complex b) + => Dispatch(a.Coefficients.Select(x => x / b).ToArray()); + + /// + /// + /// + /// + /// + /// + public static Polynomial operator *(Polynomial a, Polynomial b) + { + int cd = a.Degree + b.Degree; + if (cd > 4) + throw new Exception("TODO"); + Complex[] coefficients = new Complex[cd + 1]; + for (int i = 0; i <= a.Degree; i++) + for (int j = 0; j <= b.Degree; j++) + coefficients[i + j] += a.Coefficients[i] * b.Coefficients[j]; + return Dispatch(coefficients); + } + + /// + /// + /// + /// + public PolynomialDivisionResult LongDivision(Polynomial a) + { + if (a.Degree > Degree) + return new PolynomialDivisionResult(Zero, a); + Complex leadingC = Coefficients[Degree] / a.Coefficients[a.Degree]; + int degreeDifference = Degree - a.Degree; + Polynomial q = leadingC * new Xd(degreeDifference); + Complex[] rCoefficients = Coefficients.Select(x => x).ToArray(); + rCoefficients[Degree] = 0; + Polynomial r = Dispatch(rCoefficients); //this - q * a; + if (r.Degree < a.Degree || (r.Degree == 0 && a.Degree == 0)) + return new PolynomialDivisionResult(q, r); + PolynomialDivisionResult rDiv = r.LongDivision(a); + return new PolynomialDivisionResult(q + rDiv.Quotient, rDiv.Remainder); + } + + /// + /// + /// + /// + public Polynomial GreatestCommonDivisor(Polynomial o) + { + Polynomial g, l; + if (Degree > o.Degree) + { + g = this; + l = o; + } + else + { + g = o; + l = this; + } + + while (!l.IsZero) + { + PolynomialDivisionResult div = g.LongDivision(l); // = g.LongDivision(l); + g = l; + l = div.Remainder; + } + + return g; + } + + /// + /// + protected void Fix() + { + for (int i = 0; i < Dim; i++) + Coefficients[i] = ComplexFieldStructure.Structure.Fix(Coefficients[i]); + } + + /// + public override Complex Evaluate(Complex x) + => Enumerable + .Range(0, Degree + 1) + .Reverse() + .Select(i => Coefficients[i]) + .Aggregate((x1, x2) => x1 * x + x2); +} \ No newline at end of file diff --git a/src/Analysis/AnalyticFunctions/Polynomials/Implements/Xd.cs b/src/Analysis/AnalyticFunctions/Polynomials/Implements/Xd.cs new file mode 100644 index 0000000..2002c1c --- /dev/null +++ b/src/Analysis/AnalyticFunctions/Polynomials/Implements/Xd.cs @@ -0,0 +1,13 @@ +using System.Numerics; + +namespace Skeleton.Analysis.AnalyticFunctions.Polynomials.Implements; + +/// +/// x^d +/// +public class Xd : GeneralPolynomial +{ + /// + public Xd(int degree) : base(degree, new Complex[degree + 1]) + => Coefficients[^1] = 1; +} diff --git a/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousAddition.cs b/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousAddition.cs new file mode 100644 index 0000000..304c783 --- /dev/null +++ b/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousAddition.cs @@ -0,0 +1,30 @@ +using Skeleton.Analysis.BivariantFunctions.Real.Implements; + +namespace Skeleton.Analysis.BivariantFunctions.Real.Arithmetics; + +/// +/// +public class BivariantContinuousAddition : BivariantContinuous +{ + /// + public BivariantContinuousAddition(BivariantContinuous left, BivariantContinuous right) + { + Left = left; + Right = right; + } + + private BivariantContinuous Left { get; } + private BivariantContinuous Right { get; } + + /// + public override BivariantContinuous Dx => Left.Dx + Right.Dx; + + /// + public override BivariantContinuous Dy => Left.Dy + Right.Dy; + + /// + public override double Evaluate(double x, double y) + { + throw new NotImplementedException(); + } +} \ No newline at end of file diff --git a/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousConstant.cs b/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousConstant.cs new file mode 100644 index 0000000..16e7388 --- /dev/null +++ b/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousConstant.cs @@ -0,0 +1,36 @@ +using Skeleton.Analysis.BivariantFunctions.Real.Implements; + +namespace Skeleton.Analysis.BivariantFunctions.Real.Arithmetics; + +/// +/// +public class BivariantContinuousConstant : BivariantContinuous +{ + /// + /// + public static readonly BivariantContinuousConstant COne = new(1); + + /// + /// + public static readonly BivariantContinuousConstant CZero = new(0); + + /// + public BivariantContinuousConstant(double c) + { + Constant = c; + } + + private double Constant { get; } + + /// + public override BivariantContinuous Dx => CZero; + + /// + public override BivariantContinuous Dy => CZero; + + /// + public override double Evaluate(double x, double y) + { + return Constant; + } +} \ No newline at end of file diff --git a/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousIdentityX.cs b/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousIdentityX.cs new file mode 100644 index 0000000..5017c08 --- /dev/null +++ b/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousIdentityX.cs @@ -0,0 +1,17 @@ +using Skeleton.Analysis.BivariantFunctions.Real.Implements; + +namespace Skeleton.Analysis.BivariantFunctions.Real.Arithmetics; + +/// +/// +public class BivariantContinuousIdentityX : BivariantContinuous +{ + /// + public override BivariantContinuous Dx => BivariantContinuousConstant.COne; + + /// + public override BivariantContinuous Dy => BivariantContinuousConstant.CZero; + + /// + public override double Evaluate(double x, double y) => x; +} \ No newline at end of file diff --git a/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousIdentityY.cs b/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousIdentityY.cs new file mode 100644 index 0000000..fbdee08 --- /dev/null +++ b/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousIdentityY.cs @@ -0,0 +1,17 @@ +using Skeleton.Analysis.BivariantFunctions.Real.Implements; + +namespace Skeleton.Analysis.BivariantFunctions.Real.Arithmetics; + +/// +/// +public class BivariantContinuousIdentityY : BivariantContinuous +{ + /// + public override BivariantContinuous Dx => BivariantContinuousConstant.CZero; + + /// + public override BivariantContinuous Dy => BivariantContinuousConstant.COne; + + /// + public override double Evaluate(double x, double y) => y; +} \ No newline at end of file diff --git a/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousNegation.cs b/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousNegation.cs new file mode 100644 index 0000000..39678eb --- /dev/null +++ b/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousNegation.cs @@ -0,0 +1,22 @@ +using Skeleton.Analysis.BivariantFunctions.Real.Implements; + +namespace Skeleton.Analysis.BivariantFunctions.Real.Arithmetics; + +/// +/// +public class BivariantContinuousNegation : BivariantContinuous +{ + /// + public BivariantContinuousNegation(BivariantContinuous e) => Element = e; + + private BivariantContinuous Element { get; } + + /// + public override BivariantContinuous Dx => -Element.Dx; + + /// + public override BivariantContinuous Dy => -Element.Dy; + + /// + public override double Evaluate(double x, double y) => -Element.Evaluate(x, y); +} \ No newline at end of file diff --git a/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousProduction.cs b/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousProduction.cs new file mode 100644 index 0000000..7df4122 --- /dev/null +++ b/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousProduction.cs @@ -0,0 +1,30 @@ +using Skeleton.Analysis.BivariantFunctions.Real.Implements; + +namespace Skeleton.Analysis.BivariantFunctions.Real.Arithmetics; + +/// +/// +public class BivariantContinuousProduction : BivariantContinuous +{ + /// + public BivariantContinuousProduction(BivariantContinuous l, BivariantContinuous r) + { + Left = l; + Right = r; + } + + private BivariantContinuous Left { get; } + private BivariantContinuous Right { get; } + + /// + public override BivariantContinuous Dx => Left.Dx * Right + Left * Right.Dx; + + /// + public override BivariantContinuous Dy => Left.Dy * Right + Left * Right.Dy; + + /// + public override double Evaluate(double x, double y) + { + return Left.Evaluate(x, y) * Right.Evaluate(x, y); + } +} \ No newline at end of file diff --git a/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousX2.cs b/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousX2.cs new file mode 100644 index 0000000..ed232c5 --- /dev/null +++ b/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousX2.cs @@ -0,0 +1,15 @@ +using Skeleton.Analysis.BivariantFunctions.Real.Implements; + +namespace Skeleton.Analysis.BivariantFunctions.Real.Arithmetics; + +/// +/// +public class BivariantContinuousX2 : BivariantContinuous +{ + /// + public override BivariantContinuous Dx => 2 * X; + /// + public override BivariantContinuous Dy => Zero; + /// + public override double Evaluate(double x, double y) => x * x; +} diff --git a/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousX3.cs b/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousX3.cs new file mode 100644 index 0000000..e908e74 --- /dev/null +++ b/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousX3.cs @@ -0,0 +1,20 @@ +using Skeleton.Analysis.BivariantFunctions.Real.Implements; + +namespace Skeleton.Analysis.BivariantFunctions.Real.Arithmetics; + +/// +/// +public class BivariantContinuousX3 : BivariantContinuous +{ + /// + public override BivariantContinuous Dx => 3 * X2; + + /// + public override BivariantContinuous Dy => Zero; + + /// + public override double Evaluate(double x, double y) + { + return x * x * x; + } +} \ No newline at end of file diff --git a/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousX4.cs b/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousX4.cs new file mode 100644 index 0000000..d1a2bae --- /dev/null +++ b/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousX4.cs @@ -0,0 +1,20 @@ +using Skeleton.Analysis.BivariantFunctions.Real.Implements; + +namespace Skeleton.Analysis.BivariantFunctions.Real.Arithmetics; + +/// +/// +public class BivariantContinuousX4 : BivariantContinuous +{ + /// + public override BivariantContinuous Dx => 4 * X3; + + /// + public override BivariantContinuous Dy => Zero; + + /// + public override double Evaluate(double x, double y) + { + return x * x * x * x; + } +} \ No newline at end of file diff --git a/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousY2.cs b/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousY2.cs new file mode 100644 index 0000000..0880ade --- /dev/null +++ b/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousY2.cs @@ -0,0 +1,20 @@ +using Skeleton.Analysis.BivariantFunctions.Real.Implements; + +namespace Skeleton.Analysis.BivariantFunctions.Real.Arithmetics; + +/// +/// +public class BivariantContinuousY2 : BivariantContinuous +{ + /// + public override BivariantContinuous Dx => Zero; + + /// + public override BivariantContinuous Dy => 2 * Y; + + /// + public override double Evaluate(double x, double y) + { + return y * y; + } +} \ No newline at end of file diff --git a/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousY3.cs b/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousY3.cs new file mode 100644 index 0000000..78408eb --- /dev/null +++ b/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousY3.cs @@ -0,0 +1,20 @@ +using Skeleton.Analysis.BivariantFunctions.Real.Implements; + +namespace Skeleton.Analysis.BivariantFunctions.Real.Arithmetics; + +/// +/// +public class BivariantContinuousY3 : BivariantContinuous +{ + /// + public override BivariantContinuous Dx => Zero; + + /// + public override BivariantContinuous Dy => 3 * Y2; + + /// + public override double Evaluate(double x, double y) + { + return y * y * y; + } +} \ No newline at end of file diff --git a/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousY4.cs b/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousY4.cs new file mode 100644 index 0000000..c5ded45 --- /dev/null +++ b/src/Analysis/BivariantFunctions/Real/Arithmetics/BivariantContinuousY4.cs @@ -0,0 +1,20 @@ +using Skeleton.Analysis.BivariantFunctions.Real.Implements; + +namespace Skeleton.Analysis.BivariantFunctions.Real.Arithmetics; + +/// +/// +public class BivariantContinuousY4 : BivariantContinuous +{ + /// + public override BivariantContinuous Dx => Zero; + + /// + public override BivariantContinuous Dy => 4 * Y3; + + /// + public override double Evaluate(double x, double y) + { + return y * y * y * y; + } +} \ No newline at end of file diff --git a/src/Analysis/BivariantFunctions/Real/IBivariantContinuous.cs b/src/Analysis/BivariantFunctions/Real/IBivariantContinuous.cs new file mode 100644 index 0000000..5bcdd26 --- /dev/null +++ b/src/Analysis/BivariantFunctions/Real/IBivariantContinuous.cs @@ -0,0 +1,23 @@ +using Skeleton.Analysis.BivariantFunctions.Real.Implements; + +namespace Skeleton.Analysis.BivariantFunctions.Real; + +/// +/// +public interface IBivariantContinuous +{ + /// + /// + BivariantContinuous Dx { get; } + + /// + /// + BivariantContinuous Dy { get; } + + /// + /// + /// + /// + /// + double Evaluate(double x, double y); +} \ No newline at end of file diff --git a/src/Analysis/BivariantFunctions/Real/Implements/BivariantContinuous.cs b/src/Analysis/BivariantFunctions/Real/Implements/BivariantContinuous.cs new file mode 100644 index 0000000..3f95639 --- /dev/null +++ b/src/Analysis/BivariantFunctions/Real/Implements/BivariantContinuous.cs @@ -0,0 +1,156 @@ +using Skeleton.Analysis.BivariantFunctions.Real.Arithmetics; + +namespace Skeleton.Analysis.BivariantFunctions.Real.Implements; + +/// +/// +public abstract class BivariantContinuous : IBivariantContinuous +{ + /// + /// + public static readonly BivariantContinuous X = new BivariantContinuousIdentityX(); + + /// + /// + public static readonly BivariantContinuous Y = new BivariantContinuousIdentityY(); + + /// + /// + public static readonly BivariantContinuous Zero = BivariantContinuousConstant.CZero; + + /// + /// + public static readonly BivariantContinuous One = BivariantContinuousConstant.COne; + + /// + /// + public static readonly BivariantContinuous X2 = new BivariantContinuousX2(); + + /// + /// + public static readonly BivariantContinuous Y2 = new BivariantContinuousY2(); + + /// + /// + public static readonly BivariantContinuous X3 = new BivariantContinuousX3(); + + /// + /// + public static readonly BivariantContinuous Y3 = new BivariantContinuousY3(); + + /// + /// + public static readonly BivariantContinuous X4 = new BivariantContinuousX4(); + + /// + /// + public static readonly BivariantContinuous Y4 = new BivariantContinuousY4(); + + /// + public abstract double Evaluate(double x, double y); + + /// + public abstract BivariantContinuous Dx { get; } + + /// + public abstract BivariantContinuous Dy { get; } + + /// + /// + /// + /// + /// + public static BivariantContinuous operator +(BivariantContinuous a, BivariantContinuous b) + { + return new BivariantContinuousAddition(a, b); + } + + /// + /// + /// + /// + public static BivariantContinuous operator -(BivariantContinuous a) + { + return new BivariantContinuousNegation(a); + } + + /// + /// + /// + /// + /// + public static BivariantContinuous operator -(BivariantContinuous a, BivariantContinuous b) + { + return -b + a; + } + + /// + /// + /// + /// + /// + public static BivariantContinuous operator *(BivariantContinuous a, BivariantContinuous b) + { + return new BivariantContinuousProduction(a, b); + } + + /// + /// + /// + /// + /// + public static BivariantContinuous operator +(BivariantContinuous a, double b) + { + return a + new BivariantContinuousConstant(b); + } + + /// + /// + /// + /// + /// + public static BivariantContinuous operator +(double b, BivariantContinuous a) + { + return a + b; + } + + /// + /// + /// + /// + /// + public static BivariantContinuous operator -(BivariantContinuous a, double b) + { + return a + new BivariantContinuousConstant(-b); + } + + /// + /// + /// + /// + /// + public static BivariantContinuous operator -(double a, BivariantContinuous b) + { + return new BivariantContinuousConstant(a) - b; + } + + /// + /// + /// + /// + /// + public static BivariantContinuous operator *(BivariantContinuous a, double b) + { + return a * new BivariantContinuousConstant(b); + } + + /// + /// + /// + /// + /// + public static BivariantContinuous operator *(double b, BivariantContinuous a) + { + return a * b; + } +} \ No newline at end of file diff --git a/src/Analysis/BivariantFunctions/Real/PolynomialHelpers/PolyCHelper.cs b/src/Analysis/BivariantFunctions/Real/PolynomialHelpers/PolyCHelper.cs new file mode 100644 index 0000000..6871130 --- /dev/null +++ b/src/Analysis/BivariantFunctions/Real/PolynomialHelpers/PolyCHelper.cs @@ -0,0 +1,435 @@ +using System.Numerics; +using Skeleton.Analysis.BivariantFunctions.Real.Arithmetics; +using Skeleton.Analysis.BivariantFunctions.Real.Implements; + +namespace Skeleton.Analysis.BivariantFunctions.Real.PolynomialHelpers; + +/// +/// +public static class PolyCHelper +{ + /// + public class PolyC4Re : BivariantContinuous + { + /// + public PolyC4Re(Complex c0, Complex c1, Complex c2, Complex c3, Complex c4) + { + A0 = c0.Real; + B0 = c0.Imaginary; + A1 = c1.Real; + B1 = c1.Imaginary; + A2 = c2.Real; + B2 = c2.Imaginary; + A3 = c3.Real; + B3 = c3.Imaginary; + A4 = c4.Real; + B4 = c4.Imaginary; + } + + private double A0 { get; } + private double B0 { get; set; } + private double A1 { get; } + private double B1 { get; } + private double A2 { get; } + private double B2 { get; } + private double A3 { get; } + private double B3 { get; } + private double A4 { get; } + private double B4 { get; } + + + /// + public override BivariantContinuous Dx => + new PolyC3Re(A1, 0, 2 * A2, 2 * B2, 3 * A3, 3 * B3, 4 * A4, 4 * B4); + + /// + public override BivariantContinuous Dy => + new PolyC3Re(-B1, 0, -2 * B2, 2 * A2, -3 * B3, 3 * A3, 4 * A4, -4 * B4); + + /// + public override double Evaluate(double x, double y) + { + double x2 = x * x; + double y2 = y * y; + double x3 = x2 * x; + double y3 = y2 * y; + double x4 = x2 * x2; + double y4 = y2 * y2; + return + +A0 + + A1 * x - B1 * y + + A2 * (x2 - y2) - B2 * 2 * x * y + + A3 * (x3 - 3 * x * y2) + B3 * (y3 - 3 * x2 * y) + + A4 * (x4 - 6 * x2 * y2 + y4) - B4 * (4 * (x3 * y + x * y3)); + } + } + + /// + public class PolyC4Im : BivariantContinuous + { + /// + public PolyC4Im(Complex c0, Complex c1, Complex c2, Complex c3, Complex c4) + { + A0 = c0.Real; + B0 = c0.Imaginary; + A1 = c1.Real; + B1 = c1.Imaginary; + A2 = c2.Real; + B2 = c2.Imaginary; + A3 = c3.Real; + B3 = c3.Imaginary; + A4 = c4.Real; + B4 = c4.Imaginary; + } + + private double A0 { get; set; } + private double B0 { get; } + private double A1 { get; } + private double B1 { get; } + private double A2 { get; } + private double B2 { get; } + private double A3 { get; } + private double B3 { get; } + private double A4 { get; } + private double B4 { get; } + + /// + public override BivariantContinuous Dx => + new PolyC3Im(0, B1, 2 * A2, 2 * B2, 3 * A3, 3 * B3, 4 * A4, 4 * B4); + + /// + public override BivariantContinuous Dy => + new PolyC3Im(0, A1, -2 * B2, 2 * A2, -3 * B3, 3 * A3, -4 * B4, 4 * A4); + + /// + public override double Evaluate(double x, double y) + { + double x2 = x * x; + double y2 = y * y; + double x3 = x2 * x; + double y3 = y2 * y; + double x4 = x2 * x2; + double y4 = y2 * y2; + return + +B0 + + A1 * y + B1 * x + + A2 * 2 * x * y + B2 * (x2 - y2) + + A3 * (3 * x2 * y - y3) + B3 * (x3 - 3 * x * y2) + + A4 * (4 * (x3 * y - x * y3)) + B4 * (x4 - 6 * x2 * y2 + y4); + } + } + + + /// + public class PolyC3Re : BivariantContinuous + { + /// + public PolyC3Re(Complex c0, Complex c1, Complex c2, Complex c3) + { + A0 = c0.Real; + B0 = c0.Imaginary; + A1 = c1.Real; + B1 = c1.Imaginary; + A2 = c2.Real; + B2 = c2.Imaginary; + A3 = c3.Real; + B3 = c3.Imaginary; + } + + /// + public PolyC3Re(double a0, double b0, double a1, double b1, double a2, double b2, double a3, double b3) + { + A0 = a0; + B0 = b0; + A1 = a1; + B1 = b1; + A2 = a2; + B2 = b2; + A3 = a3; + B3 = b3; + } + + private double A0 { get; } + private double B0 { get; set; } + private double A1 { get; } + private double B1 { get; } + private double A2 { get; } + private double B2 { get; } + private double A3 { get; } + private double B3 { get; } + + /// + public override BivariantContinuous Dx => new PolyC2Re(A1, 0, 2 * A2, 2 * B2, 3 * A3, 3 * B3); + + /// + public override BivariantContinuous Dy => new PolyC2Re(-B1, 0, -2 * B2, 2 * A2, -3 * B3, 3 * A3); + + /// + public override double Evaluate(double x, double y) + { + double x2 = x * x; + double y2 = y * y; + double x3 = x2 * x; + double y3 = y2 * y; + //ddx => + // A1 +2A2 x - 2B2 y + 3A3 x2 - 6B3 xy -3A3 y2 + + //ddy => + // -B1 -2B2 x -2A2 y - 3B3 X2 - 6A3 xy + 3B3 y2 + return + +A0 + + A1 * x - B1 * y + + A2 * (x2 - y2) - B2 * 2 * x * y + + A3 * (x3 - 3 * x * y2) + B3 * (y3 - 3 * x2 * y); + } + } + + /// + public class PolyC3Im : BivariantContinuous + { + /// + public PolyC3Im(Complex c0, Complex c1, Complex c2, Complex c3) + { + A0 = c0.Real; + B0 = c0.Imaginary; + A1 = c1.Real; + B1 = c1.Imaginary; + A2 = c2.Real; + B2 = c2.Imaginary; + A3 = c3.Real; + B3 = c3.Imaginary; + } + + /// + public PolyC3Im(double a0, double b0, double a1, double b1, double a2, double b2, double a3, double b3) + { + A0 = a0; + B0 = b0; + A1 = a1; + B1 = b1; + A2 = a2; + B2 = b2; + A3 = a3; + B3 = b3; + } + + private double A0 { get; set; } + private double B0 { get; } + private double A1 { get; } + private double B1 { get; } + private double A2 { get; } + private double B2 { get; } + private double A3 { get; } + private double B3 { get; } + + + /// + public override BivariantContinuous Dx => new PolyC2Im(0, B1, 2 * A2, 2 * B2, 3 * A3, 3 * B3); + + /// + public override BivariantContinuous Dy => new PolyC2Im(0, A1, -2 * B2, 2 * A2, -3 * B3, 3 * A3); + + /// + public override double Evaluate(double x, double y) + { + double x2 = x * x; + double y2 = y * y; + double x3 = x2 * x; + double y3 = y2 * y; + //ddx => + // B1 + 2B2 x + 2A2 y + 6A3 xy + 3B3 (x2 - y2) + + //ddy => + // A1 + 2A2 x - 2B2 y - 6B3 xy + 3A3 (x2 - y2) + return + +B0 + + A1 * y + B1 * x + + A2 * 2 * x * y + B2 * (x2 - y2) + + -A3 * (y3 - 3 * x2 * y) + B3 * (x3 - 3 * x * y2); + } + } + + /// + public class PolyC2Re : BivariantContinuous + { + /// + public PolyC2Re(Complex c0, Complex c1, Complex c2) + { + A0 = c0.Real; + B0 = c0.Imaginary; + A1 = c1.Real; + B1 = c1.Imaginary; + A2 = c2.Real; + B2 = c2.Imaginary; + } + + /// + public PolyC2Re(double a0, double b0, double a1, double b1, double a2, double b2) + { + A0 = a0; + B0 = b0; + A1 = a1; + B1 = b1; + A2 = a2; + B2 = b2; + } + + private double A0 { get; } + private double B0 { get; set; } + private double A1 { get; } + private double B1 { get; } + private double A2 { get; } + private double B2 { get; } + + /// + public override BivariantContinuous Dx => new PolyC1Re(A1, 0, 2 * A2, 2 * B2); + + /// + public override BivariantContinuous Dy => new PolyC1Re(-B1, 0, -2 * B2, 2 * A2); + + /// + public override double Evaluate(double x, double y) + { + double x2 = x * x; + double y2 = y * y; + return //-B1 -2B2x - 2A2y + +A0 + + A1 * x - B1 * y + + A2 * (x2 - y2) - B2 * 2 * x * y; + } + } + + /// + public class PolyC2Im : BivariantContinuous + { + /// + public PolyC2Im(Complex c0, Complex c1, Complex c2) + { + A0 = c0.Real; + B0 = c0.Imaginary; + A1 = c1.Real; + B1 = c1.Imaginary; + A2 = c2.Real; + B2 = c2.Imaginary; + } + + /// + public PolyC2Im(double a0, double b0, double a1, double b1, double a2, double b2) + { + A0 = a0; + B0 = b0; + A1 = a1; + B1 = b1; + A2 = a2; + B2 = b2; + } + + private double A0 { get; set; } + private double B0 { get; } + private double A1 { get; } + private double B1 { get; } + private double A2 { get; } + private double B2 { get; } + + + /// + public override BivariantContinuous Dx => new PolyC1Im(0, B1, 2 * B2, 2 * A2); + + /// + public override BivariantContinuous Dy => new PolyC1Im(0, A1, -2 * B2, 2 * A2); + + /// + public override double Evaluate(double x, double y) + { + double x2 = x * x; + double y2 = y * y; + return // B1 + 2B2x + 2A2y + +B0 + + A1 * y + B1 * x + + +A2 * 2 * x * y + B2 * (x2 - y2); + } + } + + /// + public class PolyC1Re : BivariantContinuous + { + /// + public PolyC1Re(Complex c0, Complex c1) + { + A0 = c0.Real; + B0 = c0.Imaginary; + A1 = c1.Real; + B1 = c1.Imaginary; + } + + /// + public PolyC1Re(double a0, double b0, double a1, double b1) + { + A0 = a0; + B0 = b0; + A1 = a1; + B1 = b1; + } + + private double A0 { get; } + private double B0 { get; set; } + private double A1 { get; } + private double B1 { get; } + + + /// + public override BivariantContinuous Dx => new BivariantContinuousConstant(A1); + + /// + public override BivariantContinuous Dy => new BivariantContinuousConstant(-B1); + + /// + public override double Evaluate(double x, double y) + { + return + +A0 + + A1 * x - B1 * y; + } + } + + /// + public class PolyC1Im : BivariantContinuous + { + /// + public PolyC1Im(Complex c0, Complex c1) + { + A0 = c0.Real; + B0 = c0.Imaginary; + A1 = c1.Real; + B1 = c1.Imaginary; + } + + /// + public PolyC1Im(double a0, double b0, double a1, double b1) + { + A0 = a0; + B0 = b0; + A1 = a1; + B1 = b1; + } + + private double A0 { get; set; } + private double B0 { get; } + private double A1 { get; } + private double B1 { get; } + + /// + public override BivariantContinuous Dx => new BivariantContinuousConstant(B1); + + /// + public override BivariantContinuous Dy => new BivariantContinuousConstant(A1); + + /// + public override double Evaluate(double x, double y) + { + return + +B0 + + A1 * y + B1 * x; + } + } +} \ No newline at end of file diff --git a/src/Analysis/ContinuousFunctions/Arithmetics/ContinuousAddition.cs b/src/Analysis/ContinuousFunctions/Arithmetics/ContinuousAddition.cs new file mode 100644 index 0000000..a0462d5 --- /dev/null +++ b/src/Analysis/ContinuousFunctions/Arithmetics/ContinuousAddition.cs @@ -0,0 +1,27 @@ +using Skeleton.Analysis.ContinuousFunctions.Implements; + +namespace Skeleton.Analysis.ContinuousFunctions.Arithmetics; + +/// +/// +public class ContinuousAddition : Continuous +{ + /// + public ContinuousAddition(Continuous left, Continuous right) + { + Left = left; + Right = right; + } + + private Continuous Left { get; } + private Continuous Right { get; } + + /// + public override Continuous Derivative => Left.Derivative + Right.Derivative; + + /// + public override double Evaluate(double x) + { + return Left.Evaluate(x) + Right.Evaluate(x); + } +} \ No newline at end of file diff --git a/src/Analysis/ContinuousFunctions/Arithmetics/ContinuousApply.cs b/src/Analysis/ContinuousFunctions/Arithmetics/ContinuousApply.cs new file mode 100644 index 0000000..6ae0193 --- /dev/null +++ b/src/Analysis/ContinuousFunctions/Arithmetics/ContinuousApply.cs @@ -0,0 +1,27 @@ +using Skeleton.Analysis.ContinuousFunctions.Implements; + +namespace Skeleton.Analysis.ContinuousFunctions.Arithmetics; + +/// +/// +public class ContinuousApply : Continuous +{ + /// + public ContinuousApply(Continuous o, Continuous i) + { + Outer = o; + Inner = i; + } + + private Continuous Outer { get; } + private Continuous Inner { get; } + + /// + public override Continuous Derivative => Inner.Derivative * Outer.Derivative[Inner]; + + /// + public override double Evaluate(double x) + { + return Outer.Evaluate(Inner.Evaluate(x)); + } +} \ No newline at end of file diff --git a/src/Analysis/ContinuousFunctions/Arithmetics/ContinuousConstant.cs b/src/Analysis/ContinuousFunctions/Arithmetics/ContinuousConstant.cs new file mode 100644 index 0000000..ae8fdf9 --- /dev/null +++ b/src/Analysis/ContinuousFunctions/Arithmetics/ContinuousConstant.cs @@ -0,0 +1,33 @@ +using Skeleton.Analysis.ContinuousFunctions.Implements; + +namespace Skeleton.Analysis.ContinuousFunctions.Arithmetics; + +/// +/// +public class ContinuousConstant : Continuous +{ + /// + /// + public static readonly ContinuousConstant Zero = new(0); + + /// + /// + public static readonly ContinuousConstant One = new(1); + + /// + public ContinuousConstant(double c) + { + Constant = c; + } + + private double Constant { get; } + + /// + public override Continuous Derivative => Zero; + + /// + public override double Evaluate(double x) + { + return Constant; + } +} \ No newline at end of file diff --git a/src/Analysis/ContinuousFunctions/Arithmetics/ContinuousNegation.cs b/src/Analysis/ContinuousFunctions/Arithmetics/ContinuousNegation.cs new file mode 100644 index 0000000..a3cc4b5 --- /dev/null +++ b/src/Analysis/ContinuousFunctions/Arithmetics/ContinuousNegation.cs @@ -0,0 +1,25 @@ +using Skeleton.Analysis.ContinuousFunctions.Implements; + +namespace Skeleton.Analysis.ContinuousFunctions.Arithmetics; + +/// +/// +public class ContinuousNegation : Continuous +{ + /// + public ContinuousNegation(Continuous e) + { + Element = e; + } + + private Continuous Element { get; } + + /// + public override Continuous Derivative => -Element.Derivative; + + /// + public override double Evaluate(double x) + { + return -Element.Evaluate(x); + } +} \ No newline at end of file diff --git a/src/Analysis/ContinuousFunctions/Arithmetics/ContinuousProduction.cs b/src/Analysis/ContinuousFunctions/Arithmetics/ContinuousProduction.cs new file mode 100644 index 0000000..b21afed --- /dev/null +++ b/src/Analysis/ContinuousFunctions/Arithmetics/ContinuousProduction.cs @@ -0,0 +1,27 @@ +using Skeleton.Analysis.ContinuousFunctions.Implements; + +namespace Skeleton.Analysis.ContinuousFunctions.Arithmetics; + +/// +/// +public class ContinuousProduction : Continuous +{ + /// + public ContinuousProduction(Continuous left, Continuous right) + { + Left = left; + Right = right; + } + + private Continuous Left { get; } + private Continuous Right { get; } + + /// + public override Continuous Derivative => Left * Right.Derivative + Left.Derivative * Right; + + /// + public override double Evaluate(double x) + { + return Left.Evaluate(x) * Right.Evaluate(x); + } +} \ No newline at end of file diff --git a/src/Analysis/ContinuousFunctions/IContinuous.cs b/src/Analysis/ContinuousFunctions/IContinuous.cs new file mode 100644 index 0000000..0ce041f --- /dev/null +++ b/src/Analysis/ContinuousFunctions/IContinuous.cs @@ -0,0 +1,18 @@ +using Skeleton.Analysis.ContinuousFunctions.Implements; + +namespace Skeleton.Analysis.ContinuousFunctions; + +/// +/// +public interface IContinuous +{ + /// + /// + Continuous Derivative { get; } + + /// + /// + /// + /// + double Evaluate(double x); +} \ No newline at end of file diff --git a/src/Analysis/ContinuousFunctions/Implements/Continuous.cs b/src/Analysis/ContinuousFunctions/Implements/Continuous.cs new file mode 100644 index 0000000..1e8d970 --- /dev/null +++ b/src/Analysis/ContinuousFunctions/Implements/Continuous.cs @@ -0,0 +1,58 @@ +using Skeleton.Analysis.ContinuousFunctions.Arithmetics; + +namespace Skeleton.Analysis.ContinuousFunctions.Implements; + +/// +/// +public abstract class Continuous : IContinuous +{ + /// + /// + /// + public Continuous this[Continuous c] => new ContinuousApply(this, c); + + /// + public abstract Continuous Derivative { get; } + + /// + public abstract double Evaluate(double x); + + /// + /// + /// + /// + /// + public static Continuous operator +(Continuous a, Continuous b) + { + return new ContinuousAddition(a, b); + } + + /// + /// + /// + /// + public static Continuous operator -(Continuous a) + { + return new ContinuousNegation(a); + } + + /// + /// + /// + /// + /// + public static Continuous operator -(Continuous a, Continuous b) + { + return -b + a; + } + + /// + /// + /// + /// + /// + public static Continuous operator *(Continuous a, Continuous b) + { + return new ContinuousProduction(a, b); + } +} \ No newline at end of file diff --git a/src/Constants/AlgebraConstant.cs b/src/Constants/AlgebraConstant.cs new file mode 100644 index 0000000..24e8173 --- /dev/null +++ b/src/Constants/AlgebraConstant.cs @@ -0,0 +1,100 @@ +using System.Numerics; +using Skeleton.Algebra.DimensionProviders; + +namespace Skeleton.Constants; +using DiagR4 = Algebra.CategoryOf.OnField.FDiagonalMatrix; +using DiagR3 = Algebra.CategoryOf.OnField.FDiagonalMatrix; + +/// +/// constants +/// +public static class AlgebraConstant +{ + /// + /// sqrt of pi + /// + public const double SqrtPi = 1.772453850905516027298167483341145182797549456122387128213807789852911284183d; + + /// + /// erf(1) + /// + public const double Erf1 = 0.8427007877600066754197883d; + + /// + /// erfc(1) + /// + public const double Erfc1 = 1 - Erf1; + + /// + /// erf(2) + /// + public const double Erf2 = 0.9953222650189527341620691d; + + /// + /// erfc(2) + /// + public const double Erfc2 = 1 - Erf2; + + /// + /// pi * pi + /// + public const double PI2 = Math.PI * Math.PI; + + /// + /// imaginary one + /// + public static readonly Complex I = Complex.ImaginaryOne; + + /// + /// 0, 1/3, 2/3, 1 as diag + /// + public static readonly DiagR4 UniformSpectrum4 = new(0d, 1d / 3d, 2d / 3d, 1d); + + /// + /// 0, 0.5, 1 as diag + /// + public static readonly DiagR3 UniformSpectrum3 = new(0d, 1d / 2d, 1d); + + /// + /// no use + /// + /// + /// + public static int Fib(int x) + { + return x switch + { + 1 => 1, + 2 => 2, + 3 => 3, + 4 => 5, + 5 => 8, + 6 => 13, + 7 => 21, + 8 => 34, + _ => -1 + }; + } + + /// + /// no use + /// + /// + /// + public static int Ludic(int x) + { + return x switch + { + 0 => 0, + 1 => 1, + 2 => 2, + 3 => 4, + 4 => 16, + 5 => 26, + 6 => 42, + 7 => 57, + 8 => 79, + _ => -1 + }; + } +} \ No newline at end of file diff --git a/src/Constants/CC.cs b/src/Constants/CC.cs new file mode 100644 index 0000000..cb1ab98 --- /dev/null +++ b/src/Constants/CC.cs @@ -0,0 +1,76 @@ +using System.Numerics; + +namespace Skeleton.Constants; + +/// +/// +public static class CC +{ + /// + /// + public static readonly Complex I = Complex.ImaginaryOne; + + /// + /// + public static readonly double Sqrt2 = Math.Sqrt(2d); + + /// + /// + public static readonly double Sqrt3 = Math.Sqrt(3d); + + /// + /// + public static readonly double Sqrt5 = Math.Sqrt(5d); + + /// + /// + public static readonly double Sqrt6 = Math.Sqrt(6d); + + /// + /// + public static readonly double Sqrt7 = Math.Sqrt(7d); + + /// + /// + public static readonly double Sqrt10 = Math.Sqrt(10d); + + /// + /// + public static readonly double Sqrt11 = Math.Sqrt(11d); + + /// + /// + public static readonly double Sqrt13 = Math.Sqrt(13d); + + /// + /// + public static readonly double SqrtInv2 = 1d / Sqrt2; + + /// + /// + public static readonly double SqrtInv3 = 1d / Sqrt3; + + /// + /// + public static readonly double SqrtInv5 = 1d / Sqrt5; + + /// + /// + public static readonly double SqrtInv6 = 1d / Sqrt6; + + /// + /// + public static readonly double SqrtInv7 = 1d / Sqrt7; + + /// + /// + public static readonly double SqrtInv10 = 1d / Sqrt10; + + /// + /// + public static readonly double SqrtInv11 = 1d / Sqrt11; + + /// + /// + public static readonly double SqrtInv13 = 1d / Sqrt13; +} \ No newline at end of file diff --git a/src/DataStructure/CacheItem.cs b/src/DataStructure/CacheItem.cs new file mode 100644 index 0000000..45ae10f --- /dev/null +++ b/src/DataStructure/CacheItem.cs @@ -0,0 +1,82 @@ +namespace Skeleton.DataStructure; + +/// +/// store value from calculation to avoid redundant computations +/// +public abstract class CacheItem +{ + /// + /// if the item needs update + /// + public bool Expired { get; protected set; } = true; + internal HashSet References { get; set; } = new(); + internal HashSet Dependencies { get; set; } = new(); + + /// + /// tag the item and all dependencies as out dated + /// + public void Expire() + { + Expired = true; + foreach (CacheItem c in References) + if (!c.Expired) + c.Expire(); + } +} + +/// +public class CacheItem : CacheItem +{ + /// + /// construct by provide a method to update the value + /// + /// + public CacheItem(Func rec) => ProxyCalculator = rec; + private TObject? Value { get; set; } + /// + /// return cached/updated value when reference is not used to calculate another cache item + /// + public TObject Get + { + get + { + if (Expired) + { + Value = Calculator(); + Expired = false; + } + return Value!; + } + } + private Func Calculator => () => ProxyCalculator(this); + /// + /// provide trace information while calculating + /// + public Func ProxyCalculator { get; internal set; } + /// + /// trace the dependency and return cached value/updated value + /// + /// + /// + public TObject GetFrom(CacheItem d) + { + References.Add(d); + d.Dependencies.Add(this); + return Get; + } + /// + /// update the formula of item value + /// + /// + public void UpdateCalculation(Func val) + { + Expire(); + foreach (CacheItem item in Dependencies) + item.References.Remove(this); + Dependencies = new HashSet(); + ProxyCalculator = val; + } + + +} + diff --git a/src/DataStructure/IsomorphicMap.cs b/src/DataStructure/IsomorphicMap.cs new file mode 100644 index 0000000..dacdf45 --- /dev/null +++ b/src/DataStructure/IsomorphicMap.cs @@ -0,0 +1,60 @@ +namespace Skeleton.DataStructure; + +/// +/// +/// +public class IsomorphicMap + where TKey : notnull + where TValue : notnull +{ + private Dictionary F { get; set; } = new(); + private Dictionary G { get; set; } = new(); + + /// + /// + /// + /// + public TValue this[TKey key] + { + get => F[key]; + set + { + F[key] = value; + G[value] = key; + } + } + + /// + /// + /// + /// + public TKey this[TValue key] + { + get => G[key]; + set + { + F[value] = key; + G[key] = value; + } + } + + public bool ContainsKey(TKey key) => F.ContainsKey(key); + public bool ContainsKey(TValue key) => G.ContainsKey(key); + public IEnumerable Values => G.Keys; + public IEnumerable Keys => F.Keys; + + public void Remove(TKey key) + { + TValue v = F[key]; + F.Remove(key); + G.Remove(v); + } + + public void Remove(TValue val) + { + TKey key = G[val]; + F.Remove(key); + G.Remove(val); + } + +} \ No newline at end of file diff --git a/src/DataStructure/Link/Link.cs b/src/DataStructure/Link/Link.cs new file mode 100644 index 0000000..9896e55 --- /dev/null +++ b/src/DataStructure/Link/Link.cs @@ -0,0 +1,236 @@ +using System.Collections; + +namespace Skeleton.DataStructure.Link; + +/// +/// Linked List +/// +/// +public class Link : IEnumerable +{ + /// + /// singleton for default value + /// + public static readonly Link Null = new(); + + /// + /// Constructor + /// + public Link() + { + Head = new LinkHead(); + Tail = new LinkTail(); + Head.Parent = this; + Tail.Parent = this; + Head.Next = Tail; + Tail.Previous = Head; + Count = 0; + } + + /// + /// First item + /// + public LinkNode First => Head.Next; + + /// + /// Last item + /// + public LinkNode Last => Tail.Previous; + + /// + /// number of items + /// + public int Count { get; set; } + + /// + /// head of linked list + /// + private LinkHead Head { get; } + + /// + /// tail of linked list + /// + private LinkTail Tail { get; } + + /// + /// iterator in reverse direction + /// + /// + public IEnumerable Backward + { + get + { + LinkNode iter = Tail.Previous; + int current = 1; + while (!iter.IsHead) + { + if (current > Count) + throw new Exception("TODO"); + yield return iter.Value; + iter = iter.Previous; + current += 1; + } + } + } + + /// + /// + /// + /// + public IEnumerator GetEnumerator() + { + LinkNode iter = Head.Next; + int current = 1; + while (!iter.IsTail) + { + if (current > Count) + throw new Exception("TODO"); + yield return iter.Value; + iter = iter.Next; + current += 1; + } + } + + IEnumerator IEnumerable.GetEnumerator() + { + return GetEnumerator(); + } + + /// + /// add item next to head + /// + /// + public void AddFirst(TObject obj) + { + LinkNode node = new LinkNode(obj); + node.Parent = this; + Head.Next.Previous = node; + node.Previous = Head; + node.Next = Head.Next; + Head.Next = node; + Count += 1; + } + + /// + /// add item before tail + /// + /// + public void AddLast(TObject obj) + { + LinkNode node = new LinkNode(obj); + node.Parent = this; + Tail.Previous.Next = node; + node.Next = Tail; + node.Previous = Tail.Previous; + Tail.Previous = node; + Count += 1; + } + + /// + /// add item after specific node + /// + /// + /// + /// + public void AddAfter(LinkNode node, TObject obj) + { + if (node.Parent != this) + throw new Exception("TODO"); + LinkNode newNode = new LinkNode(obj); + node.Parent = this; + node.Next.Previous = newNode; + newNode.Next = node.Next; + node.Next = newNode; + newNode.Previous = node; + Count += 1; + } + + /// + /// add item before specific node + /// + /// + /// + public void AddBefore(LinkNode node, TObject obj) + { + AddAfter(node.Previous, obj); + } + + /// + /// + /// + /// + /// + public LinkNode Find(TObject value) + { + LinkNode iter = Head.Next; + int current = 1; + while (!iter.IsTail) + { + if (current > Count) + throw new Exception("TODO"); + if (iter.Value!.Equals(value)) + return iter; + iter = iter.Next; + current += 1; + } + + return LinkNode.Null; + } + + /// + /// find the node by selector + /// + /// + /// + /// + public LinkNode Find(Func selector) + { + LinkNode iter = Head.Next; + int current = 1; + while (!iter.IsTail) + { + if (current > Count) + throw new Exception("TODO"); + if (selector(iter.Value)) + return iter; + iter = iter.Next; + current += 1; + } + + return LinkNode.Null; + } + + /// + /// remove specific node + /// + /// + /// + public void Remove(LinkNode node) + { + if (!ReferenceEquals(node.Parent, this)) + throw new Exception(); + if (node.IsEnding) + throw new Exception(); + node.Previous.Next = node.Next; + node.Next.Previous = node.Previous; + node.Previous = LinkNode.Null; + node.Next = LinkNode.Null; + node.Parent = Null; + Count -= 1; + } + + /// + /// remove the first item + /// + /// + /// + public TObject RemoveFirst() + { + LinkNode res = First; + if (res.IsEnding) + throw new Exception("TODO EMPTY"); + TObject resValue = res.Value; + Remove(res); + return resValue; + } +} \ No newline at end of file diff --git a/src/DataStructure/Link/LinkHead.cs b/src/DataStructure/Link/LinkHead.cs new file mode 100644 index 0000000..fa47356 --- /dev/null +++ b/src/DataStructure/Link/LinkHead.cs @@ -0,0 +1,38 @@ +namespace Skeleton.DataStructure.Link; + +/// +/// head of linked list +/// +/// +public class LinkHead : LinkNode +{ + /// + /// no previous element + /// + /// + public override LinkNode Previous + { + get => throw new Exception("TODO"); + set => throw new Exception("TODO"); + } + + /// + /// no value + /// + /// + public new TObject Value + { + get => throw new Exception("TODO"); + set => throw new Exception("TODO"); + } + + /// + /// is head + /// + public override bool IsHead => true; + + /// + /// is ending + /// + public override bool IsEnding => true; +} \ No newline at end of file diff --git a/src/DataStructure/Link/LinkNode.cs b/src/DataStructure/Link/LinkNode.cs new file mode 100644 index 0000000..32c6130 --- /dev/null +++ b/src/DataStructure/Link/LinkNode.cs @@ -0,0 +1,164 @@ +namespace Skeleton.DataStructure.Link; + +/// +/// node of linked list +/// +/// +public class LinkNode +{ + /// + /// singleton for initialize + /// + public static readonly LinkNode Null = new(); + + /// + /// Constructor + /// + protected LinkNode() { } + + /// + /// Constructor + /// + /// + public LinkNode(TObject val) + { + Value = val; + } + + /// + /// by default, node is not the tail + /// + public virtual bool IsTail => false; + + /// + /// by default, node is not the head + /// + public virtual bool IsHead => false; + + /// + /// by default, node is not head nor tail + /// + public virtual bool IsEnding => false; + + /// + /// unexpected + /// + public bool Isolated => !IsEnding && (Next == Null || Previous == Null); + + /// + /// linked list that node belongs to + /// + public Link Parent { get; set; } = Link.Null; + + /// + /// next node + /// + public virtual LinkNode Next { get; set; } = Null; + + /// + /// previous node + /// + public virtual LinkNode Previous { get; set; } = Null; + + + private TObject? PrivateValue { get; set; } + + /// + /// item + /// + public TObject Value + { + get => PrivateValue ?? throw new Exception(); + set => PrivateValue = value; + } + + /// + /// move item forward + /// + /// + public void MoveForward() + { + if (IsEnding) + throw new Exception("TODO"); + if (Next.IsTail) + return; + LinkNode cOrder0 = Previous; + LinkNode cOrder1 = Next; + LinkNode cOrder2 = Next.Next; + cOrder0.Next = cOrder1; + cOrder1.Previous = cOrder0; + cOrder1.Next = this; + Previous = cOrder1; + Next = cOrder2; + cOrder2.Previous = this; + } + + /// + /// move item backward + /// + /// + public void MoveBackward() + { + if (IsEnding) + throw new Exception("TODO"); + if (Previous.IsHead) + return; + + LinkNode cOrder0 = Previous.Previous; + LinkNode cOrder1 = Previous; + LinkNode cOrder2 = Next; + + cOrder0.Next = this; + Previous = cOrder0; + cOrder1.Previous = this; + Next = cOrder1; + cOrder2.Previous = cOrder1; + cOrder1.Next = cOrder2; + } + + /// + /// move node after specific node + /// + /// + /// + public void MoveAfter(LinkNode target) + { + if (target.IsTail) + throw new Exception("TODO"); + if (target.Parent != Parent) + throw new Exception("TODO"); + if (!Isolated) + { + Next.Previous = Previous; + Previous.Next = Next; + } + + Previous = target; + Next = target.Next; + target.Next.Previous = this; + target.Next = this; + } + + /// + /// move item before specific node + /// + /// + /// + public void MoveBefore(LinkNode target) + { + if (target.IsHead) + throw new Exception("TODO"); + if (!ReferenceEquals(target.Parent, Parent)) + throw new Exception("TODO"); + if (!Isolated) + { + Next.Previous = Previous; + Previous.Next = Next; + } + + Previous = target.Previous; + Next = target; + target.Previous.Next = this; + target.Previous = this; + } +} \ No newline at end of file diff --git a/src/DataStructure/Link/LinkTail.cs b/src/DataStructure/Link/LinkTail.cs new file mode 100644 index 0000000..61cd2c5 --- /dev/null +++ b/src/DataStructure/Link/LinkTail.cs @@ -0,0 +1,38 @@ +namespace Skeleton.DataStructure.Link; + +/// +/// tail of linked list +/// +/// +public class LinkTail : LinkNode +{ + /// + /// no next + /// + /// + public override LinkNode Next + { + get => throw new Exception("TODO"); + set => throw new Exception("TODO"); + } + + /// + /// no value + /// + /// + public new TObject Value + { + get => throw new Exception("TODO"); + set => throw new Exception("TODO"); + } + + /// + /// is tail + /// + public override bool IsTail => true; + + /// + /// is ending + /// + public override bool IsEnding => true; +} \ No newline at end of file diff --git a/src/DataStructure/Packs/EigenPair.cs b/src/DataStructure/Packs/EigenPair.cs new file mode 100644 index 0000000..9c4a470 --- /dev/null +++ b/src/DataStructure/Packs/EigenPair.cs @@ -0,0 +1,26 @@ +using Skeleton.Algebra.Vectors; + +namespace Skeleton.DataStructure.Packs; + +/// +/// +/// +/// +public class EigenPair : Pack2 + where TVector : Vector, new() +{ + /// + public EigenPair(TScalar eigenValue, TVector[] eigenVectors) : base(eigenValue, eigenVectors) { } + + /// + /// + public TScalar EigenValue + { + get => Item1; + set => Item1 = value; + } + + /// + /// + public TVector[] EigenVectors => Item2; +} \ No newline at end of file diff --git a/src/DataStructure/Packs/LinearSystems/LinearSystem.cs b/src/DataStructure/Packs/LinearSystems/LinearSystem.cs new file mode 100644 index 0000000..c29a47a --- /dev/null +++ b/src/DataStructure/Packs/LinearSystems/LinearSystem.cs @@ -0,0 +1,37 @@ +using Skeleton.Algebra.Matrices; +using Skeleton.Algebra.Vectors; + +namespace Skeleton.DataStructure.Packs.LinearSystems; + +/// +/// system Ax = y +/// +/// +/// +/// +public class LinearSystem : Pack2 + where TVector : Vector, new() + where TMatrix : Matrix, new() + where TScalar : notnull +{ + /// + /// constructor + /// + /// + /// + public LinearSystem(TMatrix m, TVector v) + { + Item1 = m; + Item2 = v; + } + + /// + /// matrix + /// + public TMatrix A => Item1; + + /// + /// y + /// + public TVector Y => Item2; +} \ No newline at end of file diff --git a/src/DataStructure/Packs/MatrixDecompositions/GeneralUHPair.cs b/src/DataStructure/Packs/MatrixDecompositions/GeneralUHPair.cs new file mode 100644 index 0000000..8fc7088 --- /dev/null +++ b/src/DataStructure/Packs/MatrixDecompositions/GeneralUHPair.cs @@ -0,0 +1,22 @@ +namespace Skeleton.DataStructure.Packs.MatrixDecompositions; + +/// +/// +/// +/// +public class GeneralUHPair : Pack2 +{ + /// + /// + /// + /// + public GeneralUHPair(U a, H b) : base(a, b) { } + + /// + /// + public U Unitary => Item1; + + /// + /// + public H Hermitian => Item2; +} \ No newline at end of file diff --git a/src/DataStructure/Packs/MatrixDecompositions/JDPair.cs b/src/DataStructure/Packs/MatrixDecompositions/JDPair.cs new file mode 100644 index 0000000..c59507a --- /dev/null +++ b/src/DataStructure/Packs/MatrixDecompositions/JDPair.cs @@ -0,0 +1,21 @@ +namespace Skeleton.DataStructure.Packs.MatrixDecompositions; + +/// +/// +/// +public class JDPair : Pack2 +{ + /// + /// + /// + /// + public JDPair(TMat a, TMat diagonal) : base(a, diagonal) { } + + /// + /// + public TMat J => Item1; + + /// + /// + public TMat D => Item2; +} \ No newline at end of file diff --git a/src/DataStructure/Packs/MatrixDecompositions/PLUPack.cs b/src/DataStructure/Packs/MatrixDecompositions/PLUPack.cs new file mode 100644 index 0000000..4cc8308 --- /dev/null +++ b/src/DataStructure/Packs/MatrixDecompositions/PLUPack.cs @@ -0,0 +1,37 @@ +namespace Skeleton.DataStructure.Packs.MatrixDecompositions; + +/// +/// plu decomposition result +/// +/// +/// +public class PLUPack : Pack3 +{ + /// + /// + /// + /// + /// + /// + /// + public PLUPack(TMatrix p, TMatrix l, TMatrix u, TScalar det) : base(p, l, u) => Det = det; + + /// + /// permutation + /// + public TMatrix P => Item1; + /// + /// lower triangle + /// + public TMatrix L => Item2; + /// + /// upper triangle + /// + public TMatrix U => Item3; + + /// + /// det of mat + /// + public TScalar Det { get; set; } + +} \ No newline at end of file diff --git a/src/DataStructure/Packs/MatrixDecompositions/QRPair.cs b/src/DataStructure/Packs/MatrixDecompositions/QRPair.cs new file mode 100644 index 0000000..9ed33ce --- /dev/null +++ b/src/DataStructure/Packs/MatrixDecompositions/QRPair.cs @@ -0,0 +1,28 @@ +namespace Skeleton.DataStructure.Packs.MatrixDecompositions; + + +/// +/// result of qr decomposition +/// M = QR +/// Q : unitary matrix +/// R : right(upper) triangle matrix +/// +/// +public class QRPair : Pack2 +{ + /// + /// unitary matrix + /// + public TMatrix Q => Item1; + /// + /// right(upper) triangle matrix + /// + public TMatrix R => Item2; + + /// + /// Constructor + /// + /// + /// + public QRPair(TMatrix q, TMatrix r) : base(q, r) { } +} \ No newline at end of file diff --git a/src/DataStructure/Packs/Pack.cs b/src/DataStructure/Packs/Pack.cs new file mode 100644 index 0000000..620ecbb --- /dev/null +++ b/src/DataStructure/Packs/Pack.cs @@ -0,0 +1,124 @@ +namespace Skeleton.DataStructure.Packs; + +/// +/// +/// +/// +/// +/// +public class Pack4 +{ + /// + /// + public Pack4() + { + Item1 = default!; + Item2 = default!; + Item3 = default!; + Item4 = default!; + } + + /// + /// + /// + /// + /// + /// + public Pack4(T1 i1, T2 i2, T3 i3, T4 i4) + { + Item1 = i1; + Item2 = i2; + Item3 = i3; + Item4 = i4; + } + + /// + /// + protected T1 Item1 { get; set; } + + /// + /// + protected T2 Item2 { get; set; } + + /// + /// + protected T3 Item3 { get; set; } + + /// + /// + protected T4 Item4 { get; set; } +} + +/// +/// +/// +/// +/// +public class Pack3 +{ + /// + /// + protected T1 Item1; + + /// + /// + protected T2 Item2; + + /// + /// + protected T3 Item3; + + /// + /// + public Pack3() + { + Item1 = default!; + Item2 = default!; + Item3 = default!; + } + + /// + /// + /// + /// + /// + public Pack3(T1 i1, T2 i2, T3 i3) + { + Item1 = i1; + Item2 = i2; + Item3 = i3; + } +} + +/// +/// +/// +/// +public class Pack2 +{ + /// + /// + protected T1 Item1; + + /// + /// + protected T2 Item2; + + /// + /// + public Pack2() + { + Item1 = default!; + Item2 = default!; + } + + /// + /// + /// + /// + public Pack2(T1 i1, T2 i2) + { + Item1 = i1; + Item2 = i2; + } +} \ No newline at end of file diff --git a/src/DataStructure/Packs/PlainPack.cs b/src/DataStructure/Packs/PlainPack.cs new file mode 100644 index 0000000..f61f923 --- /dev/null +++ b/src/DataStructure/Packs/PlainPack.cs @@ -0,0 +1,78 @@ +namespace Skeleton.DataStructure.Packs; + +/// +/// +/// +/// +public class PlainPack2 : Pack2 +{ + /// + public PlainPack2(T1 a, T2 b) : base(a, b) { } + + /// + /// + public new T1 Item1 + { + get => base.Item1; + set => base.Item1 = value; + } + + /// + /// + public new T2 Item2 + { + get => base.Item2; + set => base.Item2 = value; + } +} + +/// +/// +/// +/// +/// +public class PlainPack3 : Pack3 +{ + /// + public PlainPack3(T1 a, T2 b, T3 c) : base(a, b, c) { } + + /// + /// + public new T1 Item1 => base.Item1; + + /// + /// + public new T2 Item2 => base.Item2; + + /// + /// + public new T3 Item3 => base.Item3; +} + +/// +/// +/// +/// +/// +/// +public class PlainPack4 : Pack4 +{ + /// + public PlainPack4(T1 a, T2 b, T3 c, T4 d) : base(a, b, c, d) { } + + /// + /// + public new T1 Item1 => base.Item1; + + /// + /// + public new T2 Item2 => base.Item2; + + /// + /// + public new T3 Item3 => base.Item3; + + /// + /// + public new T4 Item4 => base.Item4; +} \ No newline at end of file diff --git a/src/DataStructure/Packs/PolynomialDivisionResult.cs b/src/DataStructure/Packs/PolynomialDivisionResult.cs new file mode 100644 index 0000000..a6bdaa5 --- /dev/null +++ b/src/DataStructure/Packs/PolynomialDivisionResult.cs @@ -0,0 +1,19 @@ +using Skeleton.Analysis.AnalyticFunctions.Polynomials.Implements; + +namespace Skeleton.DataStructure.Packs; + +/// +/// +public class PolynomialDivisionResult : Pack2 +{ + /// + public PolynomialDivisionResult(Polynomial q, Polynomial r) : base(q, r) { } + + /// + /// + public Polynomial Quotient => Item1; + + /// + /// + public Polynomial Remainder => Item2; +} \ No newline at end of file diff --git a/src/DataStructure/Packs/QuadraticRoots.cs b/src/DataStructure/Packs/QuadraticRoots.cs new file mode 100644 index 0000000..fa6ed93 --- /dev/null +++ b/src/DataStructure/Packs/QuadraticRoots.cs @@ -0,0 +1,21 @@ +namespace Skeleton.DataStructure.Packs; + +/// +/// +/// +/// +public class QuadraticRoots : Pack2 +{ + /// + public QuadraticRoots(TScalar x1, TScalar x2) : base(x1, x2) { } + + /// + /// + /// + public TScalar X1 => Item1; + /// + /// + /// + public TScalar X2 => Item2; + +} \ No newline at end of file diff --git a/src/DataStructure/Packs/ReducedRowEchelonForm.cs b/src/DataStructure/Packs/ReducedRowEchelonForm.cs new file mode 100644 index 0000000..b556b20 --- /dev/null +++ b/src/DataStructure/Packs/ReducedRowEchelonForm.cs @@ -0,0 +1,27 @@ +namespace Skeleton.DataStructure.Packs; + +/// +/// +/// basic matrix +public class ReducedRowEchelonForm : Pack3 +{ + /// + public ReducedRowEchelonForm(M a, M b, int[] c) + { + Item1 = a; + Item2 = b; + Item3 = c; + } + + /// + /// + public M ReducedForm => Item1; + + /// + /// + public M RowOperation => Item2; + + /// + /// + public int[] Pivots => Item3; +} \ No newline at end of file diff --git a/src/Samples/C22Samples.cs b/src/Samples/C22Samples.cs new file mode 100644 index 0000000..486081c --- /dev/null +++ b/src/Samples/C22Samples.cs @@ -0,0 +1,39 @@ +using System.Numerics; +using Skeleton.Algebra.DimensionProviders; +using Skeleton.Utils.Helpers; + +namespace Skeleton.Samples; +using C22 = Algebra.CategoryOf.OnField.FMatrix; +/// +/// +public static class C22Samples +{ + /// + /// + /// + /// + public static IEnumerable C22AnySamples(this IEnumerable a) + { + (Complex, Complex)[] w = a + .RISamples() + .ShuffleZip() + .ToArray(); + foreach (((Complex, Complex), (Complex, Complex), (Complex, Complex)) info in w.ShuffleZip2()) + yield return new C22( + info.Item1.Item1, info.Item1.Item2, + info.Item2.Item1, info.Item2.Item2 + ); + } + + /// + /// + /// + /// + public static IEnumerable C22InvertibleSamples(this IEnumerable a) + { + foreach (C22 m in a.C22AnySamples()) + if (m.Det().Magnitude > 0.2d) + yield return m; + } + +} \ No newline at end of file diff --git a/src/Samples/C33Samples.cs b/src/Samples/C33Samples.cs new file mode 100644 index 0000000..6d14d6a --- /dev/null +++ b/src/Samples/C33Samples.cs @@ -0,0 +1,51 @@ +using System.Numerics; +using Skeleton.Algebra.DimensionProviders; +using Skeleton.Utils.Helpers; + +namespace Skeleton.Samples; +using C33 = Algebra.CategoryOf.OnField.FMatrix; +/// +/// +public static class C33Samples +{ + /// + /// + /// + /// + public static IEnumerable C33AnySamples(this IEnumerable a) + { + (Complex, Complex, Complex)[] w = a + .RISamples() + .ShuffleZip2() + .ToArray(); + foreach (((Complex, Complex, Complex), (Complex, Complex, Complex), (Complex, Complex, Complex)) info in + w.ShuffleZip2()) + yield return new C33( + info.Item1.Item1, info.Item1.Item2, info.Item1.Item3, + info.Item2.Item1, info.Item2.Item2, info.Item2.Item3, + info.Item3.Item1, info.Item3.Item2, info.Item3.Item3 + ); + } + + /// + /// + /// + /// + public static IEnumerable C33InvertibleSamples(this IEnumerable a) + { + foreach (C33 m in a.C33AnySamples()) + if (m.Det().Magnitude > 0.2d) + yield return m; + } + + /// + /// + public static void Test() + { + int i = 0; + foreach (C33 m in Enumerable + .Range(-10, 21).Select(x => x * 0.1d) + .C33InvertibleSamples()) + Console.WriteLine($"{i++} ^ {m.Det()}"); + } +} \ No newline at end of file diff --git a/src/Samples/C44Samples.cs b/src/Samples/C44Samples.cs new file mode 100644 index 0000000..5874665 --- /dev/null +++ b/src/Samples/C44Samples.cs @@ -0,0 +1,42 @@ +using System.Numerics; +using Skeleton.Algebra.DimensionProviders; +using Skeleton.Utils.Helpers; + +namespace Skeleton.Samples; +using C44 = Algebra.CategoryOf.OnField.FMatrix; +/// +/// +public static class C44Samples +{ + /// + /// + /// + /// + public static IEnumerable C44AnySamples(this IEnumerable a) + { + (Complex, Complex, Complex, Complex)[] w = a + .RISamples() + .ShuffleZip3() + .ToArray(); + foreach (((Complex, Complex, Complex, Complex), (Complex, Complex, Complex, Complex), (Complex, Complex, Complex + , Complex), (Complex, Complex, Complex, Complex)) info in w.ShuffleZip3()) + yield return new C44( + info.Item1.Item1, info.Item1.Item2, info.Item1.Item3, info.Item1.Item4, + info.Item2.Item1, info.Item2.Item2, info.Item2.Item3, info.Item2.Item4, + info.Item3.Item1, info.Item3.Item2, info.Item3.Item3, info.Item3.Item4, + info.Item4.Item1, info.Item4.Item2, info.Item4.Item3, info.Item4.Item4 + ); + } + + /// + /// + /// + /// + public static IEnumerable C44InvertibleSamples(this IEnumerable a) + { + foreach (C44 m in a.C44AnySamples()) + if (m.Det().Magnitude > 0.2d) + yield return m; + } + +} \ No newline at end of file diff --git a/src/Samples/ComplexSamples.cs b/src/Samples/ComplexSamples.cs new file mode 100644 index 0000000..9ad42a0 --- /dev/null +++ b/src/Samples/ComplexSamples.cs @@ -0,0 +1,48 @@ +using System.Numerics; +using Skeleton.Constants; +using Skeleton.Utils.Helpers; + +namespace Skeleton.Samples; + +/// +/// +public static class ComplexSamples +{ + /// + /// + /// + /// + public static IEnumerable PureImSamples(this IEnumerable a) + { + return a + .Shuffle() + .Select(x => x * CC.I); + } + + /// + /// Sample from real and im + /// + public static IEnumerable RISamples(this IEnumerable a) + { + return a + .ShuffleZip() + .Select(x => x.Item1 + x.Item2 * CC.I); + } + + /// + /// Samples on the unit circle + /// + public static IEnumerable UCSamples(this IEnumerable a) + { + return a + .Select(x => Complex.Exp(x * Math.PI * CC.I)); + } + + + /// + /// + public static void SampleTest() + { + foreach (Complex c in RealSamples.ZO002Samples.RISamples()) Console.WriteLine(c); + } +} \ No newline at end of file diff --git a/src/Samples/DiagC2Samples.cs b/src/Samples/DiagC2Samples.cs new file mode 100644 index 0000000..858d6fb --- /dev/null +++ b/src/Samples/DiagC2Samples.cs @@ -0,0 +1,60 @@ +using System.Numerics; +using Skeleton.Algebra.DimensionProviders; +using Skeleton.Utils.Helpers; + +namespace Skeleton.Samples; +using DiagC2 = Algebra.CategoryOf.OnField.FDiagonalMatrix; +/// +/// +public static class DiagC2Samples +{ + /// + /// sample of diag2 whose tr = 0 + /// + public static IEnumerable DiagC2Tr0Samples(this IEnumerable c) + { + return c + .Select(x => new DiagC2(x, -x)); + } + + /// + /// sample of diag2 whose det = 1 + /// + public static IEnumerable DiagC2Det1Samples(this IEnumerable c) + { + return c + .Select(x => new DiagC2(x, 1d / x)); + } + + /// + /// sample of 2x2 diag matrix that is unitary + /// + public static IEnumerable DiagC2UnitarySamples(this IEnumerable r) + { + return r + .UCSamples() + .ShuffleZip() + .Select(x => new DiagC2(x.Item1, x.Item2)); + } + + /// + /// sample of 2x2 diag matrix that is special unitary + /// + public static IEnumerable DiagC2SpecialUnitarySample(this IEnumerable r) + { + return r + .UCSamples() + .Select(x => new DiagC2(x, 1d / x)); + } + + /// + /// sample of 2x2 diag matrix with im part only + /// + public static IEnumerable DiagC2ImSample(this IEnumerable r) + { + return r + .PureImSamples() + .ShuffleZip() + .Select(x => new DiagC2(x.Item1, x.Item2)); + } +} \ No newline at end of file diff --git a/src/Samples/DiagC3Samples.cs b/src/Samples/DiagC3Samples.cs new file mode 100644 index 0000000..d781b20 --- /dev/null +++ b/src/Samples/DiagC3Samples.cs @@ -0,0 +1,63 @@ +using System.Numerics; +using Skeleton.Algebra.DimensionProviders; +using Skeleton.Utils.Helpers; + +namespace Skeleton.Samples; +using DiagC3 = Algebra.CategoryOf.OnField.FDiagonalMatrix; +/// +/// +public static class DiagC3Samples +{ + /// + /// sample of diag3 whose tr = 0 + /// + public static IEnumerable DiagC3Tr0Samples(this IEnumerable c) + { + return c + .ShuffleZip() + .Select(x => new DiagC3(x.Item1, x.Item2, -x.Item1 - x.Item2)); + } + + /// + /// sample of diag3 whose det = 1 + /// + public static IEnumerable DiagC3Det1Samples(this IEnumerable c) + { + return c + .ShuffleZip() + .Select(x => new DiagC3(x.Item1, x.Item2, 1d / (x.Item1 * x.Item2))); + } + + /// + /// sample of 3x3 diag matrix that is unitary + /// + public static IEnumerable DiagC3UnitarySamples(this IEnumerable r) + { + return r + .UCSamples() + .ShuffleZip2() + .Select(x => new DiagC3(x.Item1, x.Item2, x.Item3)); + } + + /// + /// sample of 3x3 diag matrix that is special unitary + /// + public static IEnumerable DiagC3SpecialUnitarySample(this IEnumerable r) + { + return r + .UCSamples() + .ShuffleZip() + .Select(x => new DiagC3(x.Item1, x.Item2, 1d / (x.Item1 * x.Item2))); + } + + /// + /// sample of 3x3 diag matrix with im part only + /// + public static IEnumerable DiagC3ImSample(this IEnumerable r) + { + return r + .PureImSamples() + .ShuffleZip2() + .Select(x => new DiagC3(x.Item1, x.Item2, x.Item3)); + } +} \ No newline at end of file diff --git a/src/Samples/DiagC4Samples.cs b/src/Samples/DiagC4Samples.cs new file mode 100644 index 0000000..25d959c --- /dev/null +++ b/src/Samples/DiagC4Samples.cs @@ -0,0 +1,63 @@ +using System.Numerics; +using Skeleton.Algebra.DimensionProviders; +using Skeleton.Utils.Helpers; + +namespace Skeleton.Samples; +using DiagC4 = Algebra.CategoryOf.OnField.FDiagonalMatrix; +/// +/// +public static class DiagC4Samples +{ + /// + /// sample of diag4 whose tr = 0 + /// + public static IEnumerable DiagC4Tr0Samples(this IEnumerable c) + { + return c + .ShuffleZip2() + .Select(x => new DiagC4(x.Item1, x.Item2, x.Item3, -x.Item1 - x.Item2 - x.Item3)); + } + + /// + /// sample of diag4 whose det = 1 + /// + public static IEnumerable DiagC4Det1Samples(this IEnumerable c) + { + return c + .ShuffleZip2() + .Select(x => new DiagC4(x.Item1, x.Item2, x.Item3, 1d / (x.Item1 * x.Item2 * x.Item3))); + } + + /// + /// sample of 3x3 diag matrix that is unitary + /// + public static IEnumerable DiagC4UnitarySamples(this IEnumerable r) + { + return r + .UCSamples() + .ShuffleZip3() + .Select(x => new DiagC4(x.Item1, x.Item2, x.Item3, x.Item4)); + } + + /// + /// sample of 4x4 diag matrix that is special unitary + /// + public static IEnumerable DiagC4SpecialUnitarySample(this IEnumerable r) + { + return r + .UCSamples() + .ShuffleZip2() + .Select(x => new DiagC4(x.Item1, x.Item2, x.Item3, 1d / (x.Item1 * x.Item2 * x.Item3))); + } + + /// + /// sample of 4x4 diag matrix with im part only + /// + public static IEnumerable DiagC4ImSample(this IEnumerable r) + { + return r + .PureImSamples() + .ShuffleZip3() + .Select(x => new DiagC4(x.Item1, x.Item2, x.Item3, x.Item4)); + } +} \ No newline at end of file diff --git a/src/Samples/LieSU2Samples.cs b/src/Samples/LieSU2Samples.cs new file mode 100644 index 0000000..38c6875 --- /dev/null +++ b/src/Samples/LieSU2Samples.cs @@ -0,0 +1,51 @@ +using System.Numerics; +using Skeleton.Algebra.DimensionProviders; +using Skeleton.Utils.Helpers; + +namespace Skeleton.Samples; +using LieSU2 = Algebra.CategoryOf.FSpecialLieUnitaryMatrix; +using C22 = Algebra.CategoryOf.OnField.FMatrix; + +/// +/// samples +/// +public static class LieSU2Samples +{ + /// + /// sample by log su3 + /// + /// + /// + public static IEnumerable LieSU2LogSU2Samples(this IEnumerable a) + { + return a + .SU2U2FixSamples() + .Select(x => x.Log()); + } + + /// + /// Test + /// + public static void Test() + { + C22[] w = 30 + .ScaleSamples(-1d, 1d) + .LieSU2LogSU2Samples() + .Select(x => x.AsMatrix) + .ToArray(); + int count = 0; + int hcount = 0; + int shcount = 0; + foreach ((C22, C22) x in w.ShuffleZip()) + { + C22 s = x.Item1 * x.Item2 - x.Item2 * x.Item1; + if ((s - s.Hermitian()).IsEqualApprox(C22.Zero)) hcount += 1; + + if ((s + s.Hermitian()).IsEqualApprox(C22.Zero)) shcount += 1; + + count += 1; + } + + Console.WriteLine($"{count} >> h {hcount} >> sh {shcount}"); + } +} \ No newline at end of file diff --git a/src/Samples/LieSU3Samples.cs b/src/Samples/LieSU3Samples.cs new file mode 100644 index 0000000..ec6edde --- /dev/null +++ b/src/Samples/LieSU3Samples.cs @@ -0,0 +1,51 @@ +using System.Numerics; +using Skeleton.Algebra.DimensionProviders; +using Skeleton.Utils.Helpers; + +namespace Skeleton.Samples; +using LieSU3 = Algebra.CategoryOf.FSpecialLieUnitaryMatrix; +using C33 = Algebra.CategoryOf.OnField.FMatrix; + +/// +/// samples +/// +public static class LieSU3Samples +{ + /// + /// sample by log su3 + /// + /// + /// + public static IEnumerable LieSU3LogSU3Samples(this IEnumerable a) + { + return a + .SU3U3FixSamples() + .Select(x => x.Log()); + } + + /// + /// Test + /// + public static void Test() + { + C33[] w = 30 + .ScaleSamples(-1d, 1d) + .LieSU3LogSU3Samples() + .Select(x => x.AsMatrix) + .ToArray(); + int count = 0; + int hcount = 0; + int shcount = 0; + foreach ((C33, C33) x in w.ShuffleZip()) + { + C33 s = x.Item1 * x.Item2 - x.Item2 * x.Item1; + if ((s - s.Hermitian()).IsEqualApprox(C33.Zero)) hcount += 1; + + if ((s + s.Hermitian()).IsEqualApprox(C33.Zero)) shcount += 1; + + count += 1; + } + + Console.WriteLine($"{count} >> h {hcount} >> sh {shcount}"); + } +} \ No newline at end of file diff --git a/src/Samples/LieSU4Samples.cs b/src/Samples/LieSU4Samples.cs new file mode 100644 index 0000000..726f044 --- /dev/null +++ b/src/Samples/LieSU4Samples.cs @@ -0,0 +1,58 @@ +using System.Numerics; +using Skeleton.Algebra.DimensionProviders; +using Skeleton.Utils.Helpers; + +namespace Skeleton.Samples; +using LieSU4 = Algebra.CategoryOf.FSpecialLieUnitaryMatrix; +using SU4 = Algebra.CategoryOf.FSpecialUnitaryMatrix; +using C44 = Algebra.CategoryOf.OnField.FMatrix; + +/// +/// samples +/// +public static class LieSU4Samples +{ + /// + /// sample by log su4 + /// + /// + /// + public static IEnumerable LieSU4LogSU4Samples(this IEnumerable a) + { + + return a + .SU4U4FixSamples() + .Select(x => x.Log()); + } + + /// + /// Test + /// + public static void Test() + { + SU4 a = new SU4(); + LieSU4 loga = a.Log(); + LieSU4 ww = new LieSU4(); + SU4 expww = ww.Exp(); + + C44[] w = 30 + .ScaleSamples(-1d, 1d) + .LieSU4LogSU4Samples() + .Select(x => x.AsMatrix) + .ToArray(); + int count = 0; + int hcount = 0; + int shcount = 0; + foreach ((C44, C44) x in w.ShuffleZip()) + { + C44 s = x.Item1 * x.Item2 - x.Item2 * x.Item1; + if ((s - s.Hermitian()).IsEqualApprox(C44.Zero)) hcount += 1; + + if ((s + s.Hermitian()).IsEqualApprox(C44.Zero)) shcount += 1; + + count += 1; + } + + Console.WriteLine($"{count} >> h {hcount} >> sh {shcount}"); + } +} \ No newline at end of file diff --git a/src/Samples/R2Samples.cs b/src/Samples/R2Samples.cs new file mode 100644 index 0000000..3160e65 --- /dev/null +++ b/src/Samples/R2Samples.cs @@ -0,0 +1,24 @@ +using Skeleton.Algebra.DimensionProviders; +using Skeleton.Utils.Helpers; + +namespace Skeleton.Samples; +using R2 = Algebra.CategoryOf.OnField.FVector; +/// +/// +/// +public static class R2Samples +{ + /// + /// + /// + /// + /// + public static IEnumerable R2AnySamples(this IEnumerable a) + { + (double, double)[] w = a + .ShuffleZip() + .ToArray(); + foreach ((double, double) info in w) + yield return new R2(info.Item1, info.Item2); + } +} \ No newline at end of file diff --git a/src/Samples/R3Samples.cs b/src/Samples/R3Samples.cs new file mode 100644 index 0000000..ffd4c69 --- /dev/null +++ b/src/Samples/R3Samples.cs @@ -0,0 +1,24 @@ +using Skeleton.Algebra.DimensionProviders; +using Skeleton.Utils.Helpers; + +namespace Skeleton.Samples; +using R3 = Algebra.CategoryOf.OnField.FVector; +/// +/// +/// +public static class R3Samples +{ + /// + /// + /// + /// + /// + public static IEnumerable R3AnySamples(this IEnumerable a) + { + (double, double, double)[] w = a + .ShuffleZip2() + .ToArray(); + foreach ((double, double, double) info in w) + yield return new R3(info.Item1, info.Item2, info.Item3); + } +} \ No newline at end of file diff --git a/src/Samples/R4Samples.cs b/src/Samples/R4Samples.cs new file mode 100644 index 0000000..3233266 --- /dev/null +++ b/src/Samples/R4Samples.cs @@ -0,0 +1,24 @@ +using Skeleton.Algebra.DimensionProviders; +using Skeleton.Utils.Helpers; + +namespace Skeleton.Samples; +using R4 = Algebra.CategoryOf.OnField.FVector; +/// +/// +/// +public static class R4Samples +{ + /// + /// + /// + /// + /// + public static IEnumerable R4AnySamples(this IEnumerable a) + { + (double, double, double, double)[] w = a + .ShuffleZip3() + .ToArray(); + foreach ((double, double, double, double) info in w) + yield return new R4(info.Item1, info.Item2, info.Item3, info.Item4); + } +} \ No newline at end of file diff --git a/src/Samples/RealSamples.cs b/src/Samples/RealSamples.cs new file mode 100644 index 0000000..b017967 --- /dev/null +++ b/src/Samples/RealSamples.cs @@ -0,0 +1,133 @@ +using Skeleton.Constants; + +namespace Skeleton.Samples; + +/// +/// +public static class RealSamples +{ + /// + /// first 168 primes + /// + public static IEnumerable PrimeSamples => new double[] + { + 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, + 31, 37, 41, 43, 47, 53, 59, 61, 67, + 71, 73, 79, 83, 89, 97, 101, 103, + 107, 109, 113, 127, 131, 137, 139, 149, + 151, 157, 163, 167, 173, 179, 181, 191, + 193, 197, 199, 211, 223, 227, 229, 233, + 239, 241, 251, 257, 263, 269, 271, 277, + 281, 283, 293, 307, 311, 313, 317, 331, + 337, 347, 349, 353, 359, 367, 373, 379, + 383, 389, 397, 401, 409, 419, 421, 431, + 433, 439, 443, 449, 457, 461, 463, 467, + 479, 487, 491, 499, 503, 509, 521, 523, 541, + 547, 557, 563, 569, 571, 577, 587, 593, 599, + 601, 607, 613, 617, 619, 631, 641, 643, 647, + 653, 659, 661, 673, 677, 683, 691, 701, 709, + 719, 727, 733, 739, 743, 751, 757, 761, 769, + 773, 787, 797, 809, 811, 821, 823, 827, 829, + 839, 853, 857, 859, 863, 877, 881, 883, 887, + 907, 911, 919, 929, 937, 941, 947, 953, 967, + 971, 977, 983, 991, 997 + }; + + /// + /// 100 integers from 0 to 99 + /// + public static IEnumerable Int100Samples => Enumerable + .Range(0, 100) + .Select(x => (double)x); + + /// + /// 200 integers from 0 to 199 + /// + public static IEnumerable Int200Samples => Enumerable + .Range(0, 200) + .Select(x => (double)x); + + /// + /// 500 integers from 0 to 499 + /// + public static IEnumerable Int500Samples => Enumerable + .Range(0, 500) + .Select(x => (double)x); + + /// + /// 500 uniform samples from 0 to 1 + /// + public static IEnumerable ZO002Samples => + Int500Samples + .Select(x => x * 0.002d); + + /// + /// 200 uniform samples from 0 to 1 + /// + public static IEnumerable ZO005Samples => + Int200Samples + .Select(x => x * 0.005d); + + /// + /// 100 uniform samples from 0 to 1 + /// + public static IEnumerable ZO01Samples => + Int100Samples + .Select(x => x * 0.01d); + + /// + /// 168 positive irrational numbers + /// + public static IEnumerable IRSamples => PrimeSamples + .Select(Math.Sqrt); + + /// + /// 'count' samples from 'start' to 'end'. uniformly + /// + public static IEnumerable ScaleSamples(this int count, double start, double end) + { + return Enumerable + .Range(0, count) + .Select(x => start + (end - start) * x / count); + } + + /// + /// multiple of sqrt 2 + /// + public static IEnumerable S2MSample(this int c) + { + return Enumerable + .Range(1, c) + .Select(x => x * CC.Sqrt2); + } + + /// + /// multiple of sqrt 3 + /// + public static IEnumerable S3MSample(this int c) + { + return Enumerable + .Range(1, c) + .Select(x => x * CC.Sqrt3); + } + + /// + /// multiple of pi + /// + public static IEnumerable PIMSample(this int c) + { + return Enumerable + .Range(1, c) + .Select(x => x * Math.PI); + } + + /// + /// multiple of e + /// + public static IEnumerable PESample(this int c) + { + return Enumerable + .Range(1, c) + .Select(x => x * Math.E); + } +} \ No newline at end of file diff --git a/src/Samples/SU2Samples.cs b/src/Samples/SU2Samples.cs new file mode 100644 index 0000000..a9f0187 --- /dev/null +++ b/src/Samples/SU2Samples.cs @@ -0,0 +1,82 @@ +using System.Numerics; +using Skeleton.Algebra.DimensionProviders; +using Skeleton.Algebra.Extensions; + +namespace Skeleton.Samples; +using SU2 = Algebra.CategoryOf.FSpecialUnitaryMatrix; +using C22 = Algebra.CategoryOf.OnField.FMatrix; + +/// +/// Samples of SU2 for tests +/// +public static class SU2Samples +{ + /// + /// Sample by USU* where U is sampled unitary matrix and S is sampled diag special unitary matrix + /// + /// seed + /// + public static IEnumerable SU2EigenCompositeSamples(this IEnumerable a) + { + double[] w = a.ToArray(); + return w + .U2FixSamples() + .Zip(w.DiagC2SpecialUnitarySample()) + .Select(x => new SU2(x.First * x.Second * x.First.Hermitian())); + } + + /// + /// Sample SU2 by fix random invertible matrices + /// + /// + /// + public static IEnumerable SU2U2FixSamples(this IEnumerable a) + { + return a + .C22InvertibleSamples() + .Select(x => new SU2(x)); + } + + + /// + /// test + /// + public static void SpecialUnitaryTest() + { + int count = 0; + int fcount = 0; + foreach (SU2 w in 50.ScaleSamples(-1d, 1d).SU2EigenCompositeSamples()) + { + fcount++; + if (!w.Det().IsEqualApprox(1) || !(w * w.Hermitian()).IsEqualApprox(C22.One)) + { + count++; + Console.WriteLine(w.Representation); + Console.WriteLine($"---->{w.Det()}"); + } + } + + Console.WriteLine($"{fcount}>>{count}"); + } + + /// + /// test + /// + public static void SpecialUnitaryFixTest() + { + int count = 0; + int fcount = 0; + foreach (SU2 w in 50.ScaleSamples(-1d, 1d).SU2U2FixSamples()) + { + fcount++; + if (!w.Det().IsEqualApprox(1) || !(w * w.Hermitian()).IsEqualApprox(C22.One)) + { + count++; + Console.WriteLine(w.Representation); + Console.WriteLine($"---->{w.Det()}"); + } + } + + Console.WriteLine($"{fcount}>>{count}"); + } +} \ No newline at end of file diff --git a/src/Samples/SU3Samples.cs b/src/Samples/SU3Samples.cs new file mode 100644 index 0000000..fdadd03 --- /dev/null +++ b/src/Samples/SU3Samples.cs @@ -0,0 +1,82 @@ +using System.Numerics; +using Skeleton.Algebra.DimensionProviders; +using Skeleton.Algebra.Extensions; + +namespace Skeleton.Samples; +using SU3 = Algebra.CategoryOf.FSpecialUnitaryMatrix; +using C33 = Algebra.CategoryOf.OnField.FMatrix; + +/// +/// Samples of SU3 for tests +/// +public static class SU3Samples +{ + /// + /// Sample by USU* where U is sampled unitary matrix and S is sampled diag special unitary matrix + /// + /// seed + /// + public static IEnumerable SU3EigenCompositeSamples(this IEnumerable a) + { + double[] w = a.ToArray(); + return w + .U3FixSamples() + .Zip(w.DiagC3SpecialUnitarySample()) + .Select(x => new SU3(x.First * x.Second * x.First.Hermitian())); + } + + /// + /// Sample SU3 by fix random invertible matrices + /// + /// + /// + public static IEnumerable SU3U3FixSamples(this IEnumerable a) + { + return a + .C33InvertibleSamples() + .Select(x => new SU3(x)); + } + + + /// + /// test + /// + public static void SpecialUnitaryTest() + { + int count = 0; + int fcount = 0; + foreach (SU3 w in 50.ScaleSamples(-1d, 1d).SU3EigenCompositeSamples()) + { + fcount++; + if (!w.Det().IsEqualApprox(1) || !(w * w.Hermitian()).IsEqualApprox(C33.One)) + { + count++; + Console.WriteLine(w.Representation); + Console.WriteLine($"---->{w.Det()}"); + } + } + + Console.WriteLine($"{fcount}>>{count}"); + } + + /// + /// test + /// + public static void SpecialUnitaryFixTest() + { + int count = 0; + int fcount = 0; + foreach (SU3 w in 50.ScaleSamples(-1d, 1d).SU3U3FixSamples()) + { + fcount++; + if (!w.Det().IsEqualApprox(1) || !(w * w.Hermitian()).IsEqualApprox(C33.One)) + { + count++; + Console.WriteLine(w.Representation); + Console.WriteLine($"---->{w.Det()}"); + } + } + + Console.WriteLine($"{fcount}>>{count}"); + } +} \ No newline at end of file diff --git a/src/Samples/SU4Samples.cs b/src/Samples/SU4Samples.cs new file mode 100644 index 0000000..dfb523b --- /dev/null +++ b/src/Samples/SU4Samples.cs @@ -0,0 +1,83 @@ +using System.Numerics; +using Skeleton.Algebra.DimensionProviders; +using Skeleton.Algebra.Extensions; + +namespace Skeleton.Samples; +using SU4 = Algebra.CategoryOf.FSpecialUnitaryMatrix; +using SU3 = Algebra.CategoryOf.FSpecialUnitaryMatrix; +using C33 = Algebra.CategoryOf.OnField.FMatrix; + +/// +/// Samples of SU4 for tests +/// +public static class SU4Samples +{ + /// + /// Sample by USU* where U is sampled unitary matrix and S is sampled diag special unitary matrix + /// + /// seed + /// + public static IEnumerable SU4EigenCompositeSamples(this IEnumerable a) + { + double[] w = a.ToArray(); + return w + .U4FixSamples() + .Zip(w.DiagC4SpecialUnitarySample()) + .Select(x => new SU4(x.First * x.Second * x.First.Hermitian())); + } + + /// + /// Sample SU3 by fix random invertible matrices + /// + /// + /// + public static IEnumerable SU4U4FixSamples(this IEnumerable a) + { + return a + .C44InvertibleSamples() + .Select(x => new SU4(x, true)); + } + + + /// + /// test + /// + public static void SpecialUnitaryTest() + { + int count = 0; + int fcount = 0; + foreach (SU3 w in 50.ScaleSamples(-1d, 1d).SU3EigenCompositeSamples()) + { + fcount++; + if (!w.Det().IsEqualApprox(1) || !(w * w.Hermitian()).IsEqualApprox(C33.One)) + { + count++; + Console.WriteLine(w.Representation); + Console.WriteLine($"---->{w.Det()}"); + } + } + + Console.WriteLine($"{fcount}>>{count}"); + } + + /// + /// test + /// + public static void SpecialUnitaryFixTest() + { + int count = 0; + int fcount = 0; + foreach (SU3 w in 50.ScaleSamples(-1d, 1d).SU3U3FixSamples()) + { + fcount++; + if (!w.Det().IsEqualApprox(1) || !(w * w.Hermitian()).IsEqualApprox(C33.One)) + { + count++; + Console.WriteLine(w.Representation); + Console.WriteLine($"---->{w.Det()}"); + } + } + + Console.WriteLine($"{fcount}>>{count}"); + } +} \ No newline at end of file diff --git a/src/Samples/U2Samples.cs b/src/Samples/U2Samples.cs new file mode 100644 index 0000000..64a6803 --- /dev/null +++ b/src/Samples/U2Samples.cs @@ -0,0 +1,39 @@ +using Skeleton.Algebra.DimensionProviders; +using Skeleton.Algebra.Extensions; + +namespace Skeleton.Samples; +using U2 = Algebra.CategoryOf.FUnitaryMatrix; +/// +/// +public static class U2Samples +{ + /// + /// generate random unitary matrix by gram schmidt process + /// + public static IEnumerable U2FixSamples(this IEnumerable a) + { + return a + .C22InvertibleSamples() + .Select(x => new U2(x)); + } + + /// + /// + public static void UnitaryTest() + { + int count = 0; + int fcount = 0; + foreach (U2 w in 50.ScaleSamples(-1d, 1d).U2FixSamples()) + { + fcount++; + if (!(w * w.Hermitian()).Det().Magnitude.IsEqualApprox(1)) + { + count++; + Console.WriteLine(w.Representation); + Console.WriteLine($"---->{w.Det()}"); + } + } + + Console.WriteLine($"{fcount}>>{count}"); + } +} \ No newline at end of file diff --git a/src/Samples/U3Samples.cs b/src/Samples/U3Samples.cs new file mode 100644 index 0000000..0c8d73f --- /dev/null +++ b/src/Samples/U3Samples.cs @@ -0,0 +1,39 @@ +using Skeleton.Algebra.DimensionProviders; +using Skeleton.Algebra.Extensions; + +namespace Skeleton.Samples; +using U3 = Algebra.CategoryOf.FUnitaryMatrix; +/// +/// +public static class U3Samples +{ + /// + /// generate random unitary matrix by random invertable + /// + public static IEnumerable U3FixSamples(this IEnumerable a) + { + return a + .C33InvertibleSamples() + .Select(x => new U3(x)); + } + + /// + /// + public static void UnitaryTest() + { + int count = 0; + int fcount = 0; + foreach (U3 w in 50.ScaleSamples(-1d, 1d).U3FixSamples()) + { + fcount++; + if (!(w * w.Hermitian()).Det().Magnitude.IsEqualApprox(1)) + { + count++; + Console.WriteLine(w.Representation); + Console.WriteLine($"---->{w.Det()}"); + } + } + + Console.WriteLine($"{fcount}>>{count}"); + } +} \ No newline at end of file diff --git a/src/Samples/U4Samples.cs b/src/Samples/U4Samples.cs new file mode 100644 index 0000000..bc886cd --- /dev/null +++ b/src/Samples/U4Samples.cs @@ -0,0 +1,39 @@ +using Skeleton.Algebra.DimensionProviders; +using Skeleton.Algebra.Extensions; + +namespace Skeleton.Samples; +using U4 = Algebra.CategoryOf.FUnitaryMatrix; +/// +/// +public static class U4Samples +{ + /// + /// generate random unitary matrix by gram schmidt process + /// + public static IEnumerable U4FixSamples(this IEnumerable a) + { + return a + .C44InvertibleSamples() + .Select(x => new U4(x)); + } + + /// + /// + public static void UnitaryTest() + { + int count = 0; + int fcount = 0; + foreach (U4 w in 50.ScaleSamples(-1d, 1d).U4FixSamples()) + { + fcount++; + if (!(w * w.Hermitian()).Det().Magnitude.IsEqualApprox(1)) + { + count++; + Console.WriteLine(w.Representation); + Console.WriteLine($"---->{w.Det()}"); + } + } + + Console.WriteLine($"{fcount}>>{count}"); + } +} \ No newline at end of file diff --git a/src/Utils/Approximations/ErfApproximation.cs b/src/Utils/Approximations/ErfApproximation.cs new file mode 100644 index 0000000..155a12f --- /dev/null +++ b/src/Utils/Approximations/ErfApproximation.cs @@ -0,0 +1,105 @@ +using Skeleton.Constants; +using Skeleton.Samples; + +namespace Skeleton.Utils.Approximations; + +/// +/// Approximation of Erf error function +/// +public static class ErfApproximation +{ + /// + /// erfinv coeff + /// + private static readonly double[] Ci = new[] + { + "(0,0,0,0,0,0,240,63)", + "(0,0,0,0,0,0,240,63)", + "(171,170,170,170,170,170,242,63)", + "(63,233,147,62,233,147,246,63)", + "(214,91,189,213,91,189,251,63)", + "(124,154,42,215,74,48,1,64)", + "(81,124,23,34,214,106,5,64)", + "(99,101,159,245,0,201,10,64)", + "(150,124,101,0,66,203,16,64)", + "(208,181,155,103,8,26,21,64)" + }.Select(x => x.ByteStringToDouble()).ToArray(); + + /// + /// Err error function + /// max error less than 0.0036128 + /// + /// + /// + public static double Erf(double x) + { + double p = 0.3275911d; + double a1 = 0.254829592d; + double a2 = -0.284496736d; + double a3 = 1.421413741d; + double a4 = -1.453152027d; + double a5 = 1.061405429d; + double t = 1d / (1d + p * x); + double e_x2 = Math.Exp(-x * x); + double t2 = t * t; + double t3 = t2 * t; + double t4 = t2 * t2; + double t5 = t4 * t; + return 1 - (a1 * t + a2 * t2 + a3 * t3 + a4 * t4 + a5 * t5) * e_x2; + } + + /// + /// Inverse of Erf function + /// + /// + /// + public static double ErfInv(double x) + { + if (x < -1 || x > 1) + return double.NaN; + double z = AlgebraConstant.SqrtPi * x / 2d; + double res = 0; + for (int i = 0; i < Ci.Length; i++) + res += Ci[i] / (2d * i + 1d) * Math.Pow(z, 2 * i + 1); + return res; + } + + + private static double c(int k) + { + if (k == 0) + return 1; + double res = 0; + for (int m = 0; m <= k - 1; m++) + res += c(m) * c(k - 1 - m) / ((m + 1) * (2 * m + 1)); + return res; + } + + /// + /// + public static void Test() + { + for (int m = 0; m <= 9; m++) + { + double s = c(m); + Console.WriteLine($"{s.ExactDoubleString()} {s}"); + } + } + + /// + /// + public static void TestErf() + { + double[] u = 100.ScaleSamples(-0.99d, 0.99d).ToArray(); + double[] w = u.Select(ErfInv).ToArray(); + Console.Write("["); + foreach (double x in w) + Console.Write($"{x},"); + Console.WriteLine("]"); + + Console.WriteLine(); + Console.Write("["); + foreach (double x in u) Console.Write($"{x},"); + Console.WriteLine("]"); + } +} \ No newline at end of file diff --git a/src/Utils/Helpers/AlgebraHelper.cs b/src/Utils/Helpers/AlgebraHelper.cs new file mode 100644 index 0000000..c48dbfd --- /dev/null +++ b/src/Utils/Helpers/AlgebraHelper.cs @@ -0,0 +1,152 @@ +using System.Numerics; +using Skeleton.Algebra.DimensionProviders; +using Skeleton.Algebra.Extensions; + +namespace Skeleton.Utils.Helpers; +using C22 = Algebra.CategoryOf.OnField.FMatrix; +using C33 = Algebra.CategoryOf.OnField.FMatrix; +using C44 = Algebra.CategoryOf.OnField.FMatrix; +/// +/// Tools +/// +public static class AlgebraHelper +{ + + /// + /// Fix a complex number + /// + /// + /// + public static Complex Fix(Complex val) + { + return new Complex(Fix(val.Real), Fix(val.Imaginary)); + } + + /// + /// Fix a real number + /// + /// + /// + public static double Fix(double val) + { + if (val.IsEqualApprox(0)) + return 0; + return val; + } + + /// + /// Delta function d(x) = 1 if x = 0; else d(x) = 0 + /// + /// + /// + public static Complex Delta(Complex x) + { + return x.IsEqualApprox(0) ? 1f : 0f; + } + + /// + /// if x less than 0 return 0, if x greater than 1 return 1 else return x + /// + /// + /// + public static double DoubleCut(this double x) + { + return Math.Max(0, Math.Min(1, x)); + } + + /// + /// if x less than 0 return 0, if x greater than max return max else return x + /// + /// + /// + /// + public static double DoubleCut(this double x, double max) + { + return Math.Max(0, Math.Min(max, x)); + } + + /// + /// if x less than 0 return 0, if x greater than 1 return 1 else return x + /// + /// + /// + public static float DoubleCut(this float x) + { + return Math.Max(0, Math.Min(1, x)); + } + + /// + /// l1 error + /// + /// + /// + public static double ManhattanError(this Complex a) + { + return Math.Max(Math.Abs(a.Real), Math.Abs(a.Imaginary)); + } + + + /// + /// power of 2x2 Jordan Block + /// + /// + /// + /// + public static C22 JordanBlock2Power(Complex eigenValue, Complex z) + { + if (eigenValue.IsEqualApprox(0)) + return new C22( + Delta(z), Delta(1 - z), + 0, Delta(z) + ); + return new C22( + Complex.Pow(eigenValue, z), z * Complex.Pow(eigenValue, z - 1), + 0, Complex.Pow(eigenValue, z) + ); + } + + /// + /// power of 3x3 Jordan Block + /// + /// + /// + /// + public static C33 JordanBlock3Power(Complex ev, Complex z) + { + if (ev.IsEqualApprox(0)) + return new C33( + Delta(z), Delta(1 - z), Delta(2 - z), + 0, Delta(z), Delta(1 - z), + 0, 0, Delta(z) + ); + return new C33( + Complex.Pow(ev, z), z * Complex.Pow(ev, z - 1), (z * z - z) * Complex.Pow(ev, z - 2) / 2, + 0, Complex.Pow(ev, z), z * Complex.Pow(ev, z - 1), + 0, 0, Complex.Pow(ev, z) + ); + } + + /// + /// power of 4x4 Jordan Block + /// + /// + /// + /// + public static C44 JordanBlock4Power(Complex ev, Complex z) + { + if (ev.IsEqualApprox(0)) + return new C44( + Delta(z), Delta(1 - z), Delta(2 - z), Delta(3 - z), + 0, Delta(z), Delta(1 - z), Delta(2 - z), + 0, 0, Delta(z), Delta(1 - z), + 0, 0, 0, Delta(z) + ); + return new C44( + Complex.Pow(ev, z), z * Complex.Pow(ev, z - 1), (z * z - z) * Complex.Pow(ev, z - 2) / 2, + (z * z * z - 3 * z * z + 2 * z) * Complex.Pow(ev, z - 3) / 6, + 0, Complex.Pow(ev, z), z * Complex.Pow(ev, z - 1), (z * z - z) * Complex.Pow(ev, z - 2) / 2, + 0, 0, Complex.Pow(ev, z), z * Complex.Pow(ev, z - 1), + 0, 0, 0, Complex.Pow(ev, z) + ); + } +} \ No newline at end of file diff --git a/src/Utils/Helpers/Decorators.cs b/src/Utils/Helpers/Decorators.cs new file mode 100644 index 0000000..066743a --- /dev/null +++ b/src/Utils/Helpers/Decorators.cs @@ -0,0 +1,54 @@ +using Skeleton.Algebra.Extensions; + +namespace Skeleton.Utils.Helpers; + +/// +/// +public static class Decorators +{ + /// + /// + /// + /// + /// + /// + public static T WithTol(double tol, Func f) + { + GeneralExt.SetTol(tol); + T res = f(); + GeneralExt.ResetTol(); + return res; + } + + /// + /// + /// + /// + public static void WithTol(double tol, Action a) + { + GeneralExt.SetTol(tol); + a(); + GeneralExt.ResetTol(); + } + + + /*public static M WithRealFix(Func f) + where D : IDimension + where V : class, IRealVector + where M : class, IRealMatrix + { + M res = f(); + res.Fix(); + return res; + } + + public static M WithComplexFix(Func f) + where D : IDimension + where V : class, IComplexVector + where M : class, IComplexMatrix + { + M res = f(); + res.Fix(); + return res; + }*/ +} \ No newline at end of file diff --git a/src/Utils/Helpers/LinqHelper.cs b/src/Utils/Helpers/LinqHelper.cs new file mode 100644 index 0000000..e7339ad --- /dev/null +++ b/src/Utils/Helpers/LinqHelper.cs @@ -0,0 +1,564 @@ +using System.Diagnostics; +using System.Numerics; +using Skeleton.Algebra; +using Skeleton.Algebra.Extensions; +using Skeleton.DataStructure.Packs; +using Skeleton.Utils.RandomEngines; + +namespace Skeleton.Utils.Helpers; + +/// +/// +public static class LinqHelper +{ + /// + /// + public static int ExecuteCount = 0; + + /// + /// Shuffle the enumerable + /// + public static IEnumerable Shuffle(this IEnumerable a) + { + return a.OrderBy(_ => RandomSource.Rnd.NextDouble()); + } + + /// + /// zip two shuffled copies + /// + public static IEnumerable<(T, T)> ShuffleZip(this IEnumerable a) + { + T[] w = a.ToArray(); + return w.Shuffle().Zip(w.Shuffle()); + } + + /// + /// zip three shuffled copies + /// + public static IEnumerable<(T, T, T)> ShuffleZip2(this IEnumerable a) + { + T[] w = a.ToArray(); + return w.Shuffle().Zip(w.Shuffle(), w.Shuffle()); + } + + /// + /// zip 4 shuffled copies + /// + /// + /// + /// + public static IEnumerable<(T, T, T, T)> ShuffleZip3(this IEnumerable a) + { + T[] w = a.ToArray(); + foreach ((T, (T, T, T)) x in w.Shuffle().Zip(w.ShuffleZip2())) + yield return (x.Item1, x.Item2.Item1, x.Item2.Item2, x.Item2.Item3); + } + + /// + /// Sample by specify the number of samples, with replacement + /// + public static IEnumerable Sample(this IEnumerable a, int count) + { + T[] ac = a.ToArray(); + int len = ac.Length; + for (int i = 0; i <= count; i++) + yield return ac[RandomSource.Rnd.NextInt64(0, len)]; + } + + /// + /// Sample by frequency, without replacement + /// + public static IEnumerable Sample(this IEnumerable a, double freq) + { + return a.Where(x => RandomSource.Rnd.NextDouble() < freq); + } + + /// + /// + /// + /// + /// + public static void ForEach(this IEnumerable q, Action a) + { + foreach (TIn x in q) + a(x); + } + + /// + /// + /// + /// + /// + /// + public static IEnumerable> PairBy(this IEnumerable source, Func selector) + { + HashSet ht = source.ToHashSet(); + foreach (T m1 in ht) + { + T m2 = ht.MaxBy(x => selector(m1, x))!; + yield return new PlainPack2(m1, m2); + } + } + + /// + /// + /// + /// + /// + /// + /// + public static IEnumerable> PairWithBy( + this IEnumerable source, + IEnumerable oth, + Func selector + ) + { + HashSet othC = oth.ToHashSet(); + foreach (T w1 in source) + { + T w2 = othC.MaxBy(x => selector(w1, x))!; + yield return new PlainPack2(w1, w2); + } + } + + /// + /// + /// + /// + /// + /// + /// + public static IEnumerable GroupedOrderBy + ( + this IEnumerable source, + Func key + ) + { + return source.GroupBy(key).OrderBy(g => g.Key).Select(x => x.ToArray()); + } + + /// + /// + /// + /// + /// + /// + /// + public static IEnumerable GroupedThenBy + ( + this IEnumerable source, + Func key + ) + { + foreach (TItem[] g in source) + if (g.Count() <= 1) + yield return g; + else + foreach (TItem[] gx in g.GroupedOrderBy(key)) + yield return gx; + } + + /// + /// + /// + /// + /// + /// + /// + public static IOrderedEnumerable> GroupedOrderByX(this IEnumerable source, + Func key) where TKey : notnull + { + Dictionary> groups = new Dictionary>(); + foreach (TItem item in source) + { + TKey index = key(item); + if (groups.ContainsKey(index)) + groups[index].Add(item); + + else + groups[index] = new OrderedGroup { item }; + } + + return groups + .OrderBy(pair => pair.Key) + .Select(pair => pair.Value) + .OrderBy(_ => 1); + } + + /// + /// + /// + /// + /// + /// + /// + public static IOrderedEnumerable> GroupedThenByX( + this IOrderedEnumerable> source, Func key) where TKey : notnull + { + List> res = new List>(); + foreach (OrderedGroup group in source) + if (group.Count() <= 1) + res.Add(group); + else + res.AddRange(group.GroupedOrderByX(key)); + return res.OrderBy(_ => 1); + } + + /// + /// + /// + /// + /// + /// + /// + public static IEnumerable GroupedMinByX(this IEnumerable source, Func key) + { + return source.GroupBy(key).MinBy(g => g.Key)!; + } + + /// + /// + /// + /// + /// + /// + /// + public static IEnumerable GroupedMinBy(this IEnumerable source, Func key) + { + IEnumerator iter = source.GetEnumerator(); + if (!iter.MoveNext()) + return Array.Empty(); + + TItem current = iter.Current; + TItem[] res = { current }; + if (!iter.MoveNext()) + { + iter.Dispose(); + return res; + } + + TKey minIndex = key(current); + current = iter.Current; + while (iter.MoveNext()) + { + TKey index = key(current)!; + if (!index.Equals(minIndex) && index.Equals(new[] { index, minIndex }.Min())) + { + minIndex = index; + res = new[] { current }; + } + else if (index.Equals(minIndex)) + { + Array.Resize(ref res, res.Length + 1); + res[^1] = current; + } + + current = iter.Current; + } + + TKey fIndex = key(current)!; + if (!fIndex.Equals(minIndex) && fIndex.Equals(new[] { fIndex, minIndex }.Min())) + res = new[] { current }; + if (fIndex.Equals(minIndex)) + { + Array.Resize(ref res, res.Length + 1); + res[^1] = current; + } + + iter.Dispose(); + return res; + } + + /// + /// + /// + /// + /// + /// + public static IEnumerable GroupedApproxMinBy(this IEnumerable source, Func key) + { + IEnumerator iter = source.GetEnumerator(); + if (!iter.MoveNext()) + return Array.Empty(); + TItem current = iter.Current; + TItem[] res = { current }; + if (!iter.MoveNext()) + { + iter.Dispose(); + return res; + } + + double minIndex = key(current); + current = iter.Current; + while (iter.MoveNext()) + { + double index = key(current); + if (!index.IsEqualApprox(minIndex) && index.IsEqualApprox(new[] { index, minIndex }.Min())) + { + minIndex = index; + res = new[] { current }; + } + else if (index.IsEqualApprox(minIndex)) + { + Array.Resize(ref res, res.Length + 1); + res[^1] = current; + } + + current = iter.Current; + } + + double fIndex = key(current); + if (!fIndex.IsEqualApprox(minIndex) && fIndex.IsEqualApprox(new[] { fIndex, minIndex }.Min())) + res = new[] { current }; + if (fIndex.IsEqualApprox(minIndex)) + { + Array.Resize(ref res, res.Length + 1); + res[^1] = current; + } + + iter.Dispose(); + return res; + } + + /// + /// + /// + /// + /// + /// + /// + public static IEnumerable GroupedMaxBy(this IEnumerable source, Func key) + { + IEnumerator iter = source.GetEnumerator(); + if (!iter.MoveNext()) + return Array.Empty(); + TItem current = iter.Current; + TItem[] res = { current }; + if (!iter.MoveNext()) + { + iter.Dispose(); + return res; + } + + TKey maxIndex = key(current); + current = iter.Current; + while (iter.MoveNext()) + { + TKey index = key(current)!; + if (!index.Equals(maxIndex) && index.Equals(new[] { index, maxIndex }.Min())) + { + maxIndex = index; + res = new[] { current }; + } + else if (index.Equals(maxIndex)) + { + Array.Resize(ref res, res.Length + 1); + res[^1] = current; + } + + current = iter.Current; + } + + TKey fIndex = key(current)!; + if (!fIndex.Equals(maxIndex) && fIndex.Equals(new[] { fIndex, maxIndex }.Max())) + res = new[] { current }; + if (fIndex.Equals(maxIndex)) + { + Array.Resize(ref res, res.Length + 1); + res[^1] = current; + } + + iter.Dispose(); + return res; + } + + /// + /// + /// + /// + /// + /// + public static IEnumerable GroupedApproxMaxBy(this IEnumerable source, Func key) + { + IEnumerator iter = source.GetEnumerator(); + if (!iter.MoveNext()) + return Array.Empty(); + TItem current = iter.Current; + TItem[] res = { current }; + if (!iter.MoveNext()) + { + iter.Dispose(); + return res; + } + + double maxIndex = key(current); + current = iter.Current; + while (iter.MoveNext()) + { + double index = key(current); + if (!index.IsEqualApprox(maxIndex) && index.IsEqualApprox(new[] { index, maxIndex }.Min())) + { + maxIndex = index; + res = new[] { current }; + } + else if (index.IsEqualApprox(maxIndex)) + { + Array.Resize(ref res, res.Length + 1); + res[^1] = current; + } + + current = iter.Current; + } + + double fIndex = key(current); + if (!fIndex.Equals(maxIndex) && fIndex.IsEqualApprox(new[] { fIndex, maxIndex }.Max())) + res = new[] { current }; + if (fIndex.Equals(maxIndex)) + { + Array.Resize(ref res, res.Length + 1); + res[^1] = current; + } + + iter.Dispose(); + return res; + } + + /// + /// + public static void OrderTest0() + { + int[] a = new int[10000]; + for (int i = 0; i < 10; i++) + a[i] = i; + Stopwatch f = Stopwatch.StartNew(); + + int w = a + .OrderBy(x => x % 3) + .ThenBy(x => x % 2) + .ThenBy(x => x) + .Select(x => x) + .Sum(); + Console.WriteLine(f.ElapsedMilliseconds); + } + + /// + /// + public static void OrderTest1() + { + int[] a = new int[10000]; + for (int i = 0; i < 10; i++) + a[i] = i; + Stopwatch f = Stopwatch.StartNew(); + int ep = a + .GroupedOrderBy(x => x % 3) + .GroupedThenBy(x => x % 2) + .GroupedThenBy(x => x) + .SelectMany(g => g) + .Sum(); + Console.WriteLine(f.ElapsedMilliseconds); + } + + /// + /// + public static void OrderTest2() + { + int[] a = new int[10000]; + for (int i = 0; i < 10; i++) + a[i] = i; + Stopwatch f = Stopwatch.StartNew(); + int h = a + .GroupedOrderByX(x => x % 3) + .GroupedThenByX(x => x % 2) + .GroupedThenByX(x => x) + .SelectMany(g => g) + .Sum(); + Console.WriteLine(f.ElapsedMilliseconds); + } + + /// + /// + public static void MaxMinTest() + { + int[] a = new int[1000]; + for (int i = 0; i < 1000; i++) + a[i] = i; + Stopwatch f = Stopwatch.StartNew(); + IEnumerable yy = a.GroupedMinByX(x => x % 5); + Console.WriteLine($"{yy.Count()} ---- {f.ElapsedMilliseconds}"); + f = Stopwatch.StartNew(); + IEnumerable xx = a.GroupedMinBy(x => x % 5); + Console.WriteLine($"{xx.Count()} ---- {f.ElapsedMilliseconds}"); + f = Stopwatch.StartNew(); + } + + /// + /// + public static void GroupedMinByTest() + { + HashSet x = new(); + IEnumerable s = x.GroupedMinBy(a => a - 1); + HashSet j = s.ToHashSet(); + Console.WriteLine(j.Count); + foreach (int xx in j) + Console.WriteLine(xx); + } + + /// + /// + public static void ShuffleTest() + { + IEnumerable a = Enumerable.Range(0, 10); + foreach (int x in a.Shuffle()) Console.WriteLine(x); + } + + /// + /// + public static void SampleCountTest() + { + IEnumerable a = Enumerable.Range(0, 10); + int i = 0; + foreach (int x in a.Sample(12)) + Console.WriteLine($"{i++} - {x}"); + } + + /// + /// + public static void SampleFreqTest() + { + IEnumerable a = Enumerable.Range(0, 10); + int i = 0; + foreach (int x in a.Sample(0.3d)) + Console.WriteLine($"{i++} - {x}"); + } + + /// + public class OrderedGroup : List { } + + /// + /// sum first, fix later + /// + /// + /// + /// + public static CategoryOf.FLieUnitaryMatrix LieSum + (this IEnumerable.FLieUnitaryMatrix> es) + { + CategoryOf.OnField.FMatrix res = CategoryOf.OnField.FMatrix.Zero; + CategoryOf.FLieUnitaryMatrix[] aes = es.ToArray(); + foreach (CategoryOf.FLieUnitaryMatrix lie in aes) + res += lie.AsMatrix; + return new CategoryOf.FLieUnitaryMatrix(res); + } + + /// + /// sum first, fix later + /// + /// + /// + /// + public static CategoryOf.FSpecialLieUnitaryMatrix SpLieSum + (this IEnumerable.FSpecialLieUnitaryMatrix> es) => + new(es.LieSum()); + + +} \ No newline at end of file diff --git a/src/Utils/InverseSampling/InverseSampling.cs b/src/Utils/InverseSampling/InverseSampling.cs new file mode 100644 index 0000000..a998bfc --- /dev/null +++ b/src/Utils/InverseSampling/InverseSampling.cs @@ -0,0 +1,40 @@ +using Skeleton.Utils.Approximations; + +namespace Skeleton.Utils.InverseSampling; + +/// +/// Methods for inverse sampling +/// +public static class InverseSampling +{ + /// + /// inverse of normal cdf + /// + /// 0 .. x .. 1 + /// + private static double InverseNormalCDF(double x) + { + return Math.Sqrt(2) * ErfApproximation.ErfInv(2 * x - 1); + } + + /// + /// standard normal distribution + /// + /// 0 .. x .. 1 + /// + public static double StandardNormal(double x) + { + return InverseNormalCDF(x); + } + + /// + /// + /// probability + /// mean value + /// standard deviation + /// + public static double InverseNormal(this double x, double mu, double sigma) + { + return StandardNormal(x) * sigma + mu; + } +} \ No newline at end of file diff --git a/src/Utils/RandomEngines/Normal.cs b/src/Utils/RandomEngines/Normal.cs new file mode 100644 index 0000000..07e0a4e --- /dev/null +++ b/src/Utils/RandomEngines/Normal.cs @@ -0,0 +1,35 @@ +namespace Skeleton.Utils.RandomEngines; + +/// +/// +public static class Normal +{ + /// + /// + public static double? ExNormal; + + /// + /// + /// + public static double Get() + { + if (ExNormal != null) + { + double res = ExNormal.GetValueOrDefault(); + ExNormal = null; + return res; + } + + double rsq, v1, v2; + do + { + v1 = 2d * RandomSource.Get() - 1d; + v2 = 2f * RandomSource.Get() - 1d; + rsq = v1 * v1 + v2 * v2; + } while (rsq > 1d || rsq <= 1E-10d); + + double fac = Math.Sqrt(-2d * Math.Log(rsq) / rsq); + ExNormal = fac * v1; + return fac * v2; + } +} \ No newline at end of file diff --git a/src/Utils/RandomEngines/RandomSource.cs b/src/Utils/RandomEngines/RandomSource.cs new file mode 100644 index 0000000..853f332 --- /dev/null +++ b/src/Utils/RandomEngines/RandomSource.cs @@ -0,0 +1,18 @@ +namespace Skeleton.Utils.RandomEngines; + +/// +/// +public static class RandomSource +{ + /// + /// + public static Random Rnd = new(); + + /// + /// + /// + public static double Get() + { + return Rnd.NextDouble(); + } +} \ No newline at end of file diff --git a/src/Utils/Utils.cs b/src/Utils/Utils.cs new file mode 100644 index 0000000..f4419d3 --- /dev/null +++ b/src/Utils/Utils.cs @@ -0,0 +1,326 @@ +using System.Numerics; +using Skeleton.Algebra.DimensionProviders; +using Skeleton.Algebra.Extensions; + +namespace Skeleton.Utils; +using C22 = Algebra.CategoryOf.OnField.FMatrix; +using C33 = Algebra.CategoryOf.OnField.FMatrix; +using C44 = Algebra.CategoryOf.OnField.FMatrix; + +/// +/// +public static class Utils +{ + /// + /// delta function + /// + /// + /// + public static float Delta(float x) + { + return x.IsEqualApprox(0) ? 1f : 0f; + } + + /// + /// delta function + /// + /// + /// + public static double Delta(double x) + { + return x.IsEqualApprox(0) ? 1f : 0f; + } + + /// + /// delta function + /// + /// + /// + public static Complex Delta(Complex x) + { + return x.IsEqualApprox(0) ? 1f : 0f; + } + + /// + /// power of 2x2 jordan block + /// + /// + /// + /// + public static C22 JordanPower(C22 x, Complex z) + { + if (x[1,1].IsEqualApprox(0)) + return new C22( + Delta(z), Delta(1 - z), + 0, Delta(z) + ); + return new C22( + Complex.Pow(x[1,1], z), z * Complex.Pow(x[1,1], z - 1), + 0, Complex.Pow(x[1,1], z) + ); + } + + /// + /// power of 3x3 jordan block + /// + /// + /// + /// + public static C33 JordanPower(C33 x, Complex z) + { + if (x[1,1].IsEqualApprox(0)) + return new C33( + Delta(z), Delta(1 - z), Delta(2 - z), + 0, Delta(z), Delta(1 - z), + 0, 0, Delta(z) + ); + return new C33( + Complex.Pow(x[1,1], z), z * Complex.Pow(x[1,1], z - 1), (z * z - z) * Complex.Pow(x[1,1], z - 2) / 2, + 0, Complex.Pow(x[1,1], z), z * Complex.Pow(x[1,1], z - 1), + 0, 0, Complex.Pow(x[1,1], z) + ); + } + + /// + /// + /// + /// + /// + public static C44 JordanPower4(Complex ev, Complex z) + { + if (ev.IsEqualApprox(0)) + return new C44( + Delta(z), Delta(1 - z), Delta(2 - z), Delta(3 - z), + 0, Delta(z), Delta(1 - z), Delta(2 - z), + 0, 0, Delta(z), Delta(1 - z), + 0, 0, 0, Delta(z) + ); + return new C44( + Complex.Pow(ev, z), z * Complex.Pow(ev, z - 1), (z * z - z) * Complex.Pow(ev, z - 2) / 2, + (z * z * z - 3 * z * z + 2 * z) * Complex.Pow(ev, z - 3) / 6, + 0, Complex.Pow(ev, z), z * Complex.Pow(ev, z - 1), (z * z - z) * Complex.Pow(ev, z - 2) / 2, + 0, 0, Complex.Pow(ev, z), z * Complex.Pow(ev, z - 1), + 0, 0, 0, Complex.Pow(ev, z) + ); + } + + /// + /// power of 4x4 jordan block + /// + /// + /// + /// + public static C44 JordanPower(C44 x, Complex z) + { + if (x[1,1].IsEqualApprox(0)) + return new C44( + Delta(z), Delta(1 - z), Delta(2 - z), Delta(3 - z), + 0, Delta(z), Delta(1 - z), Delta(2 - z), + 0, 0, Delta(z), Delta(1 - z), + 0, 0, 0, Delta(z) + ); + return new C44( + Complex.Pow(x[1,1], z), z * Complex.Pow(x[1,1], z - 1), (z * z - z) * Complex.Pow(x[1,1], z - 2) / 2, + (z * z * z - 3 * z * z + 2 * z) * Complex.Pow(x[1,1], z - 3) / 6, + 0, Complex.Pow(x[1,1], z), z * Complex.Pow(x[1,1], z - 1), (z * z - z) * Complex.Pow(x[1,1], z - 2) / 2, + 0, 0, Complex.Pow(x[1,1], z), z * Complex.Pow(x[1,1], z - 1), + 0, 0, 0, Complex.Pow(x[1,1], z) + ); + } + + /// + /// + /// + /// + public static double SPow(double c) + { + if ((c % 2).IsEqualApprox(0)) + return 1d; + return -1d; + } + + /// + /// integers from start to end + /// + /// + /// + /// + public static IEnumerable XRange(int start, int end) + { + for (int i = start; i <= end; i++) + yield return i; + } + + /// + /// power of complex, more accurate when order is an integer + /// + /// + /// + /// + /// + public static Complex IterPow(Complex a, int o) + { + if (o < 0) + throw new Exception(""); + if (o == 0) + return Complex.One; + Complex w = IterPow(a, o - 1); + return new Complex(a.Real * w.Real - a.Imaginary * w.Imaginary, a.Real * w.Imaginary + a.Imaginary * w.Real); + } + + /// + /// Do noting + /// + public static void Nop() { } + + /// + /// factorial of x + /// + /// + /// + /// + public static double Factorial(int x) + { + if (x < 0) + throw new Exception("TODO"); + if (x == 0) + return 1d; + return x * Factorial(x - 1); + } + + /// + /// exp of 2x2 matrix + /// + /// + /// + public static C22 C22Exp(C22 matrix) + { + Complex a = matrix[1,1]; + Complex b = matrix[1,2]; + Complex c = matrix[2,1]; + Complex d = matrix[2,2]; + + Complex h = Complex.Sqrt(a * a - 2 * a * d + 4 * b * c + d * d); + Complex s1 = Complex.Exp(0.5 * (h + a + d)); + Complex s2 = Complex.Exp(0.5 * (a + d - h)); + return new C22 + ( + (s1 * (h + a - d) - (a - d - h) * s2) / (2 * h), + (s1 - s2) * b / h, + (s1 - s2) * c / h, + (s1 * (h - a + d) - (d - a - h) * s2) / (2 * h) + ); + } + + /// + /// store double as string + /// + /// + /// + public static string ExactDoubleString(this double x) + { + return "(" + string.Join(",", BitConverter.GetBytes(x)) + ")"; + } + + /// + /// get double from string + /// + /// + /// + public static double ByteStringToDouble(this string s) + { + return BitConverter.ToDouble + ( + s + .Replace("(", "") + .Replace(")", "") + .Split(",") + .Select(x => Convert.ToByte(x)) + .ToArray() + ); + } + + + /// + /// + /// + /// + /// + /// + /// + public static bool Overlap(Vector2 pa, Vector2 sa, Vector2 pb, Vector2 sb) + { + Vector2 pa2 = pa + sa; + Vector2 pb2 = pb + sb; + if (pa.X > pb2.X || pb.X > pa2.X) + return false; + if (pa2.Y < pb.Y || pb2.Y < pa.Y) + return false; + return true; + } + + + /// + /// + /// + /// + /// + /// + public static Action OnFail(this Action test, Action fail) + => () => + { + try + { + test(); + } + catch + { + fail(); + throw; + } + }; + + /// + /// + /// + /// + /// + /// + /// + public static Action OnFail(this Action test, Action fail) + => x => + { + try + { + test(x); + } + catch + { + fail(x); + throw; + } + }; + + /// + /// + /// + /// + /// + /// + /// + /// + public static Action OnFail(this Action test, Action fail) + => (x, y) => + { + try + { + test(x, y); + } + catch (Exception e) + { + fail(x, y); + throw; + } + }; + +} \ No newline at end of file