added MMul4

JuliaBallArithmetic · Aug 26, 2024 · ad8551e · ad8551e
1 parent c738c8f
commit ad8551e
Showing 1 changed file with 99 additions and 0 deletions.
diff --git a/src/types/matrix.jl b/src/types/matrix.jl
@@ -257,6 +257,105 @@ function MMul3(A::Matrix{T}, B::Matrix{T}) where {T <: AbstractFloat}
     BallMatrix(mC, rC)
 end
 
+function MMul4(A::BallMatrix{T}, B::BallMatrix{T}) where {T <: AbstractFloat}
+    m, k = size(A)
+    mA, rA = mid(A), rad(A)
+    mB, rB = mid(B), rad(B)
+
+    C2, rC = setrounding(T, RoundUp) do
+        rC = abs.(mA) * rB + rA * (abs.(mB) + rB)
+        C2 = mA * mB + rC
+        return C2, rC
+    end
+
+    C1 = setrounding(T, RoundDown) do
+        C1 = mA * mB - rC
+        return C1
+    end
+
+    mC, rC = setrounding(T, RoundUp) do
+        mC = (C1 + C2) / 2
+        rC = mC - C1
+        return mC, rC
+    end
+
+    BallMatrix(mC, rC)
+end
+
+function MMul4(A::Matrix{T}, B::BallMatrix{T}) where {T <: AbstractFloat}
+    m, k = size(A)
+    mA = mid(A)
+    mB, rB = mid(B), rad(B)
+
+    C2, rC = setrounding(T, RoundUp) do
+        rC = abs.(mA) * rB
+        C2 = mA * mB + rC
+        return C2, rC
+    end
+
+    C1 = setrounding(T, RoundDown) do
+        C1 = mA * mB - rC
+        return C1
+    end
+
+    mC, rC = setrounding(T, RoundUp) do
+        mC = (C1 + C2) / 2
+        rC = mC - C1
+        return mC, rC
+    end
+
+    BallMatrix(mC, rC)
+end
+
+function MMul4(A::BallMatrix{T}, B::Matrix{T}) where {T <: AbstractFloat}
+    m, k = size(A)
+    mA, rA = mid(A), rad(A)
+    mB = mid(B)
+
+    C2, rC = setrounding(T, RoundUp) do
+        rC = rA * abs.(mB)
+        C2 = mA * mB + rC
+        return C2, rC
+    end
+
+    C1 = setrounding(T, RoundDown) do
+        C1 = mA * mB - rC
+        return C1
+    end
+
+    mC, rC = setrounding(T, RoundUp) do
+        mC = (C1 + C2) / 2
+        rC = mC - C1
+        return mC, rC
+    end
+
+    BallMatrix(mC, rC)
+end
+
+function MMul4(A::Matrix{T}, B::Matrix{T}) where {T <: AbstractFloat}
+    m, k = size(A)
+    mA = mid(A)
+    mB = mid(B)
+
+    C2 = setrounding(T, RoundUp) do
+        C2 = mA * mB
+        return C2
+    end
+
+    C1 = setrounding(T, RoundDown) do
+        C1 = mA * mB
+        return C1
+    end
+
+    mC, rC = setrounding(T, RoundUp) do
+        mC = (C1 + C2) / 2
+        rC = mC - C1
+        return mC, rC
+    end
+
+    BallMatrix(mC, rC)
+end
+
 # As in Revol-Theveny
 # Parallel Implementation of Interval Matrix Multiplication
 # pag. 4