Matrixmultiplikation

# 🔠 FastShift-Algorithmus in Julia, Autoren: Boris Haase & Math AI 🧠 vom 17.04.2025
using LinearAlgebra, BenchmarkTools, Printf, Random, Base.Threads
const D = 256                                        # Dachwert kleiner oder gleich G
const E = 8                                          # Exponent  etwa größer als fünf
const R = 64                                         # Richtwert  zur Prozesstrennung
const G = 1  << E                                    # Gewünschte  Obergrenze an Bits
const I = BigInt(1)                                  # Eins der  großen ganzen Zahlen
const L = (I << D ) - I                              # Obere bzw. untere Zahlengrenze
const T = AbstractMatrix{BigInt}                     # Matrixtyp großer ganzer Zahlen
const S = G  << 1                                    # Shift-Faktor als  Zweierpotenz
const O = BigInt(0)                                  # Null der  großen ganzen Zahlen

# 🔧 Hilfsfunktionen zur Bestimmung der Division mit Rest und der Matrixmultiplikation
remR(x, t, u, v) =  x < O ? -adjR(-x & u, t, v) : adjR(x & u, t, v)
adjR(x, t, v)    =  x < v ? x : x - t
divR(x, v, s)    = (x + v) >> s

function simd_mul(A::T, B::T, n::Int)::T             # Herkömmliche SIMD-Matrixmultiplikation
    C = zeros(BigInt, n, n)
    @inbounds for i in 1:n                           # Tabellenindizes werden nicht geprüft
        for j in 1:n
            z = O
            @simd ivdep for k in 1:n                 # Keine wechselseitigen Abhängigkeiten
                z += A[i, k] * B[k, j]
            end
            C[i, j] = z
        end
    end
    return C                                         # Rückgabe des Matrixprodukts
end

# 🔡 FastShift zur Matrixmultiplikation als rekursive Funktion für große ganze Zahlen
function fastshift(A::T, B::T, n::Int, s::Int)::T
    n <= R && return simd_mul(A, B, n)               # Aufruf des Standards für kleine Werte

    m = n >> 1                                       # Aufspaltung in vier Teilmatrizen
    M =    1:m                                       # Linker  Index einer Matrix
    N =  m+1:n                                       # Rechter Index einer Matrix
    v = I << (s - 1)                                 # Bezugsgröße für nächsten Nachbarn
    t = I << s                                       # Modulowert  für Division mit Rest
    r = s << 1                                       # Nächste Rekursionsstufe
    u = t -  I                                       # Bitmaske extrahiert die untere Hälfte

    R1 = (B[M, M] .<< s) .+ B[M, N]                  # Obere  Zeile der Matrixkombination
    R2 = (B[N, M] .<< s) .+ B[N, N]                  # Untere Zeile der Matrixkombination

    if n > R                                         # Aufrufen von Threads für große  Werte
        t1 =  @spawn fastshift(A[M, M], R1, m, r)
        t2 =  @spawn fastshift(A[M, N], R2, m, r)
        t3 =  @spawn fastshift(A[N, M], R1, m, r)
        S4 =         fastshift(A[N, N], R2, m, r)    # Matrix für Haupt-Thread
        S1 =  fetch(t1)                              # Ergebnisse für die übrigen Matrizen
        S2 =  fetch(t2)
        S3 =  fetch(t3)
    else
        S1 =         fastshift(A[M, M], R1, m, r)    # Verzicht auf Threads für kleine Werte
        S2 =         fastshift(A[M, N], R2, m, r)
        S3 =         fastshift(A[N, M], R1, m, r)
        S4 =         fastshift(A[N, N], R2, m, r)
    end

    C1 = divR.(S1, v, s   ) .+ divR.(S2, v, s   )    # Divisionen mit Rest wie oben vorgegeben
    C2 = remR.(S1, t, u, v) .+ remR.(S2, t, u, v)
    C3 = divR.(S3, v, s   ) .+ divR.(S4, v, s   )
    C4 = remR.(S3, t, u, v) .+ remR.(S4, t, u, v)

    return [C1 C2; C3 C4]                            # Rückgabe des Matrixprodukts
end

# 🕔 Benchmark für C = A * B mit Korrektheitsprüfung im Vergleich zum Standard von Julia
function main()
    @assert D <= G "Dachwert muss kleiner gleich Obergrenze an Bits sein"
    println("⏱️ Benchmark für das FastShift-Verfahren – m von 1 bis 10")
    println("═════════════════════════════════════════════════════════")
    println("🕸 Threads: $(nthreads())")             # Vier Threads sollen vorgegeben werden

    for m in 1:10
        n  = 1 << m                                  # Shift für Dimension als  Zweierpotenz
        A  = rand(-L:L, n, n)                        # Linke  Matrix zulässiger Zufallswerte
        B  = rand(-L:L, n, n)                        # Rechte Matrix zulässiger Zufallswerte

        println("\n❓ Test: 2^$m × 2^$m → $n × $n")  # Ausgabe der Testergebnisse
        t_fasts = @elapsed C_fasts = fastshift(A, B, n, S + m)
        @printf("🔢 Laufzeit für FastShift: %10.4f Sekunden\n", t_fasts)

        t_julia = @elapsed C_julia = A * B           # Zeitmessung des Standards
        @printf("❕ Laufzeit für Julia A*B: %10.4f Sekunden\n", t_julia)
        println("💯 Gleichheit der Ergebnisse:")     # Kontrolle der Matrixprodukte
        println(C_fasts == C_julia ? "✅ FastShift ≙ Julia" : "❌ FastShift ≠ Julia")
    end
end

# ▶️ Ausführung starten
main()