diff --git a/src/Algebra/AdditionalProperties/Extensions/ComplexMatrixExtension.cs b/src/Algebra/AdditionalProperties/Extensions/ComplexMatrixExtension.cs index f44a753..8cff58c 100644 --- a/src/Algebra/AdditionalProperties/Extensions/ComplexMatrixExtension.cs +++ b/src/Algebra/AdditionalProperties/Extensions/ComplexMatrixExtension.cs @@ -172,7 +172,6 @@ public static class ComplexMatrixExtension { for (int i = 1; i <= a.Dim; i++) a[i, i] /= (a[i, i] * Complex.Conjugate(a[i, i])); - } else { @@ -186,7 +185,5 @@ public static class ComplexMatrixExtension } } - - } diff --git a/src/Algebra/FieldStructures/ComplexFieldStructure.cs b/src/Algebra/FieldStructures/ComplexFieldStructure.cs index 74f512f..10f545c 100644 --- a/src/Algebra/FieldStructures/ComplexFieldStructure.cs +++ b/src/Algebra/FieldStructures/ComplexFieldStructure.cs @@ -1,4 +1,5 @@ using System.Numerics; +using Skeleton.Algebra.Extensions; namespace Skeleton.Algebra.FieldStructures; @@ -119,4 +120,11 @@ public class ComplexFieldStructure : FieldStructure /// public override Complex Exp(Complex x) => Complex.Exp(x); + + public override Complex Sign(Complex x) + { + if(x.IsEqualApprox(0)) + return Complex.Zero; + return x / x.Magnitude; + } } diff --git a/src/Algebra/FieldStructures/FieldStructure.cs b/src/Algebra/FieldStructures/FieldStructure.cs index 19ebe09..9b512db 100644 --- a/src/Algebra/FieldStructures/FieldStructure.cs +++ b/src/Algebra/FieldStructures/FieldStructure.cs @@ -1,4 +1,5 @@ using System.Numerics; +using System.Reflection.Metadata; using Skeleton.Algebra.Extensions; using Skeleton.DataStructure.Packs; @@ -81,6 +82,7 @@ public abstract class FieldStructure : FieldStructure /// public virtual TScalar Division(TScalar self, TScalar other) => Multiplication(self, MultiplicationInverse(other)); + /// /// conjugate of the scalar /// @@ -261,23 +263,30 @@ public abstract class FieldStructure : FieldStructure return new QuadraticRoots(x1, x2); } - internal TScalar WilkinsonShift(TScalar a1, TScalar b1, TScalar b2, TScalar a2) + /// + /// + /// + /// + /// + public virtual TScalar Sign(TScalar x) => MultiplicationUnit; + + internal TScalar WilkinsonShift(TScalar a11, TScalar a12, TScalar a21, TScalar a22) { QuadraticRoots es = QuadraticRoot( AdditionUnit, - AdditionInverse(Addition(a1, a2)), - Subtraction(Multiplication(a1, a2), Multiplication(b1, b2)) + AdditionInverse(Addition(a11, a22)), + Subtraction(Multiplication(a11, a22), Multiplication(a12, a21)) ); - TScalar k = a2; - TScalar s = Addition(Addition(Abs(a1), Abs(b1)), Addition(Abs(b2), Abs(a2))); + TScalar k = a22; + TScalar s = Addition(Addition(Abs(a11), Abs(a12)), Addition(Abs(a21), Abs(a22))); if (IsEqualApprox(s, AdditionUnit)) return k; - TScalar q = Multiplication(Division(b1, s), Division(b2, s)); + TScalar q = Multiplication(Division(a12, s), Division(a21, s)); if (IsEqualApprox(q, AdditionUnit)) return k; TScalar p = Division( - Subtraction(Division(a1, s), Division(a2, s)), + Subtraction(Division(a11, s), Division(a22, s)), Addition(MultiplicationUnit, MultiplicationUnit) ); TScalar r = SquareRoot(Addition(Multiplication(p, p), q)); @@ -289,6 +298,36 @@ public abstract class FieldStructure : FieldStructure return k; } + + internal TScalar WilkinsonShiftV2(TScalar a11, TScalar a12, TScalar a21, TScalar a22) + { + TScalar a = a11; + TScalar b = a12; + TScalar c = a22; + TScalar two = Addition(MultiplicationUnit, MultiplicationUnit); + TScalar delta = Division(Subtraction(a, c), two); + TScalar square(TScalar x) => Multiplication(x, x); + TScalar frb = Addition ( + Abs(delta), + SquareRoot( + Addition( + square(delta), + square(Abs(b)) + ) + ) + ); + return Subtraction( + c, + Division( + Multiplication( + Sign(delta), + square(Abs(b)) + ), + frb + ) + ); + } + /// /// homomorphism from complex to this field /// @@ -310,5 +349,6 @@ public abstract class FieldStructure : FieldStructure /// /// public virtual TScalar Exp(TScalar x) => MultiplicationUnit; + } diff --git a/src/Algebra/FieldStructures/RealFieldStructure.cs b/src/Algebra/FieldStructures/RealFieldStructure.cs index e214944..199c565 100644 --- a/src/Algebra/FieldStructures/RealFieldStructure.cs +++ b/src/Algebra/FieldStructures/RealFieldStructure.cs @@ -1,4 +1,5 @@ using System.Numerics; +using Skeleton.Algebra.Extensions; namespace Skeleton.Algebra.FieldStructures; @@ -70,6 +71,13 @@ public class RealFieldStructure : FieldStructure .ToArray() ); + public override double Sign(double x) + { + if (x.IsEqualApprox(0)) + return 0; + return x > 0 ? 1 : -1; + } + /// public override double MaxError(double a) => Math.Abs(a); diff --git a/src/Algebra/Matrices/Matrix.cs b/src/Algebra/Matrices/Matrix.cs index 6d7eb3d..e2e403a 100644 --- a/src/Algebra/Matrices/Matrix.cs +++ b/src/Algebra/Matrices/Matrix.cs @@ -455,12 +455,15 @@ public abstract class Matrix : Matrix @@ -1253,7 +1256,7 @@ public abstract class Matrix : Matrix(u, d); } - internal JDPair InternalUnitaryDecomposition(bool triangle = false) + public JDPair InternalUnitaryDecompositionV1(bool triangle = false) { int maxIters = 5000; double firstTol = 1E-6; @@ -1320,6 +1323,18 @@ public abstract class Matrix : Matrix InternalUnitaryDecomposition(bool triangle = false) + { + JDPair jd = InternalUnitaryDecompositionV1(triangle); + if (IsEqualApprox(jd.J * jd.D * jd.J.Hermitian())) + return jd; + throw new Exception("IDK"); + + } + + /// /// eigen values of the matrix diff --git a/src/Algebra/Matrices/NormalMatrices/FNormalMatrix.cs b/src/Algebra/Matrices/NormalMatrices/FNormalMatrix.cs index e97b7b0..2f659ff 100644 --- a/src/Algebra/Matrices/NormalMatrices/FNormalMatrix.cs +++ b/src/Algebra/Matrices/NormalMatrices/FNormalMatrix.cs @@ -27,12 +27,21 @@ public partial class CategoryOf public FNormalMatrix() { } /// - public FNormalMatrix(IEnumerable es, bool row = true) : base(es, row) => + public FNormalMatrix(IEnumerable es, bool row = true) : base(es, row) + { + if (NeedPreNormalFix()) + PreNormalConstructionFix(); this.PostConstructionFix(); + } + /// - public FNormalMatrix(params TScalar[] es) : base(es) => + public FNormalMatrix(params TScalar[] es) : base(es) + { + if(NeedPreNormalFix()) + PreNormalConstructionFix(); this.PostConstructionFix(); + } /// /// construct by unitary matrix and diagonal matrix @@ -99,6 +108,18 @@ public partial class CategoryOf /// public override IEnumerable EigenValues() => D.Diagonals; + + protected virtual bool NeedPreNormalFix() => Property.IsNormal; + + protected virtual void PreNormalConstructionFix() + { + SelfAdd(Hermitian()); + TScalar two = Structure.Addition(Structure.MultiplicationUnit, Structure.MultiplicationUnit); + for (int i = 1; i <= Dim; i++) + for (int j = 1; j <= Dim; j++) + this[i, j] = Structure.Division(this[i, j], two); + + } } } } diff --git a/src/Algebra/Matrices/NormalMatrices/FSpecialUnitaryMatrix.cs b/src/Algebra/Matrices/NormalMatrices/FSpecialUnitaryMatrix.cs index 72d6a7e..efc2e80 100644 --- a/src/Algebra/Matrices/NormalMatrices/FSpecialUnitaryMatrix.cs +++ b/src/Algebra/Matrices/NormalMatrices/FSpecialUnitaryMatrix.cs @@ -1,6 +1,8 @@ using System.Numerics; using Skeleton.Algebra.AdditionalProperties; +using Skeleton.Algebra.Extensions; using Skeleton.Algebra.FixableObjects; +using Skeleton.DataStructure.Packs.MatrixDecompositions; namespace Skeleton.Algebra; @@ -88,7 +90,7 @@ public static partial class CategoryOf /// 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 index 85a6f79..18b7e12 100644 --- a/src/Algebra/Matrices/NormalMatrices/FUnitaryMatrix.cs +++ b/src/Algebra/Matrices/NormalMatrices/FUnitaryMatrix.cs @@ -110,10 +110,31 @@ public static partial class CategoryOf return new(U, logd); } - /*/// - public override JDPair.FMatrix> EigenDecomposition() => this.UnitaryDecomposition();*/ + /* + * /// + * public override JDPair.FMatrix> EigenDecomposition() => this.UnitaryDecomposition(); + * */ internal override void UnitaryFix(OnField.FMatrix a) => a.UnitaryFix(); + protected override bool NeedPreNormalFix() => !Property.IsUnitary; + protected override void PreNormalConstructionFix() + { + if (Property.IsDiagonal) + { + for (int i = 1; i <= Dim; i++) + this[i, i] /= (this[i, i] * Complex.Conjugate(this[i, i])); + } + else + { + this[1] = this[1].Normalize; + for (int i = 2; i <= Dim; i++) + { + for (int j = 1; j < i; j++) + this[i] -= this[i].ProjectionOn(this[j]); + this[i] = this[i].Normalize; + } + } + } } } diff --git a/src/DataStructure/CacheItem.cs b/src/DataStructure/CacheItem.cs index 7c4ac24..90e1465 100644 --- a/src/DataStructure/CacheItem.cs +++ b/src/DataStructure/CacheItem.cs @@ -33,15 +33,15 @@ public class CacheItem : CacheItem /// construct by provide a method to update the value /// /// - public CacheItem(Func rec) => ProxyCalculator = rec; + public CacheItem(Func rec) => ProxyCalculator = rec; public CacheItem() => ProxyCalculator = c => default; - protected TObject? Value { get; set; } + protected TObject Value { get; set; } /// /// return cached/updated value when reference is not used to calculate another cache item /// - public virtual TObject? Get + public virtual TObject Get { get { @@ -57,14 +57,14 @@ public class CacheItem : CacheItem /// /// provide trace information while calculating /// - public Func ProxyCalculator { get; set; } + public Func ProxyCalculator { get; set; } /// /// trace the dependency and return cached value/updated value /// /// /// - public TObject? GetFrom(CacheItem d) + public TObject GetFrom(CacheItem d) { References.Add(d); d.Dependencies.Add(this); @@ -74,7 +74,7 @@ public class CacheItem : CacheItem /// update the formula of item value /// /// - public void UpdateCalculation(Func val) + public void UpdateCalculation(Func val) { Expire(); foreach (CacheItem item in Dependencies)